Fortran語言編寫的有限元結(jié)構(gòu)程序

1、算例一計(jì)算簡圖及結(jié)果輸出用平面剛架靜力計(jì)算程序下圖結(jié)構(gòu)的內(nèi)力。各桿EA , EI相同已知：EA=4.0 106KN,EI=1.6 104KN m2計(jì)算簡圖如下：(1)輸入原始數(shù)據(jù)控制參數(shù)結(jié)點(diǎn)坐標(biāo)集結(jié)點(diǎn)未知量編號單元桿端結(jié)點(diǎn)編號及單元EA、EI結(jié)點(diǎn)荷載非結(jié)點(diǎn)荷載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）輸出結(jié)果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算例二計(jì)算簡圖及結(jié)果輸出用平面剛架靜力計(jì)算程序下圖結(jié)構(gòu)的內(nèi)力。已知：桁架單元的抗拉剛度為EA=2.0 106KN，平面剛架單元的抗拉剛度為已知:EA=4.0 106KN，抗彎剛度為 El=1.84 1 04KN。計(jì)算簡圖如下：(1)輸入原始數(shù)據(jù)控制參數(shù)5,4,8 , 7 ,(NE,N

5、J,N,NW,NPJ,NPF)結(jié)點(diǎn)坐標(biāo)集結(jié)點(diǎn)未知量編號單元桿端結(jié)點(diǎn)編號及單元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）輸出結(jié)果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非結(jié)點(diǎn)荷載1. 0 , 1. 0,4. 0，-C 主程序C(一 ) 輸入原始數(shù)據(jù)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(二 ) 形成總結(jié)點(diǎn)荷載向量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(五 ) 求單元桿端內(nèi)力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(七 ) 求單元常數(shù)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(九 ) 形成單元坐標(biāo)轉(zhuǎn)換矩陣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