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

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

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,74.0,0.0,0,0,8单元杆端结点编号及单元EA、 EI 1,2,4.0E+06,1.6E+04 3,4,4.0E+06,1.6E+045,4,4.0E+06,1.6E+04结点荷载 7.0,-15.非结点荷载

2、1.0,2.0,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

3、0.160000E+057. -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.658499E-0

4、3ELEMENT 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、 : 6 EA=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+053,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

6、 N= 8 NW= 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.00

7、0000E+00ELEMENT 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

8、= 11.24152 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 (一输入原

9、始数据DIMENSION 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.TXTopen (6,file=jieguo1.dat,status=newREAD (5,* NE,NJ,N,NW,NPJ,NPFREAD (5,* (X(J,Y(J,(JN(I,J,I=1,3,J=1,NJREAD (5,* (JE(I,J,I=1,2,EA(J,

10、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,NODE,7X,X,11X,Y,12X,XX,8X,YY,8X,ZZ/ $(1X,I10,2F12.4,3I1030 FORMAT(/4X,ELEMENT,4X,NODE

11、-I,4X,NODE-J,11X,EA,13X,EI/$(1X,3I10,2E15.640 FORMAT(/7X,CODE,7X,PX-PY-PM/(1X,F10.0,F15.450 FORMAT(/4X,ELEMENT,7X,IND,10X,A,14X,Q,/$(1X,2F10.0,2F15.4C (二形成总结点荷载向量5560 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

12、L=1,6S=0.00DO 80 K=1,680 S=S-T(K,L*FO(KF(L=S75 CONTINUE DO 85 J=1,6L=JC(JP(L=P(L+F(J85 CONTINUE70 CONTINUEC (三形成整体刚度矩阵90 DO 95 I=1,N100 KB(I,J=0.0095 CONTINUE DO 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,KDCALL EJC(M,NE,NJ,JE,JN,JCDO 110 I=1,6DO 115 J=1,6S=0.

13、00DO 120 L=1,6DO 125 K=1,6125 S=S+T(L,I*KD(L,K*T(K,J120 CONTINUE KE(I,J=S115 CONTINUE110 CONTINUE DO 130 L=1,6 I=JC(L135 CONTINUE130 CONTINUE105 CONTINUEC (四解线性方程组150 KB(I,J=KB(I,J-C*KB(K,JJ145 P(I=P(I-C*P(K140 CONTINUE P(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(L

14、155 P(I=P(I/KB(I,1WRITE(6,165165 FORMAT(/7X,NODE,10X,U,14X,V,11X,CETA DO 170 I=1,NJ175 CONTINUE WRITE(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(S

15、I,CO,TCALL EJC(M,NE,NJ,JE,JN,JCDO 210 I=1,6L=JC(ID(I=0.00D(I=P(L210 CONTINUE DO 220 I=1,6 F(I=0.00DO 230 J=1,6 DO 240 K=1,6240 F(I=F(I+KD(I,J*T(J,K*D(K230 CONTINUE220L=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,

16、Q1=,F12.4,3X,M1=,F12.4 $/14X,N2=,F12.4,3X,Q2=,F12.4,3X,M2=,F12.4205 CONTINUE CLOSE(5STOP ENDC 子程序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,NJ,BL,SI,CO,JE,X,YDIMENSION JE(2,N

17、E,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 10 I=1,6DO 10 J=1,610 KD(I,J=0.00K

18、D(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=-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 CONTINUE RETURNENDC (九形成单元坐标转换矩阵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,

19、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=PF(4,IC=A/BLG=C*CB=BL-ADO 5 J=1,65 FO(J=0.00GO TO (10,20,30,40,50,60,70,IND10 S=Q*A*0.50FO(2=-S*(2.00-2.00*G+C*GFO(

20、5=-S*G*(2.00-CS=S*A/6.00FO(3=S*(6.00-8.00*C+3.00*GFO(6=-S*C*(4.00-3.00*CGO TO 10020 S=B/BLFO(2=-Q*S*S*(1.00+2.00*CFO(5=-Q*G*(1.00+2.00*SFO(3=Q*S*S*AFO(6=-Q*B*GGO TO 10030 S=B/BLFO(2=-6.00*Q*C*S/BLFO(5=-FO(2FO(3=Q*S*(2.00-3.00*SFO(6=Q*C*(2.00-3.00*CGO 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*CS=S*AFO(3=S*(2.00-3.00*C+1.200*G/1.500FO(6=-S*C*(1.00-0.800*C GO TO 100 50 FO(1=-Q*A*(1.00-0.500*C FO(4=-0.500*Q*C*A GO TO 100 60 FO(1=-Q*B/BL FO(4=-Q*C GO TO 100 70 S=B/BL FO(2=-Q*G*(3.00*S+C FO(5=-FO(2 S=S*B/BL FO(3=-Q*S*A FO(6=Q*G*B 100RETURN END