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

《Fortran语言编写的有限元结构程序》由会员分享，可在线阅读，更多相关《Fortran语言编写的有限元结构程序（11页珍藏版）》请在装配图网上搜索。

1、算例一计算简图及结果输出用平面刚架静力计算程序下图结构的内力。各杆 EA，EI 相同。已知：EA=4.0 x 106 KN,EI=1.6 x 104KN m2计算简图如下：U 土1）输入原始数据控制参数结点坐标集结点未知量编号单元杆端结点编号及单元 EA、EI结点荷载3,5,8,7,1,2(NE,NJ,N,NW,NPJ,NPF)0.0,0.0,0,0.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,-1 5.非结点荷载1

2、.0,2.0,2.0,2.0,1.0,4.0,-25.0V CETA2）输出结果NE= 3 NJ= 5 N= 8 NW= 7 NPJ= 1 NPF= 2NODEXYXXYYZZ10.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.0000ELE

3、MENTINDAQ1.2.2.0000-18.00002.1.4.0000-25.0000NODE U0.000000E+000.000000E+00 0.000000E+002 -0.221743E-02 -0.464619E-04 -0.139404E-023 -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.7119M1=-6.8477N2=

4、-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.0x 106KN，平面刚架单元的抗拉刚度为已知:EA=4.0 x106KN，抗弯刚度为 EI=1.84x 104KN m2。计算简图如下:1）输入原始数据控制参数5 , 4 , 8 , 7 (, 0NE,NJ,

5、N,NW,NPJ,NPF)结点坐标集结点未知量编号0.0,0.0,0,04.0,0.0,2,3,44.0,-3.0,5,6,08.0,0.0,7,0,8单元杆端结点编号及单元 EA、EI1,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=7NPJ= 0 NPF= 1NODEXYXXYYZZ10.00000.000000124.00000.000023434.0000-3.000056048.00000.0000708ELEME

6、NTNODE-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-03 -0.144928E-033 -0.202759E-04 -0.185

7、440E-03 0.000000E+004 -0.405518E-04 0.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.000

8、0M2=0.00005N1=-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,status=

9、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(I,J),I=1,3),J=1,N

10、J)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,I10,2F12.4,3I10)3

11、0 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)65IF(NPF.EQ.0) GO TO 90DO

12、 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,NDO 100 J=1,NW100KB(

13、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.EQ.0) GO TO 130DO

14、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)-C*P(K)140 CONTINU

15、EP(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(1),D(2),D(3)180 FOR

16、MAT(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)210 CONTINUEDO 220 I=1,6F(I)=0.00D

17、O 230 J=1,6DO 240 K=1,6240 F(I)=F(I)+KD(I,J)*T(J,K)*D(K)230 CONTINUE220 CONTINUEIF(NPF.EQ.0) GO TO 270 DO 250 I=1,NPF L=PF(1,I) IF(M.NE.L) GO TO 250 CALL 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=,F12.4,3X,M1=,

18、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) DIMENSION JE(2,NE

19、),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,

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

21、,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-ADO 5 J=1,65 FO(J)=0.00GO TO (10,20,30,40,50,60,70),I

22、ND10 S=Q*A*0.50FO(2)=-S*(2.00-2.00*G+C*G)FO(5)=-S*G*(2.00-C)S=S*A/6.00FO(3)=S*(6.00-8.00*C+3.00*G)FO(6)=-S*C*(4.00-3.00*C)GO TO 10020 S=B/BLFO(2)=-Q*S*S*(1.00+2.00*C)FO(5)=-Q*G*(1.00+2.00*S)FO(3)=Q*S*S*AFO(6)=-Q*B*GGO TO 10030 S=B/BLFO(2)=-6.00*Q*C*S/BL FO(5)=-FO(2) FO(3)=Q*S*(2.00-3.00*S) FO(6)=Q*C

23、*(2.00-3.00*C) GO TO 10040 S=Q*A*0.2500FO(2)=-S*(2.00-3.00*G+1.60*G*C)FO(5)=-S*G*(3.00-1.600*C)S=S*AFO(3)=S*(2.00-3.00*C+1.200*G)/1.500FO(6)=-S*C*(1.00-0.800*C)GO TO 10050 FO(1)=-Q*A*(1.00-0.500*C) FO(4)=-0.500*Q*C*AGO TO 10060 FO(1)=-Q*B/BLFO(4)=-Q*CGO TO 10070 S=B/BLFO(2)=-Q*G*(3.00*S+C)FO(5)=-FO(2)S=S*B/BLFO(3)=-Q*S*AFO(6)=Q*G*B100 RETURNEND