43
44
45
46#include "implicit_f.inc"
47
48
49
50
51
52
53
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
94
95
96
97
98 INTEGER NEL, NUPARAM,
100 . time,timestep,uparam(nuparam),
101 . rho(nel),rho0(nel),volume(nel),eint(nel),
102 . epspxx(nel),epspyy(nel),epspzz(nel),
103 . epspxy(nel),epspyz(nel),epspzx(nel),
104 . depsxx(nel),depsyy(nel),depszz(nel),
105 . depsxy(nel),depsyz(nel),depszx(nel),
106 . epsxx(nel) ,epsyy(nel) ,epszz(nel) ,
107 . epsxy(nel) ,epsyz(nel) ,epszx(nel) ,
108 . sigoxx(nel),sigoyy(nel),sigozz(nel),
109 . sigoxy(nel),sigoyz(nel),sigozx(nel)
110
111
112
114 . signxx(nel),signyy(nel),signzz(nel),
115 . signxy(nel),signyz(nel),signzx(nel),
116 . sigvxx(nel),sigvyy(nel),sigvzz(nel),
117 . sigvxy(nel),sigvyz(nel),sigvzx(nel),
118 . soundsp(nel),viscmax(nel)
119
120
121
123 . uvar(nel,nuvar), off(nel)
124
125
126
127 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
129 . finter
130 EXTERNAL finter
131
132
133
134
135
136
137
138
139
140
141 INTEGER I,J,K,M
143 . ak,g0,g(5),beta(5),
144 . ev,edxx,edyy,edzz,edxy,edyz,edzx,
145 . sxx,syy,szz,sxy,syz,szx,
146 . dt,sigv,sdv(30),
147 . epxx,epyy,epzz,epxy,epyz,epzx,gt,
148
149 . astas,bstas,vmisk,
150 . ssig1,ssig2
152 . eprnxx,eprnyy,eprnzz,eprnxy,eprnyz,eprnzx,
153 . evr,edrnxx,edrnyy,edrnzz,edrnxy,edrnyz,edrnzx,
154 . edrvxx,edrvyy,edrvzz,edrvxy,edrvyz,edrvzx
156 . axx,ayy,azz,axy,ayz,azx,
157 . bxx,byy,bzz,bxy,byz,bzx,
158 . aaxx,aayy,aazz,aaxy,aayz,aazx,
159 . bbxx,bbyy,bbzz,bbxy,bbyz,bbzx,
160 . ccxx,ccyy,cczz,ccxy,ccyz,cczx,
161 . evsdxx,evsdyy,evsdzz,evsdxy,evsdyz,evsdzx
162
163
164
165
166
167 ak = uparam(1)
168 g0 = two*uparam(2)
169 g(1) = two*uparam(3)
170 g(2) = two*uparam(4)
171 g(3) = two*uparam(5)
172 g(4) = two*uparam(6)
173 g(5) = two*uparam(7)
174
175 beta(1) = uparam(8)
176 beta(2) = uparam(9)
177 beta(3) = uparam(10)
178 beta(4) = uparam(11)
179 beta(5) = uparam(12)
180 astas = uparam(13)
181 bstas = uparam(14)
182 vmisk = uparam(15)
183
184
185
186 IF(time<=em20)THEN
187 DO 100 j=1,nuvar
188 DO 100 i=1,nel
189 uvar(i,j)=zero
190 100 CONTINUE
191 ENDIF
192
193 gt = g0 + g(1) + g(2) + g(3) + g(4) + g(5)
194
195 dt = timestep
196
197 DO i=1,nel
198
199 sdv(1) = uvar(i,11)
200 sdv(2) = uvar(i,12)
201 sdv(3) = uvar(i,13)
202 sdv(4) = uvar(i,14)
203 sdv(5) = uvar(i,15)
204 sdv(6) = uvar(i,16)
205 sdv(7) = uvar(i,17)
206 sdv(8) = uvar(i,18)
207 sdv(9) = uvar(i,19)
208 sdv(10) = uvar(i,20)
209 sdv(11) = uvar(i,21)
210 sdv(12) = uvar(i,22)
211 sdv(13) = uvar(i,23)
212 sdv(14) = uvar(i,24)
213 sdv(15) = uvar(i,25)
214 sdv(16) = uvar(i,26)
215 sdv(17) = uvar(i,27)
216 sdv(18) = uvar(i,28)
217 sdv(19) = uvar(i,29)
218 sdv(20) = uvar(i,30)
219 sdv(21) = uvar(i,31)
220 sdv(22) = uvar(i,32)
221 sdv(23) = uvar(i,33)
222 sdv(24) = uvar(i,34)
223 sdv(25) = uvar(i,35)
224 sdv(26) = uvar(i,36)
225 sdv(27) = uvar(i,37)
226 sdv(28) = uvar(i,38)
227 sdv(29) = uvar(i,39)
228 sdv(30) = uvar(i,40)
229
230 epxx = epsxx(i)
231 epyy = epsyy(i)
232 epzz = epszz(i)
233 epxy = epsxy(i)
234 epyz = epsyz(i)
235 epzx = epszx(i)
236
237 ev = (epxx + epyy + epzz) * third
238 edxx = epxx - ev
239 edyy = epyy - ev
240 edzz = epzz - ev
241 edxy = epxy * half
242 edyz = epyz * half
243 edzx = epzx * half
244
245 eprnxx = epspxx(i)
246 eprnyy = epspyy(i)
247 eprnzz = epspzz(i)
248 eprnxy = epspxy(i)
249 eprnyz = epspyz(i)
250 eprnzx = epspzx(i)
251
252 evr = (eprnxx + eprnyy + eprnzz) * third
253 edrnxx = eprnxx - evr
254 edrnyy = eprnyy - evr
255 edrnzz = eprnzz - evr
256 edrnxy = eprnxy * half
257 edrnyz = eprnyz * half
258 edrnzx = eprnzx * half
259
260 edrvxx = uvar(i,5)
261 edrvyy = uvar(i,6)
262 edrvzz = uvar(i,7)
263 edrvxy = uvar(i,8)
264 edrvyz = uvar(i,9)
265 edrvzx = uvar(i,10)
266
267 axx = edrvxx
268
269 bxx = two * (edrnxx - edrvxx) /
max(dt,em20)
270 ayy = edrvyy
271 byy = two * (edrnyy - edrvyy) /
max(dt,em20)
272 azz = edrvzz
273 bzz = two * (edrnzz - edrvzz) /
max(dt,em20)
274 axy = edrvxy
275 bxy = two * (edrnxy - edrvxy) /
max(dt,em20)
276 ayz = edrvyz
277 byz = two * (edrnyz - edrvyz) /
max
278 azx = edrvzx
279 bzx = two * (edrnzx - edrvzx) /
max(dt,em20)
280
281
282
283 j = 1
284 k = (j -1)*6
285 aaxx = g(j)/beta(j) * (axx - bxx/beta(j))
286 bbxx = g(j)/beta(j) * bxx
287 ccxx = sdv(k+1) - aaxx
288 evsdxx = aaxx + bbxx*dt + ccxx*exp(-beta(j)*dt)
289 aayy = g(j)/beta(j) * (ayy - byy/beta(j))
290 bbyy = g(j)/beta(j) * byy
291 ccyy = sdv(k+2) - aayy
292 evsdyy = aayy + bbyy*dt + ccyy*exp(-beta(j)*dt)
293 aazz = g(j)/beta(j) * (azz - bzz/beta(j))
294 bbzz = g(j)/beta(j) * bzz
295 cczz = sdv(k+3) - aazz
296 evsdzz = aazz + bbzz*dt + cczz*exp(-beta(j)*dt)
297 aaxy = g(j)/beta(j) * (axy - bxy/beta(j))
298 bbxy = g(j)/beta(j) * bxy
299 ccxy = sdv(k+4) - aaxy
300 evsdxy = aaxy + bbxy*dt + ccxy*exp(-beta(j)*dt)
301 aayz = g(j)/beta(j) * (ayz - byz/beta(j))
302 bbyz = g(j)/beta(j) * byz
303 ccyz = sdv(k+5) - aayz
304 evsdyz = aayz + bbyz*dt + ccyz*exp(-beta(j)*dt)
305 aazx = g(j)/beta(j) * (azx - bzx/beta(j))
306 bbzx = g(j)/beta(j) * bzx
307 cczx = sdv(k+6) - aazx
308 evsdzx = aazx + bbzx*dt + cczx*exp(-beta(j)*dt)
309 uvar(i,(k + 11)) = evsdxx
310 uvar(i,(k + 12)) = evsdyy
311 uvar(i,(k + 13)) = evsdzz
312 uvar(i,(k + 14)) = evsdxy
313 uvar(i,(k + 15)) = evsdyz
314 uvar(i,(k + 16)) = evsdzx
315
316 j = 2
317 k = (j -1)*6
318 aaxx = g(j)/beta(j) * (axx - bxx/beta(j))
319 bbxx = g(j)/beta(j) * bxx
320 ccxx = sdv(k+1) - aaxx
321 evsdxx = aaxx + bbxx*dt + ccxx*exp(-beta(j)*dt)
322 aayy = g(j)/beta(j) * (ayy - byy/beta(j))
323 bbyy = g(j)/beta(j) * byy
324 ccyy = sdv(k+2) - aayy
325 evsdyy = aayy + bbyy*dt + ccyy*exp(-beta(j)*dt)
326 aazz = g(j)/beta(j) * (azz - bzz/beta(j))
327 bbzz = g(j)/beta(j) * bzz
328 cczz = sdv(k+3) - aazz
329 evsdzz = aazz + bbzz*dt + cczz*exp(-beta(j)*dt)
330 aaxy = g(j)/beta(j) * (axy - bxy/beta(j))
331 bbxy = g(j)/beta(j) * bxy
332 ccxy = sdv(k+4) - aaxy
333 evsdxy = aaxy + bbxy*dt + ccxy*exp(-beta(j)*dt)
334 aayz = g(j)/beta(j) * (ayz - byz/beta(j))
335 bbyz = g(j)/beta(j) * byz
336 ccyz = sdv(k+5) - aayz
337 evsdyz = aayz + bbyz*dt + ccyz*exp(-beta(j)*dt)
338 aazx = g(j)/beta(j) * (azx - bzx/beta(j))
339 bbzx = g(j)/beta(j) * bzx
340 cczx = sdv(k+6) - aazx
341 evsdzx = aazx + bbzx*dt + cczx*exp(-beta(j)*dt)
342 uvar(i,(k + 11)) = evsdxx
343 uvar(i,(k + 12)) = evsdyy
344 uvar(i,(k + 13)) = evsdzz
345 uvar(i,(k + 14)) = evsdxy
346 uvar(i,(k + 15)) = evsdyz
347 uvar(i,(k + 16)) = evsdzx
348
349 j = 3
350 k = (j -1)*6
351 aaxx = g(j)/beta(j) * (axx - bxx/beta(j))
352 bbxx = g(j)/beta(j) * bxx
353 ccxx = sdv(k+1) - aaxx
354 evsdxx = aaxx + bbxx*dt + ccxx*exp(-beta(j)*dt)
355 aayy = g(j)/beta(j) * (ayy - byy/beta(j))
356 bbyy = g(j)/beta(j) * byy
357 ccyy = sdv(k+2) - aayy
358 evsdyy = aayy + bbyy*dt + ccyy*exp(-beta(j)*dt)
359 aazz = g(j)/beta(j) * (azz - bzz/beta(j))
360 bbzz = g(j)/beta(j) * bzz
361 cczz = sdv(k+3) - aazz
362 evsdzz = aazz + bbzz*dt + cczz*exp(-beta(j)*dt)
363 aaxy = g(j)/beta(j) * (axy - bxy/beta(j))
364 bbxy = g(j)/beta(j) * bxy
365 ccxy = sdv(k+4) - aaxy
366 evsdxy = aaxy + bbxy*dt + ccxy*exp(-beta(j)*dt)
367 aayz = g(j)/beta(j) * (ayz - byz/beta(j))
368 bbyz = g(j)/beta(j) * byz
369 ccyz = sdv(k+5) - aayz
370 evsdyz = aayz + bbyz*dt + ccyz*exp(-beta(j)*dt)
371 aazx = g(j)/beta(j) * (azx - bzx/beta(j))
372 bbzx = g(j)/beta(j) * bzx
373 cczx = sdv(k+6) - aazx
374 evsdzx = aazx + bbzx*dt + cczx*exp(-beta(j)*dt)
375 uvar(i,(k + 11)) = evsdxx
376 uvar(i,(k + 12)) = evsdyy
377 uvar(i,(k + 13)) = evsdzz
378 uvar(i,(k + 14)) = evsdxy
379 uvar(i,(k + 15)) = evsdyz
380 uvar(i,(k + 16)) = evsdzx
381
382 j = 4
383 k = (j -1)*6
384 aaxx = g(j)/beta(j) * (axx - bxx/beta(j))
385 bbxx = g(j)/beta(j) * bxx
386 ccxx = sdv(k+1) - aaxx
387 evsdxx = aaxx + bbxx*dt + ccxx*exp(-beta(j)*dt)
388 aayy = g(j)/beta(j) * (ayy - byy/beta(j))
389 bbyy = g(j)/beta(j) * byy
390 ccyy = sdv(k+2) - aayy
391 evsdyy = aayy + bbyy*dt + ccyy*exp(-beta(j)*dt)
392 aazz = g(j)/beta(j) * (azz - bzz/beta(j))
393 bbzz = g(j)/beta(j) * bzz
394 cczz = sdv(k+3) - aazz
395 evsdzz = aazz + bbzz*dt + cczz*exp(-beta(j)*dt)
396 aaxy = g(j)/beta(j) * (axy - bxy/beta(j))
397 bbxy = g(j)/beta(j) * bxy
398 ccxy = sdv(k+4) - aaxy
399 evsdxy = aaxy + bbxy*dt + ccxy*exp(-beta(j)*dt)
400 aayz = g(j)/beta(j) * (ayz - byz/beta(j))
401 bbyz = g(j)/beta(j) * byz
402 ccyz = sdv(k+5) - aayz
403 evsdyz = aayz + bbyz*dt + ccyz*exp(-beta(j)*dt)
404 aazx = g(j)/beta(j) * (azx - bzx/beta(j))
405 bbzx = g(j)/beta(j) * bzx
406 cczx = sdv(k+6) - aazx
407 evsdzx = aazx + bbzx*dt + cczx*exp(-beta(j)*dt)
408 uvar(i,(k + 11)) = evsdxx
409 uvar(i,(k + 12)) = evsdyy
410 uvar(i,(k + 13)) = evsdzz
411 uvar(i,(k + 14)) = evsdxy
412 uvar(i,(k + 15)) = evsdyz
413 uvar(i,(k + 16)) = evsdzx
414
415 j = 5
416 k = (j -1)*6
417 aaxx = g(j)/beta(j) * (axx - bxx/beta(j))
418 bbxx = g(j)/beta(j) * bxx
419 ccxx = sdv(k+1) - aaxx
420 evsdxx = aaxx + bbxx*dt + ccxx*exp(-beta(j)*dt)
421 aayy = g(j)/beta(j) * (ayy - byy/beta(j))
422 bbyy = g(j)/beta(j) * byy
423 ccyy = sdv(k+2) - aayy
424 evsdyy = aayy + bbyy*dt + ccyy*exp(-beta(j)*dt)
425 aazz = g(j)/beta(j) * (azz - bzz/beta(j))
426 bbzz = g(j)/beta(j) * bzz
427 cczz = sdv(k+3) - aazz
428 evsdzz = aazz + bbzz*dt + cczz*exp(-beta(j)*dt)
429 aaxy = g(j)/beta(j) * (axy - bxy/beta(j))
430 bbxy = g(j)/beta(j) * bxy
431 ccxy = sdv(k+4) - aaxy
432 evsdxy = aaxy + bbxy*dt + ccxy*exp(-beta(j)*dt)
433 aayz = g(j)/beta(j) * (ayz - byz/beta(j))
434 bbyz = g(j)/beta(j) * byz
435 ccyz = sdv(k+5) - aayz
436 evsdyz = aayz + bbyz*dt + ccyz*exp(-beta(j)*dt)
437 aazx = g(j)/beta(j) * (azx - bzx/beta(j))
438 bbzx = g(j)/beta(j) * bzx
439 cczx = sdv(k+6) - aazx
440 evsdzx = aazx + bbzx*dt + cczx*exp(-beta(j)*dt)
441 uvar(i,(k + 11)) = evsdxx
442 uvar(i,(k + 12)) = evsdyy
443 uvar(i,(k + 13)) = evsdzz
444 uvar(i,(k + 14)) = evsdxy
445 uvar(i,(k + 15)) = evsdyz
446 uvar(i,(k + 16)) = evsdzx
447
448 sigv = (sigoxx(i) + sigoyy(i) + sigozz(i))* third
449 sigv = sigv + (3. * ak) *
450 . (depsxx(i) + depsyy(i) + depszz(i)) * third
451
452
453 sxx = g0 * edxx
454 syy = g0 * edyy
455 szz = g0 * edzz
456 sxy = g0 * edxy
457 syz = g0 * edyz
458 szx = g0 * edzx
459
460 sxx = sxx + uvar(i,11)
461 syy = syy + uvar(i,12)
462 szz = szz + uvar(i,13)
463 sxy = sxy + uvar(i,14)
464 syz = syz + uvar(i,15)
465 szx = szx + uvar(i,16)
466
467 sxx = sxx + uvar(i,17)
468 syy = syy + uvar(i,18)
469 szz = szz + uvar(i,19)
470 sxy = sxy + uvar(i,20)
471 syz = syz + uvar(i,21)
472 szx = szx + uvar(i,22)
473
474 sxx = sxx + uvar(i,23)
475 syy = syy + uvar(i,24)
476 szz = szz + uvar(i,25)
477 sxy = sxy + uvar(i,26)
478 syz = syz + uvar(i,27)
479 szx = szx + uvar(i,28)
480
481 sxx = sxx + uvar(i,29)
482 syy = syy + uvar(i,30)
483 szz = szz + uvar(i,31)
484 sxy = sxy + uvar(i,32)
485 syz = syz + uvar(i,33)
486 szx = szx + uvar(i,34)
487
488 sxx = sxx + uvar(i,35)
489 syy = syy + uvar(i,36)
490 szz = szz + uvar(i,37)
491 sxy = sxy + uvar(i,38)
492 syz = syz + uvar(i,39)
493 szx = szx + uvar(i,40)
494
495 uvar(i,5) = two * edrnxx - edrvxx
496 uvar(i,6) = two * edrnyy - edrvyy
497 uvar(i,7) = two * edrnzz - edrvzz
498 uvar(i,8) = two * edrnxy - edrvxy
499 uvar(i,9) = two * edrnyz - edrvyz
500 uvar(i,10) = two * edrnzx - edrvzx
501
502 signxx(i) = sxx + sigv
503 signyy(i) = syy + sigv
504 signzz(i) = szz + sigv
505 signxy(i) = sxy
506 signyz(i) = syz
507 signzx(i) = szx
508
509 soundsp(i) = sqrt((ak/rho(i)) +
510 . ((two*two*gt)/(rho(i)*three)))
511 viscmax(i) = zero
512
513
514
515
516
517 ssig1=signxx(i)+signyy(i)+signzz(i)
518 ssig2=half*(sxx**2+syy**2+szz**2)+sxy**2+syz**2+szx**2
519 ssig2=three*ssig2
520
521 IF(ssig2>zero) THEN
522 uvar(i,1)=exp( half*log(ssig2) )/vmisk
523 ELSE
524 uvar(i,1) = zero
525 ENDIF
526
527 IF( (ssig1**2.+astas*2.*ssig2)>zero ) THEN
528 uvar(i,2)=( ssig1+exp( half*log(ssig1**2+astas*2*ssig2) ) )/bstas
529 ELSE
530 uvar(i,2) = ssig1/bstas
531 ENDIF
532 uvar(i,3)=
max(uvar(i,3),uvar(i,1))
533 uvar(i,4)=
max(uvar(i,4),uvar(i,2))
534
535 ENDDO
536 RETURN