Fortran语言编写的有限元结构程序.docx
《Fortran语言编写的有限元结构程序.docx》由会员分享，可在线阅读，更多相关《Fortran语言编写的有限元结构程序.docx（19页珍藏版）》请在冰豆网上搜索。
Fortran语言编写的有限元结构程序
算例一计算简图及结果输出
用平面刚架静力计算程序下图结构的内力。
各杆
EA，EI相同。
已知：
EA=4.0106KN,EI=1.6104KNm2
计算简图如下：
.
（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
1,2,4.0E+06,1.6E+04
单元杆端结点编号及单元
EA、EI
3,4,4.0E+06,1.6E+04
5,4,4.0E+06,1.6E+04
结点荷载
7.0,-15.
非结点荷载
1.0,2.0,2.0,-
2.0,1.0,4.0,-25.0
（2）输出结果
NE=
3NJ=
5N=
8
NW=
7NPJ=
1NPF=
2
NODE
X
Y
XX
YY
ZZ
1
0.0000
0.0000
0
0
0
2
0.0000
4.0000
1
2
3
3
0.0000
4.0000
1
2
4
4
4.0000
4.0000
5
6
7
5
4.0000
0.0000
0
0
8
ELEMENT
NODE-I
NODE-J
EA
EI
1
1
2
0.400000E+07
0.160000E+05
2
3
4
0.400000E+07
0.160000E+05
3
5
4
0.400000E+07
0.160000E+05
CODE
PX-PY-PM
7.
-15.0000
ELEMENT
IND
A
Q
1.
2.
2.0000
-18.0000
2.
1.
4.0000
-25.0000
NODE
U
V
CETA
1
0.000000E+00
0.000000E+00
0.000000E+00
2
-0.221743E-02
-0.464619E-04
-0.139404E-02
3
-0.221743E-02
-0.464619E-04
0.357876E-02
4
-0.222472E-02-0.535381E-04-0.298554E-02
5
0.000000E+00
0.000000E+00
0.658499E-03
ELEMENT
N
Q
M
1
N1=
46.4619
Q1=
10.7119
M1=
-6.8477
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
算例二计算简图及结果输出
用平面刚架静力计算程序下图结构的内力。
已知：
桁架单元的抗拉刚度为
EA=2.0
106KN
，平面刚架单元的抗拉刚度为已知：
EA=4.0
106KN
，抗弯刚度为
EI=1.84
104KN
m2。
计算简图如下：
（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
1,2,6.0E+06,1.84E+05
2,4,6.0E+06,1.84E+05
单元杆端结点编号及单元EA、EI3,1,2.0E+06,0.0
3,2,2.0E+06,0.0
3,4,2.0E+06,0.0
非结点荷载
1.0,1.0,4.0,-
（2）输出结果
NE=
5NJ=
4N=
8NW=
7NPJ=
0
NPF=
1
NODE
X
Y
XX
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
ELEMENT
NODE-I
NODE-J
EA
EI
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
ELEMENT
IND
A
Q
1.
1.
4.0000
-20.0000
NODE
U
V
CETA
1
0.000000E+00
0.000000E+000.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
ELEMENT
N
Q
M
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
N1=
-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.0000
N2=
38.0173
Q2=
0.0000
M2=
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,NDO100J=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）
120CONTINUEKE（I,J）=S
115CONTINUE
110CONTINUEDO130L=1,6
135CONTINUE
130CONTINUE
105CONTINUE
C（四）解线性方程组
N1=N-1DO140K=1,N1IM=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）
140CONTINUEP（N）=P（N）/KB（N,1）
DO155K=1,N1
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
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）
210CONTINUEDO220I=1,6
240F（I）=F（I）+KD（I,J）*T（J,K）*D（K）
230CONTINUE
220CONTINUEIF（NPF.EQ.0）GOTO270
DO250I=1,NPFL=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）
205CONTINUECLOSE（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
30KD（J,I）=KD（I,J）
20CONTINUE
RETURN
END
C（九）形成单元坐标转换矩阵
10T（I,J）=0.00T（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.50FO
（2）=-S*（2.00-2.00*G+C*G）
20S=B/BLFO
（2）=-Q*S*S*（1.00+2.00*C）
30S=B/BLFO
（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/BLFO
（2）=-Q*G*（3.00*S+C）
100RETURNEND
升级会员 