OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
m2law.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| m2law ../engine/source/materials/mat/mat002/m2law.F
25!||--- called by ------------------------------------------------------
26!|| mmain ../engine/source/materials/mat_share/mmain.F90
27!||--- calls -----------------------------------------------------
28!|| m2iter_imp ../engine/source/materials/mat/mat002/m2iter_imp.F
29!|| mdtsph ../engine/source/materials/mat_share/mdtsph.F
30!|| mqviscb ../engine/source/materials/mat_share/mqviscb.F
31!|| mstrain_rate ../engine/source/materials/mat_share/mstrain_rate.F
32!||--- uses -----------------------------------------------------
33!|| glob_therm_mod ../common_source/modules/mat_elem/glob_therm_mod.F90
34!||====================================================================
35 SUBROUTINE m2law(
36 1 PM, OFF, SIG, EINT,
37 2 RHO, QOLD, EPXE, EPSD,
38 3 VOL, STIFN, DT2T, NELTST,
39 4 ITYPTST, OFFG, GEO, PID,
40 5 AMU, VOL_AVG, MUMAX, MAT,
41 6 NGL, SSP, DVOL, AIRE,
42 7 VNEW, VD2, DELTAX, VIS,
43 8 D1, D2, D3, D4,
44 9 D5, D6, PNEW, PSH,
45 A QNEW, SSP_EQ, SOLD1, SOLD2,
46 B SOLD3, SOLD4, SOLD5, SOLD6,
47 C SIGY, DEFP, DPLA,
48 D EPSP, TSTAR, ETSE, MSSA,
49 E DMELS, TEMPEL, SIGBAK, AL_IMP,
50 F SIGNOR, CONDE, DTEL, G_DT,
51 G NEL, IPM, RHOREF, RHOSP,
52 H IPG, DMG, ITY, JTUR,
53 I JTHE, JSPH, ISMSTR, JSMS,
54 J PLAP, NPG , IEOS , DPDM ,
55 K FHEAT, glob_therm, JLAG)
56C-----------------------------------------------
57C M o d u l e s
58C-----------------------------------------------
59 use glob_therm_mod
60C-----------------------------------------------
61C I m p l i c i t T y p e s
62C-----------------------------------------------
63#include "implicit_f.inc"
64#include "comlock.inc"
65C-----------------------------------------------
66C G l o b a l P a r a m e t e r s
67C-----------------------------------------------
68#include "mvsiz_p.inc"
69C-----------------------------------------------
70C C o m m o n B l o c k s
71C-----------------------------------------------
72#include "com08_c.inc"
73#include "param_c.inc"
74#include "impl1_c.inc"
75#include "units_c.inc"
76C-----------------------------------------------
77C D u m m y A r g u m e n t s
78C-----------------------------------------------
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
88C
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
92 my_real :: DT2T
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 ! Heat due to plastic work for Heat Equation with lagrangian framework
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
108C-----------------------------------------------
109C L o c a l V a r i a b l e s
110C-----------------------------------------------
111 INTEGER :: I, J, ICC,IRTY,ISRATE,VP,IDEV,MX,IBID,NINDX
112 INTEGER :: INDX(NEL)
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
124C-----------------------------------------------
125C B o d y
126C-----------------------------------------------
127 facq0 = one
128C
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 ! Zerilli
151 rhocpi = pm(53,mx)
152 t_ref = pm(54,mx)
153 z4 = pm(52,mx)
154 ELSE ! Johnson-Cook
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
167C
168C------------------------------------------
169C ECROUISSAGE CINE
170C------------------------------------------
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
181C
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)
188c -- deviatoric elastic stress predictor
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)
195C
196c -- von Mises stress at elastic predictor
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
201C
202c -- storing elastic predictors for stiffness computation
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
213C-------------
214C STRAIN RATE (JOHNSON-COOK, ZERILLI-ARMSTRONG)
215C-------------
216 idev = vp - 2
217 CALL mstrain_rate(nel ,israte ,asrate ,epsd ,idev ,
218 . d1, d2, d3, d4, d5, d6)
219
220 epsp(1:nel) = epsd(1:nel)
221 epd(1:nel) = one
222c
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 ! VP
235!
236 IF (irty == 0) THEN ! J-C
237 DO i=1,nel
238 mt=max(em15,z3)
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
252 mt = max(em15,z3)
253 epd(1:nel) = one - tstar(1:nel)**mt
254 ENDIF
255C-------------
256C CRITERE
257C-------------
258 IF (fisokin == zero) THEN ! isotropic hardening
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 ! kinematic hardening
291!
292 DO i=1,nel
293 beta = one-fisokin
294 ! SIGY is used for hourglass stress compute--
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 ! I=1,NEL
326 END IF ! (FISOKIN == 0 ) THEN
327!---------------------------------------------------------
328 IF (impl_s==0) THEN
329 DO i=1,nel
330 scale= min(one,ak(i)/ max(aj2(i),em15))
331C plastic strain increment.
332 dpla(i)=(one-scale)*aj2(i)/max(three*g(i)+qh(i),em15)
333C actual yield stress.
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 ! implicit : nonlinear hardening requires iterations in radial return
346 CALL m2iter_imp(sig, epxe, aj2, g,
347 1 ca, cb, cn, epd,
348 2 sigmx, epmx, dpla, ak,
349 3 qh, sigy, fisokin, nel)
350!
351 END IF ! IMPL_S
352c-----------------------------------------
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)
361C
362 hkin = two_third*fisokin*qh(i)
363 alpha = hkin/max(two*g(i)+hkin,em15)
364C ..updates back stresses
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
371C ..gets stresses from shifted stresses and back stresses
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 !(FISOKIN > 0 ) THEN
380c-----------------------------------------
381!
382 bidmvsiz(1:mvsiz) = zero
383 IF (jsph==0) THEN
384 CALL mqviscb(
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
398 CALL mdtsph(
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
409C
410 dta = half*dt1
411C
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 ! add pressure to the deviatoric stress
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 ! if EOS is used, material law calculates only deviatoric stress tensor
446 ! sound speed depends on pressure derivative over volume change
447 ! calculated in EOS
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
465C
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! implicit
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
492C
493 DO i = 1,nel
494 IF (dpla(i)>zero) THEN
495
496c ..... Von Mises stress at the elastic predictor (point B)
497c ..... SIGEXX, etc. are deviatoric stresses
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)
504c ..... NORM_1 normalizes deviatoric stresses, includes consistent
505c ..... stiffness matrix parameter beta, von Mises at B, and two_pmi
506 norm_1=g3*vm_1*sqrt(scale/g3h)
507c ..... Deviatoric stresses "normalized"
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
515c ..... Parameter alpha of consistent matrix
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 ! end implicit
524!---------------------------------------------------------
525! update and filter plastic strain rate for VP=1
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
534C----------------------------------------------
535C TEMPERATURE (Heating due to plastic work)
536C----------------------------------------------
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 ! temperature and internal energy must be incremented consistantly
545 ! internal energy is incremented later in parent subroutine (mmain)
546 ! with total energy deformation which already includes plastic work
547 ! Edef = 0.5 * VOL * sum ( sig.eps_dot , i=1..6)
548 ! so internal energy and temperature remain consistant
549 ENDDO
550 END IF
551!---------------------------------------------------------
552 ! Printout element deletion
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,tt
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
569 END
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
m2law
subroutine m2law(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, 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, plap, npg, ieos, dpdm, fheat, glob_therm, jlag)
Definition m2law.F:56
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