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

1、 N2= -46.4619 Q2= 7.2881 M2= 0.0000 2 N1= 7.2881 Q1= 46.4619 M1= 0.0000 N2= -7.2881 Q2= 53.5381 M2= 14.1523 3 N1= 53.5381 Q1= 7.2881 M1= 0.0000 N2= -53.5381 Q2= -7.2881 M2= -29.1523算例二计算简图及结果输出桁架单元的抗拉刚度为，平面刚架单元的抗拉刚度为已知：，抗弯刚度为。控制参数 （NE,NJ,N,NW,NPJ,NPF） NE= 5 NJ= 4 N= 8 NW= 7 NPJ= 0 NPF= 1 NODE X Y XX

2、 YY ZZ 1 0.0000 0.0000 0 0 1 2 4.0000 0.0000 2 3 4 3 4.0000 -3.0000 5 6 0 4 8.0000 0.0000 7 0 8 1 1 2 0.600000E+07 0.184000E+06 2 2 4 0.600000E+07 0.184000E+06 3 3 1 0.200000E+07 0.000000E+00 4 3 2 0.200000E+07 0.000000E+00 5 3 4 0.200000E+07 0.000000E+00 1. 1. 4.0000 -20.0000 1 0.000000E+00 0.00000

3、0E+00 0.312593E-03 2 -0.202759E-04 -0.253871E-03 -0.144928E-03 3 -0.202759E-04 -0.185440E-03 0.000000E+00 4 -0.405518E-04 0.000000E+00 -0.227378E-04 1 N1= 30.4138 Q1= 37.1896 M1= 0.0000 N2= -30.4138 Q2= 42.8104 M2= 11.2415 2 N1= 30.4138 Q1= 2.8104 M1= -11.2415 N2= -30.4138 Q2= -2.8104 M2= 0.0000 3 N

4、1= -38.0173 Q1= 0.0000 M1= 0.0000 N2= 38.0173 Q2= 0.0000 M2= 0.0000 4 N1= 45.6207 Q1= 0.0000 M1= 0.0000 N2= -45.6207 Q2= 0.0000 M2= 0.0000 5 N1= -38.0173 Q1= 0.0000 M1= 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

5、）,T（6,6）,P（300）,KB（200,20）,F（6）,FO（6）, $D（6）,BL,SI,CO,S,C OPEN （5,FILE=RPF1.TXT） open （6,file=jieguo1.dat,status=new READ （5,*） NE,NJ,N,NW,NPJ,NPF READ （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

6、.0） READ （5,*） （PF（I,J）,I=1,4）,J=1,NPF） WRITE （6,10） NE,NJ,N,NW,NPJ,NPF WRITE （6,20） （J,X（J）,Y（J）,（JN（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

7、,NJ=N=NW=,I5,2X, $NPJ=NPF=I5）20 FORMAT（/7X,NODE,7X,X,11X,Y,12X,XX,8X,YYZZ/ $（1X,I10,2F12.4,3I10）30 FORMAT（/4X,ELEMENT,4X,NODE-INODE-JEA,13X,EI $（1X,3I10,2E15.6）40 FORMAT（/7X,CODEPX-PY-PM/（1X,F10.0,F15.4）50 FORMAT（/4X,IND,10X,A,14X,Q,/ $（1X,2F10.0,2F15.4）C （二）形成总结点荷载向量 DO 55 I=1,N55 P（I）=0.00 IF（NPJ.

8、EQ.0） GO TO 65 DO 60 I=1,NPJ L=PJ（1,I）60 P（L）=PJ（2,I）65 IF（NPF.EQ.0） GO TO 90 DO 70 I=1,NPF M=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,6 S=0.00 DO 80 K=1,680 S=S-T（K,L）*FO（K） F（L）=S75 CONTINUE DO 85 J=1,6 L=JC（J） IF（L.E

9、Q.0） GO TO 85 P（L）=P（L）+F（J）85 CONTINUE70 CONTINUEC （三）形成整体刚度矩阵90 DO 95 I=1,N DO 100 J=1,NW100 KB（I,J）=0.0095 CONTINUE DO 105 M=1,NE CALL ESM（M,NE,BL,EA,EI,KD） DO 110 I=1,6 DO 115 J=1,6 DO 120 L=1,6 DO 125 K=1,6125 S=S+T（L,I）*KD（L,K）*T（K,J）120 CONTINUE KE（I,J）=S115 CONTINUE110 CONTINUE DO 130 L=1,6

10、I=JC（L） IF（I.EQ.0） GO TO 130 DO 135 K=1,6 J=JC（K） IF（J.EQ.0.OR.J.LT.I）GO TO 135 JJ=J-I+1 KB（I,JJ）=KB（I,JJ）+KE（L,K）135 CONTINUE130 CONTINUE105 CONTINUEC （四）解线性方程组 N1=N-1 DO 140 K=1,N1 IM=K+NW-1 IF（N.LT.IM） IM=N I1=K+1 DO 145 I=I1,IM L=I-K+1 C=KB（K,L）/KB（K,1） JM=NW-L+1 DO 150 J=1,JM JJ=J+I-K150 KB（I,J

11、）=KB（I,J）-C*KB（K,JJ）145 P（I）=P（I）-C*P（K）140 CONTINUE P（N）=P（N）/KB（N,1） DO 155 K=1,N1 I=N-K JM=K+1 IF（NW.LT.JM） JM=NW DO 160 J=2,JM L=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,UVCETA DO 170 I=1,NJ DO 175 J=1,3 D（J）=0.00 L=JN（J,I） IF（L.EQ.0） GO TO 175 D（J）=P（L）175

12、 CONTINUE WRITE（6,180） I,D（1）,D（2）,D（3）180 FORMAT（1X,I10,3E15.6）170 CONTINUEC （五）求单元杆端内力 WRITE（6,200）200 FORMAT（/4X,N,17X,M） DO 205 M=1,NE DO 210 I=1,6 L=JC（I） D（I）=0.00 IF（L.EQ.0） GO TO 210 D（I）=P（L）210 CONTINUE DO 220 I=1,6 F（I）=0.00 DO 230 J=1,6 DO 240 K=1,6240 F（I）=F（I）+KD（I,J）*T（J,K）*D（K）230 CO

13、NTINUE220 CONTINUE IF（NPF.EQ.0） GO TO 270 DO 250 I=1,NPF L=PF（1,I） IF（M.NE.L） GO TO 250 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=M1=,F12.4 $/14X,N2=Q2=M2=,F12.4）205 CONTINUE CLOSE（5） STOP ENDC 子程序C （六）形成单元定位向量 SUBROUTINE EJC（M,NE

14、,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,3 JC（I）=JN（I,J1）10 JC（I+3）=JN（I,J2） RETURNC （七）求单元常数 SUBROUTINE SCL（M,NE,NJ,BL,SI,CO,JE,X,Y） DIMENSION JE（2,NE）,X（NJ）,Y（NJ） REAL*8 BL,SI,CO,DX,DY DX=X（J2）-X（J1） DY=Y（J2）-Y（J1） BL=DSQRT（DX*DX+DY*DY） SI=DY/BL CO=DX/BLC （八）

15、形成单元刚度矩阵 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）/BL G1=2.00*EI（M）/BL G2=3.00*G1/BL G3=2.00*G2/BL DO 10 I=1,6 DO 10 J=1,610 KD（I,J）=0.00 KD（1,1）=G KD（1,4）=-G KD（4,4）=G KD（2,2）=G3 KD（5,5）=G3 KD（2,5）=-G3 KD（2,3）=-G2 KD（2,6）=-G2 KD（3,5）=G2 KD（5,6）=

16、G2 KD（3,3）=2.00*G1 KD（6,6）=2.00*G1 KD（3,6）=G1 DO 20 I=1,5 I1=I+1 DO 30 J=I1,630 KD（J,I）=KD（I,J）20 CONTINUE RETURN C （九）形成单元坐标转换矩阵 SUBROUTINE CTM（SI,CO,T） REAL*8 T（6,6）,SI,CO10 T（I,J）=0.00 T（1,1）=CO T（1,2）=SI T（2,1）=-SI T（2,2）=CO T（3,3）=1.00 DO 20 I=1,3 DO 20 J=1,320 T（I+3,J+3）=T（I,J）C （十）形成单元固端力 SUB

17、ROUTINE EFX（I,NPF,BL,PF,FO） DIMENSION PF（4,NPF） REAL*8 FO（6）,A,B,C,G,Q,S,BL IND=PF（2,I） A=PF（3,I） Q=PF（4,I） C=A/BL G=C*C B=BL-A DO 5 J=1,65 FO（J）=0.00 GO TO （10,20,30,40,50,60,70）,IND10 S=Q*A*0.50 FO（2）=-S*（2.00-2.00*G+C*G） FO（5）=-S*G*（2.00-C） S=S*A/6.00 FO（3）=S*（6.00-8.00*C+3.00*G） FO（6）=-S*C*（4.00

18、-3.00*C） GO TO 10020 S=B/BL FO（2）=-Q*S*S*（1.00+2.00*C） FO（5）=-Q*G*（1.00+2.00*S） FO（3）=Q*S*S*A FO（6）=-Q*B*G30 S=B/BL FO（2）=-6.00*Q*C*S/BL FO（5）=-FO（2） FO（3）=Q*S*（2.00-3.00*S） FO（6）=Q*C*（2.00-3.00*C）40 S=Q*A*0.2500 FO（2）=-S*（2.00-3.00*G+1.60*G*C） FO（5）=-S*G*（3.00-1.600*C） S=S*A FO（3）=S*（2.00-3.00*C+1.200*G）/1.500 FO（6）=-S*C*（1.00-0.800*C）50 FO（1）=-Q*A*（1.00-0.500*C） FO（4）=-0.500*Q*C*A60 FO（1）=-Q*B/BL FO（4）=-Q*C70 S=B/BL FO（2）=-Q*G*（3.00*S+C） S=S*B/BL FO（3）=-Q*S*A FO（6）=Q*G*B100 RETURN