35 1 NEL ,NGL ,NUPARAM ,NUVAR ,NVARTMP ,NUMTABL ,
36 2 UPARAM ,UVAR ,VARTMP ,ITABLE ,TABLE ,JTHE ,
37 3 TIME ,TIMESTEP,OFF ,RHO ,PLA ,DPLA ,
38 4 SOUNDSP ,SIGY ,ET ,TEMP ,EPSD ,DPDM ,
39 5 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
40 6 SIGOXX ,SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
41 7 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
42 8 INLOC ,DPLANL ,IEOS ,JLAG ,FHEAT ,VOLUME ,
43 9 SEQ ,LPLANL ,PLA_NL ,LEPSDNL ,DPDT_NL )
52#include "implicit_f.inc"
60 INTEGER ,
INTENT(IN) :: IEOS
61 INTEGER ,
INTENT(IN) :: JTHE
62 INTEGER ,
INTENT(IN) :: JLAG
63 INTEGER NEL,NUPARAM,NUVAR,NVARTMP,NUMTABL,INLOC
64 INTEGER ,
DIMENSION(NUMTABL),
INTENT(IN) :: ITABLE
65 INTEGER ,
DIMENSION(NEL) ,
INTENT(IN) :: NGL
67 my_real :: TIME,TIMESTEP
68 my_real,
DIMENSION(NUPARAM) ,
INTENT(IN) :: UPARAM
69 my_real,
DIMENSION(NEL) ,
INTENT(IN) :: volume
70 my_real,
DIMENSION(NEL) ,
INTENT(IN) :: rho,dpdm,off,
71 . depsxx,depsyy,depszz,depsxy,depsyz,depszx,
72 . sigoxx,sigoyy,sigozz,sigoxy,sigoyz,sigozx,dplanl
73 my_real ,
DIMENSION(NEL) ,
INTENT(OUT) :: soundsp,sigy,et,
74 . signxx,signyy,signzz,signxy,signyz,signzx
75 my_real ,
DIMENSION(NEL) ,
INTENT(INOUT) :: pla,dpla,epsd
76 my_real ,
DIMENSION(NEL,NUVAR) ,
INTENT(INOUT) :: uvar
77 INTEGER ,
DIMENSION(NEL,NVARTMP) ,
INTENT(INOUT) :: VARTMP
78 my_real ,
DIMENSION(NEL) ,
INTENT(INOUT) :: temp
79 my_real ,
DIMENSION(NEL) ,
INTENT(INOUT) :: fheat
80 my_real ,
DIMENSION(NEL) ,
INTENT(INOUT) :: seq
81 INTEGER,
INTENT(IN) :: LPLANL, LEPSDNL
82 my_real,
DIMENSION(NEL*LPLANL) :: pla_nl
83 my_real,
DIMENSION(NEL*LEPSDNL) :: dpdt_nl
85 TYPE(
ttable),
DIMENSION(NTABLE) :: TABLE
89 INTEGER I,II,NINDX,ITER,NITER,ISMOOTH,
90 . func_yld,func_temp,func_eta,ndim_yld,ndim_temp,ndim_eta
91 INTEGER ,
DIMENSION(NEL) :: INDEX
93 my_real :: young,bulk,lame,g
94 . xrate,xscale,yscale,dtinv,j2,q2,dphi_dlam
95 my_real :: asrate,ddep,dfdsig2,sig_dfdsig
97 my_real,
DIMENSION(NEL) ::svm
98 . sxx,syy,szz,sxy,syz,szx,sigm,
99 . fact_eta,dydx,hardp,hardr,yld_i,hardp_i,hardr_i,dxdyv,dlam,phi,
100 . ftherm,tfac,pla0,normxx,normyy,normzz,normxy,
101 . normyz,normzx,dpxx,dpyy,dpzz,dpxy,dpyz,dpzx,dpla_dlam
102 my_real,
DIMENSION(NEL,3) :: xvec_eta
103 my_real,
DIMENSION(NEL,4) :: xvec
104 INTEGER,
DIMENSION(NEL,3) :: IPOS_ETA
105 INTEGER,
DIMENSION(NEL,2) ::
123 ismooth = nint(uparam(6))
135 dtinv = one /
max(em20, timestep)
136 alpha = asrate*timestep
139 IF (jthe /= 0) eta = zero
143 func_temp = itable(2)
145 ndim_yld = table(func_yld)%NDIM
146 IF (func_temp > 0)
THEN
147 ndim_temp = table(func_temp)%NDIM
149 IF (func_eta > 0)
THEN
150 ndim_eta = table(func_eta)%NDIM
164 IF (time == zero) temp(1:nel) = tini
168 IF (func_eta > 0)
THEN
170 xvec_eta(1:nel,1) = epsd(1:nel) * xrate
172 xvec_eta(1:nel,1) = dpdt_nl(1:nel
174 ipos_eta(1:nel,1) = 1
175 IF (ndim_eta > 1)
THEN
176 xvec_eta(1:nel,2) = temp(1:nel)
177 ipos_eta(1:nel,2) = vartmp(1:nel,4)
179 IF (ndim_eta > 2)
THEN
181 xvec_eta(1:nel,3) = pla(1:nel)
183 xvec_eta(1:nel,3) = pla_nl(1:nel)
185 ipos_eta(1:nel,3) = vartmp(1:nel,5)
187 CALL table_vinterp(table(func_eta),nel,nel,ipos_eta,xvec_eta,
189 IF (ndim_eta > 1) vartmp(1:nel,4) = ipos_eta(1:nel,2)
190 IF (ndim_eta > 2) vartmp(1:nel,5) = ipos_eta(1:nel,3)
192 ftherm(i) =
min(eta*fact_eta(i), one)
195 ftherm(1:nel) =
min(eta, one)
205 ldav = (depsxx(i) + depsyy(i) + depszz(i)) * lame
206 signxx(i) = sigoxx(i) + depsxx(i)*g2 + ldav
207 signyy(i) = sigoyy(i) + depsyy(i)*g2 + ldav
208 signzz(i) = sigozz(i) + depszz(i)*g2 + ldav
210 signyz(i) = sigoyz(i) + depsyz(i)*g
211 signzx(i) = sigozx(i) + depszx(i)*g
213 sigm(i) = (signxx(i) + signyy(i) + signzz(i)) * third
215 sxx(i) = signxx(i) - sigm(i)
216 syy(i) = signyy(i) - sigm(i)
217 szz(i) = signzz(i) - sigm(i)
222 j2 = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) +
223 . sxy(i)**2 + syz(i)**2 + szx(i)**2
224 svm(i) = sqrt(three*j2)
228 xvec(1:nel,1) = pla(1:nel)
229 xvec(1:nel,2) = epsd(1:nel) * xscale
230 ipos(1:nel,1) = vartmp(1:nel,1)
233 . xvec,yld,hardp,hardr)
234 yld(1:nel) = yld(1:nel)*yscale
235 hardp(1:nel) = hardp(1:nel)*yscale
236 vartmp(1:nel,1) = ipos(1:nel,1)
238 IF (func_temp > 0)
THEN
240 ipos(1:nel,1) = vartmp(1:nel,2)
241 ipos(1:nel,2) = vartmp(1:nel,3)
242 CALL table_vinterp(table(func_temp),nel,nel,ipos,xvec,yld_tref,dydx)
243 vartmp(1:nel,2) = ipos(1:nel,1)
244 vartmp(1:nel,3) = ipos(1:nel,2)
245 xvec(1:nel,2) = temp(1:nel)
246 CALL table_vinterp(table(func_temp),nel,nel,ipos,xvec,yld_temp,dydx)
247 tfac(1:nel) = yld_temp(1:nel) / yld_tref(1:nel)
248 yld(1:nel) = yld(1:nel) * tfac(1:nel)
249 hardp(1:nel) = hardp(1:nel) * tfac(1:nel)
257 phi(1:nel) = svm(1:nel) - yld(1:nel)
260 IF (phi(i) >= zero .AND. off(i) == one)
THEN
277#include "vectorize.inc"
294 normxx(i) = three_half*sxx(i)/(
max(svm(i),em20))
295 normyy(i) = three_half*syy(i)/(
max(svm(i),em20))
296 normzz(i) = three_half*szz(i)/(
max(svm(i),em20))
297 normxy(i) = three*sxy(i)/(
max(svm(i),em20))
298 normyz(i) = three*syz(i)/(
max(svm(i),em20))
299 normzx(i) = three*szx(i)/(
max(svm(i),em20))
306 dfdsig2 = normxx(i) * normxx(i) * g2
307 . + normyy(i) * normyy(i) * g2
308 . + normzz(i) * normzz(i) * g2
309 . + normxy(i) * normxy(i) * g
310 . + normyz(i) * normyz(i) * g
311 . + normzx(i) * normzx(i) * g
315 sig_dfdsig = signxx(i) * normxx(i)
316 . + signyy(i) * normyy(i)
317 . + signzz(i) * normzz(i)
318 . + signxy(i) * normxy(i)
320 . + signzx(i) * normzx(i)
321 dpla_dlam(i) = sig_dfdsig / yld(i)
325 dphi_dlam = - dfdsig2 + hardp(i)*dpla_dlam(i)
330 dlam(i) = - phi(i) / dphi_dlam
335 ddep = dpla_dlam(i)*dlam(i)
337 dpla(i) =
max(dpla(i) + ddep,zero)
339 pla(i) = pla0(i) + dpla(i)
341 dpxx(i) = dlam(i)*normxx(i)
342 dpyy(i) = dlam(i)*normyy(i)
343 dpzz(i) = dlam(i)*normzz(i)
344 dpxy(i) = dlam(i)*normxy(i)
345 dpyz(i) = dlam(i)*normyz(i)
346 dpzx(i) = dlam(i)*normzx(i)
352 xvec(1:nel,1:2) = zero
358 ipos(ii,1) = vartmp(i,1)
361 CALL table2d_vinterp_log(table(func_yld),ismooth,nel,nindx,ipos,xvec,yld_i,hardp_i,hardr_i)
364 vartmp(i,1) = ipos(ii,1)
365 hardp(i) = hardp_i(ii)*yscale*tfac(i)
366 yld(i) = yld_i(ii)*yscale*tfac(i)
371#include "vectorize.inc"
376 signxx(i) = signxx(i) - dpxx(i)*g2
377 signyy(i) = signyy(i) - dpyy(i)*g2
378 signzz(i) = signzz(i) - dpzz(i)*g2
379 signxy(i) = signxy(i) - dpxy(i)*g
380 signyz(i) = signyz(i) - dpyz(i)*g
381 signzx(i) = signzx(i) - dpzx(i)*g
384 sigm(i) = (signxx(i) + signyy(i) + signzz(i)) * third
387 sxx(i) = signxx(i) - sigm(i)
388 syy(i) = signyy(i) - sigm(i)
389 szz(i) = signzz(i) - sigm(i)
395 j2 = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) +
396 . sxy(i)**2 + syz(i)**2 + szx(i)**2
397 svm(i) = sqrt(three*j2)
400 phi(i) = svm(i) - yld(i)
414 et(i) = hardp(i) / (hardp(i) + young)
418 IF (eta > zero .AND. inloc == 0)
THEN
419 IF (jthe /= 0 .AND. jlag /= 0)
THEN
422 fheat(i) = fheat(i) + ftherm(i)*yld(i)*dpla(i) * volume(i)
424 ELSE IF (cp > zero)
THEN
427 temp(i) = temp(i) + ftherm(i)*yld(i)*dpla(i) / (cp*rho(i))
436 IF (inloc > 0 .AND. eta > zero)
THEN
437 IF (jthe /= 0 .AND. jlag /= 0)
THEN
439 IF (off(i) == one)
THEN
440 fheat(i) = fheat(i) +
441 . ftherm(i)*yld(i)*
max(dplanl(i),zero)*volume(i)
446 IF (off(i) == one)
THEN
447 temp(i) = temp(i) + ftherm(i)*yld(i)*
max(dplanl(i),zero)/(cp*rho(i))
459 signxx(1:nel) = signxx(1:nel) - sigm(1:nel)
460 signyy(1:nel) = signyy(1:nel) - sigm(1:nel)
461 signzz(1:nel) = signzz(1:nel) - sigm(1:nel)
462 signxy(1:nel) = signxy(1:nel)
463 signyz(1:nel) = signyz(1:nel)
464 signzx(1:nel) = signzx(1:nel)
465 soundsp(1:nel)= sqrt((dpdm(1:nel) + four_over_3*g) / rho(1:nel))
467 soundsp(1:nel)= sqrt((bulk + four_over_3*g) / rho(1:nel))
478 epsd(i) =
alpha*dpla(i)*dtinv + alphi*epsd(i)
subroutine sigeps109(nel, ngl, nuparam, nuvar, nvartmp, numtabl, uparam, uvar, vartmp, itable, table, jthe, time, timestep, off, rho, pla, dpla, soundsp, sigy, et, temp, epsd, dpdm, depsxx, depsyy, depszz, depsxy, depsyz, depszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, inloc, dplanl, ieos, jlag, fheat, volume, seq, lplanl, pla_nl, lepsdnl, dpdt_nl)