33 2 TF ,TIME ,UPARAM ,RHO0 ,
35 4 EPSPXX ,EPSPYY ,EPSPXY ,EPSPYZ ,EPSPZX ,
36 5 DEPSXX ,DEPSYY ,DEPSXY ,DEPSYZ ,DEPSZX ,
37 5 DEPBXX ,DEPBYY ,DEPBXY ,
38 6 EPSXX ,EPSYY ,EPSXY ,EPSYZ ,EPSZX ,
39 7 SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
40 7 MOMOXX ,MOMOYY ,MOMOXY ,
41 8 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
42 8 MOMNXX ,MOMNYY ,MOMNXY ,
43 9 SIGVXX ,SIGVYY ,SIGVXY ,SIGVYZ ,SIGVZX ,
44 A SOUNDSP,VISCMAX,THK ,PLA ,UVAR ,
45 B OFF ,NGL ,IPM ,MAT ,ETSE ,
46 C GS ,SIGY ,SHF ,SEQ_OUTPUT,EPSP )
50#include "implicit_f.inc"
60 INTEGER NEL, NUVAR, NGL(NEL), MAT(NEL),IPM(NPROPMI,*)
63 . AREA(NEL),RHO0(NEL),EINT(NEL,2),
65 . EPSPXX(NEL),EPSPYY(NEL),
66 . EPSPXY(NEL),EPSPYZ(NEL),EPSPZX(NEL),
67 . DEPSXX(),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)
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)
88 my_real uvar(nel,nuvar), off(nel),thk(nel)
97 INTEGER I,J,,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,C,FAC,,,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,
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),
119 . a_1,a_2,a_3,amn_xy(mvsiz),jac_i(3),jac_2(3),
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
136 ifunc(j) = ipm(10+j,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)
149 IF(tf(npf(ifunc(1)+1)-1)==zero)
THEN
150 epsmax=tf(npf(ifunc(1)+1)-2)
155 nnu1 = nu / (one - nu)
156 epsr1 =uparam(iadbuf+ns+2*nrate+2)
157 epsr2 =uparam(iadbuf+ns+2*nrate+3)
159 opte = uparam(iadbuf+ns+2*nrate + 10)
160 einf = uparam(iadbuf+ns+2*nrate+ 11)
161 ce = uparam(iadbuf+ns+2*nrate+ 12)
164 e(i) = uparam(iadbuf+2)
165 a1(i) = uparam(iadbuf+3)
166 a2(i) = uparam(iadbuf+4)
167 g(i) = uparam(iadbuf+5)
170 c1=thk0(i)*one_over_12
177 ifunce = uparam(iadbuf+ns+2*nrate+ 9)
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)
185 a1(i) = e(i)/(one - nu*nu)
192 ELSEIF ( ce /= zero)
THEN
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)
199 a1(i) = e(i)/(one - nu*nu)
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(i)*depbxx(i)+am1(i)*depbyy(i)
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)
230 soundsp(i) = sqrt(a1(i)/rho0(i))
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) ) )
244 epst = half*( epsxx(i)+epsyy(i)
245 . + sqrt( (epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
246 . + epsxy(i)*epsxy(i) ) )
247 fail(i) =
max(zero,
min(one,(epsr2-epst)/(epsr2-epsr1)))
259 IF(epsp(i)>=uparam(iadbuf+ns+j)) jj(i) = j
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)
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)
277 CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
278 CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,y2)
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)
289 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
290 uvar(i,j1) = ipos1(i)
291 uvar(i,j2) = ipos2(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)
302 qtier(i)=four_over_3*gama2(i)
303 h(i) =
max(zero,h(i))
305 s1=a01*(signxx(i)*signxx(i)+cm(i)*momnxx(i)*momnxx(i))
306 s2=a02*(signyy(i)*signyy(i)+cm(i)*momnyy(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))
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))
319 dezz =-(depsxx(i)+depsyy(i))*nnu1
320 thk(i) = thk(i)+thk(i)*dezz*off(i)
327 IF(svm(i)>yld(i).AND.off(i)==one)
THEN
340 nu4(i) = half -nu3(i)
342 dpla_j(i)=(svm(i)-yld(i))/(g3(i)*qtier(i)+h(i))
343 etse(i)= h(i)/(h(i)+e(i))
347 s3=sqrt(nu2(i)*(a01-a02)*(a01-a02)+s12*s12)
348 IF (abs(s1)<em20)
THEN
351 q12(i)=-(a01-a02+s3)/s1
353 IF (abs(s2)<em20)
THEN
356 q21(i)=(a01-a02+s3)/s2
358 jq(i)=one/(one - q12(i)*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
392 yld_i =yld(i)+h(i)*dpla_i
393 dr =a1(i)*dpla_i/yld_i
401 jac(k)=one+(djac(k)+c1)*xp(k)*half
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
409 fm(k)=one +(dfm(k)*xp(k)+a)*half
411 fnm(k)=one+dfnm(k)*xp(k)*half
416 sfn=a*(anxx(i)-an_xy(i))+b*(anyy(i)-an_xy(i))+c*anxy(i)
420 sfm=a*(amxx(i)-am_xy(i))+b*(amyy(i)-am_xy(i))+c*amxy(i)
424 sfnm=a*(anmxx(i)-anm_xy(i))+b*(anmyy(i)-amn_xy(i))
426 c1=abs(sfnm)/
max(sfn,sfm)
429 ELSEIF (sfnm<zero)
THEN
434 f =sfn+sfm+s(i)*sfnm-yld_i*yld_i
438 aa1=two*jac_i(1)*djac(1)*s12
441 aa2=two*jac_i(2)*djac(2)*s12
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))
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)
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
478 jac(k)=1.+c1*xp(k)+a*xp(k)
480 jac_2(k)= jac_i(k)*jac_i(k)
481 fnm(k)=one-a*xp(k)*half
483 pn(k)=jac_i(k)+qtier(i)*a
485 pnm1(k)=-sqr4_3*gama(i)*a
486 pnm2(k)=pnm1(k)*one_over_12
491 sfnm=a*(anmxx(i)-anm_xy(i))+b*(anmyy(i)-amn_xy(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))
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)
511 IF(pla(i)>epsmax.AND.off(i)==one)off(i)=four_over_5