36 1 NEL, NUPARAM, NUVAR, MFUNC,
37 2 KFUNC, NPF, TF, TIME,
38 3 TIMESTEP,UPARAM, RHO0, RHO,
39 4 OFFG, RHOREF, NVARTMP, VARTMP,
41 7 DEPSYY, DEPSZZ, DEPSXY, DEPSYZ,
42 8 DEPSZX, EPSXX, EPSYY, EPSZZ,
43 9 EPSXY, EPSYZ, EPSZX, SIG0XX,
44 A SIG0YY, SIG0ZZ, SIG0XY, SIG0YZ,
45 B SIG0ZX, SIGNXX, SIGNYY, SIGNZZ,
46 C SIGNXY, SIGNYZ, SIGNZX,
47 E SOUNDSP, VISCMAX, UVAR,
49 G MAT, EPSP, ET, ISMSTR,
50 H IHET, JSMS, MATPARAM)
59#include
"implicit_f.inc"
102#include "param_c.inc"
104#include "impl1_c.inc"
106 INTEGER ,
INTENT(IN) :: NVARTMP
107 INTEGER ,
INTENT(INOUT) :: VARTMP(NEL,NVARTMP)
108 INTEGER,
INTENT(IN) :: ISMSTR, IHET, JSMS
109 INTEGER NEL, NUPARAM, NUVAR,
110 . NGL(NEL),MAT(),IPM(NPROPMI,*)
112 . TIME,TIMESTEP,UPARAM(*),
113 . RHO(NEL),RHO0(NEL),
114 . DEPSXX(NEL),DEPSYY(NEL),DEPSZZ(NEL),
115 . (NEL),DEPSYZ(NEL),DEPSZX(NEL),
116 . EPSXX(NEL) ,EPSYY(NEL) ,EPSZZ(NEL) ,
117 . EPSXY(NEL) ,EPSYZ(NEL) ,EPSZX(NEL) ,
118 . SIG0XX(NEL),(NEL),SIG0ZZ(NEL),
119 . SIG0XY(NEL),SIG0YZ(NEL),SIG0ZX(NEL),
121 TYPE(matparam_struct_) ,
INTENT(IN) :: MATPARAM
126 . signxx(nel),signyy(nel),signzz(nel),
127 . signxy(nel),signyz(nel),signzx(nel),
128 . soundsp(nel),viscmax(nel),et(nel),
134 . uvar(nel,nuvar), off(nel), pla(nel), offg(nel), rhoref(mvsiz)
138 INTEGER NPF(*), MFUNC, (MFUNC)
153 INTEGER :: I,II,J,J1,J2,IADBUF,NFUNC,
154 . nindex_epssm,nindex_epst,nindex_mu,nload,nunload,mx,
155 . iload0,idamage,itens,numtabl
156 INTEGER,
DIMENSION(100) :: IFUNC
157 INTEGER,
DIMENSION(NEL) :: JJ,IE_CST,ILOAD
158 INTEGER,
DIMENSION(NEL) :: INDEX_EPSSM,INDEX_EPST,INDEX_MU
159 INTEGER,
DIMENSION(NEL) :: IPOS1,IADP1,ILENP1
160 INTEGER,
DIMENSION(NEL,2) :: IPOS
161 my_real :: YLD_EMAX,R,FAC,YP1,YP2,EPSP0,EPSPMAX,E,DELTA,SVM2,SVM,
162 . AA,E0,NU,EMAX,EPSSMAX,DSIG,P,EXPO,HYS,DE,YFACC,MU,ALPHA
163 my_real,
DIMENSION(NEL) :: ALPHA_1,AA1,AA2,EPST,YLDMIN,YLDMAX,YLDELAS,
164 . G,SIG0,EPS0,EPSS,TAB_LOC1,TAB_LOC2,,DFC
165 my_real,
DIMENSION(NEL,2) :: XVEC
170 ifunc(j)=ipm(10+j,mx)
174 e0 = uparam(iadbuf + 2)
175 aa = uparam(iadbuf + 3)
176 epssmax = uparam(iadbuf + 4)
177 g(:) = uparam(iadbuf + 5)
178 nu = uparam(iadbuf + 6)
179 epsp0 = uparam(iadbuf + 9)
180 nload = nint(uparam(iadbuf + 7))
181 nunload = nint(uparam(iadbuf + 8))
182 idamage = nint(uparam(iadbuf + 2*nfunc + 9))
183 itens = nint(uparam(iadbuf + 2*nfunc + 13))
184 yfacc = uparam(iadbuf + 2*nfunc + 8)
185 expo = uparam(iadbuf + 2*nfunc + 10)
186 hys = uparam(iadbuf + 2*nfunc + 11)
187 emax = uparam(iadbuf + 2*nfunc + 12)
188 yld_emax= uparam(iadbuf + 2*nfunc + 15)
189 numtabl = matparam%NTABLE
194 epst(i) = epsxx(i)**2 + epsyy(i)**2 + epszz(i)**2
195 . + half*(epsxy(i)**2+epsyz(i)**2+epszx(i)**2)
196 epst(i) = sqrt(epst(i))
208 ipos(i,1) = vartmp(i,1)
215 vartmp(1:nel,1)=ipos(1:nel,1)
218 IF (epst(i) >= epssmax)
THEN
219 nindex_epssm = nindex_epssm + 1
220 index_epssm(nindex_epssm) = i
221 yldelas(i) = yld_emax + emax*(epst(i) - epssmax)
223 nindex_epst = nindex_epst + 1
224 index_epst(nindex_epst) = i
232 epspmax = matparam%TABLE(1)%X(2)%VALUES(nload)
233#include "vectorize.inc"
235 ipos(i,1) = vartmp(i,2)
236 ipos(i,2) = vartmp(i,5)
237 xvec(i,1) =
min(epst(i), epssmax)
245 ELSE IF (nload == 1)
THEN
247 ipos(:,1) = vartmp(:,2)
248 ipos(:,2) = vartmp(:,5)
250 xvec(i,1) =
min(epst(i), epssmax)
257 vartmp(1:nel,2) = ipos(1:nel,1)
258 vartmp(1:nel,5) = ipos(1:nel,2)
264 IF (nindex_epssm > 0)
THEN
265#include "vectorize.inc"
268 yldmax(i) = yldmax(i) + emax*(epst(i) - epssmax)
274 IF (nunload > 0)
THEN
275 IF (nunload > 1)
THEN
276 epspmax = matparam%TABLE(1)%X(2)%VALUES(nunload)
277#include "vectorize.inc"
279 ipos(i,1) = vartmp(i,3)
280 ipos(i,2) = vartmp(i,6)
281 xvec(i,1) =
min(epst(i), epssmax)
287 ELSE IF (nunload == 1)
THEN
288 ipos(:,1) = vartmp(:,3)
289 ipos(:,2) = vartmp(:,6)
298 vartmp(1:nel,3)=ipos(1:nel,1)
299 vartmp(1:nel,6)=ipos(1:nel,2)
305 IF (nindex_epssm > 0)
THEN
306#include "vectorize.inc"
309 yldmin(i) = yldmin(i) + emax*(epst(i) - epssmax)
319 alpha = one/
max(em20,uvar(i,10))
320 sig0xx(i) = alpha*sig0xx(i)
321 sig0yy(i) = alpha*sig0yy(i)
322 sig0zz(i) = alpha*sig0zz(i)
323 sig0xy(i) = alpha*sig0xy(i)
324 sig0zx(i) = alpha*sig0zx(i)
325 sig0yz(i) = alpha*sig0yz(i)
330 iload0 = int(uvar(i,5))
332 delta = epst(i) - uvar(i,4)
335 IF (delta >= zero )
THEN
338 ELSEIF (delta < zero)
THEN
341 IF (idamage /= 0 )
THEN
346 epss(i) = epst(i) - yld(i)/ e
347 epss(i) =
max(zero, epss(i))
348 de = aa*(epss(i) - uvar(i,1))
349 IF (iload(i) == 1)
THEN
350 e = e +
max(de, zero)
351 IF(iload0 == -1) e= uvar(i,3)
352 uvar(i,1) =
max(uvar(i,1), epss(i))
354 e = e +
min(de ,zero)
355 IF(iload0 == 1) e= uvar(i,3)
356 uvar(i,1) =
min(epss(i),uvar(i,1))
362 aa1(i) = e*(one-nu)/(one + nu)/(one - two*nu)
363 aa2(i) = aa1(i)*nu/(one - nu)
364 g(i) = half*e/(one + nu)
366 signxx(i)= aa1(i)*depsxx(i) + aa2(i)*(depsyy(i) + depszz(i))
367 signyy(i)= aa1(i)*depsyy(i) + aa2(i)*(depsxx(i) + depszz(i))
368 signzz(i)= aa1(i)*depszz(i) + aa2(i)*(depsxx(i) + depsyy(i))
369 signxy(i)= g(i) *depsxy(i)
370 signyz(i)= g(i) *depsyz(i)
371 signzx(i)= g(i) *depszx(i)
372 dsig = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
373 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
376 signxx(i)=sig0xx(i) + aa1(i)*depsxx(i)
377 . + aa2(i)*(depsyy(i) + depszz(i))
378 signyy(i)=sig0yy(i) + aa1(i)*depsyy(i)
379 . + aa2(i)*(depsxx(i) + depszz(i))
380 signzz(i)=sig0zz(i) + aa1(i)*depszz(i)
381 . + aa2(i)*(depsxx(i) + depsyy(i))
382 signxy(i)=sig0xy(i) + g(i) *depsxy(i)
383 signyz(i)=sig0yz(i) + g(i) *depsyz(i)
384 signzx(i)=sig0zx(i) + g(i) *depszx(i)
386 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
387 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
390 IF (idamage == 0 )
THEN
391 IF(delta >= zero )
THEN
392 yld(i) =
min(yldmax(i),svm)
393 ELSEIF(delta < zero)
THEN
394 yld(i) =
max(yldmin(i),svm)
395 yld(i) =
min(yld(i),uvar(i,6))
397 IF (dsig > yldmin(i) .AND. dsig > svm)
THEN
404 uvar(i,8) = uvar(i,8) + half*(yldelas(i) + uvar(i,9))*delta
405 uvar(i,8) =
max(zero, uvar(i,8))
406 uvar(i,2) =
max(uvar(i,2) , uvar(i,8))
408 uvar(i,9) = yldelas(i)
414 IF (idtmins/=2.OR.jsms==0)
THEN
420 soundsp(i) = sqrt(aa1(i)/rho0(i))
424 ELSEIF (ismstr==1 .OR. ismstr==11)
THEN
430 soundsp(i) = sqrt(aa1(i)/rhoref(i))
436 IF (abs(offg(i)) <= one)
THEN
441 soundsp(i) = sqrt(aa1(i
448 soundsp(i) = sqrt(aa1(i)/rhoref(i
458 IF (idamage == 0 )
THEN
460 signxx(i)= aa1(i)*epsxx(i) + aa2(i)*(epsyy(i) + epszz(i))
461 signyy(i)= aa1(i)*epsyy(i) + aa2(i)*(epsxx(i) + epszz(i))
462 signzz(i)= aa1(i)*epszz(i) + aa2(i)*(epsxx(i) + epsyy(i))
463 signxy(i)= g(i) *epsxy(i)
464 signyz(i)= g(i) *epsyz(i)
465 signzx(i)= g(i) *epszx(i
467 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
468 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
470 r = yld(i)/
max(em20,svm)
472 signxx(i)=signxx(i)*r
473 signyy(i)=signyy(i)*r
474 signzz(i)=signzz(i)*r
475 signxy(i)=signxy(i)*r
476 signyz(i)=signyz(i)*r
477 signzx(i)=signzx(i)*r
485 IF (ie_cst(i) == 1)
THEN
486 iload0 = int(uvar(i,5))
487 IF(iload0 /= iload(i))
THEN
499 ELSEIF (idamage == 1)
THEN
503 signxx(i)= aa1(i)*epsxx(i) + aa2(i)*(epsyy(i) + epszz(i))
504 signyy(i)= aa1(i)*epsyy(i) + aa2(i)*(epsxx(i) + epszz(i))
505 signzz(i)= aa1(i)*epszz(i) + aa2(i)*(epsxx(i) + epsyy(i))
506 signxy(i)= g(i) *epsxy(i)
507 signyz(i)= g(i) *epsyz(i)
508 signzx(i)= g(i) *epszx(i)
510 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
511 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
513 r = (yld(i)/
max(em20,svm))
515 signxx(i)=signxx(i)*r
516 signyy(i)=signyy(i)*r
517 signzz(i)=signzz(i)*r
518 signxy(i)=signxy(i)*r
519 signyz(i)=signyz(i)*r
520 signzx(i)=signzx(i)*r
522 IF (iload(i) == -1)
THEN
523 r = (yldmin(i)/
max(em20,yldelas(i)))
524 p = third*(signxx(i) + signyy(i) + signzz(i))
525 signxx(i)=(signxx(i) - p)*r + p
526 signyy(i)=(signyy(i) - p)*r + p
527 signzz(i)=(signzz(i) - p)*r + p
528 signxy(i)=signxy(i)*r
529 signyz(i)=signyz(i)*r
542 ELSEIF(idamage == 2 )
THEN
546 signyy(i)= aa1(i)*epsyy(i) + aa2(i)*(epsxx(i) + epszz(i))
548 signxy(i)= g(i) *epsxy(i)
549 signyz(i)= g(i) *epsyz(i)
552 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
553 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
555 r = yld(i)/
max(em20,svm)
557 signxx(i)=signxx(i)*r
558 signyy(i)=signyy(i)*r
559 signzz(i)=signzz(i)*r
560 signxy(i)=signxy(i)*r
561 signyz(i)=signyz(i)*r
562 signzx(i)=signzx(i)*r
564 IF (iload(i) == -1)
THEN
565 r = yldmin(i)/
max(em20,yldelas(i))
566 signxx(i)=signxx(i)*r
567 signyy(i)=signyy(i)*r
568 signzz(i)=signzz(i)*r
569 signxy(i)=signxy(i)*r
570 signyz(i)=signyz(i)*r
571 signzx(i)=signzx(i)*r
575 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
576 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
584 ELSEIF (idamage == 3)
THEN
587 signxx(i)= aa1(i)*epsxx(i) + aa2(i)*(epsyy(i) + epszz(i))
588 signyy(i)= aa1(i)*epsyy(i) + aa2(i)*(epsxx(i) + epszz(i))
589 signzz(i)= aa1(i)*epszz(i) + aa2(i)*(epsxx(i) + epsyy(i))
590 signxy(i)= g(i) *epsxy(i)
591 signyz(i)= g(i) *epsyz(i)
592 signzx(i)= g(i) *epszx(i)
594 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
595 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
597 r = (yld(i)/
max(em20,svm))
598 signxx(i)=signxx(i)*r
599 signyy(i)=signyy(i)*r
600 signzz(i)=signzz(i)*r
601 signxy(i)=signxy(i)*r
602 signyz(i)=signyz(i)*r
603 signzx(i)=signzx(i)*r
606 IF (iload(i)==-1 .AND.uvar(i,2)/=zero)
THEN
607 r = one - (uvar(i,8)/uvar(i,2))**expo
608 r = one - (one - hys)*r
609 p = third*(signxx(i) + signyy(i) + signzz(i))
610 signxx(i)=(signxx(i) - p)*r + p
611 signyy(i)=(signyy(i) - p)*r + p
612 signzz(i)=(signzz(i) - p)*r + p
613 signxy(i)=signxy(i)*r
614 signyz(i)=signyz(i)*r
615 signzx(i)=signzx(i)*r
625 ELSEIF (idamage == 4)
THEN
628 signxx(i)= aa1(i)*epsxx(i) + aa2(i)*(epsyy(i) + epszz(i))
629 signyy(i)= aa1(i)*epsyy(i) + aa2(i)*(epsxx(i) + epszz(i))
630 signzz(i)= aa1(i)*epszz(i) + aa2(i)*(epsxx(i) + epsyy(i))
631 signxy(i)= g(i) *epsxy(i)
632 signyz(i)= g(i) *epsyz(i)
633 signzx(i)= g(i) *epszx(i)
635 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
636 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
638 r = (yld(i)/
max(em20,svm))
640 signxx(i)=signxx(i)*r
641 signyy(i)=signyy(i)*r
642 signzz(i)=signzz(i)*r
643 signxy(i)=signxy(i)*r
644 signyz(i)=signyz(i)*r
645 signzx(i)=signzx(i)*r
647 IF (iload(i)==-1 .AND.uvar(i,2)/= zero)
THEN
648 r = one - (uvar(i,8)/uvar
649 r = one - (one - hys)*r
650 signxx(i)=signxx(i)*r
651 signyy(i)=signyy(i)*r
652 signzz(i)=signzz(i)*r
653 signxy(i)=signxy(i)*r
654 signyz(i)=signyz(i)*r
655 signzx(i)=signzx(i)*r
670 mu = 1 - rho(i)/rho0(i)
672 nindex_mu = nindex_mu + 1
673 index_mu(nindex_mu) = i
674 ipos1(nindex_mu) = vartmp(i,4)
675 iadp1(nindex_mu) = npf(ifunc(nfunc)) / 2 + 1
676 ilenp1(nindex_mu) = npf(ifunc(nfunc)+1) / 2 - iadp1(nindex_mu) - ipos1(nindex_mu)
677 tab_loc1(nindex_mu) = mu
681 IF ( nindex_mu > 0 )
THEN
683 CALL vinter2(tf,iadp1,ipos1,ilenp1,nindex_mu,tab_loc1,dfc,yld)
687 alpha_1(i) = yfacc*yld(ii)
688 alpha_1(i) =
max(zero, alpha_1(i))
689 vartmp(i,4) = ipos1(ii)
694 signxx(i) = alpha_1(i)*signxx(i)
695 signyy(i) = alpha_1(i)*signyy(i)
696 signzz(i) = alpha_1(i)*signzz(i)
697 signxy(i) = alpha_1(i)*signxy(i)
698 signzx(i) = alpha_1(i)*signzx(i)
699 signyz(i) = alpha_1(i)*signyz(i)
700 uvar(i,10) = alpha_1(i)
705 IF (impl_s > 0 .OR. ihet > 1)
THEN
707 et(i) =
min(one , et(i))
subroutine sigeps70(nel, nuparam, nuvar, mfunc, kfunc, npf, tf, time, timestep, uparam, rho0, rho, offg, rhoref, nvartmp, vartmp, rhosp, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, sig0xx, sig0yy, sig0zz, sig0xy, sig0yz, sig0zx, signxx, signyy, signzz, signxy, signyz, signzx, soundsp, viscmax, uvar, off, ngl, ipm, mat, epsp, et, ismstr, ihet, jsms, matparam)
subroutine mulaw(lft, llt, nft, mtn, jcvt, pm, off, sig, eint, rho, vol, strain, gama, uvar, bufmat, tf, npf, imat, ngl, nuvar, nvartmp, vartmp, geo, pid, epsd, wxx, wyy, wzz, jsph, ssp, voln, vis, d1, d2, d3, d4, d5, d6, dvol, sold1, sold2, sold3, sold4, sold5, sold6, rx, ry, rz, sx, sy, sz, tx, ty, tz, ismstr, mfxx, mfxy, mfxz, mfyx, mfyy, mfyz, mfzx, mfzy, mfzz, ipm, isorth, nel, matparam)