44
45
46
47#include "implicit_f.inc"
48
49
50
51#include "mvsiz_p.inc"
52
53
54
55 INTEGER NEL,NPT,IPT,IFLA,NUPARAM, NUVAR,NGL(NEL),ISRATE
56
58 . time , timestep , uparam(nuparam),
59 . rho0(nel), eint(nel,2),
60 . epspxx(nel), epspyy(nel),shf(nel),
61 . epspxy(nel), epspyz(nel), epspzx(nel),
62 . depsxx(nel), depsyy(nel),
63 . depsxy(nel), depsyz(nel), depszx(nel),
64 . epsxx(nel), epsyy(nel),
65 . epsxy(nel), epsyz(nel), epszx(nel),
66 . sigoxx(nel), sigoyy(nel),
67 . sigoxy(nel), sigoyz(nel), sigozx(nel),
68 . thkly(nel) ,thk(nel),
area(nel),
69 . asrate
70
71
72
74 . signxx(nel), signyy(nel),
75 . signxy(nel), signyz(nel), signzx(nel),
76 . sigvxx(nel), sigvyy(nel),
77 . sigvxy(nel), sigvyz(nel), sigvzx(nel),
78 . soundsp(nel), viscmax(nel),etse(nel)
79
80
81
82 my_real uvar(nel,nuvar), off(nel)
83 my_real ,
INTENT(INOUT) :: epsd(nel)
84
85
86
87 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
89 . finter,tf(*)
90 EXTERNAL finter
91
92
93
94 INTEGER I,KF,IFLAG,ICORRECT
95
97 . poisson,e,e1,e2,g02,
98 . gt2,bulkt,relvexp,
99 . vmu2,vbulk,fac,
100 . c1,c2,c3,pmin,dpdmu,epspc,
101 . p0,phi,gama0, var,epsp,
102 . gama(mvsiz),amu(mvsiz),
103 . sm(mvsiz),em(mvsiz),dedtm(mvsiz),
104 . g2(mvsiz),bulk(mvsiz),
105 . dsxx(mvsiz),dsyy(mvsiz),dsxy(mvsiz),
106 . dexx(mvsiz),deyy(mvsiz),dezz,
107 . dexy(mvsiz),epspzz(mvsiz),epszz(mvsiz),
108 . dedtxx(mvsiz),dedtyy(mvsiz),dedtxy(mvsiz),
109 . dsdtxx(mvsiz),dsdtyy(mvsiz),dsdtxy(mvsiz),
110 . dpdro(mvsiz),p(mvsiz),pdot(mvsiz),
111 . mg2(mvsiz),pg2(mvsiz),mk(mvsiz),pk(mvsiz),
112 . sigair(mvsiz),relvol(mvsiz),enew(mvsiz),
113 . rho(mvsiz)
114
116 . ssm,epsm,epspm,small,cshear,midstep,dt05
117
118 small = em3
119 cshear= zep426667
120
121 IF (time == zero) THEN
122 DO i=1,nel
123 uvar(i,3) =
area(i)*thk(i)
124 ENDDO
125 ENDIF
126
127
128 poisson= uparam(2)
129 e = uparam(1)
130 g02 = half*e/(1+poisson)
131 e1 = uparam(3)
132 e2 = uparam(4)
133 gt2 = uparam(5)/(one+uparam(6))
134 bulkt = uparam(5)/(one-two*uparam(6))
135 vmu2 = two*uparam(7)
136 vbulk = three*uparam(8)+vmu2
137
138 p0 = uparam(9)
139 phi = uparam(10)
140 gama0 = uparam(11)
141
142 c1 = uparam(12)
143 c2 = uparam(13)
144 c3 = uparam(14)
145 iflag = uparam(15)
146 pmin = uparam(16)
147 relvexp = uparam(17)
148 fac = uparam(18)
149
150 kf = ifunc(1)
151
152 icorrect=0
153 IF(nuparam>=19) icorrect=nint(uparam(19))
154
155 IF(icorrect == 0)THEN
156 dt05=half*timestep
157 ELSE
158 dt05=zero
159 END IF
160
161 DO i=1,nel
162 epsxx(i)=epsxx(i)-dt05*epspxx(i)
163 epsyy(i)=epsyy(i)-dt05*epspyy(i)
164 epsxy(i)=epsxy(i)-dt05*epspxy(i)
165 epsyz(i)=epsyz(i)-dt05*epspyz(i)
166 epszx(i)=epszx(i)-dt05*epspzx(i)
167 END DO
168
169 DO 200 i=1,nel
170
171 epspm = epspxx(i)+epspyy(i)
172 epspc = epspxx(i)-epspyy(i)
173 epsp = sqrt(half*(epspc*epspc+epspxx(i)*epspxx(i)
174 1 + epspyy(i)*epspyy(i)) +three_half* (epspxy(i)*epspxy(i)
175 2 + epspyz(i)*epspyz(i)+epspzx(i)*epspzx(i)))
176 IF (israte > 0) THEN
177 epsd(i) = asrate*epsp + (one - asrate)*uvar(i,4)
178 uvar(i,4) = epsd(i)
179 ELSE
180 epsd(i) = epsp
181 ENDIF
182
183
184 rho(i)=uvar(i,3)*rho0(i)/
area(i)/thk(i)
185 relvol(i)=rho0(i)/rho(i)
186 enew(i)=(
max(e,e1*epsp+e2))/(exp(relvexp*log(relvol(i))))
187 g2(i)=enew(i)/(one+poisson)
188 bulk(i)=enew(i)/(one-two*poisson)
189 mk(i) = (bulk(i)+bulkt)/vbulk
190 pk(i) = bulk(i)*bulkt/vbulk
191 mg2(i) = (g2(i)+gt2)/vmu2
192 pg2(i) = g2(i)*gt2/vmu2
193 gama(i) = (relvol(i)-one+gama0)
194 IF(one+gama(i)-phi<=small) gama(i)=-(one-phi-small)
195 ssm=sigoxx(i)+sigoyy(i)
196 sm(i)=third*ssm
197 epsm =epsxx(i)+epsyy(i)
198 var=two*pg2(i)+c3*pk(i)
199 epspzz(i)=((g2(i)-c1*bulk(i))*epspm-(mg2(i)-c2*mk(i))*ssm
200 1 +(pg2(i)-c3*pk(i))*epsm-var*uvar(i,2))
201 2 /(g2(i)+c1*bulk(i)+g2(i)+var*timestep)
202 dezz=epspzz(i)*timestep
203 epszz(i)=uvar(i,2)+dezz
204 uvar(i,2)=epszz(i)
205 thk(i) = thk(i) + dezz*thkly(i)*off(i)
206 em(i)=third*(epsm+epszz(i))
207 dedtm(i)=third*(epspm+epspzz(i))
208 sigair(i)=
max(zero,-(p0*gama(i))/(1+gama(i)-phi))
209 200 CONTINUE
210
211 IF(icorrect/=0)THEN
212 DO 210 i=1,nel
213
214 dsxx(i)=sigoxx(i)-sm(i)
215 dsyy(i)=sigoyy(i)-sm(i)
216 dsxy(i)=sigoxy(i)
217
218
219 dexx(i)=epsxx(i)-em(i)
220 deyy(i)=epsyy(i)-em(i)
221 dexy(i)=epsxy(i)
222
223
224 dedtxx(i)=epspxx(i)-dedtm(i)
225 dedtyy(i)=epspyy(i)-dedtm(i)
226 dedtxy(i)=epspxy(i)
227 210 CONTINUE
228 ELSE
229 DO i=1,nel
230
231 dsxx(i)=sigoxx(i)-sm(i)
232 dsyy(i)=sigoyy(i)-sm(i)
233 dsxy(i)=sigoxy(i)
234
235
236 dexx(i)=epsxx(i)-em(i)
237 deyy(i)=epsyy(i)-em(i)
238 dexy(i)=epsxy(i)*half
239
240
241 dedtxx(i)=epspxx(i)-dedtm(i)
242 dedtyy(i)=epspyy(i)-dedtm(i)
243 dedtxy(i)=epspxy(i)*half
244 END DO
245 END IF
246
247
248
249 DO 250 i=1,nel
250 dsdtxx(i)=g2(i)*dedtxx(i)-mg2(i)*dsxx(i)+pg2(i)*dexx(i)
251 dsdtyy(i)=g2(i)*dedtyy(i)-mg2(i)*dsyy(i)+pg2(i)*deyy(i)
252 dsdtxy(i)=g2(i)*dedtxy(i)-mg2(i)*dsxy(i)+pg2(i)*dexy(i)
253 midstep=one/(one+mg2(i)*dt05)
254 dsdtxx(i)=dsdtxx(i)*midstep
255 dsdtyy(i)=dsdtyy(i)*midstep
256 dsdtxy(i)=dsdtxy(i)*midstep
257 250 CONTINUE
258
259
260 DO 260 i=1,nel
261 sm(i)=sm(i)+uvar(i,1)
262 IF(kf/=0) THEN
263 amu(i)=rho(i)/rho0(i)-one
264 IF(iflag == 0) THEN
265 p(i)=-fac*finter(kf,amu(i),npf,tf,dpdmu)
266 ELSE
267 pdot(i)=( c1*bulk(i)*dedtm(i)
268 & -c2*mk(i)*sm(i)
269 & +c3*pk(i)*em(i))
270 & /(one+c2*dt05*mk(i))
271 p(i)=sm(i)+pdot(i)*timestep
272 & -fac*finter(kf,amu(i),npf,tf,dpdmu)
273 ENDIF
274 ELSE
275 pdot(i)=( c1*bulk(i)*dedtm(i)
276 & -c2*mk(i)*sm(i)
277 & +c3*pk(i)*em(i))
278 & /(one+c2*dt05*mk(i))
279 p(i)=sm(i)+pdot(i)*timestep
280 ENDIF
281 IF(p(i)<=pmin) p(i)=pmin
282 p(i)=p(i)-sigair(i)
283 260 CONTINUE
284
285
286 DO 270 i=1,nel
287 dpdro(i)= g2(i)/(one-poisson)
288 & +p0*(one-phi)/(one+gama(i)-phi)**2
289 270 CONTINUE
290
291 DO 280 i=1,nel
292 soundsp(i)=sqrt(dpdro(i)/rho(i))
293 280 CONTINUE
294
295
296
297 IF(icorrect/=0)THEN
298 DO 290 i=1,nel
299 signxx(i)=dsxx(i)+dsdtxx(i)*timestep+p(i)
300 signyy(i)=dsyy(i)+dsdtyy(i)*timestep+p(i)
301 signxy(i)=dsxy(i)+dsdtxy(i)*timestep*0.5
302 signyz(i)=sigoyz(i)+g02*depsyz(i)*cshear
303 signzx(i)=sigozx(i)+g02*depszx(i)*cshear
304 viscmax(i)= zero
305 uvar(i,1)=sigair(i)
306 290 CONTINUE
307 ELSE
308 DO i=1,nel
309 signxx(i)=dsxx(i)+dsdtxx(i)*timestep+p(i)
310 signyy(i)=dsyy(i)+dsdtyy(i)*timestep+p(i)
311 signxy(i)=dsxy(i)+dsdtxy(i)*timestep
312 signyz(i)=sigoyz(i)+g02*depsyz(i)*cshear
313 signzx(i)=sigozx(i)+g02*depszx(i)*cshear
314 viscmax(i)= zero
315 uvar(i,1)=sigair(i)
316 END DO
317 END IF
318
319 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)