37 1 NEL ,NUPARAM,NUVAR ,MFUNC ,KFUNC ,NPF ,
38 2 TF ,TIME ,TIMESTEP,UPARAM ,RHO0 ,RHO ,
39 3 VOLUME ,EINT ,L_DMG ,DMG ,
40 4 EPSPXX ,EPSPYY ,EPSPZZ ,EPSPXY ,EPSPYZ ,EPSPZX ,
41 5 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
42 6 EPSXX ,EPSYY ,EPSZZ ,EPSXY ,EPSYZ ,EPSZX ,
43 7 SIGOXX ,SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
44 8 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
45 9 SIGVXX ,SIGVYY ,SIGVZZ ,SIGVXY ,SIGVYZ ,SIGVZX ,
46 A SOUNDSP,VISCMAX,UVAR ,OFF ,NGL ,IPT ,
47 B IPM ,MAT ,EPSP ,IPLA ,SIGY ,DEFP ,
55#include "implicit_f.inc"
104#include "param_c.inc"
105#include "scr17_c.inc"
106#include "units_c.inc"
107#include "com08_c.inc"
109 INTEGER NEL, NUPARAM, NUVAR,IPT,
110 . NGL(NEL),MAT(NEL),IPM(NPROPMI,*),IPLA, ITER, IFLAG
112 . TIME,TIMESTEP,UPARAM(*),
113 . RHO(NEL),RHO0(NEL),VOLUME(NEL),EINT(NEL),
114 . EPSPXX(NEL),EPSPYY(NEL),EPSPZZ(NEL),
115 . EPSPXY(NEL),EPSPYZ(NEL),EPSPZX(NEL),
116 . (NEL),DEPSYY(NEL),DEPSZZ(NEL),
117 . DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
118 . EPSXX(NEL) ,EPSYY(NEL) ,EPSZZ(NEL) ,
119 . EPSXY(NEL) ,EPSYZ(NEL) ,EPSZX(NEL) ,
120 . SIGOXX(NEL),SIGOYY(NEL),SIGOZZ(NEL),
121 . sigoxy(nel),sigoyz(nel),sigozx(nel),
125 INTEGER,
INTENT(IN) :: L_DMG
130 . SIGNXX(NEL),SIGNYY(NEL),SIGNZZ(NEL),
131 . SIGNXY(NEL),SIGNYZ(NEL),SIGNZX(NEL),
132 . SIGVXX(NEL),SIGVYY(NEL),SIGVZZ(NEL),
133 . SIGVXY(NEL),SIGVYZ(NEL),SIGVZX(NEL),
134 . SOUNDSP(NEL),VISCMAX(NEL),SIGY(NEL),DEFP(NEL)
139 . uvar(nel,nuvar), off(nel)
140 my_real,
DIMENSION(NEL,L_DMG),
INTENT(INOUT) :: dmg
144 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
152 INTEGER I,J,IADBUFV,J1,J2,JJ(MVSIZ),NFUNC,
153 . NRATE(MVSIZ),IPOS1(MVSIZ),IPOS2(MVSIZ),IAD1(MVSIZ),
154 . ILEN1(MVSIZ),IAD2(MVSIZ),ILEN2(MVSIZ),
155 . NFUNCV(MVSIZ),PFUN(MVSIZ),
156 . iposp(mvsiz),iadp(mvsiz),ilenp(mvsiz),nfuncm,nratem,
157 . ipflag,iparam,npar,iyeild_tab,iadbuf,
158 , itab(100),ntable,nxh,ndim,ipos,nxk
160 . f(mvsiz), sigm(mvsiz), epsm(mvsiz), epsp1(mvsiz),
161 . lamda(mvsiz), fstar(mvsiz),p0(mvsiz),pn(mvsiz),
163 . visp(mvsiz), q1, q2, q3,
166 . sigxx(mvsiz),sigyy(mvsiz),sigxy(mvsiz),
167 . sigzx(mvsiz),sigyz(mvsiz), sigzz(mvsiz),
168 . sganxx(mvsiz), sganyy(mvsiz), sganzz(mvsiz),
169 . sganxy(mvsiz), sganxz(mvsiz), sganyz(mvsiz)
171 . e,nu,p,dav,vm,r,fac,epst,ep1,ep2,ep3,ep4,ep5,ep6,
172 . e1,e2,e3,e4,e5,e6,c,cc,d,y,yp,e42,e52,e62,epst2,
173 . dp11,dp22,dp33,dp12,dp13,dp23,a1,a2,
174 . d11,d22,d33,d12,d13,d23,dcrf,dcrm, dsepp,dcd,lam1,
175 . c1(mvsiz),c2(mvsiz) ,g(mvsiz),df,coh,sih,va,crit,conv,
176 . fg(mvsiz),fn1(mvsiz),yld,c11,p1, var, dtinv, sqr22,
177 . sxx,syy,szz, sigm1, sigm2, va1, vm1, va11,a21,yp1,yp2,
178 . df1, df2, xfac,yfac,dx2,xx(2),dft,yy
182 ntable = ipm(226,mat(1))
186 itab(1)=ipm(226+1,mat(1))
187 iadbuf = ipm(7,mat(1))-1
188 xfac = uparam(iadbuf + 22 )
189 yfac = uparam(iadbuf + 23 )
192 iadbufv = ipm(7,mat(1))-1
193 e = uparam(iadbufv+1)
194 nu = uparam(iadbufv+2)
195 iflag = uparam(iadbufv+18)
196 et = uparam(iadbufv+4)
197 n = uparam(iadbufv+5)
198 csd = one/uparam(iadbufv+6)
199 q1 = uparam(iadbufv+8)
200 q2 = uparam(iadbufv+9)
201 q3 = uparam(iadbufv+10)
202 sn = uparam(iadbufv+11)
203 IF(sn==zero)sn = ep20
204 epsn = uparam(iadbufv+12)
205 fi = uparam(iadbufv+13)
206 fn = uparam(iadbufv+14)
207 fc = uparam(iadbufv+15)
208 ff = uparam(iadbufv+16)
209 fu = uparam(iadbufv+17)
212 yeild0(i)= uparam(iadbufv+3)
213 visp(i) = uparam(iadbufv+7)
215 g(i) = half*e/(one + nu)
216 c11 = e/three/(one- two*nu)
217 soundsp(i) = sqrt((c11 + four_over_3*g(i))/rho0(i))
219 c1(i) = e*(one-nu) /((one + nu)*(one - two*nu))
220 c2(i) = c1(i)*nu/(one - nu)
227 IF(n<one)epsm(i)=em20
235 IF(iyeild_tab > 0)
THEN
241 uvar(i,2) = yeild0(i)
247 sqr22= one/sqrt(two*pi)
260 signxx(i)=sigoxx(i) + c1(i)*depsxx(i)
261 . + c2(i)*(depsyy(i) + depszz(i))
262 signyy(i)=sigoyy(i) + c1(i)*depsyy(i)
263 . + c2(i)*(depsxx(i) + depszz(i))
264 signzz(i)=sigozz(i) + c1(i)*depszz(i)
265 . + c2(i)*(depsxx(i) + depsyy(i))
266 signxy(i)=sigoxy(i) + g(i)*depsxy(i)
267 signyz(i)=sigoyz(i) + g(i)*depsyz(i)
268 signzx(i)=sigozx(i) + g(i)*depszx(i)
269 IF(off(i)==zero)
THEN
279 dtinv=timestep/
max(timestep**2,em20)
294 df = (fu- fc)/(ff-fc)
295 fstar(i)= fc + df*(f(i)-fc)
298 pn(i)= (signxx(i)+signyy(i)+signzz(i))*third
299 sxx = signxx(i)-pn(i)
300 syy = signyy(i)-pn(i)
301 szz = signzz(i)-pn(i)
302 vm = half*(sxx**2 + syy**2 + szz**2)
303 . + (signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
305 sigm1 = one/
max(sigm(i),em20)
306 var = three_half * q2 * pn(i)* sigm1
308 coh = half * (var + one/
max(em20,var))
309 sih = half * (var - one/
max(em20,var))
310 va = one + q3 * fstar(i)**2 - two * q1 * fstar(i) * coh
315 IF(vm < yld .OR. sigy(i) == zero )
THEN
334 a21 = sigm1 /(one - f(i))
336 a1 = fn*exp(-half*((epsm(i)-epsn)/sn)**2)/sn
339 vm1 = one/
max(vm,em20)
340 va1 = one/
max(va,em20)
341 va11 =half*q1*q2*fstar(i)*sih*va1
342 d11 = half*(two*signxx(i) - signyy(i)-signzz(i))*vm1
344 d22 = half*(two*signyy(i) - signxx(i)-signzz
346 d33 = half*(two*signzz(i) - signxx(i)-signyy(i))*vm1
348 d12 = three*signxy(i)*vm1
349 d13 = three*signzx(i)*vm1
353 . two*(d12*signxy(i) + d13*signzx(i) + d23*signyz(i))
356 dcrf=-sigm(i)*(q3*fstar(i
357 dcrm= - va - three*va11*pn(i)*sigm1
359 dsepp = et*(one+ ( epsp(i)*csd)**visp(i))
361 dsepp = et * n * epsm(i) ** (n - 1)*
362 . (one + ( epsp(i)*csd)**visp(i))
365 dcd = c1(i)*(d11**2 + d22**2 + d33**2) +
366 . two*c2(i)*(d11*d22 + d11*d33 + d22*d33) +
369 lam1= dcd - dcrm* dsepp*a2 -
370 . dcrf*((one - f(i))*(d11+d22+d33) + a1*a2)
371 IF(lam1/=zero) lamda(i) =
max(zero,(vm - yld))/lam1
373 dp11 = dp11 + lamda(i) * d11
374 dp22 = dp22 + lamda(i) * d22
375 dp33 = dp33 + lamda(i) * d33
376 dp12 = dp12 + lamda(i) * d12
377 dp13 = dp13 + lamda(i) * d13
378 dp23 = dp23 + lamda(i) * d23
380 signxx(i)= sigxx(i) - c1(i)*dp11 - c2(i)*(dp22 + dp33)
381 signyy(i)= sigyy(i) - c1(i)*dp22 - c2(i)*(dp11 + dp33)
382 signzz(i)= sigzz(i) - c1(i)*dp33 - c2(i)*(dp22 + dp11)
383 signxy(i)= sigxy(i) - two*g(i)*dp12
384 signzx(i)= sigzx(i) - two*g(i)*dp13
385 signyz(i)= sigyz(i) - two*g(i)*dp23
388 epsp1(i)= signxx(i)*dp11 + signyy(i)*dp22 + signzz(i)*dp33
389 . + two*(signxy(i)*dp12 + signzx(i)*dp13 + signyz(i)*dp23)
390 epsp1(i)=epsp1(i)*a21
391 epsp1(i) =
max(zero, epsp1(i))
392 pn(i)= (signxx(i)+signyy(i)+signzz(i))*third
393 sxx = signxx(i)-pn(i)
394 syy = signyy(i)-pn(i)
395 szz = signzz(i)-pn(i)
396 vm = half*(sxx**2 + syy**2 + szz**2)
397 . + (signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
399 var = three_half * q2*pn(i) *sigm1
401 coh = half * ( var + one/
max(em20,var))
402 sih = half * ( var - one/
max(em20,var))
403 va = one + q3*fstar(i)**2 - two*q1*fstar(i)*coh
404 va = sqrt(
max(zero,va))
409 fg(i) = fg(i) + (one-f(i))*(dp11+dp22+dp33)
411 fn1(i) = fn1(i) + a1*epsp1(i)
412 epsm(i) = epsm(i) + epsp1(i)
413 f(i) = fi + fg(i) + fn1(i)
414 IF(q1==zero.and.q2==zero.and.q3==zero)f(i)=zero
415 IF(f(i)<zero) f(i)=zero
416 sigm(i) = (yeild0(i) + et * epsm(i) ** n)*
417 . (one+ (epsp(i)*csd) ** visp(i))
425 fstar(i) = fc+ (fu-fc)/(ff-fc)*(f(i)-fc)
436 IF(off(i)==one.AND.fstar(i)>=ff)
THEN
438 WRITE(iout, 1000) ngl(i)
439 WRITE(istdo,1100) ngl(i),tt
442 IF(off(i)<em01)off(i)=zero
443 off(i)=off(i)*four_over_5
457 fstar(i) = fc + (fu - fc)/(ff-fc)*(f(i)-fc)
458 df = (fu - fc)/(ff - fc)
459 fstar(i) = fc + df*(f(i)-fc)
461 pn(i)= (signxx(i)+signyy(i)+signzz(i))*third
462 sxx = signxx(i)-pn(i)
463 syy = signyy(i)-pn(i)
464 szz = signzz(i)-pn(i)
465 vm = half*(sxx**2 + syy**2 + szz**2)
466 . + (signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
468 sigm1 = one/
max(em20,sigm(i))
470 var = three_half * q2 * pn(i) * sigm1
472 coh = half*(var + one/
max(em20,var))
473 sih = half*(var - one/
max(em20,var))
475 va = one +q3*fstar(i)**2 - two*q1*fstar(i)
477 va= one +q3*fstar(i)**2 - two*q1*fstar(i)*coh
481 crit = vm * sigm2 - va
483 sigy(i) = yld*sqrt(
max(zero, va))
484 IF(crit < zero .OR. sigy(i) == zero)
THEN
499 a21 = sigm1 /(one - f(i))
501 a1 = fn*exp(-half*((epsm(i)-epsn)/sn)**2)/sn
511 d11 = (two*signxx(i) - signyy(i) - signzz(i))*sigm2
512 d22 = (two*signyy(i) - signxx(i) - signzz(i))*sigm2
513 d33 = (two*signzz(i) - signxx(i) - signyy(i))*sigm2
514 d12 = six*signxy(i)*sigm2
515 d13 = six*signzx(i)*sigm2
516 d23 = six*signyz(i)*sigm2
517 dcrm = -two * vm * sigm2 * sigm1
518 dcrf = two * q1 * df - two * q3 * df * fstar(i)
520 va1 = q1 * q2 * fstar(i) * sih*sigm1
521 d11 = (two*signxx(i) - signyy(i) - signzz(i))*sigm2 + va1
522 d22 = (two*signyy(i) - signxx(i) - signzz(i))*sigm2 + va1
523 d33 = (two*signzz(i) - signxx(i) - signyy(i))*sigm2 + va1
524 d12 = six * signxy(i)*sigm2
525 d13 = six * signzx(i)*sigm2
526 d23 = six * signyz(i)*sigm2
527 dcrm =-two* vm*sigm1*sigm2 - three * va1 * pn(i)*sigm1
528 dcrf = two * q1*df*coh - two*q3*df*fstar(i)
531 a2 = d11 * signxx(i) + d22 * signyy(i) + d33 * signzz(i) +
532 . two*(d12 * signxy(i) + d13 * signzx(i) + d23 * signyz(i) )
536 dsepp = et*(one+ ( epsp(i)*csd)**visp(i))
538 dsepp = et * n * epsm(i) ** (n - 1)*
539 . (one+ ( epsp(i)*csd)**visp(i))
541 dcd = c1(i) * (d11**2 + d22**2 + d33**2) +
542 . two*c2(i)*(d11*d22 + d11*d33 + d22*d33) +
543 . two*g(i) *(d12**2 + d13**2 + d23**2)
544 lam1 = dcd - dcrm * dsepp * a2 -
545 . dcrf * ((one - f(i)) * (d11 + d22 + d33) + a1*a2)
546 IF(lam1/=zero) lamda(i) =
max(zero,(vm * sigm2 - va))/lam1
550 dp11 = dp11 + lamda(i) * d11
551 dp22 = dp22 + lamda(i) * d22
552 dp33 = dp33 + lamda(i) * d33
553 dp12 = dp12 + lamda(i) * d12
554 dp13 = dp13 + lamda(i) * d13
555 dp23 = dp23 + lamda(i) * d23
557 signxx(i)= sigxx(i) - c1(i) * dp11 - c2(i) * (dp22 + dp33)
558 signyy(i)= sigyy(i) - c1(i) * dp22 - c2(i) * (dp11 + dp33)
559 signzz(i)= sigzz(i) - c1(i) * dp33 - c2(i) * (dp22 + dp11)
560 signxy(i)= sigxy(i) - two*g(i
561 signzx(i)= sigzx(i) - two*g(i) * dp13
562 signyz(i)= sigyz(i) - two*g(i) * dp23
566 epsp1(i) = signxx(i)*dp11 + signyy
567 . two*(signxy(i)*dp12 + signzx(i)*dp13 + signyz(i)*dp23)
568 epsp1(i) = epsp1(i)*a21
569 epsp1(i) =
max(zero, epsp1(i))
570 pn(i)= (signxx(i)+signyy(i)+signzz(i))*third
571 sxx = signxx(i)-pn(i)
572 syy = signyy(i)-pn(i)
573 szz = signzz(i)-pn(i)
574 vm = half*(sxx**2 + syy**2 + szz**2)
575 . + (signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
577 var = three_half * q2 * pn(i)*sigm1
579 coh = half * (var + one/
max(em20,var))
580 sih = half * (var - one/
max(em20,var))
582 va = one + q3*fstar(i)**2 - two*q1*fstar(i)
584 va = one + q3*fstar(i)**2 - two*q1*fstar(i)*coh
589 fg(i) = fg(i) + (one-f(i))*(dp11 + dp22 +dp33)
590 sigy(i) = sigm(i)*sqrt(
max(zero,va))
591 fn1(i) = fn1(i) + a1 * epsp1(i)
592 epsm(i) = epsm(i) + epsp1(i)
593 f(i) = fi + fg(i) + fn1(i)
594 IF(f(i)<zero) f(i)=zero
595 sigm(i)= (yeild0(i) + et * epsm(i)**n)*(one +
596 . (epsp(i)*csd)**visp(i))
604 fstar(i) = fc+ (fu-fc)*(f(i)-fc)/(ff-fc)
617 IF(off(i)==one.AND.fstar(i)>=ff)
THEN
619 WRITE(iout, 1000) ngl(i)
620 WRITE(istdo,1100) ngl(i),tt
623 IF(off(i)<em01)off(i)=zero
624 off(i)=off(i)*four_over_5
627 ELSEIF(iflag==2)
THEN
638 df = (fu - fc)/(ff - fc)
639 fstar(i) = fc + df*(f(i)-fc)
642 pn(i)= (signxx(i)+signyy(i)+signzz(i))*third
643 sxx = signxx(i)-pn(i)
644 syy = signyy(i)-pn(i)
645 szz = signzz(i)-pn(i)
646 vm = half*(sxx**2 + syy**2 + szz**2)
647 . + (signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
650 sigm1 = one/
max(em20,sigm(i))
652 var = half * q2 * pn(i) * sigm1
654 coh = half*(var + one/
max(em20,var))
655 sih = half*(var - one/
max(em20,var))
657 va = one +q3*fstar(i)**2 - two*q1*fstar(i)
659 va= one+q3*fstar(i)**2 - two*q1*fstar(i)*coh
662 crit = vm * sigm2 - va
663 sigy(i) = yld*sqrt(
max(zero, va))
664 IF(crit < zero .OR. sigy(i) == zero)
THEN
678 a21 = sigm1 /(one - f(i))
680 a1 = fn*exp(-half*((epsm(i)-epsn)/sn)**2)/sn
692 d11 = (two*signxx(i) - signyy(i) - signzz(i))*sigm2
693 d22 = (two*signyy(i) - signxx(i) - signzz(i))*sigm2
694 d33 = (two*signzz(i) - signxx(i) - signyy(i))*sigm2
695 d12 = six*signxy(i)*sigm2
696 d13 = six*signzx(i)*sigm2
697 d23 = six*signyz(i)*sigm2
698 dcrm = -two * vm * sigm2 * sigm1
699 dcrf = two * q1 * df - two * q3 * df * fstar(i)
701 va1 = q1 * q2 * fstar(i) * sih*sigm1
702 d11 = (two*signxx(i) - signyy(i) - signzz(i))*sigm2 + va1
703 d22 = (two*signyy(i) - signxx(i) - signzz(i))*sigm2 + va1
704 d33 = (two*signzz(i) - signxx(i) - signyy(i))*sigm2 + va1
705 d12 = six * signxy(i)*sigm2
706 d13 = six * signzx(i)*sigm2
707 d23 = six * signyz(i)*sigm2
708 dcrm =-two * vm*sigm1*sigm2 - three * va1 * pn(i)*sigm1
709 dcrf = two * q1*df*coh - two*q3*df*fstar(i)
712 a2 = d11 * signxx(i) + d22 * signyy(i) + d33 * signzz(i) +
713 . two*(d12 * signxy(i) + d13 * signzx(i) + d23 * signyz(i) )
717 dsepp = et*(one + ( epsp(i)*csd)**visp(i))
719 dsepp = et * n * epsm(i) ** (n - one)*
720 . (one + ( epsp(i)*csd)**visp(i))
722 dcd = c1(i) * (d11**2 + d22**2 + d33**2) +
723 . two*c2(i)*(d11*d22 + d11*d33 + d22*d33) +
724 . two*g(i) *(d12**2 + d13**2 + d23**2)
725 lam1 = dcd - dcrm * dsepp * a2 -
726 . dcrf * ((one - f(i)) * (d11 + d22 + d33) + a1*a2)
727 IF(lam1/=zero) lamda(i) =
max(zero,(vm * sigm2 - va))/lam1
731 dp11 = dp11 + lamda(i) * d11
732 dp22 = dp22 + lamda(i) * d22
733 dp33 = dp33 + lamda(i) * d33
734 dp12 = dp12 + lamda(i) * d12
735 dp13 = dp13 + lamda(i) * d13
736 dp23 = dp23 + lamda(i) * d23
738 signxx(i)= sigxx(i) - c1(i) * dp11 - c2(i) * (dp22 + dp33)
739 signyy(i)= sigyy(i) - c1(i) * dp22 - c2(i) * (dp11 + dp33)
740 signzz(i)= sigzz(i) - c1(i) * dp33 - c2(i) * (dp22 + dp11)
741 signxy(i)= sigxy(i) - two*g(i) * dp12
742 signzx(i)= sigzx(i) - two*g(i) * dp13
743 signyz(i)= sigyz(i) - two*g(i) * dp23
747 epsp1(i) = signxx(i)*dp11 + signyy(i)*dp22 + signzz(i)*dp33+
748 . two*(signxy(i)*dp12 + signzx(i)*dp13 + signyz(i)*dp23)
749 epsp1(i) = epsp1(i)*a21
750 epsp1(i) =
max(zero, epsp1(i))
752 pn(i)= (signxx(i)+signyy(i)+signzz(i))*third
753 sxx = signxx(i)-pn(i)
754 syy = signyy(i)-pn(i)
755 szz = signzz(i)-pn(i)
756 vm = half*(sxx**2 + syy**2 + szz**2)
757 . + (signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
759 var = three_half * q2 * pn(i)*sigm1
761 coh = half * (var + one/
max(em20,var))
762 sih = half * (var - one/
max(em20,var))
764 va = one + q3*fstar(i)**2 - two*q1*fstar(i)
766 va = one + q3*fstar(i)**2 - two*q1*fstar(i)*coh
771 fg(i) = fg(i) + (one -f(i))*(dp11 + dp22 +dp33)
772 sigy(i) = sigm(i)*sqrt(
max(zero,va))
774 pn(i)= (signxx(i)+signyy(i)+signzz(i))*third
775 IF(pn(i)>=zero) fn1(i)= fn1(i) + a1 * epsp1(i)
776 epsm(i) = epsm(i) + epsp1(i)
777 f(i) = fi + fg(i) + fn1(i)
778 IF(f(i)<zero) f(i)=zero
779 sigm(i)= (yeild0(i) + et * epsm(i)**n)*(one+
780 . (epsp(i)*csd)**visp(i))
788 fstar(i) = fc+ (fu-fc)*(f(i)-fc)/(ff-fc)
802 IF(off(i)==one.AND.fstar(i)>=ff)
THEN
804 WRITE(iout, 1000) ngl(i)
805 WRITE(istdo,1100) ngl(i),tt
809 IF(off(i)<em01)off(i)=zero
810 off(i)=off(i)*four_over_5
826 df = (fu- fc)/(ff-fc)
827 fstar(i)= fc + df*(f(i)-fc)
830 pn(i)= (signxx(i)+signyy(i)+signzz(i))*third
831 sxx = signxx(i)-pn(i)
832 syy = signyy(i)-pn(i)
833 szz = signzz(i)-pn(i)
834 vm = half*(sxx**2 + syy**2 + szz**2)
835 . + (signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
837 sigm1 = one/
max(sigm(i),em20)
838 var = three_half * q2 * pn(i)* sigm1
840 coh = half * (var + one/
max(em20,var))
841 sih = half * (var - one/
max(em20,var))
842 va = one + q3 * fstar(i)**2 - two *q1*fstar(i)*coh
843 va= sqrt(
max(zero,va))
848 IF(vm < yld .OR. sigy(i) == zero)
THEN
869 a21 = sigm1 /(one - f(i))
871 a1 = fn*exp(-half*((epsm(i)-epsn)/sn)**2)/sn
875 vm1 = one/
max(vm,em20)
876 va1 = one/
max(va,em20)
877 va11 =half*q1*q2*fstar(i)*sih*va1
878 d11 = half*(two*signxx(i)
879 . - signyy(i)-signzz(i))*vm1 + va11
880 d22 = half*(two*signyy(i)
881 . - signxx(i)-signzz(i))*vm1 + va11
882 d33 = half*(two*signzz(i)
883 . - signxx(i)-signyy(i))*vm1 + va11
884 d12 = three*signxy(i)*vm1
885 d13 = three*signzx(i)*vm1
886 d23 = three*signyz(i)*vm1
888 a2 = d11*signxx(i) + d22*signyy(i) + d33*signzz(i) +
889 . two*(d12*signxy(i) + d13*signzx(i) + d23*signyz(i))
892 dcrf=-sigm(i)*(q3*fstar(i)*df-q1*coh*df)*va1
893 dcrm= - va - three*va11*pn(i)*sigm1
895 dsepp = et*(one + ( epsp(i)*csd)**visp(i))
897 dsepp = et * n * epsm(i) ** (n - one)*
898 . (one + ( epsp(i)*csd)**visp(i))
900 dcd = c1(i)*(d11**2 + d22**2 + d33**2) +
901 . two*c2(i)*(d11*d22 + d11*d33 + d22*d33) +
902 . two*g(i)*(d12**2 + d13**2 + d23**2)
904 lam1= dcd - dcrm* dsepp*a2 -
905 . dcrf*((one-f(i))*(d11+d22+d33) + a1*a2)
906 IF(lam1/=zero) lamda(i) =
max(zero,(vm - yld))/lam1
908 dp11 = dp11 + lamda(i) * d11
909 dp22 = dp22 + lamda(i) * d22
910 dp33 = dp33 + lamda(i) * d33
911 dp12 = dp12 + lamda(i) * d12
912 dp13 = dp13 + lamda(i) * d13
913 dp23 = dp23 + lamda(i) * d23
915 signxx(i)= sigxx(i) - c1(i)*dp11 - c2(i)*(dp22 + dp33)
916 signyy(i)= sigyy(i) - c1(i)*dp22 - c2(i)*(dp11 + dp33)
917 signzz(i)= sigzz(i) - c1(i)*dp33 - c2(i)*(dp22 + dp11)
918 signxy(i)= sigxy(i) - two*g(i)*dp12
919 signzx(i)= sigzx(i) - two*g(i)*dp13
920 signyz(i)= sigyz(i) - two*g(i)*dp23
922 epsp1(i)= signxx(i)*dp11 + signyy(i)*dp22 + signzz(i)*dp33
923 . + two*(signxy(i)*dp12 + signzx(i)*dp13 + signyz(i)*dp23)
924 epsp1(i)=epsp1(i)*a21
925 epsp1(i) =
max(zero, epsp1(i))
926 pn(i)= (signxx(i)+signyy(i)+signzz(i))*third
927 sxx = signxx(i)-pn(i)
928 syy = signyy(i)-pn(i)
929 szz = signzz(i)-pn(i)
930 vm = half*(sxx**2 + syy**2 + szz**2)
931 . + (signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
933 var = three_half * q2*pn(i) *sigm1
935 coh = half * ( var + one/
max(em20,var))
936 sih = half * ( var - one/
max(em20,var))
937 va = one + q3*fstar(i)**2 - two*q1*fstar(i)*coh
938 va = sqrt(
max(zero,va))
942 fg(i) = fg(i) + (one-f(i))*(dp11+dp22+dp33)
943 pn(i) = (sigxx(i) + sigyy(i) + sigzz(i))*third
946 IF(pn(i)>=zero) fn1(i) = fn1(i) + a1*epsp1(i)
947 epsm(i) = epsm(i) + epsp1(i)
948 f(i) = fi + fg(i) + fn1(i)
950 . q2==zero.and.q3==zero) f(i)= zero
951 IF(f(i)<zero) f(i)=zero
952 sigm(i) = (yeild0(i) + et * epsm(i) ** n)*
953 . (one + (epsp(i)*csd) ** visp(i))
961 fstar(i) = fc+ (fu-fc)/(ff-fc)*(f(i)-fc)
972 IF(off(i)==one.AND.fstar(i)>=ff)
THEN
974 WRITE(iout, 1000) ngl(i)
975 WRITE(istdo,1100) ngl(i),tt
978 IF(off(i)<em01)off(i)=zero
979 off(i)=off(i)*four_over_5
988 IF(iyeild_tab > 0)
THEN
989 iadbuf = ipm(7,mat(1))-1
990 xfac = uparam(iadbuf + 22 )
991 yfac = uparam(iadbuf + 23 )
993 ndim= table(itab(1))%NDIM
995 nxk=
SIZE(table(itab(1))%X(2)%VALUES)
997 dx2 = table(itab(1))%X(2)%VALUES(j)*xfac- epsp(i)
998 IF(dx2 >= zero .OR. j == nxk)
THEN
1000 r=(table(itab(1))%X(2)%VALUES(j)*xfac-epsp(i))/
1001 . (table(itab(1))%X(2)%VALUES(j)*xfac
1002 . -table(itab(1))%X(2)%VALUES(j-1)*xfac)
1011 xx(2) = epsp(i)*xfac
1018 1000
FORMAT(1x,
'RUPTURE OF SOLID ELEMENT NUMBER ',i10)
1019 1100
FORMAT(1x,
'RUPTURE OF SOLID ELEMENT NUMBER ',i10,
1020 .
' AT TIME :',g20.12)