Fortran语言编写的有限元结构程序.docx
《Fortran语言编写的有限元结构程序.docx》由会员分享，可在线阅读，更多相关《Fortran语言编写的有限元结构程序.docx（14页珍藏版）》请在冰豆网上搜索。
Fortran语言编写的有限元结构程序
算例一计算简图及结果输出
用平面刚架静力计算程序下图结构的内力。
各杆EA，EI相同。
已知：
计算简图如下：
（1）输入原始数据
控制参数3,5,8,7,1,2（NE,NJ,N,NW,NPJ,NPF）
结点坐标集结点未知量编号
单元杆端结点编号及单元EA、EI
结点荷载
非结点荷载
（2）输出结果
NE=3NJ=5N=8NW=7NPJ=1NPF=2
NODEXYXXYYZZ
10.00000.0000000
20.00004.0000123
30.00004.0000124
44.00004.0000567
54.00000.0000008
ELEMENTNODE-INODE-JEAEI
1120.400000E+070.160000E+05
2340.400000E+070.160000E+05
3540.400000E+070.160000E+05
CODEPX-PY-PM
7.-15.0000
ELEMENTINDAQ
1.2.2.0000-18.0000
2.1.4.0000-25.0000
NODEUVCETA
10.000000E+000.000000E+000.000000E+00
2-0.221743E-02-0.464619E-04-0.139404E-02
3-0.221743E-02-0.464619E-040.357876E-02
4-0.222472E-02-0.535381E-04-0.298554E-02
50.000000E+000.000000E+000.658499E-03
ELEMENTNQM
1N1=46.4619Q1=10.7119M1=-6.8477
N2=-46.4619Q2=7.2881M2=0.0000
2N1=7.2881Q1=46.4619M1=0.0000
N2=-7.2881Q2=53.5381M2=14.1523
3N1=53.5381Q1=7.2881M1=0.0000
N2=-53.5381Q2=-7.2881M2=-29.1523
算例二计算简图及结果输出
用平面刚架静力计算程序下图结构的内力。
已知：
桁架单元的抗拉刚度为
，平面刚架单元的抗拉刚度为已知：
，抗弯刚度为
。
计算简图如下：
（1）输入原始数据
控制参数
（NE,NJ,N,NW,NPJ,NPF）
结点坐标集结点未知量编号
单元杆端结点编号及单元EA、EI
非结点荷载
（2）输出结果
NE=5NJ=4N=8NW=7NPJ=0NPF=1
NODEXYXXYYZZ
10.00000.0000001
24.00000.0000234
34.0000-3.0000560
48.00000.0000708
ELEMENTNODE-INODE-JEAEI
1120.600000E+070.184000E+06
2240.600000E+070.184000E+06
3310.200000E+070.000000E+00
4320.200000E+070.000000E+00
5340.200000E+070.000000E+00
ELEMENTINDAQ
1.1.4.0000-20.0000
NODEUVCETA
10.000000E+000.000000E+000.312593E-03
2-0.202759E-04-0.253871E-03-0.144928E-03
3-0.202759E-04-0.185440E-030.000000E+00
4-0.405518E-040.000000E+00-0.227378E-04
ELEMENTNQM
1N1=30.4138Q1=37.1896M1=0.0000
N2=-30.4138Q2=42.8104M2=11.2415
2N1=30.4138Q1=2.8104M1=-11.2415
N2=-30.4138Q2=-2.8104M2=0.0000
3N1=-38.0173Q1=0.0000M1=0.0000
N2=38.0173Q2=0.0000M2=0.0000
4N1=45.6207Q1=0.0000M1=0.0000
N2=-45.6207Q2=0.0000M2=0.0000
5N1=-38.0173Q1=0.0000M1=0.0000
N2=38.0173Q2=0.0000M2=0.0000
C主程序
C
（一）输入原始数据
DIMENSIONJE（2,100）,JN（3,100）,JC（6）,EA（100）,EI（100）,X（100）,
$Y（100）,PJ（2,50）,PF（4,100）
REAL*8KE（6,6）,KD（6,6）,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.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）
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,3I10））
30FORMAT（/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
（二）形成总结点荷载向量
DO55I=1,N
55P（I）=0.00
IF（NPJ.EQ.0）GOTO65
DO60I=1,NPJ
L=PJ（1,I）
60P（L）=PJ（2,I）
65IF（NPF.EQ.0）GOTO90
DO70I=1,NPF
M=PF（1,I）
CALLSCL（M,NE,NJ,BL,SI,CO,JE,X,Y）
CALLEFX（I,NPF,BL,PF,FO）
CALLCTM（SI,CO,T）
CALLEJC（M,NE,NJ,JE,JN,JC）
DO75L=1,6
S=0.00
DO80K=1,6
80S=S-T（K,L）*FO（K）
F（L）=S
75CONTINUE
DO85J=1,6
L=JC（J）
IF（L.EQ.0）GOTO85
P（L）=P（L）+F（J）
85CONTINUE
70CONTINUE
C（三）形成整体刚度矩阵
90DO95I=1,N
DO100J=1,NW
100KB（I,J）=0.00
95CONTINUE
DO105M=1,NE
CALLSCL（M,NE,NJ,BL,SI,CO,JE,X,Y）
CALLCTM（SI,CO,T）
CALLESM（M,NE,BL,EA,EI,KD）
CALLEJC（M,NE,NJ,JE,JN,JC）
DO110I=1,6
DO115J=1,6
S=0.00
DO120L=1,6
DO125K=1,6
125S=S+T（L,I）*KD（L,K）*T（K,J）
120CONTINUE
KE（I,J）=S
115CONTINUE
110CONTINUE
DO130L=1,6
I=JC（L）
IF（I.EQ.0）GOTO130
DO135K=1,6
J=JC（K）
IF（J.EQ.0.OR.J.LT.I）GOTO135
JJ=J-I+1
KB（I,JJ）=KB（I,JJ）+KE（L,K）
135CONTINUE
130CONTINUE
105CONTINUE
C（四）解线性方程组
N1=N-1
DO140K=1,N1
IM=K+NW-1
IF（N.LT.IM）IM=N
I1=K+1
DO145I=I1,IM
L=I-K+1
C=KB（K,L）/KB（K,1）
JM=NW-L+1
DO150J=1,JM
JJ=J+I-K
150KB（I,J）=KB（I,J）-C*KB（K,JJ）
145P（I）=P（I）-C*P（K）
140CONTINUE
P（N）=P（N）/KB（N,1）
DO155K=1,N1
I=N-K
JM=K+1
IF（NW.LT.JM）JM=NW
DO160J=2,JM
L=J+I-1
160P（I）=P（I）-KB（I,J）*P（L）
155P（I）=P（I）/KB（I,1）
WRITE（6,165）
165FORMAT（/7X,'NODE',10X,'U',14X,'V',11X,'CETA'）
DO170I=1,NJ
DO175J=1,3
D（J）=0.00
L=JN（J,I）
IF（L.EQ.0）GOTO175
D（J）=P（L）
175CONTINUE
WRITE（6,180）I,D
（1）,D
（2）,D（3）
180FORMAT（1X,I10,3E15.6）
170CONTINUE
C（五）求单元杆端内力
WRITE（6,200）
200FORMAT（/4X,'ELEMENT',13X,'N',17X,'Q',17X,'M'）
DO205M=1,NE
CALLSCL（M,NE,NJ,BL,SI,CO,JE,X,Y）
CALLESM（M,NE,BL,EA,EI,KD）
CALLCTM（SI,CO,T）
CALLEJC（M,NE,NJ,JE,JN,JC）
DO210I=1,6
L=JC（I）
D（I）=0.00
IF（L.EQ.0）GOTO210
D（I）=P（L）
210CONTINUE
DO220I=1,6
F（I）=0.00
DO230J=1,6
DO240K=1,6
240F（I）=F（I）+KD（I,J）*T（J,K）*D（K）
230CONTINUE
220CONTINUE
IF（NPF.EQ.0）GOTO270
DO250I=1,NPF
L=PF（1,I）
IF（M.NE.L）GOTO250
CALLEFX（I,NPF,BL,PF,FO）
DO260J=1,6
260F（J）=F（J）+FO（J）
250CONTINUE
270WRITE（6,280）M,（F（I）,I=1,6）
280FORMAT（/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）
STOP
END
C子程序
C（六）形成单元定位向量
SUBROUTINEEJC（M,NE,NJ,JE,JN,JC）
DIMENSIONJE（2,NE）,JN（3,NJ）,JC（6）
J1=JE（1,M）
J2=JE（2,M）
DO10I=1,3
JC（I）=JN（I,J1）
10JC（I+3）=JN（I,J2）
RETURN
END
C（七）求单元常数
SUBROUTINESCL（M,NE,NJ,BL,SI,CO,JE,X,Y）
DIMENSIONJE（2,NE）,X（NJ）,Y（NJ）
REAL*8BL,SI,CO,DX,DY
J1=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/BL
CO=DX/BL
RETURN
END
C（八）形成单元刚度矩阵
SUBROUTINEESM（M,NE,BL,EA,EI,KD）
DIMENSIONEA（NE）,EI（NE）
REAL*8KD（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
DO10I=1,6
DO10J=1,6
10KD（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）=G2
KD（3,3）=2.00*G1
KD（6,6）=2.00*G1
KD（3,6）=G1
DO20I=1,5
I1=I+1
DO30J=I1,6
30KD（J,I）=KD（I,J）
20CONTINUE
RETURN
END
C（九）形成单元坐标转换矩阵
SUBROUTINECTM（SI,CO,T）
REAL*8T（6,6）,SI,CO
DO10I=1,6
DO10J=1,6
10T（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
DO20I=1,3
DO20J=1,3
20T（I+3,J+3）=T（I,J）
RETURN
END
C（十）形成单元固端力
SUBROUTINEEFX（I,NPF,BL,PF,FO）
DIMENSIONPF（4,NPF）
REAL*8FO（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
DO5J=1,6
5FO（J）=0.00
GOTO（10,20,30,40,50,60,70）,IND
10S=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-3.00*C）
GOTO100
20S=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*G
GOTO100
30S=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）
GOTO100
40S=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）
GOTO100
50FO
（1）=-Q*A*（1.00-0.500*C）
FO（4）=-0.500*Q*C*A
GOTO100
60FO
（1）=-Q*B/BL
FO（4）=-Q*C
GOTO100
70S=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
升级会员 