47
48
49
50#include "implicit_f.inc"
51
52
53
54#include "mvsiz_p.inc"
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
105
106#include "param_c.inc"
107#include "com01_c.inc"
108
109
110
111
112 INTEGER NEL0, NUPARAM, NUVAR, NPT0, IPT,IFLAG(*),
113 . NGL(NEL0),MAT(NEL0),ISRATE,IPM(NPROPMI,*),INLOC
115 my_real ,
DIMENSION(NEL0) ,
INTENT(IN) :: epsd_pg
116 my_real ,
DIMENSION(NEL0) ,
INTENT(INOUT) :: epsp
117 my_real time,timestep,uparam(*),
118 .
area(nel0),rho0(nel0),eint(nel0,2),
119 . thkly(nel0),pla(nel0),
120 . epspxx(nel0),epspyy(nel0),
121 . epspxy(nel0),epspyz(nel0),epspzx(nel0),
122 . depsxx(nel0),depsyy(nel0),
123 . depsxy(nel0),depsyz(nel0),depszx(nel0),
124 . epsxx(nel0) ,epsyy(nel0) ,
125 . epsxy(nel0) ,epsyz(nel0) ,epszx(nel0) ,
126 . sigoxx(nel0),sigoyy(nel0),
127 . sigoxy(nel0),sigoyz(nel0),sigozx(nel0),
128 . gs(*) ,shf(nel0), dplanl(nel0)
129 my_real,
DIMENSION(NEL0),
INTENT(IN) :: loff
130
131
132
134 . signxx(nel0),signyy(nel0),
135 . signxy(nel0),signyz(nel0),signzx(nel0),
136 . sigvxx(nel0),sigvyy(nel0),
137 . sigvxy(nel0),sigvyz(nel0),sigvzx(nel0),
138 . soundsp(nel0),viscmax(nel0),etse(nel0),
139 . dpla_i(nel0)
140
141
142
143 my_real uvar(nel0,nuvar), off(nel0),thk(nel0),yld(nel0)
144
145
146
147 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
149 EXTERNAL finter
150
151
152
153
154
155
156
157
158
159
160 INTEGER I,J,NRATE(MVSIZ),J1(MVSIZ),J2(MVSIZ),N,NINDX,NMAX,IADBUF,NFUNC,
161 . IAD1(MVSIZ),IPOS1(MVSIZ),ILEN1(MVSIZ),
162 . IAD2(MVSIZ),IPOS2(MVSIZ),ILEN2(MVSIZ),
163 . IAD3(MVSIZ),IPOS3(MVSIZ),ILEN3(MVSIZ),
164 . IAD4(MVSIZ),IPOS4(MVSIZ),ILEN4(MVSIZ),
165 . JJ(MVSIZ),INDEX(MVSIZ),IFUNC(MVSIZ,100),NRATEM,
166 . NRATE1,IFUNC1(100),IADBUFV(MVSIZ),NFUNCV(MVSIZ),
167 . , J3(MVSIZ), J4(MVSIZ), NN, INDEX1,MOPTE,MEINF,MCE,IFUNCE
169 . e(mvsiz),a1(mvsiz),a2(mvsiz),g(mvsiz),g3(mvsiz),
170 . dydx1(mvsiz),dydx2(mvsiz),rate(mvsiz,4),svm(mvsiz),
171 . y1(mvsiz),y2(mvsiz),dr(mvsiz),
172 . yfac(mvsiz,4),nnu1(mvsiz),nu1(mvsiz),
173 . nu2(mvsiz),nu3(mvsiz),nu4(mvsiz),nu5(mvsiz),nu6(mvsiz),
174 . aa(mvsiz),bb(mvsiz),dpla_j(mvsiz),
175 . pp(mvsiz),qq(mvsiz),fail(mvsiz),svmo(mvsiz),h(mvsiz),
176 . epsmax(mvsiz),epsr1(mvsiz),epsr2(mvsiz),fisokin(mvsiz),
177 . sigexx(mvsiz),sigeyy(mvsiz),sigexy(mvsiz),
178 . hk(mvsiz),nu(mvsiz),y3(mvsiz),y4(mvsiz),
179 . dydx3(mvsiz),dydx4(mvsiz),einf(mvsiz),
180 . ce(mvsiz),escale(mvsiz),e0(mvsiz),dydxe(mvsiz)
182 . r,umr,a,b,c,amu,s11,s22,s12,p,p2,fac,dezz,
183 . sigz,s1,s2,s3,dpla,vm2,epst,nnu2,ce1, einf1,
184 . err,f,df,pla_i,q2,yld_i,sigpxx,sigpyy,sigpxy,
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, x(mvsiz)
190 . me, ma1, ma2, mg, mnu,
191 . mepsmax,mepsr1,mepsr2,mfisokin,
192 . y11(mvsiz), y22(mvsiz), y33(mvsiz), y44(mvsiz), y(mvsiz),yp(mvsiz),
193 . x1(mvsiz), x2(mvsiz), x3(mvsiz), x4(mvsiz)
194 INTEGER (MVSIZ+1), IC, MNRATE, OPTE1,OPTE(MVSIZ)
195C
196 DATA nmax/3/
197
198
199
200 nfunc = ipm(10,mat(1))
201 DO j=1,nfunc
202 ifunc1(j)=ipm(10+j,mat(1))
203 ENDDO
204
205 iadbuf = ipm(7,mat(1)) - 1
206 e1 = uparam(iadbuf+2)
207 a11 = uparam(iadbuf+3)
208 a21 = uparam(iadbuf+4)
209 g1 = uparam(iadbuf+5)
210 g31 = 3.*g1
211 nux = uparam(iadbuf+6)
212 nrate1 = nint(uparam(iadbuf+1))
213 epsmax1=uparam(iadbuf+6+2*nrate1+1)
214
215 opte1 = uparam(iadbuf+2*nrate1 + 18)
216 einf1 = uparam
217 ce1 = uparam(iadbuf+2*nrate1 + 20)
218
219 DO i=1,nel0
220 e(i) = e1
221 a1(i) = a11
222 a2(i) = a21
223 g(i) = g1
224 g3(i) = g31
225 ENDDO
226 IF (opte1 == 1)THEN
227 ifunce = uparam(iadbuf+2*nrate1 + 17)
228 DO i=1,nel0
229 IF(pla(i) > zero)THEN
230 escale(i) = finter(kfunc(ifunce),pla(i),npf,tf,dydxe(i))
231 ENDIF
232 ENDDO
233 DO i=1,nel0
234 IF(pla(i) > zero)THEN
235 e(i) = escale(i)* e1
236 a1(i) = e(i)/(one - nux*nux)
237 a2(i) = nux*a1(i)
238 g(i) = half*e(i)/(one+nux)
239 g3(i) = three*g(i)
240 gs(i) = g(i) * shf(i)
241 ENDIF
242 ENDDO
243 ELSEIF ( ce1 /= zero) THEN
244 DO i=1,nel0
245 IF(pla(i) > zero)THEN
246 e(i) = e1-(e1-einf1)*(one-exp(-ce1*pla(i)))
247 a1(i) = e(i)/(one - nux*nux)
248 a2(i) = nux*a1(i)
249 g(i) = half*e(i)/(one+nux)
250 g3(i) = three*g(i)
251 gs(i) = g(i) * shf(i)
252 ENDIF
253 ENDDO
254 ENDIF
255 IF(epsmax1==0.)THEN
256 IF(tf(npf(ifunc1(1)+1)-1)==zero)THEN
257 epsmax1=tf(npf(ifunc1(1)+1)-2)
258 ELSE
259 epsmax1= ep30
260 ENDIF
261 ENDIF
262
263 nnu11 = nux / (one - nux)
264 nnu2 = nnu11*nnu11
265 nu11 = one/(one - nux)
266 nu21 = one/(one + nux)
267 nu31 = one - nnu11
268 nu41 = one + nnu2 + nnu11
269 nu51 = one + nnu2 - two*nnu11
270 nu61 = half - nnu2 + half*nnu11
271
272 epsr11 =uparam(iadbuf+6+2*nrate1+2)
273 epsr21 =uparam(iadbuf+6+2*nrate1+3)
274 fisokin1=uparam(iadbuf+6+2*nrate1+8)
275
276 IF (isigi==0) THEN
277 IF(time==zero)THEN
278 DO i=1,nel0
279 uvar(i,1)=zero
280 uvar(i,2)=zero
281 uvar(i,3)=zero
282 uvar(i,4)=zero
283 ENDDO
284 DO j=1,nrate1
285 DO i=1,nel0
286 uvar(i,j+4)=0
287 ENDDO
288 ENDDO
289 ENDIF
290 ENDIF
291
292
293
294 DO i=1,nel0
295
296
297
298 ENDDO
299
300
301 DO i=1,nel0
302
303 signxx(i)=sigoxx(i) - uvar(i,2) +a1(i)*depsxx(i)+a2(i)*depsyy(i)
304 signyy(i)=sigoyy(i) - uvar(i,3) +a2(i)*depsxx(i)+a1(i)*depsyy(i)
305 signxy(i)=sigoxy(i) - uvar(i,4) +g(i) *depsxy(i)
306 signyz(i)=sigoyz(i)+gs(i) *depsyz(i)
307 signzx(i)=sigozx(i)+gs(i) *depszx(i)
308 sigexx(i) = signxx(i)
309 sigeyy(i) = signyy(i)
310 sigexy(i) = signxy(i)
311
312 soundsp(i) = sqrt(a1(i)/rho0(i))
313 viscmax(i) = zero
314 etse(i) = one
315
316
317
318
319 IF (israte==0) THEN
320 epsp(i) = half*( abs(epspxx(i)+epspyy(i))
321 . + sqrt( (epspxx(i)-epspyy(i))*(epspxx(i)-epspyy(i))
322 . + epspxy(i)*epspxy(i) ) )
323 ELSE
324 epsp(i) = asrate*epsd_pg(i) + (one-asrate)*epsp(i)
325 ENDIF
326
327
328
329
330 epst = half*( epsxx(i)+epsyy(i)
331 . + sqrt( (epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
332 . + epsxy(i)*epsxy(i) ) )
333 fail(i) =
max(zero,
min(one,(epsr21-epst)/(epsr21-epsr11)))
334
335 ENDDO
336
337
338
339 DO i=1,nel0
340 jj(i) = 1
341 ENDDO
342 iadbuf = ipm(7,mat(1)) - 1
343
344 DO j=2,nrate1-1
345 DO i=1,nel0
346 IF(epsp(i)>=uparam(iadbuf+6+j)) jj(i) = j
347 ENDDO
348 ENDDO
349#include "vectorize.inc"
350 DO i=1,nel0
351 IF(jj(i)==1)THEN
352 j1(i) = jj(i)
353 j2(i) = j1(i)+1
354 j3(i) = j2(i)+1
355 j4(i) = j3(i)+1
356 ELSEIF(jj(i)==(nrate1-1))THEN
357 j1(i) = jj(i) - 2
358 j2(i) = j1(i)+1
359 j3(i) = j2(i)+1
360 j4(i) = j3(i)+1
361 ELSE
362 j1(i) = jj(i) - 1
363 j2(i) = j1(i)+1
364 j3(i) = j2(i)+1
365 j4(i) = j3(i)+1
366 ENDIF
367 ENDDO
368#include "vectorize.inc"
369 DO i=1,nel0
370 rate(i,1)=uparam(iadbuf + 6 + j1(i) )
371 rate(i,2)=uparam(iadbuf + 6 + j2(i) )
372 rate(i,3)=uparam(iadbuf + 6 + j3(i) )
373 rate(i,4)=uparam(iadbuf + 6 + j4(i) )
374 yfac(i,1)=uparam(iadbuf + 6 + nrate1 + j1(i) )
375 yfac(i,2)=uparam(iadbuf + 6 + nrate1 + j2(i) )
376 yfac(i,3)=uparam(iadbuf + 6 + nrate1 + j3(i) )
377 yfac(i,4)=uparam(iadbuf + 6 + nrate1 + j4(i) )
378 ENDDO
379
380#include "vectorize.inc"
381 DO i=1,nel0
382 IF(jj(i)==1)THEN
383 j1(i) = jj(i)
384 j2(i) = j1(i)+1
385 j3(i) = j2(i)+1
386 j4(i) = j3(i)+1
387 ELSEIF(jj(i)==(nrate1-1))THEN
388 j1(i) = jj(i) - 2
389 j2(i) = j1(i)+1
390 j3(i) = j2(i)+1
391 j4(i) = j3(i)+1
392 ELSE
393 j1(i) = jj(i) - 1
394 j2(i) = j1(i)+1
395 j3(i) = j2(i)+1
396 j4(i) = j3(i)+1
397 ENDIF
398 ENDDO
399#include "vectorize.inc"
400 DO i=1,nel0
401 ipos1(i) = nint(uvar(i,j1(i)+4))
402 iad1(i) = npf(ifunc1(j1(i))) / 2 + 1
403 ilen1(i) = npf(ifunc1(j1(i))+1) / 2 - iad1(i) - ipos1(i)
404 ipos2(i) = nint(uvar(i,j2(i)+4))
405 iad2(i) = npf(ifunc1(j2(i))) / 2 + 1
406 ilen2(i) = npf(ifunc1(j2(i))+1) / 2 - iad2(i) - ipos2(i)
407 ipos3(i) = nint(uvar(i,j3(i)+4))
408 iad3(i) = npf(ifunc1(j3(i))) / 2 + 1
409 ilen3(i) = npf(ifunc1(j3(i))+1) / 2 - iad3(i) - ipos3(i)
410 ipos4(i) = nint(uvar(i,j4(i)+4))
411 iad4(i) = npf(ifunc1(j4(i))) / 2 + 1
412 ilen4(i) = npf(ifunc1(j4(i))+1) / 2 - iad4(i) - ipos4(i)
413 ENDDO
414
415 CALL vinter(tf,iad1,ipos1,ilen1,nel0,pla,dydx1,y1)
416 CALL vinter(tf,iad2,ipos2,ilen2,nel0,pla,dydx2,y2)
417 CALL vinter(tf,iad3,ipos3,ilen3,nel0,pla,dydx3,y3)
418 CALL vinter(tf,iad4,ipos4,ilen4,nel0,pla,dydx4,y4)
419
420
421 IF (fisokin1==0.) THEN
422
423#include "vectorize.inc"
424 DO i=1,nel0
425 IF(jj(i)==1)THEN
426 j1(i) = jj(i)
427 j2(i) = j1(i)+1
428 j3(i) = j2(i)+1
429 j4(i) = j3(i)+1
430 dydx1(i)= dydx1(i)*yfac(i,1)
431 dydx2(i) = dydx2(i)*yfac(i,2)
432 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
433 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
434 ELSEIF(jj(i)==(nrate1-1))THEN
435 j1(i) = jj(i) - 2
436 j2(i) = j1(i)+1
437 j3(i) = j2(i)+1
438 j4(i) = j3(i)+1
439 dydx3(i) = dydx3(i)*yfac(i,3)
440 dydx4(i) = dydx4(i)*yfac(i,4)
441 fac = (epsp(i) - rate(i,3))/(rate(i,4) - rate(i,3))
442 h(i) = fail(i)*(dydx3(i) + fac*(dydx4(i)-dydx3(i)))
443 ELSE
444 j1(i) = jj(i) - 1
445 j2(i) = j1(i)+1
446 j3(i) = j2(i)+1
447 j4(i) = j3(i)+1
448 dydx2(i) = dydx2(i)*yfac(i,2)
449 dydx3(i) = dydx3(i)*yfac(i,3)
450 fac = (epsp(i) - rate(i,2))/(rate(i,3) - rate(i,2))
451 h(i) = fail(i)*(dydx2(i) + fac*(dydx3(i)-dydx2(i)))
452 ENDIF
453 ENDDO
454 DO i=1,nel0
455 uvar(i,j1(i)+4) = ipos1(i)
456 uvar(i,j2(i)+4) = ipos2(i)
457 uvar(i,j3(i)+4) = ipos3(i)
458 uvar(i,j4(i)+4) = ipos4(i)
459 ENDDO
460 x1(1:nel0) = rate(1:nel0,1)
461 x2(1:nel0) = rate(1:nel0,2)
462 x3(1:nel0) = rate(1:nel0,3)
463 x4(1:nel0) = rate(1:nel0,4)
464 x(1:nel0) = epsp(1:nel0)
465 y11(1:nel0) = y1(1:nel0)*yfac(1:nel0,1)
466 y22(1:nel0) = y2(1:nel0)*yfac(1:nel0,2)
467 y33(1:nel0) = y3(1:nel0)*yfac(1:nel0,3)
468 y44(1:nel0) = y4(1:nel0)*yfac(1:nel0,4)
470 . x,y,yp,jj, nrate1)
471 yld(1:nel0) = fail(1:nel0)*y(1:nel0)
472 yld(1:nel0) =
max(yld(1:nel0),em20)
473
474 ELSEIF (fisokin1==one) THEN
475#include "vectorize.inc"
476 DO i=1,nel0
477 IF(jj(i)==1)THEN
478 j1(i) = jj(i)
479 j2(i) = j1(i)+1
480 j3(i) = j2(i)+1
481 j4(i) = j3(i)+1
482 dydx1(i)= dydx1(i)*yfac(i,1)
483 dydx2(i) = dydx2(i)*yfac(i,2)
484 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
485 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
486 ELSEIF(jj(i)==(nrate1-1))THEN
487 j1(i) = jj(i) - 2
488 j2(i) = j1(i)+1
489 j3(i) = j2(i)+1
490 j4(i) = j3(i)+1
491 dydx3(i) = dydx3(i)*yfac(i,3)
492 dydx4(i) = dydx4(i)*yfac(i,4)
493 fac = (epsp(i) - rate(i,3))/(rate(i,4) - rate(i,3))
494 h(i) = fail(i)*(dydx3(i) + fac*(dydx4(i)-dydx3(i)))
495 ELSE
496 j1(i) = jj(i) - 1
497 j2(i) = j1(i)+1
498 j3(i) = j2(i)+1
499 j4(i) = j3(i)+1
500 dydx2(i) = dydx2(i)*yfac(i,2)
501 dydx3(i) = dydx3(i)*yfac(i,3)
502 fac = (epsp(i) - rate(i,2))/(rate(i,3) - rate(i,2))
503 h(i) = fail(i)*(dydx2(i) + fac*(dydx3(i)-dydx2(i)))
504 ENDIF
505 ENDDO
506 DO i=1,nel0
507 uvar(i,j1(i)+4) = ipos1(i)
508 uvar(i,j2(i)+4) = ipos2(i)
509 uvar(i,j3(i)+4) = ipos3(i)
510 uvar(i,j4(i)+4) = ipos4(i)
511 ENDDO
512
513 y11(1:nel0)=tf(npf(ifunc1(j1(1:nel0)))+1)*yfac(1:nel0,1)
514 y22(1:nel0)=tf(npf(ifunc1(j2(1:nel0)))+1)*yfac(1:nel0,2)
515 y33(1:nel0)=tf(npf(ifunc1(j3(1:nel0)))+1)*yfac(1:nel0,3)
516 y44(1:nel0)=tf(npf(ifunc1(j4(1:nel0)))+1)*yfac(1:nel0,4)
517 x1(1:nel0) = rate(1:nel0,1)
518 x2(1:nel0) = rate(1:nel0,2)
519 x3(1:nel0) = rate(1:nel0,3)
520 x4(1:nel0) = rate(1:nel0,4)
521 x(1:nel0) = epsp(1:nel0)
523 . x,y,yp,jj,nrate1)
524 yld(1:nel0) = fail(1:nel0)*y(1:nel0)
525 ELSE
526#include "vectorize.inc"
527 DO i=1,nel0
528 IF(jj(i)==1)THEN
529 j1(i) = jj(i)
530 j2(i) = j1(i)+1
531 j3(i) = j2(i)+1
532 j4(i) = j3(i)+1
533 dydx1(i)= dydx1(i)*yfac(i,1)
534 dydx2(i) = dydx2(i)*yfac(i,2)
535 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
536 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
537 ELSEIF(jj(i)==(nrate1-1))THEN
538 j1(i) = jj(i) - 2
539 j2(i) = j1(i)+1
540 j3(i) = j2(i)+1
541 j4(i) = j3(i)+1
542 dydx3(i) = dydx3(i)*yfac(i,3)
543 dydx4(i) = dydx4(i)*yfac(i,4)
544 fac = (epsp(i) - rate(i,3))/(rate(i,4) - rate(i,3))
545 h(i) = fail(i)*(dydx3(i) + fac*(dydx4(i)-dydx3(i)))
546 ELSE
547 j1(i) = jj(i) - 1
548 j2(i) = j1(i)+1
549 j3(i) = j2(i)+1
550 j4(i) = j3(i)+1
551 dydx2(i) = dydx2(i)*yfac(i,2)
552 dydx3(i) = dydx3(i)*yfac(i,3)
553 fac = (epsp(i) - rate(i,2))/(rate(i,3) - rate(i,2))
554 h(i) = fail(i)*(dydx2(i) + fac*(dydx3(i)-dydx2(i)))
555 ENDIF
556 ENDDO
557 x1(1:nel0) = rate(1:nel0,1)
558 x2(1:nel0) = rate(1:nel0,2)
559 x3(1:nel0) = rate(1:nel0,3)
560 x4(1:nel0) = rate(1:nel0,4)
561 x(1:nel0) = epsp(1:nel0)
562 y11(1:nel0) = y1(1:nel0)*yfac(1:nel0,1)
563 y22(1:nel0) = y2(1:nel0)*yfac(1:nel0,2)
564 y33(1:nel0) = y3(1:nel0)*yfac(1:nel0,3)
565 y44(1:nel0) = y4(1:nel0)*yfac(1:nel0,4)
567 . x,y,yp,jj,nrate1)
568 yld(1:nel0) = fail(1:nel0)*y(1:nel0)
569 yld(1:nel0) =
max(yld(1:nel0),em20)
570
571 DO i=1,nel0
572 uvar(i,j1(i)+4) = ipos1(i)
573 uvar(i,j2(i)+4) = ipos2(i)
574 uvar(i,j3(i)+4) = ipos3(i)
575 uvar(i,j4(i)+4) = ipos4(i)
576 ENDDO
577
578 y11(1:nel0)=tf(npf(ifunc1(j1(1:nel0)))+1)*yfac(1:nel0,1)
579 y22(1:nel0)=tf(npf(ifunc1(j2(1:nel0)))+1)*yfac(1:nel0,2)
580 y33(1:nel0)=tf(npf(ifunc1(j3(1:nel0)))+1)*yfac(1:nel0,3)
581 y44(1:nel0)=tf(npf(ifunc1(j4(1:nel0)))+1)*yfac(1:nel0,4)
583 . x,y,yp,jj,nrate1)
584 yld(1:nel0) = (one -fisokin1) * yld(1:nel0) + fisokin1 * (fail(1:nel0)*y(1:nel0))
585 ENDIF
586
587
588
589
590 IF(iflag(1)==0)THEN
591
592 DO i=1,nel0
593 svm(i)=sqrt(signxx(i)*signxx(i)
594 . +signyy(i)*signyy(i)
595 . -signxx(i)*signyy(i)
596 . + three*signxy(i)*signxy(i))
597 r =
min(one,yld(i)/
max(em20,svm(i)))
598 signxx(i)=signxx(i)*r
599 signyy(i)=signyy(i)*r
600 signxy(i)=signxy(i)*r
601 umr = one - r
602 dpla_i(i) = off(i)*svm(i)*umr/(g3(i)+h(i))
603 pla(i) = pla(i) + dpla_i(i)
604 s1=half*(signxx(i)+signyy(i))
605 IF (inloc == 0) THEN
606 dezz = dpla_i(i) * half*(signxx(i)+signyy(i)) / yld(i)
607 dezz=-(depsxx(i)+depsyy(i))*nnu11-nu31*dezz
608 thk(i) = thk(i) + dezz*thkly(i)*off(i)
609 ENDIF
610 IF(r<one) etse(i)= h(i)/(h(i)+e(i))
611 ENDDO
612
613 ELSEIF(iflag(1)==1)THEN
614
615
616
617 DO i=1,nel0
618 h(i) =
max(zero,h(i))
619 s1=signxx(i)+signyy(i)
620 s2=signxx(i)-signyy(i)
621 s3=signxy(i)
622 aa(i)=fourth*s1*s1
623 bb(i)=three_over_4*s2*s2+three*s3*s3
624 svm(i)=sqrt(aa(i)+bb(i))
625 IF (inloc == 0) THEN
626 dezz = -(depsxx(i)+depsyy(i))*nnu11
627 thk(i) = thk(i) + dezz*thkly(i)*off(i)
628 ENDIF
629 ENDDO
630
631
632
633 nindx=0
634 DO i=1,nel0
635 IF(svm(i)>yld(i).AND.off(i)==one) THEN
636 nindx=nindx+1
637 index(nindx)=i
638 ENDIF
639 ENDDO
640
641 IF(nindx/=0) THEN
642
643
644
645#include "vectorize.inc"
646 DO j=1,nindx
647 i=index(j)
648 dpla_j(i)=(svm(i)-yld(i))/(g3(i)+h(i))
649 etse(i)= h(i)/(h(i)+e(i))
650 hk(i) = h(i)*(one-fisokin1)
651 ENDDO
652
653 DO n=1,nmax
654#include "vectorize.inc"
655 DO j=1,nindx
656 i=index(j)
657 dpla_i(i)=dpla_j(i)
658 yld_i =yld(i)+hk(i)*dpla_i(i)
659 dr(i) =half*e(i)*dpla_i(i)/yld_i
660 pp(i) =one/(one + dr(i)*nu11)
661 qq(i) =one/(one + three*dr(i)*nu21)
662 p2 =pp(i)*pp(i)
663 q2 =qq(i)*qq(i)
664 f =aa(i)*p2+bb(i)*q2-yld_i*yld_i
665 df =-(aa(i)*nu11*p2*pp(i)+ three*bb(i)*nu21*q2*qq(i))
666 . *(e(i)- two*dr(i)*hk(i))/yld_i
667 . - two*hk(i)*yld_i
668 df = sign(
max(abs(df),em20),df)
669 IF(dpla_i(i)>zero) THEN
670 dpla_j(i)=
max(zero,dpla_i(i)-f/df)
671 ELSE
672 dpla_j(i)=zero
673 ENDIF
674 ENDDO
675 ENDDO
676
677
678
679#include "vectorize.inc"
680 DO j=1,nindx
681 i=index(j)
682 pla(i) = pla(i) + dpla_i(i)
683 s1=(signxx(i)+signyy(i))*pp(i)
684 s2=(signxx(i)-signyy(i))*qq(i)
685 signxx(i)=half*(s1+s2)
686 signyy(i)=half*(s1-s2)
687 signxy(i)=signxy(i)*qq(i)
688 IF (inloc == 0) THEN
689 dezz = - nu31*dr(i)*s1/e(i)
690 thk(i) = thk(i) + dezz*thkly(i)*off(i)
691 ENDIF
692 ENDDO
693 ENDIF
694
695 ELSEIF(iflag(1)==2)THEN
696
697
698
699
700 DO i=1,nel0
701 p = -(signxx(i)+signyy(i))*third
702 s11 = signxx(i)+p
703 s22 = signyy(i)+p
704
705 s12 = signxy(i)
706
707 p2 = p*p
708 vm2= three*(s12*s12 - s11*s22)
709 a = p2*nu41 + vm2
710 vm2= three*p2 + vm2
711 b = p2*nu61
712 c = p2*nu51 - yld(i)*yld(i)
713 r =
min(one,(-b + sqrt(
max(zero,b*b-a*c)))/
max(a ,em20))
714 signxx(i) = s11*r - p
715 signyy(i) = s22*r - p
716 signxy(i) = s12*r
717 umr = one - r
718 sigz = nnu11*p*umr
719 signxx(i) = signxx(i) + sigz
720 signyy(i) = signyy(i) + sigz
721 svm(i)=sqrt(vm2)
722 dpla_i(i) = off(i)*svm(i)*umr/(g3(i)+h(i))
723 pla(i) = pla(i) + dpla_i(i)
724 IF (inloc == 0) THEN
725 dezz = dpla_i(i) * half*(signxx(i)+signyy(i)) / yld(i)
726 dezz=-(depsxx(i)+depsyy(i))*nnu11-nu31*dezz
727 thk(i) = thk(i) + dezz*thkly(i)*off(i)
728 ENDIF
729 IF(r<one) etse(i)= h(i)/(h(i)+e(i))
730 ENDDO
731 ENDIF
732
733 DO i=1,nel0
734 IF(pla(i)>epsmax1.AND.off(i)==one)off(i)=four_over_5
735 ENDDO
736
737
738
739
740 IF (fisokin1/=zero) THEN
741 DO i=1,nel0
742 dsxx = sigexx(i) - signxx(i)
743 dsyy = sigeyy(i) - signyy(i)
744 dsxy = sigexy(i) - signxy(i)
745 dexx = (dsxx - nux*dsyy)
746 deyy = (dsyy - nux*dsxx)
747 dexy = two*(one + nux)*dsxy
748 alpha = fisokin1*h(i)/(e(i)+h(i))/three
749 sigpxx =
alpha*(four*dexx + two*deyy)
750 sigpyy =
alpha*(four*deyy + two*dexx)
752 signxx(i) = signxx(i) + uvar(i,2)
753 signyy(i) = signyy(i) + uvar(i,3)
754 signxy(i) = signxy(i) + uvar(i,4)
755 uvar(i,2) = uvar(i,2) + sigpxx
756 uvar(i,3) = uvar(i,3) + sigpyy
757 uvar(i,4) = uvar(i,4) + sigpxy
758 ENDDO
759 ENDIF
760
761
762
763 IF (inloc > 0) THEN
764 DO i = 1,nel0
765 IF (loff(i) == one) THEN
766 svm(i) = sqrt(signxx(i)*signxx(i)
767 . + signyy(i)*signyy(i)
768 . - signxx(i)*signyy(i)
769 . + three*signxy(i)*signxy(i))
770 dezz =
max(dplanl(i),zero)*half*(signxx(i)+signyy(i))/
max(svm(i),em20)
771 dezz = -nux*((signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/e1) - dezz
772 thk(i) = thk(i) + dezz*thkly(i)*off(i)
773 ENDIF
774 ENDDO
775 ENDIF
776
777 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)
subroutine inter_rat_v(nel0, x0, x1, x2, x3, y0, y1, y2, y3, x, y, yp, jj, n)
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)