45
46
47
48#include "implicit_f.inc"
49
50
51
52#include "mvsiz_p.inc"
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
99
100
101
102
103
104#include "param_c.inc"
105#include "com01_c.inc"
106
107
108
109
110 INTEGER NEL0, NUPARAM, NUVAR, NPT0, IPT,IFLAG(*),
111 . NGL(NEL0),(NEL0),ISRATE,IPM(NPROPMI,*)
113 my_real ,
DIMENSION(NEL0) ,
INTENT(IN) :: epsd_pg
114 my_real ,
DIMENSION(NEL0) ,
INTENT(INOUT) :: epsp
116 . time,timestep,
117 .
area(nel0),rho0(nel0),eint(nel0,2),
118 . thkly(nel0),pla(nel0),
119 . epspxx(nel0),epspyy(nel0),
120 . epspxy(nel0),epspyz(nel0),epspzx(nel0),
121 . depsxx(nel0),depsyy(nel0),
122 . depsxy(nel0),depsyz(nel0),depszx(nel0),
123 . epsxx(nel0) ,epsyy(nel0) ,
124 . epsxy(nel0) ,epsyz(nel0) ,epszx(nel0) ,
125 . sigoxx(nel0),sigoyy(nel0),
126 . sigoxy(nel0),sigoyz(nel0),sigozx(nel0),gs(*)
127 my_real ,
DIMENSION(NUPARAM) :: uparam
128
129
130
132 . signxx(nel0),signyy(nel0),
133 . signxy(nel0),signyz(nel0),signzx(nel0),
134 . sigvxx(nel0),sigvyy(nel0)
135 . sigvxy(nel0),sigvyz(nel0),sigvzx(nel0),
136 . soundsp(nel0),viscmax(nel0),etse(nel0),
137 . dpla_i(nel0)
138
139
140
142 . uvar(nel0,nuvar), off(nel0),thk(nel0),yld(nel0)
143
144
145
146 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
148 . finter ,tf(*)
149 EXTERNAL finter
150
151
152
153
154
155
156
157
158
159
160 INTEGER I,J,NRATE(MVSIZ),J1,J2,N,NINDX,NMAX,
161 . IAD1(MVSIZ),IPOS1(MVSIZ),ILEN1(MVSIZ),
162 . IAD2(MVSIZ),IPOS2(MVSIZ),ILEN2(MVSIZ),
163 . JJ(MVSIZ),INDEX(MVSIZ),NRATEM,
164 . NRATE1
166 . e(mvsiz),a1(mvsiz),a2(mvsiz),g(mvsiz),
167 . dydx1(mvsiz),dydx2(mvsiz),rate(mvsiz),svm(mvsiz),
168 . y1(mvsiz),y2(mvsiz),dr(mvsiz),
169 . yfac(mvsiz,2),
170 . aa(mvsiz),bb(mvsiz),dpla_j(mvsiz),
171 . pp(mvsiz),qq(mvsiz),fail(mvsiz),svmo(mvsiz),h(mvsiz),
172 . epsmax(mvsiz),epsr1(mvsiz),epsr2(mvsiz),fisokin(mvsiz),
173 . sigexx(mvsiz),sigeyy(mvsiz),sigexy(mvsiz),
174 . hk(mvsiz),nu(mvsiz),nnu1(mvsiz),epsf(mvsiz),
175 . svm2(mvsiz),yld2(mvsiz)
177 . r,umr,a,b,c,amu,s11,s22,s12,p,p2,fac,dezz,
178 . sigz,s1,s2,s3,dpla,vm2,epst(mvsiz),nnu2,
179 . err,f,df,pla_i,q2,yld_i,sigpxx,sigpyy,sigpxy,
181 . e1,a11,a21,g1,nnu11,nu11,nu21,nu31,nu41,nu51,nu61,
182 . epsmax1,epsr11,epsr21,fisokin1,g31,
183 . dsxx,dsyy,dsxy,dexx,deyy,dexy,nux, epsf1, mepsf
185 . y11(mvsiz),y21(mvsiz), me, ma1, ma2, mg, mnu,
186 . mepsmax,mepsr1,mepsr2,mfisokin
187 INTEGER JST(MVSIZ+1), IC, MNRATE
188
189 DATA nmax/3/
190
191
192
193 e1 = uparam(2)
194 a11 = uparam(3)
195 a21 = uparam(4)
196 g1 = uparam(5)
197 g31 = three*g1
198 nux = uparam(6)
199 nnu11 = nux / (1. - nux)
200 nrate1 = nint(uparam(1))
201 epsmax1 =uparam(7+2*nrate1)
202 epsr11 =uparam(8+2*nrate1)
203 epsr21 =uparam(9+2*nrate1)
204 fisokin1=uparam(14+2*nrate1)
205 epsf1= uparam(15+2*nrate1)
206 IF(epsmax1==zero)THEN
207 IF(tf(npf(ipm(11,mat(1))+1)-1)==zero)THEN
208 epsmax1=tf(npf(ipm(11,mat(1))+1)-2)
209 ELSE
210 epsmax1= ep30
211 ENDIF
212 ENDIF
213
214 IF (isigi==0) THEN
215 IF(time==zero)THEN
216 DO i=1,nel0
217 uvar(i,1)=zero
218 uvar(i,2)=zero
219 uvar(i,3)=zero
220 uvar(i,4)=zero
221 DO j=1,nrate1
222 uvar(i,j+4)=zero
223 ENDDO
224 ENDDO
225 ENDIF
226 ENDIF
227
228
229
230 DO i=1,nel0
231
232
233
234
235 signxx(i)=sigoxx(i) - uvar(i,2) +a11*depsxx(i)+a21*depsyy(i)
236 signyy(i)=sigoyy(i) - uvar(i,3) +a21*depsxx(i)+a11*depsyy(i)
237 signxy(i)=sigoxy(i) - uvar(i,4) +g1 *depsxy(i)
238 signyz(i)=sigoyz(i)+gs(i) *depsyz(i)
239 signzx(i)=sigozx(i)+gs(i) *depszx(i)
240 sigexx(i) = signxx(i)
241 sigeyy(i) = signyy(i)
242 sigexy(i) = signxy(i)
243
244 soundsp(i) = sqrt(a11/rho0(i))
245 viscmax(i) = zero
246 etse(i) = one
247
248
249
250
251 IF (israte == 0) THEN
252 epsp(i) = half*( abs(epspxx(i)+epspyy(i))
253 . + sqrt( (epspxx(i)-epspyy(i))*(epspxx(i)-epspyy(i))
254 . + epspxy(i)*epspxy(i) ) )
255 ELSE
256 epsp(i) = asrate*epsd_pg(i) + (one-asrate)*epsp(i)
257 END IF
258
259
260
261
262 epst(i) = half*( epsxx(i)+epsyy(i)
263 . + sqrt( (epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
264 . + epsxy(i)*epsxy(i) ) )
265 jj(i) = 1
266 ENDDO
267
268
269
270
271 DO j=2,nrate1-1
272 DO i=1,nel0
273 IF(epsp(i)>=uparam(6+j)) jj(i) = j
274 ENDDO
275 ENDDO
276
277 DO i=1,nel0
278 fac=uparam(6+jj(i))
279 rate(i)=(epsp(i) - fac)/(uparam(7+jj(i)) - fac)
280 yfac(i,1)=uparam(6+nrate1+jj(i))
281 yfac(i,2)=uparam(6+nrate1+jj(i)+1)
282 ENDDO
283
284 DO i=1,nel0
285 j1 = jj(i)
286 j2 = j1+1
287 ipos1(i) = nint(uvar(i,j1+4))
288 iad1(i) = npf(ipm(10+j1,mat(1))) / 2 + 1
289 ilen1(i) = npf(ipm(10+j1,mat(1))+1) / 2 - iad1(i)-ipos1(i)
290 ipos2(i) = nint(uvar(i,j2+4))
291 iad2(i) = npf(ipm(10+j2,mat(1))) / 2 + 1
292 ilen2(i) = npf(ipm(10+j2,mat(1))+1) / 2 - iad2(i)-ipos2(i)
293 END DO
294
295 CALL vinter(tf,iad1,ipos1,ilen1,nel0,pla,dydx1,y1)
296 CALL vinter(tf,iad2,ipos2,ilen2,nel0,pla,dydx2,y2)
297
298 IF (fisokin1==zero) THEN
299 DO i=1,nel0
300 j1 = jj(i)
301 j2 = j1+1
302 y1(i)=y1(i)*yfac(i,1)
303 y2(i)=y2(i)*yfac(i,2)
304 fac = rate(i)
305 yld(i) = y1(i) + fac*(y2(i)-y1(i))
306 yld(i) =
max(yld(i),em20)
307 dydx1(i)=dydx1(i)*yfac(i,1)
308 dydx2(i)=dydx2(i)*yfac(i,2)
309 h(i) =dydx1(i) + fac*(dydx2(i)-dydx1(i))
310 uvar(i,j1+4) = ipos1(i)
311 uvar(i,j2+4) = ipos2(i)
312 ENDDO
313 ELSEIF (fisokin1==1.) THEN
314
315
316
317 DO i=1,nel0
318 j1 = jj(i)
319 j2 = j1+1
320 fac = rate(i)
321 dydx1(i)=dydx1(i)*yfac(i,1)
322 dydx2(i)=dydx2(i)*yfac(i,2)
323 h(i) = dydx1(i) + fac*(dydx2(i)-dydx1(i))
324 uvar(i,j1+4) = ipos1(i)
325 uvar(i,j2+4) = ipos2(i)
326
327 y1(i)=tf(npf(ipm(10+j1,mat(1)))+1)
328 y2(i)=tf(npf(ipm(10+j2,mat(1)))+1)
329 y1(i)=y1(i)*yfac(i,1)
330 y2(i)=y2(i)*yfac(i,2)
331 yld(i) = y1(i) + fac*(y2(i)-y1(i))
332 ENDDO
333
334 ELSE
335
336
337
338 DO i=1,nel0
339 j1 = jj(i)
340 j2 = j1+1
341 y1(i)=y1(i)*yfac(i,1)
342 y2(i)=y2(i)*yfac(i,2)
343 fac = rate(i)
344 yld(i) = y1(i) + fac*(y2(i)-y1(i))
345 yld(i) =
max(yld(i),em20)
346 dydx1(i)=dydx1(i)*yfac(i,1)
347 dydx2(i)=dydx2(i)*yfac(i,2)
348 h(i) = dydx1(i) + fac*(dydx2(i)-dydx1(i))
349 uvar(i,j1+4) = ipos1(i)
350 uvar(i,j2+4) = ipos2(i)
351
352 y1(i)=tf(npf(ipm(10+j1,mat(1)))+1)
353 y2(i)=tf(npf(ipm(10+j2,mat(1)))+1)
354 y1(i)=y1(i)*yfac(i,1)
355 y2(i)=y2(i)*yfac(i,2)
356 yld(i) = (1.-fisokin1) * yld(i) +
357 . fisokin1 * (y1(i) + fac*(y2(i)-y1(i)))
358 ENDDO
359 ENDIF
360
361
362
363
364 IF(iflag(1)==0)THEN
365 nu31 = one-nnu11
366
367 DO i=1,nel0
368 svm(i)=sqrt(signxx(i)*signxx(i)
369 . +signyy(i)*signyy(i)
370 . -signxx(i)*signyy(i)
371 . +three*signxy(i)*signxy(i))
372 r =
min(one,yld(i)/
max(em20,svm(i)))
373 signxx(i)=signxx(i)*r
374 signyy(i)=signyy(i)*r
375 signxy(i)=signxy(i)*r
376 umr = one - r
377 dpla_i(i) = off(i)*svm(i)*umr/(g31+h(i))
378 pla(i) = pla(i) + dpla_i(i)
379 s1=half*(signxx(i)+signyy(i))
380 dezz = dpla_i(i) * half*(signxx(i)+signyy(i)) / yld(i)
381 dezz=-(depsxx(i)+depsyy(i))*nnu11-nu31*dezz
382 thk(i) = thk(i) + dezz*thkly(i)*off(i)
383 IF(r<one) etse(i)= h(i)/(h(i)+e1)
384 ENDDO
385
386 ELSEIF(iflag(1)==1)THEN
387
388
389
390 DO i=1,nel0
391 h(i) =
max(zero,h(i))
392 s1=signxx(i)+signyy(i)
393 s2=signxx(i)-signyy(i)
394 s3=signxy(i)
395 aa(i)=fourth*s1*s1
396 bb(i)=three_over_4*s2*s2+3.*s3*s3
397 svm(i)=sqrt(aa(i)+bb(i))
398 dezz = -(depsxx(i)+depsyy(i))*nnu11
399 thk(i) = thk(i) + dezz*thkly(i)*off(i)
400 ENDDO
401
402
403
404 nindx=0
405 DO i=1,nel0
406 IF(svm(i)>yld(i).AND.off(i)==one) THEN
407 nindx=nindx+1
408 index(nindx)=i
409 ENDIF
410 ENDDO
411
412 IF(nindx/=0) THEN
413
414
415
416 DO j=1,nindx
417 i=index(j)
418 dpla_j(i)=(svm(i)-yld(i))/(g31+h(i))
419 etse(i)= h(i)/(h(i)+e1)
420 hk(i) = h(i)*(one-fisokin1)
421 ENDDO
422
423 nu11 = one/(one-nux)
424 nu21 = one/(one+nux)
425 nu31 = one-nnu11
426 DO n=1,nmax
427#include "vectorize.inc"
428 DO j=1,nindx
429 i=index(j)
430 dpla_i(i)=dpla_j(i)
431 yld_i =yld(i)+hk(i)*dpla_i(i)
432 dr(i) =half*e1*dpla_i(i)/yld_i
433 pp(i) =one/(one+dr(i)*nu11)
434 qq(i) =one/(one+three*dr(i)*nu21)
435 p2 =pp(i)*pp(i)
436 q2 =qq(i)*qq(i)
437 f =aa(i)*p2+bb(i)*q2-yld_i*yld_i
438 df =-(aa(i)*nu11*p2*pp(i)+three*bb(i)*nu21*q2*qq(i))
439 . *(e1-two*dr(i)*hk(i))/yld_i
440 . -two*hk(i)*yld_i
441 df = sign(
max(abs(df),em20),df)
442 IF(dpla_i(i)>zero) THEN
443
444
445
446
447 dpla_j(i)=
max(zero,dpla_i(i)-f/df)
448 ELSE
449 dpla_j(i)=zero
450 ENDIF
451 ENDDO
452 ENDDO
453
454
455
456#include "vectorize.inc"
457 DO j=1,nindx
458 i=index(j)
459 pla(i) = pla(i) + dpla_i(i)
460 s1=(signxx(i)+signyy(i))*pp(i)
461 s2=(signxx(i)-signyy(i))*qq(i)
462 signxx(i)=half*(s1+s2)
463 signyy(i)=half*(s1-s2)
464 signxy(i)=signxy(i)*qq(i)
465 dezz = - nu31*dr(i)*s1/e1
466 thk(i) = thk(i) + dezz*thkly(i)*off(i)
467 ENDDO
468 ENDIF
469
470 ELSEIF(iflag(1)==2)THEN
471
472
473
474 DO i=1,nel0
475 h(i) =
max(zero,h(i))
476 svm2(i)=signxx(i)*signxx(i)
477 . +signyy(i)*signyy(i)
478 . -signxx(i)*signyy(i)
479 . +three*signxy(i)*signxy(i)
480 svm(i)=sqrt(svm2(i))
481 dezz = -(depsxx(i)+depsyy(i))*nnu11
482 thk(i) = thk(i) + dezz*thkly(i)*off(i)
483 ENDDO
484
485
486
487 nindx=0
488 DO i=1,nel0
489 yld2(i)=yld(i)*yld(i)
490 IF(svm2(i)>yld2(i).AND.off(i)==one) THEN
491 nindx=nindx+1
492 index(nindx)=i
493 ENDIF
494 ENDDO
495
496 IF(nindx/=0) THEN
497
498
499
500 nu31 = one-nnu11
501 DO j=1,nindx
502 i=index(j)
503 a=(svm2(i)-yld2(i))
504 . /(five*svm2(i)+three*(-signxx(i)*signyy(i)+signxy(i)*signxy(i)))
505 s1=(one-two*a)*signxx(i)+ a*signyy(i)
506 s2= a*signxx(i)+(one-two*a)*signyy(i)
507 s3=(one-three*a)*signxy(i)
508 signxx(i)=s1
509 signyy(i)=s2
510 signxy(i)=s3
511 dpla_i(i) = off(i)*(svm(i)-yld(i))/(g31+h(i))
512
513 hk(i) = h(i)*(one-fisokin1)
514 yld(i) =yld(i)+hk(i)*dpla_i(i)
515 END DO
516
517 DO j=1,nindx
518 i=index(j)
519 svm(i)=sqrt( signxx(i)*signxx(i)
520 . +signyy(i)*signyy(i)
521 . -signxx(i)*signyy(i)
522 . +three*signxy(i)*signxy(i))
523 r =
min(one,yld(i)/
max(em20,svm(i)))
524 signxx(i)=signxx(i)*r
525 signyy(i)=signyy(i)*r
526 signxy(i)=signxy(i)*r
527 pla(i) = pla(i) + dpla_i(i)
528 dezz = dpla_i(i) * half*(signxx(i)+signyy(i)) / yld(i)
529 dezz = -nu31*dezz
530 thk(i) = thk(i) + dezz*thkly(i)*off(i)
531 etse(i)= h(i)/(h(i)+e1)
532 END DO
533 END IF
534 ENDIF
535
536 DO i=1,nel0
537 IF(off(i)==one)THEN
538 IF(pla(i)>epsmax1.OR.epst(i)>epsf1)off(i)=four_over_5
539 END IF
540 ENDDO
541
542
543
544
545 IF (fisokin1/=zero) THEN
546 DO i=1,nel0
547 dsxx = sigexx(i) - signxx(i)
548 dsyy = sigeyy(i) - signyy(i)
549 dsxy = sigexy(i) - signxy(i)
550 dexx = (dsxx - nux*dsyy)
551 deyy = (dsyy - nux*dsxx)
552 dexy = two*(one+nux)*dsxy
553 alpha = fisokin1*h(i)/(e1+h(i))/three
554 sigpxx =
alpha*(four*dexx+two*deyy)
555 sigpyy =
alpha*(four*deyy+two*dexx)
557
558 signxx(i) = signxx(i) + uvar(i,2)
559 signyy(i) = signyy(i) + uvar(i,3)
560 signxy(i) = signxy(i) + uvar(i,4)
561 uvar(i,2) = uvar(i,2) + sigpxx
562 uvar(i,3) = uvar(i,3) + sigpyy
563 uvar(i,4) = uvar(i,4) + sigpxy
564 ENDDO
565 ENDIF
566
567
568 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)