42
43
44
45
46
47
48#include "implicit_f.inc"
49#include "comlock.inc"
50
51
52
53#include "mvsiz_p.inc"
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93#include "param_c.inc"
94#include "scr05_c.inc"
95
96
97 INTEGER ,INTENT(IN) :: ,NUPARAM,NUVAR
98 INTEGER ,INTENT(IN) :: JTHE
99 INTEGER NGL(NEL),MAT(NEL),IPM(NPROPMI,*)
101 . time,timestep,uparam(nuparam),
102 . rho(nel),rho0(nel),volume(nel),eint(nel),
103 . epspxx(nel),epspyy(nel),epspzz(nel),
104 . epspxy(nel),epspyz(nel),epspzx(nel),
105 . depsxx(nel),depsyy(nel),depszz(nel),
106 . depsxy(nel),depsyz(nel),depszx(nel),
107 . epsxx(nel) ,epsyy(nel) ,epszz(nel) ,
108 . epsxy(nel) ,epsyz(nel) ,epszx(nel) ,
109 . sigoxx(nel),sigoyy(nel),sigozz(nel),
110 . sigoxy(nel),sigoyz(nel),sigozx(nel)
111 my_real,
DIMENSION(NEL) ,
INTENT(INOUT)
112
113
114
116 . signxx(nel),signyy(nel),signzz(nel),
117 . signxy(nel),signyz(nel),signzx(nel),
118 . sigvxx(nel),sigvyy(nel),sigvzz(nel),
119 . sigvxy(nel),sigvyz(nel),sigvzx(nel),
120 . soundsp(nel),etse(nel),viscmax(nel)
121
122
123
124 my_real,
DIMENSION(NEL) ,
INTENT(INOUT) :: off
125 my_real,
DIMENSION(NEL,NUVAR) ,
INTENT(INOUT) ::uvar
126
127
128
129 INTEGER I,J,J1,J2,I1,I2,KK,IADBUF,EFLAG,ISMSTR,IFUNC(100)
131 . e,emart,nu,g,g2,wave,sqdt,a,b,c,fct,fctp, dftr,unmxn,db,
132 .
alpha,epsl,aa1,fm,dfmsa,dfmas,uxx,nn,beta,gm,km,
133 . cb,cc,caas,cbas,pold,dgt,dkt,cp,tini,inve, n1,n2,n3,
134 . k,p,sxx,cas,csa,tsas,tfas, tssa,tfsa,rv_pui,dfs,
135 . syy,szz,sxy,syz,szx,fass,fsas,fasf,fsaf,rsas,rfas,
136 . sv,fs,fs0,yld_ass,yld_asf,yld_sas,yld_saf,rssa,rfsa,
137 . dfm, fss,dsxx,dsyy,dsxy,dsyz,dszx,dszz,var,h,et,
138 . pm,delta,x1,x2,test,test2,ftest,gnew,knew,betan,
139 . nx2,ny2 ,nz2,nxy2,nyz2,nzx2,ne,dnx,dny,dnz,dnxy,dnyz,dnzx,
140 . nxx(nel),nyy(nel),nzz(nel),nxy(nel),nyz(nel),nzx(nel),
141 . e1(nel),e2(nel),e3(nel),e4(nel),e5(nel),e6(nel),trde(nel),
142 . de1(nel),de2(nel),de3(nel),gt(nel),kt(nel),
143 . ee1(nel),ee2(nel),ee3(nel),pp(nel),nne(nel),det(nel),
144 . sigxx(nel), sigyy(nel), sigzz(nel)
146 . ev(mvsiz,3)
148 . av(mvsiz,6),evv(mvsiz,3),dirprv(mvsiz,3,3)
149
150
151
152
153
154 e = uparam(1)
155 nu = uparam(2)
156 g = uparam(3)
157 k = uparam(4)
158 aa1 = uparam(5)
159 yld_ass = uparam(6)
160 yld_asf = uparam(7)
161 yld_sas = uparam(8)
162 yld_saf = uparam(9)
164 epsl = uparam(11)/(sqrt(two_third)+
alpha)
165 emart = uparam(14)
166 eflag = int(uparam(15))
167 gm = uparam(16)
168 km = uparam(17)
169 g2 = two*g
171 sqdt = sqrt(two/three)
172 cas = uparam(18)
173 csa = uparam(19)
174 tsas = uparam(20)
175 tfas = uparam(21)
176 tssa = uparam(22)
177 tfsa = uparam(23)
178 cp = uparam(24)
179 tini = uparam(25)
180
181
182 DO i=1,nel
183 av(i,1)=epsxx(i)
184 av(i,2)=epsyy(i)
185 av(i,3)=epszz(i)
186 av(i,4)=epsxy(i) * half
187 av(i,5)=epsyz(i) * half
188 av(i,6)=epszx(i) * half
189 ENDDO
190
191
192 IF (iresp==1) THEN
194 ELSE
196 ENDIF
197
198 IF(ismstr== 0.OR. ismstr==2.OR. ismstr==4) THEN
199 DO i=1,nel
200
201 ev(i,1)=exp(evv(i,1))
202 ev(i,2)=exp(evv(i,2))
203 ev(i,3)=exp(evv(i,3))
204 ENDDO
205 ELSEIF(ismstr==10.OR.ismstr==12) THEN
206 DO i =1,nel
207 ev(i,1)=sqrt(evv(i,1)+ one)
208 ev(i,2)=sqrt(evv(i,2)+ one)
209 ev(i,3)=sqrt(evv(i,3)+ one)
210 ENDDO
211 ELSE
212
213 DO i=1,nel
214 ev(i,1)=evv(i,1)+ one
215 ev(i,2)=evv(i,2)+ one
216 ev(i,3)=evv(i,3)+ one
217 ENDDO
218 ENDIF
219 DO i =1,nel
220
221 det(i) =ev(i,1)*ev(i,2)*ev(i,3)
222 IF(det(i)/=zero) THEN
223 trde(i) = log(det(i))
224 rv_pui = exp((-third)*trde(i))
225 ELSE
226 rv_pui = zero
227 trde(i)= zero
228 ENDIF
229 ee1(i) = log(ev(i,1)*rv_pui)
230 ee2(i) = log(ev(i,2)*rv_pui)
231 ee3(i) = log(ev(i,3)*rv_pui)
232 ENDDO
233
234
235
236 IF (jthe == 0) THEN
237 DO i=1,nel
238 temp(i) = tini + eint(i) / rho0(i)/cp/
max(em15,volume(i))
239 ENDDO
240 ENDIF
241
242
243 IF (eflag > zero)THEN
244
245 rsas = yld_ass *(sqdt+
alpha)-cas*tsas
246 rfas = yld_asf *(sqdt+
alpha)-cas*tfas
247 rssa = yld_sas *(sqdt+
alpha)-csa*tssa
248 rfsa = yld_saf *(sqdt+
alpha)-csa*tfsa
249 DO i = 1,nel
250
251 fm = uvar(i,1)
252 gt(i) = g + fm * (gm - g)
253 kt(i) = k + fm * (km - k)
254
255 p = kt(i) * (trde(i) - three*
alpha*epsl*fm)
256
257 ne = sqrt( ee1(i)**2 + ee2(i)**2 + ee3(i)**2)
258 nxx(i) =ee1(i)/
max(ne,em20)
259 nyy(i) =ee2(i)/
max(ne,em20)
260 nzz(i) =ee3(i)/
max(ne,em20)
261
262 sxx= two*gt(i)*(ee1(i) -epsl*fm*nxx(i))
263 syy= two*gt(i)*(ee2(i) -epsl*fm*nyy(i))
264 szz= two*gt(i)*(ee3(i) -epsl*fm*nzz(i))
265
266 sv = sqrt( sxx*sxx + syy*syy + szz*szz )
267
268 soundsp(i) = sqrt(aa1/rho0(i))
269
270 viscmax(i) = zero
271
272
273
274 dfmsa = zero
275 dfmas = zero
276
277 fs = sv + three*
alpha*p - cas*temp(i)
278
279
280 fass = fs - rsas
281 fasf = fs - rfas
282 fs0 = uvar(i,2)
283 beta = epsl*(two*gt(i)+nine*kt(i)*
alpha*
alpha)
284 IF((fs - fs0) > zero .AND. fass > zero.AND. fasf < zero .AND. fm < one )THEN
285
286 db = (two * (gm-g) +nine*
alpha*
alpha*(km-k)) *epsl
287 unmxn = one - fm
288 dftr = two*ne*(gm-g) + three*
alpha*trde(i)*(km-k)
289 dfmas =
min(one, -(fs-fs0)*(one-fm)/(fasf-beta*(one-fm) ) )
290 a = unmxn *db
291 b = (rfas-fs+unmxn*(beta-dftr))
292 c = unmxn*(fs0 - fs)
293 DO kk = 1,3
294 fct = dfmas*dfmas *a+ dfmas* b +c
295 fctp = two*dfmas *a+ b
296 dfmas = dfmas - fct / fctp
297 ENDDO
298 dfmas =
min(one,dfmas )
299 ENDIF
300
301
302
303 fs = sv + three*
alpha*p - csa*temp(i)
304 fsas = fs -rssa
305 fsaf = fs -rfsa
306 fs0 = uvar(i,3)
307
308 IF((fs - fs0) < zero .AND. fsas < zero.AND. fsaf > zero )THEN
309
310 db = (two * (gm-g) +nine*
alpha*
alpha*(km-k))*epsl
311 dftr = two * (gm-g)*ne+ three*
alpha*(km-k)*trde(i)
312 dfmsa = zero
313 a = fm *db
314 b = -(rfsa-fs+fm*(dftr-beta))
315 c = -fm*(fs - fs0)
316 DO kk = 1,3
317 fct = dfmsa*dfmsa *a+ dfmsa* b +c
318 fctp = two*dfmsa *a+ b
319 dfmsa = dfmsa - fct / fctp
320 ENDDO
321 dfmsa =
max(-one , dfmsa )
322 ENDIF
323
324
325 dfm = dfmas + dfmsa
326 IF(dfm < zero .AND. fm == zero) dfm = zero
327
328
329
330
331 dgt = dfm * (gm - g)
332 dkt = dfm * (km - k)
333
334 sxx = sxx -two*gt(i) * epsl*nxx(i)*dfm
335 . + two*dgt* (ee1(i)-epsl*nxx(i)*dfm)
336 syy = syy -two*gt(i)* epsl*dfm*nyy(i)
337 . + two*dgt* (ee2(i)-epsl*nyy(i)*dfm)
338 szz = szz -two*gt(i)* epsl*dfm*nzz(i)
339 . + two*dgt* (ee3(i)-epsl*nzz(i)*dfm)
340
341 p = p - kt(i)*epsl*three*
alpha*dfm
342 . + dkt *(trde(i) -epsl*three*
alpha*dfm)
343
344 sigxx(i)= sxx + p
345 sigyy(i)= syy + p
346 sigzz(i)= szz + p
347 sv = sqrt( sxx*sxx + syy*syy + szz*szz )
348 fs = sv + three*
alpha*p
349
350 IF(det(i)/=zero)THEN
351 inve = one/det(i)
352 ELSE
353 inve = zero
354 ENDIF
355 sigxx(i)= sigxx(i) *inve
356 sigyy(i)= sigyy(i) *inve
357 sigzz(i)= sigzz(i) *inve
358
359 signxx(i) = dirprv(i,1,1)*dirprv(i,1,1)*sigxx(i)
360 . + dirprv(i,1,2)*dirprv(i,1,2)*sigyy(i)
361 . + dirprv(i,1,3)*dirprv(i,1,3)*sigzz(i)
362 signyy(i) = dirprv(i,2,2)*dirprv(i,2,2)*sigyy(i)
363 . + dirprv(i,2,3)*dirprv(i,2,3)*sigzz(i)
364 . + dirprv(i,2,1)*dirprv(i,2,1)*sigxx(i)
365 signzz(i) = dirprv(i,3,3)*dirprv(i,3,3)*sigzz(i)
366 . + dirprv(i,3,1)*dirprv(i,3,1)*sigxx(i)
367 . + dirprv(i,3,2)*dirprv(i,3,2)*sigyy(i)
368 signxy(i) = dirprv(i,1,1)*dirprv(i,2,1)*sigxx(i)
369 . + dirprv(i,1,2)*dirprv(i,2,2)*sigyy(i)
370 . + dirprv(i,1,3)*dirprv(i,2,3)*sigzz(i)
371 signyz(i) = dirprv(i,2,2)*dirprv(i,3,2)*sigyy(i)
372 . + dirprv(i,2,3)*dirprv(i,3,3)*sigzz(i)
373 . + dirprv(i,2,1)*dirprv(i,3,1)*sigxx(i)
374 signzx(i) = dirprv(i,3,3)*dirprv(i,1,3)*sigzz(i)
375 . + dirprv(i,3,1)*dirprv(i,1,1)*sigxx(i)
376 . + dirprv(i,3,2)*dirprv(i,1,2)*sigyy(i)
377
378 uvar(i,1) = uvar(i,1) + dfm
379 uvar(i,1) =
max(zero, uvar(i,1))
380 uvar(i,1) =
min(one, uvar(i,1))
381 uvar(i,2) = fs- cas*temp(i)
382 uvar(i,3) = fs- csa*temp(i)
383 IF (dfmas /= zero) dfs = abs(uvar(i,2) - fs0)
384 IF (dfmsa /= zero) dfs = abs(uvar(i,3) - fs0)
385 IF (dfs /= zero .AND. epsl /= zero .AND. dfm/= zero) THEN
386 h = dfs/epsl/dfm
387 etse(i) = h *(one+nu) /( (e + uvar(i,1)*(emart-e)) + h)
388 ELSE
389 etse(i) = one
390 ENDIF
391 uvar(i,10) = epsxx(i)
392 ENDDO
393
394 ELSE
395 rsas = yld_ass *(sqdt+
alpha)-cas*tsas
396 rfas = yld_asf *(sqdt+
alpha)-cas*tfas
397 rssa = yld_sas *(sqdt+
alpha)-csa*tssa
398 rfsa = yld_saf *(sqdt+
alpha)-csa*tfsa
399 DO i = 1,nel
400
401 fm = uvar(i,1)
402
403 p = k * (trde(i) - three*
alpha*epsl*fm)
404
405 ne = sqrt( ee1(i)**2 + ee2(i)**2 + ee3(i)**2)
406 nxx(i) =ee1(i)/
max(ne,em20)
407 nyy(i) =ee2(i)/
max(ne,em20)
408 nzz(i) =ee3(i)/
max(ne,em20)
409
410 sxx= g2*(ee1(i) -epsl*fm*nxx(i))
411 syy= g2*(ee2(i) -epsl*fm*nyy(i))
412 szz= g2*(ee3(i) -epsl*fm*nzz(i))
413
414 sv = sqrt( sxx*sxx + syy*syy + szz*szz )
415
416 soundsp(i) = sqrt(aa1/rho0(i))
417 viscmax(i) = zero
418
419
420
421 dfmsa = zero
422 dfmas = zero
423
424
425 fs = sv + three*
alpha*p - cas*temp(i)
426
427
428
429 fass = fs -rsas
430 fasf = fs -rfas
431 fs0 = uvar(i,2)
432 IF((fs - fs0) > zero .AND. fass > zero.AND. fasf < zero.AND. fm < one ) then
433 dfmas =
min(one, -(fs-fs0)*(one-fm)/(fasf-beta*(one-fm) ) )
434 ENDIF
435
436
437 fs = sv + three*
alpha*p - csa*temp(i)
438 fsas = fs - rssa
439 fsaf = fs - rfsa
440 fs0 = uvar(i,3)
441 IF((fs - fs0) < zero .AND. fsas < zero .AND. fsaf > zero .AND. fm > zero ) THEN
442 dfmsa =
max(-one , fm * (fs - fs0)/ (fsaf+beta*fm) )
443 ENDIF
444
445
446 dfm = dfmas + dfmsa
447 IF(dfm < zero .AND. fm == zero) dfm = zero
448
449
450 sxx = sxx -g2* epsl*dfm*nxx(i)
451 syy = syy -g2* epsl*dfm*nyy(i)
452 szz = szz -g2* epsl*dfm*nzz(i)
453 p = p - k*epsl*three*
alpha*dfm
454
455 sigxx(i)= sxx + p
456 sigyy(i)= syy + p
457 sigzz(i)= szz + p
458 sv = sqrt( sxx*sxx + syy*syy + szz*szz )
459 fs = sv + three*
alpha*p
460
461 IF(det(i)/=zero)THEN
462 inve = one/det(i)
463 ELSE
464 inve = zero
465 ENDIF
466 sigxx(i)= sigxx(i) *inve
467 sigyy(i)= sigyy(i) *inve
468 sigzz(i)= sigzz(i) *inve
469
470 signxx(i) = dirprv(i,1,1)*dirprv(i,1,1)*sigxx(i)
471 . + dirprv(i,1,2)*dirprv(i,1,2)*sigyy(i)
472 . + dirprv(i,1,3)*dirprv(i,1,3)*sigzz(i)
473 signyy(i) = dirprv(i,2,2)*dirprv(i,2,2)*sigyy(i)
474 . + dirprv(i,2,3)*dirprv(i,2,3)*sigzz(i)
475 . + dirprv(i,2,1)*dirprv(i,2,1)*sigxx(i)
476 signzz(i) = dirprv(i,3,3)*dirprv(i,3,3)*sigzz(i)
477 . + dirprv(i,3,1)*dirprv(i,3,1)*sigxx(i)
478 . + dirprv(i,3,2)*dirprv(i,3,2)*sigyy(i)
479 signxy(i) = dirprv(i,1,1)*dirprv(i,2,1)*sigxx(i)
480 . + dirprv(i,1,2)*dirprv(i,2,2)*sigyy(i
481 . + dirprv(i,1,3)*dirprv(i,2,3)*sigzz(i)
482 signyz(i) = dirprv(i,2,2)*dirprv(i,3,2)*sigyy(i)
483 . + dirprv(i,2,3)*dirprv(i,3,3)*sigzz(i)
484 . + dirprv(i,2,1)*dirprv(i,3,1)*sigxx(i)
485 signzx(i) = dirprv(i,3,3)*dirprv(i,1,3)*sigzz(i)
486 . + dirprv(i,3,1)*dirprv(i,1,1)*sigxx(i)
487 . + dirprv(i,3,2)*dirprv(i,1,2)*sigyy(i)
488 uvar(i,1) = uvar(i,1) + dfm
489 uvar(i,1) =
max(zero, uvar(i,1))
490 uvar(i,1) =
min(one, uvar(i,1))
491 uvar(i,2) = fs- cas*temp(i)
492 uvar(i,3) = fs- csa*temp(i)
493 dfs = zero
494 IF (dfmas /= zero) dfs = abs(uvar(i,2) - fs0)
495 IF (dfmsa /= zero) dfs = abs(uvar(i,3) - fs0)
496 IF (dfs /= zero .AND. epsl /= zero.AND.dfm/= zero) THEN
497 h = dfs/epsl/dfm
498 etse(i) = h * (one+nu)/(e + h)
499 ELSE
500 etse(i) = one
501 ENDIF
502 uvar(i,4) = epsl*dfm
503 uvar(i,7) = uvar(i,7)+epsl*dfm
504 uvar(i,10) = epsxx(i)
505 ENDDO
506
507 ENDIF
508 RETURN
subroutine valpvecdp_v(sig, val, vec, nel)
subroutine valpvec_v(sig, val, vec, nel)