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