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) :: JCVT
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(MVSIZ),
83
85 . ax(mvsiz), ay(mvsiz), az(mvsiz), bx(mvsiz), by(mvsiz),
86 . bz(mvsiz), cx(mvsiz), cy(mvsiz),cz(mvsiz),
87 . wvec(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 . sigmy(mvsiz),
alpha(mvsiz), efib(mvsiz),
91 . epsft(mvsiz), epsfc(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(mvsiz),
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
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(i,2)>zero) THEN
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 DO i=1,nel
487 ds1(i)=t1(i)-so1(i)
488 ds2(i)=t2(i)-so2(i)
489 ds3(i)=t3(i)-so3(i)
490 ds4(i)=t4(i)-so4(i)
491 ds5(i)=t5(i)-so5(i)
492 ds6(i)=t6(i)-so6(i)
493 ENDDO
494
495 DO i=1,nel
496 lamda(i)=(dp1(i)*ds1(i)+dp2(i)*ds2(i)+dp3(i)*ds3(i)
497 . +dp4(i)*ds4(i)+dp5(i)*ds5(i)+dp6(i)*ds6(i))*coef(i)
498 ENDDO
499
500 DO i=1,nel
501 cnn(i)=cn(i)-one
502 ENDDO
503
504 DO i=1,nel
505 plas(i)=one
506 IF (pla(i)>zero) plas(i)=pla(i)**cnn(i)
507 ENDDO
508
509 DO 208 i=1,nel
510 IF (lamda(i)==zero) GO TO 208
511 lamda(i)=lamda(i)*coef(i)/
512 . (dp1(i)*(d11(i)*dp1(i)+d12(i)*dp2(i)+d13(i)*dp3(i))+
513 . dp2(i)*(d12(i)*dp1(i)+d22(i)*dp2(i)+d23(i)*dp3(i))+
514 . dp3(i)*(d13(i)*dp1(i)+d23(i)*dp2(i)+d33(i)*dp3(i))+
515 . two*dp4(i)*g12(i)*dp4(i)+
516 . two*dp5(i)*g23(i)*dp5(i)+
517 . two*dp6(i)*g31(i)*dp6(i)+
518 . (so1(i)*dp1(i)+so2(i)*dp2(i)+so3(i)*dp3(i)+
519 . two*so4(i)*dp4(i)+2.*so5(i)*dp5(i)+2.*so6(i)*dp6(i))
520 . *cn(i)*cb(i)*plas(i))
521 208 CONTINUE
522
523 3104 FORMAT(' 208 LAMDA ',e11.4)
524
525 DO i=1,nel
526 dp1(i)=lamda(i)*dp1(i)
527 dp2(i)=lamda(i)*dp2(i)
528 dp3(i)=lamda(i)*dp3(i)
529 dp4(i)=lamda(i)*dp4(i)
530 dp5(i)=lamda(i)*dp5(i)
531 dp6(i)=lamda(i)*dp6(i)
532 ENDDO
533
534 DO i=1,nel
535 epe(i,1)=epe(i,1)-dp1(i)
536 epe(i,2)=epe(i,2)-dp2(i)
537 epe(i,3)=epe(i,3)-dp3(i)
538 ENDDO
539
540 DO i=1,nel
541 t1(i)=t1(i)-d11(i)*dp1(i)-d12(i)*dp2(i)-d13(i)*dp3(i)
542 t2(i)=t2(i)-d12(i)*dp1(i)-d22(i)*dp2(i)-d23(i)*dp3(i)
543 t3(i)=t3(i)-d13(i)*dp1(i)-d23(i)*dp2(i)-d33(i)*dp3(i)
544 t4(i)=t4(i)-g12(i)*dp4(i)*two
545 t5(i)=t5(i)-g23(i)*dp5(i)*two
546 t6(i)=t6(i)-g31(i)*dp6(i)*two
547 ENDDO
548
549 mx = mat(1)
550 wplaref = pm(98 ,mx)
551 DO i=1,nel
552 dwpla = half*
553 . (dp1(i)*(t1(i)+so1(i))+
554 . dp2(i)*(t2(i)+so2(i))+
555 . dp3(i)*(t3(i)+so3(i))+
556 . two*dp4(i)*(t4(i)+so4(i))+
557 . two*dp5(i)*(t5(i)+so5(i))+
558 . two*dp6(i)*(t6(i)+so6(i)))
559 dwpla =
max(dwpla ,zero) / wplaref
560 pla(i) = pla(i) + dwpla
561 pla(i)=
max(pla(i),zero)
562 ENDDO
563
564
565
566
567
568
569 CALL m14ftg(sig,ax ,ay ,az ,bx ,by ,
570 2 bz ,cx ,cy ,cz ,t1 ,t2 ,
571 3 t3 ,t4 ,t5 ,t6 ,nel)
572
573 DO i=1,nel
574 sig(i,1)=sig(i,1)*off(i)
575 sig(i,2)=sig(i,2)*off(i)
576 sig(i,3)=sig(i,3)*off(i)
577 sig(i,4)=sig(i,4)*off(i)
578 sig(i,5)=sig(i,5)*off(i)
579 sig(i,6)=sig(i,6)*off(i)
580 ENDDO
581
582
583
584 dt5=half*dt1
585 DO i=1,nel
586 eint(i)=eint(i)+dt5*vnew(i)*
587 . ( d1(i)*(sold1(i)+sig(i,1))
588 . + d2(i)*(sold2(i)+sig(i,2))
589 . + d3(i)*(sold3(i)+sig(i,3))
590 . + d4(i)*(sold4(i)+sig(i,4))
591 . + d5(i)*(sold5(i)+sig(i,5))
592 . + d6(i)*(sold6(i)+sig(i,6)))
593 eint(i)=eint(i)/vol(i)
594 dam(i,5)=kd1(i)*1000 + kd2(i)*100 + kd3(i)*10 + kd4(i) + 10000
595 ENDDO
596
597 DO i=1,nel
598 sigym =
max(em20,half*(one/f11(i)+one/f22(i)))
599 sigy(i)=sigmy(i)*sqrt(sigym)
600 defp(i)=pla(i)
601 ENDDO
602
603 1000 FORMAT(1x,'RUPTURE OF SOLID ELEMENT NUMBER ',i10)
604
605 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)