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, ,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
81 . soundsp(nel), viscmax(nel), et(nel)
82
83
84
86 . uvar(nel,nuvar), off(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,II
97
99 . g,rbulk,aa,bb,cc,p,trace,c(5),
100 . t1(mvsiz), t2(mvsiz),t3(mvsiz),av(mvsiz,6),ev(mvsiz,3),
101 . evv(mvsiz,3),rv(mvsiz),rvd,l1di1lam1,l2di1lam2,l3di1lam3,
102 . dirprv(mvsiz,3,3),evd(3),c0(3),d(3),invr,eti,et1,et2,et3,
103 . ev2(3),clam(3),lam_2(3),lam_4(3),amax,clp
104 my_real ,
DIMENSION(NEL,3) :: evm,cii
105 my_real ,
DIMENSION(NEL) :: gtmax,rkmax,trac3
106
107
108 g = uparam(1)
109 rbulk = uparam(2)
110 c0(1) = uparam(4)
111 c0(2) = uparam(5)
112 c0(3) = uparam(6)
113
114 d(1) = uparam(7)
115 d(2) = uparam(8)
116 d(3) = uparam(9)
117
118
119
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) 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
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 DO i=1,nel
178
179 ev2(1) = ev(i,1)*ev(i,1)
180 ev2(2) = ev(i,2)*ev(i,2)
181 ev2(3) = ev(i,3)*ev(i,3)
182
183 trace = ev2(1) + ev2(2) + ev2(3)
184 aa = (trace - three)
185 bb = aa*aa
186 cc = (c0(1) + two*c0(2)*aa + three*c0(3)*bb)
187 et1 = ev2(1)*cc + ev2(1)**2*(two
188 et2 = ev2(2)*cc + ev2(2)**2*(two*c0(2) + six*c0(3)*aa)
189 et3 = ev2(3)*cc + ev2(3)**2*(two*c0(2) + six*c0(3)*aa)
190 et(i)=
max(et1,et2,et3)
191 et(i) = four*et(i)
192 et(i)=
max(one,et(i)/
max(em20,g))
193 ENDDO
194 ENDIF
195
196
197
198 DO i = 1,nel
199
200 rv(i) = ev(i,1)*ev(i,2)*ev(i,3)
201 ENDDO
202
203 IF(iexpan > 0.AND.(ismstr==10.OR.ismstr==11.OR.ismstr==12)) THEN
204 DO i=1,nel
205 rv(i) = rv(i) -epsth3(i)
206 ENDDO
207 ENDIF
208
209 DO i = 1,nel
210 invr = one/(
max(rv(i), em20))
211
212
213
214
215 rvd = zero
216 IF(rv(i) > zero) rvd = exp((-third)*log(rv(i)))
217
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 = evd(1)*evd(1) + evd(2)*evd(2) + evd(3)*evd(3)
224 l1di1lam1 = two*(evd(1)**2 - third*trace)
225 l2di1lam2 = two*(evd(2)**2 - third*trace)
226 l3di1lam3 = two*(evd(3)**2 - third*trace)
227
228 aa = (trace - three)
229 bb = aa*aa
230 cc = (c0(1) + two*c0(2)*aa + three*c0(3)*bb)*invr
231 trac3(i) = aa
232
233
234 t1(i) = l1di1lam1*cc
235 t2(i) = l2di1lam2*cc
236 t3(i) = l3di1lam3*cc
237
238 p = rbulk*(rv(i) - one) + four*d(2)*(rv(i) - one)**3 + six*d(3)*(rv(i) - one)**5
239
240 t1(i) = t1(i) + p
241 t2(i) = t2(i) + p
242 t3(i) = t3(i) + p
243 ENDDO
244
245
246 gtmax(1:nel
247 rkmax(1:nel) = rbulk
248 cii(1:nel,1:3) = zero
249 DO ii = 1,3
250 clp = four*ii*c0(ii)
251 DO i=1,nel
252 IF (abs(trac3(i))<em10) cycle
253 lam_2(1:3) = evm(i,1:3)*evm(i,1:3)
254 lam_4(1:3) = lam_2(1:3)*lam_2(1:3)
255 aa = one_over_9*ii*trac3(i)**ii
256 bb = zero
257 cc = zero
258 IF (ii>1) bb =third*(3-ii)*trac3(i)**(ii-1)
259 IF (ii>2) cc =(ii-1)*trac3(i)**(ii-
260 cii(i,1:3) = cii(i,1:3) +clp*(aa+bb*lam_2(1:3) +cc*lam_4(1:3))
261 ENDDO
262 ENDDO
263 DO i=1,nel
264 amax=
max(cii(i,1),cii(i,2),cii
265
266 eti =
max(one,amax*0.81)
267 gtmax(i) = g*eti
268
269 rkmax(i) = rbulk+twelve*d(2)*(rv(i) - one)**2 + 30*d(3)*(rv(i) - one)**4
270 rkmax(i) =
max(rbulk,rkmax(i))
272 ENDDO
273
274 DO i=1,nel
275 signxx(i) = dirprv(i,1,1)*dirprv(i,1,1)*t1(i)
276 . + dirprv(i,1,2)*dirprv(i,1,2)*t2(i)
277 .
278
279 signyy(i) = dirprv(i,2,2)*dirprv(i,2,2)*t2(i)
280 . + dirprv(i,2,3)*dirprv(i,2,3)*t3(i)
281 . + dirprv(i,2,1)*dirprv(i,2,1)*t1(i)
282
283 signzz(i) = dirprv(i,3,3)*dirprv(i,3,3)*t3(i)
284 . + dirprv(i,3,1)*dirprv(i,3,1)*t1(i)
285 . + dirprv(i,3,2)*dirprv(i,3,2)*t2(i)
286
287 signxy(i) = dirprv(i,1,1)*dirprv(i,2,1)*t1(i)
288 . + dirprv(i,1,2)*dirprv(i,2,2)*t2(i)
289 . + dirprv(i,1,3)*dirprv(i,2,3)*t3(i)
290
291 signyz(i) = dirprv(i,2,2)*dirprv(i,3,2)*t2(i)
292 . + dirprv(i,2,3)*dirprv(i,3,3)*t3(i)
293 . + dirprv(i,2,1)*dirprv(i,3,1)*t1(i)
294
295 signzx(i) = dirprv(i,3,3)*dirprv(i,1,3)*t3(i)
296 . + dirprv(i,3,1)*dirprv(i,1,1)*t1(i)
297 . + dirprv(i,3,2)*dirprv(i,1,2)*t2
298
299
300 soundsp(i)=sqrt((four_over_3*gtmax(i) + rkmax(i))/rho(i))
301
302 viscmax(i) = zero
303 ENDDO
304
305 RETURN
subroutine valpvecdp_v(sig, val, vec, nel)
subroutine valpvec_v(sig, val, vec, nel)