38 1 NEL , NUPARAM, NUVAR , NFUNC , IFUNC , NPF ,
39 2 TF , TIME , TIMESTEP, UPARAM, RHO0 , RHO ,
40 3 VNEW , EINT , MUOLD ,
41 4 EPSPXX , EPSPYY , EPSPZZ , EPSPXY, EPSPYZ, EPSPZX,
42 5 DEPSXX , DEPSYY , DEPSZZ , DEPSXY, DEPSYZ, DEPSZX,
43 6 EPSXX , EPSYY , EPSZZ , EPSXY , EPSYZ , EPSZX ,
44 7 SIGOXX , SIGOYY , SIGOZZ , SIGOXY, SIGOYZ, SIGOZX,
45 8 SIGNXX , SIGNYY , SIGNZZ , SIGNXY, SIGNYZ, SIGNZX,
46 9 SIGVXX , SIGVYY , SIGVZZ , SIGVXY, SIGVYZ, SIGVZX,
47 A SOUNDSP, VISCMAX, UVAR , OFF , DVOL , VOL0 ,
48 B PM , IPM , MAT , PSH , BUFMAT,
49 C VAREOS ,NVAREOS ,MAT_PARAM,NVARTMP_EOS, VARTMP_EOS)
106#include "implicit_f.inc"
110#include "com04_c.inc"
111#include "param_c.inc"
115 INTEGER NEL, NUPARAM, NUVAR, IPM(NPROPMI,NUMGEO), MAT(NEL)
116 INTEGER,
INTENT(IN) ::NVAREOS
118 . TIME , TIMESTEP , UPARAM(NUPARAM),
119 . RHO (NEL), RHO0 (NEL), (NEL), EINT(NEL),
120 . EPSPXX(NEL), EPSPYY(NEL), EPSPZZ(NEL),
121 . EPSPXY(NEL), EPSPYZ(NEL), EPSPZX(NEL),
122 . DEPSXX(NEL), DEPSYY(NEL), DEPSZZ(NEL),
123 . DEPSXY(NEL), DEPSYZ(NEL), DEPSZX(NEL),
124 . EPSXX (NEL), EPSYY (NEL), EPSZZ (NEL),
125 . EPSXY (NEL), EPSYZ (NEL), EPSZX (NEL),
126 . SIGOXX(NEL), SIGOYY(NEL), SIGOZZ(NEL),
127 . sigoxy(nel), sigoyz(nel), sigozx(nel),
128 . dvol(nel) , vol0(nel) , pm(npropm,nummat),
129 . psh(nel) , bufmat(*) , muold(nel)
131 my_real,
INTENT(INOUT) :: vareos(nvareos*nel)
132 TYPE() ,
DIMENSION(NUMMAT),
INTENT(IN) :: MAT_PARAM
133 INTEGER,
INTENT(IN) :: NVARTMP_EOS
134 INTEGER,
INTENT(INOUT) :: VARTMP_EOS(NEL,NVARTMP_EOS)
139 . SIGNXX (NEL), SIGNYY (NEL), SIGNZZ(NEL),
140 . SIGNXY (NEL), (NEL), SIGNZX(NEL),
141 . SIGVXX (NEL), SIGVYY (NEL), SIGVZZ(NEL),
142 . SIGVXY (NEL), SIGVYZ (NEL), SIGVZX(NEL),
143 . soundsp(nel), viscmax(nel)
147 my_real uvar(nel,nuvar), off(nel)
151 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
152 my_real FINTER,FINTTE,TF(*),FINT2V
153 EXTERNAL FINTER,FINTTE
157 INTEGER I, ITER, MATS(NEL)
161 . dfs(nel), mus(nel), mu2s(nel), dpdm(nel),
162 . theta(nel), dpde(nel)
166 . ps,pe,alphae,alphap,en,beta,eta,xi,unsurn,
167 . bulk0,c0,dalpdpe,aa,hh,hh2,g,g2,dav,
168 . bulk1, crit, tol, pnew1, dfda, dgda,
171 . rhos(nel),rhos0(nel),espes(nel),
172 . vnews(nel),dvols(nel),pnew(nel),
173 . alphayld(nel),pyld(nel),alpha0(nel),
174 . alphaold(nel),dpdael(nel),eintn(nel),pold(nel)
176 . sigold(nel,6), pnewg(nel)
177 DOUBLE PRECISION ALPHA1(NEL),ALPHA(NEL),DALPHA
185 iflag1=nint(uparam(7))
186 iflag2=nint(uparam(8))
187 itemax=nint(uparam(17))
196 dgdalpha= (g-gs)/(alphae-one)
202 eostyp = mat_param(mxs)%IEOS
205 alphaold(i)= uvar(i,3)
206 alphayld(i)= uvar(i,4)
207 pold(i)=-(sigoxx(i)+sigoyy(i)+sigozz(i))*third
208 pyld(i)=ps-xi*(alphayld(i)-one)**unsurn
211 mu(i) =rho(i)/rho0(i)-one
214 sigold(i,1:3)=-pold(i)
221 IF(alphaold(i) == one)
THEN
227 hh =aa*(alphaold(i)-one)+one
229 dpdael(i)=bulk0/(alphaold(i)*(one-alphaold(i)/hh2))
230 bulk1=bulk0*(one+mu(i))/alphae
232 alpha0(i)=(pold(i)+bulk0-alphaold(i)*dpdael(i))/(bulk1-dpdael(i))
233 ELSEIF(iflag1 == 2)
THEN
234 alpha0(i)= bulk0/(bulk1-pold(i))
236 alpha0(i)=
max(alpha0(i),one)
252 rhos(i) =rho(i)*alpha(i)
253 dvols(i)=dvol(i)-vnew(i)*(one-alphaold
254 dvols(i)=dvols(i)/alphaold(i)
261 mus(i) = (one + mu(i))*alpha(i)/alphae - one
262 dfs(i) = one/(one+mus(i))
263 eintn(i)=eint(i)-half*pold(i)*dvols(i)
268 CALL eosmain(0 ,nel ,eostyp ,pm ,off ,eintn ,
269 . rhos ,rhos0 ,mus ,mu2s ,espes ,
270 . dvols ,dfs ,vnews ,mats ,psh ,
271 . pnew ,dpdm ,dpde ,theta ,
272 . bufmat ,sigold ,muold ,75 ,
273 . npf ,tf ,vareos ,nvareos,mat_param(mxs),
274 . bid ,nvartmp_eos, vartmp_eos)
275 CALL eosmain(1 ,nel ,eostyp ,pm ,off ,eintn ,
276 . rhos ,rhos0 ,mus ,mu2s ,espes ,
277 . dvols ,dfs ,vnews ,mats ,psh ,
278 . pnew ,dpdm ,dpde ,theta ,
279 . bufmat ,sigold ,muold ,75 ,
280 . npf ,tf ,vareos ,nvareos,mat_param(mxs),
281 . bid ,nvartmp_eos, vartmp_eos)
284 pnew(i)=pnew(i)/alpha(i)
290 IF (iter <= itemax)
THEN
297 dfda=dpdm(i)*(one+mu(i))/alphae + dpde(i)*eint(i)/vol0(i)
299 pnewg(i)=pold(i)+dpdael(i)*(alpha(i)-alphaold(i))
302 dalpha=(pnew(i)-pnewg(i))/(dfda-dgda)
303 ELSEIF(iflag1 == 2)
THEN
304 dalpha=(pnew(i)-pnewg(i))/((dfda-pnew(i))/alpha(i)-dgda)
306 alpha1(i)=alpha(i)-dalpha
307 alpha1(i)=
max(alpha1(i),one)
313 beta=abs(alpha(i)-alpha1(i))
330 CALL ancmsg(msgid=137,anmode=aninfo,i1=itemax)
339 IF(pnew(i) > pyld(i))
THEN
344 IF(iplas==0)
GO TO 300
349 IF(ipla(i)==1) alpha(i)=alpha0(i)
357 rhos(i) =rho(i)*alpha(i)
358 dvols(i)=dvol(i)-vnew(i)*(one-alphaold(i)/alpha(i))
359 dvols(i)=dvols(i)/alphaold(i)
360 vnews(i)=vnew(i)/alpha(i)
366 mus(i) = (one + mu(i))*alpha(i)/alphae - one
367 dfs(i) = one/(one+mus(i))
368 eintn(i)=eint(i)-half*pold(i)*dvols(i)
373 CALL eosmain(0 ,nel ,eostyp ,pm ,off ,eintn ,
374 . rhos ,rhos0 ,mus ,mu2s ,espes ,
375 . dvols ,dfs ,vnews ,mats ,psh ,
376 . pnew ,dpdm ,dpde ,theta ,
377 . bufmat,sigold,muold ,75 ,
378 . npf ,tf ,vareos ,nvareos,mat_param(mxs),
379 . bid, nvartmp_eos, vartmp_eos)
380 CALL eosmain(1 ,nel ,eostyp ,pm ,off ,eintn ,
381 . rhos ,rhos0 ,mus ,mu2s ,espes ,
382 . dvols ,dfs ,vnews ,mats ,psh ,
383 . pnew ,dpdm ,dpde ,theta ,
384 . bufmat,sigold,muold ,75 ,
385 . npf ,tf ,vareos ,nvareos,mat_param(mxs),
386 . bid, nvartmp_eos, vartmp_eos)
389 pnew(i)=pnew(i)/alpha(i)
395 IF (iter <= itemax)
THEN
401 dfda=dpdm(i)*(one+mu(i))/alphae + dpde(i)*eint(i)/vol0(i)
402 dgda=-xi*unsurn*(
max(alpha(i)-one,em20))**(unsurn-one)
403 pnewg(i)=ps-xi*(
max(alpha(i)-one,em20))**unsurn
406 dalpha=(pnew(i)-pnewg(i))/(dfda-dgda)
407 ELSEIF(iflag1 == 2)
THEN
408 dalpha=(pnew(i)-pnewg(i))/((dfda-pnew(i))/alpha(i)-dgda)
410 alpha1(i)=alpha(i)-dalpha
411 alpha1(i)=
max(alpha1(i),one)
417 beta=abs(alpha(i)-alpha1(i))
434 CALL ancmsg(msgid=137,anmode=aninfo,i1=itemax)
442 IF(alpha(i) > one)
THEN
443 IF(ipla(i) == 1)
THEN
444 IF(pnew(i) < ps)
THEN
445 alphayld(i)= one+ eta*(ps-pnew(i))**en
457 dav = (depsxx(i)+depsyy(i)+depszz(i))*third
458 g =
max(zero,gs+dgdalpha*(alpha(i)-one))
460 signxx(i)=sigoxx(i)+pold(i)+g2*(depsxx(i)-dav)
461 signyy(i)=sigoyy(i)+pold(i)+g2*(depsyy(i)-dav)
462 signzz(i)=sigozz(i)+pold(i)+g2*(depszz(i)-dav)
463 signxy(i)=sigoxy(i)+g*depsxy(i)
464 signyz(i)=sigoyz(i)+g*depsyz(i)
465 signzx(i)=sigozx(i)+g*depszx(i)
474 uvar(i,4)=alphayld(i)
475 signxx(i)=signxx(i)-pnew(i)
476 signyy(i)=signyy(i)-pnew(i)
477 signzz(i)=signzz(i)-pnew(i)
483 IF(alpha(i) > one)
THEN
484 soundsp(i)=sqrt(bulk0/rhos0(1))*(aa*(alpha(i)-one)+one)
486 soundsp(i)=sqrt(bulk0/rhos0(1))
subroutine sigeps75(nel, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, rho0, rho, vnew, eint, muold, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigvxx, sigvyy, sigvzz, sigvxy, sigvyz, sigvzx, soundsp, viscmax, uvar, off, dvol, vol0, pm, ipm, mat, psh, bufmat, vareos, nvareos, mat_param, nvartmp_eos, vartmp_eos)