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 . (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 NPF(*), MFUNC, KFUNC(MFUNC)
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),NRATEM,
169 . NRATE1,IFUNC1(100),NFUNCV(MVSIZ),
170 . NFUNCM,IPERF
171 INTEGER JST(MVSIZ+1),,MNRATE
173 . mx,mx2,mx3,mx4,me,ma1,ma2,mg,mnu,
174 . mepsmax,mepsr1,mepsr2,mfisokin
176 . dydx1(mvsiz),dydx2(mvsiz),rate(mvsiz,2),svm(mvsiz),
177 . y1(mvsiz),y2(mvsiz),dr(mvsiz),
178 . yfac(mvsiz,2),
179 . aa(mvsiz),bb(mvsiz),dpla_j(mvsiz),
180 . pp(mvsiz),qq(mvsiz),fail(mvsiz),svmo(mvsiz),h(mvsiz),
181 . sigexx(mvsiz),sigeyy(mvsiz),sigexy(mvsiz),
182 . hk(mvsiz)
184 . r,umr,a,b,c,amu,s11,s22,s12,p,p2,fac,dezz,
185 . sigz,s1,s2,s3,dpla,vm2,epst,nnu2,
186 . err,f,df,pla_i,q2,yld_i,sigpxx,sigpyy,sigpxy,
alpha,
187 . e1,a11,a21,g1,g31,nnu11,nu11,nu21,nu31,nu41,nu51,nu61,
188 . epsmax1,epsr11,epsr21,fisokin1,
189 . dsxx,dsyy,dsxy,dexx,deyy,dexy,nux
190
191
193 . ng_e,ng_a1,ng_a2,ng_g,ng_nu,ng_g3,
194 . ng_espmax,ng_epsr1,ng_epsr2,ng_fisokin,
195 . ng_nrates(mvsiz,100),
196 . ng_rate1,ng_rate2,ng_yfac1,ng_yfac2,
197 . ng_ipos1,ng_ipos2,ng_iad1,ng_iad2,ng_ilen1,ng_ilen2,
198 . ng_eps1,ng_eps2,ng_y11,ng_y12,ng_epsmax,
199 . y11(mvsiz),y12(mvsiz)
200 INTEGER NG_NRATE, NG_CRVVAL,EPSMAX0LIST(MVSIZ)
201 integer
202 . k,l,m,jjs(mvsiz+1,100+1),jcnt,
203 . isqrt,isqrtlist(mvsiz),
204 . diff(mvsiz),imix,imixv(mvsiz),tempic,jend(mvsiz),
205 . ifuncm(mvsiz,100),nfuncmv(mvsiz),jlast(mvsiz+1),
206 . iepsmax0,jjindex(mvsiz),jjstart,jjend
207
208 DATA nmax/3/,iperf/0/
209
210
211
212 nfunc = ipm(10,mat(1))
213 DO j=1,nfunc
214 ifunc1(j)=ipm(10+j,mat(1))
215 ENDDO
216
217 iadbuf = ipm(7,mat(1))-1
218 e1 = uparam(iadbuf+2)
219 a11 = uparam(iadbuf+3)
220 a21 = uparam(iadbuf+4)
221 g1 = uparam(iadbuf+5)
222 g31 = three*g1
223 nux = uparam(iadbuf+6)
224 nrate1 = nint(uparam(iadbuf+1))
225 epsmax1=uparam(iadbuf+6+2*nrate1+1)
226 IF(epsmax1==zero)THEN
227 IF(tf(npf(ifunc1(1)+1)-1)==zero)THEN
228 epsmax1=tf(npf(ifunc1(1)+1)-2)
229 ELSE
230 epsmax1= ep30
231 ENDIF
232 ENDIF
233
234 nnu11 = nux / (one - nux)
235 nnu2 = nnu11*nnu11
236 nu11 = one/(one-nux)
237 nu21 = one/(one+nux)
238 nu31 = one-nnu11
239 nu41 = one + nnu2 + nnu11
240 nu51 = one + nnu2 - two*nnu11
241 nu61 = half - nnu2 + half*nnu11
242
243 epsr11 =uparam(iadbuf+6+2*nrate1+2)
244 epsr21 =uparam(iadbuf+6+2*nrate1+3)
245 fisokin1=uparam(iadbuf+6+2*nrate1+8)
246
247 IF (isigi==0) THEN
248 IF(time==zero)THEN
249 DO i=1,nel0
250 uvar(i,1)=zero
251 uvar(i,2)=zero
252 uvar(i,3)=zero
253 uvar(i,4)=zero
254 DO j=1,nrate1
255 uvar(i,j+4)=zero
256 ENDDO
257 ENDDO
258 ENDIF
259 ENDIF
260
261 DO i=1,nel0
262 dpla_i(i) =zero
263 ENDDO
264
265
266 DO i=1,nel0
267
268 signxx(i)=sigoxx(i) - uvar(i,2) +a11*depsxx(i)+a21*depsyy(i)
269 signyy(i)=sigoyy(i) - uvar(i,3) +a21*depsxx(i)+a11*depsyy(i)
270 signxy(i)=sigoxy(i) - uvar(i,4) +g1 *depsxy(i)
271 signyz(i)=sigoyz(i)+gs(i) *depsyz(i)
272 signzx(i)=sigozx(i)+gs(i) *depszx(i)
273 sigexx(i) = signxx(i)
274 sigeyy(i) = signyy(i)
275 sigexy(i) = signxy(i)
276
277 soundsp(i) = sqrt(a11/rho0(i))
278 viscmax(i) = zero
279 etse(i) = one
280
281
282
283 IF (israte==0) THEN
284 epsp(i) = half*( abs(epspxx(i)+epspyy(i))
285 . + sqrt( (epspxx(i)-epspyy(i))*(epspxx(i)-epspyy(i))
286 . + epspxy(i)*epspxy(i) ) )
287 ELSE
288 epsp(i) = asrate*epsd_pg(i) + (one-asrate)*epsp(i)
289 ENDIF
290
291
292
293 epst = half*( epsxx(i)+epsyy(i)
294 . + sqrt( (epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
295 . + epsxy(i)*epsxy(i) ) )
296 fail(i) =
max(zero,
min(one,(epsr21-epst)/(epsr21-epsr11)))
297
298 ENDDO
299
300
301
302 DO i=1,nel0
303 jj(i) = 1
304 ENDDO
305 iadbuf = ipm(7,mat(1))-1
306
307 DO j=2,nrate1-1
308 DO i=1,nel0
309 IF(epsp(i)>=uparam(iadbuf+6+j)) jj(i) = j
310 ENDDO
311 ENDDO
312 DO i=1,nel0
313 rate(i,1)=uparam(iadbuf+6+jj(i))
314 rate(i,2)=uparam(iadbuf+6+jj(i)+1)
315 yfac(i,1)=uparam(iadbuf+6+nrate1+jj(i))
316 yfac(i,2)=uparam(iadbuf+6+nrate1+jj(i)+1)
317 ENDDO
318
319 DO i=1,nel0
320 j1 = jj(i)
321 j2 = j1+1
322 ipos1(i) = nint(uvar(i,j1+4))
323 iad1(i) = npf(ifunc1(j1)) / 2 + 1
324 ilen1(i) = npf(ifunc1(j1)+1) / 2 - iad1(i) - ipos1(i)
325 ipos2(i) = nint(uvar(i,j2+4))
326 iad2(i) = npf(ifunc1(j2)) / 2 + 1
327 ilen2(i) = npf(ifunc1(j2)+1) / 2 - iad2(i) - ipos2(i)
328 ENDDO
329
330 CALL vinter(tf,iad1,ipos1,ilen1,nel0,pla,dydx1,y1)
331 CALL vinter(tf,iad2,ipos2,ilen2,nel0,pla,dydx2,y2)
332
333 IF (fisokin1==zero) THEN
334 DO i=1,nel0
335 j1 = jj(i)
336 j2 = j1+1
337 y1(i)=y1(i)*yfac(i,1)
338 y2(i)=y2(i)*yfac(i,2)
339 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
340 yld(i) = fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
341 yld(i) =
max(yld(i),em20)
342 dydx1(i)=dydx1(i)*yfac(i,1)
343 dydx2(i)=dydx2(i)*yfac(i,2)
344 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
345 uvar(i,j1+4) = ipos1(i)
346 uvar(i,j2+4) = ipos2(i)
347 ENDDO
348 ELSEIF (fisokin1==1.) THEN
349 DO i=1,nel0
350 j1 = jj(i)
351 j2 = j1+1
352 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
353 dydx1(i)=dydx1(i)*yfac(i,1)
354 dydx2(i)=dydx2(i)*yfac(i,2)
355 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
356 uvar(i,j1+4) = ipos1(i)
357 uvar(i,j2+4) = ipos2(i)
358
359 y1(i)=tf(npf(ifunc1(j1))+1)
360 y2(i)=tf(npf(ifunc1(j2))+1)
361 y1(i)=y1(i)*yfac(i,1)
362 y2(i)=y2(i)*yfac(i,2)
363 yld(i) = fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
364 ENDDO
365 ELSE
366 DO i=1,nel0
367 j1 = jj(i)
368 j2 = j1+1
369 y1(i)=y1(i)*yfac(i,1)
370 y2(i)=y2(i)*yfac(i,2)
371 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
372 yld(i) = fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
373 yld(i) =
max(yld(i),em20)
374 dydx1(i)=dydx1(i)*yfac(i,1)
375 dydx2(i)=dydx2(i)*yfac(i,2)
376 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
377 uvar(i,j1+4) = ipos1(i)
378 uvar(i,j2+4) = ipos2(i)
379
380 y1(i)=tf(npf(ifunc1(j1))+1)
381 y2(i)=tf(npf(ifunc1(j2))+1)
382 y1(i)=y1(i)*yfac(i,1)
383 y2(i)=y2(i)*yfac(i,2)
384 yld(i) = (1.-fisokin1) * yld(i) +
385 . fisokin1 * (fail(i)*(y1(i) + fac*(y2(i)-y1(i))))
386 ENDDO
387 ENDIF
388
389
390
391
392 IF(iflag(1)==0)THEN
393
394 DO i=1,nel0
395 svm(i)=sqrt(signxx(i)*signxx(i)
396 . +signyy(i)*signyy(i)
397 .
398 . +three*signxy(i)*signxy(i))
399 r =
min(one,yld(i)/
max(em20,svm(i)))
400 signxx(i)=signxx(i)*r
401 signyy(i)=signyy(i)*r
402 signxy(i)=signxy(i)*r
403 umr = one - r
404 dpla_i(i) = off(i)*svm(i)*umr/(g31+h(i))
405 pla(i) = pla(i) + dpla_i(i)
406 s1=half*(signxx(i)+signyy(i))
407 dezz = dpla_i(i) * half*(signxx(i)+signyy(i)) / yld(i)
408 dezz=-(depsxx(i)+depsyy(i))*nnu11-nu31*dezz
409 thk(i) = thk(i) + dezz*thkly(i)*off(i)
410 IF(r<one) etse(i)= h(i)/(h(i)+e1)
411 ENDDO
412
413 ELSEIF(iflag(1)==1)THEN
414
415
416
417 DO i=1,nel0
418 h(i) =
max(zero,h(i))
419 s1=signxx(i)+signyy(i)
420 s2=signxx(i)-signyy(i)
421 s3=signxy(i)
422 aa(i)=fourth*s1*s1
423 bb(i)=three_over_4*s2*s2+3.*s3*s3
424 svm(i)=sqrt(aa(i)+bb(i))
425 dezz = -(depsxx(i)+depsyy(i))*nnu11
426 thk(i) = thk(i) + dezz*thkly(i)*off(i)
427 ENDDO
428
429
430
431 nindx=0
432 DO i=1,nel0
433 IF(svm(i)>yld(i).AND.off(i)==one) THEN
434 nindx=nindx+1
435 index(nindx)=i
436 ENDIF
437 ENDDO
438
439 IF(nindx/=0) THEN
440
441
442
443 DO j=1,nindx
444 i=index(j)
445 dpla_j(i)=(svm(i)-yld(i))/(g31+h(i))
446 etse(i)= h(i)/(h(i)+e1)
447 hk(i) = h(i)*(one-fisokin1)
448 ENDDO
449
450 DO n=1,nmax
451#include "vectorize.inc"
452 DO j=1,nindx
453 i=index(j)
454 dpla_i(i)=dpla_j(i)
455 yld_i =yld(i)+hk(i)*dpla_i(i)
456 dr(i) =half*e1*dpla_i(i)/yld_i
457 pp(i) =one/(one+dr(i)*nu11)
458 qq(i) =one/(one+three*dr(i)*nu21)
459 p2 =pp(i)*pp(i)
460 q2 =qq(i)*qq(i)
461 f =aa(i)*p2+bb(i)*q2-yld_i*yld_i
462 df =-(aa(i)*nu11*p2*pp(i)+three*bb(i)*nu21*q2*qq(i))
463 . *(e1-two*dr(i)*hk(i))/yld_i
464 . -two*hk(i)*yld_i
465 df = sign(
max(abs(df),em20),df)
466 IF(dpla_i(i)>zero) THEN
467 dpla_j(i)=
max(zero,dpla_i(i)-f/df)
468 ELSE
469 dpla_j(i)=zero
470 ENDIF
471 ENDDO
472 ENDDO
473
474
475
476#include "vectorize.inc"
477 DO j=1,nindx
478 i=index(j)
479 pla(i) = pla(i) + dpla_i(i)
480 s1=(signxx(i)+signyy(i))*pp(i)
481 s2=(signxx(i)-signyy(i))*qq(i)
482 signxx(i)=half*(s1+s2)
483 signyy(i)=half*(s1-s2)
484 signxy(i)=signxy(i)*qq(i)
485 dezz = - nu31*dr(i)*s1/e1
486 thk(i) = thk(i) + dezz*thkly(i)*off(i)
487 ENDDO
488 ENDIF
489
490 ELSEIF(iflag(1)==2)THEN
491
492
493
494
495 DO i=1,nel0
496 p = -(signxx(i)+signyy(i))*third
497 s11 = signxx(i)+p
498 s22 = signyy(i)+p
499
500 s12 = signxy(i)
501
502 p2 = p*p
503 vm2= three*(s12*s12 - s11*s22)
504 a = p2*nu41 + vm2
505 vm2= three*p2 + vm2
506 b = p2*nu61
507 c = p2*nu51 - yld(i)*yld(i)
508 r =
min(one,(-b + sqrt(
max(zero,b*b-a*c)))/
max(a ,em20))
509 signxx(i) = s11*r - p
510 signyy(i) = s22*r - p
511 signxy(i) = s12*r
512
513
514 umr = one - r
515 sigz = nnu11*p*umr
516 signxx(i) = signxx(i) + sigz
517 signyy(i) = signyy(i) + sigz
518 svm(i)=sqrt(vm2)
519 dpla_i(i) = off(i)*svm(i)*umr/(g31+h(i))
520 pla(i) = pla(i) + dpla_i(i)
521 dezz = dpla_i(i) * half*(signxx(i)+signyy(i)) / yld(i)
522 dezz=-(depsxx(i)+depsyy(i))*nnu11-nu31*dezz
523 thk(i) = thk(i) + dezz*thkly(i)*off(i)
524 IF(r<one) etse(i)= h(i)/(h(i)+e1)
525 ENDDO
526 ENDIF
527
528 DO i=1,nel0
529 IF(pla(i)>epsmax1.AND.off(i)==one)off(i)=four_over_5
530 ENDDO
531
532
533
534
535 IF (fisokin1/=zero) THEN
536 DO i=1,nel0
537 dsxx = sigexx(i) - signxx(i)
538 dsyy = sigeyy(i) - signyy(i)
539 dsxy = sigexy(i) - signxy(i)
540 dexx = (dsxx - nux*dsyy)
541 deyy = (dsyy - nux*dsxx)
542 dexy = two*(one+nux)*dsxy
543 alpha = fisokin1*h(i)/(e1+h(i))/three
544 sigpxx =
alpha*(four*dexx+two*deyy)
545 sigpyy =
alpha*(four*deyy+two*dexx)
547 signxx(i) = signxx(i) + uvar(i,2)
548 signyy(i) = signyy(i) + uvar(i,3)
549 signxy(i) = signxy(i) + uvar(i,4)
550 uvar(i,2) = uvar(i,2) + sigpxx
551 uvar(i,3) = uvar(i,3) + sigpyy
552 uvar(i,4) = uvar(i,4) + sigpxy
553 ENDDO
554 ENDIF
555
556 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)