45
46
47
48#include "implicit_f.inc"
49
50
51
52#include "mvsiz_p.inc"
53
54
55
56#include "scr05_c.inc"
57#include "impl1_c.inc"
58
59
60
61 INTEGER NEL, NUPARAM, NUVAR,ISMSTR,IHET,IEXPAN
62
64 . time , timestep , uparam(nuparam),
65 . rho(nel), rho0(nel), volume(nel), eint(nel),
66 . epspxx(nel), epspyy(nel), epspzz(nel),
67 . epspxy(nel), epspyz(nel), epspzx(nel),
68 . depsxx(nel), depsyy(nel), depszz(nel),
69 . depsxy(nel), depsyz(nel), depszx(nel),
70 . epsxx(nel), epsyy(nel), epszz(nel),
71 . epsxy(nel), epsyz(nel), epszx(nel),
72 . sigoxx(nel), sigoyy(nel), sigozz(nel),
73 . sigoxy(nel), sigoyz(nel), sigozx(nel), offg(nel),
74 . epsth3(nel)
75
76
77
79 . signxx(nel), signyy(nel), signzz(nel),
80 . signxy(nel), signyz(nel), signzx(nel),
81 . sigvxx(nel), sigvyy(nel), sigvzz(nel),
82 . sigvxy(nel), sigvyz(nel), sigvzx(nel),
83 . soundsp(nel), viscmax(nel), et(nel)
84
85
86
88 . uvar(nel,nuvar), off(nel)
89
90
91
92 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
93 my_real finter,fintte,tf(*),fint2v
94 EXTERNAL finter,fintte
95
96
97
98 INTEGER I,J,K,II,NORDRE,NPRONY,ITYPE,IVISC
100 . mu(100),al(100),taux(100), gama(100),gamainf,nu,
101 . gmax,s(mvsiz,3),sd(mvsiz,3),
102 . ec(mvsiz,3),hd(mvsiz,3,100),
103 . ev(mvsiz,3),av(mvsiz,6),dirprv(mvsiz,3,3),
104 . rv(mvsiz),p(mvsiz),rbulk,
105 . evv(mvsiz,3),ssp(mvsiz,3),fac, beta(100),
106 . dtinv,visc,fac1,
107 . e1,e2,e3,e4,e5,e6,epsp, vm,dav,visc1,sigxx(mvsiz),
108 . sigyy(mvsiz),sigzz(mvsiz),sigxy(mvsiz),sigzx(mvsiz),
109 . sigyz(mvsiz),dvisc,a11,a12,a13,a22,a21,a23,
110 . a31,a32,a33,sdg0(mvsiz,6),sdg(mvsiz,6),hg(mvsiz,6),hg0(6),h(mvsiz,3,100)
111
113 . a(3,3),dpra(3,3),eigen(3),amax,eti(mvsiz),efac,rv_pui,rv_pui_tab(mvsiz),
114 . cj(mvsiz),lam_al(3),rvl,lamin,cimax(mvsiz),gvis,cmax,cmax0,
115 . d11,d12,d13,d21,d22,d23,d31,d32,d33,mu2
116 my_real ,
DIMENSION(NEL,3) :: ai,bi
117
118
119
120 gmax = zero
121 nu = uparam(1)
122 nordre = nint(uparam(2))
123 nprony = nint(uparam(3))
124 gamainf = uparam(4)
125 rbulk = uparam(5)
126 dvisc = uparam(6)
127
128 DO i= 1,nordre
129 mu(i) = uparam(6 + i)
130 al(i) = uparam(6 + nordre + i)
131 beta(i) = uparam(2*nordre + 2*nprony + 7 + i)
132
133 gmax = gmax + mu(i)
134 ENDDO
135
136 ivisc = uparam(2*nordre + 2*nprony + 7 )
137 IF(nprony /= zero)THEN
138 DO i=1,nprony
139 gama(i) = uparam(6 + 2*nordre + i)
140 taux(i) = uparam(6 + 2*nordre + nprony + i)
141 ENDDO
142 ENDIF
143
144
145 gmax = two*gmax
146
147 IF(time==zero)THEN
148 DO j = 1,nuvar
149 DO i = 1, nel
150 uvar(i,j) = zero
151 ENDDO
152 ENDDO
153 ENDIF
154 DO i=1,nel
155 av(i,1)=epsxx(i)
156 av(i,2)=epsyy(i)
157 av(i,3)=epszz(i)
158 av(i,4)=epsxy(i)/2
159 av(i,5)=epsyz(i)/2
160 av(i,6)=epszx(i)/2
161 ENDDO
162
163
164 IF (iresp==1) THEN
166 ELSE
168 ENDIF
169
170
171
172
173 IF(ismstr==0.OR.ismstr==2.OR.ismstr==4) THEN
174 DO i=1,nel
175
176 ev(i,1)=exp(evv(i,1))
177 ev(i,2)=exp(evv(i,2))
178 ev(i,3)=exp(evv(i,3))
179 ENDDO
180 ELSEIF(ismstr==10.OR.ismstr==12) THEN
181 DO i =1,nel
182 IF(offg(i)<=one) THEN
183 ev(i,1)=sqrt(evv(i,1)+ one)
184 ev(i,2)=sqrt(evv(i,2)+ one)
185 ev(i,3)=sqrt(evv(i,3)+ one)
186 ELSE
187 ev(i,1)=evv(i,1)+ one
188 ev(i,2)=evv(i,2)+ one
189 ev(i,3)=evv(i,3)+ one
190 END IF
191 ENDDO
192 ELSE
193
194 DO i=1,nel
195 ev(i,1)=evv(i,1)+ one
196 ev(i,2)=evv(i,2)+ one
197 ev(i,3)=evv(i,3)+ one
198 ENDDO
199 ENDIF
200
201 IF (impl_s > 0 .OR. ihet > 1) THEN
202 DO j =1,3
203 eti(1:nel)=zero
204 DO k= 1,nordre
205
206 DO i=1,nel
207 amax = ev(i,j)
208 IF(amax>zero) THEN
209 IF(mu(k)/=zero) THEN
210 IF((al(k)-one)/=zero) THEN
211 eti(i) = eti(i) + mu(k)*exp((al(k)-one)*log(amax))
212 ELSE
213 eti(i) = eti(i) + mu(k)
214 ENDIF
215 ELSE
216 eti(i) = eti(i) + zero
217 ENDIF
218 ELSE
219 eti(i) = zero
220 ENDIF
221 ENDDO
222 ENDDO
223 et(1:nel)=
max(eti(1:nel),et(1:nel))
224 ENDDO
225 efac=two
226 DO i=1,nel
227 et(i)=
max(one,efac*et(i)/gmax)
228 ENDDO
229 ENDIF
230
231
232 rv(1:nel) = ev(1:nel,1)*ev(1:nel,2)*ev(1:nel,3)
233 ec(1:nel,1:3) = ev(1:nel,1:3)**2
234
235
236 IF(iexpan > 0.AND.(ismstr==10.OR.ismstr==11.OR.ismstr==12)) rv(1:nel)= rv(1:nel)-epsth3(1:nel)
237
238
239 p(1:nel) = zero
240 DO ii = 1,nordre
241 fac = two*mu(ii)/(al(ii))
242 IF(fac/=zero) THEN
243 IF(al(ii)/=zero) THEN
244 DO i=1,nel
245 fac1 = fac / rv(i)
246 IF(ev(i,1)>zero) p(i) = p(i) +
247 . fac1*third*( exp(al(ii)*log(ev(i,1))) )
248 IF(ev(i,2)>zero) p(i) = p(i) +
249 . fac1*third*( exp(al(ii)*log(ev(i,2))) )
250 IF(ev(i,3)>zero) p(i) = p(i) +
251 . fac1*third*( exp(al(ii)*log(ev(i,3))) )
252
253 IF(rv(i)>zero) p(i) = p(i)
254 . - fac1*exp((-beta(ii)*al(ii))*log(rv(i)))
255 ENDDO
256 ELSE
257 DO i=1,nel
258 fac1 = fac / rv(i)
259 IF(ev(i,1)/=zero) p(i) = p(i) + fac1*third
260 IF(ev(i,2)/=zero) p(i) = p(i) + fac1*third
261 IF(ev(i,3)/=zero) p(i) = p(i) + fac1*third
262 IF(rv(i)>zero) p(i) = p(i) - fac1
263 ENDDO
264 ENDIF
265 ELSE
266 p(1:nel) = zero
267 ENDIF
268 ENDDO
269
270 DO j=1,3
271 s(1:nel,j)= zero
272 DO ii= 1, nordre
273 fac = two*mu(ii)/(al(ii))
274 IF(fac/=zero) THEN
275 DO i=1,nel
276 fac1 = fac/(ev(i,j)**2)
277 IF(al(ii)/=zero) THEN
278 IF(ev(i,j)>zero) s(i,j)= s(i,j) + fac1*exp(al(ii)*log(ev(i,j)))
279 IF(rv(i)>zero) s(i,j)= s(i,j) - fac1*exp((-beta(ii)*al(ii))*log(rv(i)))
280 ELSE
281 IF(ev(i,j)==zero) THEN
282 s(i,j)= s(i,j) - fac1
283 ENDIF
284 ENDIF
285 ENDDO
286 ENDIF
287 ENDDO
288 ENDDO
289
290 DO i=1,nel
291
292
293
294
295 IF(rv(i)>zero) THEN
296 rv_pui_tab(i) = exp( two_third*log(rv(i)) )
297 ELSE
298 rv_pui_tab(i) = zero
299 ENDIF
300 ENDDO
301 DO j=1,3
302 ssp(1:nel,j) = p(1:nel)/(ev(1:nel,j)**2)
303 sd(1:nel,j) = s(1:nel,j) - ssp(1:nel,j)*rv(1:nel)
304 sd(1:nel,j) = rv_pui_tab(1:nel)*sd(1:nel,j)
305 ENDDO
306
307 IF(ivisc == 2) sd(1:nel,1:3) = s(1:nel,1:3)
308
309 gvis = zero
310
311 IF(ivisc > 0)THEN
312 gvis = gmax
313
314 sdg0(1:nel,1:6) = uvar(1:nel,1:6)
315 DO i=1,nel
316
317
318
319 sdg(i,1) = dirprv(i,1,1)*dirprv(i,1,1)*sd(i,1)
320 . + dirprv(i,1,2)*dirprv(i,1,2)*sd(i,2)
321 . + dirprv(i,1,3)*dirprv(i,1,3)*sd(i,3)
322
323 sdg(i,2) = dirprv(i,2,2)*dirprv(i,2,2)*sd(i,2)
324 . + dirprv(i,2,3)*dirprv(i,2,3)*sd(i,3)
325 . + dirprv(i,2,1)*dirprv(i,2,1)*sd(i,1)
326
327 sdg(i,3) = dirprv(i,3,3)*dirprv(i,3,3)*sd(i,3)
328 . + dirprv(i,3,1)*dirprv(i,3,1)*sd(i,1)
329 . + dirprv(i,3,2)*dirprv(i,3,2)*sd(i,2)
330
331 sdg(i,4) = dirprv(i,1,1)*dirprv(i,2,1)*sd(i,1)
332 . + dirprv(i,1,2)*dirprv(i,2,2)*sd(i,2)
333 . + dirprv(i,1,3)*dirprv(i,2,3)*sd(i,3)
334
335 sdg(i,5) = dirprv(i,2,2)*dirprv(i,3,2)*sd(i,2)
336 . + dirprv(i,2,3)*dirprv(i,3,3)*sd(i,3)
337 . + dirprv(i,2,1)*dirprv(i,3,1)*sd(i,1)
338
339 sdg(i,6) = dirprv(i,3,3)*dirprv(i,1,3)*sd(i,3)
340 . + dirprv(i,3,1)*dirprv(i,1,1)*sd(i,1)
341 . + dirprv(i,3,2)*dirprv(i,1,2)*sd(i,2)
342
343 ENDDO
344 DO j=1,6
345 uvar(1:nel,j) = sdg(1:nel,j)
346 ENDDO
347
348 DO ii= 1,nprony
349 fac= -timestep/taux(ii)
350 DO j= 1,6
351 DO i=1,nel
352 hg0(j) = uvar(i, 6 + (ii - 1)*6 + j)
353 hg(i,j) = exp(fac)*hg0(j) + exp(half*fac)*(sdg(i,j)-sdg0(i,j))
354 uvar(i, 6 + (ii - 1)*6 + j)= hg(i,j)
355 ENDDO
356 ENDDO
357
358
359
360 DO i=1,nel
361 a11 = dirprv(i,1,1)*hg(i,1) + dirprv(i,2,1)*hg(i,4)
362 . + dirprv(i,3,1)*hg(i,6)
363 a12 = dirprv(i,1,2)*hg(i,1) + dirprv(i,2,2)*hg(i,4)
364 . + dirprv(i,3,2)*hg(i,6)
365 a13 = dirprv(i,1,3)*hg(i,1) + dirprv(i,2,3)*hg(i,4)
366 . + dirprv(i,3,3)*hg(i,6)
367
368
369 a21 = dirprv(i,1,1)*hg(i,4) + dirprv(i,2,1)*hg(i,2)
370 . + dirprv(i,3,1)*hg(i,5)
371 a22 = dirprv(i,1,2)*hg(i,4) + dirprv(i,2,2)*hg(i,2)
372 . + dirprv(i,3,2)*hg(i,5)
373 a23 = dirprv(i,1,3)*hg(i,4) + dirprv(i,2,3)*hg(i,2)
374 . + dirprv(i,3,3)*hg(i,5)
375
376 a31 = dirprv(i,1,1)*hg(i,6) + dirprv(i,2,1)*hg(i,5)
377 . + dirprv(i,3,1)*hg(i,3)
378 a32 = dirprv(i,1,2)*hg(i,6) + dirprv(i,2,2)*hg(i,5)
379 . + dirprv(i,3,2)*hg(i,3)
380 a33 = dirprv(i,1,3)*hg(i,6) + dirprv(i,2,3)*hg(i,5)
381 . + dirprv(i,3,3)*hg(i,3)
382
383 h(i,1,ii) = dirprv(i,1,1)*a11 + dirprv(i,2,1)*a21
384 . + dirprv(i,3,1)*a31
385 h(i,2,ii) = dirprv(i,1,2)*a12 + dirprv(i,2,2)*a22
386 . + dirprv(i,3,2)*a32
387 h(i,3,ii) = dirprv(i,1,3)*a13 + dirprv(i,2,3)*a23
388 . + dirprv(i,3,3)*a33
389 ENDDO
390 DO j=1,3
391 DO i=1,nel
392
393 fac = third*(h(i,1,ii)*ec(i,1) + h(i,2,ii)*ec(i,2)
394 . + h(i,3,ii)*ec(i,3))
395 hd(i,j,ii) = h(i,j,ii) - fac/
max(em20,ec(i,j))
396 ENDDO
397 ENDDO
398 ENDDO
399 ENDIF
400
401 IF(ivisc == 1) THEN
402 DO i =1,nel
403
404
405 ! else pui = 0
406 IF(rv(i)>zero) THEN
407 rv_pui_tab(i) = exp( (-two_third)*log(rv(i)) )
408 ELSE
409 rv_pui_tab(i) = zero
410 ENDIF
411 ENDDO
412
413 DO j=1,3
414 DO i =1,nel
415 fac= gamainf*rv_pui_tab(i)
416 s(i,j) = ssp(i,j)*rv(i) + fac*sd(i,j)
417 ENDDO
418 DO ii=1,nprony
419 DO i =1,nel
420 s(i,j) = s(i,j) + gama(ii)*rv_pui_tab(i)*hd(i,j,ii)
421 ENDDO
422 ENDDO
423
424 DO i =1,nel
425 s(i,j) = (s(i,j)*ev(i,j)**2)/rv(i)
426 ENDDO
427 ENDDO
428 ELSEIF(ivisc == 2) THEN
429 DO j=1,3
430 DO i =1,nel
431 s(i,j) = gamainf*s(i,j)
432 ENDDO
433 DO ii=1,nprony
434 DO i =1,nel
435 s(i,j) = s(i,j) + gama(ii)*h(i,j,ii)
436 ENDDO
437 ENDDO
438
439 DO i =1,nel
440 s(i,j) = (s(i,j)*ev(i,j)**2)/rv(i)
441 ENDDO
442 ENDDO
443 ELSE
444
445 DO j=1,3
446 DO i =1,nel
447 s(i,j) = s(i,j)*ev(i,j)**2/rv(i)
448 ENDDO
449 ENDDO
450 ENDIF
451
452 cmax0 = two_third*gmax+ rbulk
453 ai(1:nel,1:3) = zero
454 bi(1:nel,1:3) = zero
455 cj(1:nel) = zero
456 DO ii = 1,nordre
457 IF(al(ii)/=zero) THEN
458 fac = one/al(ii)
459 fac1 = beta(ii)+fac
460 mu2 = two*mu(ii)
461 DO i=1,nel
462 lam_al(1:3) = exp(al(ii)*log(ev(i,1:3)))
463 rvl = fac1*exp(-beta(ii)*al(ii)*log(rv(i)))
464 cj(i) = cj(i) + mu2*rvl
465 ai(i,1:3) = ai(i,1:3) + mu2*lam_al(1:3)
466 bi(i,1:3) = bi(i,1:3) + fac*mu2*lam_al(1:3)
467 ENDDO
468 ENDIF
469 ENDDO
470
471 DO i=1,nel
472 d11 = ai(i,1)-bi(i,1)+cj(i)
473 d22 = ai(i,2)-bi(i,2)+cj(i)
474 d33 = ai(i,3)-bi(i,3)+cj(i)
475
476 fac = one/rv(i)
477 cmax = fac*
max(d11,d22,d33)
478 cimax(i) = two_third*gvis +
max(cmax,cmax0)
479 et(i)= one
480 ENDDO
481
482 DO i=1,nel
483 signxx(i) = dirprv(i,1,1)*dirprv(i,1,1)*s(i,1)
484 . + dirprv(i,1,2)*dirprv(i,1,2)*s(i,2)
485 . + dirprv(i,1,3)*dirprv(i,1,3)*s(i,3)
486
487 signyy(i) = dirprv(i,2,2)*dirprv(i,2,2)*s(i,2)
488 . + dirprv(i,2,3)*dirprv(i,2,3)*s(i,3)
489 . + dirprv(i,2,1)*dirprv(i,2,1)*s(i,1)
490
491 signzz(i) = dirprv(i,3,3)*dirprv(i,3,3)*s(i,3)
492 . + dirprv(i,3,1)*dirprv(i,3,1)*s(i,1)
493 . + dirprv(i,3,2)*dirprv(i,3,2)*s(i,2)
494
495 signxy(i) = dirprv(i,1,1)*dirprv(i,2,1)*s(i,1)
496 . + dirprv(i,1,2)*dirprv(i,2,2)*s(i,2)
497 . + dirprv(i,1,3)*dirprv(i,2,3)*s(i,3)
498
499 signyz(i) = dirprv(i,2,2)*dirprv(i,3,2)*s(i,2)
500 . + dirprv(i,2,3)*dirprv(i,3,3)*s(i,3)
501 . + dirprv(i,2,1)*dirprv(i,3,1)*s(i,1)
502
503 signzx(i) = dirprv(i,3,3)*dirprv(i,1,3)*s(i,3)
504 . + dirprv(i,3,1)*dirprv(i,1,1)*s(i,1)
505 . + dirprv(i,3,2)*dirprv(i,1,2)*s(i,2)
506
507
508
509
510 soundsp(i)=sqrt(cimax(i)/rho(i))
511
512
513 viscmax(i) = zero
514 ENDDO
515 RETURN
subroutine valpvecdp_v(sig, val, vec, nel)
subroutine valpvec_v(sig, val, vec, nel)