47
48
49
50#include "implicit_f.inc"
51
52
53
54#include "mvsiz_p.inc"
55#include "param_c.inc"
56#include "com01_c.inc"
57
58
59
60 INTEGER NEL, NUVAR, NGL(NEL), MAT(NEL),IPM(NPROPMI,*)
62 . time,uparam(*),
63 .
area(nel),rho0(nel),eint(nel,2),
64 . thk0(nel),pla(nel),
65 . epspxx(nel),epspyy(nel),
66 . epspxy(nel),epspyz(nel),epspzx(nel),
67 . depsxx(nel),depsyy(nel),depsxy(nel),
68 . depbxx(nel),depbyy(nel),depbxy(nel),
69 . depsyz(nel),depszx(nel),
70 . epsxx(nel) ,epsyy(nel) ,
71 . epsxy(nel) ,epsyz(nel) ,epszx(nel) ,
72 . sigoxx(nel),sigoyy(nel),sigoxy(nel),
73 . momoxx(nel),momoyy(nel),momoxy(nel),
74 . sigoyz(nel),sigozx(nel),gs(nel),shf(nel),seq_output(nel)
75
76
77
79 . signxx(nel),signyy(nel),signxy(nel),
80 . momnxx(nel),momnyy(nel),momnxy(nel),
81 . signyz(nel),signzx(nel),epsp(nel),
82 . sigvxx(nel),sigvyy(nel),sigy(nel),
83 . sigvxy(nel),sigvyz(nel),sigvzx(nel),
84 . soundsp(nel),viscmax(nel),etse(nel)
85
86
87
88 my_real uvar(nel,nuvar), off(nel),thk(nel)
89
90
91
92 INTEGER NPF(*)
94
95
96
97 INTEGER I,J,K,J1,J2,N,NINDX,NMAX,IADBUF,NFUNC,NS,
98 . NRATE,IAD1(MVSIZ),IPOS1(MVSIZ),ILEN1(MVSIZ),
99 . IAD2(MVSIZ),IPOS2(MVSIZ),ILEN2(MVSIZ),OPTE,
100 . JJ(MVSIZ),INDEX(MVSIZ),IFUNC(100), IFUNCE,MX,ISRATE
102 . e(mvsiz),nu,a,b,c,fac,dezz,s1,s2,s12,s3,
103 . dpla,epst,a1(mvsiz),a2(mvsiz),g(mvsiz),g3(mvsiz),
104 . dydx1(mvsiz),dydx2(mvsiz),rate(mvsiz,2),svm(mvsiz),
105 . yld(mvsiz),y1(mvsiz),y2(mvsiz),dr,f1,f2,
106 . yfac(mvsiz,2),nnu1,nu1(mvsiz),
107 . nu2(mvsiz),nu3(mvsiz),nu4(mvsiz),nu5(mvsiz),dpla_i,
108 . fail(mvsiz),epsmax,epsr1,epsr2,
109 . err,f,df,yld_i,dpla_j(mvsiz),
110 . aa1,aa2,aa3,fn(3),fm(3),fnm(3),sn1,sn2,sm1,sm2,
111 . pn(3),p_m(3),dfn(3),dfm(3),dfnm(3),djac(3),jac(3),c1,
112 . xp(3),pnm1(3),pnm2(3),jq(mvsiz),jq2(mvsiz),s(mvsiz),
113 . gm(mvsiz),cm(mvsiz),qtier(mvsiz),cnm(mvsiz),
114 . num1(mvsiz),num2(mvsiz),an(mvsiz),bn(mvsiz),
115 . am1(mvsiz),am2(mvsiz),gama(mvsiz),gama2(mvsiz),
116 . h(mvsiz),einf,ce,dydxe(mvsiz),
117 . escale(mvsiz)
119 . a_1,a_2,a_3,amn_xy(mvsiz),jac_i(3),jac_2(3),
120 . a01,a02,a03,a12,
121 . anxx(mvsiz),anyy(mvsiz),anxy(mvsiz),an_xy(mvsiz),
122 . amxx(mvsiz),amyy(mvsiz),amxy(mvsiz),am_xy(mvsiz),
123 . anmxx(mvsiz),anmyy(mvsiz),anmxy(mvsiz),anm_xy(mvsiz),
124 . sn11(mvsiz),sn22(mvsiz),sm11(mvsiz),sm22(mvsiz),
125 . b_1(mvsiz),b_2(mvsiz),b_3(mvsiz),q12(mvsiz),q21(mvsiz),
126 . sfn,sfm,sfnm,dsfn,dsfm,dsfnm,tol
127
128 DATA nmax/3/,ns/10/
129 tol =em6
130
131
132
133 mx = mat(1)
134 nfunc = ipm(10,mx)
135 DO j=1,nfunc
136 ifunc(j) = ipm(10+j,mx)
137 ENDDO
138 iadbuf = ipm(7,mx)-1
139 israte = ipm(3,mx)
140 nu = uparam(iadbuf+6)
141 a01 = uparam(iadbuf+7)
142 a02 = uparam(iadbuf+8)
143 a03 = uparam(iadbuf+9)
144 a12 = uparam(iadbuf+10)
145 nrate = nint(uparam(iadbuf+1))
146 epsmax=uparam(iadbuf+ns+2*nrate+1)
147
148 IF(epsmax==zero)THEN
149 IF(tf(npf(ifunc(1)+1)-1)==zero)THEN
150 epsmax=tf(npf(ifunc(1)+1)-2)
151 ELSE
152 epsmax= ep30
153 ENDIF
154 ENDIF
155 nnu1 = nu / (one - nu)
156 epsr1 =uparam(iadbuf+ns+2*nrate+2)
157 epsr2 =uparam(iadbuf+ns+2*nrate+3)
158
159 opte = uparam(iadbuf+ns+2*nrate + 10)
160 einf = uparam(iadbuf+ns+2*nrate+ 11)
161 ce = uparam(iadbuf+ns+2*nrate+ 12)
162
163 DO i=1,nel
164 e(i) = uparam(iadbuf+2)
165 a1(i) = uparam(iadbuf+3)
166 a2(i) = uparam(iadbuf+4)
167 g(i) = uparam(iadbuf+5)
168 g3(i) = three*g(i)
169
170 c1=thk0(i)*one_over_12
171 am1(i) = a1(i)*c1
172 am2(i) = a2(i)*c1
173 gm(i) = g(i)*c1
174 ENDDO
175
176 IF (opte == 1)THEN
177 ifunce = uparam(iadbuf+ns+2*nrate+ 9)
178 DO i=1,nel
179 IF(pla(i) > zero)THEN
180 escale(i) = finter(ifunc(ifunce),pla(i),npf,tf,dydxe(i))
181 e(i) = escale(i)* e(i)
182 g(i) = half*e(i)/(one+nu)
183 gs(i) = g(i)*shf(i)
184 g3(i) = three*g(i)
185 a1(i) = e(i)/(one - nu*nu)
186 a2(i) = nu*a1(i)
187 am1(i) = a1(i)*c1
188 am2(i) = a2(i)*c1
189 gm(i) = g(i)*c1
190 ENDIF
191 ENDDO
192 ELSEIF ( ce /= zero) THEN
193 DO i=1,nel
194 IF(pla(i) > zero)THEN
195 e(i) = e(i)-(e(i)-einf)*(one-exp(-ce*pla(i)))
196 g(i) = half*e(i)/(one+nu)
197 gs(i) = g(i)*shf(i)
198 g3(i) = three*g(i)
199 a1(i) = e(i)/(one - nu*nu)
200 a2(i) = nu*a1(i)
201 am1(i) = a1(i)*c1
202 am2(i) = a2(i)*c1
203 gm(i) = g(i)*c1
204 ENDIF
205 ENDDO
206 ENDIF
207
208 IF (isigi==0) THEN
209 IF(time==0.0)THEN
210 DO i=1,nel
211 uvar(i,1)=zero
212 uvar(i,2)=zero
213 DO j=1,nrate
214 uvar(i,j+2)=zero
215 ENDDO
216 ENDDO
217 ENDIF
218 ENDIF
219
220 DO i=1,nel
221 signxx(i)=sigoxx(i)+a1(i)*depsxx(i)+a2(i)*depsyy(i)
222 signyy(i)=sigoyy(i)+a2(i)*depsxx(i)+a1(i)*depsyy(i)
223 signxy(i)=sigoxy(i)+g(i) *depsxy(i)
224 momnxx(i)=momoxx(i)+am1(i)*depbxx(i)+am2(i)*depbyy(i)
225 momnyy(i)=momoyy(i)+am2
226 momnxy(i)=momoxy(i)+gm(i) *depbxy(i)
227 signyz(i)=sigoyz(i)+gs(i) *depsyz(i)
228 signzx(i)=sigozx(i)+gs(i) *depszx(i)
229
230 soundsp(i) = sqrt(a1(i)/rho0(i))
231 viscmax(i) = zero
232 etse(i) = one
233
234
235
236 IF (israte == 0) THEN
237 epsp(i) = half*( abs(epspxx(i)+epspyy(i))
238 . + sqrt( (epspxx(i)-epspyy(i))*(epspxx(i)-epspyy(i))
239 . + epspxy(i)*epspxy(i) ) )
240 ENDIF
241
242
243
244 epst = half*( epsxx(i)+epsyy(i)
245 . + sqrt( (epsxx(i
246 . + epsxy(i)*epsxy(i) ) )
247 fail(i)
248
249 ENDDO
250
251
252
253 DO i=1,nel
254 jj(i) = 1
255 ENDDO
256 DO i=1,nel
257 iadbuf = ipm(7,mat(i))-1
258 DO j=2,nrate-1
259 IF(epsp(i)>=uparam(iadbuf+ns+j)) jj(i) = j
260 ENDDO
261 rate(i,1)=uparam(iadbuf+ns+jj(i))
262 rate(i,2)=uparam(iadbuf+ns+jj(i)+1)
263 yfac(i,1)=uparam(iadbuf+ns+nrate+jj(i))
264 yfac(i,2)=uparam(iadbuf+ns+nrate+jj(i)+1)
265 ENDDO
266 DO i=1,nel
267 j1 = jj(i)
268 j2 = j1+1
269 ipos1(i) = nint(uvar(i,j1))
270 iad1(i) = npf(ifunc(j1)) / 2 + 1
271 ilen1(i) = npf(ifunc(j1)+1) / 2 - iad1(i) - ipos1(i)
272 ipos2(i) = nint(uvar(i,j2))
273 iad2(i) = npf(ifunc(j2)) / 2 + 1
274 ilen2(i) = npf(ifunc(j2)+1) / 2 - iad2(i) - ipos2(i)
275 ENDDO
276
277 CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
278 CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,y2)
279
280 DO i=1,nel
281 j1 = jj(i)
282 j2 = j1+1
283 y1(i)=y1(i)*yfac(i,1)
284 y2(i)=y2(i)*yfac(i,2)
285 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
286 yld(i) = fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
287 yld(i) =
max(yld(i),em20)
288 sigy(i) = yld(i)
289 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
290 uvar(i,j1) = ipos1(i)
291 uvar(i,j2) = ipos2(i)
292 ENDDO
293
294
295
296 DO i=1,nel
297 c1=pla(i)*e(i)
298 gama(i)=three_half*(c1+yld(i))/(three_half*c1+yld(i))
299 gama2(i)=gama(i)*gama(i)
300 cm(i)= sixteen*gama2(i)
301 cnm(i)= sqr16_3*gama(i)
302 qtier(i)=four_over_3*gama2(i)
303 h(i) =
max(zero,h(i))
304
305 s1=a01*(signxx(i)*signxx(i)+cm(i)*momnxx(i)*momnxx
306 s2=a02*(signyy(i)*signyy(i)+cm(i)*momnyy(i)*momnyy(i))
307 s3=a03*(signyy(i)*signxx(i)+cm(i)*momnxx(i)*momnyy(i))
308 anxy(i)=a12*signxy(i)*signxy(i)
309 amxy(i)=a12*momnxy(i)*momnxy(i)*cm(i)
310 f1=s1+s2-s3+anxy(i)+amxy(i)
311 s1=a01*(signxx(i)*momnxx(i))
312 s2=a02*(signyy(i)*momnyy(i))
313 s3=a03*(signxx(i)*momnyy(i)+signyy(i)*momnxx(i))*half
314 anmxy(i)=a12*signxy(i)*momnxy(i)
315 f2=cnm(i)*(s1+s2-s3+anmxy(i))
316 anmxy(i)=anmxy(i)*cnm(i)
317 svm(i)=sqrt(f1+abs(f2))
318
319 dezz =-(depsxx(i)+depsyy(i))*nnu1
320 thk(i) = thk(i)+thk(i)*dezz*off(i)
321 ENDDO
322
323
324
325 nindx=0
326 DO i=1,nel
327 IF(svm(i)>yld(i).AND.off(i)==one) THEN
328 nindx=nindx+1
329 index(nindx)=i
330 ENDIF
331 ENDDO
332 IF(nindx==0) RETURN
333
334
335
336 DO j=1,nindx
337 i=index(j)
338 nu2(i) = 1.-nu*nu
339 nu3(i) = nu*half
340 nu4(i) = half -nu3(i)
341 nu5(i) = one - nnu1
342 dpla_j(i)=(svm(i)-yld(i))/(g3(i)*qtier(i)+h(i))
343 etse(i)= h(i)/(h(i)+e(i))
344 s1=a01*nu*two - a03
345 s2=a02*nu*two - a03
346 s12=a03-nu*(a01+a02)
347 s3=sqrt(nu2(i)*(a01-a02)*(a01-a02)+s12*s12)
348 IF (abs(s1)<em20) THEN
349 q12(i)=zero
350 ELSE
351 q12(i)=-(a01-a02+s3)/s1
352 ENDIF
353 IF (abs(s2)<em20) THEN
354 q21(i)=zero
355 ELSE
356 q21(i)=(a01-a02+s3)/s2
357 ENDIF
358 jq(i)=one/(one - q12(i)*q21(i))
359 jq2(i)=jq(i)*jq(i)
360 a=a01*q12(i)
361 b=a02*q21(i)
362 a_1=(a01+a03*q21(i)+b*q21(i))*jq2(i)
363 a_2=(a02+a03*q12(i)+a*q12(i))*jq2(i)
364 a_3=(a+b)*jq2(i)*2.0+a03*(jq2(i)*2.0-jq(i))
365 sn11(i)=signxx(i)+signyy(i)*q12(i)
366 sn22(i)=q21(i)*signxx(i)+signyy(i)
367 sm11(i)=momnxx(i)+momnyy(i)*q12(i)
368 sm22(i)=q21(i)*momnxx(i)+momnyy(i)
369 anxx(i)=a_1*sn11(i)*sn11(i)
370 anyy(i)=a_2*sn22(i)*sn22(i)
371 an_xy(i)=a_3*sn11(i)*sn22(i)
372 amxx(i)=a_1*sm11(i)*sm11(i)*cm(i)
373 amyy(i)=a_2*sm22(i)*sm22(i)*cm(i)
374 am_xy(i)=a_3*sm11(i)*sm22(i)*cm(i)
375 anmxx(i)=a_1*sn11(i)*sm11(i)*cnm(i)
376 anmyy(i)=a_2*sn22(i)*sm22(i)*cnm(i)
377 anm_xy(i)=a_3*sn11(i)*sm22(i)*cnm(i)*half
378 amn_xy(i)=a_3*sm11(i)*sn22(i)*cnm(i)*half
379 a=a03*nu3(i)
380 b=s3*jq(i)
381 b_1(i)=a02-a-b
382 b_2(i)=a01-a+b
383 b_3(i)=a12*nu4(i)
384 ENDDO
385
386
387
388 DO n=1,nmax
389 DO j=1,nindx
390 i=index(j)
391 dpla_i=dpla_j(i)
392 yld_i =yld(i)+h(i)*dpla_i
393 dr =a1(i)*dpla_i/yld_i
394 xp(1) =b_1(i)*dr
395 xp(2) =b_2(i)*dr
396 xp(3) =b_3(i)*dr
397 c1=one + qtier(i)
398 b=twop444*gama2(i)
399 DO k=1,3
400 djac(k)=c1+b*xp(k)
401 jac(k)=one
402 jac_i(k)=one/jac(k)
403 jac_2(k)= jac_i(k)*jac_i(k)
404 a=xp(k)*zep444*gama2(i)
405 dfn(k)=fivep5+sixteenp5*a
406 fn(k)=one+(dfn(k)+fivep5)*a*half
407 a=onep8333*xp(k)
408 dfm(k)=onep8333+a
409 fm(k)=one +(dfm(k)*xp(k)+a)*half
410 dfnm(k)=-b*xp(k)
411 fnm(k)=one+dfnm(k)*xp(k)*half
412 ENDDO
413 a=jac_2(1)*fn(1)
414 b=jac_2(2)*fn(2)
415 c=jac_2(3)*fn(3)
416 sfn=a*(anxx(i)-an_xy(i))+b*(anyy(i)-an_xy(i))+c*anxy(i)
417 a=jac_2(1)*fm(1)
418 b=jac_2(2)*fm(2)
419 c=jac_2(3)*fm(3)
420 sfm=a*(amxx(i)-am_xy(i))+b*(amyy(i)-am_xy(i))+c*amxy(i)
421 a=jac_2(1)*fnm(1)
422 b=jac_2(2)*fnm(2)
423 c=jac_2(3)*fnm(3)
424 sfnm=a*(anmxx(i)-anm_xy(i))+b*(anmyy(i)-amn_xy(i))
425 . +c*anmxy(i)
426 c1=abs(sfnm)/
max(sfn,sfm)
427 IF (c1<tol) THEN
428 s(i)=zero
429 ELSEIF (sfnm<zero) THEN
430 s(i)=-one
431 ELSE
432 s(i)=one
433 ENDIF
434 f =sfn+sfm+s(i)*sfnm-yld_i*yld_i
435 c1=zep444*gama2(i)
436 s12=jac_2(1)*b_1(i)
437 s1=s12*c1
438 aa1=two*jac_i(1)*djac(1)*s12
439 s12=jac_2(2)*b_2(i)
440 s2=s12*c1
441 aa2=two*jac_i(2)*djac(2)*s12
442 s12=jac_2(3)*b_3(i)
443 s3=s12*c1
444 aa3=two*jac_i(3)*djac(3)*s12
445 a=s1*dfn(1)-aa1*fn(1)
446 b=s2*dfn(2)-aa2*fn(2)
447 c=s3*dfn(3)-aa3*fn(3)
448 dsfn=a*(anxx(i)-an_xy(i))+b*(anyy(i)-an_xy(i))+c*anxy(i)
449 a=s1*dfm(1)-aa1*fm(1)
450 b=s2*dfm(2)-aa2*fm(2)
451 c=s3*dfm(3)-aa3*fm(3)
452 dsfm=a*(amxx(i)-am_xy(i))+b*(amyy(i)-am_xy(i))+c*amxy(i)
453 a=s1*dfnm(1)-aa1*fnm(1)
454 b=s2*dfnm(2)-aa2*fnm(2)
455 c=s3*dfnm(3)-aa3*fnm(3)
456 dsfnm=a*(anmxx(i)-anm_xy(i))+b*(anmyy(i)-amn_xy(i))
457 . +c*anmxy(i)
458 df =(dsfn+dsfm+s(i)*dsfnm)*
459 . (a1(i)-dr*h(i))/yld_i-two*h(i)*yld_i
460 dpla_j(i)=
max(zero,dpla_i-f/df)
461 ENDDO
462 ENDDO
463
464
465
466 DO j=1,nindx
467 i=index(j)
468 pla(i) = pla(i) + dpla_j(i)
469 yld_i =yld(i)+h(i)*dpla_j(i)
470 dr =a1(i)*dpla_j(i)/yld_i
471 xp(1) =b_1(i)*dr
472 xp(2) =b_2(i)*dr
473 xp(3) =b_3(i)*dr
474 c1=one+qtier(i)
475 b=twop444*gama2(i)
476 DO k=1,3
477 a=b*xp(k)
478 jac(k)=1.+c1*xp(k)+a*xp(k)
479 jac_i(k)=one/jac(k)
480 jac_2(k)= jac_i(k)*jac_i(k)
481 fnm(k)=one-a*xp(k)*half
482 a=xp(k)*jac_i(k)
483 pn(k)=jac_i(k)+qtier(i)*a
484 p_m(k)=jac_i(k)+a
485 pnm1(k)=-sqr4_3*gama(i)*a
486 pnm2(k)=pnm1(k)*one_over_12
487 ENDDO
488 a=jac_2(1)*fnm(1)
489 b=jac_2(2)*fnm(2)
490 c=jac_2(3)*fnm(3)
491 sfnm=a*(anmxx(i)-anm_xy(i))+b*(anmyy(i)-amn_xy(i))
492 . +c*anmxy(i)
493 sn1=sn11(i)*pn(1)+sm11(i)*pnm1(1)*s(i)
494 sn2=sn22(i)*pn(2)+sm22(i)*pnm1(2)*s(i)
495 s3=signxy(i)*pn(3)+momnxy(i)*pnm1(3)*s(i)
496 sm1=sm11(i)*p_m(1)+sn11(i)*pnm2(1)*s(i)
497 sm2=sm22(i)*p_m(2)+sn22(i)*pnm2(2)*s(i)
498 momnxy(i)=momnxy(i)*p_m(3)+signxy(i)*pnm2(3)*s(i)
499 signxx(i)=jq(i)*(sn1-sn2*q12(i))
500 signyy(i)=jq(i)*(sn2-sn1*q21(i))
501 signxy(i)=s3
502 momnxx(i)=jq(i)*(sm1-sm2*q12(i))
503 momnyy(i)=jq(i)*(sm2-sm1*q21(i))
504 s1=a01*signxx(i)+a02*signyy(i)
505 . -a03*(signxx(i)+signyy(i))*half
506 dezz = - nu5(i)*dpla_j(i)*s1/yld_i
507 thk(i) = thk(i)+thk(i)* dezz*off(i)
508 ENDDO
509
510 DO i=1,nel
511 IF(pla(i)>epsmax.AND.off(i)==one)off(i)=four_over_5
512 ENDDO
513
514 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)