空间桁架结构程序设计

1、空间桁架静力分析空间桁架静力分析程序及算例1、变量及数组说明输入数据控制数据NF单个节点的自由度数NP结构离散节点的总数NE结构离散单兀的总数NM结构中单兀/、同的特征数类的总数NR结构受约束节点的总数NCF结构受外荷载作用的节点总数1ND一个单元的节点总数几何数据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,

2、IN)单元的弹性模量AE(2,IN)单兀的横截画面积NAE(NE)单元特征类信息存储数组荷载数据PF(4,NCF)外荷载信息数组PF(1,NCF)存放外荷载作用的节点号PF(2,NCF)WX的外«PF(3,NCF)WY的外«PF(4,NCF)存放Z方向的外荷载输出数据位移DIST(NPF)节点位移数组DIST(NF*I-2)存放I节点X方向的位移DIST(NF*I-1)存放I节点Y方向的位移DIST(NF*I)存放I节点Z方向的位移力SG(NE)单元内力数组SM(NE)单元截面应力数组FL(NF*NR)支座反力数组FL(NF*I-2)存放受约束的I节点X方向的反力FL(NF

3、*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结构矩阵方程的方程总数（去掉约束）NNMNNM=NN+1V(NN)已知节点荷载

4、列阵数组，回代完成后为存放结构位移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,*)

5、*1WRITE(12,*)'*Program for Analysis of Space TrussesWRITE(12,*)'*School of Civil Engineering CSUWRITE(12,*)'*2012.6.25 Designed By MuZhaoxiangWRITE(12,*)*1-1 -WRITE(12,*)''WRITE(12,*)'*1WRITE(12,*)'*The Input Data*WRITE(12,*)'*1WRITE(12,100)READ(11,*)NF,NP,NE,NM,NR,NC

6、F,NDWRITE(12,110)NF,NP,NE,NM,NR,NCF,ND100FORMAT(6X,'TheGeneralInformation72X,'NF',5X,'NP',5X,'NE',5X,'NM',5X,'NR',&5X,'NCF',5X,'ND')110FORMAT(2X,I2,6I7)NPF=NF*NPNDF=ND*NF空间桁架静力分析CALLANALYSE(NF,NP,NE,NM,NR,NCF,ND,NPF,NDF)END!*!总计算程序SUBRO

7、UTINEANALYSE(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

8、)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

9、),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)ENDIF120FORMAT(/6X,'TheInformationofJoints'/2x,'Joint',5X,'X',5X,'Y',5X,'Z')121FORMAT(1X,I4,3F8.1)130FORMAT(/6X,'TheInformationofMembers'/2x,'Member',2X,

10、9;START',4X,'END',6X,'NAE')131FORMAT(1X,I4,3I8)140FORMAT(/6X,'TheInformationofSUPPORTS'/2x,'Joint',5X,'S')141FORMAT(1X,I4,F8.3)150FORMAT(/6X,'TheInformationofSections'/4x,'E0',8X,'A0')151FORMAT(1X,1PE8.2,F8.4)160FORMAT(/6X,'TheL

11、oadingatJoints'/2x,'Joint',5X,'FX',5X,'FY',7X,'FZ')161FORMAT(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(

12、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*!矩阵转置子程序SUBROUTINEMAT(M,N,A,B)DIMENSIONA(M,N),B(N,M)DOI=1,MDOJ=1,NB(J,I)=A(I,J)ENDDOENDDORETURNEND!单元刚度矩阵的形成SUBROUTINEFKE(NP,NE,NM,IE,X,Y,Z,ME,NAE,AE,AKE)DIMENS

13、IONX(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!单元坐标转换矩阵SUBROUTINEFT(IE,NP,NE,X,Y,Z,

14、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!生成单元联系数组LMTSUBROUTINEFLMT(NP,NE,

15、NN,NNM,NR,RR,ND,NF,NDF,ME,IT,LMT)DIMENSIONIT(NF,NP),LMT(NDF,NE),ME(ND,NE),RR(2,NR)NN=0;NNM=0;IT=0;LMT=0N=0DOI=1,NPC=0DOK=1,NRKR=RR(1,K)IF(KR.EQ.I)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"是为了防止四舍五入是出现错误

16、C=C-LIF(L.EQ.0)THENN=N+1IT(J,I)=NELSEIT(J,I)=0ENDIFENDDOENDDONN=NNNM=NN+1DOIE=1,NEDOI=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

17、=1,NDFIF(LMT(J,IE).EQ.IP)THENDOK=1,NDFIF(LMT(K,IE).GT.0.AND.LMT(K,IE).LE.IG)IG=LMT(K,IE)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)

18、,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,6DO220J=1,6NI=LMT(I,IE)NJ=LMT(J,IE)IF(NJ-NI).GE.0.AND.NI*NJ.GT.0)THENIJ=MA

19、XA(NJ)+NJ-NICKK(IJ)=CKK(IJ)+TAK(I,J)ENDIF220CONTINUE10CONTINUERETURNEND!生成荷载矩阵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,NCFIF(J.EQ.PF(1,K).AND.IT(I,J).NE.0)THENV(IT(I,J)=PF(I+1,K)ENDIFENDDOENDDOENDDODOK=1,NCFDOI=1,NPIF(I.EQ.PF(1,K

20、)THENPP(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)IF(NN.EQ.1)RETURNDO200N=1,NNKN=MAXA(N)KL=KN+1KU=MAXA(N+1)-1KH=KU-KLIF(KH)304,240,210210K=N-KHIC=0KLT=KUDO260J=1,KHK

21、LT=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)IF(ABS(C).LT.1.0E+07)GOTO290WRITE(IOUT,2010)N,CSTOP290B=B+C*A(KK)300A(KK)=CA(KN)=A(KN)-B304IF(A(KN)310,310,200

22、310IF(ISH.EQ.0)GOTO320IF(A(KN).EQ.0.0)A(KN)=-1.0E-16GOTO200320WRITE(IOUT,2000)N,A(KN)STOP200CONTINUERETURN2000FORMAT(/,'Stop-stiffnessmatrixnotpositivedefinite',/,'nopositive&pivotforequation',I4,&/,'pivot=',E20.10)2010FORMAT(/,'Stop-strumsequencecheckfailed+becau

23、seofmultiplier&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,410410K=NC=0.0DO420KK=KL,KUK=K-1420C=C+A(KK)*V(K,IP)V(N,IP)=V(N,

24、IP)-C400CONTINUEDO480N=1,NNK=MAXA(N)480V(N,IP)=V(N,IP)/A(K)IF(NN.EQ.1)RETURNN=NNDO500L=2,NNKL=MAXA(N)+1KU=MAXA(N+1)-1IF(KU-KL)500,510,510510K=NDO520KK=KL,KUK=K-1520V(K,IP)=V(K,IP)-A(KK)*V(N,IP)500N=N-1ENDDORETURNEND!求解杆件内力、支反力和位移SUBROUTINEDISPLS(NP,NE,NPF,NM,NN,IT,FTOOL,DIST,AE,NAE,X,Y,Z,PP,FF,SG,SM

25、,ME,&NR,RR,NF)DIMENSIONIT(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)IF(LAB.EQ.0)THENDIST(3*(I-1)+J)=0.0ELS

26、EIF(LAB.GT.0.AND.LAB.LE.NN)THENDIST(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

27、,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)IF(LAB.EQ.0)THENK=K+1FL(K)=FF2(3*(I-1)+J)ENDIFENDDOENDDOWRITE(12,*)' 'WRITE(12,*)'*WRITE(12,

28、*)'*WRITE(12,*)'*WRITE(12,600)，*The Results of Calculation*'，*WRITE(12,610)(I,DIST(3*I-2)*1000,DIST(3*I-1)*1000,&DIST(3*I)*1000, I=1,NP)600610620WRITE(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(6

29、X,'The Joint Displacement'/2x,'Joint',6X,'X(mm)',8X,'Y(mm)',6X,'Z(mm)')FORMAT(1X,I4,2X,1P3E12.2)FORMAT(/6X,' The Terminal Forces'/2x,'Member', 6X,'FN(kN)',6X,'(MPa)')- 9 -空间桁架静力分析630FORMAT(3X,I4,2X,F8.2,6X,F8.2)640FORMAT(/6X,

30、9;TheBearingForce72x,'Joint',8X,'X',8X,'Y',8X,'Z')650FORMAT(2X,I4,2X,3F10.2)RETURNEND3、算例以下图所示空间桁架为例：圆形桁架穹项，其几何尺寸如图(a)所示，整体坐标系原点取在拱顶，集中荷载P作用于拱顶，各杆截面面积A和弹性模量E都相同(取E=210GPa,A=0.04m2);各杆件及结点编号如图(b)所示(a)空间桁架立面图(尺寸：m,荷载：kN)(b)空间桁架平面图(结点和杆件编号)(1)输入数据说明第一部分：控制数据TheGeneralInf

31、ormationNFNPNENMNRNCFND313241612第二部分：节点坐标数据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-201012.5-2-21.6511-12.5-2-21.6512-25-2013000第三部分：单元信息数据TheInformationofMembersMEMBERSTARTENDNAE11712181328142915391631017410184111

32、9511110512111612112671137811489115910116101111711121187121197131208131219131221013123111312412131第四部分：约束信息TheInformationofSupportsJOINTS10.11120.11130.11140.11150.11160.111第五部分：单元截面信息TheInformationofSectionsE0(KPa)A0(m2)210E60.04第六部分：外荷载信息TheLoadingatJointsJOINTFX(KN)FY(KN)FZ(KN)130-5000(2)计算结果第一部分：

33、节点位移结果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-

34、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.175-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.03

35、21.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第三部分：支反力结果JOINTFX(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-1

36、02543.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-2000011712 1813 2814 2915 3916 31017 41018 41119 511110 512111 612112 67113 78114891-13 -空间桁架静力分析1591011610111171112118712119713120813121913122101312311131241213110.11120.11130.11140.11150.11160.111210E

37、60.04130-5000输出文件.txt*ProgramforAnalysisofSpaceTrussesSchoolofCivilEngineeringCSU2012.6.25DesignedByMuZhaoxiang*TheInputData*TheGeneralInformation- 15 -NFNPNENMNR1324NCF1ND2TheInformationofJointsJointXY1.0-10.050.0243.3-10.025.0343.3-10.0-25.04.0-10.0-50.05-43.3-10.0-25.06-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