空间桁架结构程序设计Fortran

1、空间桁架静力分析程序及算例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)单元的弹性模

2、量AE(2,IN)单元的横截面面积NAE(NE)单元特征类信息存储数组荷载数据PF(4,NCF)外荷载信息数组PF(1,NCF)存放外荷载作用的节点号PF(2,NCF)存放X方向的外荷载PF(3,NCF)存放Y方向的外荷载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*I-1)存放受约束

3、的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、空间桁架结构有限元分析程序源代码!主程序（读入文件，调用总计算程序，输出结果）CHARACTER IDFUT*20,OUTFUT*20WRITE(*,*) Input Data File name:READ (*,*)IDFUTOPEN (11,FILE=IDFUT,STATUS=OLD)WRITE(*,*) Output File name:READ (*,*)OUTFUTOPEN(12,FILE=OUTFUT,STATUS=UNKNOWN)WRITE(12,*)*WRITE(12,*)* Program for Analysis of S

5、pace Trusses *WRITE(12,*)* School of Civil Engineering CSU *WRITE(12,*)* 2012.6.25 Designed By MuZhaoxiang *WRITE(12,*)*WRITE(12,*) WRITE(12,*)* WRITE(12,*)*The Input Data* WRITE(12,*)*WRITE(12,100)READ(11,*)NF,NP,NE,NM,NR,NCF,NDWRITE(12,110)NF,NP,NE,NM,NR,NCF,ND100 FORMAT(6X,The General Information

6、/2X,NF,5X,NP,5X,NE,5X,NM,5X,NR,& 5X,NCF,5X,ND)110 FORMAT(2X,I2,6I7)NPF=NF*NPNDF=ND*NFCALL ANALYSE(NF,NP,NE,NM,NR,NCF,ND,NPF,NDF)END!*!总计算程序SUBROUTINE ANALYSE(NF,NP,NE,NM,NR,NCF,ND,NPF,NDF) DIMENSION X(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(

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

8、),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)ENDIF120 FORMAT(/6X,The Information of Joints/2x,Joint,5X,X,5X,Y,5X,Z)

9、121 FORMAT(1X,I4,3F8.1)130 FORMAT(/6X,The Information of Members/2x,Member,2X,START,4X,END,6X,NAE)131 FORMAT(1X,I4,3I8)140 FORMAT(/6X,The Information of SUPPORTS/2x,Joint,5X,S)141 FORMAT(1X,I4,F8.3)150 FORMAT(/6X,The Information of Sections/4x,E0,8X,A0)151 FORMAT(1X,1PE8.2,F8.4)160 FORMAT(/6X,The Lo

10、ading at Joints/2x,Joint,5X,FX,5X,FY,7X,FZ)161 FORMAT(1X,I4,3F8.2)CALL FLMT(NP,NE,NN,NNM,NR,RR,ND,NF,NDF,ME,IT,LMT)CALL FMAXA(NNM,NE,LMT,MAXA,NWK,NPF,NDF)CALL LP(V,PP,IT,PF,NN,NCF,NF,NP,NPF)CALL CONKB(NP,NE,NM,NWK,ME,X,Y,Z,AE,NAE,LMT,MAXA,CKK,NNM)ISH=1CALL LDLT(CKK,MAXA,NN,ISH,IOUT,NWK,NNM)CALL REBA

11、CK(CKK,V,MAXA,NN,NWK,NNM)CALL DISPLS(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) DO I=1,M DO J=1,N B(J,I)=A(I,J) END DO END DORETURNEND!单元刚度矩阵的形成SUBROUTINE FKE(NP,NE,NM,IE,X,Y,Z,ME,NAE,AE,AKE) DIMENSION X(NP),Y(NP),Z(NP

12、),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/BL AKE(1,1)=C AKE(1,2)=-C AKE(2,1)=-C AKE(2,2)=CRETURNEND!单元坐标转换矩阵 SUBROUTINE FT(IE,NP,NE,X,Y,Z,ME,T

13、) DIMENSION X(NP),Y(NP),Z(NP),ME(2,NE),T(2,6) T=0 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) CX=(X2-X1)/BL CY=(Y2-Y1)/BL CZ=(Z2-Z1)/BL T(1,1)=CX;T(2,4)=CX T(1,2)=CY;T(2,5)=CY T(1,3)=CZ;T(2,6)=CZ RETURNEND!生成单元联系数组LMT SUBROUTINE

14、 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=0 N=0 DO I=1,NP C=0 DO K=1,NR KR=RR(1,K) IF(KR.EQ.I) C=RR(2,K) ENDDO NC=C !NC=0,提取了整数部分 C=C-NC !C=0.*,例如C=0.111 DO J=1,NF C=C*10.0 !例如C=1.21 L=C+0.1 !提取C整数部分,例如L=1,即提取了约束RR(2,K)十分位 !上

15、的数字,这里+0.1是为了防止四舍五入是出现错误 C=C-L IF(L.EQ.0)THEN N=N+1 IT(J,I)=N ELSE IT(J,I)=0 ENDIF ENDDO ENDDO NN=N NNM=NN+1 DO IE=1,NE DO I=1,ND NI=ME(I,IE) DO J=1,NF LMT(I-1)*NF+J,IE)=IT(J,NI) ENDDO ENDDO ENDDO RETURNEND!二维总刚中对角线元地址数组SUBROUTINE FMAXA(NNM,NE,LMT,MAXA,NWK,NPF,NDF) DIMENSION MAXA(NPF),LMT(NDF,NE) MA

16、XA=0;NWK=0 MAXA(1)=1 DO I=2,NNM IP=I-1 IG=IP DO IE=1,NE DO J=1,NDF IF(LMT(J,IE).EQ.IP) THEN DO K=1,NDF IF(LMT(K,IE).GT.0.AND.LMT(K,IE).LE.IG) IG=LMT(K,IE) ENDDO END IF ENDDO ENDDO MAXA(I)= MAXA(I-1)+IP-IG+1 ENDDO NWK= MAXA(NNM)-1 RETURNEND!生成一维存储结构总刚度矩阵SUBROUTINE CONKB(NP,NE,NM,NWK,ME,X,Y,Z,AE,NAE,L

17、MT,MAXA,CKK,NNM) DIMENSION CKK(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=0 DO 10 IE=1,NE TAK=0 CALL FKE(NP,NE,NM,IE,X,Y,Z,ME,NAE,AE,AKE) CALL FT(IE,NP,NE,X,Y,Z,ME,T) CALL MAT(2,6,T,TT) AK=MATMUL(TT,AKE) TAK=MATMUL(AK,T) !总体坐标系

18、下的单元刚度矩阵 DO 220 I=1,6 DO 220 J=1,6 NI=LMT(I,IE) NJ=LMT(J,IE) IF(NJ-NI).GE.0.AND.NI*NJ.GT.0) THEN IJ=MAXA(NJ)+NJ-NI CKK(IJ)=CKK(IJ)+TAK(I,J) ENDIF220 CONTINUE10 CONTINUE RETURN END!生成荷载矩阵SUBROUTINE LP(V,PP,IT,PF,NN,NCF,NF,NP,NPF) DIMENSION V(NN),PP(NPF),IT(NF,NP),PF(4,NCF) V=0 PP=0 DO I=1,NF DO J=1,N

19、P DO K=1,NCF IF(J.EQ.PF(1,K).AND.IT(I,J).NE.0)THENV(IT(I,J)=PF(I+1,K)ENDIF ENDDOENDDO ENDDO DO K=1,NCF DO I=1,NP IF(I.EQ.PF(1,K)THEN PP(NF*(I-1)+1)=PF(2,K) PP(NF*(I-1)+2)=PF(3,K) PP(NF*(I-1)+3)=PF(4,K) ENDIF ENDDO ENDDO RETURNEND!对一维结构总刚度矩阵进行矩阵分解（LDLT） SUBROUTINE LDLT(A,MAXA,NN,ISH,IOUT,NWK,NNM) DIM

20、ENSION A(NWK),MAXA(NNM) IF(NN.EQ.1) RETURN DO 200 N=1,NN KN=MAXA(N) KL=KN+1 KU=MAXA(N+1)-1 KH=KU-KL IF(KH)304,240,210210 K=N-KH IC=0 KLT=KU DO 260 J=1,KH KLT=KLT-1 IC=IC+1 KI=MAXA(K) ND=MAXA(K+1)-KI-1 IF(ND) 260,260,270270 KK=MIN0(IC,ND) C=0.0 DO 280 L=1,KK280 C=C+A(KI+L)*A(KLT+L) A(KLT)=A(KLT)-C260

21、 K=K+1240 K=N B=0.0 DO 300 KK=KL,KU K=K-1 KI=MAXA(K) C=A(KK)/A(KI) IF(ABS(C).LT.1.0E+07) GOTO 290 WRITE(IOUT,2010) N,C STOP290 B=B+C*A(KK)300 A(KK)=C A(KN)=A(KN)-B304 IF(A(KN) 310,310,200310 IF(ISH.EQ.0) GOTO 320 IF(A(KN).EQ.0.0) A(KN)=-1.0E-16 GOTO 200320 WRITE(IOUT,2000) N,A(KN) STOP200 CONTINUE R

22、ETURN2000FORMAT(/, Stop-stiffness matrix not positive definite,/,no positive& pivot for equation,I4,& /, pivot =,E20.10)2010FORMAT(/, Stop-strum sequence check failed + because of multiplier& growth for column & number,I4,/, Multiplier = ,E20.8) END!回代，求得节点位移 SUBROUTINE REBACK(A,V,MAXA,NN,NWK,NNM) D

23、IMENSION A(NWK),V(NN,1),MAXA(NNM) NIP=1 DO IP=1,NIP DO 400 N=1,NN KL=MAXA(N)+1 KU=MAXA(N+1)-1 IF(KU-KL) 400,410,410410 K=N C=0.0 DO 420 KK=KL,KU K=K-1420 C=C+A(KK)*V(K,IP) V(N,IP)=V(N,IP)-C400 CONTINUE DO 480 N=1,NN K=MAXA(N)480 V(N,IP)=V(N,IP)/A(K) IF(NN.EQ.1)RETURN N=NN DO 500 L=2,NN KL=MAXA(N)+1

24、KU=MAXA(N+1)-1 IF(KU-KL) 500,510,510510 K=N DO 520 KK=KL,KU K=K-1520 V(K,IP)=V(K,IP)-A(KK)*V(N,IP)500 N=N-1 ENDDO RETURN END!求解杆件内力、支反力和位移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),

25、 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=0 DO I=1,NP DO J=1,NF LAB=IT(J,I) IF(LAB.EQ.0) THEN DIST(3*(I-1)+J)=0.0 ELSEIF(LAB.GT.0.AND.LAB.LE.NN) THEN DIST(3*(I-1)+J)=FTOOL(LAB) ENDIF ENDDO ENDDO DO IE=1

26、,NE N1=ME(1,IE);N2=ME(2,IE) DO J=1,NF UE(J)=DIST(3*(N1-1)+J) UE(3+J)=DIST(3*(N2-1)+J) ENDDO CALL FT(IE,NP,NE,X,Y,Z,ME,T) CALL FKE(NP,NE,NM,IE,X,Y,Z,ME,NAE,AE,AKE) U=MATMUL(T,UE) FE1=MATMUL(AKE,U) CALL MAT(2,6,T,TT) FE=MATMUL(TT,FE1) DO J=1,NF FF(3*(N1-1)+J)=FF(3*(N1-1)+J)+FE(J) FF(3*(N2-1)+J)=FF(3*(

27、N2-1)+J)+FE(3+J) ENDDO ISW=NAE(IE) AO=AE(2,ISW) SG(IE)=FE1(2) SM(IE)=FE1(2)/AO DO I=1,NPF FF2(I)=FF(I)-PP(I) ENDDO ENDDODO I=1,NPDO J=1,NFLAB=IT(J,I)IF(LAB.EQ.0)THENK=K+1FL(K)=FF2(3*(I-1)+J)ENDIFENDDOENDDO WRITE(12,*) WRITE(12,*)* WRITE(12,*)*The Results of Calculation* WRITE(12,*)* WRITE(12,600)WRI

28、TE(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)600 FORMAT(6X,The Joint Displacement/2x,Joint,6X,X(mm),8X,Y(mm),6X,Z(mm)610 FORMAT(1X,I4,2X,1P

29、3E12.2)620 FORMAT(/6X,The Terminal Forces/2x,Member, 6X,FN(kN),6X,(MPa)630 FORMAT(3X,I4,2X,F8.2,6X,F8.2)640 FORMAT(/6X,The Bearing Force/2x,Joint,8X,X,8X,Y,8X,Z)650 FORMAT(2X,I4,2X,3F10.2) RETURN END3、算例以下图所示空间桁架为例：圆形桁架穹项，其几何尺寸如图（a）所示，整体坐标系原点取在拱顶，集中荷载P作用于拱顶，各杆截面面积A和弹性模量E都相同（取E=210GPa，A=0.04m2）；各杆件及结

30、点编号如图（b）所示。（a）空间桁架立面图（尺寸：m，荷载：kN）（b）空间桁架平面图（结点和杆件编号）（1）输入数据说明第一部分：控制数据The General InformationNFNPNENMNRNCFND313241612第二部分：节点坐标数据The Information of JointsJOINTX(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-20130

31、00第三部分：单元信息数据The Information of MembersMEMBERSTARTENDNAE117121813281429153916310174101841119511110512111612112671137811489115910116101111711121187121197131208131219131221013123111312412131第四部分：约束信息The Information of SupportsJOINTS10.11120.11130.11140.11150.11160.111第五部分：单元截面信息The Information of Sect

32、ionsE0(KPa)A0(m2)210E60.04第六部分：外荷载信息The Loading at JointsJOINTFX(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+

33、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.175-166.66-4.176-166.66-4.177-166.66-4.178-166.66-4.179-

34、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第三部分：支反力结果JOINTFX(KN)FY(KN)FZ(KN)1083.33-295.32-255.7383.33-147.653-255.7

35、383.33147.654083.33295.35255.7383.33147.656255.7383.33-147.65（3）输入及输出文件 输入文件.txt3 13 241 6 1 20 -10 5043.3 -10 2543.3 -10 -250 -10 -50-43.3 -10 -25-43.3 -10 25-12.5 -2 21.6512.5 -221.6525 -2 012.5 -2-21.65-12.5 -2 -21.65-25 -2 0000117121813281429153916310174101841119511110512111612112671137811489115

36、91011610111171112118712119713120813121913122101312311131241213110.11120.11130.11140.11150.11160.111210E60.04130-5000 输出文件.txt * * Program for Analysis of Space Trusses * * School of Civil Engineering CSU * * 2012.6.25 Designed By MuZhaoxiang * * * *The Input Data* * The General Information NF NP NE

37、NM NR NCF ND 3 13 24 1 6 1 2 The Information of Joints Joint X Y Z 1 .0 -10.0 50.0 2 43.3 -10.0 25.0 3 43.3 -10.0 -25.0 4 .0 -10.0 -50.0 5 -43.3 -10.0 -25.0 6 -43.3 -10.0 25.0 7 -12.5 -2.0 21.6 8 12.5 -2.0 21.6 9 25.0 -2.0 .0 10 12.5 -2.0 -21.6 11 -12.5 -2.0 -21.6 12 -25.0 -2.0 .0 13 .0 .0 .0 The Information of Members Member START END NAE 1 1 7 1 2 1 8 1 3 2 8 1 4 2 9 1 5 3 9 1 6 3 10 1 7 4 10 1 8 4 11 1 9 5 11 1 10 5 12 1 11 6 12 1 12 6 7 1 13 7 8 1 14 8 9 1 15 9 10 1 16 10 11 1 17 11 12 1 18 7 12 1 19 7 13 1 20 8 13 1 21 9 13 1 22 10 13 1 23 11 13 1 24 12