-
Notifications
You must be signed in to change notification settings - Fork 0
/
stress-pro.f90
530 lines (520 loc) · 11.7 KB
/
stress-pro.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
!此程序是对原计算平面应力问题的改进，细化了网格（5x5）
!本程序计算弹性力学平面应力问题（带厚度），在F90上调试通过
!包含主程序STRESS和四个子程序ELSTMX(KK)、MODIFY、DCMPBD和SLVBD
!主程序主要进行数据输入、结果输出、最大半带宽的计算、组装总体刚度矩阵以及计算单元内
!的应力、应变值
!子程序ELSTMX(KK)计算刚度矩阵系数，生成第KK单元的刚度矩阵K
!子程序MODIFY输入载荷节点的载荷值，位移边界节点处的位移值，对总体刚度矩阵、位移数组
!和节点力数组进行修改
!子程序DEMPBD用高斯消元法将对称等带宽总刚度矩阵的上半部分化为上三角阵
!子程序SLVBD对子程序DEMPBD得到的上三角矩阵进行回代计算，得到节点位移向量DD
!
PROGRAM STRESS
COMMON /ELMATX/ESM(6,6),X(3),Y(3),D(3,3)
COMMON / GRAD/ B(3,6),AR2
COMMON /MTL/ EM,PR,TH
COMMON /AV/A(8500),JGF,JGSM,NP,NBW,JEND
COMMON /PRO/IRATION,IROWNUM,ICOLNUM,NN,NE
DIMENSION NS(6),U(6),STRA(3),STRE(3)
DATA IN/60/,IO/61/
!
!将单元以及组成单元的节点数组按照动态数组分配方法
REAL,ALLOCATABLE::XC(:),YC(:)
INTEGER,ALLOCATABLE::NEL(:,:)
! INTEGER :: RATION
! INTEGER :: ROWNUM
! INTEGER :: COLNUM
INTEGER::A1
INTEGER::B1
INTEGER::C1
INTEGER::A0
INTEGER::B0
INTEGER::C0
INTEGER::FLAG
!
!ESM(6,6)――单元刚度矩阵，X(3),Y(3)――单元节点坐标，D(3,3)材料性质矩阵
!B(3,6)――几何矩阵，AR2――三角形面积的二倍，NP――自由度总数，NBW――最大半带宽
!A(8500)――存储节点位移向量，节点力向量和总体刚度矩阵的数组A
!JGF，JEND为计数单元
!JGF=NP-节点位移向量在数组A中的位置，JGSM=JGF+NP-节点力向量在数组A中的位置
!JEND=JGSM+NP*NBW――刚度矩阵在数组A中的位置，数组A总长度
!NS(6)――一个单元节点的自由度编号数组，U(6)――一个单元的节点自由度
!IN/60/,IO/61/一个输入输出文件设备号，STRA(3),STRE(3)――存储单元的应变、应力
!
!---------------------------程序输入段-------------------------------------
! TITLE――存储计算内容标题
! NN――节点总数，NE――单元总数，EM――杨氏模量，PR――泊松比，TH――板的厚度，N――单元编号
! NEL(I)――节点的X轴坐标，YC(I)――节点的Y轴坐标
! NEL(N,I)――组成第N个三角形的第I个节点的编号(I=1,2,3)
!
OPEN(60,FILE='INPUT.DAT',STATUS='UNKNOWN')
OPEN(61,FILE='OUTPUT.DAT',STATUS='UNKNOWN')
READ(IN,1) TITLE
1 FORMAT(900A)
WRITE(IO,1) TITLE
WRITE(IO,*)
!
! 输入计算模型的节点总数NN和单元总数NE
! READ(IN,*)NN
! READ(IN,*)NE
IRATION=5
NN=(1+2*IRATION+1)*(2*IRATION+1)/2
NE=((2*(2*IRATION)-1)+1)*(2*IRATION)/2
WRITE(IO,2)NN,NE
2 FORMAT(/,'输入数据:',/,'节点数=',I3,5X,'单元数=',I4,/)
NP=2*NN
ALLOCATE(NEL(1:NE,1:3),XC(NN),YC(NN))
!
! 输入材料的杨氏模量EM，泊松比PR，平板厚度TH，节点坐标XC(I),YC(I)和组成单元的
! 节点NEL(N,I)
! 组成单元的节点的编号都按逆时针顺序输入
!
READ(IN,*)EM
READ(IN,*)PR
READ(IN,*)TH
!这里的节点坐标采用程序自动输入，需要输入所划分区域的最外侧顶点的坐标
!利用一个SETCOOR函数实现
XLENGTH=2./(2*IRATION)
YLENGTH=2./(2*IRATION)
IROWNUM=2*IRATION+1
ICOLNUM=2*IRATION+1
DO 456 I=1,IROWNUM
DO 345 J=1,I
XC(I*(I-1)/2+J)=(J-1)*XLENGTH
YC(I*(I-1)/2+J)=2.0-(I-1)*YLENGTH
345 CONTINUE
456 CONTINUE
!
! 输出材料的性质和计算模型拓扑数据便于检查时对照
!
WRITE(IO,3)
3 FORMAT('材料常数为:')
WRITE(IO,4)EM,PR,TH
4 FORMAT('弹性模量为:',E12.5,5X,'泊松比为:',E12.5,5X,'厚度为:',E12.5)
WRITE(IO,*)
DO 6 I=1,NN
WRITE(IO,5)I,XC(I),YC(I)
5 FORMAT('节点',I3,'坐标为:X=',F8.3,3X,'Y=',F8.3)
6 CONTINUE
WRITE(IO,*)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!获得NEL数组的值，NEL数组存储单元与其节点的对应关系
A1=1
B1=2
C1=3
A0=1
B0=2
C0=3
NEL(1,1)=A1
NEL(1,2)=B1
NEL(1,3)=C1
M=2
FLAG=1
I=1
GOTO 91
90 IF(MOD(FLAG,2).EQ.1) THEN
B1=C1
C1=A1+1
NEL(M,1)=A1
NEL(M,2)=B1
NEL(M,3)=C1
FLAG=FLAG+1
M=M+1
GOTO 91
ENDIF
IF(MOD(FLAG,2).EQ.0) THEN
A1=C1
C1=B1+1
NEL(M,1)=A1
NEL(M,2)=B1
NEL(M,3)=C1
M=M+1
FLAG=FLAG+1
GOTO 91
ENDIF
91 IF(C1.LT.((I+1)**2+(I+1))/2) THEN
GOTO 90
ENDIF
IF(C1.GE.((I+1)**2+(I+1))/2) THEN
IF((I+1).LT.IROWNUM) THEN
A0=B0
B0=B0+I+1
C0=B0+1
A1=A0
B1=B0
C1=C0
NEL(M,1)=A1
NEL(M,2)=B1
NEL(M,3)=C1
FLAG=1
I=I+1
M=M+1
GOTO 91
ENDIF
ENDIF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
WRITE(IO,*)
!
! ---------------------------计算开始-------------------------------
! 计算最大半带宽，B=MAX(D+1)*F
! D是一个单元各节点遍点号之差的最大值，F是一个节点的自由度数
!
INBW=0
NBW=0
DO 20 KK=1,NE
DO 25 I=1,3
25 NS(I)=NEL(KK,I)
DO 21 I=1,2
IJ=I+1
DO 21 J=IJ,3
NB=IABS(NS(I)-NS(J)) !节点号之差的绝对值
IF(NB.LE.NBW) GOTO 21
INBW=KK
NBW=NB
21 CONTINUE
20 CONTINUE
NBW=(NBW+1)*2 !平面问题节点自由度F=2,NBW此时为最大半带宽
!
!数组A中数据的安排:A(1,2,3....NP|NP+1.....2NP|2NP+1.....JEND)
! 节点位移向量 节点力向量 总体刚度矩阵
!初始化数组A
JGF=NP !JGF=2*NN
JGSM=JGF+NP !JGSM=4*NN
JEND=JGSM+NP*NBW !JEND=4*NN+2*NN*12
! NP为K中KNN的N，NBW为最大半带宽
! |XXXX|
! |XXXX|
! |XXXX|
! |XXX*|
! |XX**|
! |X***|
!用等带宽二维数组方法存储的刚度矩阵中包含*表示的位置，但在实际程序中未使用
!
JL=JEND-JGF
DO 24 I=1,JEND
24 A(I)=0.0
GOTO 30
!
! 生成材料性质矩阵
! E为杨氏模量
! V为泊松比
!
30 R=EM/(1.-PR**2)
D(1,1)=R
D(2,2)=D(1,1)
D(3,3)=R*(1.-PR)/2.
D(1,2)=PR*R
D(2,1)=D(1,2)
D(1,3)=0.0
D(3,1)=0.0
D(2,3)=0.0
D(3,2)=0.0
!
! 单元矩阵循环的开始
!
KK=1
! KK为单元号
! 节点自由度的生成，节点坐标的局部变化
!
32 DO 31 I=1,3
J=NEL(KK,I) !元件KK的三个节点编号
NS(2*I-1)=J*2-1
NS(2*I)=J*2 !元件KK各个自由度在总刚度矩阵中的方程号
X(I)=XC(J)
31 Y(I)=YC(J) !元件KK的节点坐标值
!
! 调用子程序ELSTMX计算单元刚度矩阵ESM(6,6)并输出
!
CALL ELSTMX(KK)
!
! 单元刚度矩阵组装成总体刚度矩阵
DO 33 I=1,6
II=NS(I)
DO 34 J=1,6
JJ=NS(J)+1-II
IF(JJ.LE.0)GOTO 34
J1=JGSM+(JJ-1)*NP+II-(JJ-1)*(JJ-2)/2
A(J1)=A(J1)+ESM(I,J)
34 CONTINUE
33 CONTINUE
KK=KK+1
IF(KK.LE.NE)GOTO 32
!
! 调用子程序MODIFY输入载荷节点处的载荷值，位移边界节点处的位移值，对总体刚度矩阵
! 位移数组和节点力数组进行相应的修改
CALL MODIFY
WRITE(IO,35)
35 FORMAT(/,'计算结果为:',/)
WRITE(IO,36)NBW
36 FORMAT('最大半带宽为',I3)
WRITE(IO,37)JEND
37 FORMAT('总数组大小为:',I4)
!
!调用子程序DEMPBD，用高斯消元法将对称等带宽总刚度矩阵化为上三角矩阵
CALL DCMPBD
!
! 调用子程序SLVBD对子程序DEMPBD得到的上三角矩阵进行回代计算，得到节点位移向量DD
CALL SLVBD
!
! ---------------------------输出各个节点的位移向量----------------------------
!
WRITE(IO,*)
DO 45 I=1,NP/2
WRITE(IO,43)I,A(2*I-1),A(2*I)
43 FORMAT('节点号',I3,5X,'X方向的位移UX=',E12.5,5X,'Y方向的位移UY=',E12.5)
45 CONTINUE
!
! ---------------------------附加计算-----------------------------------------
! 计算节点处的应变
!
DO 96 KK=1,NE
! 生成节点的自由度
! 节点坐标局部变化
DO 51 I=1,3
J=NEL(KK,I)
NS(2*I-1)=2*J-1
NS(2*I)=2*J
X(I)=XC(J)
51 Y(I)=YC(J)
!
! 单元节点位移局部变化
!
65 DO 73 I=1,6,2
NS1=NS(I)
NS2=NS(I+1)
U(I)=A(NS1) !节点X方向的位移量
73 U(I+1)=A(NS2) !节点Y方向的位移量
!
!计算单元应变
!
CALL ELSTMX(KK)
DO 52 I=1,3
STRA(I)=0.0
DO 52 K=1,6
52 STRA(I)=STRA(I)+B(I,K)*U(K)/AR2
!
! 计算单元应力值
!
DO 58 I=1,3
STRE(I)=0
DO 58 K=1,3
58 STRE(I)=STRE(I)+D(I,K)*(STRA(K))
!
! 计算主应力S1，S2，TM
!
AA=(STRE(1)+STRE(2))/2.
AB=SQRT((ABS(STRE(1)-STRE(2))/2.)**2+STRE(3)**2)
S1=AA+AB
S2=AA-AB
TM=AB
!
! 计算主应力方向与X轴的夹角
!
IF(ABS(STRE(1)-STRE(2)).LT.0.0001) GOTO 93
AC=ATAN2(2*STRE(3),(STRE(1)-STRE(2)))
THM=(180/3.1415926*AC)/2
GOTO 94
93 THM=90
!
! ------------------------输出单元应变和应力的计算结果--------------------------
!
94 WRITE(IO,57)KK
57 FORMAT(/,'单元',I4)
WRITE(IO,95)STRA(1),STRA(2),STRA(3)
95 FORMAT('EPTOX=',E12.5,2X,'EPTOY=',E12.5,2X,'EPTOXY=',E12.5)
WRITE(IO,97)STRE(1),STRE(2),STRE(3)
97 FORMAT('SX=',E12.5,5X,'SY=',E12.5,5X,'SXY=',E12.5)
WRITE(IO,98)S1,S2,TM
98 FORMAT('S1=',E12.5,5X,'S2=',E12.5,5X,'TMAX=',E12.5)
WRITE(IO,99)THM
99 FORMAT(44X,'ANGEL=',F8.2,'°')
96 CONTINUE
CLOSE(60)
CLOSE(61)
STOP
END
!
! 计算刚度系数矩阵，计算第KK单元的刚度矩阵K
!
SUBROUTINE ELSTMX(KK)
COMMON/MTL/EM,PR,TH
COMMON/GRAD/B(3,6),AR2
COMMON/ELMATX/ESM(6,6),X(3),Y(3),D(3,3)
DIMENSION C(6,3)
DATA IO/61/ IN/60/
!
! 计算矩阵B
!
DO 20 I=1,3
DO 20 J=1,6
20 B(I,J)=0.0
B(1,1)=Y(2)-Y(3)
B(1,3)=Y(3)-Y(1)
B(1,5)=Y(1)-Y(2)
B(2,2)=X(3)-X(2)
B(2,4)=X(1)-X(3)
B(2,6)=X(2)-X(1)
B(3,1)=B(2,2)
B(3,2)=B(1,1)
B(3,3)=B(2,4)
B(3,4)=B(1,3)
B(3,5)=B(2,6)
B(3,6)=B(1,5)
AR2=X(2)*Y(3)+X(3)*Y(1)+X(1)*Y(2)-X(2)*Y(1)-X(3)*Y(2)-X(1)*Y(3)
!
! AR2=2*DLT
!
! 计算矩阵C=(B3*6)*(D3*3)
DO 22 I=1,6
DO 22 J=1,3
C(I,J)=0.0
DO 22 K=1,3
22 C(I,J)=C(I,J)+B(K,I)*D(K,J)
!
! 计算矩阵ESM=DTDB=CB
!
DO 27 I=1,6
DO 27 J=1,6
SUM=0.0
DO 28 K=1,3
28 SUM=SUM+C(I,K)*B(K,J)
ESM(I,J)=SUM*TH/(2.*AR2)
27 CONTINUE
RETURN
END
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!输入载荷节点处的载荷值，位移边界节点处的位移值，对总体刚度矩阵、位移数组和节点力
! 数组进行修改
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE MODIFY
COMMON/AV/A(8500),JGF,JGSM,NP,NBW,JEND
COMMON /PRO/IRATION,IROWNUM,ICOLNUM,NN,NE
DATA IN/60/,IO/61/
INTEGER :: ROW
INTEGER :: COL
INTEGER :: STA
!
!输入节点的集中载荷，放到A数组中节点力向量的相应位置
!IB为施加外载荷节点的自由度(X,Y),BV为该自由度上的载荷值
202 READ(IN,*)IB
IF(IB.LE.0) THEN
ROW=0
COL=0
GOTO 208
ENDIF
READ(IN,*)BV
IF(MOD(IB,2).EQ.1) GOTO 204
WRITE(IO,203) IB/2,BV
203 FORMAT('节点',I3,'载荷为:PY=',F11.3)
GOTO 206
204 WRITE(IO,205) IB/2+1,BV
205 FORMAT('节点',I3,'载荷为:PX=',F8.3)
206 A(JGF+IB)=A(JGF+IB)+BV
GOTO 202
!
!输入位移边界节点处的位移值，放到A数组中节点位移向量的相对位置
!IB为节点位移的自由度，BV为位移值
!IB为基奇数时代表x，偶数代表y
208 IF(ROW.LT.IROWNUM) GOTO 11
!下面的一横排的y方向约束
ROW=IROWNUM
COL=COL+1
STA=ROW*(ROW-1)/2.0+COL
IB=2*STA
GOTO 12
!左侧的一竖排x方向约束为0
11 COL=1
ROW=ROW+1
STA=ROW*(ROW-1)/2.0+COL
IB=2*STA-1
COL=0
12 IF(COL.GT.ICOLNUM) RETURN
BV=0.0
IF(MOD(IB,2).EQ.1) GOTO 214
WRITE(IO,213) (IB+1)/2,BV
213 FORMAT('节点',I3,'位移约束为:V=',F8.3)
GOTO 209
214 WRITE(IO,215)(IB+1)/2,BV
215 FORMAT('节点',I3,'位移约束为:U=',F8.3)
209 K=IB-1
DO 211 J=2,NBW
M=IB+J-1
IF(M.GT.NP) GOTO 210
IJ=JGSM+(J-1)*NP+IB-(J-1)*(J-2)/2
A(JGF+M)=A(JGF+M)-A(IJ)*BV
A(IJ)=0.0
210 IF(K.LE.0) GOTO 211
KJ=JGSM+(J-1)*NP+K-(J-1)*(J-2)/2
A(JGF+k)=A(JGF+K)-A(KJ)*BV
A(KJ)=0.0
K=K-1
211 CONTINUE
A(JGF+IB)=A(JGSM+IB)*BV
221 CONTINUE
GOTO 208
END
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! GUASS 消元子程序
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE DCMPBD
COMMON /AV/A(8500),JGF,JGSM,NP,NBW,JEND
!
!用高斯消元法将对称等带宽刚度矩阵化为上三角矩阵
!
NP1=NP-1
DO 226 I=1,NP1
MJ=I+NBW-1
IF(MJ.GT.NP) MJ=NP
NJ=I+1
MK=NBW
IF((NP-I+1).LT.NBW) MK=NP-I+1
ND=0
DO 225 J=NJ,MJ
MK=MK-1
ND=ND+1
NL=ND+1
DO 225 K=1,MK
NK=ND+k
JK=JGSM+(K-1)*NP+J-(K-1)*(K-2)/2
INL=JGSM+(NL-1)*NP+I-(NL-1)*(NL-2)/2
INK=JGSM+(NK-1)*NP+I-(NK-1)*(NK-2)/2
II=JGSM+I
225 A(JK)=A(JK)-A(INL)*A(INK)/A(II)
226 CONTINUE
RETURN
END
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 回代求解子程序
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE SLVBD
COMMON /AV/A(8500),JGF,JGSM,NP,NBW,JEND
DATA IN/60/,IO/61/
NP1=NP-1
!
! 对子程序DEMPBD得到的上三角矩阵进行回代计算，得到节点位移向量DD
!
DO 250 I=1,NP1
MJ=I+NBW-1
IF(MJ.GT.NP)MJ=NP
NJ=I+1
L=1
DO 250 J=NJ,MJ
L=L+1
IL=JGSM+(L-1)*NP+I-(L-1)*(L-2)/2
250 A(JGF+J)=A(JGF+J)-A(IL)*A(JGF+I)/A(JGSM+I)
!
!求解节点自由度的回代计算，从下往上迭代
!
A(NP)=A(JGF+NP)/A(JGSM+NP)
DO 252 K=1,NP1
I=NP-K
MJ=NBW
IF((I+NBW-1).GT.NP) MJ=NP-I+1
SUM=0.0
DO 251 J=2,MJ
N=I+J-1
IJ=JGSM+(J-1)*NP+I-(J-1)*(J-2)/2
251 SUM=SUM+A(IJ)*A(N)
252 A(I)=(A(JGF+I)-SUM)/A(JGSM+I)
RETURN
END