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

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

2、15.1. 0,2. 0 , 2. 0 ,-2. 0 , 1. 0,4. 0,-2 5. 07 NPJ= 1 NPF= 22）输出结果NODEXYXXYY10.00000.000000020.00004.000012330.00004.000012444.00004.000056754.00000.0000008ELEMENTNODE-INODE-JEAEI112 0.400000E+070.160000E+05234 0.400000E+070.160000E+05354 0.400000E+070.160000E+05CODEPX-PY-PM7.-15.0000ELEMENTINDAQ1.

3、2.2.0000-18.00002.1.4.0000-25.0000NODEUVCETA10.000000E+000.000000E+000.000000E+002-0.221743E-02-0.464619E-04-0.139404E-02NE= 3 NJ= 5 N= 8 NW=ZZ3 -0.221743E-02 -0.464619E-04 0.357876E-024 -0.222472E-02 -0.535381E-04 -0.298554E-025 0.000000E+00 0.000000E+00 0.658499E-03ELEMENTNQM1N1=46.4619Q1=10.7119M

4、1=-6.8477N2=-46.4619Q2=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.0 106KN，平面刚架单元的抗拉刚度为已知:EA=4.0 106KN，抗弯刚度为 El=1.84 1 04KN。计算简图如下：(1)输入原始数据控制参数5,4,8 , 7 ,(NE,N

5、J,N,NW,NPJ,NPF)结点坐标集结点未知量编号单元杆端结点编号及单元EA、EI0.0,0.0,0,04.0,0.0,2,3,44.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+053,1,2.0E+06,0.03,2,2.0E+06,0.03,4,2.0E+06,0.02）输出结果NE= 5 NJ= 4 N= 8 NW= 7 NPJ= 0 NPF= 1NODEXYXXYYZZ10.00000.000000124.00000.000023434.0000-3.000056048.00000.0000708EL

6、EMENTNODE-INODE-JEAEI112 0.600000E+070.184000E+06224 0.600000E+070.184000E+06331 0.200000E+070.000000E+00432 0.200000E+070.000000E+00534 0.200000E+070.000000E+00ELEMENTINDAQ1.1.4.0000-20.0000NODEUVCETA10.000000E+000.000000E+000.312593E-032-0.202759E-04-0.253871E-030.144928E-033-0.202759E-04-0.185440

7、E-030.000000E+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.0

8、0005N1=-38.0173Q1=0.0000M1=0.0000N2=38.0173Q2=0.0000M2=0.0000非结点荷载1. 0 , 1. 0,4. 0，-C 主程序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=j

9、ieguo1.dat,status=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)IF (NPJ.NE.0) READ (5,*) (PJ(I,J),I=1,2),J=1,NPJ)IF (NPF.NE.0) READ (5,*) (PF(I,J),I=1,4),J=1,NPF)WRITE (6,10) NE,NJ,N,NW,NPJ,NPFWRITE (6,20) (J,X(J),Y(J),(JN

10、(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)IF (NPJ.NE.0) WRITE (6,40) (PJ(I,J),I=1,2),J=1,NPJ)IF (NPF.NE.0) WRITE (6,50) (PF(I,J),I=1,4),J=1,NPF)10 FORMAT(/6X,NE=,I5,2X,NJ=,I5,2X,N=,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

11、,I10,2F12.4,3I10)30 FORMAT(/4X,ELEMENT,4X,NODE-I,4X,NODE-J,11X,EA,13X,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.00IF(NPJ.EQ.0) GO TO 65DO 60 I=1,NPJL=PJ(1,I)60 P(L)=PJ(2,I)65 IF(

12、NPF.EQ.0) GO TO 90DO 70 I=1,NPFM=PF(1,I)CALL SCL(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)=S75 CONTINUEDO 85 J=1,6L=JC(J)IF(L.EQ.0) GO TO 85P(L)=P(L)+F(J)85CONTINUE70CONTINUEC（三）形成整体刚度矩阵90DO 95 I=1,N

13、DO 100 J=1,NW100KB(I,J)=0.0095CONTINUEDO 105 M=1,NECALL SCL(M,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,6125S=S+T(L,I)*KD(L,K)*T(K,J)120CONTINUEKE(I,J)=S115CONTINUE110CONTINUEDO 130 L=1,6I=JC(L)IF(I

14、.EQ.0) GO TO 130DO 135 K=1,6J=JC(K)IF(J.EQ.0.OR.J.LT.I)GO TO 135JJ=J-I+1KB(I,JJ)=KB(I,JJ)+KE(L,K)135CONTINUE130CONTINUE105CONTINUEC（四）解线性方程组N1=N-1DO 140 K=1,N1IM=K+NW-1IF(N.LT.IM) IM=NI1=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

15、)-C*P(K)140 CONTINUEP(N)=P(N)/KB(N,1)DO 155 K=1,N1I=N-KJM=K+1IF(NW.LT.JM) JM=NWDO 160 J=2,JML=J+I-1160 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,CETA)DO 170 I=1,NJDO 175 J=1,3D(J)=0.00L=JN(J,I)IF(L.EQ.0) GO TO 175D(J)=P(L)175 CONTINUEWRITE(6,180) I,D(

16、1),D(2),D(3)180 FORMAT(1X,I10,3E15.6)170 CONTINUEC(五 ) 求单元杆端内力WRITE(6,200)200 FORMAT(/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)CALL EJC(M,NE,NJ,JE,JN,JC)DO 210 I=1,6L=JC(I)D(I)=0.00IF(L.EQ.0) GO TO 210D(I)=P(L)210CONTINUE24

17、0DO 220 I=1,6F(I)=0.00DO 230 J=1,6DO 240 K=1,6F(I)=F(I)+KD(I,J)*T(J,K)*D(K)230CONTINUE220CONTINUEIF(NPF.EQ.0) GO TO 270DO 250 I=1,NPFL=PF(1,I)IF(M.NE.L) GO TO 250CALL EFX(I,NPF,BL,PF,FO)DO 260 J=1,6260 F(J)=F(J)+FO(J)250 CONTINUE270 WRITE(6,280) M,(F(I),I=1,6)280 FORMAT(/1X,I10,3X,N1=,F12.4,3X,Q1=,F

18、12.4,3X,M1=,F12.4 $/14X,N2=,F12.4,3X,Q2=,F12.4,3X,M2=,F12.4)205 CONTINUECLOSE(5)STOPENDC 子程序C(六 ) 形成单元定位向量SUBROUTINE EJC(M,NE,NJ,JE,JN,JC)DIMENSION JE(2,NE),JN(3,NJ),JC(6)J1=JE(1,M)J2=JE(2,M)DO 10 I=1,3JC(I)=JN(I,J1)10 JC(I+3)=JN(I,J2)RETURNENDC(七 ) 求单元常数SUBROUTINE SCL(M,NE,NJ,BL,SI,CO,JE,X,Y)DIMENS

19、ION 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,G3 G=EA(M)/BLG1=2.00*EI(M)/BLG2=3.00*G1/BLG3=2.00*G2/BLDO 10 I

20、=1,6DO 10 J=1,610 KD(I,J)=0.00KD(1,1)=GKD(1,4)=-GKD(4,4)=GKD(2,2)=G3KD(5,5)=G3KD(2,5)=-G3KD(2,3)=-G2KD(2,6)=-G2KD(3,5)=G2KD(5,6)=G2KD(3,3)=2.00*G1KD(6,6)=2.00*G1KD(3,6)=G1DO 20 I=1,5I1=I+1DO 30 J=I1,630 KD(J,I)=KD(I,J)20 CONTINUERETURNENDC(九 ) 形成单元坐标转换矩阵SUBROUTINE CTM(SI,CO,T)REAL*8 T(6,6),SI,CODO 10 I=1,6DO 10 J=1,610 T(I,J)=0.00T(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(十 ) 形成单元固端力SUBROUTINE EFX(I,NPF,BL,PF,FO)DIMENSION PF(4,NPF)REAL*8 FO(6),A,B,C,G ,Q,S,BLIND=PF(2,I)A=PF(3,I)Q=PF(4,I)C=A/BLG=C*CB=BL