33 1 NEL ,NUPARAM ,NUVAR ,NFUNC ,IFUNC ,
34 2 NPF ,TF ,TIME ,TIMESTEP ,UPARAM ,
35 3 AREA ,THKLY ,SOUNDSP ,VISCMAX ,UVAR ,
36 4 DEPSXX ,DEPSYY ,DEPSXY ,DEPSYZ ,DEPSZX ,
37 5 EPSXX ,EPSYY ,EPSXY ,EPSYZ ,EPSZX ,
38 6 SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
39 7 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
40 8 SIGVXX ,SIGVYY ,SIGVXY ,TAN_PHI ,OFFGG ,
41 9 RHO0 ,ETSE ,SHF ,ALDT ,NSENSOR ,
42 A SENSOR_TAB,NIPARAM ,IPARAM)
50#include "implicit_f.inc"
94#include "tabsiz_c.inc"
98 INTEGER ,
INTENT(IN) :: NEL
99 INTEGER ,
INTENT(IN) :: NFUNC
100 INTEGER ,
INTENT(IN) :: NUPARAM
101 INTEGER ,
INTENT(IN) :: NIPARAM
102 INTEGER ,
INTENT(IN) :: NUVAR
103 INTEGER ,
INTENT(IN) :: NSENSOR
104 my_real ,
INTENT(IN) :: TIME,TIMESTEP
105 INTEGER ,
DIMENSION(NFUNC) ,
INTENT(IN) :: IFUNC
106 INTEGER ,
DIMENSION(SNPC) ,
INTENT(IN) :: NPF
107 my_real ,
DIMENSION(STF) ,
INTENT(IN) :: TF
108 INTEGER ,
DIMENSION(NIPARAM) ,
INTENT(IN) :: IPARAM
109 my_real ,
DIMENSION(NUPARAM) ,
INTENT(IN) :: uparam
110 my_real ,
DIMENSION(NEL) ,
INTENT(IN) ::
area,rho0,thkly,offgg,
111 . depsxx,depsyy,depsxy,depsyz,depszx,epsxx,epsyy,epsxy,epsyz,epszx,
112 . sigoxx,sigoyy,sigoxy,sigoyz,sigozx,shf,aldt
113 TYPE () ,
DIMENSION(NSENSOR) ,
INTENT(IN) :: SENSOR_TAB
117 my_real ,
DIMENSION(NEL) ,
INTENT(OUT) :: tan_phi,etse,
118 . signxx,signyy,signxy,signyz,signzx,
119 . sigvxx,sigvyy,sigvxy,soundsp,viscmax
123 my_real ,
DIMENSION(NEL,NUVAR) ,
INTENT(INOUT) :: uvar
127 INTEGER :: I,ITER,NITER,ISENS,IFUNCC,IFUNCT,IFUNCS,IFUNF1,IFUNF2
128 INTEGER ,
DIMENSION(NEL) :: IADF1,IADF2,IAD1,IAD2,,ILEN1,ILEN2,ILEN3,
129 . ILENF1,ILENF2,IPOS,IPOS1,IPOS2
130 my_real :: KMAX,GMAX,GFROT,GSH,STIFF
131,HC0,HT0,VISCE,VISCG,LC0,,HDC,HDT,
132 . TFROT,SIGG,EC2,ET2,SINN,TAN2,DAMP,V1,V2,DTINV,,ETT,TRACE,
133 . zerostress,dsig,lmin,tstart,phi1,phi2,j11,j12,j21,j22,det
134 my_real ,
DIMENSION(3) :: yfac
135 my_real ,
DIMENSION(NEL) :: sigc,sigt,ec,et,lc,lt,fn,angle,kg,gxy,
136 . dtang,tfold,cvisc,cvist,flexc,flext,flexf,hc,ht,dc,dt,dcc,dtt,
137 . fc,ft,fpc,fpt,dflxc,dflxt,epsf,epsfc,epsft,xc,xt,beta,ph01,ph02
171 zerostress = uparam(8)
189 IF (zerostress > zero .and. isens > 0)
THEN
190 tstart = sensor_tab(isens)%TSTART
194 dtinv = timestep/
max(timestep**2,em20)
199#include "vectorize.inc"
202 mass = rho0(i)*
area(i)*thkly(i)*fourth
203 cvisc(i) = sqrt(mass*kflex1)*dtinv*third
204 cvist(i) = sqrt(mass*kflex2)*dtinv*third
208#include "vectorize.inc"
210 etc = uvar(i,4) + depsxx(i)
211 ett = uvar(i,5) + depsyy(i)
214 ec(i) = exp(etc) - one
215 et(i) = exp(ett) - one
216 lc(i) = lc0 * (one + ec(i))
217 lt(i) = lt0 * (one + et(i))
226#include "vectorize.inc"
229 iad1(i) = npf(ifuncc) / 2 + 1
230 iad2(i) = npf(ifunct) / 2 + 1
231 iad3(i) = npf(ifuncs) / 2 + 1
232 ilen1(i) = npf(ifuncc+1) / 2 - iad1(i) - ipos(i)
233 ilen2(i) = npf(ifunct+1) / 2 - iad2(i) - ipos(i)
234 ilen3(i) = npf(ifuncs+1) / 2 - iad3(i) - ipos(i)
238 iadf1(i) = npf(ifunf1) / 2 + 1
239 ilenf1(i) = npf(ifunf1+1) / 2 - iadf1(i) - ipos(i)
244 iadf2(i) = npf(ifunf2) / 2 + 1
245 ilenf2(i) = npf(ifunf2+1) / 2 - iadf2(i) - ipos(i)
251 epsfc(1:nel) = uvar(1:nel,7)
252 epsft(1:nel) = uvar(1:nel,8)
256 epsf(1:nel) =(hc0 * epsfc(1:nel) + ht0 * epsft(1:nel)) / h0
257 flexf(1:nel) = kflex * epsf(1:nel)
261 hc(i) = hc0 * (epsfc(i) + one)
262 ht(i) = ht0 * (epsft(i) + one)
263 dc(i) = sqrt(lc(i)**2 + hc(i)**2)
264 dt(i) = sqrt(lt(i)**2 + ht(i)**2)
271 CALL vinter2(tf,iad1,ipos1,ilen1,nel,dcc,fpc,fc)
272 CALL vinter2(tf,iad2,ipos2,ilen2,nel,dtt,fpt,ft)
278 flexc(1:nel) = flexc(1:nel) * kflex1
279 dflxc(1:nel) = dflxc(1:nel) * kflex1
282 flexc(i) = kflex1 * log(epsfc(i) + one)
289 CALL vinter2(tf,iadf2,ipos,ilenf2,nel,epsft,dflxt,flext)
290 flext(1:nel) = flext(1:nel) * kflex2
291 dflxt(1:nel) = dflxt(1:nel) * kflex2
294 flext(i) = kflex2 * log(epsft(i) + one)
295 dflxt(i) = kflex2 / (epsft(i) + one)
299#include "vectorize.inc"
301 dyc = epsfc(i) - uvar(i,7)
302 dyt = epsft(i) - uvar(i,8)
305 phi1 = fc(i) * hdc + flexc(i) + flexf(i) + cvisc(i)*dyc
306 phi2 = ft(i) * hdt + flext(i) + flexf(i) + cvist(i)*dyt
307 j12 = kflex * ht0 / h0
308 j21 = kflex * hc0 / h0
309 j11 = j12 + fpc(i)*hdc*hc0 + dflxc(i) + cvisc(i)
310 j22 = j21 + fpt(i)*hdt*ht0 + dflxt(i) + cvist(i)
311 det = j11 * j22 - j12 * j21
313 epsfc(i) = epsfc(i) - beta(i) * (j22 * phi1 - j12 * phi2) / det
314 epsft(i) = epsft(i) + beta(i) * (j12 * phi1 - j11 * phi2) / det
316 epsfc(i) =
max(epsfc(i), em04 - one)
317 epsft(i) =
max(epsft(i), em04 - one)
319 IF (abs(phi1) > ph01(i) .and. abs(phi2) > ph02(i))
THEN
322 beta(i) = beta(i) * half
323 beta(i) =
max(beta(i), em02)
331#include "vectorize.inc"
333 sigc(i) = fc(i) * lc(i) / dc(i)
334 sigt(i) = ft(i) * lt(i) / dt(i)
351#include "vectorize.inc"
353 trace = exp(epsxx(i) + epsyy(i))
354 ec2 =
max(trace / (ec(i) + one), em6)
355 et2 =
max(trace / (et(i) + one), em6)
357 signxx(i) = sigc(i) / ec2
358 signyy(i) = sigt(i) / et2
363#include "vectorize.inc"
365 tan_phi(i)= depsxy(i)
366 angle(i) = atan(tan_phi(i))*hundred80/pi
370 CALL vinter2(tf,iad3,ipos,ilen3,nel,angle,gxy,signxy)
372#include "vectorize.inc"
375 kg(i) = gxy(i) * tan2 * yfac(3)
377 damp = sqrt(rho0(i)*
area(i)*thkly(i)*half)
378 v1 = visce*damp*sqrt(kmax)
379 v2 = visce*damp*sqrt(kmax)
380 sigvxx(i) = dtinv*(depsxx(i))*v1
381 sigvyy(i) = dtinv*(depsyy(i))*v2
383 IF (fn(i) > zero)
THEN
384 tfrot = two_third*viscg*fn(i)*(hc0+ht0)/(lc(i)+lt(i))
385 dtang(i) = depsxy(i) - tan_phi(i)
386 sigg = tfold(i) + gfrot*dtang(i)
387 IF (abs(sigg) > tfrot)
THEN
388 sigvxy(i) = sign(tfrot,sigg)
394 sinn = tan_phi(i) / sqrt(one + tan2)
395 stiff = kmax*(one+sinn) + gmax
396 lmin =
min(dc(i)/dc0,dt(i)/dt0)*uvar(i,14)
397 soundsp(i) = sqrt(stiff/(rho0(i)))*aldt(i)/lmin
398 viscmax(i) =
max(v1,v2)
402#include "vectorize.inc"
404 uvar(i,3) = signxy(i) + sigvxy(i)
405 tan_phi(i)= depsxy(i)
406 uvar(i,6) = depsxy(i)
407 uvar(i,9) = sigvxy(i)
409 signyz(i) = sigoyz(i) + gsh * depsyz(i)
410 signzx(i) = sigozx(i) + gsh * depszx(i)
415 IF (zerostress > zero)
THEN
416 IF (time <= tstart)
THEN
417#include "vectorize.inc"
419 uvar(i,11) = signxx(i)
420 uvar(i,12) = signyy(i)
421 uvar(i,13) = signxy(i)
427#include "vectorize.inc"
429 dsig = signxx(i) - sigoxx(i) - uvar(i,11)
430 IF((uvar(i,11) > zero).AND.(dsig < zero))
THEN
431 uvar(i,11) =
max(zero,uvar(i,11)+zerostress*dsig)
432 ELSEIF((uvar(i,11) < zero).AND.(dsig > zero))
THEN
433 uvar(i,11) =
min(zero,uvar(i,11)+zerostress*dsig)
435 dsig = signyy(i) - sigoyy(i) - uvar(i,12)
436 IF((uvar(i,12) > zero).AND.(dsig < zero))
THEN
437 uvar(i,12) =
max(zero,uvar(i,12)+zerostress*dsig)
438 ELSEIF((uvar(i,12) < zero).AND.(dsig > zero))
THEN
439 uvar(i,12) =
min(zero,uvar(i,12)+zerostress*dsig)
441 dsig = signxy(i) - sigoxy(i) - uvar(i,13)
442 IF((uvar(i,13) > zero).AND.(dsig < zero))
THEN
443 uvar(i,13) =
max(zero,uvar(i,13)+zerostress*dsig)
444 ELSEIF((uvar(i,13) < zero).AND.(dsig > zero))
THEN
445 uvar(i,13) =
min(zero,uvar(i,13)+zerostress*dsig)
447 signxx(i) = signxx(i) - uvar(i,11)
448 signyy(i) = signyy(i) - uvar(i,12)
449 signxy(i) = signxy(i) - uvar(i,13)
455 soundsp(i) = sqrt(kmax*two/(rho0(i)))
subroutine sigeps158c(nel, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, area, thkly, soundsp, viscmax, uvar, depsxx, depsyy, depsxy, depsyz, depszx, epsxx, epsyy, epsxy, epsyz, epszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, sigvxx, sigvyy, sigvxy, tan_phi, offgg, rho0, etse, shf, aldt, nsensor, sensor_tab, niparam, iparam)