44
45
46
47#include "implicit_f.inc"
48
49
50
51#include "mvsiz_p.inc"
52
53
54
55 INTEGER ,NPT,IPT,,NUPARAM, ,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,J,,IFLAG,ICORRECT
95
97 . poisson,e,e1,e2,g02,bulk0,
98 . gt2,bulkt,lamda,lamdat,relvexp,
99 . vmu2,vlamda,vbulk,fac,
100 . c1,c2,c3,pmin,dpdmu,epspc,
101 . p0,phi,gama0, var,epsp,
102 . dpdgama(mvsiz),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,epin,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 epsp = asrate*epsp + (one - asrate)*epsd(i)
178 ENDIF
179 epsd(i) = epsp
180
181
182 rho(i)=uvar(i,3)*rho0(i)/
area(i)/thk(i)
183 relvol(i)=rho0(i)/rho(i)
184 enew(i)=(
max(e,e1*epsp+e2))/(exp(relvexp*log(relvol(i))))
185 g2(i)=enew(i)/(one+poisson)
186 bulk(i)=enew(i)/(one-two*poisson)
187 mk(i) = (bulk(i)+bulkt)/vbulk
188 pk(i) = bulk(i)*bulkt/vbulk
189 mg2(i) = (g2(i)+gt2)/vmu2
190 pg2(i) = g2(i)*gt2/vmu2
191 gama(i) = (relvol(i)-one+gama0)
192 IF(one+gama(i)-phi<=small) gama(i)=-(one-phi-small)
193 ssm=sigoxx(i)+sigoyy(i)
194 sm(i)=third*ssm
195 epsm =epsxx(i)+epsyy(i)
196 var=two*pg2(i)+c3*pk(i)
197 epspzz(i)=((g2(i)-c1*bulk(i))*epspm-(mg2(i)-c2*mk(i))*ssm
198 1 +(pg2(i)-c3*pk(i))*epsm-var*uvar(i,2))
199 2 /(g2(i)+c1*bulk(i)+g2(i)+var*timestep)
200 dezz=epspzz(i)*timestep
201 epszz(i)=uvar(i,2)+dezz
202 uvar(i,2)=epszz(i)
203 thk(i) = thk(i) + dezz*thkly(i)*off(i)
204 em(i)=third*(epsm+epszz(i))
205 dedtm(i)=third*(epspm+epspzz(i))
206 sigair(i)=
max(zero,-(p0*gama(i))/(1+gama(i)-phi))
207 200 CONTINUE
208
209 IF(icorrect/=0)THEN
210 DO 210 i=1,nel
211
212 dsxx(i)=sigoxx(i)-sm(i)
213 dsyy(i)=sigoyy(i)-sm(i)
214 dsxy(i)=sigoxy(i)
215
216
217 dexx(i)=epsxx(i)-em(i)
218 deyy(i)=epsyy(i)-em(i)
219 dexy(i)=epsxy(i)
220
221
222 dedtxx(i)=epspxx(i)-dedtm(i)
223 dedtyy(i)=epspyy(i)-dedtm(i)
224 dedtxy(i)=epspxy(i)
225 210 CONTINUE
226 ELSE
227 DO i=1,nel
228
229 dsxx(i)=sigoxx(i)-sm(i)
230 dsyy(i)=sigoyy(i)-sm(i)
231 dsxy(i)=sigoxy(i)
232
233
234 dexx(i)=epsxx(i)-em(i)
235 deyy(i)=epsyy(i)-em(i)
236 dexy(i)=epsxy(i)*half
237
238
239 dedtxx(i)=epspxx(i)-dedtm(i)
240 dedtyy(i)=epspyy(i)-dedtm(i)
241 dedtxy(i)=epspxy(i)*half
242 END DO
243 END IF
244
245
246
247 DO 250 i=1,nel
248 dsdtxx(i)=g2(i)*dedtxx(i)-mg2(i)*dsxx(i)+pg2(i)*dexx(i)
249 dsdtyy(i)=g2(i)*dedtyy(i)-mg2(i)*dsyy(i)+pg2(i)*deyy(i)
250 dsdtxy(i)=g2(i)*dedtxy(i)-mg2(i)*dsxy(i)+pg2(i)*dexy(i)
251 midstep=one/(one+mg2(i)*dt05)
252 dsdtxx(i)=dsdtxx(i)*midstep
253 dsdtyy(i)=dsdtyy(i)*midstep
254 dsdtxy(i)=dsdtxy(i)*midstep
255 250 CONTINUE
256
257
258 DO 260 i=1,nel
259 sm(i)=sm(i)+uvar(i,1)
260 IF(kf/=0) THEN
261 amu(i)=rho(i)/rho0(i)-one
262 IF(iflag == 0) THEN
263 p(i)=-fac*finter(kf,amu(i),npf
264 ELSE
265 pdot(i)=( c1*bulk(i)*dedtm(i)
266 & -c2*mk(i)*sm(i)
267 & +c3*pk(i)*em(i))
268 & /(one+c2*dt05*mk(i))
269 p(i)=sm(i)+pdot(i)*timestep
270 & -fac*finter(kf,amu(i),npf,tf,dpdmu)
271 ENDIF
272 ELSE
273 pdot(i)=( c1*bulk(i)*dedtm(i)
274 & -c2*mk(i)*sm(i)
275 & +c3*pk(i)*em(i))
276 & /(one+c2*dt05*mk(i))
277 p(i)=sm(i)+pdot(i)*timestep
278 ENDIF
279 IF(p(i)<=pmin) p(i)=pmin
280 p(i)=p(i)-sigair(i)
281 260 CONTINUE
282
283
284 DO 270 i=1,nel
285 dpdro(i)= g2(i)/(one-poisson)
286 & +p0*(one-phi)/(one+gama(i)-phi)**2
287 270 CONTINUE
288
289 DO 280 i=1,nel
290 soundsp(i)=sqrt(dpdro(i)/rho(i))
291 280 CONTINUE
292
293
294
295 IF(icorrect/=0)THEN
296 DO 290 i=1,nel
297 signxx(i)=dsxx(i)+dsdtxx(i)*timestep+p(i)
298 signyy(i)=dsyy(i)+dsdtyy(i)*timestep+p(i)
299 signxy(i)=dsxy(i)+dsdtxy(i)*timestep*0.5
300 signyz(i)=sigoyz(i)+g02*depsyz(i)*cshear
301 signzx(i)=sigozx(i)+g02*depszx(i)*cshear
302 viscmax(i)= zero
303 uvar(i,1)=sigair(i)
304 290 CONTINUE
305 ELSE
306 DO i=1,nel
307 signxx(i)=dsxx(i)+dsdtxx(i)*timestep+p(i)
308 signyy(i)=dsyy(i)+dsdtyy(i)*timestep+p(i)
309 signxy(i)=dsxy(i)+dsdtxy(i)*timestep
310 signyz(i)=sigoyz(i)+g02*depsyz(i)*cshear
311 signzx(i)=sigozx(i)+g02*depszx(i)*cshear
312 viscmax(i)= zero
313 uvar(i,1)=sigair(i)
314 END DO
315 END IF
316
317 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)