Fortran 讲解平面有限元源程序Word下载.docx
《Fortran 讲解平面有限元源程序Word下载.docx》由会员分享，可在线阅读，更多相关《Fortran 讲解平面有限元源程序Word下载.docx（23页珍藏版）》请在冰豆网上搜索。
*格式说明，似乎1X前多个/
WRITE（6,50）EK
*输出EK
50
FORMAT（1X,6E11.4）
*同上
60
RETURN
*返回操作系统
END
*程序结束
EK（I1-1,J1-1）=A*（B（I）*B（J）+A2*B（I+3）*B（J+3））
假如I=2,J=3
I1=2*I=4,J1=2*J=6，这样上面的语句就相当于：
EK（3,5）=A*（B
（2）*B（3）+A2*B（5）*B（6））
呵呵，就这意思。
一维、二维数组你总该知道吧？
B
（2）就是元素在数组中第二个，EK（3，5）第三行第五列那个元素（行列式）。
有限元离不开行列式，甚至需要用到张量（N维数组）。
FORTRAN在计算上得天独厚。
平面弹性力学有限元源程序（FORTRAN）
1.$DEBUG
2.
PROGRAMPLANE
3.
IMPLICITREAL*8（A-H,O-Z）,INTEGER（I-N）
4.
ALLOCATABLE:
:
IJK（:
:
）,XY（:
）,BCA（:
）,SK（:
）,STR（:
）,MB（:
）,ZB（:
）,B（:
）
5.
DELD（:
）,TOD（:
）,DELST（:
）,TOST（:
）,DELSUP（:
）,TOTSUP（:
6.
DIMENSIONEK（6,6）
7.
CHARACTERPN*40,FN*12
8.
9.
WRITE（*,'
（A）'
）'
本程序为计算平面问题的有限元程序'
10.
特点：
（1）采用三结点三角形单元；
'
11.
（2）采用等带宽存贮技术；
12.
（3）采用高斯消元法解线性方程组。
13.
（/A）'
输入计算问题名（PN）:
14.
READ（*,'
）PN
15.
CALLFNAME（PN,'
.DAT'
FN）
16.
（2A）'
输入数据文件名为：
FN
17.
OPEN（5,FILE=FN,STATUS='
OLD'
18.
.OUT'
19.
（/2A）'
结果输出数据文件名为：
'
20.
OPEN（6,FILE=FN,STATUS='
UNKNOWN'
21.
.OU1'
22.
参数输出数据文件名为：
23.
OPEN（7,FILE=FN,STATUS='
24.
25.
READ（5,*）NG,NE,MC,NX,NB,EO,VO,DENSITY,T
26.
WRITE（6,120）NG,NE,MC,NX,NB
27.
WRITE（6,130）EO,VO,DENSITY,T
28.
READ（5,*）NWA,NWE,NWK,NWP
29.
NT=2*NG
30.
ALLOCATE（IJK（3,NE）,XY（2,NG）,BCA（7,NE）,STR（3,NE）,MB（2,NB）,ZB（NB）,B（NT））
31.
ALLOCATE（DELD（2,NG,NX）,TOD（2,NG）,DELST（3,NE,NX）,TOST（3,NE）,DELSUP（NB,NX）,TOTSUP（NB））
32.
CALLCLEAR（2,NG,TOD）
33.
CALLCLEAR（3,NE,TOST）
34.
CALLCLEAR1（NB,TOTSUP）
35.
IF（NG.EQ.0）THEN
36.
STOP000
37.
ENDIF
38.
CALLINPUT（NG,NE,NB,IJK,XY,MB,ZB）
39.
CALLCALND（NG,NE,IJK,ND）
40.
ALLOCATE（SK（NT,ND））
41.
IF（MC.EQ.0）THEN
42.
E=EO
43.
V=VO
44.
ELSE
45.
E=EO/（1-VO*VO）
46.
V=VO/（1-VO）
47.
48.
IP=0
49.
NX1=NX
50.
A1=E/（1-V*V）/4.0
51.
A2=0.5*（1-V）
52.
CALLABC（NG,NE,IJK,XY,BCA,NWA）
53.
CALLCLEAR（NT,ND,SK）
54.
DO100K=1,NE
55.
CALLKE（K,NE,T,V,A1,A2,IJK,BCA,EK,NWE）
56.
CALLSUMK（K,NE,NT,ND,IJK,EK,SK）
57.100
CONTINUE
58.
CALLCHECK（NT,ND,SK）
59.110
60.
IP=IP+1
61.
WRITE（6,'
（///1X,A,I4）'
）"
荷载工况="
IP
62.
READ（5,*）NF,NP,NM
63.
DOI=1,NT
64.
B（I）=0.0
65.
ENDDO
66.
IF（NF.GT.0）THEN
67.
CALLPF（NF,NT,B）
68.
69.
IF（NP.GT.0）THEN
70.
CALLPP（NP,NG,NT,T,XY,B）
71.
72.
IF（NM.GT.0）THEN
73.
CALLPG（NM,NE,NT,DENSITY,T,IJK,BCA,B）
74.
75.
DOI=1,NB
76.
I1=2*（MB（1,I）-1）+2-MB（2,I）
77.
DELSUP（I,IP）=-B（I1）
78.
79.
IF（NF.NE.0.OR.NP.NE.0.OR.NM.NE.0）THEN
80.
CALLDBC（NB,ND,NT,NX,NX1,MB,ZB,SK,B）
81.
CALLGAUSS（NT,ND,NX1,NX,SK,B,NWK,NWP）
82.
CALLSTRESS（NE,NT,V,A1,A2,IJK,BCA,B,STR）
83.
CALLTRES（IP,NG,NE,NX,NT,B,STR,DELD,TOD,DELST,TOST）
84.
CALLSUPT（IP,NG,NE,NB,NX,T,V,A1,A2,IJK,MB,BCA,DELD,DELSUP,TOTSUP）
85.
CALLOUTPUT（IP,NG,NE,NB,NX,MB,DELD,TOD,DELST,TOST,DELSUP,TOTSUP）
86.
87.
WRITE（*,'
（2x,A,I4,A）'
IP,"
没有荷载!
"
88.
89.
90.
NX1=NX1-1
91.
IF（NX1.GT.0）GOTO110
92.120
FORMAT（1X,'
结点总数='
I3,1X,'
单元总数='
问题类型='
&
93.
&
I1,1X,'
荷载工况数='
I2,1X,'
支承位移数='
I2）
94.130
FORMAT（/1X,'
弹性模量='
E10.3,1X,'
泊松比='
F5.3,1X,'
密度='
E10.3,&
95.
1X,'
厚度='
F6.3）
96.
END
97.
98.
99.
SUBROUTINEGAUSS（NT,ND,NX1,NX,SK,B,NWK,NWP）
100.
101.
DIMENSIONSK（NT,ND）,A（NT,NT）,B（NT）
102.
CALLCLEAR（NT,NT,A）
103.
104.
DOJ=1,ND
105.
IF（（I+J-1）.LE.NT）THEN
106.
A（I,I+J-1）=SK（I,J）
107.
ENDIF
108.
ENDDO
109.
110.
IF（NWK.EQ.1.AND.NX1.EQ.NX）THEN
111.
WRITE（7,'
结构总刚（列出右上三角元素）"
112.
113.
WRITE（7,100）（A（I,J）,J=1,NT）
114.
115.
116.
IF（NWP.EQ.1）THEN
117.
（1X,A,I3,A）'
第"
NX-NX1+1,"
荷载工况的荷载列阵"
118.
119.
（1x,E11.4）'
）B（I）
120.
121.
122.
DO50K=1,NT-1
123.
DO60I=K+1,NT
124.
C1=A（K,I）/A（K,K）
125.
DO70J=I,NT
126.
A（I,J）=A（I,J）-C1*A（K,J）
127.70
CONTINUE
128.
B（I）=B（I）-C1*B（K）
129.60
130.50
131.
B（NT）=B（NT）/A（NT,NT）
132.
DO80I=NT-1,1,-1
133.
DO90J=I+1,NT
134.
B（I）=B（I）-A（I,J）*B（J）
135.90
136.
B（I）=B（I）/A（I,I）
137.80
138.100
FORMAT（1x,4000E10.3）
139.
140.
141.
142.
SUBROUTINECHECK（NT,ND,SK）
143.
144.
DIMENSIONSK（NT,ND）
145.
M=0
146.
147.
IF（SK（I,1）.LE.0.0）THEN
148.
M=M+1
149.
150.
151.
IF（M.GT.0）THEN
152.
WRITE（*,*）"
总刚矩阵的对角元素为负值!
!
153.
STOP2222
154.
155.
156.
157.
SUBROUTINESTRESS（NE,NT,V,A1,A2,IJK,BCA,B,STR）
158.
159.
DIMENSIOND1（6）,D2（3）,B（NT）,S（3,6）,IJK（3,NE）,BCA（7,NE）,STR（3,NE）
160.
CALLCLEAR（3,NE,STR）
161.
DO10K=1,NE
162.
DOI=1,3
163.
II=IJK（I,K）
164.
D1（2*（I-1）+1）=B（2*（II-1）+1）
165.
D1（2*（I-1）+2）=B（2*（II-1）+2）
166.
167.
CALLCLEAR（3,6,S）
168.
C1=2*A1/BCA（7,K）
169.
DO20I=1,3
170.
S（1,2*（I-1）+1）=C1*BCA（I,K）
171.
S（1,2*（I-1）+2）=C1*V*BCA（I+3,K）
172.
S（2,2*（I-1）+1）=C1*V*BCA（I,K）
173.
S（2,2*（I-1）+2）=C1*BCA（I+3,K）
174.
S（3,2*（I-1）+1）=C1*A2*BCA（I+3,K）
175.
S（3,2*（I-1）+2）=C1*A2*BCA（I,K）
176.20
177.
CALLMATMUL（3,6,1,S,D1,D2）
178.
179.
STR（I,K）=D2（I）
180.
181.10
182.
183.
184.
SUBROUTINEOUTPUT（IP,NG,NE,NB,NX,MB,DELD,TOD,DELST,TOST,DELSUP,TOTSUP）
185.
186.
DIMENSIONMB（2,NB）,DELD（2,NG,NX）,TOD（2,NG）
187.
DIMENSIONDELST（3,NE,NX）,TOST（3,NE）,DELSUP（NB,NX）,TOTSUP（NB）
188.
CHARACTERVE*6
189.
WRITE（6,20）IP
190.
WRITE（6,30）
191.
DOI=1,NG
192.
WRITE（6,40）I,DELD（1,I,IP）,TOD（1,I）,DELD（2,I,IP）,T
升级会员 