56
57
58
59 use glob_therm_mod
60
61
62
63#include "implicit_f.inc"
64#include "comlock.inc"
65
66
67
68#include "mvsiz_p.inc"
69
70
71
72#include "com08_c.inc"
73#include "param_c.inc"
74#include "impl1_c.inc"
75#include "units_c.inc"
76
77
78
79 INTEGER, INTENT(IN) :: NPG
80 INTEGER, INTENT(IN) :: ISMSTR
81 INTEGER, INTENT(IN) :: JSMS
82 INTEGER, INTENT(IN) :: ITY
83 INTEGER, INTENT(IN) :: JTUR
84 INTEGER, INTENT(IN) :: JTHE
85 INTEGER, INTENT(IN) :: JSPH
86 INTEGER, INTENT(IN) :: IEOS
87 INTEGER, INTENT(IN) :: JLAG
88
89 INTEGER NELTST,ITYPTST,PID(NEL),G_DT,NEL,IPG
90 INTEGER MAT(NEL),NGL(NEL), IPM(NPROPMI,*)
91 my_real,
DIMENSION(NEL) ,
INTENT(INOUT) :: plap
93 my_real pm(npropm,*), off(nel), sig(nel,6), eint(nel), rho(nel), qold(nel),
94 . epxe(nel), epsd(nel), vol(nel), stifn(nel), offg(nel),geo(npropg,*),
95 . mumax(nel),amu(nel),vol_avg(nel)
96 my_real vnew(nel), vd2(nel), deltax(nel), ssp(nel), aire(nel), vis(nel),
97 . psh(nel), pnew(nel),qnew(nel) ,ssp_eq(nel), dvol(nel),
98 . d1(nel), d2(nel), d3(nel), d4(nel), d5(nel), d6(nel),
99 . sold1(mvsiz), sold2(mvsiz), sold3(mvsiz),
100 . sold4(mvsiz), sold5(mvsiz), sold6(mvsiz),
101 . tstar(mvsiz), dpla(mvsiz), epsp(mvsiz), dmg(nel)
102 my_real,
DIMENSION(MVSIZ) ,
INTENT(IN) :: dpdm
103 my_real,
DIMENSION(NEL) ,
INTENT(INOUT) :: fheat
104 my_real sigy(nel) ,defp(nel),etse(nel), mssa(nel), dmels(nel), tempel(nel),
105 . sigbak(nel,6),al_imp(nel),signor(mvsiz,6),conde(nel),
dtel(nel),
106 . rhoref(nel) ,rhosp(nel)
107 type (glob_therm_) ,intent(inout) :: glob_therm
108
109
110
111 INTEGER :: I, , ICC,IRTY,ISRATE,VP,IDEV,MX,IBID,NINDX
112 INTEGER :: INDX()
113 my_real rho0, plap1,pold,pshift,bulk,p, epmx,cc,cn,epdr, cmx,
114 . z3,z4,rhocp,rhocpi,t_ref,t_melt,asrate,
115 . e1, e2, e3, e4, e5,e6, einc, g0,g1,g2,g3,g3h,
116 . scale, bid1, bid2, bid3, dta,mt,
117 . dsxx,dsyy,dszz,dsxy,dsyz,dszx,
alpha,hkin,
118 . facq0,fisokin,beta,vm,vm_1,norm_1,ca0,cb,sigm0
119 my_real epd(mvsiz),g(mvsiz), ak(mvsiz),pc(mvsiz), qh(mvsiz),
120 . aj2(mvsiz), dav(mvsiz),ca(mvsiz),sigmx(mvsiz),
121 . sigexx(mvsiz),sigeyy(mvsiz),sigezz(mvsiz),
122 . sigexy(mvsiz),sigeyz(mvsiz),sigezx(mvsiz)
123 my_real,
DIMENSION(MVSIZ) :: bidmvsiz
124
125
126
127 facq0 = one
128
129 mx = mat(1)
130 rho0 = pm( 1,mx)
131 bulk = pm(32,mx)
132 cn = pm(40,mx)
133 cc = pm(43,mx)
134 icc = nint(pm(49,mx))
135 epmx = pm(41,mx)
136 epdr = pm(44,mx)
137 irty = nint(pm(50,mx))
138 z3 = pm(51,mx)
139 cmx = pm(45,mx)
140 israte = ipm(3,mx)
141 asrate =
min(one,pm(9,mx)*dt1)
142 fisokin = pm(55,mx)
143 g0 = pm(22,mx)
144 ca0 = pm(38,mx)
145 cb = pm(39,mx)
146 sigm0 = pm(42,mx)
147 vp = ipm(255,mx)
148 pshift = pm(88,mx)
149
150 IF (irty == 1) THEN
151 rhocpi = pm(53,mx)
152 t_ref = pm(54,mx)
153 z4 = pm(52,mx)
154 ELSE
155 rhocp = pm(69,mx)
156 IF (rhocp > zero) rhocpi = one / rhocp
157 t_ref = pm(79,mx)
158 t_melt = pm(80,mx)
159 END IF
160
161 DO i=1,nel
162 g(i) = g0*off(i)
163 ca(i) = ca0
164 sigmx(i)= sigm0
165 etse(i) = one
166 ENDDO
167
168
169
170
171 IF (fisokin > zero) THEN
172 DO i=1,nel
173 sig(i,1)=sig(i,1)-sigbak(i,1)
174 sig(i,2)=sig(i,2)-sigbak(i,2)
175 sig(i,3)=sig(i,3)-sigbak(i,3)
176 sig(i,4)=sig(i,4)-sigbak(i,4)
177 sig(i,5)=sig(i,5)-sigbak(i,5)
178 sig(i,6)=sig(i,6)-sigbak(i,6)
179 ENDDO
180 ENDIF
181
182 DO i=1,nel
183 p =-third*(sig(i,1)+sig(i,2)+sig(i,3))
184 dav(i)=-third*(d1(i)+d2(i)+d3(i))
185 g1=dt1*g(i)
186 g2=two*g1
187 ssp(i)=sqrt((onep333*g(i)+bulk)/rho0)
188
189 sig(i,1)=sig(i,1)+p+g2*(d1(i)+dav(i))
190 sig(i,2)=sig(i,2)+p+g2*(d2(i)+dav(i))
191 sig(i,3)=sig(i,3)+p+g2*(d3(i)+dav(i))
192 sig(i,4)=sig(i,4)+g1*d4(i)
193 sig(i,5)=sig(i,5)+g1*d5(i)
194 sig(i,6)=sig(i,6)+g1*d6(i)
195
196
197 aj2(i)=half*(sig(i,1)**2+sig(i,2)**2+sig(i,3)**2)
198 . + sig(i,4)**2+sig(i,5)**2+sig(i,6)**2
199 aj2(i)=sqrt(three*aj2(i))
200 ENDDO
201
202
203 IF (impl_s>0 .OR. fisokin>zero) THEN
204 DO i=1,nel
205 sigexx(i) = sig(i,1)
206 sigeyy(i) = sig(i,2)
207 sigezz(i) = sig(i,3)
208 sigexy(i) = sig(i,4)
209 sigeyz(i) = sig(i,5)
210 sigezx(i) = sig(i,6)
211 ENDDO
212 ENDIF
213
214
215
216 idev = vp - 2
218 . d1, d2, d3, d4, d5, d6)
219
220 epsp(1:nel) = epsd(1:nel)
221 epd(1:nel) = one
222
223 IF (cc /= zero) THEN
224 IF(vp == 1)THEN
225 DO i=1,nel
226 epd(i)=
max(plap(i),epdr)
227 epd(i)= log(epd(i)/epdr)
228 ENDDO
229 ELSE
230 DO i=1,nel
231 epd(i)=
max(epsd(i),em15)
232 epd(i)= log(epd(i)/epdr)
233 ENDDO
234 ENDIF
235
236 IF (irty == 0) THEN
237 DO i=1,nel
239 epd(i)=
max(zero,epd(i))
240 epd(i)= (one +cc * epd(i))*(one -tstar(i)**mt)
241 IF(icc==1) sigmx(i)= sigm0*epd(i)
242 ENDDO
243 ELSEIF (irty==1) THEN
244 DO i=1,nel
245 epd(i)= cc*exp((-z3+z4 * epd(i))*tempel(i))
246 IF (icc==1) sigmx(i) = sigm0 + epd(i)
247 ca(i) = ca0 + epd(i)
248 epd(i)=one
249 ENDDO
250 END IF
251 ELSEIF (irty == 0) THEN
253 epd(1:nel) = one - tstar(1:nel)**mt
254 ENDIF
255
256
257
258 IF (fisokin == zero) THEN
259 IF(cn==one) THEN
260 ak(1:nel)= ca(1:nel)+cb*epxe(1:nel)
261 qh(1:nel)= cb*epd(1:nel)
262 ELSE
263 DO i=1,nel
264 IF(epxe(i)>zero) THEN
265 ak(i)=ca(i)+cb*epxe(i)**cn
266 IF(cn>one) THEN
267 qh(i) = (cb*cn*epxe(i)**(cn - one))*epd(i)
268 ELSE
269 qh(i) = (cb*cn/epxe(i)**(one -cn))*epd(i)
270 ENDIF
271 ELSE
272 ak(i)=ca(i)
273 qh(i)=zero
274 ENDIF
275 ENDDO
276 ENDIF
277 DO i=1,nel
278 ak(i)=ak(i)*epd(i)
279 IF(sigmx(i)<ak(i))THEN
280 ak(i)=sigmx(i)
281 qh(i)=zero
282 ENDIF
283 sigy(i) = ak(i)
284 IF(epxe(i)>epmx)THEN
285 ak(i)=zero
286 qh(i)=zero
287 ENDIF
288 ENDDO
289
290 ELSE
291
292 DO i=1,nel
293 beta = one-fisokin
294
295 IF(cn==one) THEN
296 sigy(i) = ca(i)+cb*epxe(i)
297 ak(i)= ca(i)+beta*cb*epxe(i)
298 qh(i)= cb*epd(i)
299 ELSE
300 IF(epxe(i)>zero) THEN
301 sigy(i)=ca(i)+cb*epxe(i)**cn
302 ak(i)=ca(i)+beta*cb*epxe(i)**cn
303 IF(cn>one) THEN
304 qh(i)= (cb*cn*epxe(i)**(cn - one))*epd(i)
305 ELSE
306 qh(i)= (cb*cn/epxe(i)**(one -cn))*epd(i)
307 ENDIF
308 ELSE
309 ak(i)=ca(i)
310 sigy(i)=ca(i)
311 qh(i)=zero
312 ENDIF
313 ENDIF
314 ak(i)=ak(i)*epd(i)
315 sigy(i)=sigy(i)*epd(i)
316 IF(sigmx(i)<ak(i))THEN
317 ak(i)=sigmx(i)
318 qh(i)=zero
319 ENDIF
320 sigy(i)=
min(sigy(i),sigmx(i))
321 IF(epxe(i)>epmx)THEN
322 ak(i)=zero
323 qh(i)=zero
324 ENDIF
325 END DO
326 END IF
327
328 IF (impl_s==0) THEN
329 DO i=1,nel
330 scale=
min(one,ak(i)/
max(aj2(i),em15))
331
332 dpla(i)=(one-scale)*aj2(i)/
max(three*g(i)+qh(i),em15)
333
334 ak(i)=ak(i)+(one - fisokin)*dpla(i)*qh(i)
335 scale=
min(one,ak(i)/
max(aj2(i),em15))
336 sig(i,1)=scale*sig(i,1)
337 sig(i,2)=scale*sig(i,2)
338 sig(i,3)=scale*sig(i,3)
339 sig(i,4)=scale*sig(i,4)
340 sig(i,5)=scale*sig(i,5)
341 sig(i,6)=scale*sig(i,6)
342 epxe(i)=epxe(i)+dpla(i)
343 ENDDO
344
345 ELSE
347 1 ca, cb, cn, epd,
348 2 sigmx, epmx, dpla, ak,
349 3 qh, sigy, fisokin, nel)
350
351 END IF
352
353 IF (fisokin > zero) THEN
354 DO i=1,nel
355 dsxx = sigexx(i) - sig(i,1)
356 dsyy = sigeyy(i) - sig(i,2)
357 dszz = sigezz(i) - sig(i,3)
358 dsxy = sigexy(i) - sig(i,4)
359 dsyz = sigeyz(i) - sig(i,5)
360 dszx = sigezx(i) - sig(i,6)
361
362 hkin = two_third*fisokin*qh(i)
363 alpha = hkin/
max(two*g(i)+hkin,em15)
364
365 sigbak(i,1) = sigbak(i,1) +
alpha*dsxx
366 sigbak(i,2) = sigbak(i,2) +
alpha*dsyy
367 sigbak(i,3) = sigbak(i,3) +
alpha*dszz
368 sigbak(i,4) = sigbak(i,4) +
alpha*dsxy
369 sigbak(i,5) = sigbak(i,5) +
alpha*dsyz
370 sigbak(i,6) = sigbak(i,6) +
alpha*dszx
371
372 sig(i,1)=sig(i,1) + sigbak(i,1)
373 sig(i,2)=sig(i,2) + sigbak(i,2)
374 sig(i,3)=sig(i,3) + sigbak(i,3)
375 sig(i,4)=sig(i,4) + sigbak(i,4)
376 sig(i,5)=sig(i,5) + sigbak(i,5)
377 sig(i,6)=sig(i,6) + sigbak(i,6)
378 ENDDO
379 END IF
380
381
382 bidmvsiz(1:mvsiz) = zero
383 IF (jsph==0) THEN
385 1 pm, off, rho, bidmvsiz,
386 2 bidmvsiz,ssp, bidmvsiz,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, bidmvsiz,
400 2 bidmvsiz,bidmvsiz,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 dta = half*dt1
411
412 nindx = 0
413 indx(1:nel) = 0
414 DO i=1,nel
415 IF ((epxe(i) > epmx).AND.(dmg(i)==zero)) THEN
416 nindx = nindx + 1
417 indx(nindx) = i
418 dmg(i) = one
419 ENDIF
420 ENDDO
421
422 IF (ieos == 0) THEN
423 DO i=1,nel
424 pnew(i) = bulk*amu(i)
425 sig(i,1) =(sig(i,1) - pnew(i))*off(i)
426 sig(i,2) =(sig(i,2) - pnew(i))*off(i)
427 sig(i,3) =(sig(i,3) - pnew(i))*off(i)
428 sig(i,4) = sig(i,4) * off(i)
429 sig(i,5) = sig(i,5) * off(i)
430 sig(i,6) = sig(i,6) * off(i)
431 ENDDO
432
433 DO i=1,nel
434 e1 = d1(i)*(sold1(i)+sig(i,1))
435 e2 = d2(i)*(sold2(i)+sig(i,2))
436 e3 = d3(i)*(sold3(i)+sig(i,3))
437 e4 = d4(i)*(sold4(i)+sig(i,4))
438 e5 = d5(i)*(sold5(i)+sig(i,5))
439 e6 = d6(i)*(sold6(i)+sig(i,6))
440 einc = vol_avg(i)*(e1+e2+e3+e4+e5+e6)*dta - half*dvol(i)*(qold(i)+qnew(i))
441 eint(i) = (eint(i)+einc*off(i)) /
max(em15,vol(i))
442 ENDDO
443
444 ELSE
445
446
447
448 DO i = 1, nel
449 ssp(i) = sqrt((dpdm(i) + four*g(i)/three)/rho0)
450
451 pold = (sold1(i)+sold2(i)+sold3(i)) * third
452 sold1(i) = sold1(i) - pold
453 sold2(i) = sold2(i) - pold
454 sold3(i) = sold3(i) - pold
455 e1 = d1(i) * (sold1(i)+sig(i,1))
456 e2 = d2(i) * (sold2(i)+sig(i,2))
457 e3 = d3(i) * (sold3(i)+sig(i,3))
458 e4 = d4(i) * (sold4(i)+sig(i,4))
459 e5 = d5(i) * (sold5(i)+sig(i,5))
460 e6 = d6(i) * (sold6(i)+sig(i,6))
461 einc = vol_avg(i) * (e1+e2+e3+e4+e5+e6) * dta
462 eint(i) = eint(i) + (einc+half*dvol(i)*(pold-pshift-qold(i)-qnew(i)))*off(i)
463 ENDDO
464 END IF
465
466 DO i=1,nel
467 qold(i) = qnew(i)
468 defp(i) = epxe(i)
469 sigy(i) =
max(sigy(i),ak(i))
470 ENDDO
471 DO i=1,nel
472 IF (dpla(i) > 0) etse(i)= half*qh(i)*off(i)/
max(g(i),em15)
473 ENDDO
474
475
476
477 IF (impl_s>0) THEN
478 IF (ikt==0) RETURN
479 IF (fisokin == zero) THEN
480 DO i=1,nel
481 IF (dpla(i)>0)THEN
482 IF(cn==one) THEN
483 qh(i)= cb*epd(i)
484 ELSEIF(cn>one) THEN
485 qh(i)= (cb*cn*epxe(i)**(cn - one))*epd(i)
486 ELSE
487 qh(i)= (cb*cn/epxe(i)**(one -cn))*epd(i)
488 ENDIF
489 ENDIF
490 ENDDO
491 ENDIF
492
493 DO i = 1,nel
494 IF (dpla(i)>zero) THEN
495
496
497
498 vm =half*(sigexx(i)**2+sigeyy(i)**2+sigezz(i)**2)
499 1 +sigexy(i)**2+sigeyz(i)**2+sigezx(i)**2
500 vm_1 =one/sqrt(three*vm)
501 g3 = three*g(i)
502 g3h =
max(g3+qh(i),em15)
503 scale =
max(zero,one-g3h*dpla(i)*vm_1)
504
505
506 norm_1=g3*vm_1*sqrt(scale/g3h)
507
508 signor(i,1)=sigexx(i)*norm_1
509 signor(i,2)=sigeyy(i)*norm_1
510 signor(i,3)=sigezz(i)*norm_1
511 signor(i,4)=sigexy(i)*norm_1
512 signor(i,5)=sigeyz(i)*norm_1
513 signor(i,6)=sigezx(i)*norm_1
514
515
516 al_imp(i)= one - g3*dpla(i)*vm_1
517 ELSE
518 al_imp(i)=one
519 ENDIF
520 ENDDO
521 ENDIF
522
523
524
525
526
527
528 IF (vp== 1) THEN
529 DO i=1,nel
530 plap1 = dpla(i)/
max(em20,dt1)
531 plap(i) = asrate * plap1 + (one - asrate) * plap(i)
532 ENDDO
533 ENDIF
534
535
536
537 IF (jthe < 0 .AND. jlag /= 0) THEN
538 DO i=1,nel
539 fheat(i) = fheat(i) + sigy(i)*dpla(i)*vol(i)
540 ENDDO
541 ELSEIF(rhocp > zero)THEN
542 DO i=1,nel
543 tempel(i) = tempel(i) + sigy(i)*dpla(i) / rhocp
544
545
546
547
548
549 ENDDO
550 END IF
551
552
553 IF (nindx > 0) THEN
554 DO j=1,nindx
555#include "lockon.inc"
556 WRITE(iout, 1000) ngl(indx(j)),ipg
557 WRITE(istdo,1100) ngl(indx(j)),ipg
558#include "lockoff.inc"
559 ENDDO
560 ENDIF
561
562 1000 FORMAT(1x,'EXCEEDED EPS_MAX ON SOLID ELEMENT NUMBER ',i10,
563 . ': DEVIATORIC STRESS SET TO 0 ON INTEGRATION POINT ',i5 )
564 1100 FORMAT(1x,'EXCEEDED EPS_MAX ON SOLID ELEMENT NUMBER ',i10,
565 . ': DEVIATORIC STRESS SET TO 0 ON INTEGRATION POINT ',i5 ,
566 . ' AT TIME :',g11.4)
567
568 RETURN
subroutine dtel(ssp, pm, geo, pid, mat, rho0, vis, deltax, aire, vol, dtx)
subroutine m2iter_imp(sig, epxe, aj2, g, ca, cb, cn, epd, sigmx, epmx, dpla1, ak, qh, sigy, fisokin, nel)
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)
subroutine mstrain_rate(nel, israte, asrate, epsd, idev, ep1, ep2, ep3, ep4, ep5, ep6)