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!|| matparam_def_mod ../common_source/modules/mat_elem/matparam_def_mod.F90
35!|| prop_param_mod ../common_source/modules/mat_elem/prop_param_mod.F90
36!||====================================================================
37 SUBROUTINE m2law(mat_param ,
38 1 PM, OFF, SIG, EINT,
39 2 RHO, QOLD, EPXE, EPSD,
40 3 VOL, STIFN, DT2T, NELTST,
41 4 ITYPTST, OFFG, GEO, PID,
42 5 AMU, VOL_AVG, MUMAX, MAT,
43 6 NGL, SSP, DVOL, AIRE,
44 7 VNEW, VD2, DELTAX, VIS,
45 8 D1, D2, D3, D4,
46 9 D5, D6, PNEW, PSH,
47 A QNEW, SSP_EQ, SOLD1, SOLD2,
48 B SOLD3, SOLD4, SOLD5, SOLD6,
49 C SIGY, DEFP, DPLA, JLAG,
50 D EPSP, TSTAR, ETSE, MSSA,
51 E DMELS, TEMPEL, SIGBAK, AL_IMP,
52 F SIGNOR, CONDE, DTEL, G_DT,
53 G NEL, IPM, RHOREF, RHOSP,
54 H IPG, DMG, ITY, JTUR,
55 I JTHE, JSPH, ISMSTR, JSMS,
56 J NPG , DPDM , FHEAT, glob_therm,
57 H NUMGEO, IGEO)
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
579 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(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)
Definition m2law.F:58
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