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

1、Fortran语言编写的有限元结构程序算例一计算简图及结果输出用平面刚架静力计算程序下图结构的内力。各杆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

2、+04结点荷载7. 0,-15.非结点荷载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

3、-PM7.-15.0000ELEMENTINDAQ1.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=

4、-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.0106 KN，平面刚架单元的抗拉刚度为已知：EA=4.0106KN，抗弯刚度为EI=1.84104KNm 2 。计算简图如下：（ 1）输入原始数据控制参数5,4,8,7,0（NE,NJ,N

5、,NW,NPJ,NPF）0.0,0.0,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 、EI 3,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

6、.00000.0000708ELEMENTNODE-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.202759

7、E-04-0.185440E-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

8、=0.0000M2=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,

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

10、,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)10FORMAT(/6X,NE=,I5,2X,NJ=,I5,2X,N=,I5,2X,NW=,I5,2X, $NPJ=,I5,2X,NPF=I5)20FORMAT(/7X,NODE,7X,X,11X,Y,12X,XX,8X,YY,8X,ZZ/ $(1X,I10,2F12.4

11、,3I10)30 FORMAT(/4X,ELEMENT,4X,NODE-I,4X,NODE-J,11X,EA,13X,EI/ $(1X,3I10,2E15.6)40FORMAT(/7X,CODE,7X,PX-PY-PM/(1X,F10.0,F15.4)50FORMAT(/4X,ELEMENT,7X,IND,10X,A,14X,Q,/ $(1X,2F10.0,2F15.4)C(二)形成总结点荷载向量DO 55 I=1,N55P(I)=0.00IF(NPJ.EQ.0) GO TO 65DO 60 I=1,NPJL=PJ(1,I)60 P(L)=PJ(2,I)65 IF(NPF.EQ.0) GO T

12、O 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,680S=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)85 CONTINUE70 CONTINUEC(三)形成整体刚度矩阵90DO 95 I=1,N DO 100 J=1,

13、NW100 KB(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)120CONTINUE KE(I,J)=S115CONTINUE110CONTINUE DO 130 L=1,6135CONTINUE130CONTI

14、NUE105CONTINUEC(四)解线性方程组N1=N-1 DO 140 K=1,N1 IM=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-K150KB(I,J)=KB(I,J)-C*KB(K,JJ)145P(I)=P(I)-C*P(K)140CONTINUE P(N)=P(N)/KB(N,1)DO 155 K=1,N1160P(I)=P(I)-KB(I,J)*P(L)155P(I)=P(I)/KB(I,1)WRITE(6,165)165FORMAT

15、(/7X,NODE,10X,U,14X,V,11X,CETA) DO 170 I=1,NJ175CONTINUEWRITE(6,180) I,D(1),D(2),D(3)180FORMAT(1X,I10,3E15.6)170CONTINUEC(五)求单元杆端内力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)D

16、O 210 I=1,6L=JC(I)D(I)=0.00IF(L.EQ.0) GO TO 210D(I)=P(L)210CONTINUE DO 220 I=1,6240F(I)=F(I)+KD(I,J)*T(J,K)*D(K)230CONTINUE220CONTINUE IF(NPF.EQ.0) GO TO 270DO 250 I=1,NPF L=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)2

17、80FORMAT(/1X,I10,3X,N1=,F12.4,3X,Q1=,F12.4,3X,M1=,F12.4 $/14X,N2=,F12.4,3X,Q2=,F12.4,3X,M2=,F12.4)205CONTINUE CLOSE(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

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

19、3.00*G1/BLG3=2.00*G2/BLDO 10 I=1,6DO 10 J=1,610KD(I,J)=0.00KD(1,1)=G30KD(J,I)=KD(I,J)20CONTINUERETURNENDC(九)形成单元坐标转换矩阵10T(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,320T(I+3,J+3)=T(I,J)RETURNENDC(十)形成单元固端力SUBROUTINE EFX(I,NPF,BL,PF,FO)DIMENSION PF(4,NPF)REAL*8 FO(

20、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,65FO(J)=0.00GO TO (10,20,30,40,50,60,70),IND10S=Q*A*0.50 FO(2)=-S*(2.00-2.00*G+C*G)20S=B/BL FO(2)=-Q*S*S*(1.00+2.00*C)30S=B/BL FO(2)=-6.00*Q*C*S/BLFO(5)=-FO(2)FO(3)=Q*S*(2.00-3.00*S)FO(6)=Q*C*(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 10070S=B/BL FO(2)=-Q*G*(3.00*S+C)100RETURN END