OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
m2law.F File Reference
#include "implicit_f.inc"
#include "comlock.inc"
#include "mvsiz_p.inc"
#include "com08_c.inc"
#include "param_c.inc"
#include "impl1_c.inc"
#include "units_c.inc"
#include "lockon.inc"
#include "lockoff.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine m2law (mat_param, pm, off, sig, eint, rho, qold, epxe, epsd, vol, stifn, dt2t, neltst, ityptst, offg, geo, pid, amu, vol_avg, mumax, mat, ngl, ssp, dvol, aire, vnew, vd2, deltax, vis, d1, d2, d3, d4, d5, d6, pnew, psh, qnew, ssp_eq, sold1, sold2, sold3, sold4, sold5, sold6, sigy, defp, dpla, jlag, epsp, tstar, etse, mssa, dmels, tempel, sigbak, al_imp, signor, conde, dtel, g_dt, nel, ipm, rhoref, rhosp, ipg, dmg, ity, jtur, jthe, jsph, ismstr, jsms, npg, dpdm, fheat, glob_therm, numgeo, igeo)

Function/Subroutine Documentation

◆ m2law()

subroutine m2law ( type (matparam_struct_), intent(in) mat_param,
pm,
off,
sig,
eint,
rho,
qold,
epxe,
epsd,
vol,
stifn,
dt2t,
integer neltst,
integer ityptst,
offg,
geo,
integer, dimension(nel) pid,
amu,
vol_avg,
mumax,
integer, dimension(nel) mat,
integer, dimension(nel) ngl,
ssp,
dvol,
aire,
vnew,
vd2,
deltax,
vis,
d1,
d2,
d3,
d4,
d5,
d6,
pnew,
psh,
qnew,
ssp_eq,
sold1,
sold2,
sold3,
sold4,
sold5,
sold6,
sigy,
defp,
dpla,
integer, intent(in) jlag,
epsp,
tstar,
etse,
mssa,
dmels,
tempel,
sigbak,
al_imp,
signor,
conde,
dtel,
integer g_dt,
integer nel,
integer, dimension(npropmi,*) ipm,
rhoref,
rhosp,
integer ipg,
dmg,
integer, intent(in) ity,
integer, intent(in) jtur,
integer, intent(in) jthe,
integer, intent(in) jsph,
integer, intent(in) ismstr,
integer, intent(in) jsms,
integer, intent(in) npg,
intent(in) dpdm,
intent(inout) fheat,
type (glob_therm_), intent(inout) glob_therm,
integer, intent(in) numgeo,
integer, dimension(n_var_igeo,numgeo), intent(in) igeo )

Definition at line 37 of file m2law.F.

58C-----------------------------------------------
59C M o d u l e s
60C-----------------------------------------------
61 use glob_therm_mod
62 use matparam_def_mod
63 USE prop_param_mod , only : n_var_igeo
64C-----------------------------------------------
65C I m p l i c i t T y p e s
66C-----------------------------------------------
67#include "implicit_f.inc"
68#include "comlock.inc"
69C-----------------------------------------------
70C G l o b a l P a r a m e t e r s
71C-----------------------------------------------
72#include "mvsiz_p.inc"
73C-----------------------------------------------
74C C o m m o n B l o c k s
75C-----------------------------------------------
76#include "com08_c.inc"
77#include "param_c.inc"
78#include "impl1_c.inc"
79#include "units_c.inc"
80C-----------------------------------------------
81C D u m m y A r g u m e n t s
82C-----------------------------------------------
83 INTEGER, INTENT(IN) :: NPG
84 INTEGER, INTENT(IN) :: ISMSTR
85 INTEGER, INTENT(IN) :: JSMS
86 INTEGER, INTENT(IN) :: ITY
87 INTEGER, INTENT(IN) :: JTUR
88 INTEGER, INTENT(IN) :: JTHE
89 INTEGER, INTENT(IN) :: JSPH
90 INTEGER, INTENT(IN) :: JLAG
91 INTEGER,INTENT(IN) :: NUMGEO
92C
93 INTEGER NELTST,ITYPTST,PID(NEL),G_DT,NEL,IPG
94 INTEGER MAT(NEL),NGL(NEL), IPM(NPROPMI,*)
95 my_real :: dt2t
96 my_real pm(npropm,*), off(nel), sig(nel,6), eint(nel), rho(nel), qold(nel),
97 . epxe(nel), epsd(nel), vol(nel), stifn(nel), offg(nel),geo(npropg,*),
98 . mumax(nel),amu(nel),vol_avg(nel)
99 my_real vnew(nel), vd2(nel), deltax(nel), ssp(nel), aire(nel), vis(nel),
100 . psh(nel), pnew(nel),qnew(nel) ,ssp_eq(nel), dvol(nel),
101 . d1(nel), d2(nel), d3(nel), d4(nel), d5(nel), d6(nel),
102 . sold1(mvsiz), sold2(mvsiz), sold3(mvsiz),
103 . sold4(mvsiz), sold5(mvsiz), sold6(mvsiz),
104 . tstar(mvsiz), dpla(mvsiz), epsp(mvsiz), dmg(nel)
105 my_real, DIMENSION(MVSIZ) ,INTENT(IN) :: dpdm
106 my_real, DIMENSION(NEL) ,INTENT(INOUT) :: fheat ! Heat due to plastic work for Heat Equation with lagrangian framework
107 my_real sigy(nel) ,defp(nel),etse(nel), mssa(nel), dmels(nel), tempel(nel),
108 . sigbak(nel,6),al_imp(nel),signor(mvsiz,6),conde(nel),dtel(nel),
109 . rhoref(nel) ,rhosp(nel)
110 type (glob_therm_) ,intent(inout) :: glob_therm
111 type (matparam_struct_) ,intent(in) :: mat_param
112 integer,dimension(n_var_igeo,numgeo),intent(in) :: igeo
113C-----------------------------------------------
114C L o c a l V a r i a b l e s
115C-----------------------------------------------
116 INTEGER :: I, J, ICC,IFORM,ISRATE,VP,IDEV,NINDX
117 INTEGER :: INDX(NEL)
118 my_real rho0, plap,pold,pshift,bulk,p, epmx,cc,cn,epdr,
119 . z3,z4,rhocp,rhocpi,t_ref,t_melt,asrate,
120 . e1, e2, e3, e4, e5,e6, einc, g0,g1,g2,g3,g3h,
121 . scale, dta,mt,
122 . dsxx,dsyy,dszz,dsxy,dsyz,dszx,alpha,hkin,
123 . facq0,fisokin,beta,vm,vm_1,norm_1,ca0,cb,sigm0
124 my_real epd(mvsiz),g(mvsiz), ak(mvsiz), qh(mvsiz),
125 . aj2(mvsiz), dav(mvsiz),ca(mvsiz),sigmx(mvsiz),
126 . sigexx(mvsiz),sigeyy(mvsiz),sigezz(mvsiz),
127 . sigexy(mvsiz),sigeyz(mvsiz),sigezx(mvsiz)
128 my_real, DIMENSION(MVSIZ) :: bidmvsiz
129C-----------------------------------------------
130C B o d y
131C-----------------------------------------------
132 facq0 = one
133C
134 iform = mat_param%iparam(1)
135 icc = mat_param%iparam(2)
136 vp = mat_param%iparam(3)
137 israte = mat_param%iparam(4)
138
139 rho0 = mat_param%rho
140 bulk = mat_param%bulk
141 g0 = mat_param%shear
142
143
144 ca0 = mat_param%uparam(1)
145 cb = mat_param%uparam(2)
146 cn = mat_param%uparam(3)
147 epmx = mat_param%uparam(4)
148 sigm0 = mat_param%uparam(5)
149 cc = mat_param%uparam(6)
150 epdr = mat_param%uparam(7)
151 fisokin = mat_param%uparam(8)
152 asrate = mat_param%uparam(9)
153 z3 = mat_param%uparam(10)
154 asrate = min(one, asrate*dt1)
155
156 t_ref = mat_param%therm%tref
157 t_melt = mat_param%therm%tmelt
158 rhocp = mat_param%therm%rhocp
159 IF (rhocp > zero) THEN
160 rhocpi = one / rhocp
161 ELSE
162 rhocpi = zero
163 END IF
164 pshift = pm(88,mat(1)) ! from EOS
165!
166 IF (iform == 1) THEN ! Zerilli
167 z4 = mat_param%uparam(11)
168 END IF
169!
170 DO i=1,nel
171 g(i) = g0*off(i)
172 ca(i) = ca0
173 sigmx(i)= sigm0
174 etse(i) = one
175 ENDDO
176C
177C------------------------------------------
178C ECROUISSAGE CINE
179C------------------------------------------
180 IF (fisokin > zero) THEN
181 DO i=1,nel
182 sig(i,1)=sig(i,1)-sigbak(i,1)
183 sig(i,2)=sig(i,2)-sigbak(i,2)
184 sig(i,3)=sig(i,3)-sigbak(i,3)
185 sig(i,4)=sig(i,4)-sigbak(i,4)
186 sig(i,5)=sig(i,5)-sigbak(i,5)
187 sig(i,6)=sig(i,6)-sigbak(i,6)
188 ENDDO
189 ENDIF
190C
191 DO i=1,nel
192 p =-third*(sig(i,1)+sig(i,2)+sig(i,3))
193 dav(i)=-third*(d1(i)+d2(i)+d3(i))
194 g1=dt1*g(i)
195 g2=two*g1
196 ssp(i)=sqrt((onep333*g(i)+bulk)/rho0)
197c -- deviatoric elastic stress predictor
198 sig(i,1)=sig(i,1)+p+g2*(d1(i)+dav(i))
199 sig(i,2)=sig(i,2)+p+g2*(d2(i)+dav(i))
200 sig(i,3)=sig(i,3)+p+g2*(d3(i)+dav(i))
201 sig(i,4)=sig(i,4)+g1*d4(i)
202 sig(i,5)=sig(i,5)+g1*d5(i)
203 sig(i,6)=sig(i,6)+g1*d6(i)
204C
205c -- von Mises stress at elastic predictor
206 aj2(i)=half*(sig(i,1)**2+sig(i,2)**2+sig(i,3)**2)
207 . + sig(i,4)**2+sig(i,5)**2+sig(i,6)**2
208 aj2(i)=sqrt(three*aj2(i))
209 ENDDO
210C
211c -- storing elastic predictors for stiffness computation
212 IF (impl_s>0 .OR. fisokin>zero) THEN
213 DO i=1,nel
214 sigexx(i) = sig(i,1)
215 sigeyy(i) = sig(i,2)
216 sigezz(i) = sig(i,3)
217 sigexy(i) = sig(i,4)
218 sigeyz(i) = sig(i,5)
219 sigezx(i) = sig(i,6)
220 ENDDO
221 ENDIF
222C-------------
223C STRAIN RATE (JOHNSON-COOK, ZERILLI-ARMSTRONG)
224C-------------
225 idev = vp - 2
226 CALL mstrain_rate(nel ,israte ,asrate ,epsd ,idev ,
227 . d1, d2, d3, d4, d5, d6)
228
229 epsp(1:nel) = epsd(1:nel)
230 epd(1:nel) = one
231c
232 IF (cc /= zero) THEN
233 IF (vp == 1) THEN
234 DO i=1,nel
235 epd(i)= max(epsd(i),epdr)
236 epd(i)= log(epd(i)/epdr)
237 ENDDO
238 ELSE
239 DO i=1,nel
240 epd(i)= max(epsd(i),em15)
241 epd(i)= log(epd(i)/epdr)
242 ENDDO
243 ENDIF ! VP
244!
245 IF (iform == 0) THEN ! J-C
246 DO i=1,nel
247 mt=max(em15,z3)
248 epd(i)= max(zero,epd(i))
249 epd(i)= (one +cc * epd(i))*(one -tstar(i)**mt)
250 IF (icc==1) sigmx(i) = sigm0*epd(i)
251 ENDDO
252 ELSEIF (iform==1) THEN
253 DO i=1,nel
254 epd(i)= cc*exp((-z3+z4 * epd(i))*tempel(i))
255 IF (icc==1) sigmx(i) = sigm0 + epd(i)
256 ca(i) = ca0 + epd(i)
257 epd(i)=one
258 ENDDO
259 END IF
260 ELSEIF (iform == 0) THEN
261 mt = max(em15,z3)
262 epd(1:nel) = one - tstar(1:nel)**mt
263 IF (icc==1) sigmx(1:nel) = sigm0*epd(1:nel)
264 ENDIF
265C-------------
266C CRITERE
267C-------------
268 IF (fisokin == zero) THEN ! isotropic hardening
269 IF(cn==one) THEN
270 ak(1:nel)= ca(1:nel)+cb*epxe(1:nel)
271 qh(1:nel)= cb*epd(1:nel)
272 ELSE
273 DO i=1,nel
274 IF(epxe(i)>zero) THEN
275 ak(i)=ca(i)+cb*epxe(i)**cn
276 IF(cn>one) THEN
277 qh(i) = (cb*cn*epxe(i)**(cn - one))*epd(i)
278 ELSE
279 qh(i) = (cb*cn/epxe(i)**(one -cn))*epd(i)
280 ENDIF
281 ELSE
282 ak(i)=ca(i)
283 qh(i)=zero
284 ENDIF
285 ENDDO
286 ENDIF
287 DO i=1,nel
288 ak(i)=ak(i)*epd(i)
289 IF(sigmx(i)<ak(i))THEN
290 ak(i)=sigmx(i)
291 qh(i)=zero
292 ENDIF
293 sigy(i) = ak(i)
294 IF(epxe(i)>epmx)THEN
295 ak(i)=zero
296 qh(i)=zero
297 ENDIF
298 ENDDO
299!
300 ELSE ! kinematic hardening
301!
302 DO i=1,nel
303 beta = one-fisokin
304 ! SIGY is used for hourglass stress compute--
305 IF(cn==one) THEN
306 sigy(i) = ca(i)+cb*epxe(i)
307 ak(i)= ca(i)+beta*cb*epxe(i)
308 qh(i)= cb*epd(i)
309 ELSE
310 IF(epxe(i)>zero) THEN
311 sigy(i)=ca(i)+cb*epxe(i)**cn
312 ak(i)=ca(i)+beta*cb*epxe(i)**cn
313 IF(cn>one) THEN
314 qh(i)= (cb*cn*epxe(i)**(cn - one))*epd(i)
315 ELSE
316 qh(i)= (cb*cn/epxe(i)**(one -cn))*epd(i)
317 ENDIF
318 ELSE
319 ak(i)=ca(i)
320 sigy(i)=ca(i)
321 qh(i)=zero
322 ENDIF
323 ENDIF
324 ak(i)=ak(i)*epd(i)
325 sigy(i)=sigy(i)*epd(i)
326 IF(sigmx(i)<ak(i))THEN
327 ak(i)=sigmx(i)
328 qh(i)=zero
329 ENDIF
330 sigy(i)=min(sigy(i),sigmx(i))
331 IF(epxe(i)>epmx)THEN
332 ak(i)=zero
333 qh(i)=zero
334 ENDIF
335 END DO ! I=1,NEL
336 END IF ! (FISOKIN == 0 ) THEN
337!---------------------------------------------------------
338 IF (impl_s==0) THEN
339 DO i=1,nel
340 scale= min(one,ak(i)/ max(aj2(i),em15))
341C plastic strain increment.
342 dpla(i)=(one-scale)*aj2(i)/max(three*g(i)+qh(i),em15)
343C actual yield stress.
344 ak(i)=ak(i)+(one - fisokin)*dpla(i)*qh(i)
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)+dpla(i)
353 ENDDO
354!
355 ELSE ! implicit : nonlinear hardening requires iterations in radial return
356 CALL m2iter_imp(sig, epxe, aj2, g,
357 1 ca, cb, cn, epd,
358 2 sigmx, epmx, dpla, ak,
359 3 qh, sigy, fisokin, nel)
360!
361 END IF ! IMPL_S
362c-----------------------------------------
363 IF (fisokin > zero) THEN
364 DO i=1,nel
365 dsxx = sigexx(i) - sig(i,1)
366 dsyy = sigeyy(i) - sig(i,2)
367 dszz = sigezz(i) - sig(i,3)
368 dsxy = sigexy(i) - sig(i,4)
369 dsyz = sigeyz(i) - sig(i,5)
370 dszx = sigezx(i) - sig(i,6)
371C
372 hkin = two_third*fisokin*qh(i)
373 alpha = hkin/max(two*g(i)+hkin,em15)
374C ..updates back stresses
375 sigbak(i,1) = sigbak(i,1) + alpha*dsxx
376 sigbak(i,2) = sigbak(i,2) + alpha*dsyy
377 sigbak(i,3) = sigbak(i,3) + alpha*dszz
378 sigbak(i,4) = sigbak(i,4) + alpha*dsxy
379 sigbak(i,5) = sigbak(i,5) + alpha*dsyz
380 sigbak(i,6) = sigbak(i,6) + alpha*dszx
381C ..gets stresses from shifted stresses and back stresses
382 sig(i,1)=sig(i,1) + sigbak(i,1)
383 sig(i,2)=sig(i,2) + sigbak(i,2)
384 sig(i,3)=sig(i,3) + sigbak(i,3)
385 sig(i,4)=sig(i,4) + sigbak(i,4)
386 sig(i,5)=sig(i,5) + sigbak(i,5)
387 sig(i,6)=sig(i,6) + sigbak(i,6)
388 ENDDO
389 END IF !(FISOKIN > 0 ) THEN
390c-----------------------------------------
391!
392 bidmvsiz(1:mvsiz) = zero
393 IF (jsph==0) THEN
394 CALL mqviscb(
395 1 pm, off, rho, bidmvsiz,
396 2 bidmvsiz,ssp, bidmvsiz,stifn,
397 3 dt2t, neltst, ityptst, aire,
398 4 offg, geo, pid, vnew,
399 5 vd2, deltax, vis, d1,
400 6 d2, d3, pnew, psh,
401 7 mat, ngl, qnew, ssp_eq,
402 8 vol, mssa, dmels, igeo,
403 9 facq0, conde, dtel, g_dt,
404 a ipm, rhoref, rhosp, nel,
405 b ity, ismstr, jtur, jthe,
406 c jsms, npg , glob_therm)
407 ELSE
408 CALL mdtsph(
409 1 pm, off, rho, bidmvsiz,
410 2 bidmvsiz,bidmvsiz,stifn, dt2t,
411 3 neltst, ityptst, offg, geo,
412 4 pid, mumax, ssp, vnew,
413 5 vd2, deltax, vis, d1,
414 6 d2, d3, pnew, psh,
415 7 mat, ngl, qnew, ssp_eq,
416 8 g_dt, dtel, nel, ity,
417 9 jtur, jthe)
418 ENDIF
419C
420 dta = half*dt1
421C
422 nindx = 0
423 indx(1:nel) = 0
424 DO i=1,nel
425 IF ((epxe(i) > epmx).AND.(dmg(i)==zero)) THEN
426 nindx = nindx + 1
427 indx(nindx) = i
428 dmg(i) = one
429 ENDIF
430 ENDDO
431
432 IF (mat_param%IEOS == 0) THEN ! add pressure to the deviatoric stress
433 DO i=1,nel
434 pnew(i) = bulk*amu(i)
435 sig(i,1) =(sig(i,1) - pnew(i))*off(i)
436 sig(i,2) =(sig(i,2) - pnew(i))*off(i)
437 sig(i,3) =(sig(i,3) - pnew(i))*off(i)
438 sig(i,4) = sig(i,4) * off(i)
439 sig(i,5) = sig(i,5) * off(i)
440 sig(i,6) = sig(i,6) * off(i)
441 ENDDO
442!
443 DO i=1,nel
444 e1 = d1(i)*(sold1(i)+sig(i,1))
445 e2 = d2(i)*(sold2(i)+sig(i,2))
446 e3 = d3(i)*(sold3(i)+sig(i,3))
447 e4 = d4(i)*(sold4(i)+sig(i,4))
448 e5 = d5(i)*(sold5(i)+sig(i,5))
449 e6 = d6(i)*(sold6(i)+sig(i,6))
450 einc = vol_avg(i)*(e1+e2+e3+e4+e5+e6)*dta - half*dvol(i)*(qold(i)+qnew(i))
451 eint(i) = (eint(i)+einc*off(i)) / max(em15,vol(i))
452 ENDDO
453
454 ELSE
455 ! if EOS is used, material law calculates only deviatoric stress tensor
456 ! sound speed depends on pressure derivative over volume change
457 ! calculated in EOS
458 DO i = 1, nel
459 ssp(i) = sqrt((dpdm(i) + four*g(i)/three)/rho0)
460!
461 pold = (sold1(i)+sold2(i)+sold3(i)) * third
462 sold1(i) = sold1(i) - pold
463 sold2(i) = sold2(i) - pold
464 sold3(i) = sold3(i) - pold
465 e1 = d1(i) * (sold1(i)+sig(i,1))
466 e2 = d2(i) * (sold2(i)+sig(i,2))
467 e3 = d3(i) * (sold3(i)+sig(i,3))
468 e4 = d4(i) * (sold4(i)+sig(i,4))
469 e5 = d5(i) * (sold5(i)+sig(i,5))
470 e6 = d6(i) * (sold6(i)+sig(i,6))
471 einc = vol_avg(i) * (e1+e2+e3+e4+e5+e6) * dta
472 eint(i) = eint(i) + (einc+half*dvol(i)*(pold-pshift-qold(i)-qnew(i)))*off(i)
473 ENDDO
474 END IF
475C
476 DO i=1,nel
477 qold(i) = qnew(i)
478 defp(i) = epxe(i)
479 sigy(i) = max(sigy(i),ak(i))
480 ENDDO
481 DO i=1,nel
482 IF (dpla(i) > 0) etse(i)= half*qh(i)*off(i)/max(g(i),em15)
483 ENDDO
484!----------------------------------------------------------------
485! implicit
486!----------------------------------------------------------------
487 IF (impl_s>0) THEN
488 IF (ikt==0) RETURN
489 IF (fisokin == zero) THEN
490 DO i=1,nel
491 IF (dpla(i)>0)THEN
492 IF(cn==one) THEN
493 qh(i)= cb*epd(i)
494 ELSEIF(cn>one) THEN
495 qh(i)= (cb*cn*epxe(i)**(cn - one))*epd(i)
496 ELSE
497 qh(i)= (cb*cn/epxe(i)**(one -cn))*epd(i)
498 ENDIF
499 ENDIF
500 ENDDO
501 ENDIF
502C
503 DO i = 1,nel
504 IF (dpla(i)>zero) THEN
505
506c ..... Von Mises stress at the elastic predictor (point B)
507c ..... SIGEXX, etc. are deviatoric stresses
508 vm =half*(sigexx(i)**2+sigeyy(i)**2+sigezz(i)**2)
509 1 +sigexy(i)**2+sigeyz(i)**2+sigezx(i)**2
510 vm_1 =one/sqrt(three*vm)
511 g3 = three*g(i)
512 g3h = max(g3+qh(i),em15)
513 scale = max(zero,one-g3h*dpla(i)*vm_1)
514c ..... NORM_1 normalizes deviatoric stresses, includes consistent
515c ..... stiffness matrix parameter beta, von Mises at B, and two_pmi
516 norm_1=g3*vm_1*sqrt(scale/g3h)
517c ..... Deviatoric stresses "normalized"
518 signor(i,1)=sigexx(i)*norm_1
519 signor(i,2)=sigeyy(i)*norm_1
520 signor(i,3)=sigezz(i)*norm_1
521 signor(i,4)=sigexy(i)*norm_1
522 signor(i,5)=sigeyz(i)*norm_1
523 signor(i,6)=sigezx(i)*norm_1
524
525c ..... Parameter alpha of consistent matrix
526 al_imp(i)= one - g3*dpla(i)*vm_1
527 ELSE
528 al_imp(i)=one
529 ENDIF
530 ENDDO
531 ENDIF
532!---------------------------------------------------------
533 ! end implicit
534!---------------------------------------------------------
535! update and filter plastic strain rate for VP=1
536!---------------------------------------------------------
537!
538 IF (vp== 1) THEN
539 DO i=1,nel
540 plap = dpla(i) / max(em20,dt1)
541 epsd(i) = asrate * plap + (one - asrate) * epsd(i)
542 ENDDO
543 ENDIF
544C----------------------------------------------
545C TEMPERATURE (Heating due to plastic work)
546C----------------------------------------------
547 IF (jthe < 0 .AND. jlag /= 0) THEN
548 DO i=1,nel
549 fheat(i) = fheat(i) + sigy(i)*dpla(i)*vol(i)
550 ENDDO
551 ELSEIF (rhocp > zero) THEN
552 DO i=1,nel
553 tempel(i) = tempel(i) + sigy(i)*dpla(i) * rhocpi
554 ! temperature and internal energy must be incremented consistantly
555 ! internal energy is incremented later in parent subroutine (mmain)
556 ! with total energy deformation which already includes plastic work
557 ! Edef = 0.5 * VOL * sum ( sig.eps_dot , i=1..6)
558 ! so internal energy and temperature remain consistant
559 ENDDO
560 END IF
561!---------------------------------------------------------
562 ! Printout element deletion
563 IF (nindx > 0) THEN
564 DO j=1,nindx
565#include "lockon.inc"
566 WRITE(iout, 1000) ngl(indx(j)),ipg
567 WRITE(istdo,1100) ngl(indx(j)),ipg,tt
568#include "lockoff.inc"
569 ENDDO
570 ENDIF
571!---------------------------------------------------------------
572 1000 FORMAT(1x,'EXCEEDED EPS_MAX ON SOLID ELEMENT NUMBER ',i10,
573 . ': DEVIATORIC STRESS SET TO 0 ON INTEGRATION POINT ',i5 )
574 1100 FORMAT(1x,'EXCEEDED EPS_MAX ON SOLID ELEMENT NUMBER ',i10,
575 . ': DEVIATORIC STRESS SET TO 0 ON INTEGRATION POINT ',i5 ,
576 . ' AT TIME :',g11.4)
577!-----------
578 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
dtel
subroutine dtel(ssp, pm, geo, pid, mat, rho0, vis, deltax, aire, vol, dtx)
Definition dtel.F:46
alpha
#define alpha
Definition eval.h:35
m2iter_imp
subroutine m2iter_imp(sig, epxe, aj2, g, ca, cb, cn, epd, sigmx, epmx, dpla1, ak, qh, sigy, fisokin, nel)
Definition m2iter_imp.F:33
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
mdtsph
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)
Definition mdtsph.F:47
mqviscb
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)
Definition mqviscb.F:56
mstrain_rate
subroutine mstrain_rate(nel, israte, asrate, epsd, idev, ep1, ep2, ep3, ep4, ep5, ep6)
Definition mstrain_rate.F:33