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,ISRATE, IPT,(*),
111 . NGL(NEL0),MAT(NEL0),IPM(NPROPMI,*)
113 my_real ,
DIMENSION(NEL0) ,
INTENT(IN) :: epsd_pg
114 my_real ,
DIMENSION(NEL0) ,
INTENT(INOUT) :: epsp
115 my_real time,timestep(nel0),uparam(*),
116 .
area(nel0),rho0(nel0),eint(nel0,2),
117 . thkly(nel0),pla(nel0),
118 . epspxx(nel0),epspyy(nel0),
119 . epspxy(nel0),epspyz(nel0),epspzx(nel0),
120 . depsxx(nel0),depsyy(nel0),
121 . depsxy(nel0),depsyz(nel0),depszx(nel0),
122 . epsxx(nel0) ,epsyy(nel0) ,
123 . epsxy(nel0) ,epsyz(nel0) ,epszx(nel0) ,
124 . sigoxx(nel0),sigoyy(nel0),
125 . sigoxy(nel0),sigoyz(nel0),sigozx(nel0),gs(*)
126
127
128
130 . signxx(nel0),signyy(nel0),
131 . signxy(nel0),signyz(nel0),signzx(nel0),
132 . sigvxx(nel0),sigvyy(nel0),
133 . sigvxy(nel0),sigvyz(nel0),sigvzx(nel0),
134 . soundsp(nel0),viscmax(nel0),etse(nel0)
135
136
137
138 my_real uvar(nel0,nuvar), off(nel0),thk(nel0),yld(nel0)
139
140
141
142 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
144 EXTERNAL finter
145
146
147
148
149
150
151
152
153
154
155 INTEGER I,J,NRATE(MVSIZ),J1,J2,N,NINDX,NMAX,IADBUF,NFUNC,
156 . IAD1(MVSIZ),IPOS1(MVSIZ),(MVSIZ),
157 . (MVSIZ),IPOS2(MVSIZ),ILEN2(MVSIZ),
158 . JJ(MVSIZ),INDEX(MVSIZ),IFUNC(MVSIZ,100),NRATEM,
159 . NRATE1,IFUNC1(100),IADBUFV(MVSIZ),NFUNCV(MVSIZ),
160 . NFUNCM, NRATEP(MVSIZ), NRATEN(MVSIZ),I1, I2,
161 . IAD1P(MVSIZ),(MVSIZ),ILEN1P(MVSIZ),
162 . IAD2P(MVSIZ),IPOS2P(MVSIZ),ILEN2P(MVSIZ),
163 . IAD1N(MVSIZ),IPOS1N(MVSIZ),ILEN1N(MVSIZ),
164 . IAD2N(MVSIZ),IPOS2N(MVSIZ),ILEN2N(MVSIZ),II(MVSIZ),
165 . NRATEN1, NRATEP1
167 . e(mvsiz),a1(mvsiz),a2(mvsiz),g(mvsiz),g3(mvsiz),
168 . dydx1(mvsiz),dydx2(mvsiz),rate(mvsiz,2),svm(mvsiz),
169 . y1(mvsiz),y2(mvsiz),dr(mvsiz),epsmax(mvsiz),
170 . yfac(mvsiz,2),nnu1(mvsiz),nu1(mvsiz),
171 . nu2(mvsiz),nu3(mvsiz),nu4(mvsiz),nu5(mvsiz),nu6(mvsiz),
172 . aa(mvsiz),bb(mvsiz),dpla_i(mvsiz),dpla_j(mvsiz),
173 . pp(mvsiz),qq(mvsiz),fail(mvsiz),svmo(mvsiz),h(mvsiz),
174 . epsr1(mvsiz),epsr2(mvsiz),fisokin(mvsiz),
175 . sigexx(mvsiz),sigeyy(mvsiz),sigexy(mvsiz),sigezx(mvsiz),
176 . sigeyz(mvsiz),yn1(mvsiz),yn2(mvsiz),yp1(mvsiz),yp2(mvsiz),
177 . hk(mvsiz),nu(mvsiz), visc(mvsiz
178 . c2(mvsiz),ratep(mvsiz,2),raten(mvsiz,2), yldmax(mvsiz),
179 . yldmin(mvsiz), epst(mvsiz), yldelas(mvsiz),nnu(mvsiz)
181 . r,umr,a,b,c,amu,s11,s22,s12,p,p2,fac
182 . sigz,s1,s2,s3,dpla,vm2,nnu2,
183 . err,f,df,pla_i,q2,yld_i,sigpxx,sigpyy,sigpxy,
185 . e1,a11,a21,g1,g31,nnu11,nu11,nu21,nu31,nu41,nu51,nu61,
186 . visc1, viscv1
187 . dsxx,dsyy,dsxy,dexx,deyy,dexy,nux,c11, c21
189 . y11(mvsiz),y21(mvsiz), me, ma1, ma2, mg, mnu,svm_old,
190 . dev, dtinv, lamda, deps33, svm1
191 INTEGER JST(MVSIZ+1), IC, MNRATE, ITER, MX
192
193
194
195
196
197
198 iflag(1) = 2
199 mx = mat(1)
200 nfunc = ipm(10,mx)
201 DO i= 1, nel0
202
203 DO j=1,nfunc
204 ifunc(i,j)=ipm(10+j,mx)
205 ENDDO
206 ENDDO
207
208 mx = mat(1)
209 iadbuf = ipm(7,mx)-1
210 DO i=1,nel0
211 e(i) = uparam(iadbuf+2)
212 a1(i) = uparam(iadbuf+3)
213 a2(i) = uparam(iadbuf+4)
214 g(i) = uparam(iadbuf+5)
215 g3(i) = three*g(i)
216 nu(i) = uparam(iadbuf+6)
217 visc(i) = uparam(iadbuf+7)
218 viscv(i) = uparam(iadbuf+8)
219 nratep(i) = uparam(iadbuf+9)
220 nraten(i) = uparam(iadbuf+10)
221 epsmax(i) = uparam(iadbuf+11)
222 nrate(i) =nraten(i) + nratep(i) +1
223
224
225 nnu(i) = nu(i) / (one - nu(i))
226 ENDDO
227 IF (isigi==0) THEN
228 IF(time==zero)THEN
229 DO i=1,nel0
230 DO j=1,nrate(i)
231 uvar(i,j)=0
232 ENDDO
233 ENDDO
234 ENDIF
235 ENDIF
236
237
238
239 DO i=1,nel0
240 signxx(i)=sigoxx(i) + a1(i)*depsxx(i)+a2(i)*depsyy(i)
241 signyy(i)=sigoyy(i) + a2(i)*depsxx(i)+a1(i)*depsyy(i)
242 signxy(i)=sigoxy(i) + g(i) *depsxy(i)
243 signyz(i)=sigoyz(i) + gs(i) *depsyz(i)
244 signzx(i)=sigozx(i) + gs(i) *depszx(i)
245
246 visc(i) = visc(i)*rho0(i)
247 viscv(i) = viscv(i)*rho0(i)
248
249 c11 = one/
max(em20,three*visc(i) + four*viscv(i))
250 c1(i) = four* visc(i)*(three*viscv(i)
251 . + visc(i))*c11
252 c2(i) = two*viscv(i)*(three*visc(i)
253 . - two*viscv(i))*c11
254
255 dtinv = timestep(i)**2/
max(em20, timestep(i))
256
257 sigvxx(i)= c2(i)*(epspxx(i)+epspyy(i))
258 . + two*viscv(i)*epspxx(i)
259 sigvyy(i)= c2(i)*(epspxx(i)+epspyy(i))
260 . + two*viscv(i)*epspyy(i)
261 sigvxy(i)= viscv(i)*epspxy(i)
262 sigvyz(i)= viscv(i)*epspyz(i)
263 sigvzx(i)= viscv(i)*epspzx(i)
264
265 soundsp(i) = sqrt(a1(i)/rho0(i))
266 c11 = c2(i) + 2*viscv(i)
267
268 viscmax(i) =
max(c2(i), c11)
269 viscmax(i) = viscmax(i)/
270 . (onep414*rho0(i)*soundsp(i)*sqrt(
area(i)))
271 etse(i) = one
272
273
274
275 IF (israte == 0) THEN
276 epsp(i) = half*(epspxx(i)+epspyy(i)
277 . + sqrt( (epspxx(i)-epspyy(i))*(epspxx(i)-epspyy(i))
278 . + epspxy(i)*epspxy(i) ) )
279 ELSE
280 epsp(i) = asrate*epsd_pg(i) + (one-asrate)*epsp(i)
281 END IF
282
283
284
285
286 epst(i) = half*(epsxx(i)+epsyy(i)
287 . + sqrt( (epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
288 . + epsxy(i)*epsxy(i) ) )
289 IF(off(i)==one.AND.epst(i)>epsmax(i))off(i) =four_over_5
290 ENDDO
291
292
293
294 DO i=1,nel0
295 jj(i) = 1
296 ii(i) = 1 + nratep(i)
297 ENDDO
298
299 DO i=1,nel0
300 DO j=2,nratep(i)-1
301 iadbuf = ipm(7,mat(i)) - 1
302 IF(abs(epsp(i))>=uparam(iadbuf+ 10 +j)) jj(i) = j
303 ENDDO
304 ENDDO
305
306 DO i=1,nel0
307 DO j=2,nraten(i)-1
308 iadbuf = ipm(7,mat(i))-1
309 IF(abs(epsp(i))>=
310 . abs(uparam(iadbuf+10 + nratep(i) + j))) ii(i) = nratep(i) + j
311 ENDDO
312 ENDDO
313 DO i=1,nel0
314 iadbuf = ipm(7,mat(i))-1
315 ratep(i,1)=uparam(iadbuf+10+jj(i))
316 ratep(i,2)=uparam(iadbuf+10+jj(i)+1)
317 raten(i,1)=uparam(iadbuf+10+ii(i))
318 raten(i,2)=uparam(iadbuf+10+ii(i)+1)
319 ENDDO
320
321 DO i=1,nel0
322
323 j1 = jj(i) + 1
324 j2 = j1+1
325 i1 = ii(i) +1
326 i2 = i1 +1
327 ipos1p(i) = nint(uvar(i,j1))
328 iad1p(i) = npf(ifunc(i,j1)) / 2 + 1
329 ilen1p(i) = npf(ifunc(i,j1)+1) / 2 - iad1p(i) - ipos1p(i)
330 ipos2p(i) = nint(uvar(i,j2))
331 iad2p(i) = npf(ifunc(i,j2)) / 2 + 1
332 ilen2p(i) = npf(ifunc(i,j2)+1) / 2 - iad2p(i) - ipos2p(i)
333
334 ipos1n(i) = nint(uvar(i,i1))
335 iad1n(i) = npf(ifunc(i,i1)) / 2 + 1
336 ilen1n(i) = npf(ifunc(i,i1)+1) / 2 - iad1n(i) - ipos1n(i)
337 ipos2n(i) = nint(uvar(i,i2))
338 iad2n(i) = npf(ifunc(i,i2)) / 2 + 1
339 ilen2n(i) = npf(ifunc(i,i2)+1) / 2 - iad2n(i) - ipos2n(i)
340
341 ipos1(i) = nint(uvar(i,1))
342 iad1(i) = npf(ifunc(i,1)) / 2 + 1
343 ilen1(i) = npf(ifunc(i,1)+1) / 2 - iad1(i) - ipos1(i)
344 ENDDO
345
346 CALL vinter(tf,iad1p,ipos1p,ilen1p,nel0,epst,dydx1,yp1)
347 CALL vinter(tf,iad2p,ipos2p,ilen2p,nel0,epst,dydx2,yp2)
348
349 CALL vinter(tf,iad1n,ipos1n,ilen1n,nel0,epst,dydx1,yn1)
350 CALL vinter(tf,iad2n,ipos2n,ilen2n,nel0,epst,dydx2,yn2)
351
352 CALL vinter(tf,iad1,ipos1,ilen1,nel0,epst,dydx1,y1)
353
354
355
356
357 DO i=1,nel0
358 j1 = jj(i) +1
359 j2 = j1+1
360 fac = (epsp(i) - ratep(i,1))/(ratep(i,2) - ratep(i,1))
361 yldmax(i) = yp1(i) + fac*(yp2(i)-yp1(i))
362 yldmax(i) =
max(yldmax(i),em20)
363 uvar(i,j1) = ipos1p(i)
364 uvar(i,j2) = ipos2p(i)
365
366 i1 = ii(i) +1
367 i2 = i1+1
368 fac = (-epsp(i) - raten(i,2))/(raten(i,1) - raten(i,2))
369 yldmin(i) = yn2(i) + fac*(yn1(i)-yn2(i))
370 yldmin(i) =
max(yldmin(i),em20)
371 uvar(i,i1) = ipos1n(i)
372 uvar(i,i2) = ipos2n(i)
373
374 yldelas(i)= y1(i)
375 uvar(i,1) = ipos1(i)
376 ENDDO
377
378
379
380
381
382 IF(iflag(1)==2)THEN
383
384 DO i=1,nel0
385 IF(off(i)==1) THEN
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
391 IF(svm(i)/=zero) THEN
392 IF(svm(i)>=yldelas(i))THEN
393 DO iter = 1, 10
394 r =
min(one,yldelas(i)/
max(em20,svm(i)))
395 dev = (signxx(i) + signyy(i))*third
396 signxx(i)=(signxx(i)-dev)*r
397 signyy(i)=(signyy(i)-dev)*r
398 signxy(i)=signxy(i)*r
399 signyz(i)=signyz(i)*r
400 signzx(i)=signzx(i)*r
401
402 signxx(i)=(signxx(i)+ dev)
403 signyy(i)=(signyy(i)+ dev)
404 svm(i)=sqrt(signxx(i)*signxx(i)
405 . +signyy(i)*signyy(i)
406 . -signxx(i)*signyy(i)
407 . + three*signxy(i)*signxy(i))
408 ENDDO
409 r = one
410 DO iter = 1,10
411 sigexx(i)=signxx(i) + sigvxx(i)
412 sigeyy(i)=signyy(i) + sigvyy(i)
413 sigexy(i)=signxy(i) + sigvxy(i)
414 sigeyz(i)=signyz(i) + sigvyz(i)
415 sigezx(i)=signzx(i) + sigvzx(i)
416 svm(i)=sqrt(sigexx(i)*sigexx(i)
417 . +sigeyy(i)*sigeyy(i)
418 . -sigexx(i)*sigeyy(i)
419 . + three*sigexy(i)*sigexy(i))
420 IF(svm(i)>=yldmax(i)) THEN
421 r =
min(one,yldmax(i)/
max(em20,svm(i)))
422 ELSEIF(svm(i)<=yldmin(i))THEN
423 r =
max(one,yldmin(i)/
max(em20,svm(i)))
424 ENDIF
425 dev = (sigexx(i) + sigeyy(i))*third
426 sigexx(i)=(sigexx(i)-dev)*r
427 sigeyy(i)=(sigeyy(i)-dev)*r
428 sigexy(i)=sigexy(i)*r
429 sigeyz(i)=sigeyz(i)*r
430 sigezx(i)=sigezx(i)*r
431 sigexx(i)=sigexx(i)+ dev
432 sigeyy(i)=sigeyy(i)+ dev
433 svm(i)=sqrt(sigexx(i
434 . +sigeyy(i)*sigeyy(i)
435 . -sigexx(i)*sigeyy(i)
436 . + three*sigexy(i)*sigexy(i))
437
438 IF(r/=1.) THEN
439 sigvxx(i)=sigexx(i) - signxx(i)
440 sigvyy(i)=sigeyy(i) - signyy(i)
441 sigvxy(i)=sigexy(i) - signxy(i)
442 sigvyz(i)=sigeyz(i) - signyz(i)
443 sigvzx(i)=sigezx(i) - signzx(i)
444 ELSE
445 GO TO 200
446
447 ENDIF
448 ENDDO
449 ELSE
450 ENDIF
451 ENDIF
452 ENDIF
453 200 ENDDO
454
455
456 DO i=1,nel0
457 lamda = nu(i) * e(i) /((one + nu(i))*(one - two* nu(i)))
458 lamda = lamda + dtinv*third*(three*visc(i) - two*viscv(i))
459 c11= lamda + two*g(i) + two*dtinv*viscv(i)
460 dezz = -( lamda/
max(em20,c11))*(depsxx(i) + depsyy(i))
461 thk(i) = thk(i) + dezz*thkly(i)*off(i)
462 ENDDO
463
464 ELSEIF(iflag(1)==1)THEN
465
466
467
468
469
470 ELSEIF(iflag(1)==0)THEN
471
472 ENDIF
473
474 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)