51
52
53
54 use glob_therm_mod
55
56
57
58#include "implicit_f.inc"
59#include "comlock.inc"
60
61
62
63#include "mvsiz_p.inc"
64
65
66
67#include "com08_c.inc"
68#include "param_c.inc"
69#include "scr03_c.inc"
70#include "scr17_c.inc"
71#include "units_c.inc"
72
73
74
75 INTEGER, INTENT(IN) :: ISMSTR
76 INTEGER, INTENT(IN) :: JSMS
77 INTEGER, INTENT(IN) :: ITY
78 INTEGER, INTENT(IN) :: JTUR
79 INTEGER, INTENT(IN) :: JTHE
80 INTEGER, INTENT(IN) :: NFT
81 INTEGER, INTENT(IN) :: JSPH,NPG
82 INTEGER NELTST,ITYPTST ,PID(*), G_DT
83 INTEGER MAT(*),NGL(*),IPLA,NEL,IPM(NPROPMI,*)
85 . dt2t
87 . pdam
89 . pm(npropm,*), off(*), sig(nel,6), eint(*), rho(*), qold(*),
90 . epxe(*), epd(*), vol(*),stifn(*), offg(*),geo(npropg,*) ,
91 . mumax(*), sigy(*), defp(*), dpla(mvsiz),
92 . amu(*), vol_avg(*)
94 . vnew(*), vd2(*), deltax(*), ssp(*), aire(*), vis(*),
95 . psh(*), pnew(1), qnew(*) ,ssp_eq(*),
96 . d1(*), d2(*), d3(*), d4(*), d5(*), d6(*),
97 . sold1(mvsiz), sold2(mvsiz), sold3(mvsiz),
98 . sold4(mvsiz), sold5(mvsiz), sold6(mvsiz),
99 . mssa(*), dmels(*), conde(*),
dtel(*),
100 . rhoref(*) ,rhosp(*)
101 type (glob_therm_) ,intent(inout) :: glob_therm
102
103
104
105 INTEGER ICC(MVSIZ), I, MX, II,,ICC_1
107 . rho0(mvsiz),
108 . g(mvsiz),ak(mvsiz),
109 . qh(mvsiz), c1(mvsiz),
110 . p(mvsiz), epmx(mvsiz), thetl(mvsiz),
111 . ca(mvsiz), cb(mvsiz), cc(mvsiz), cn(mvsiz),
112 . epdr(mvsiz),
113 . dvol(mvsiz), sigmx(mvsiz),
114 . ce(mvsiz),epsl(mvsiz), ql(mvsiz), yldl(mvsiz), aj2(mvsiz),
115 . g0(mvsiz),
116 . e1, e2, e3, e4, e5, e6, g1, g2,
117 . epsp, hl, depsl, alpe, dav, scale, bid1, bid2,
118 . bid3, einc, dta,facq0,
119 . rho0_1,c1_1,ca_1,cb_1,cn_1,
120 . epmx_1,sigmx_1,cc_1,epdr_1,epsl_1,
121 . yldl_1,ql_1
122
123 facq0 = one
124 IF(ipla==0)THEN
125 mx = mat(1)
126 rho0_1 =pm( 1,mx)
127 c1_1 =pm(32,mx)
128 ca_1 =pm(38,mx)
129 cb_1 =pm(39,mx)
130 cn_1 =pm(40,mx)
131 epmx_1 =pm(41,mx)
132 sigmx_1=pm(42,mx)
133 cc_1 =pm(43,mx)
134 epdr_1 =pm(44,mx)
135 epsl_1 =pm(45,mx)
136 yldl_1 =pm(47,mx)
137 ql_1 =pm(48,mx)
138 icc_1 =nint(pm(49,mx))
139 DO 10 i=1,nel
140 rho0(i) =rho0_1
141 g0(i) =pm(22,mx)*off(i)
142 c1(i) =c1_1
143 ca(i) =ca_1
144 cb(i) =cb_1
145 cn(i) =cn_1
146 epmx(i) =epmx_1
147 sigmx(i)=sigmx_1
148 cc(i) =cc_1
149 epdr(i) =epdr_1
150 epsl(i) =epsl_1
151 yldl(i) =yldl_1
152 ql(i) =ql_1
153 icc(i) =icc_1
154 10 CONTINUE
155 ELSE
156
157 mx = mat(1)
158 rho0_1 =pm( 1,mx)
159 c1_1 =pm(32,mx)
160 ca_1 =pm(38,mx)
161 cb_1 =pm(39,mx)
162 cn_1 =pm(40,mx)
163 epmx_1 =pm(41,mx)
164 sigmx_1=pm(42,mx)
165 cc_1 =pm(43,mx)
166 epdr_1 =pm(44,mx)
167 epsl_1 =pm(45,mx)
168 ql_1 =pm(46,mx)
169 yldl_1 =pm(47,mx)
170 icc_1 =nint(pm(49,mx))
171 DO 11 i=1,nel
172 rho0(i) =rho0_1
173 g0(i) =pm(22,mx)*off(i)
174 c1(i) =c1_1
175 ca(i) =ca_1
176 cb(i) =cb_1
177 cn(i) =cn_1
178 epmx(i) =epmx_1
179 sigmx(i)=sigmx_1
180 cc(i) =cc_1
181 epdr(i) =epdr_1
182 epsl(i) =epsl_1
183 ql(i) =ql_1
184 yldl(i) =yldl_1
185 icc(i) =icc_1
186 11 CONTINUE
187 ENDIF
188
189 DO 15 i=1,nel
190 epd(i)=off(i) *
max( abs(d1(i)), abs(d2(i)), abs(d3(i)), half*abs(d4(i)),half*abs(d5(i)),half*abs(d6(i)))
191 epsp =
max(epd(i),epdr(i))
192 ce(i) = one +cc(i) * log(epsp/epdr(i))
193 15 CONTINUE
194
195 IF(ipla/=2)THEN
196 DO 20 i=1,nel
197 IF(cn(i)==one) THEN
198 ak(i)= ca(i)+cb(i)*epxe(i)
199 qh(i)= cb(i)*ce(i)
200 ELSE
201 IF(epxe(i)>zero) THEN
202 ak(i)=ca(i)+cb(i)*epxe(i)**cn(i)
203 IF(cn(i)>one) THEN
204 qh(i)= (cb(i)*cn(i)*epxe(i)**(cn(i)-one))*ce(i)
205 ELSE
206 qh(i)= (cb(i)*cn(i)/epxe(i)**(one - cn(i)))*ce(i)
207 ENDIF
208 ELSE
209 ak(i)=ca(i)
210 qh(i)=zero
211 ENDIF
212 ENDIF
213 sigy(i) = ak(i)
214 IF(epxe(i)>epmx(i))ak(i)=zero
215 IF(epxe(i)>epsl(i))qh(i)=ql(i)
216 20 CONTINUE
217 ELSE
218 DO i=1,nel
219 IF(cn(i)==one) THEN
220 ak(i)= ca(i)+cb(i)*epxe(i)
221 ELSE
222 IF(epxe(i)>zero) THEN
223 ak(i)=ca(i)+cb(i)*epxe(i)**cn(i)
224 ELSE
225 ak(i)=ca(i)
226 ENDIF
227 ENDIF
228 sigy(i) = ak(i)
229 IF(epxe(i)>epmx(i))ak(i)=zero
230 ENDDO
231 ENDIF
232
233 IF(ipla==0)THEN
234 DO i=1,nel
235 ak(i) =
min(ak(i),sigmx(i))
236 hl = three*g0(i)*ql(i)/(three*g0(i)+ql(i))
237 depsl =
max(zero,epxe(i)-epsl(i))
238 ak(i) =
min(ak(i),yldl(i)+hl*depsl)
239 ak(i) =
max(ak(i),zero)
240 alpe =
min(one,ak(i)/
max(ak(i)+three*g0(i)*depsl,em15))
241 ak(i) = ak(i)*ce(i)
242 IF(icc(i)==2) ak(i) =
min(ak(i),sigmx(i))
243 g(i) = alpe*g0(i)
244 c1(i) = pdam*alpe*c1(i) + (1.-pdam)*c1(i)
245 p(i) =-third*(sig(i,1)+sig(i,2)+sig(i,3))
246 g1=dt1*g(i)
247 g2=two*g1
248 ssp(i)=sqrt((onep333*g(i)+c1(i))/rho0(i))
249
250
251
252 dav =-third*(d1(i)+d2(i)+d3(i))
253 sig(i,1)=sig(i,1)+p(i)+g2*(d1(i)+dav)
254 sig(i,2)=sig(i,2)+p(i)+g2*(d2(i)+dav)
255 sig(i,3)=sig(i,3)+p(i)+g2*(d3(i)+dav)
256 sig(i,4)=sig(i,4)+g1*d4(i)
257 sig(i,5)=sig(i,5)+g1*d5(i)
258 sig(i,6)=sig(i,6)+g1*d6(i)
259 aj2(i)=half*(sig(i,1)**2+sig(i,2)**2+sig(i,3)**2)
260 1 +sig(i,4)**2+sig(i,5)**2+sig(i,6)**2
261 aj2(i)=sqrt(three*aj2(i))
262 ENDDO
263 ELSEIF(ipla==2)THEN
264 DO i=1,nel
265 ak(i) =
min(ak(i),sigmx(i))
266 depsl =
max(zero,epxe(i)-epsl(i))
267 ak(i) =
min(ak(i),yldl(i)+ql(i)*depsl)
268 ak(i) =
max(ak(i),zero)
269 alpe =
min(one,ak(i)/
max(ak(i) + three*g0(i)*depsl,em15))
270 ak(i) = ak(i)*ce(i)
271 IF(icc(i)==2) ak(i) =
min(ak(i),sigmx(i))
272 g(i) = alpe*g0(i)
273 c1(i) = pdam*alpe*c1(i) + (one - pdam)*c1(i)
274 p(i) =-third*(sig(i,1)+sig(i,2)+sig(i,3))
275 g1=dt1*g(i)
276 g2=two*g1
277 ssp(i)=sqrt((onep333*g(i)+c1(i))/rho0(i))
278
279
280
281 IF(epxe(i)>epsl(i).AND.three*g(i)<em15)THEN
282 sig(i,1)=zero
283 sig(i,2)=zero
284 sig(i,3)=zero
285 sig(i,4)=zero
286 sig(i,5)=zero
287 sig(i,6)=zero
288 aj2(i) =zero
289 ELSE
290 dav =-third*(d1(i)+d2(i)+d3(i))
291 sig(i,1)=sig(i,1)+p(i)+g2*(d1(i)+dav)
292 sig(i,2)=sig(i,2)+p(i)+g2*(d2(i)+dav)
293 sig(i,3)=sig(i,3)+p(i)+g2*(d3(i)+dav)
294 sig(i,4)=sig(i,4)+g1*d4(i)
295 sig(i,5)=sig(i,5)+g1*d5(i)
296 sig(i,6)=sig(i,6)+g1*d6(i)
297 aj2(i)=half*(sig(i,1)**2+sig(i,2)**2+sig(i,3)**2)+
298 . sig(i,4)**2+sig(i,5)**2+sig(i,6)**2
299 aj2(i)=sqrt(three*aj2(i))
300 ENDIF
301 ENDDO
302 ELSE
303 DO i=1,nel
304 ak(i) =
min(ak(i),sigmx(i))
305 depsl =
max(zero,epxe(i)-epsl(i))
306 ak(i) =
min(ak(i),yldl(i)+ql(i)*depsl)
307 ak(i) =
max(ak(i),zero)
308 alpe =
min(one,ak(i)/
max(ak(i) + three*g0(i)*depsl,em15))
309 ak(i) = ak(i)*ce(i)
310 IF(icc(i)==2) ak(i) =
min(ak(i),sigmx(i))
311 g(i) = alpe*g0(i)
312 c1(i) = pdam*alpe*c1(i) + (one - pdam)*c1(i)
313 p(i) =-third*(sig(i,1)+sig(i,2)+sig(i,3))
314 g1=dt1*g(i)
315 g2=two*g1
316 ssp(i)=sqrt((onep333*g(i)+c1(i))/rho0(i))
317
318
319
320 IF(epxe(i)>epsl(i).AND.g(i)*(three*g0(i)+qh(i))<em15)THEN
321 sig(i,1)=zero
322 sig(i,2)=zero
323 sig(i,3)=zero
324 sig(i,4)=zero
325 sig(i,5)=zero
326 sig(i,6)=zero
327 aj2(i) =zero
328 ELSE
329 dav =-third*(d1(i)+d2(i)+d3(i))
330 sig(i,1)=sig(i,1)+p(i)+g2*(d1(i)+dav)
331 sig(i,2)=sig(i,2)+p(i)+g2*(d2(i)+dav)
332 sig(i,3)=sig(i,3)+p(i)+g2*(d3(i)+dav)
333 sig(i,4)=sig(i,4)+g1*d4(i)
334 sig(i,5)=sig(i,5)+g1*d5(i)
335 sig(i,6)=sig(i,6)+g1*d6(i)
336 aj2(i)=half*(sig(i,1)**2+sig(i,2)**2+sig(i,3)**2)+
337 . sig(i,4)**2+sig(i,5)**2+sig(i,6)**2
338 aj2(i)=sqrt(three*aj2(i))
339 ENDIF
340 ENDDO
341 ENDIF
342
343 IF(ipla==0)THEN
344 DO i=1,nel
345 scale=
min(one,ak(i) /
max(aj2(i),em15))
346 sig(i,1)=scale*sig(i,1)
347 sig(i,2)=scale*sig(i,2)
348 sig(i,3)=scale*sig(i,3)
349 sig(i,4)=scale*sig(i,4)
350 sig(i,5)=scale*sig(i,5)
351 sig(i,6)=scale*sig(i,6)
352 epxe(i)=epxe(i)+(one -scale)*aj2(i)/
max(three*g(i)+qh(i),em15)
353 dpla(i) = (one -scale)*aj2(i)/
max(three*g(i)+qh(i),em15)
354 ENDDO
355 ELSEIF(ipla==2)THEN
356 DO i=1,nel
357 scale=
min(one,ak(i) /
max(aj2(i),em15))
358 sig(i,1)=scale*sig(i,1)
359 sig(i,2)=scale*sig(i,2)
360 sig(i,3)=scale*sig(i,3)
361 sig(i,4)=scale*sig(i,4)
362 sig(i,5)=scale*sig(i,5)
363 sig(i,6)=scale*sig(i,6)
364 epxe(i)=epxe(i)+(one -scale)*aj2(i)/
max(three*g(i),em15)
365 dpla(i) = (one -scale)*aj2(i)/
max(three*g(i),em15)
366 ENDDO
367 ELSEIF(ipla==1)THEN
368 DO i=1,nel
369 scale=
min(one,ak(i) /
max(aj2(i),em15))
370 dpla(i)=(one -scale)*aj2(i)*g0(i)/
max(g(i)*(three*g0(i)+qh(i)),em15)
371 ak(i)=
max(zero,ak(i)+dpla(i)*qh(i))
372 scale=
min(one,ak(i)/
max(aj2(i),em15))
373 sig(i,1)=scale*sig(i,1)
374 sig(i,2)=scale*sig(i,2)
375 sig(i,3)=scale*sig(i,3)
376 sig(i,4)=scale*sig(i,4)
377 sig(i,5)=scale*sig(i,5)
378 sig(i,6)=scale*sig(i,6)
379 epxe(i)=epxe(i)+dpla(i)
380 ENDDO
381 ENDIF
382
383 IF (jsph==0)THEN
385 1 pm, off, rho, bid1,
386 2 bid2, ssp, bid3, stifn,
387 3 dt2t, neltst, ityptst, aire,
388 4 offg, geo, pid, vnew,
389 5 vd2, deltax, vis, d1,
390 6 d2, d3, pnew, psh,
391 7 mat, ngl, qnew, ssp_eq,
392 8 vol, mssa, dmels, ibid,
393 9 facq0, conde,
dtel, g_dt,
394 a ipm, rhoref, rhosp, nel,
395 b ity, ismstr, jtur, jthe,
396 c jsms, npg , glob_therm)
397 ELSE
399 1 pm, off, rho, bid1,
400 2 bid2, bid3, stifn, dt2t,
401 3 neltst, ityptst, offg, geo,
402 4 pid, mumax, ssp, vnew,
403 5 vd2, deltax, vis, d1,
404 6 d2, d3, pnew, psh,
405 7 mat, ngl, qnew, ssp_eq,
406 8 g_dt,
dtel, nel, ity,
407 9 jtur, jthe)
408 ENDIF
409
410 DO 120 i=1,nel
411 pnew(i)=c1(i)*amu(i)
412 120 CONTINUE
413
414
415
416
417 IF(pdam==one.AND.codvers>=42)THEN
418 DO i=1,nel
419 IF(off(i)<em01) off(i)=zero
420 IF(off(i)<one) off(i)=off(i)*four_over_5
421 ENDDO
422
423 DO i=1,nel
424 IF(epmx(i)/=zero.AND.off(i)>=one.AND. epxe(i)>epmx(i))THEN
425 off(i)=off(i)*four_over_5
426 ii=i+nft
427#include "lockon.inc"
428 WRITE(iout,1000) ngl(i)
429#include "lockoff.inc"
430 ENDIF
431 ENDDO
432 ENDIF
433
434 dta = half*dt1
435
436 DO i=1,nel
437 sig(i,1)=(sig(i,1)-pnew(i))*off(i)
438 sig(i,2)=(sig(i,2)-pnew(i))*off(i)
439 sig(i,3)=(sig(i,3)-pnew(i))*off(i)
440 sig(i,4)= sig(i,4)*off(i)
441 sig(i,5)= sig(i,5)*off(i)
442 sig(i,6)= sig(i,6)*off(i)
443 e1=d1(i)*(sold1(i)+sig(i,1))
444 e2=d2(i)*(sold2(i)+sig(i,2))
445 e3=d3(i)*(sold3(i)+sig(i,3))
446 e4=d4(i)*(sold4(i)+sig(i,4))
447 e5=d5(i)*(sold5(i)+sig(i,5))
448 e6=d6(i)*(sold6(i)+sig(i,6))
449 einc= vol_avg(i)*(e1+e2+e3+e4+e5+e6)*dta - half*dvol(i)*(qold(i)+qnew(i))
450 eint(i)=(eint(i)+einc*off(i)) /
max(em15,vol(i))
451 qold(i)=qnew(i)
452 ENDDO
453
454 DO i=1,nel
455 defp(i)=epxe(i)
456 sigy(i)=
max(sigy(i),ak(i))
457 ENDDO
458
459 1000 FORMAT(1x,'RUPTURE OF SOLID ELEMENT NUMBER ',i10)
460 RETURN
subroutine dtel(ssp, pm, geo, pid, mat, rho0, vis, deltax, aire, vol, dtx)
subroutine mdtsph(pm, off, rho, rk, t, re, sti, dt2t, neltst, ityptst, offg, geo, pid, mumax, ssp, vol, vd2, deltax, vis, d1, d2, d3, pnew, psh, mat, ngl, qvis, ssp_eq, g_dt, dtsph, nel, ity, jtur, jthe)
subroutine mqviscb(pm, off, rho, rk, temp, ssp, re, sti, dt2t, neltst, ityptst, aire, offg, geo, pid, vol, vd2, deltax, vis, d1, d2, d3, pnew, psh, mat, ngl, qvis, ssp_eq, vol0, mssa, dmels, igeo, facq0, conde, dtel, g_dt, ipm, rhoref, rhosp, nel, ity, ismstr, jtur, jthe, jsms, npg, glob_therm)
1.15.0