Fortran语言编写的有限元结构程序

1、算例一计算简图及结果输出用平面刚架静力计算程序下图结构的内力。各杆EA ，EI 相同。已知： EA=4.0 106 KN,EI=1.6 104 KN m2计算简图如下：.（ 1）输入原始数据控制参数3,5,8,7,1,2（ NE,NJ,N,NW,NPJ,NPF）0.0,0. 0,0,00.0,4. 0,1,2,3结点坐标集结点未知量编号0.0,4. 0,1,2,44.0,4. 0,5,6,74.0,0. 0,0,0,81,2,4.0E+06,1.6E+04单元杆端结点编号及单元EA、EI3,4,4.0E+06,1.6E+045,4,4.0E+06,1.6E+04结点荷载7. 0,-15.非结点

2、荷载1. 0,2. 0,2. 0,-2. 0,1. 0,4. 0,-25. 0（ 2）输出结果NE=3 NJ=5 N=8NW=7 NPJ=1 NPF=2NODEXYXXYYZZ10.00000.000000020.00004.000012330.00004.000012444.00004.000056754.00000.0000008ELEMENTNODE-INODE-JEAEI1120.400000E+070.160000E+052340.400000E+070.160000E+053540.400000E+070.160000E+05CODEPX-PY-PM7.-15.0000ELEMEN

3、TINDAQ1.2.2.0000-18.00002.1.4.0000-25.0000NODEUVCETA10.000000E+000.000000E+000.000000E+002-0.221743E-02-0.464619E-04-0.139404E-023-0.221743E-02-0.464619E-040.357876E-024-0.222472E-02 -0.535381E-04 -0.298554E-0250.000000E+000.000000E+000.658499E-03ELEMENTNQM1N1=46.4619Q1=10.7119M1=-6.8477N2=-46.4619Q

4、2=7.2881M2=0.00002N1=7.2881Q1=46.4619M1=0.0000N2=-7.2881Q2=53.5381M2=14.15233N1=53.5381Q1=7.2881M1=0.0000N2=-53.5381Q2=-7.2881M2=-29.1523算例二计算简图及结果输出用平面刚架静力计算程序下图结构的内力。已知：桁架单元的抗拉刚度为EA=2.0106 KN，平面刚架单元的抗拉刚度为已知：EA=4.0106KN，抗弯刚度为EI=1.84104KNm 2 。计算简图如下：（ 1）输入原始数据控制参数5,4,8,7,0（NE,NJ,N,NW,NPJ,NPF）0.0,0.0

5、,0,04.0,0.0,2,3,4结点坐标集结点未知量编号4.0,-3.0,5,6,08.0,0.0,7,0,81,2,6.0E+06,1.84E+052,4,6.0E+06,1.84E+05单元杆端结点编号及单元EA 、EI3,1,2.0E+06,0.03,2,2.0E+06,0.03,4,2.0E+06,0.0非结点荷载1. 0,1. 0,4. 0,-（ 2）输出结果NE=5 NJ=4 N=8 NW=7 NPJ=0NPF=1NODEXYXXYYZZ10.00000.000000124.00000.000023434.0000-3.000056048.00000.0000708ELEMENT

6、NODE-INODE-JEAEI1120.600000E+070.184000E+062240.600000E+070.184000E+063310.200000E+070.000000E+004320.200000E+070.000000E+005340.200000E+070.000000E+00ELEMENTINDAQ1.1.4.0000-20.0000NODEUVCETA10.000000E+000.000000E+00 0.312593E-032-0.202759E-04-0.253871E-03-0.144928E-033-0.202759E-04-0.185440E-030.00

7、0000E+004-0.405518E-040.000000E+00-0.227378E-04ELEMENTNQM1N1=30.4138Q1=37.1896M1=0.0000N2=-30.4138Q2=42.8104M2=11.24152N1=30.4138Q1=2.8104M1=-11.2415N2=-30.4138Q2=-2.8104M2=0.00003N1=-38.0173Q1=0.0000M1=0.0000N2=38.0173Q2=0.0000M2=0.00004N1=45.6207Q1=0.0000M1=0.0000N2=-45.6207Q2=0.0000M2=0.00005N1=-

8、38.0173Q1=0.0000M1=0.0000N2=38.0173Q2=0.0000M2=0.0000C主程序C (一)输入原始数据DIMENSION JE(2,100),JN(3,100),JC(6),EA(100),EI(100),X(100), $Y(100),PJ(2,50),PF(4,100)REAL*8 KE(6,6),KD(6,6),T(6,6),P(300),KB(200,20),F(6),FO(6), $D(6),BL,SI,CO,S,COPEN (5,FILE='RPF1.TXT')open (6,file='jieguo1.dat',s

9、tatus='new')READ (5,*) NE,NJ,N,NW,NPJ,NPFREAD (5,*) (X(J),Y(J),(JN(I,J),I=1,3),J=1,NJ)READ (5,*) (JE(I,J),I=1,2),EA(J),EI(J),J=1,NE)WRITE (6,20) (J,X(J),Y(J),(JN(I,J),I=1,3),J=1,NJ)WRITE (6,30) (J,(JE(I,J),I=1,2),EA(J),EI(J),J=1,NE)10 FORMAT(/6X,'NE=',I5,2X,'NJ=',I5,2X,'N

10、=',I5,2X,'NW=',I5,2X, $'NPJ=',I5,2X,'NPF='I5)20 FORMAT(/7X,'NODE',7X,'X',11X,'Y',12X,'XX',8X,'YY',8X,'ZZ'/ $(1X,I10,2F12.4,3I10)30 FORMAT(/4X,'ELEMENT',4X,'NODE-I',4X,'NODE-J',11X,'EA',13X,

11、9;EI'/ $(1X,3I10,2E15.6)40 FORMAT(/7X,'CODE',7X,'PX-PY-PM'/(1X,F10.0,F15.4)50 FORMAT(/4X,'ELEMENT',7X,'IND',10X,'A',14X,'Q',/ $(1X,2F10.0,2F15.4)C(二)形成总结点荷载向量DO 55 I=1,N55 P(I)=0.00DO 60 I=1,NPJL=PJ(1,I)60P(L)=PJ(2,I)65DO 70 I=1,NPFM=PF(1,I)CALL SC

12、L(M,NE,NJ,BL,SI,CO,JE,X,Y)CALL EFX(I,NPF,BL,PF,FO)CALL CTM(SI,CO,T)CALL EJC(M,NE,NJ,JE,JN,JC)DO 75 L=1,6S=0.00DO 80 K=1,680 S=S-T(K,L)*FO(K) F(L)=S75CONTINUEDO 85 J=1,6L=JC(J)P(L)=P(L)+F(J)85CONTINUE70CONTINUEC (三)形成整体刚度矩阵90DO 95 I=1,N DO 100 J=1,NW100 KB(I,J)=0.0095 CONTINUEDO 105 M=1,NECALL SCL(M,

13、NE,NJ,BL,SI,CO,JE,X,Y)CALL CTM(SI,CO,T)CALL ESM(M,NE,BL,EA,EI,KD)CALL EJC(M,NE,NJ,JE,JN,JC)DO 110 I=1,6DO 115 J=1,6S=0.00DO 120 L=1,6DO 125 K=1,6125 S=S+T(L,I)*KD(L,K)*T(K,J)120 CONTINUE KE(I,J)=S115 CONTINUE110 CONTINUE DO 130 L=1,6135 CONTINUE130 CONTINUE105 CONTINUEC(四)解线性方程组N1=N-1 DO 140 K=1,N1

14、IM=K+NW-1I1=K+1DO 145 I=I1,IML=I-K+1C=KB(K,L)/KB(K,1)JM=NW-L+1DO 150 J=1,JMJJ=J+I-K150 KB(I,J)=KB(I,J)-C*KB(K,JJ)145 P(I)=P(I)-C*P(K)140 CONTINUE P(N)=P(N)/KB(N,1)DO 155 K=1,N1160 P(I)=P(I)-KB(I,J)*P(L)155 P(I)=P(I)/KB(I,1)WRITE(6,165)165 FORMAT(/7X,'NODE',10X,'U',14X,'V',11X

15、,'CETA') DO 170 I=1,NJ175 CONTINUEWRITE(6,180) I,D(1),D(2),D(3)180 FORMAT(1X,I10,3E15.6)170 CONTINUEC(五)求单元杆端内力WRITE(6,200)200FORMAT(/4X,'ELEMENT',13X,'N',17X,'Q',17X,'M')DO 205 M=1,NECALL SCL(M,NE,NJ,BL,SI,CO,JE,X,Y)CALL ESM(M,NE,BL,EA,EI,KD)CALL CTM(SI,CO,T)

16、CALL EJC(M,NE,NJ,JE,JN,JC)DO 210 I=1,6L=JC(I)D(I)=0.00D(I)=P(L)210 CONTINUE DO 220 I=1,6240 F(I)=F(I)+KD(I,J)*T(J,K)*D(K)230 CONTINUE220DO 250 I=1,NPF L=PF(1,I)CALL EFX(I,NPF,BL,PF,FO)DO 260 J=1,6260F(J)=F(J)+FO(J)250CONTINUE270WRITE(6,280) M,(F(I),I=1,6)280 FORMAT(/1X,I10,3X,'N1=',F12.4,3X,

17、'Q1=',F12.4,3X,'M1=',F12.4 $/14X,'N2=',F12.4,3X,'Q2=',F12.4,3X,'M2=',F12.4)205 CONTINUE CLOSE(5)STOPENDC子程序C (六)形成单元定位向量SUBROUTINE EJC(M,NE,NJ,JE,JN,JC)DIMENSIONJE(2,NE),JN(3,NJ),JC(6)J1=JE(1,M)J2=JE(2,M)DO 10 I=1,3JC(I)=JN(I,J1)10JC(I+3)=JN(I,J2)RETURNENDC (七

18、)求单元常数SUBROUTINE SCL(M,NE,NJ,BL,SI,CO,JE,X,Y)DIMENSION JE(2,NE),X(NJ),Y(NJ)REAL*8 BL,SI,CO,DX,DYJ1=JE(1,M)J2=JE(2,M)DX=X(J2)-X(J1)DY=Y(J2)-Y(J1)BL=DSQRT(DX*DX+DY*DY)SI=DY/BLCO=DX/BLRETURNENDC (八)形成单元刚度矩阵SUBROUTINE ESM(M,NE,BL,EA,EI,KD)DIMENSION EA(NE),EI(NE)REAL*8 KD(6,6),BL,S,G ,G1,G2,G3G=EA(M)/BLG1=2.00*EI(M)/BLG2=3.00*G1/BLG3=2.00*G2/BLDO 10 I=1,6DO 10 J=1,610 KD(I,J)=0.00KD(1,1)=G30 KD(J,I)=KD(I,J)20 CONTINUERETURNENDC (九)形成单元坐标转换矩阵10 T(I,J)=0.00 T(1,1)=COT(1,2)=SIT(2,1)=-SIT(2,2)=COT(3,3)=1.00DO 20 I=1,3DO 20 J=1,320 T(I+3,J+3)=T(I,J)RETURNENDC (十)形成单元固