44
45
46
47#include "implicit_f.inc"
48#include "comlock.inc"
49
50
51
52#include "mvsiz_p.inc"
53
54
55
56#include "com01_c.inc"
57#include "param_c.inc"
58#include "com08_c.inc"
59#include "scr17_c.inc"
60#include "units_c.inc"
61
62
63
64 INTEGER, INTENT(IN) ::
65 INTEGER, INTENT(IN) :: JSPH
66 INTEGER MAT(MVSIZ),NGL(MVSIZ),NEL
67
69 . pm(npropm,*), off(*), sig(nel,6), eint(*), pla(*), sigf(*),
70 . rx(*), ry(*), rz(*), sx(*), sy(*), sz(*),
71 . epsf(*), dam(nel,5), epe(nel,3), epc(nel,3), a(mvsiz,6), vol(*),
72 . vnew(mvsiz), dvol(mvsiz), ssp(mvsiz),
73 . d1(mvsiz), d2(mvsiz), d3(mvsiz),d4(mvsiz),d5(mvsiz),d6(mvsiz),
74 . sold1(mvsiz), sold2(mvsiz), sold3(mvsiz), sold4(mvsiz),
75 . sold5(mvsiz), sold6(mvsiz),sigy(*),defp(*),seq_output(*),
76 . tsaiwu(mvsiz)
77
78
79
80 INTEGER KD1(MVSIZ), KD2(MVSIZ), KD3(MVSIZ),
81 . KD4(MVSIZ),ICC(MVSIZ),ICC_1,
82 . I, IDAM, KDX, MX, NINDX, INDEX(), J
83
85 . ax(mvsiz), ay(mvsiz), az(mvsiz), bx(mvsiz), by(mvsiz),
86 . bz(mvsiz), cx(mvsiz), cy(mvsiz),cz(mvsiz),
87 . wvec(mvsiz),s1(mvsiz), s2(mvsiz), s3(mvsiz),
88 . t1(mvsiz), t2(mvsiz), t3(mvsiz),t4(mvsiz),t5(mvsiz),t6(mvsiz),
89 . e1(mvsiz), e2(mvsiz), e3(mvsiz),e4(mvsiz),e5(mvsiz),e6(mvsiz),
90 . hard(mvsiz), sigmy(mvsiz),
alpha(mvsiz), efib(mvsiz),
91 . epsft(mvsiz), epsfc(mvsiz), sigeq(mvsiz), p(mvsiz),
92 . d11(mvsiz), d12(mvsiz), d13(mvsiz), d22(mvsiz),
93 . d23(mvsiz), d33(mvsiz), g12(mvsiz), g23(mvsiz), g31(mvsiz),
94 . c11(mvsiz), c22(mvsiz), c33(mvsiz), c12(mvsiz), c23(mvsiz),
95 . c13(mvsiz), f11(mvsiz), f22(mvsiz), f44(mvsiz), f55
96 . f12(mvsiz), f23(mvsiz), f1(mvsiz), f2(mvsiz), f4(mvsiz),
97 . f5(mvsiz), delta(mvsiz), so1(mvsiz),
98 . so2(mvsiz), so3(mvsiz), so4(mvsiz), so5(mvsiz), so6(mvsiz),
99 . ds1(mvsiz), ds2(mvsiz), ds3(mvsiz), ds4(mvsiz), ds5(mvsiz),
100 . ds6(mvsiz), dp1(mvsiz), dp2(mvsiz), dp3(mvsiz), dp4(mvsiz),
101 . dp5(mvsiz), dp6(mvsiz), lamda(mvsiz), coef(mvsiz), plas(mvsiz),
102 . cn(mvsiz), cb(mvsiz), cnn(mvsiz), sigt1(mvsiz), sigt2(mvsiz),
103 . sigt3(mvsiz),cc(mvsiz),epdr(mvsiz)
105 . fib, ca, sigmx, epsp, dt5 ,sigym,
106 . cb_1,cn_1,cc_1,epdr_1,
107 . f1_1,f2_1,f4_1,f5_1,f11_1,
108 . f22_1,f44_1,f55_1,f12_1,f23_1,
109 . c11_1,c22_1,c33_1,c12_1,c23_1,
110 . c13_1,d11_1,d12_1,d13_1,d22_1,
111 . d23_1,d33_1,g12_1,g23_1,g31_1,
112 . alpha_1,wplaref,dwpla
113
114 IF (ncycle==0 .AND. irun==1) THEN
115 DO i=1,nel
116 kd1(i) = 0
117 kd2(i) = 0
118 kd3(i) = 0
119 kd4(i) = 0
120 ENDDO
121 ELSE
122 DO i=1,nel
123 idam = int(dam(i,5))-10000
124 kd1(i) = idam/1000
125 kdx = idam - kd1(i)*1000
126 kd2(i) = kdx/100
127 kdx = kdx - kd2(i)*100
128 kd3(i) = kdx/10
129 kd4(i) = kdx - kd3(i)*10
130 ENDDO
131 ENDIF
132
133
134
135 fib = zero
136 mx = mat(1)
137 alpha_1 = pm(39,mx)
138
139 DO i=1,nel
142 ENDDO
143
144
145
146
148 1 pm, a, rx, ry,
149 2 rz, sx, sy, sz,
150 3 ax, ay, az, bx,
151 4 by, bz, cx, cy,
152 5 cz, nel, jcvt, jsph)
153 CALL m14gtf(sig,ax ,ay ,az ,bx ,by,
154 2 bz ,cx ,cy ,cz ,d1 ,d2,
155 3 d3 ,d4 ,d5 ,d6 ,t1 ,t2,
156 4 t3 ,t4 ,t5 ,t6 ,e1 ,e2,
157 5 e3 ,e4 ,e5 ,e6 ,nel)
158
159
160
161 nindx=0
162 mx =mat(1)
163 cb_1 =pm(26,mx)
164 cn_1 =pm(27,mx)
165 cc_1 =pm(50,mx)
166 epdr_1 =pm(51,mx)
167 icc_1 =nint(pm(52,mx))
168
169 DO i=1,nel
170 ca =pm(25,mx)
171 cb(i) =cb_1
172 cn(i) =cn_1
173 cc(i) =cc_1
174 epdr(i) =epdr_1
175 icc(i) =icc_1
176 epsp =
max(abs(e1(i)),abs(e2(i)),abs(e3(i)),abs(e4(i)),abs(e5(i)))
177 IF (epsp>epdr(i) .AND. cc(i)/=zero) THEN
178 epsp=one + cc(i) * log(epdr(i)/epsp)
179 ELSE
180 epsp=one
181 ENDIF
182 IF (icc(i)==1) THEN
183 sigmx =pm(28,mx)*epsp
184 ELSEIF (icc(i)==2) THEN
185 sigmx =pm(28,mx)
186 ELSEIF (icc(i)==3) THEN
187 sigmx =pm(28,mx)*epsp
188
189 ELSEIF (icc(i)==4) THEN
190 sigmx =pm(28,mx)
191
192 ENDIF
193 cb(i) =cb(i)*epsp
194 ca =ca*epsp
195 sigmy(i)=
min(sigmx,ca+cb(i)*pla(i)**cn(i))
196 IF (sigmy(i)==sigmx .AND. off(i)==one) THEN
197 off(i)=zep99
198 kd4(i)=2
199 nindx=nindx+1
200 index(nindx)=i
201 ENDIF
202 sigt1(i)=pm(33,mx)
203 sigt2(i)=pm(34,mx)
204 sigt3(i)=pm(35,mx)
205 ssp(i) =pm(49,mx)
206 delta(i)=pm(78,mx)
207 ENDDO
208
209 IF (nindx/=0) THEN
210 DO j=1,nindx
211 i=index(j)
212#include "lockon.inc"
213 WRITE(iout,1000) ngl(i)
214#include "lockoff.inc"
215 END DO
216 END IF
217
218 DO i=1,nel
219 e1(i)=e1(i)*dt1
220 e2(i)=e2(i)*dt1
221 e3(i)=e3(i)*dt1
222 e4(i)=e4(i)*dt1
223 e5(i)=e5(i)*dt1
224 e6(i)=e6(i)*dt1
225 ENDDO
226
227 DO i=1,nel
228 epe(i,1)=epe(i,1)+e1(i)
229 epe(i,2)=epe(i,2)+e2(i)
230 epe(i,3)=epe(i,3)+e3(i)
231 ENDDO
232
233 mx =mat(1)
234 f1_1 =pm(59,mx)
235 f2_1 =pm(60,mx)
236 f4_1 =pm(61,mx)
237 f5_1 =pm(62,mx)
238 f11_1 =pm(63,mx)
239 f22_1 =pm(64,mx)
240 f44_1 =pm(65,mx)
241 f55_1 =pm(66,mx)
242 f12_1 =pm(67,mx)
243 f23_1 =pm(68,mx)
244
245 c11_1 =pm(69,mx)
246 c22_1 =pm(73,mx)
247 c33_1 =pm(74,mx)
248 c12_1 =pm(75,mx)
249 c23_1 =pm(76,mx)
250 c13_1 =pm(77,mx)
251
252 d11_1 =pm(40,mx)
253 d12_1 =pm(41,mx)
254 d13_1 =pm(42,mx)
255 d22_1 =pm(43,mx)
256 d23_1 =pm(44,mx)
257 d33_1 =pm(45,mx)
258 g12_1 =pm(46,mx)
259 g23_1 =pm(47,mx)
260 g31_1 =pm(48,mx)
261
262 DO i=1,nel
263 f1(i) = f1_1
264 f2(i) = f2_1
265 f4(i) = f4_1
266 f5(i) = f5_1
267 f11(i) = f11_1
268 f22(i) = f22_1
269 f44(i) = f44_1
270 f55(i) = f55_1
271 f12(i) = f12_1
272 f23(i) = f23_1
273
274 c11(i) = c11_1
275 c22(i) = c22_1
276 c33(i) = c33_1
277 c12(i) = c12_1
278 c23(i) = c23_1
279 c13(i) = c13_1
280
281 d11(i) = d11_1
282 d12(i) = d12_1
283 d13(i) = d13_1
284 d22(i) = d22_1
285 d23(i) = d23_1
286 d33(i) = d33_1
287 g12(i) = g12_1
288 g23(i) = g23_1
289 g31(i) = g31_1
290 ENDDO
291
292
293
294 DO i=1,nel
295 so1(i)=t1(i)
296 so2(i)=t2(i)
297 so3(i)=t3(i)
298 so4(i)=t4(i)
299 so5(i)=t5(i)
300 so6(i)=t6(i)
301 ENDDO
302
303 DO i=1,nel
304 t1(i)=t1(i)+d11(i)*e1(i)+d12(i)*e2(i)+d13(i)*e3(i)
305 t2(i)=t2(i)+d12(i)*e1(i)+d22(i)*e2(i)+d23(i)*e3(i)
306 t3(i)=t3(i)+d13(i)*e1(i)+d23(i)*e2(i)+d33(i)*e3(i)
307 t4(i)=t4(i)+g12(i)*e4(i)
308 t5(i)=t5(i)+g23(i)*e5(i)
309 t6(i)=t6(i)+g31(i)*e6(i)
310 ENDDO
311
312
313
314 IF (fib>0) THEN
315 mx = mat(1)
316 DO i=1,nel
317 efib(i) = pm(36,mx)
318 epsft(i)= pm(37,mx)
319 epsfc(i)= pm(38,mx)
320 epsf(i) = epsf(i)+ e1(i)
321 sigf(i) = efib(i)*epsf(i)
322 ENDDO
323
324
325
326 ENDIF
327
328
329
330 DO i=1,nel
331 IF (off(i)<em01) off(i)=zero
332 IF (off(i)<one) off(i)=off(i)*four_over_5
333 ENDDO
334
335
336
337 DO i=1,nel
338
339
340
341 wvec(i)=(one-dam(i,1))*sigt1(i)
342 IF (t1(i)>wvec(i)) THEN
343
344 IF (epc(i,1)==zero) THEN
345 epc(i,1) =
max(epe(i,1),zero)
346 ELSE
347 epc(i,1) =
max(epc(i,1)+e1(i),zero)
348 ENDIF
349
350 t1(i)=wvec(i)
351 t2(i)=t2(i)-d12(i)*e1(i)*dam(i,1)
352 t3(i)=t3(i)-d13(i)*e1(i)*dam(i,1)
353 IF (kd1(i)==0) THEN
354 kd1(i)=1
355
356
357 ENDIF
358 dam(i,1)=
min((dam(i,1)+delta(i)),one)
359 IF (dam(i,1)>=1. .AND. kd1(i)/=2) THEN
360 kd1(i)=2
361
362
363 ENDIF
364 ENDIF
365
366 IF( e1(i)<zero .AND. dam(i,1)>zero)
367 . epc(i,1)=
max(epc(i,1)+e1(i),zero)
368
369
370
371 wvec(i)=(one-dam(i,2))*sigt2(i)
372 IF (t2(i)>wvec(i)) THEN
373
374 IF (epc(i,2)==zero) THEN
375 epc(i,2)=
max(epe(i,2),zero)
376 ELSE
377 epc(i,2)=
max(epc(i,2)+e2(i),zero)
378 ENDIF
379
380 t1(i)=t1(i)-d12(i)*e2(i)*dam(i,2)
381 t2(i)=wvec(i)
382 t3(i)=t3(i)-d23(i)*e2(i)*dam(i,2)
383 IF (kd2(i)==0) THEN
384 kd2(i)=1
385
386
387 ENDIF
388 dam(i,2)=
min((dam(i,2)+delta(i)),one)
389 IF (dam(i,2)>=one.AND. kd2(i)/=2) THEN
390 kd2(i)=2
391
392
393 ENDIF
394 ENDIF
395
396 IF (e2(i)<zero .AND. dam(i,2)>zero)
397 . epc(i,2)=
max(epc(i,2)+e2(i),zero)
398
399
400
401 wvec(i)=(one-dam(i,3))*sigt3(i)
402 IF (t3(i)>wvec(i)) THEN
403
404 IF (epc(i,3)==zero) THEN
405 epc(i,3)=
max(epe(i,3),zero)
406 ELSE
407 epc(i,3)=
max(epc(i,3)+e3(i),zero)
408 ENDIF
409
410 t1(i)=t1(i)-d13(i)*e3(i)*dam(i,3)
411 t2(i)=t2(i)-d23(i)*e3(i)*dam(i,3)
412 t3(i)=wvec(i)
413 IF (kd3(i)==0) THEN
414 kd3(i)=1
415
416
417 ENDIF
418 dam(i,3)=
min((dam(i,3)+delta(i)),one)
419 IF (dam(i,3)>=1. .AND. kd3(i)/=2) THEN
420 kd3(i)=2
421
422
423 ENDIF
424 ENDIF
425
426 IF (e3(i)<zero .AND. dam(i,3)>zero)
427 . epc(i,3)=
max(epc(i,3)+e3(i),zero)
428
429 ENDDO
430
431
432
433 DO i=1,nel
434 IF (t1(i)<zero.AND.epc(i,1)>zero) THEN
435 t1(i)=zero
436 t2(i)=t2(i)-d12(i)*e1(i)*dam(i,1)
437 t3(i)=t3(i)-d13(i)*e1(i)*dam(i,1)
438 ENDIF
439 IF (t2(i)<zero.AND.epc(iTHEN
440 t1(i)=t1(i)-d12(i)*e2(i)*dam(i,2)
441 t2(i)=zero
442 t3(i)=t3(i)-d23(i)*e2(i)*dam(i,2)
443 ENDIF
444 IF (t3(i)<zero.AND.epc(i,3)>zero) THEN
445 t3(i)=zero
446 t1(i)=t1(i)-d13(i)*e3(i)*dam(i,3)
447 t2(i)=t2(i)-d23(i)*e3(i)*dam(i,3)
448 ENDIF
449 ENDDO
450
451
452
453 DO i=1,nel
454 wvec(i)=f1(i)*t1(i)+f2(i)*(t2(i)+t3(i))+
455 . f11(i)*t1(i)*t1(i)+f22(i)*(t2(i)*t2(i)+t3(i)*t3(i))+
456 . f55(i)*t5(i)*t5(i)+f44(i)*(t4(i)*t4(i)+t6(i)*t6(i))+
457 . two*f12(i)*(t1(i)*t2(i)+t1(i)*t3(i))+2*f23(i)*t2(i)*t3(i)
458 dam(i,4)=wvec(i)
459 tsaiwu(i)=
max(
min(wvec(i)/sigmy(i),one),tsaiwu(i))
460
461 ENDDO
462
463 DO i=1,nel
464 coef(i)=zero
465 IF (wvec(i)>sigmy(i).AND.off(i)==one) THEN
466 coef(i)=one
467 IF (kd4(i)==0) THEN
468 kd4(i)=1
469
470
471 ENDIF
472 ENDIF
473 ENDDO
474
475 DO i=1,nel
476 dp1(i)=f1(i)+two*f11(i)*so1(i)+two*f12(i)*(so2(i)+so3(i))
477 dp2(i)=f2(i)+two*f22(i)*so2(i)+two*f12(i)*so1(i)
478 . + so3(i)*f23(i)*two
479 dp3(i)=f2(i)+two*f22(i)*so3(i)+two*f12(i)*so1(i)
480 . + so2(i)*f23(i)*two
481 dp4(i)=two*f44(i)*so4(i)
482 dp5(i)=two*f55(i)*so5(i)
483 dp6(i)=two*f44(i)*so6(i)
484 ENDDO
485
486
487 3101 FORMAT(/' 205 WVEC SIGMY T1-T6 ',2e11.4/6e11.4)
488
489
490 3102 FORMAT(' F1 F2 F11 F22 F12 F23 F44 F55',3e11.4/5e11.4)
491
492 DO i=1,nel
493 ds1(i)=t1(i)-so1(i)
494 ds2(i)=t2(i)-so2(i)
495 ds3(i)=t3(i)-so3(i)
496 ds4(i)=t4(i)-so4(i)
497 ds5(i)=t5(i)-so5(i)
498 ds6(i)=t6(i)-so6(i)
499 ENDDO
500
501 DO i=1,nel
502 lamda(i)=(dp1(i)*ds1(i)+dp2(i)*ds2(i)+dp3(i)*ds3(i)
503 . +dp4(i)*ds4(i)+dp5(i)*ds5(i)+dp6(i)*ds6(i))*coef(i)
504 ENDDO
505
506
507 3103 FORMAT(' 207 LAMDA DS1-DS6 ',e11.4/6e11.4)
508
509 DO i=1,nel
510 cnn(i)=cn(i)-one
511 ENDDO
512
513 DO i=1,nel
514 plas(i)=one
515 IF (pla(i)>zero) plas(i)=pla(i)**cnn(i)
516 ENDDO
517
518 DO 208 i=1,nel
519 IF (lamda(i)==zero) GO TO 208
520 lamda(i)=lamda(i)*coef(i)/
521 . (dp1(i)*(d11(i)*dp1(i)+d12(i)*dp2(i)+d13(i)*dp3(i))+
522 . dp2(i)*(d12(i)*dp1(i)+d22(i)*dp2(i)+d23(i)*dp3(i))+
523 . dp3(i)*(d13(i)*dp1(i)+d23(i)*dp2(i)+d33(i)*dp3(i))+
524 . two*dp4(i)*g12(i)*dp4(i)+
525 . two*dp5(i)*g23(i)*dp5(i)+
526 . two*dp6(i)*g31(i)*dp6(i)+
527 . (so1(i)*dp1(i)+so2(i)*dp2(i)+so3(i)*dp3(i)+
528 . two*so4(i)*dp4(i)+2.*so5(i)*dp5(i)+2.*so6(i)*dp6(i))
529 . *cn(i)*cb(i)*plas(i))
530 208 CONTINUE
531
532 3104 FORMAT(' 208 LAMDA ',e11.4)
533
534 DO i=1,nel
535 dp1(i)=lamda(i)*dp1(i)
536 dp2(i)=lamda(i)*dp2(i)
537 dp3(i)=lamda(i)*dp3(i)
538 dp4(i)=lamda(i)*dp4(i)
539 dp5(i)=lamda(i)*dp5(i)
540 dp6(i)=lamda(i)*dp6(i)
541 ENDDO
542
543 3105 FORMAT(' 209 PLA DP1-DP6',e11.4/6e11.4)
544
545 DO i=1,nel
546 epe(i,1)=epe(i,1)-dp1(i)
547 epe(i,2)=epe(i,2)-dp2(i)
548 epe(i,3)=epe(i,3)-dp3(i)
549 ENDDO
550
551 DO i=1,nel
552 t1(i)=t1(i)-d11(i)*dp1(i)-d12(i)*dp2(i)-d13(i)*dp3(i)
553 t2(i)=t2(i)-d12(i)*dp1(i)-d22(i)*dp2(i)-d23(i)*dp3(i)
554 t3(i)=t3(i)-d13(i)*dp1(i)-d23(i)*dp2(i)-d33(i)*dp3(i)
555 t4(i)=t4(i)-g12(i)*dp4(i)*two
556 t5(i)=t5(i)-g23(i)*dp5(i)*two
557 t6(i)=t6(i)-g31(i)*dp6(i)*two
558 ENDDO
559
560 mx = mat(1)
561 wplaref = pm(98 ,mx)
562 DO i=1,nel
563 dwpla = half*
564 . (dp1(i)*(t1(i)+so1(i))+
565 . dp2(i)*(t2(i)+so2(i))+
566 . dp3(i)*(t3(i)+so3(i))+
567 . two*dp4(i)*(t4(i)+so4(i))+
568 . two*dp5(i)*(t5(i)+so5(i))+
569 . two*dp6(i)*(t6(i)+so6(i)))
570 dwpla =
max(dwpla ,zero) / wplaref
571 pla(i) = pla(i) + dwpla
572 pla(i)=
max(pla(i),zero)
573 ENDDO
574
575
576
577
578
579
580 CALL m14ftg(sig,ax ,ay ,az ,bx ,by ,
581 2 bz ,cx ,cy ,cz ,t1 ,t2 ,
582 3 t3 ,t4 ,t5 ,t6 ,nel)
583
584 DO i=1,nel
585 sig(i,1)=sig(i,1)*off(i)
586 sig(i,2)=sig(i,2)*off(i)
587 sig(i,3)=sig(i,3)*off(i)
588 sig(i,4)=sig(i,4)*off(i)
589 sig(i,5)=sig(i,5)*off(i)
590 sig(i,6)=sig(i,6)*off(i)
591 ENDDO
592
593
594
595 dt5=half*dt1
596 DO i=1,nel
597 eint(i)=eint(i)+dt5*vnew(i)*
598 . ( d1(i)*(sold1(i)+sig(i,1))
599 . + d2(i)*(sold2(i)+sig(i,2))
600 . + d3(i)*(sold3(i)+sig(i,3))
601 . + d4(i)*(sold4(i)+sig(i,4))
602 . + d5(i)*(sold5(i)+sig(i,5))
603 . + d6(i)*(sold6(i)+sig(i,6)))
604 eint(i)=eint(i)/vol(i)
605 dam(i,5)=kd1(i)*1000 + kd2(i)*100 + kd3(i)*10 + kd4(i) + 10000
606 ENDDO
607
608 DO i=1,nel
609 sigym =
max(em20,half*(one/f11(i)+one/f22(i)))
610 sigy(i)=sigmy(i)*sqrt(sigym)
611 defp(i)=pla(i)
612 ENDDO
613
614 1000 FORMAT(1x,'RUPTURE OF SOLID ELEMENT NUMBER ',i10)
615 6001 FORMAT(i6,' FAIL(14) E(',i5,') M(',a3,') T(',6e10.3,')')
616
617 RETURN
subroutine m14ama(pm, a, rx, ry, rz, sx, sy, sz, ax, ay, az, bx, by, bz, cx, cy, cz, nel, jcvt, jsph)
subroutine m14ftg(sig, ax, ay, az, bx, by, bz, cx, cy, cz, t1, t2, t3, t4, t5, t6, nel)
subroutine m14gtf(sig, ax, ay, az, bx, by, bz, cx, cy, cz, d1, d2, d3, d4, d5, d6, t1, t2, t3, t4, t5, t6, e1, e2, e3, e4, e5, e6, nel)