44
45
46
47#include "implicit_f.inc"
48
49
50
51#include "mvsiz_p.inc"
52
53
54
55#include "scr05_c.inc"
56#include "impl1_c.inc"
57
58
59
60 INTEGER NEL, NUPARAM, NUVAR,ISMSTR,NGL(*),IHET,IEXPAN
62 . time , timestep , uparam(nuparam),
63 . rho(nel), rho0(nel), volume(nel), eint(nel),
64 . epspxx(nel), epspyy(nel), epspzz(nel),
65 . epspxy(nel), epspyz(nel), epspzx(nel),
66 . depsxx(nel), depsyy(nel), depszz(nel),
67 . depsxy(nel), depsyz(nel), depszx(nel),
68 . epsxx(nel), epsyy(nel), epszz(nel),
69 . epsxy(nel), epsyz(nel), epszx(nel),
70 . sigoxx(nel), sigoyy(nel), sigozz(nel),
71 . sigoxy(nel), sigoyz(nel), sigozx(nel),offg(nel),
72 . epsth3(nel)
73
74
75
77 . signxx(nel), signyy(nel), signzz(nel),
78 . signxy(nel), signyz(nel), signzx(nel),
79 . sigvxx(nel), sigvyy(nel), sigvzz(nel),
80 . sigvxy(nel), sigvyz(nel), sigvzx(nel),
81 . soundsp(nel), viscmax(nel), et(nel)
82
83
84
86 . uvar(nel,nuvar), off(nel) ,ww(nel)
87
88
89
90 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
91 my_real finter,fintte,tf(*),fint2v
92 EXTERNAL finter,fintte
93
94
95
96 INTEGER I,J,MULLINS,II
97
99 . mu,d,lam,g,rbulk,aa,bb,cc,p,trace(mvsiz),c(5),beta,
100 . t1(mvsiz), t2(mvsiz),t3(mvsiz),av(mvsiz,6),ev(mvsiz,3),
101 . evv(mvsiz,3),rv(mvsiz),rvd,di1lam1(mvsiz),di1lam2(mvsiz),di1lam3(mvsiz),
102 . dirprv(mvsiz,3,3),evd(3),eti(mvsiz)
104 . clam(5),lam_2(3),lam_4(3),amax,clp
105 my_real ,
DIMENSION(NEL,3) :: evm,cii
106 my_real ,
DIMENSION(NEL) :: gtmax,rkmax,rv_1
107
108
109 mu = uparam(1)
110 d = uparam(2)
111 lam = uparam(3)
112 g = uparam(4)
113 rbulk = uparam(5)
114 c(1) = uparam(6)
115 c(2) = uparam(7)
116 c(3) = uparam(8)
117 c(4) = uparam(9)
118 c(5) = uparam(10)
119 mullins = int(uparam(13))
120
121 IF(time==zero)THEN
122 DO j = 1,nuvar
123 DO i = 1, nel
124 uvar(i,j) = zero
125 ENDDO
126 ENDDO
127 ENDIF
128
129 DO i=1,nel
130 av(i,1)=epsxx(i)
131 av(i,2)=epsyy(i)
132 av(i,3)=epszz(i)
133 av(i,4)=epsxy(i)/2
134 av(i,5)=epsyz(i)/2
135 av(i,6)=epszx(i)/2
136 ENDDO
137
138
139 IF (iresp==1) THEN
141 ELSE
143 ENDIF
144
145
146
147
148 IF(ismstr==0.OR.ismstr==2.OR.ismstr==4) THEN
149 DO i=1,nel
150
151 ev(i,1)=exp(evv(i,1))
152 ev(i,2)=exp(evv(i,2))
153 ev(i,3)=exp(evv(i,3))
154 ENDDO
155 ELSEIF(ismstr==10.OR.ismstr==12) THEN
156 DO i =1,nel
157 IF(offg(i)<=one) THEN
158 ev(i,1)=sqrt(evv(i,1)+ one)
159 ev(i,2)=sqrt(evv(i,2)+ one)
160 ev(i,3)=sqrt(evv(i,3)+ one)
161 ELSE
162 ev(i,1)=evv(i,1)+ one
163 ev(i,2)=evv(i,2)+ one
164 ev(i,3)=evv(i,3)+ one
165 END IF
166 ENDDO
167 ELSE
168
169 DO i=1,nel
170 ev(i,1)=evv(i,1)+ one
171 ev(i,2)=evv(i,2)+ one
172 ev(i,3)=evv(i,3)+ one
173 ENDDO
174 ENDIF
175
176 IF (impl_s > 0 .OR. ihet > 1) THEN
177 eti(1:nel)=zero
178 DO j= 1,5
179 bb = one/lam**(2*j - 2)
180 aa = j*(j-1)*c(j)
181 DO i=1,nel
182 trace(i) = ev(i,1)*ev(i,1)+ ev(i,2)*ev(i,2) + ev(i,3)*ev(i,3)
183 cc = aa*bb*trace(i)**(j-2)
184 eti(i) = eti(i) + cc
185 ENDDO
186 ENDDO
187 DO i=1,nel
188 eti(i) = mu*eti(i)
189 et(i)= eti(i)
190 et(i)=
max(one,et(i)/g)
191 ENDDO
192 ENDIF
193
194
195
196 DO i = 1,nel
197
198 rv(i) = ev(i,1)*ev(i,2)*ev(i,3)
199 ENDDO
200
201 IF(iexpan > 0.AND.(ismstr==10.OR.ismstr==11.OR.ismstr==12)) THEN
202 DO i=1,nel
203 rv(i) = rv(i) -epsth3(i)
204 ENDDO
205 ENDIF
206
207
208
209 DO i = 1,nel
210
211
212
213 IF(rv(i)>zero) THEN
214 rvd = exp((-third)*log(rv(i)))
215 ELSE
216 rvd = zero
217 ENDIF
218 evd(1) = ev(i,1)*rvd
219 evd(2) = ev(i,2)*rvd
220 evd(3) = ev(i,3)*rvd
221 evm(i,1:3) = evd(1:3)
222
223 trace(i) = evd(1)**2 + evd(2)**2 + evd(3)**2
224 di1lam1(i) = two*evd(1)**2 - two_third*trace(i)
225 di1lam2(i) = two*evd(2)**2 - two_third*trace(i)
226 di1lam3(i) = two*evd(3)**2 - two_third*trace(i)
227 t1(i) = zero
228 t2(i) = zero
229 t3(i) = zero
230 ENDDO
231 DO j=1,5
232 bb = one/lam**(2*j - 2)
233 aa = j*c(j)
234 clam(j) = bb*c(j)
235 DO i=1,nel
236 cc = aa*bb*trace(i)**(j-1)
237 t1(i) = t1(i) + cc*di1lam1(i)
238 t2(i) = t2(i) + cc*di1lam2(i)
239 t3(i) = t3(i) + cc*di1lam3(i)
240 ENDDO
241 ENDDO
242
243 DO i=1,nel
244 rv_1(i) = one/rv(i)
245 p = rv_1(i)*rbulk*(rv(i)**2 - one)/two
246 t1(i) = (mu*t1(i) + p)*rv_1(i)
247 t2(i) = (mu*t2(i) + p)*rv_1(i)
248 t3(i) = (mu*t3(i) + p)*rv_1(i)
249 ENDDO
250
251 IF(mullins == 1)THEN
252 DO i = 1,nel
253 beta = one/lam**2
254 ww(i) = mu *(c(1)*(trace(i)- three)+c(2) * beta *(trace(i)**2- nine)
255 . + c(3) * beta**2 *(trace(i)**3- three**3)
256 . + c(4) * beta**3 *(trace(i)**4- three**4)
257 . + c(5) * beta**4 *(trace(i)**5- three**5) )
258
259 ENDDO
260 ENDIF
261
262 gtmax(1:nel) = g
263 rkmax(1:nel) = rbulk/two
264 cii(1:nel,1:3) = zero
265 DO ii = 1,5
266 clp = four*ii*clam(ii)
267 DO i=1,nel
268 lam_2(1:3) = evm(i,1:3)**2
269 lam_4(1:3) = lam_2(1:3)**2
270 aa = one_over_9*ii*trace(i)**ii
271 bb = zero
272 cc = zero
273 IF (ii>1) bb =third*(3-ii)*trace(i)**(ii-1)
274 IF (ii>2) cc =(ii-1)*trace(i)**(ii-2)
275 cii(i,1:3) = cii(i,1:3) +clp*(aa+bb*lam_2(1:3) +cc*lam_4(1:3))
276 ENDDO
277 ENDDO
278 DO i=1,nel
279 amax=
max(cii(i,1),cii(i,2),cii(i,3))
280
281 eti(i) =
max(one,amax*0.81)
282 gtmax(i) = g*eti(i)
283 rkmax(i) = rkmax(i)*(one+rv_1(i)*rv_1(i))
284 rkmax(i) =
max(rbulk,rkmax(i))
285
287 ENDDO
288
289
290 DO i=1,nel
291 signxx(i) = dirprv(i,1,1)*dirprv(i,1,1)*t1(i)
292 . + dirprv(i,1,2)*dirprv(i,1,2)*t2(i)
293 . + dirprv(i,1,3)*dirprv(i,1,3)*t3(i)
294
295 signyy(i) = dirprv(i,2,2)*dirprv(i,2,2)*t2(i)
296 . + dirprv(i,2,3)*dirprv(i,2,3)*t3(i)
297 . + dirprv(i,2,1)*dirprv(i,2,1)*t1(i)
298
299 signzz(i) = dirprv(i,3,3)*dirprv(i,3,3)*t3(i)
300 . + dirprv(i,3,1)*dirprv(i,3,1)*t1(i)
301 . + dirprv(i,3,2)*dirprv(i,3,2)*t2(i)
302
303 signxy(i) = dirprv(i,1,1)*dirprv(i,2,1)*t1(i)
304 . + dirprv(i,1,2)*dirprv(i,2,2)*t2(i)
305 . + dirprv(i,1,3)*dirprv(i,2,3)*t3(i)
306
307 signyz(i) = dirprv(i,2,2)*dirprv(i,3,2)*t2(i)
308 . + dirprv(i,2,3)*dirprv(i,3,3)*t3(i)
309 . + dirprv(i,2,1)*dirprv(i,3,1)*t1(i)
310
311 signzx(i) = dirprv(i,3,3)*dirprv(i,1,3)*t3(i)
312 . + dirprv(i,3,1)*dirprv(i,1,1)*t1(i)
313 . + dirprv(i,3,2)*dirprv(i,1,2)*t2(i)
314
315
316 soundsp(i)=sqrt((four_over_3*gtmax(i) + rkmax(i))/rho(i))
317
318 viscmax(i) = zero
319 ENDDO
320
321 RETURN
subroutine valpvecdp_v(sig, val, vec, nel)
subroutine valpvec_v(sig, val, vec, nel)