空间桁架结构程序设计(Fortran)

实用文案空间桁架静力分析程序及算例1、变量及数组说明输入数据NF 单个节点的自由度数NP 结构离散节点的总数NE 结构离散单元的总数控制NM 结构中单元不同的特征数类的总数数据NR 结构受约束节点的总数NCF 结构受外荷载作用的节点总数ND 一个单元的节点总数X(NP) 节点X坐标数组Y(NP) 节点Y坐标数组Z(NP) 节点Z坐标数组单元节点信息存储矩阵几何ME(ND,NE) ME(1,NE)存储杆件始端节点号数据ME(2,NE)储存杆件末端节点号结构约束信息矩阵RR(2,NR) RR(1,NR)存放受有约束的节点号RR(2,NR)存放节点位移约束情况单元 单元特征数类数组特征 AE(2,IN) AE(1,IN)单元的弹性模量数据 AE(2,IN)单元的横截面面积标准文档实用文案NAE(NE) 单元特征类信息存储数组外荷载信息数组PF(1,NCF)存放外荷载作用的节点号荷载PF(4,NCF) PF(2,NCF)存放X方向的外荷载数据PF(3,NCF)存放Y方向的外荷载PF(4,NCF)存放Z方向的外荷载输出数据节点位移数组DIST(NF*I-2) 存放I节点X方向的位移位移 DIST(NPF)DIST(NF*I-1) 存放I节点Y方向的位移DIST(NF*I) 存放I节点Z方向的位移SG(NE) 单元内力数组SM(NE) 单元截面应力数组支座反力数组力FL(NF*I-2)存放受约束的 I节点X方向的反力FL(NF*NR)FL(NF*I-1)存放受约束的 I节点Y方向的反力FL(NF*I)存放受约束的 I节点Z方向的反力中间变量NPF=NF*NP 二维总刚度矩阵的最大行数NDF=ND*NF 一个单元的自由度总数（ 2*3=6）IN 单元特征类总数标准文档实用文案AKE(2,2) 单元在局部坐标系中的刚度局矩阵BL 杆件单元长度T(2,6) 坐标转换矩阵TAK(6,6) 单元在总体坐标系中的刚度矩阵IT(NF,NP) 节点联系数组LMT(NDF,NE) 单元联系数组MAXA(NPF) 结构二维总刚度矩阵主对角元地址数组NWK 结构一维总刚度矩阵的总容量CKK(NWK) 结构一维总刚度矩阵NN 结构矩阵方程的方程总数（去掉约束）NNM NNM=NN+1V(NN) 已知节点荷载列阵数组，回代完成后为存放结构位移PP(NPF) 所有节点荷载列阵数组2、空间桁架结构有限元分析程序源代码主程序（读入文件，调用总计算程序，输出结果）CHARACTERIDFUT*20,OUTFUT*20WRITE(*,*)'InputDataFilename:'READ(*,*)IDFUTOPEN(11,FILE=IDFUT,STATUS='OLD')WRITE(*,*)'OutputFilename:'READ(*,*)OUTFUTOPEN(12,FILE=OUTFUT,STATUS='UNKNOWN')WRITE(12,*)'*****************************************'标准文档实用文案WRITE(12,*)'* ProgramforAnalysisofSpaceTrusses *'WRITE(12,*)'* SchoolofCivilEngineeringCSU *'WRITE(12,*)'* *'WRITE(12,*)'*****************************************'WRITE(12,*)''WRITE(12,*)'*****************************************'WRITE(12,*)'*************TheInputData****************'WRITE(12,*)'*****************************************'WRITE(12,100)READ(11,*)NF,NP,NE,NM,NR,NCF,NDWRITE(12,110)NF,NP,NE,NM,NR,NCF,NDFORMAT(6X,'TheGeneralInformation'/2X,'NF',5X,'NP',5X,'NE',5X,'NM',5X,'NR',&5X,'NCF',5X,'ND')FORMAT(2X,I2,6I7)NPF=NF*NPNDF=ND*NFCALLANALYSE(NF,NP,NE,NM,NR,NCF,ND,NPF,NDF)END!********************************************************************总计算程序SUBROUTINEANALYSE(NF,NP,NE,NM,NR,NCF,ND,NPF,NDF)DIMENSIONX(NP),Y(NP),Z(NP),MM(NE),ME(ND,NE),IT(NF,NP),RR(ND,NR),NAE(NE),&标准文档实用文案AE(1,2),PF(4,NCF),LMT(NDF,NE),MAXA(NPF),CKK(1000),V(NPF),DIST(NPF),&PP(NPF),FF(NPF),SG(NE),SM(NE)READ(11,*)(X(I),Y(I),Z(I),I=1,NP)READ(11,*)(MM(I),ME(1,I),ME(2,I),NAE(I),I=1,NE)READ(11,*)(RR(1,J),RR(2,J),J=1,NR)READ(11,*)(AE(1,J),J=1,2)WRITE(12,120)WRITE(12,121)(I,X(I),Y(I),Z(I),I=1,NP)WRITE(12,130)WRITE(12,131)(MM(I),ME(1,I),ME(2,I),NAE(I),I=1,NE)WRITE(12,140)WRITE(12,141)(INT(RR(1,J)),RR(2,J),J=1,NR)WRITE(12,150)WRITE(12,151)(AE(1,J),J=1,2)IF(NCF/=0)THENREAD(11,*)((PF(I,J),I=1,4),J=1,NCF)WRITE(12,160)WRITE(12,161)(INT(PF(1,J)),PF(2,J),PF(3,J),PF(4,J),J=1,NCF)ENDIFFORMAT(/6X,'TheInformationofJoints'/2x,'Joint',5X,'X',5X,'Y',5X,'Z')FORMAT(1X,I4,3F8.1)130 FORMAT(/6X,'The Information ofMembers'/2x,'Member',2X,'START',4X,'END',6X,'NAE')标准文档实用文案131 FORMAT(1X,I4,3I8)FORMAT(/6X,'TheInformationofSUPPORTS'/2x,'Joint',5X,'S')FORMAT(1X,I4,F8.3)FORMAT(/6X,'TheInformationofSections'/4x,'E0',8X,'A0')FORMAT(1X,1PE8.2,F8.4)FORMAT(/6X,'TheLoadingatJoints'/2x,'Joint',5X,'FX',5X,'FY',7X,'FZ')FORMAT(1X,I4,3F8.2)CALLFLMT(NP,NE,NN,NNM,NR,RR,ND,NF,NDF,ME,IT,LMT)CALLFMAXA(NNM,NE,LMT,MAXA,NWK,NPF,NDF)CALLLP(V,PP,IT,PF,NN,NCF,NF,NP,NPF)CALLCONKB(NP,NE,NM,NWK,ME,X,Y,Z,AE,NAE,LMT,MAXA,CKK,NNM)ISH=1CALLLDLT(CKK,MAXA,NN,ISH,IOUT,NWK,NNM)CALLREBACK(CKK,V,MAXA,NN,NWK,NNM)CALLDISPLS(NP,NE,NPF,NM,NN,IT,V,DIST,AE,NAE,X,Y,Z,PP,FF,SG,SM,ME,NR,RR,NF)END!********************************************************************矩阵转置子程序SUBROUTINE MAT(M,N,A,B)DIMENSION A(M,N),B(N,M)DOI=1,MDOJ=1,NB(J,I)=A(I,J)标准文档实用文案ENDDOENDDORETURNEND单元刚度矩阵的形成SUBROUTINE FKE(NP,NE,NM,IE,X,Y,Z,ME,NAE,AE,AKE)DIMENSION X(NP),Y(NP),Z(NP),ME(2,NE),NAE(NE),AE(2,NM),AKE(2,2)N1=ME(1,IE)N2=ME(2,IE)X1=X(N1);Y1=Y(N1);Z1=Z(N1)X2=X(N2);Y2=Y(N2);Z2=Z(N2)BL=SQRT((X2-X1)**2+(Y2-Y1)**2+(Z2-Z1)**2)NMI=NAE(IE)E0=AE(1,NMI);A0=AE(2,NMI)C=E0*A0/BLAKE(1,1)=CAKE(1,2)=-CAKE(2,1)=-CAKE(2,2)=CRETURNEND标准文档实用文案单元坐标转换矩阵SUBROUTINE FT(IE,NP,NE,X,Y,Z,ME,T)DIMENSIONX(NP),Y(NP),Z(NP),ME(2,NE),T(2,6)T=0N1=ME(1,IE);N2=ME(2,IE)X1=X(N1);Y1=Y(N1);Z1=Z(N1)X2=X(N2);Y2=Y(N2);Z2=Z(N2)BL=SQRT((X2-X1)**2+(Y2-Y1)**2+(Z2-Z1)**2)CX=(X2-X1)/BLCY=(Y2-Y1)/BLCZ=(Z2-Z1)/BLT(1,1)=CX;T(2,4)=CXT(1,2)=CY;T(2,5)=CYT(1,3)=CZ;T(2,6)=CZRETURNEND生成单元联系数组LMTSUBROUTINE FLMT(NP,NE,NN,NNM,NR,RR,ND,NF,NDF,ME,IT,LMT)DIMENSION IT(NF,NP),LMT(NDF,NE),ME(ND,NE),RR(2,NR)NN=0;NNM=0;IT=0;LMT=0N=0DOI=1,NP标准文档实用文案C=0DOK=1,NRKR=RR(1,K) C=RR(2,K)ENDDONC=C !NC=0,提取了整数部分C=C-NC !C=0.***, 例如C=0.111DOJ=1,NFC=C*10.0 !例如C=1.21L=C+0.1 !提取C整数部分,例如L=1,即提取了约束 RR(2,K)十分位上的数字,这里"+0.1"是为了防止四舍五入是出现错误C=C-LN=N+1IT(J,I)=NELSEIT(J,I)=0ENDIFENDDOENDDONN=NNNM=NN+1DOIE=1,NE标准文档实用文案DOI=1,NDNI=ME(I,IE)DOJ=1,NFLMT((I-1)*NF+J,IE)=IT(J,NI)ENDDOENDDOENDDORETURNEND二维总刚中对角线元地址数组SUBROUTINEFMAXA(NNM,NE,LMT,MAXA,NWK,NPF,NDF)DIMENSIONMAXA(NPF),LMT(NDF,NE)MAXA=0;NWK=0MAXA(1)=1DOI=2,NNMIP=I-1IG=IPDOIE=1,NEDOJ=1,NDFIF(LMT(J,IE).EQ.IP)THENDOK=1,NDF标准文档实用文案ENDDOENDIFENDDOENDDOMAXA(I)=MAXA(I-1)+IP-IG+1ENDDONWK=MAXA(NNM)-1RETURNEND生成一维存储结构总刚度矩阵SUBROUTINECONKB(NP,NE,NM,NWK,ME,X,Y,Z,AE,NAE,LMT,MAXA,CKK,NNM)DIMENSIONCKK(NWK),X(NP),Y(NP),Z(NP),AE(2,NM),NAE(NE),LMT(6,NE),ME(2,NE),&MAXA(NNM),AK(6,2),AKE(2,2),T(2,6),TT(6,2),TAK(6,6)CKK=0DO10IE=1,NETAK=0CALLFKE(NP,NE,NM,IE,X,Y,Z,ME,NAE,AE,AKE)CALLFT(IE,NP,NE,X,Y,Z,ME,T)CALLMAT(2,6,T,TT)AK=MATMUL(TT,AKE)TAK=MATMUL(AK,T) !总体坐标系下的单元刚度矩阵DO220I=1,6标准文档实用文案DO220J=1,6NI=LMT(I,IE)NJ=LMT(J,IE)IJ=MAXA(NJ)+NJ-NICKK(IJ)=CKK(IJ)+TAK(I,J)ENDIFCONTINUECONTINUERETURNEND生成荷载矩阵SUBROUTINELP(V,PP,IT,PF,NN,NCF,NF,NP,NPF)DIMENSIONV(NN),PP(NPF),IT(NF,NP),PF(4,NCF)V=0PP=0DOI=1,NFDOJ=1,NPDOK=1,NCFV(IT(I,J))=PF(I+1,K)ENDIFENDDO标准文档实用文案ENDDOENDDODOK=1,NCFDOI=1,NPPP(NF*(I-1)+1)=PF(2,K)PP(NF*(I-1)+2)=PF(3,K)PP(NF*(I-1)+3)=PF(4,K)ENDIFENDDOENDDORETURNEND!对一维结构总刚度矩阵进行矩阵分解（ LDLT）SUBROUTINELDLT(A,MAXA,NN,ISH,IOUT,NWK,NNM)DIMENSIONA(NWK),MAXA(NNM)DO200N=1,NNKN=MAXA(N)KL=KN+1KU=MAXA(N+1)-1KH=KU-KL标准文档实用文案IF(KH)304,240,210210K=N-KHIC=0KLT=KUDO260J=1,KHKLT=KLT-1IC=IC+1KI=MAXA(K)ND=MAXA(K+1)-KI-1IF(ND)260,260,270270KK=MIN0(IC,ND)C=0.0DO280L=1,KK280C=C+A(KI+L)*A(KLT+L)A(KLT)=A(KLT)-C260K=K+1240K=NB=0.0DO300KK=KL,KUK=K-1KI=MAXA(K)C=A(KK)/A(KI)标准文档实用文案WRITE(IOUT,2010)N,CSTOPB=B+C*A(KK)A(KK)=CA(KN)=A(KN)-BIF(A(KN))310,310,200WRITE(IOUT,2000)N,A(KN)STOPCONTINUERETURN2000 FORMAT(//,'Stop-stiffnessmatrixnotpositivedefinite',//,'nopositive&pivotforequation',I4,&//,'pivot=',E20.10)2010 FORMAT(//,'Stop-strumsequencecheckfailed+becauseofmultiplier&growthforcolumn&number',I4,//,'Multiplier=',E20.8)END回代，求得节点位移SUBROUTINEREBACK(A,V,MAXA,NN,NWK,NNM)DIMENSIONA(NWK),V(NN,1),MAXA(NNM)标准文档实用文案NIP=1DOIP=1,NIPDO400N=1,NNKL=MAXA(N)+1KU=MAXA(N+1)-1IF(KU-KL)400,410,410410 K=NC=0.0DO420KK=KL,KUK=K-1420 C=C+A(KK)*V(K,IP)V(N,IP)=V(N,IP)-CCONTINUEDO480N=1,NNK=MAXA(N)480 V(N,IP)=V(N,IP)/A(K)N=NNDO500L=2,NNKL=MAXA(N)+1KU=MAXA(N+1)-1IF(KU-KL)500,510,510510 K=N标准文档实用文案DO520KK=KL,KUK=K-1520 V(K,IP)=V(K,IP)-A(KK)*V(N,IP)500 N=N-1ENDDORETURNEND求解杆件内力、支反力和位移SUBROUTINE DISPLS(NP,NE,NPF,NM,NN,IT,FTOOL,DIST,AE,NAE,X,Y,Z,PP,FF,SG,SM,ME,&NR,RR,NF)DIMENSION IT(3,NP),DIST(NPF),FTOOL(NPF),X(NP),Y(NP),Z(NP),T(2,6),TT(6,2),AE(2,NM),&ME(2,NE),NAE(NE),UE(6),U(2),AKE(2,2),FE1(2),FE(6),FF(NPF),PP(NPF),&SG(NE),SM(NE),FF2(NPF),RR(2,NR),FL(3*NR)SG=0;SM=0;FF=0;FF2=0DOI=1,NPDOJ=1,NFLAB=IT(J,I)DIST(3*(I-1)+J)=0.0标准文档实用文案DIST(3*(I-1)+J)=FTOOL(LAB)ENDIFENDDOENDDODOIE=1,NEN1=ME(1,IE);N2=ME(2,IE)DOJ=1,NFUE(J)=DIST(3*(N1-1)+J)UE(3+J)=DIST(3*(N2-1)+J)ENDDOCALLFT(IE,NP,NE,X,Y,Z,ME,T)CALLFKE(NP,NE,NM,IE,X,Y,Z,ME,NAE,AE,AKE)U=MATMUL(T,UE)FE1=MATMUL(AKE,U)CALLMAT(2,6,T,TT)FE=MATMUL(TT,FE1)DOJ=1,NFFF(3*(N1-1)+J)=FF(3*(N1-1)+J)+FE(J)FF(3*(N2-1)+J)=FF(3*(N2-1)+J)+FE(3+J)ENDDOISW=NAE(IE)AO=AE(2,ISW)SG(IE)=FE1(2)标准文档实用文案SM(IE)=FE1(2)/AODOI=1,NPFFF2(I)=FF(I)-PP(I)ENDDOENDDODOI=1,NPDOJ=1,NFLAB=IT(J,I)K=K+1FL(K)=FF2(3*(I-1)+J)ENDIFENDDOENDDOWRITE(12,*)''WRITE(12,*)'****************************************'WRITE(12,*)'*********TheResultsofCalculation**********'WRITE(12,*)'****************************************'WRITE(12,600)WRITE(12,610)(I,DIST(3*I-2)*1000,DIST(3*I-1)*1000,&DIST(3*I)*1000,I=1,NP)WRITE(12,620)WRITE(12,630)(IE,SG(IE),SM(IE)/1000,IE=1,NE)标准文档实用文案WRITE(12,640)WRITE(12,650)(INT(RR(1,I)),FL(3*I-2),FL(3*I-1),FL(3*I),I=1,NR)FORMAT(6X,'TheJointDisplacement'/2x,'Joint',6X,'X(mm)',8X,'Y(mm)',6X,'Z(mm)')FORMAT(1X,I4,2X,1P3E12.2)620 FORMAT(//6X,'TheTerminalForces'/2x,'Member',6X,'FN(kN)',6X,' σ(MPa)')630 FORMAT(3X,I4,2X,F8.2,6X,F8.2)FORMAT(//6X,'TheBearingForce'/2x,'Joint',8X,'X',8X,'Y',8X,'Z')FORMAT(2X,I4,2X,3F10.2)RETURNEND3、算例以下图所示空间桁架为例：圆形桁架穹项，其几何尺寸如图（ a）所示，整体坐标系原点取在拱顶，集中荷载 P作用于拱顶，各杆截面面积 A和弹性模量E都相同（取E=210GPa，A=0.04m2）；各杆件及结点编号如图（ b）所示。（a）空间桁架立面图（尺寸： m，荷载：kN）标准文档实用文案（b）空间桁架平面图（结点和杆件编号）（1）输入数据说明第一部分：控制数据TheGeneralInformationNFNPNENMNRNCFND313241612第二部分：节点坐标数据TheInformationofJointsJOINTX(m)Y(m)Z(m)10-1050243.3-1025343.3-10-2540-10-505-43.3-10-256-43.3-10257-12.5-221.65812.5-221.65925-20标准文档实用文案1012.5-2-21.6511-12.5-2-21.6512-25-2013000第三部分：单元信息数据TheInformationofMembersMEMBERSTARTENDNAE11712181328142915391631017410184111951111051211161211267113781148911591011610111标准文档实用文案1711121187121197131208131219131221013123111312412131第四部分：约束信息 TheInformationofSupportsJOINT S10.11120.11130.11140.11150.11160.111第五部分：单元截面信息TheInformationofSectionsE0(KPa)A0(m2)210E60.04第六部分：外荷载信息TheLoadingatJointsJOINTFX(KN)FY(KN)FZ(KN)130-5000标准文档实用文案（2）计算结果第一部分：节点位移结果JOINTX(mm)Y(mm)Z(mm)10.00E+000.00E+000.00E+0020.00E+000.00E+000.00E+0030.00E+000.00E+000.00E+0040.00E+000.00E+000.00E+0050.00E+000.00E+000.00E+0060.00E+000.00E+000.00E+007-1.27E+003.25E+002.19E+0081.27E+003.25E+002.19E+0092.53E+003.25E+00-5.25E-08101.27E+003.25E+00-2.19E+0011-1.27E+003.25E+00-2.19E+0012-2.53E+003.25E+001.93E-0713-3.38E-07-6.75E+015.95E-07第二部分：单元内力及应力结果MEMBERFN(kN)(MPa)1-166.66-4.172-166.66-4.173-166.66-4.174-166.65-4.17标准文档实用文案5-166.66-4.176-166.66-4.177-166.66-4.178-166.66-4.179-166.66-4.1710-166.65-4.1711-166.66-4.1712-166.66-4.1713851.0321.2814851.0121.2815851.0121.2816851.0321.2817851.0121.2818851.0121.2819-1044.98-26.1220-1044.98-26.1221-1044.98-26.1222-1044.98-26.1223-1044.98-26.1224-1044.98-26.12第三部分：支反力结果JOINT FX(KN) FY(KN) FZ(KN)标准文档实用文案1083.33-295.32-255.7383.33-147.653-255.7383.33147.654083.33295.35255.7383.33147.656255.7383.33-147.65（3）输入及输出文件①输入文件.txt3132416120-105043.3-102543.3-10-250-10-50-43.3-10-25-43.3-1025-12.5-221.6512.5-221.6525-2012.5-2-21.65-12.5-2-21.65-25-20000117121813281429153916310174101841115111512161211267113781标准文档实用文案14891159101161011117111211871211971312081312191312210131231113124121310.1110.1110.1110.1110.1110.111210E60.0413 0 -500 0②输出文件.txt****************************************** ProgramforAnalysisofSpaceTrusses ** SchoolofCivilEngineeringCSU ** ************************************************************************************************TheInputData*********************************************************TheGeneralInformationNF NP NE NM NR NCF ND3 13 24 1 6 1 2TheInformationofJointsJointXYZ1.0-10.050.0243.3-10.025.0343.3-10.0-25.04.0-10.0-50.05-43.3-10.0-25.0标准文档实用文案6-43.3-10.025.07-12.5-2.021.6812.5-2.021.6925.0-2.0.01012.5-2.0-21.611-12.5-2.0-21.612-25.0-2.0.013.0.0.0TheInformationofMembersMemberSTARTENDNAE11712181328142915391631017410184111951111