1 C 这是一个采用平面四边形八节点等参单元的平面有限元分析程序。
2 PROGRAM PLANEFEM
3 IMPLICIT REAL*8(A-H,O-Z)
4 CHARACTER*80 LINECHAR,NEWLINECHAR
5 COMMON A(30000),L(4000)
6 COMMON /SOL/NPOIN,NELEM,NTYPE,NMATS
7 OPEN(5,FILE='FEMDATA',STATUS='OLD')
8 OPEN(6,FILE='FEMOUT',STATUS='UNKNOWN')
9 C 以下进行的是从数据文件FEMDATA中读入数据文件主标题信息。
10 READ(5,5000) LINECHAR
11 LOCATECHAR=INDEX(LINECHAR,'输入')
12 IF(LOCATECHAR.NE.0) THEN
13 LINECHAR(LOCATECHAR:LOCATECHAR+3)='输出'
14 ENDIF
15 WRITE(6,5100) LINECHAR
16 5000 FORMAT(A)
17 5100 FORMAT(80('*')/A/80('*')/)
18 C 以下进行的是先读入一行字符，然后从字符中读入网格单元总数NELEM。
19 READ(5,5000) LINECHAR
20 WRITE(6,5000) LINECHAR
21 LOCATECHAR=INDEX(LINECHAR,'=')
22 LOCATECHAR=LOCATECHAR+1
23 NEWLINECHAR=LINECHAR(LOCATECHAR:80)
24 READ(NEWLINECHAR,5200) NELEM
25 5200 FORMAT(I5)
26 C 以下进行的是先读入一行字符，然后从字符中读入单元节点总数NPOIN。
27 READ(5,5000) LINECHAR
28 WRITE(6,5000) LINECHAR
29 LOCATECHAR=INDEX(LINECHAR,'=')
30 LOCATECHAR=LOCATECHAR+1
31 NEWLINECHAR=LINECHAR(LOCATECHAR:80)
32 READ(NEWLINECHAR,5200) NPOIN
33 C 以下进行的是先读入一行字符，然后从字符中读入问题类型编号NTYPE。
34 READ(5,5000) LINECHAR
35 WRITE(6,5000) LINECHAR
36 LOCATECHAR=INDEX(LINECHAR,'=')
37 LOCATECHAR=LOCATECHAR+1
38 NEWLINECHAR=LINECHAR(LOCATECHAR:80)
39 READ(NEWLINECHAR,5200) NTYPE
40 C 以下进行的是先读入一行字符，然后从字符中读入材料类型总数NMATS。
41 READ(5,5000) LINECHAR
42 WRITE(6,5000) LINECHAR
43 LOCATECHAR=INDEX(LINECHAR,'=')
44 LOCATECHAR=LOCATECHAR+1
45 NEWLINECHAR=LINECHAR(LOCATECHAR:80)
46 READ(NEWLINECHAR,5200) NMATS
47 C 以下进行的是根据以输入的控制参数确定虚拟数组在实数组中的起始位置。
48 N1=1
49 N2=N1+NPOIN*2
50 N3=N2+NMATS*3
51 N4=N3+NPOIN*2
52 NN2=N1+NELEM*8
53 NN3=NN2+NPOIN
54 NN4=NN3+NELEM
55 NN5=NN4+NELEM
56 C 以下进行的是调用子程序进行有限元分析计算。
57 CALL INPUT(A(N1),A(N2),L(N1),L(NN2),L(NN3))
58 CALL STIF(A(N1),A(N2),L(N1),L(NN3))
59 CALL LOAD(A(N1),L(N1),L(NN4))
60 CALL ASEM(A(N3),A(N4),L(N1),L(NN2),L(NN4),L(NN5))
61 CALL SOLVE(A(N3),A(N4),L(NN5))
62 CALL STRE(A(N1),A(N2),A(N3),L(N1),L(NN3))
63 WRITE(*,5300)
64 5300 FORMAT(////////////////)
65 STOP ' *** 程序正常运行结束 ***'
66 END
67 C*****************************************************************************
68 C 这是一个进行数据输入并进行数据检验的子程序。 *
69 C*****************************************************************************
70 SUBROUTINE INPUT(COORD,PROPS,LNODS,IFPRE,MATNO)
71 IMPLICIT REAL*8(A-H,O-Z)
72 CHARACTER*80 LINECHAR,NEWLINECHAR,CHARCOMP
73 LOGICAL NOTREADCHAR
74 DIMENSION COORD(2,1),PROPS(3,1),LNODS(8,1),IFPRE(1),MATNO(1)
75 DIMENSION NEROR(6)
76 COMMON/SOL/NPOIN,NELEM,NTYPE,NMATS
77 CHARCOMP='aa(这个字符常量在ASCII序列中介于数字和字符之间)'
78 NNODE=8
79 C 以下进行的是对材料类型编号赋初值。
80 I=1
81 DO WHILE (I.LE.NELEM)
82 MATNO(I)=1
83 I=I+1
84 ENDDO
85 C 以下进行的是输入材料类型编号信息标题。
86 READ(5,5000) LINECHAR
87 WRITE(6,5000) LINECHAR
88 5000 FORMAT(A)
89 C 以下进行的是输入材料类型编号。
90 DO WHILE (NMATS.GT.1)
91 NUMBERREAD=1
92 DO WHILE (NUMBERREAD.LE.NMATS)
93 READ(5,*) I,MATNO(I)
94 NUMBERREAD=NUMBERREAD+1
95 ENDDO
96 ENDDO
97 C 以下进行的是输出材料类型编号。
98 I=1
99 DO WHILE (I.LE.NELEM)
100 WRITE(6,6000) I,MATNO(I)
101 I=I+1
102 ENDDO
103 6000 FORMAT(2X,'单元编号为：',I4,10X,'材料类型编号为：'I6)
104 C 以下进行的是输入单元节点总体编号标题。
105 READ(5,5000) LINECHAR
106 LOCATECHAR=INDEX(LINECHAR,'输入')
107 IF(LOCATECHAR.NE.0) THEN
108 LINECHAR(LOCATECHAR:LOCATECHAR+3)='输出'
109 ENDIF
110 WRITE(6,5000) LINECHAR
111 C 以下进行的是输入单元节点总体编号。
112 NOTREADCHAR=.TRUE.
113 DO WHILE (NOTREADCHAR)
114 READ(5,5000) LINECHAR
115 IF (LINECHAR(1:2).LT.CHARCOMP(1:2)) THEN
116 READ(LINECHAR,6100) IELEM,(LNODS(I,IELEM),I=1,NNODE)
117 J=1
118 DO WHILE (J.LE.10)
119 KELEM=IELEM+J
120 IF (KELEM.GT.NELEM) GOTO 100
121 K=1
122 DO WHILE (K.LE.8)
123 LNODS(K,KELEM)=LNODS(K,IELEM)+J
124 K=K+1
125 ENDDO
126 100 J=J+1
127 ENDDO
128 ELSE
129 NOTREADCHAR=.FALSE.
130 END IF
131 ENDDO
132 6100 FORMAT(9I5)
133 C 以下进行的是输出单元节点总体编号。
134 IELEM=1
135 DO WHILE (IELEM.LE.NELEM)
136 WRITE(6,6200) IELEM,(LNODS(I,IELEM),I=1,NNODE)
137 IELEM=IELEM+1
138 ENDDO
139 6200 FORMAT(2X,'单元编号为：',I4,' 节点总体编号为：',8I5)
140 C 以下进行的是输出单元节点的直角坐标信息标题。
141 LOCATECHAR=INDEX(LINECHAR,'输入')
142 IF(LOCATECHAR.NE.0) THEN
143 LINECHAR(LOCATECHAR:LOCATECHAR+3)='输出'
144 ENDIF
145 WRITE(6,5000) LINECHAR
146 C 以下进行的是输入单元节点的直角坐标信息。
147 NOTREADCHAR=.TRUE.
148 DO WHILE (NOTREADCHAR)
149 READ(5,5000) LINECHAR
150 IF (LINECHAR(1:2).LT.CHARCOMP(1:2)) THEN
151 READ(LINECHAR,6300) IPOIN,(COORD(J,IPOIN),J=1,2)
152 ELSE
153 NOTREADCHAR=.FALSE.
154 END IF
155 ENDDO
156 6300 FORMAT(I5,2F15.5)
157 C 以下进行的是输出单元节点的极坐标信息标题。
158 LOCATECHAR=INDEX(LINECHAR,'输入')
159 IF(LOCATECHAR.NE.0) THEN
160 LINECHAR(LOCATECHAR:LOCATECHAR+3)='输出'
161 ENDIF
162 WRITE(6,5000) LINECHAR
163 C 以下进行的是输入单元节点的极坐标信息。
164 PI=DATAN(1.D0)*4.D0
165 NOTREADCHAR=.TRUE.
166 DO WHILE (NOTREADCHAR)
167 READ(5,5000) LINECHAR
168 IF (LINECHAR(1:2).LT.CHARCOMP(1:2)) THEN
169 READ(LINECHAR,6300) IPOIN,R,ANGLE
170 ANGLE=ANGLE/180.D0*PI
171 COORD(1,IPOIN)=R*DCOS(ANGLE)
172 COORD(2,IPOIN)=R*DSIN(ANGLE)
173 NUMBERREAD=NUMBERREAD+1
174 ELSE
175 NOTREADCHAR=.FALSE.
176 END IF
177 ENDDO
178 C 以下进行的是根据输入信息计算出单元节点的坐标。
179 I=1
180 DO WHILE (I.LE.NELEM)
181 INODE=1
182 DO WHILE (INODE.LE.NNODE)
183 NODST =LNODS(INODE,I)
184 IGASH=INODE+2
185 IF(IGASH.GT.NNODE) IGASH=1
186 NDOFN=LNODS(IGASH,I)
187 MIDPT=INODE+1
188 NODMD=LNODS(MIDPT,I)
189 J=1
190 DO WHILE (J.LE.2)
191 IF(DABS(COORD(J,NODMD)).LT.1.E-6) THEN
192 COORD(J,NODMD)=(COORD(J,NDOFN)+COORD(J,NODST))/2.0
193 ENDIF
194 J=J+1
195 ENDDO
196 INODE=INODE+2
197 ENDDO
198 I=I+1
199 ENDDO
200 C 以下进行的是输出单元节点的坐标信息。
201 WRITE(6,6400)
202 I=1
203 DO WHILE (I.LE.NPOIN)
204 WRITE(6,6500) I,(COORD(J,I),J=1,2)
205 I=I+1
206 ENDDO
207 6400 FORMAT(4X,'节点编号',10X,'Ｘ坐标',16X,'Ｙ坐标')
208 6500 FORMAT(2X,I10,F16.5,6X,F16.5)
209 C 以下进行的是输出单元节点约束代码信息标题。
210 LOCATECHAR=INDEX(LINECHAR,'输入')
211 IF(LOCATECHAR.NE.0) THEN
212 LINECHAR(LOCATECHAR:LOCATECHAR+3)='输出'
213 ENDIF
214 WRITE(6,5000) LINECHAR
215 C 以下进行的是输入单元节点约束代码信息。
216 NOTREADCHAR=.TRUE.
217 DO WHILE (NOTREADCHAR)
218 READ(5,5000) LINECHAR
219 IF (LINECHAR(1:2).LT.CHARCOMP(1:2)) THEN
220 READ(LINECHAR,6600) I,IFPRE(I)
221 ELSE
222 NOTREADCHAR=.FALSE.
223 END IF
224 ENDDO
225 6600 FORMAT(2I5)
226 C 以下进行的是输出单元节点约束代码信息。
227 LINECHAR='节点约束代码信息为一个两位数的阿拉伯数字，其中，'
228 WRITE(6,5000) LINECHAR
229 LINECHAR='十位数表示Ｘ方向的约束状况，'
230 NEWLINECHAR='个位数表示Ｙ方向的约束状况，'
231 LOCATECHAR=LEN_TRIM(LINECHAR)
232 LOCATECHAR=LOCATECHAR+1
233 LINECHAR(LOCATECHAR:)=NEWLINECHAR(1:LOCATECHAR)
234 WRITE(6,5000) LINECHAR
235 LINECHAR='具体的，0表示自由，1表示固定，2表示给定非零位移。'
236 WRITE(6,5000) LINECHAR
237 I=1
238 DO WHILE (I.LE.NPOIN)
239 WRITE(6,6700) I,IFPRE(I)
240 I=I+1
241 ENDDO
242 6700 FORMAT(2X,'节点编号为：',I4,13X,'约束代码为：',I4)
243 C 以下进行的是输入材料的特性参数信息标题。
244 LINECHAR='以下输出的是材料特性参数（格式为：'
245 NEWLINECHAR='材料类型编号NMATS,弹性模量Ｅ,泊松比μ,厚度t）'
246 LOCATECHAR=LEN_TRIM(LINECHAR)
247 LOCATECHAR=LOCATECHAR+1
248 LINECHAR(LOCATECHAR:)=NEWLINECHAR(1:46)
249 WRITE(6,5000) LINECHAR
250 C 以下进行的是输入材料的特性参数。
251 NUMBERREAD=1
252 DO WHILE (NUMBERREAD.LE.NMATS)
253 READ(5,*) I,(PROPS(J,I),J=1,3)
254 WRITE(6,6800) I
255 WRITE(6,6900) (PROPS(J,I),J=1,3)
256 NUMBERREAD=NUMBERREAD+1
257 ENDDO
258 6800 FORMAT(2X,'材料类型编号为：',I2,' 的特性参数为：')
259 6900 FORMAT(2X,1PE13.3,1PE13.3,1PE13.3)
260 C 以下进行的是对输入和计算出的数据进行检验。
261 IERROR=1
262 DO WHILE (IERROR.LE.6)
263 NEROR(IERROR)=0
264 IERROR=IERROR+1
265 ENDDO
266 C 以下进行的是对输入和计算出的单元节点坐标进行检验。
267 IPOIN=1
268 DO WHILE (IPOIN.LE.NPOIN)
269 KPOIN=IPOIN-1
270 JPOIN=1
271 DO WHILE (JPOIN.LE.KPOIN)
272 J=1
273 DO WHILE (J.LE.2)
274 IF(COORD(J,IPOIN).NE.COORD(J,JPOIN)) GOTO 200
275 J=J+1
276 ENDDO
277 NEROR(1)=NEROR(1)+1
278 WRITE(6,7000) IPOIN,JPOIN
279 JPOIN=JPOIN+1
280 ENDDO
281 200 CONTINUE
282 IPOIN=IPOIN+1
283 ENDDO
284 7000 FORMAT(2X,'坐标错误！节点编号为',I4,' 的坐标等于',I4,' 的坐标')
285 C 以下进行的是对输入和计算出的材料类型编号进行检验。
286 I=1
287 DO WHILE (I.LE.NELEM)
288 IF(MATNO(I).LE.0.OR.MATNO(I).GT.NMATS) THEN
289 NEROR(2)=NEROR(2)+1
290 END IF
291 I=I+1
292 ENDDO
293 IF(NEROR(2).NE.0) THEN
294 WRITE(6,7100) NEROR(2)
295 END IF
296 7100 FORMAT(2X,' 材料类型编号错误,错误个数为',I2,' 个')
297 C 以下进行的是对输入和计算出的单元节点总体编号进行检验。
298 I=1
299 DO WHILE (I.LE.NELEM)
300 J=1
301 DO WHILE (J.LE.8)
302 IF(LNODS(J,I).LE.0.OR.LNODS(J,I).GT.NPOIN) THEN
303 NEROR(3)=NEROR(3)+1
304 END IF
305 J=J+1
306 ENDDO
307 I=I+1
308 ENDDO
309 IF(NEROR(3).NE.0) THEN
310 WRITE(6,7200) NEROR(3)
311 END IF
312 7200 FORMAT(2X,'单元节点总体编号错误,错误个数为',I5,' 个')
313 C 以下进行的是对输入和计算出的单元节点编号进行检验。
314 I=1
315 DO WHILE (I.LE.NPOIN)
316 KSTAR=0
317 IELEM=1
318 DO WHILE (IELEM.LE.NELEM)
319 KZERO=0
320 J=1
321 DO WHILE(J.LE.8)
322 IF(LNODS(J,IELEM).EQ.I) THEN
323 KZERO=KZERO+1
324 KSTAR=IELEM
325 END IF
326 J=J+1
327 ENDDO
328 IF(KZERO.GT.1) THEN
329 NEROR(4)=NEROR(4)+1
330 WRITE(6,7300) I,IELEM
331 END IF
332 IELEM=IELEM+1
333 ENDDO
334 IF(KSTAR.EQ.0)THEN
335 NEROR(5)=NEROR(5)+1
336 WRITE(6,7400) I
337 END IF
338 I=I+1
339 ENDDO
340 7300 FORMAT(2X,'节点编号',I5,' 在单元',I5,' 中重复。')
341 7400 FORMAT(2X,'总体节点编号',I5,' 在单元中从未出现')
342 C 以下进行的是对输入的单元节点约束代码进行检验。
343 I=1
344 DO WHILE (I.LE.NPOIN)
345 IF(IFPRE(I).LT.0) THEN
346 NEROR(6)=NEROR(6)+1
347 END IF
348 IF(IFPRE(I).GT.22) THEN
349 NEROR(6)=NEROR(6)+1
350 END IF
351 I=I+1
352 ENDDO
353 IF(NEROR(6).NE.0) THEN
354 WRITE(6,7500) NEROR(6)
355 END IF
356 7500 FORMAT(2X,'总体节点编号为',I5,' 的节点约束代码出现错误。')
357 C 以下进行的是判断出错类型。
358 KEROR=0
359 IERROR=1
360 DO WHILE (IERROR.LE.6)
361 IF(NEROR(IERROR).NE.0) THEN
362 WRITE(6,7600) IERROR
363 KEROR=KEROR+1
364 END IF
365 IERROR=IERROR+1
366 ENDDO
367 7600 FORMAT(2X,'出现第 ',I1,' 类错误。')
368 C 以下进行的是判断是否出错,如出错,则中断程序运行,否则返回调用主程序。
369 7700 FORMAT(2X,'*****在输入数据文件中发现错误。程序运行中断*****')
370 IF(KEROR.NE.0) THEN
371 WRITE(6,7700)
372 STOP '*** 程序运行被中断,这是由于输入数据出现错误引起的。***'
373 END IF
374 RETURN
375 END
376 C*****************************************************************************
377 C 这是一个形成单元刚度矩阵的子程序。 *
378 C*****************************************************************************
379 SUBROUTINE STIF(COORD,PROPS,LNODS,MATNO)
380 IMPLICIT REAL*8(A-H,O-Z)
381 DIMENSION COORD(2,1),PROPS(3,1),LNODS(8,1),MATNO(1)
382 DIMENSION POSGP(3),WEIGP(3),ESTIF(136),ELCOD(2,8),DBMAT(3)
383 COMMON/SOL/NPOIN,NELEM,NTYPE,NMATS
384 COMMON/GCOD/GPCOD(2,9)
385 COMMON/BMDMP/BMATX(3,16),DMATX(3,3)
386 OPEN(7,FILE='ESTIF',FORM='UNFORMATTED')
387 C 以下进行的是高斯积分点位置赋值。
388 POSGP(1)=-DSQRT(0.6D0)
389 POSGP(2)=0.0
390 POSGP(3)=DSQRT(0.6D0)
391 C 以下进行的是加权系数赋值。
392 WEIGP(1)=5.0/9.0
393 WEIGP(2)=8.0/9.0
394 WEIGP(3)=5.0/9.0
395 C 以下进行的是应力分量个数赋值。
396 NSTRE=3
397 C 以下进行的是形成弹性矩阵[D]。
398 DO 90 IELEM=1,NELEM
399 WRITE(6,615) IELEM
400 615 FORMAT(2X,'网格单元编号为：',I5)
401 LPROP=MATNO(IELEM)
402 DO 10 I=1,8
403 LNODE=LNODS(I,IELEM)
404 DO 10 J=1,2
405 10 ELCOD(J,I)=COORD(J,LNODE)
406 YOUNG=PROPS(1,LPROP)
407 POISS=PROPS(2,LPROP)
408 THICK=PROPS(3,LPROP)
409 CALL MODPS(YOUNG,POISS)
410 C 以下进行的是存放[K]的一维数组赋初值零。
411 KEVAB=0
412 DO 20 I=1,16
413 DO 20 J=1,I
414 KEVAB=KEVAB+1
415 ESTIF(KEVAB)=0.0
416 20 CONTINUE
417 C 以下进行的是计算积分点处形函数及其导数之值以及对整体坐标的偏导数之值。
418 KGASP=0
419 DO 80 IGAUS=1,3
420 DO 80 JGAUS=1,3
421 KGASP=KGASP+1
422 EXISP=POSGP(IGAUS)
423 ETASP=POSGP(JGAUS)
424 WRITE(6,705) KGASP
425 705 FORMAT(6X,'高斯积分点编号为：',I5)
426 CALL SFR2(EXISP,ETASP)
427 CALL JACOB2(IELEM,DJACB,KGASP,ELCOD)
428 DVOLU=DJACB*WEIGP(IGAUS)*WEIGP(JGAUS)*THICK
429 CALL BMATPS
430 C 以下进行的是弹性矩阵[D]赋初值零。
431 KEVAB=0
432 DO 70 IEVAB=1,16
433 DO 40 ISTRE=1,NSTRE
434 DBMAT(ISTRE)=0.0
435 DO 30 JSTRE=1,NSTRE
436 DBMAT(ISTRE)=DBMAT(ISTRE)+BMATX(JSTRE,IEVAB)*DMATX(JSTRE,ISTRE)
437 30 CONTINUE
438 40 CONTINUE
439 DO 60 JEVAB=1,IEVAB
440 KEVAB=KEVAB+1
441 BTDBM=0.0
442 DO 50 ISTRE=1,NSTRE
443 BTDBM=BTDBM+DBMAT(ISTRE)*BMATX(ISTRE,JEVAB)
444 50 CONTINUE
445 ESTIF(KEVAB)=ESTIF(KEVAB)+BTDBM*DVOLU
446 60 CONTINUE
447 70 CONTINUE
448 80 CONTINUE
449 WRITE(7) ESTIF
450 90 CONTINUE
451 REWIND 7
452 RETURN
453 END
454 C*****************************************************************************
455 C 这是一个形成弹性矩阵的子程序。 *
456 C*****************************************************************************
457 SUBROUTINE MODPS(Y,P)
458 IMPLICIT REAL*8(A-H,O-Z)
459 COMMON/SOL/NPOIN,NELEM,NTYPE,NMATS
460 COMMON/BMDMP/BMATX(3,16),DMATX(3,3)
461 DO 100 I=1,3
462 DO 100 J=1,3
463 DMATX(I,J)=0.0
464 100 CONTINUE
465 C 如果NTYPE=1,表示为平面应力问题。
466 IF(NTYPE.EQ.1) THEN
467 CONST=Y/(1.0-P*P)
468 DMATX(1,1)=CONST
469 DMATX(2,2)=CONST
470 DMATX(1,2)=CONST*P
471 DMATX(2,1)=CONST*P
472 DMATX(3,3)=CONST*(1.0-P)/2.0
473 ENDIF
474 C 如果NTYPE=2,表示为平面应变问题。
475 IF(NTYPE.EQ.2) THEN
476 CONST=Y*(1.0-P)/((1.0+P)*(1.0-2.0*P))
477 DMATX(1,1)=CONST
478 DMATX(2,2)=CONST
479 DMATX(1,2)=CONST*P/(1.0-P)
480 DMATX(2,1)=CONST*P/(1.0-P)
481 DMATX(3,3)=Y/(2.0*(1.0+P))
482 ENDIF
483 RETURN
484 END
485 C*****************************************************************************
486 C 这是一个计算形函数当前积分点及形函数对局部坐标的导数值的子程序。 *
487 C*****************************************************************************
488 SUBROUTINE SFR2(S,T)
489 IMPLICIT REAL*8(A-H,O-Z)
490 COMMON/SD/SHAPE(8),DERIV(2,8)
491 C 以下进行的是为简化表达式定义一些变量。
492 S2=S*2.0
493 T2=T*2.0
494 SS=S*S
495 TT=T*T
496 ST=S*T
497 STT=S*T*T
498 SST=S*S*T
499 ST2=S*T*2.0
500 C 以下进行的是给形函数赋值。
501 SHAPE(1)=(-1.0+ST+SS+TT-SST-STT)/4.0
502 SHAPE(2)=(1.0-T-SS+SST)/2.0
503 SHAPE(3)=(-1.0-ST+SS+TT-SST+STT)/4.0
504 SHAPE(4)=(1.0+S-TT-STT)/2.0
505 SHAPE(5)=(-1.0+ST+SS+TT+SST+STT)/4.0
506 SHAPE(6)=(1.0+T-SS-SST)/2.0
507 SHAPE(7)=(-1.0-ST+SS+TT+SST-STT)/4.0
508 SHAPE(8)=(1.0-S-TT+STT)/2.0
509 C 以下进行的是给形函数对局部坐标的导数赋值。
510 DERIV(1,1)=(T+S2-ST2-TT)/4.0
511 DERIV(1,2)=-S+ST
512 DERIV(1,3)=(-T+S2-ST2+TT)/4.0
513 DERIV(1,4)=(1.0-TT)/2.0
514 DERIV(1,5)=(T+S2+ST2+TT)/4.0
515 DERIV(1,6)=-S-ST
516 DERIV(1,7)=(-T+S2+ST2-TT)/4.0
517 DERIV(1,8)=(-1.0+TT)/2.0
518 DERIV(2,1)=(S+T2-SS-ST2)/4.0
519 DERIV(2,2)=(-1.0+SS)/2.0
520 DERIV(2,3)=(-S+T2-SS+ST2)/4.0
521 DERIV(2,4)=-T-ST
522 DERIV(2,5)=(S+T2+SS+ST2)/4.0
523 DERIV(2,6)=(1.0-SS)/2.0
524 DERIV(2,7)=(-S+T2+SS-ST2)/4.0
525 DERIV(2,8)=-T+ST
526 RETURN
527 END
528 C*****************************************************************************
529 C 这是一个形成雅可比矩阵的子程序。 *
530 C*****************************************************************************
531 SUBROUTINE JACOB2(IELEM,DJACB,KGASP,ELCOD)
532 IMPLICIT REAL*8(A-H,O-Z)
533 DIMENSION XJACM(2,2),XJACI(2,2),ELCOD(2,8)
534 COMMON/SD/SHAPE(8),DERIV(2,8)
535 COMMON/CAR/CARTD(2,8)
536 COMMON/GCOD/GPCOD(2,9)
537 C 以下进行的是计算当前积分点的整体坐标。
538 DO 100 I=1,2
539 GPCOD(I,KGASP)=0.0
540 DO 100 J=1,8
541 GPCOD(I,KGASP)=GPCOD(I,KGASP)+ELCOD(I,J)*SHAPE(J)
542 100 CONTINUE
543 C 以下进行的是形成雅可比矩阵[J]的各元素。
544 DO 200 I=1,2
545 DO 200 J=1,2
546 XJACM(I,J)=0.0
547 DO 200 K=1,8
548 XJACM(I,J)=XJACM(I,J)+DERIV(I,K)*ELCOD(J,K)
549 200 CONTINUE
550 C 以下进行的是计算雅可比行列式┃Ｊ┃的值。
551 DJACB=XJACM(1,1)*XJACM(2,2)-XJACM(1,2)*XJACM(2,1)
552 WRITE(6,6000) DJACB
553 6000 FORMAT(14X,'雅可比行列式┃Ｊ┃的值为：',F12.5)
554 IF(DJACB.LT.1.E-6)THEN
555 WRITE(6,6100) IELEM
556 6100 FORMAT('单元编号为：',I4,' 的雅可比行列式的值小于或等于零！')
557 STOP '****** 程序运行被中断于子程序JACOB2 ******'
558 END IF
559 C 以下进行的是形成雅可比矩阵[J]的逆矩阵的各元素。
560 XJACI(1,1)=XJACM(2,2)/DJACB
561 XJACI(2,2)=XJACM(1,1)/DJACB
562 XJACI(1,2)=-XJACM(1,2)/DJACB
563 XJACI(2,1)=-XJACM(2,1)/DJACB
564 C 以下进行的是计算形函数对整体坐标的导数。
565 DO 300 I=1,2
566 DO 300 K=1,8
567 CARTD(I,K)=0.0
568 DO 300 J=1,2
569 CARTD(I,K)=CARTD(I,K)+XJACI(I,J)*DERIV(J,K)
570 300 CONTINUE
571 RETURN
572 END
573 C*****************************************************************************
574 C 这是一个形成应变矩阵的子程序。 *
575 C*****************************************************************************
576 SUBROUTINE BMATPS
577 IMPLICIT REAL*8(A-H,O-Z)
578 COMMON/CAR/CARTD(2,8)
579 COMMON/SD/SHAPE(8),DERIV(2,8)
580 COMMON/BMDMP/BMATX(3,16),DMATX(3,3)
581 COMMON/SOL/NPOIN,NELEM,NTYPE,NMATS
582 COMMON/GCOD/GPCOD(2,9)
583 C 以下进行的是应变矩阵[B]赋初值零。
584 DO 100 I=1,16
585 DO 100 J=1,3
586 100 BMATX(J,I)=0.0
587 C 以下进行的是形成应变矩阵[B]。
588 NGASH=0
589 DO 200 I=1,8
590 MGASH=NGASH+1
591 NGASH=MGASH+1
592 BMATX(1,MGASH)=CARTD(1,I)
593 BMATX(1,NGASH)=0.0
594 BMATX(2,MGASH)=0.0
595 BMATX(2,NGASH)=CARTD(2,I)
596 BMATX(3,MGASH)=CARTD(2,I)
597 BMATX(3,NGASH)=CARTD(1,I)
598 200 CONTINUE
599 RETURN
600 END
601 C*****************************************************************************
602 C 这是一个计算等效节点载荷的子程序。 *
603 C*****************************************************************************
604 SUBROUTINE LOAD(COORD,LNODS,NLOAD)
605 IMPLICIT REAL*8(A-H,O-Z)
606 CHARACTER*80 LINECHAR,NEWLINECHAR
607 DIMENSION COORD(2,1),LNODS(8,1),NLOAD(1)
608 DIMENSION ELCOD(2,3),ELOAD(16),POSGP(3),WEIGP(3),NOPRS(3)
609 DIMENSION PRESS(3,2),PGASH(2),DGASH(2)
610 COMMON/SOL/NPOIN,NELEM,NTYPE,NMATS
611 COMMON/SD/SHAPE(8),DERIV(2,8)
612 OPEN(8,FILE='ELOAD',FORM='UNFORMATTED')
613 C 以下进行的是高斯积分点位置赋值。
614 POSGP(1)=-DSQRT(0.6D0)
615 POSGP(2)=0.0
616 POSGP(3)=DSQRT(0.6D0)
617 C 以下进行的是加权系数赋值。
618 WEIGP(1)=5.0/9.0
619 WEIGP(2)=8.0/9.0
620 WEIGP(3)=5.0/9.0
621 C 以下进行的是输入两行单元受载信息标题。
622 READ(5,5000) LINECHAR
623 LOCATECHAR=INDEX(LINECHAR,'输入')
624 IF(LOCATECHAR.NE.0) THEN
625 LINECHAR(LOCATECHAR:LOCATECHAR+3)='输出'
626 ENDIF
627 WRITE(6,5000) LINECHAR
628 READ(5,5000) LINECHAR
629 LOCATECHAR=INDEX(LINECHAR,'输入')
630 IF(LOCATECHAR.NE.0) THEN
631 LINECHAR(LOCATECHAR:LOCATECHAR+3)='输出'
632 ENDIF
633 WRITE(6,5000) LINECHAR
634 LOCATECHAR=INDEX(LINECHAR,'=')
635 LOCATECHAR=LOCATECHAR+1
636 NEWLINECHAR=LINECHAR(LOCATECHAR:80)
637 READ(NEWLINECHAR,5200) NUMBER
638 5000 FORMAT(A)
639 5200 FORMAT(I5)
640 C 以下进行的是单元受载边数赋初值零。
641 DO 10 I=1,NELEM
642 10 NLOAD(I)=0
643 C 以下进行的是输入受载网格单元的编号及其受载边数。
644 755 FORMAT(2X,'网格单元编号为：',I4,', 单元受载边数为：',I3)
645 DO 20 I=1,NUMBER
646 READ(5,*) KEL,NEDGE
647 NLOAD(KEL)=NEDGE
648 WRITE(6,755) KEL,NLOAD(KEL)
649 20 CONTINUE
650 DO 200 IELEM=1,NELEM
651 NEDGE=NLOAD(IELEM)
652 C 如果当前单元受载边数为零,则转移到循环的终端语句进行下一个网格单元。
653 IF(NEDGE.EQ.0) GO TO 200
654 C 以下进行的是输入单元受载信息标题。
655 READ(5,5000) LINECHAR
656 LOCATECHAR=INDEX(LINECHAR,'输入')
657 IF(LOCATECHAR.NE.0) THEN
658 LINECHAR(LOCATECHAR:LOCATECHAR+3)='输出'
659 ENDIF
660 WRITE(6,5000) LINECHAR
661 C 以下进行的是单元等效节点载荷列阵赋初值零。
662 DO 50 I=1,16
663 ELOAD(I)=0.0
664 50 CONTINUE
665 C 以下进行的是输入受载边三个节点的总体编号及当前受载边均布载荷σ和τ。
666 DO 160 IEDGE=1,NEDGE
667 READ(5,*) (NOPRS(I),I=1,3)
668 READ(5,*) SGMAR,TAU
669 C 如果均布载荷为零,则输入三个节点的法向载荷σ(1,2,3)及切向载荷τ(1,2,3)。
670 IF(ABS(SGMAR).LT.1.D-8.AND.ABS(TAU).LT.1.D-8) THEN
671 READ(5,*) ((PRESS(J,I),J=1,3),I=1,2)
672 ELSE
673 C 如果均布载荷不为零,则给三个节点赋值。
674 DO 60 I=1,3
675 PRESS(I,1)=SGMAR
676 PRESS(I,2)=TAU
677 60 CONTINUE
678 END IF
679 C 给当前受载边三个节点坐标赋值。
680 DO 100 I=1,3
681 LNODE=NOPRS(I)
682 DO 100 J=1,2
683 ELCOD(J,I)=COORD(J,LNODE)
684 100 CONTINUE
685 T=-1.0
686 DO 150 IGAUS=1,3
687 S=POSGP(IGAUS)
688 CALL SFR2(S,T)
689 DO 110 J=1,2
690 PGASH(J)=0.0
691 DGASH(J)=0.0
692 C 以下进行的是计算当前积分点的σ和τ及偏导数在当前积分点的值。
693 DO 110 I=1,3
694 PGASH(J)=PGASH(J)+PRESS(I,J)*SHAPE(I)
695 DGASH(J)=DGASH(J)+ELCOD(J,I)*DERIV(1,I)
696 110 CONTINUE
697 DVOLU=WEIGP(IGAUS)
698 PXCOM=DGASH(1)*PGASH(2)-DGASH(2)*PGASH(1)
699 PYCOM=DGASH(1)*PGASH(1)+DGASH(2)*PGASH(2)
700 C 查找当前受载边的第I个节点是当前单元的第几个节点,查找结果为第I个节点。
701 DO 120 I=1,8
702 NLOCA=LNODS(I,IELEM)
703 IF(NLOCA.EQ.NOPRS(1)) GO TO 130
704 120 CONTINUE
705 130 J=I+2
706 KOUNT=0
707 DO 140 K=I,J
708 KOUNT=KOUNT+1
709 N=(K-1)*2+1
710 M=N+1
711 C 如果单元节点局部编号大于8,则应等于1。
712 IF(K.GT.8)THEN
713 N=1
714 M=2
715 END IF
716 C 以下进行的是累加得到当前单元的等效节点载荷。
717 ELOAD(N)=ELOAD(N)+SHAPE(KOUNT)*PXCOM*DVOLU
718 ELOAD(M)=ELOAD(M)+SHAPE(KOUNT)*PYCOM*DVOLU
719 140 CONTINUE
720 150 CONTINUE
721 160 CONTINUE
722 WRITE(6,795) ELOAD
723 795 FORMAT(2X,'单元等效载荷列阵为：'/4(1P4E15.5/))
724 WRITE(8)ELOAD
725 200 CONTINUE
726 RETURN
727 END
728 C*****************************************************************************
729 C 这是一个形成整体刚度矩阵和整体载荷列阵的单元组集子程序。 *
730 C*****************************************************************************
731 SUBROUTINE ASEM(V,A,LNODS,IFPRE,NLOAD,MAXA)
732 IMPLICIT REAL*8(A-H,O-Z)
733 DIMENSION V(1),A(1)
734 DIMENSION LNODS(8,1),IFPRE(1),NLOAD(1),MAXA(1)
735 DIMENSION NCODF(2),ESTIF(136),ELOAD(16)
736 COMMON/SOL/NPOIN,NELEM,NTYPE,NMATS
737 MAXA(1)=1
738 DO 100 IPOIN=1,NPOIN
739 KMIN=NPOIN+1
740 DO 40 IELEM=1,NELEM
741 DO 10 I=1,8
742 LNODE=LNODS(I,IELEM)
743 IF(LNODE.EQ.IPOIN)GO TO 20
744 10 CONTINUE
745 GO TO 40
746 20 CONTINUE
747 DO 30 I=1,8
748 LNODE=LNODS(I,IELEM)
749 IF(LNODE.GE.KMIN)GO TO 30
750 KMIN=LNODE
751 30 CONTINUE
752 40 CONTINUE
753 KHIGH=(IPOIN-KMIN)*2
754 DO 60 J=1,2
755 KDOFN=(IPOIN-1)*2+J+1
756 MAXA(KDOFN)=MAXA(KDOFN-1)+KHIGH+J
757 60 CONTINUE
758 100 CONTINUE
759 NN=NPOIN*2
760 NNM=NN+1
761 NWK=MAXA(NNM)-1
762 WRITE(6,735) NWK
763 735 FORMAT(2X,'剖面元素之和=',I7)
764 DO 120 IWK=1,NWK
765 A(IWK)=0.0
766 120 CONTINUE
767 DO 130 IN=1,NN
768 130 V(IN)=0.0
769 DLARG=1.0D+15
770 REWIND 8
771 DO 300 IELEM=1,NELEM
772 WRITE(6,135) IELEM
773 135 FORMAT(2X,'网格单元编号为：',I5)
774 READ(7) ESTIF
775 IF(NLOAD(IELEM).GT.0)THEN
776 READ(8) ELOAD
777 END IF
778 DO 280 INODE=1,8
779 LNODI=LNODS(INODE,IELEM)
780 ICODF=IFPRE(LNODI)
781 NCODF(1)=ICODF/10
782 NCODF(2)=ICODF-NCODF(1)*10
783 DO 260 IDOFN=1,2
784 IEVAB=(INODE-1)*2+IDOFN
785 ICOLN=(LNODI-1)*2+IDOFN
786 DO 240 JNODE=1,8
787 LNODJ=LNODS(JNODE,IELEM)
788 DO 220 JDOFN=1,2
789 JEVAB=(JNODE-1)*2+JDOFN
790 JROW=(LNODJ-1)*2+JDOFN
791 IF(JROW.GT.ICOLN)GO TO 220
792 IF(JEVAB.GT.IEVAB)THEN
793 KEVAB=JEVAB*(JEVAB-1)/2+IEVAB
794 ELSE
795 KEVAB=IEVAB*(IEVAB-1)/2+JEVAB
796 END IF
797 LDIS=MAXA(ICOLN)+(ICOLN-JROW)
798 A(LDIS)=A(LDIS)+ESTIF(KEVAB)
799 220 CONTINUE
800 240 CONTINUE
801 IF(NLOAD(IELEM).GT.0)THEN
802 V(ICOLN)=V(ICOLN)+ELOAD(IEVAB)
803 END IF
804 IF(NCODF(IDOFN).EQ.0)GO TO 260
805 KEVAB=IEVAB*(IEVAB+1)/2
806 LDIS=MAXA(ICOLN)
807 A(LDIS)=A(LDIS)+DLARG*ESTIF(KEVAB)
808 260 CONTINUE
809 280 CONTINUE
810 300 CONTINUE
811 C*******
812 C ENTER GIVED DISPLACEMENTS
813 C*******
814 DO 400 I=1,NPOIN
815 ICODF=IFPRE(I)
816 NCODF(1)=ICODF/10
817 NCODF(2)=ICODF-NCODF(1)*10
818 DO 320 J=1,2
819 IF(NCODF(J).EQ.2)THEN
820 C*******
821 READ(5,*) X
822 C*******
823 ICOLN=(I-1)*2+J
824 LDIS=MAXA(ICOLN)
825 V(ICOLN)=A(LDIS)*X
826 END IF
827 320 CONTINUE
828 400 CONTINUE
829 RETURN
830 END
831 C*****************************************************************************
832 C 这是一个求解结构平衡方程组的子程序。 *
833 C*****************************************************************************
834 SUBROUTINE SOLVE(V,A,MAXA)
835 IMPLICIT REAL*8(A-H,O-Z)
836 DIMENSION V(1),A(1),MAXA(1)
837 COMMON/SOL/NPOIN,NELEM,NTYPE,NMATS
838 K=0
839 NN=NPOIN*2
840 NNM=NN+1
841 NWK=MAXA(NNM)-1
842 DO 640 N=1,NN
843 KN=MAXA(N)
844 KL=KN+1
845 KU=MAXA(N+1)-1
846 KH=KU-KL
847 IF(KH)610,590,550
848 550 K=N-KH
849 IC=0
850 KLT=KU
851 DO 580 J=1,KH
852 IC=IC+1
853 KLT=KLT-1
854 KI=MAXA(K)
855 ND=MAXA(K+1)-KI-1
856 IF(ND)580,580,560
857 560 KK=MIN(IC,ND)
858 C=0
859 DO 570 L=1,KK
860 M1=KI+L
861 M2=KLT+L
862 570 C=C+A(M1)*A(M2)
863 A(KLT)=A(KLT)-C
864 580 K=K+1
865 590 K=N
866 B=0.0
867 DO 600 KK=KL,KU
868 K=K-1
869 KI=MAXA(K)
870 C=A(KK)/A(KI)
871 B=B+C*A(KK)
872 600 A(KK)=C
873 A(KN)=A(KN)-B
874 610 IF(A(KN)) 620,620,640
875 620 WRITE(6,2000) N,A(KN)
876 STOP
877 640 CONTINUE
878 C********
879 C REDUCE RIGHT--HAND--SIDE LOAD VECTOR
880 C********
881 650 CONTINUE
882 DO 680 N=1,NN
883 KL=MAXA(N)+1
884 KU=MAXA(N+1)-1
885 KUL=KU-KL
886 IF(KUL)680,660,660
887 660 K=N
888 C=0.0
889 DO 670 KK=KL,KU
890 K=K-1
891 670 C=C+A(KK)*V(K)
892 V(N)=V(N)-C
893 680 CONTINUE
894 C*******
895 C BACK-SUBSTITUTE
896 C********
897 DO 700 N=1,NN
898 K=MAXA(N)
899 V(N)=V(N)/A(K)
900 700 CONTINUE
901 IF(NN.EQ.1)GO TO 740
902 N=NN
903 DO 730 L=2,NN
904 KL=MAXA(N)+1
905 KU=MAXA(N+1)-1
906 KUL=KU-KL
907 IF(KUL)730,710,710
908 710 K=N
909 DO 720 KK=KL,KU
910 K=K-1
911 720 V(K)=V(K)-A(KK)*V(N)
912 730 N=N-1
913 740 CONTINUE
914 2000 FORMAT(//10X,'***** 停止运行！系数矩阵非正定！*****'//'NON
915 + POSITIVE PIVOT FOR EQUATION'I4//' PIVOT=',E20.12)
916 C********
917 C OUTPUT DISPLACEMENT
918 C********
919 WRITE(6,903)
920 903 FORMAT(//80('=')/36X,'移动输出'/80('=')//)
921 WRITE(6,905)
922 905 FORMAT(32X,'Ｘ方向',10X,'Ｙ方向')
923 M2=0
924 DO 1000 I=1,NPOIN
925 M1=M2+1
926 M2=M1+1
927 1000 WRITE(6,915)I,(V(J),J=M1,M2)
928 915 FORMAT(2X,'高斯积分点编号为：',I5,1PE16.5,1PE16.5)
929 RETURN
930 END
931 C*****************************************************************************
932 C 这是一个计算单元应力的子程序。 *
933 C*****************************************************************************
934 SUBROUTINE STRE(COORD,PROPS,V,LNODS,MATNO)
935 IMPLICIT REAL*8(A-H,O-Z)
936 CHARACTER*80 LINECHAR
937 DIMENSION COORD(2,1),PROPS(3,1),V(1),LNODS(8,1),MATNO(1)
938 DIMENSION POSGP(3),ELCOD(2,8)
939 DIMENSION ELDIS(16),BE(3),S(3),SP(2)
940 COMMON/SOL/NPOIN,NELEM,NTYPE,NMATS
941 COMMON/BMDMP/BMATX(3,16),DMATX(3,3)
942 COMMON/GCOD/GPCOD(2,9)
943 POSGP(1)=-DSQRT(0.6D0)
944 POSGP(2)=0.0
945 POSGP(3)=DSQRT(0.6D0)
946 WRITE(6,616)
947 616 FORMAT(80('=')/26X,'以下是高斯积分点的应力输出'/80('=')/)
948 DO 90 IELEM=1,NELEM
949 WRITE(6,636) IELEM
950 636 FORMAT(2X,'网格单元编号为：',I5)
951 LPROP=MATNO(IELEM)
952 DO 20 I=1,8
953 LNODE=LNODS(I,IELEM)
954 LDOFN=LNODE*2-2
955 DO 10 J=1,2
956 IEVAB=I*2-2+J
957 MDOFN=LDOFN+J
958 ELDIS(IEVAB)=V(MDOFN)
959 ELCOD(J,I)=COORD(J,LNODE)
960 10 CONTINUE
961 20 CONTINUE
962 YOUNG=PROPS(1,LPROP)
963 POISS=PROPS(2,LPROP)
964 THICK=PROPS(3,LPROP)
965 CALL MODPS(YOUNG,POISS)
966 KGASP=0
967 DO 80 IGAUS=1,3
968 DO 80 JGAUS=1,3
969 KGASP=KGASP+1
970 EXISP=POSGP(IGAUS)
971 ETASP=POSGP(JGAUS)
972 CALL SFR2(EXISP,ETASP)
973 CALL JACOB2(IELEM,DJACB,KGASP,ELCOD)
974 CALL BMATPS
975 DO 40 I=1,3
976 BE(I)=0.0
977 DO 30 J=1,16
978 BE(I)=BE(I)+BMATX(I,J)*ELDIS(J)
979 30 CONTINUE
980 40 CONTINUE
981 DO 60 L=1,3
982 S(L)=0.0
983 DO 50 J=1,3
984 S(L)=S(L)+DMATX(L,J)*BE(J)
985 50 CONTINUE
986 60 CONTINUE
987 WRITE(6,755) KGASP
988 755 FORMAT(/2X,'高斯积分点的编号为：',I2)
989 LINECHAR='高斯积分点的坐标为：'
990 WRITE(6,775) LINECHAR,GPCOD(1,KGASP),GPCOD(2,KGASP)
991 775 FORMAT(2X,A26,' X=',F9.3,5X,'Y=',F9.3)
992 C *****
993 XA=(S(1)+S(2))*0.5
994 XI=(S(1)-S(2))*0.5
995 XE=S(3)
996 X0=DSQRT(XI*XI+XE*XE)
997 SP(1)=XA+X0
998 SP(2)=XA-X0
999 C ******
1000 WRITE(6,815)
1001 815 FORMAT(14X,'正应力σx',6X,'正应力σy',7X,'剪应力τ')
1002 WRITE(6,835) S
1003 835 FORMAT(10X,1PE13.4,2X,1PE13.4,2X,1PE13.4)
1004 WRITE(6,855)
1005 855 FORMAT(14X,'主应力σ1',6X,'主应力σ2')
1006 WRITE(6,875) SP
1007 875 FORMAT(10X,1PE13.4,2X,1PE13.4)
1008 80 CONTINUE
1009 90 CONTINUE
1010 RETURN
1011 END
1012
