32 1 NEL0 ,NUPARAM,NUVAR ,NFUNC ,IFUNC ,
34 2 TF ,TIME ,TIMESTEP,UPARAM ,RHO0 ,
36 4 EPSPXX ,EPSPYY ,EPSPZZ ,EPSPXY ,EPSPYZ ,EPSPZX ,
37 5 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
38 6 EPSXX ,EPSYY ,EPSZZ ,EPSXY ,EPSYZ ,EPSZX ,
39 7 SIGOXX ,SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
40 8 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
41 9 SIGVXX ,SIGVYY ,SIGVZZ ,SIGVXY ,SIGVYZ ,SIGVZX ,
42 A SOUNDSP,VISCMAX,UVAR ,OFF ,SIGY ,DEFP ,
47#include "implicit_f.inc"
100 INTEGER NEL0, NUPARAM, NUVAR
102 . TIME,TIMESTEP,UPARAM(NUPARAM),
103 . RHO(NEL0),RHO0(NEL0),VOLUME(NEL0),EINT(NEL0),
104 . EPSPXX(NEL0),EPSPYY(NEL0),EPSPZZ(NEL0),
105 . EPSPXY(NEL0),EPSPYZ(NEL0)
108 . EPSXX(NEL0) ,EPSYY(NEL0) ,EPSZZ(NEL0) ,
109 . (NEL0) ,EPSYZ(NEL0) ,EPSZX(NEL0) ,
110 . SIGOXX(NEL0),SIGOYY(NEL0),SIGOZZ(NEL0),
111 . SIGOXY(NEL0),SIGOYZ(NEL0),SIGOZX(NEL0),
117 . signxx(nel0),signyy(nel0),signzz(nel0),
118 . signxy(nel0),signyz(nel0),signzx(nel0),
119 . sigvxx(nel0),sigvyy(nel0),sigvzz(nel0),
120 . sigvxy(nel0),sigvyz(nel0),sigvzx(nel0),
121 . soundsp(nel0),viscmax(nel0),sigy(nel0),defp(nel0)
125 my_real uvar(nel0,nuvar), off(nel0)
129 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
130 my_real FINTER ,TF(*)
144 . g,ca,cb,cn,epsm,sigm,cc,cd,cm,eps0,
146 . ch1,ch2,ch3,qh1,qh2,svm,
148 . eprnxx,eprnyy,eprnzz,eprnxy,eprnyz,eprnzx,
149 . evr,edrnxx,edrnyy,edrnzz,edrnxy,edrnyz,edrnzx,
196 beta = timestep*two*pi*cutfre
206 press = - (sigoxx(i) + sigoyy(i) + sigozz(i)) * third
215 evr = - (eprnxx + eprnyy + eprnzz) * third
216 edrnxx = eprnxx + evr
217 edrnyy = eprnyy + evr
218 edrnzz = eprnzz + evr
223 signxx(i)=sigoxx(i) + press + dtg2*edrnxx
224 signyy(i)=sigoyy(i) + press + dtg2*edrnyy
225 signzz(i)=sigozz(i) + press + dtg2*edrnzz
226 signxy(i)=sigoxy(i) + dtg2*edrnxy
227 signyz(i)=sigoyz(i) + dtg2*edrnyz
228 signzx(i)=sigozx(i) + dtg2*edrnzx
230 soundsp(i) = sqrt(c14g3/rho0(i))
239 uvar(i,2) = half*(edrnxx*edrnxx+edrnyy*edrnyy+edrnzz*edrnzz)
240 . +edrnxy*edrnxy+edrnyz*edrnyz+edrnzx*edrnzx
241 uvar(i,2) = off(i) * sqrt(three*uvar(i,2)) / three_half
243 uvar(i,2) = beta*uvar(i,2) + (one - beta)*uvar(i,5)
244 uvar(i,5) = uvar(i,2)
251 IF(uvar(i,1)<=zero)
THEN
253 ELSEIF(uvar(i,1)>epsm)
THEN
256 ch1=ca+cb*uvar(i,1)**cn
258 IF(uvar(i,2)<=eps0)
THEN
260 ELSEIF(uvar(i,1)<=zero)
THEN
261 ch2=cc*log(uvar(i,2)/eps0)
263 ch2=(cc-cd*uvar(i,1)**cm)*log(uvar(i,2)/eps0)
265 IF(uvar(i,2)<=zero)
THEN
271 uvar(i,3)=
min(sigm+ch3,ch1+ch2+ch3)
273 IF(uvar(i,1)>epsm) uvar(i,3)=zero
279 IF(uvar(i,1)>zero. and .cn>=one)
THEN
280 qh1= cb*cn*uvar(i,1)**(cn- one)
281 ELSEIF(uvar(i,1)>zero. and .cn<one)
THEN
282 qh1= cb*cn*uvar(i,1)**(one - cn)
286 IF(uvar(i,1)<=zero. or .uvar(i,2)<=eps0)
THEN
289 qh2=cd*cm*uvar(i,1)**(cm- one)*log(uvar(i,2)/eps0)
291 qh2=cd*cm*uvar(i,1)**(one - cm)*log(uvar(i,2)/eps0)
299 svm = half*(signxx(i)*signxx(i)
300 . +signyy(i)*signyy(i)
301 . +signzz(i)*signzz(i)
302 . +signxy(i)*signxy(i)
303 . +signyz(i)*signyz(i)
304 . +signzx(i)*signzx(i))
305 svm = sqrt(three * svm)
306 r =
min(one,uvar(i,3)/
max(em20,svm))
307 signxx(i)=signxx(i)*r
308 signyy(i)=signyy(i)*r
309 signzz(i)=signzz(i)*r
310 signxy(i)=signxy(i)*r
311 signyz(i)=signyz(i)*r
312 signzx(i)=signzx(i)*r
314 dpla_i = off(i)*svm*umr/(g3+uvar(i,4))
315 uvar(i,1) = uvar(i,1) + dpla_i
324 signxx(i) = (signxx(i) - press) * off(i)
325 signyy(i) = (signyy(i) - press) * off(i)
326 signzz(i) = (signzz(i) - press) * off(i)
327 signxy(i) = signxy(i) * off(i)
328 signyz(i) = signyz(i) * off(i)
329 signzx(i) = signzx(i) * off(i)
subroutine sigeps45(nel0, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, rho0, rho, volume, eint, 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, sigy, defp, amu)