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

1、算例一计算简图及结果输出用平面刚架静力计算程序下图结构的内力。各杆EA,EI相同。已知:642EA=4.010KN,EI=1.610KN m计算简图如下: (1输入原始数据控制参数3,5,8,7,1,2(NE,NJ,N,NW,NPJ,NPF结点坐标集结点未知量编号0.0,0.0,0,0 0.0,4.0,1,2,3 0.0,4.0,1,2,4 4.0,4.0,5,6,7 4.0,0.0,0,0,8单元杆端结点编号及单元EA、EI 1,2,4.0E+06,1.6E+04 3,4,4.0E+06,1.6E+04 5,4,4.0E+06,1.6E+04结点荷载7.0,-15.非结点荷载1.0,2.0,

2、2.0,-2.0,1.0,4.0,-25.0(2输出结果NE= 3 NJ= 5 N= 8 NW= 7 NPJ= 1 NPF= 2 NODE X Y XX YY ZZ1 0.0000 0.0000 0 0 02 0.0000 4.0000 1 2 33 0.0000 4.0000 1 2 44 4.0000 4.000056 75 4.0000 0.0000 0 0 8ELEMENT NODE-I NODE-J EA EI1 12 0.400000E+07 0.160000E+052 3 4 0.400000E+07 0.160000E+053 54 0.400000E+07 0.160000E

3、+05CODE PX-PY-PM7. -15.0000ELEMENT IND A Q1. 2. 2.0000 -18.00002. 1. 4.0000 -25.0000NODE U V CETA1 0.000000E+00 0.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.65849

4、9E-03ELEMENT N Q M1 N1= 46.4619 Q1= 10.7119 M1= -6.8477N2= -46.4619 Q2= 7.2881 M2= 0.00002 N1= 7.2881 Q1= 46.4619 M1= 0.0000N2= -7.2881 Q2= 53.5381 M2= 14.15233 N1= 53.5381 Q1= 7.2881 M1= 0.0000N2= -53.5381 Q2= -7.2881 M2= -29.1523算例二计算简图及结果输出用平面刚架静力计算程序下图结构的内力。已知:桁架单元的抗拉刚度为6EA=2.010KN,平面刚架单元的抗拉刚度为已

5、知: 6EA=4.010KN,抗弯刚度为42EI=1.8410KN m。计算简图如下: (1输入原始数据控制参数5,4,8,7,0(NE,NJ,N,NW,NPJ,NPF结点坐标集结点未知量编号0.0,0.0,0,0 4.0,0.0,2,3,4 4.0,-3.0,5,6,0 8.0,0.0,7,0,8单元杆端结点编号及单元EA、EI 1,2,6.0E+06,1.84E+05 2,4,6.0E+06,1.84E+05 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=

6、7 NPJ= 0 NPF= 1 NODE X Y XX YY ZZ1 0.0000 0.0000 0 0 12 4.0000 0.0000 23 43 4.0000 -3.0000 5 6 04 8.0000 0.0000 7 0 8ELEMENT NODE-I NODE-J EA EI1 12 0.600000E+07 0.184000E+062 2 4 0.600000E+07 0.184000E+063 3 1 0.200000E+07 0.000000E+004 3 2 0.200000E+07 0.000000E+005 3 4 0.200000E+07 0.000000E+00EL

7、EMENT IND A Q1. 1. 4.0000 -20.0000NODE U V CETA1 0.000000E+00 0.000000E+00 0.312593E-032 -0.202759E-04 -0.253871E-03 -0.144928E-033 -0.202759E-04 -0.185440E-03 0.000000E+004 -0.405518E-04 0.000000E+00 -0.227378E-04ELEMENT N Q M1 N1= 30.4138 Q1= 37.1896 M1= 0.0000N2= -30.4138 Q2= 42.8104 M2= 11.24152

8、 N1= 30.4138 Q1= 2.8104 M1= -11.2415N2= -30.4138 Q2= -2.8104 M2= 0.00003 N1= -38.0173 Q1= 0.0000 M1= 0.0000N2= 38.0173 Q2= 0.0000 M2= 0.00004 N1= 45.6207 Q1= 0.0000 M1= 0.0000N2= -45.6207 Q2= 0.0000 M2= 0.00005 N1= -38.0173 Q1= 0.0000 M1= 0.0000N2= 38.0173 Q2= 0.0000 M2= 0.0000C 主程序C (一输入原始数据DIMENSI

9、ON JE(2,100,JN(3,100,JC(6,EA(100,EI(100,X(100,$Y(100,PJ(2,50,PF(4,100REAL*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='new'READ (5,* NE,NJ,N,NW,NPJ,NPFREAD (5,* (X(J,Y(J,(JN(I,J,I=1,3,J=1,NJREAD (5,*

10、 (JE(I,J,I=1,2,EA(J,EI(J,J=1,NEWRITE (6,10 NE,NJ,N,NW,NPJ,NPFWRITE (6,20 (J,X(J,Y(J,(JN(I,J,I=1,3,J=1,NJWRITE (6,30 (J,(JE(I,J,I=1,2,EA(J,EI(J,J=1,NE10 FORMAT(/6X,'NE=',I5,2X,'NJ=',I5,2X,'N=',I5,2X,'NW=',I5,2X,$'NPJ=',I5,2X,'NPF='I520 FORMAT(/7X,'N

11、ODE',7X,'X',11X,'Y',12X,'XX',8X,'YY',8X,'ZZ'/$(1X,I10,2F12.4,3I1030 FORMAT(/4X,'ELEMENT',4X,'NODE-I',4X,'NODE-J',11X,'EA',13X,'EI'/$(1X,3I10,2E15.640 FORMAT(/7X,'CODE',7X,'PX-PY-PM'/(1X,F10.0,F15.450

12、FORMAT(/4X,'ELEMENT',7X,'IND',10X,'A',14X,'Q',/$(1X,2F10.0,2F15.4C (二形成总结点荷载向量DO 55 I=1,N55 P(I=0.00DO 60 I=1,NPJL=PJ(1,I60 P(L=PJ(2,IDO 70 I=1,NPFM=PF(1,ICALL SCL(M,NE,NJ,BL,SI,CO,JE,X,YCALL EFX(I,NPF,BL,PF,FOCALL CTM(SI,CO,TCALL EJC(M,NE,NJ,JE,JN,JCDO 75 L=1,6S=0.00

13、DO 80 K=1,680 S=S-T(K,L*FO(KF(L=S75 CONTINUEDO 85 J=1,6L=JC(JP(L=P(L+F(J85 CONTINUE70 CONTINUEC (三形成整体刚度矩阵90 DO 95 I=1,NDO 100 J=1,NW100 KB(I,J=0.0095 CONTINUEDO 105 M=1,NECALL SCL(M,NE,NJ,BL,SI,CO,JE,X,Y CALL CTM(SI,CO,TCALL ESM(M,NE,BL,EA,EI,KDCALL EJC(M,NE,NJ,JE,JN,JCDO 110 I=1,6DO 115 J=1,6S=0.0

14、0DO 120 L=1,6DO 125 K=1,6125 S=S+T(L,I*KD(L,K*T(K,J120 CONTINUEKE(I,J=S115 CONTINUE110 CONTINUEDO 130 L=1,6I=JC(LDO 135 K=1,6J=JC(KJJ=J-I+1KB(I,JJ=KB(I,JJ+KE(L,K135 CONTINUE130 CONTINUE105 CONTINUEC (四解线性方程组N1=N-1DO 140 K=1,N1IM=K+NW-1I1=K+1DO 145 I=I1,IML=I-K+1C=KB(K,L/KB(K,1JM=NW-L+1DO 150 J=1,JMJ

15、J=J+I-K150 KB(I,J=KB(I,J-C*KB(K,JJ145 P(I=P(I-C*P(K140 CONTINUEP(N=P(N/KB(N,1DO 155 K=1,N1I=N-KJM=K+1DO 160 J=2,JML=J+I-1160 P(I=P(I-KB(I,J*P(L155 P(I=P(I/KB(I,1WRITE(6,165165 FORMAT(/7X,'NODE',10X,'U',14X,'V',11X,'CETA' DO 170 I=1,NJDO 175 J=1,3D(J=0.00L=JN(J,ID(J=P(

16、L175 CONTINUEWRITE(6,180 I,D(1,D(2,D(3180 FORMAT(1X,I10,3E15.6170 CONTINUEC (五求单元杆端内力WRITE(6,200200 FORMAT(/4X,'ELEMENT',13X,'N',17X,'Q',17X,'M' DO 205 M=1,NECALL SCL(M,NE,NJ,BL,SI,CO,JE,X,YCALL ESM(M,NE,BL,EA,EI,KDCALL CTM(SI,CO,TCALL EJC(M,NE,NJ,JE,JN,JCDO 210 I=1,6

17、L=JC(ID(I=0.00D(I=P(L210 CONTINUEDO 220 I=1,6F(I=0.00DO 230 J=1,6DO 240 K=1,6240 F(I=F(I+KD(I,J*T(J,K*D(K230 CONTINUE220 CONTINUEDO 250 I=1,NPFL=PF(1,ICALL EFX(I,NPF,BL,PF,FODO 260 J=1,6260 F(J=F(J+FO(J250 CONTINUE270 WRITE(6,280 M,(F(I,I=1,6280 FORMAT(/1X,I10,3X,'N1=',F12.4,3X,'Q1='

18、,F12.4,3X,'M1=',F12.4 $/14X,'N2=',F12.4,3X,'Q2=',F12.4,3X,'M2=',F12.4205 CONTINUECLOSE(5STOPENDC 子程序C (六形成单元定位向量SUBROUTINE EJC(M,NE,NJ,JE,JN,JCDIMENSION JE(2,NE,JN(3,NJ,JC(6J1=JE(1,MJ2=JE(2,MDO 10 I=1,3JC(I=JN(I,J110 JC(I+3=JN(I,J2RETURNENDC (七求单元常数SUBROUTINE SCL(M,NE

19、,NJ,BL,SI,CO,JE,X,YDIMENSION JE(2,NE,X(NJ,Y(NJREAL*8 BL,SI,CO,DX,DYJ1=JE(1,MJ2=JE(2,MDX=X(J2-X(J1DY=Y(J2-Y(J1BL=DSQRT(DX*DX+DY*DYSI=DY/BLCO=DX/BLRETURNENDC (八形成单元刚度矩阵SUBROUTINE ESM(M,NE,BL,EA,EI,KD DIMENSION EA(NE,EI(NEREAL*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

20、 10 I=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,J20 CONTINUERETURNENDC (九形成单元坐标转换矩阵SUBROUTINE CTM(SI,CO,TREAL*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,JRETURNENDC (十形成单元固端力SUBROUTINE EFX(I,NPF,BL,PF,FO DIMENSION PF(4,NPFREAL*8 FO(6,A,B,C,G,Q,S,BLIND=PF(2,IA=PF(3,IQ=