OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat104c_ldam_newton.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!|| mat104c_ldam_newton ../engine/source/materials/mat/mat104/mat104c_ldam_newton.F
25!||--- called by ------------------------------------------------------
26!|| sigeps104c ../engine/source/materials/mat/mat104/sigeps104c.F
27!||====================================================================
28 SUBROUTINE mat104c_ldam_newton(
29 1 NEL ,NGL ,NUPARAM ,NUVAR ,
30 2 TIME ,TIMESTEP,UPARAM ,UVAR ,JTHE ,OFF ,
31 3 GS ,RHO ,PLA ,DPLA ,EPSD ,SOUNDSP ,
32 4 DEPSXX ,DEPSYY ,DEPSXY ,DEPSYZ ,DEPSZX ,
33 5 SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
34 6 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,THKLY ,
35 7 THK ,SIGY ,ET ,TEMPEL ,DPLA_NL ,DMG ,
36 8 TEMP ,SEQ ,PLA_NL ,L_PLANL ,PLAP_NL ,L_EPSDNL)
37 !=======================================================================
38 ! Implicit types
39 !=======================================================================
40#include "implicit_f.inc"
41 !=======================================================================
42 ! Common
43 !=======================================================================
44#include "com01_c.inc"
45 !=======================================================================
46 ! Dummy arguments
47 !=======================================================================
48 INTEGER NEL,NUPARAM,NUVAR,JTHE
49 INTEGER ,DIMENSION(NEL), INTENT(IN) :: NGL
50 INTEGER, INTENT(IN) :: L_PLANL,L_EPSDNL
51 my_real
52 . TIME,TIMESTEP
53 my_real,DIMENSION(NUPARAM), INTENT(IN) ::
54 . UPARAM
55 my_real,DIMENSION(NEL), INTENT(IN) ::
56 . rho,tempel,
57 . depsxx,depsyy,depsxy,depsyz,depszx,
58 . sigoxx,sigoyy,sigoxy,sigoyz,sigozx,
59 . gs,thkly,dpla_nl
60 my_real, DIMENSION(NEL*L_PLANL), INTENT(IN) ::
61 . pla_nl
62 my_real, DIMENSION(NEL*L_EPSDNL), INTENT(IN) ::
63 . plap_nl
64c
65 my_real ,DIMENSION(NEL), INTENT(OUT) ::
66 . soundsp,sigy,et,
67 . signxx,signyy,signxy,signyz,signzx
68c
69 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
70 . pla,dpla,epsd,off,thk,temp,seq
71 my_real ,DIMENSION(NEL,6), INTENT(INOUT) ::
72 . dmg
73 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) ::
74 . uvar
75 !=======================================================================
76 ! Local Variables
77 !=======================================================================
78 INTEGER I,II,IGURSON,ITER,NITER,NINDX,INDEX(NEL)
79c
80 my_real ::
81 . young,bulk,lam,g,g2,nu,cdr,kdr,hard,yld0,qvoce,bvoce,jcc,
82 . epsp0,mtemp,tini,tref,eta,cp,dpis,dpad,asrate,afiltr,hkhi,
83 . q1,q2,ed,an,epn,kw,fr,fc,f0,a11,a12,nnu,dti
84 my_real ::
85 . h,ldav,trdfds,sigvm,omega,
86 . dtemp,fcosh,fsinh,dpdt,dlam,ddep
87 my_real ::
88 . dsdrdj2,dsdrdj3,
89 . dj3dsxx,dj3dsyy,dj3dsxy,dj3dszz,
90 . dj2dsxx,dj2dsyy,dj2dsxy,normsig,
91 . dfdsxx,dfdsyy,dfdsxy,dyld_dpla_nl,
92 . normxx,normyy,normxy,normzz,
93 . sdpla,dphi_dtrsig,sig_dfdsig,dfdsig2,sdv_dfdsig,
94 . dphi_dsig,dphi_dyld,dphi_dfdr,df_dfs,dfs_dft,dphi_dft,
95 . dphi_dfs,dfn_dlam,dfsh_dlam,dfg_dlam,dft_dlam,
96 . dfn,dfsh,dfg,dft,dyld_dpla,dyld_dtemp,dtemp_dlam
97C
98 my_real, DIMENSION(NEL) ::
99 . dsigxx,dsigyy,dsigxy,trsig,trdep,
100 . sxx,syy,sxy,szz,sigm,j2,j3,sigdr,yld,weitemp,
101 . hardp,fhard,frate,ftherm,dtherm,fdr,phi0,triax,
102 . fdam,phi,ft,fs,fg,fn,fsh,dpla_dlam,dphi_dlam,dezz,etat,
103 . sigdr2,yld2i,dpxx,dpyy,dpzz,dpxy,dlam_nl
104 .
105 !=======================================================================
106 ! DRUCKER - VOCE - JOHNSON-COOK MATERIAL LAW WITH GURSON DAMAGE
107 ! USING LOCAL DAMAGE OR NON-LOCAL MICROMORPHIC METHOD
108 !=======================================================================
109c
110 !=======================================================================
111 ! - INITIALISATION OF COMPUTATION ON TIME STEP
112 !=======================================================================
113 ! Recovering model parameters
114 ! Elastic parameters
115 young = uparam(1) ! Young modulus
116 bulk = uparam(2) ! Bulk modulus
117 g = uparam(3) ! Shear modulus
118 g2 = uparam(4) ! 2*Shear modulus
119 lam = uparam(5) ! Lambda Hooke parameter
120 nu = uparam(6) ! Poisson ration
121 nnu = uparam(7) ! NU/(1-NU)
122 a11 = uparam(9) ! Diagonal term, elastic matrix in plane stress
123 a12 = uparam(10) ! Non-diagonal term, elastic matrix in plane stress
124 ! Plastic criterion and hardening parameters [Drucker, 1948]
125 cdr = uparam(12) ! Drucker coefficient
126 kdr = uparam(13) ! Drucker 1/K coefficient
127 tini = uparam(14) ! Initial temperature
128 hard = uparam(15) ! Linear hardening
129 yld0 = uparam(16) ! Initial yield stress
130 qvoce = uparam(17) ! 1st Voce parameter
131 bvoce = uparam(18) ! 2nd Voce parameter
132 ! Strain-rate dependence parameters
133 jcc = uparam(20) ! Johnson-Cook strain rate coefficient
134 epsp0 = uparam(21) ! Johnson-Cook reference strain rate
135 ! Thermal softening and self-heating parameters
136 mtemp = uparam(22) ! Thermal softening slope
137 tref = uparam(23) ! Reference temperature
138 eta = uparam(24) ! Taylor-Quinney coefficient
139 cp = uparam(25) ! Thermal mass capacity
140 dpis = uparam(26) ! Isothermal plastic strain rate
141 dpad = uparam(27) ! Adiabatic plastic strain rate
142 ! Plastic strain-rate filtering parameters
143 asrate = uparam(28) ! Plastic strain rate filtering frequency
144 afiltr = asrate*timestep/(asrate*timestep + one)
145 dti = one / max(timestep, em20)
146c
147 ! Gurson damage model parameters parameters
148 igurson = nint(uparam(30)) ! Gurson switch flag:
149 ! = 0 => no damage model
150 ! = 1 => local damage model
151 ! = 2 => non local (Forest - micromorphic) damage model
152 ! = 3 => non local (Peerlings) damage model
153 q1 = uparam(31) ! Gurson yield criterion 1st parameter
154 q2 = uparam(32) ! Gurson yield criterion 2nd parameter
155 ed = uparam(34) ! Plastic strain trigger for nucleation
156 an = uparam(35) ! Nucleation rate
157 kw = uparam(36) ! Shear coefficient (Nahshon & Hutchinson)
158 fr = uparam(37) ! Failure void volume fraction
159 fc = uparam(38) ! Critical void volume fraction
160 f0 = uparam(39) ! Initial void volume fraction
161 hkhi = uparam(40) ! Micromorphic penalty parameter
162c
163 ! Recovering internal variables
164 DO i=1,nel
165 ! If the element is failing
166 IF (off(i) < em03) off(i) = zero
167 IF (off(i) < one) off(i) = off(i)*four_over_5
168 ! Damage variables
169 fg(i) = dmg(i,2) ! Growth damage
170 fn(i) = dmg(i,3) ! Nucleation damage
171 fsh(i) = dmg(i,4) ! Shear damage
172 ft(i) = dmg(i,5) ! Total damage
173 fs(i) = dmg(i,6) ! Effective damage
174 ! Standard inputs
175 dpla(i) = zero ! Initialization of the plastic strain increment
176 et(i) = one ! Initialization of hourglass coefficient
177 hardp(i) = zero ! Initialization of hardening modulus
178 dezz(i) = zero ! Initialization of the transverse strain
179 ENDDO
180c
181 ! Initialization of damage, temperature and self-heating weight factor
182 IF (time == zero) THEN
183 temp(1:nel) = tini
184 IF (isigi == 0) THEN
185 dmg(1:nel,5) = f0
186 ft(1:nel) = f0
187 dmg(1:nel,1) = f0/fr
188 IF (f0<fc) THEN
189 dmg(1:nel,6) = f0
190 ELSE
191 dmg(1:nel,6) = fc + (one/q1-fc)*(f0-fc)/(fr-fc)
192 ENDIF
193 fs(1:nel) = dmg(1:nel,6)
194 ENDIF
195 ENDIF
196 IF (cp > zero) THEN
197 IF (jthe == 0) THEN
198 DO i=1,nel
199 ! Computation of the weight factor
200 IF (epsd(i) < dpis) THEN
201 weitemp(i) = zero
202 ELSEIF (epsd(i) > dpad) THEN
203 weitemp(i) = one
204 ELSE
205 weitemp(i) = ((epsd(i)-dpis)**2 )
206 . * (three*dpad - two*epsd(i) - dpis)
207 . / ((dpad-dpis)**3)
208 ENDIF
209 ENDDO
210 ELSE
211 temp(1:nel) = tempel(1:nel)
212 ENDIF
213 ENDIF
214c
215 ! Computation of the initial yield stress
216 DO i = 1,nel
217 ! a) - Hardening law
218 fhard(i) = yld0 + hard*pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
219 IF (igurson == 2) THEN
220 IF (pla_nl(i) - pla(i) < zero) THEN
221 fhard(i) = fhard(i) - hkhi*(pla_nl(i) - pla(i))
222 ENDIF
223 ENDIF
224 ! b) - Correction factor for strain-rate dependence (Johnson-Cook)
225 frate(i) = one
226 IF (epsd(i) > epsp0) frate(i) = frate(i) + jcc*log(epsd(i)/epsp0)
227 ! c) - Correction factor for thermal softening
228 ftherm(i) = one
229 IF (cp > zero) ftherm(i) = one - mtemp * (temp(i) - tref)
230 ! d) - Computation of the yield function and value check
231 yld(i) = fhard(i)*frate(i)*ftherm(i)
232 ! e) - Checking values
233 yld(i) = max(em10, yld(i))
234 ENDDO
235c
236 !========================================================================
237 ! - COMPUTATION OF TRIAL VALUES
238 !========================================================================
239 DO i=1,nel
240c
241 ! Computation of the trial stress tensor
242 signxx(i) = sigoxx(i) + (a11*depsxx(i) + a12*depsyy(i))
243 signyy(i) = sigoyy(i) + (a11*depsyy(i) + a12*depsxx(i))
244 signxy(i) = sigoxy(i) + (depsxy(i)*g)
245 signyz(i) = sigoyz(i) + (depsyz(i)*gs(i))
246 signzx(i) = sigozx(i) + (depszx(i)*gs(i))
247 ! Computation of the trace of the trial stress tensor
248 trsig(i) = signxx(i) + signyy(i)
249 sigm(i) = -trsig(i) * third
250 ! Computation of the deviatoric trial stress tensor
251 sxx(i) = signxx(i) + sigm(i)
252 syy(i) = signyy(i) + sigm(i)
253 szz(i) = sigm(i)
254 sxy(i) = signxy(i)
255 dezz(i) = -nnu*depsxx(i) - nnu*depsyy(i)
256 ! Second deviatoric invariant
257 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) + sxy(i)**2
258 ! Third deviatoric invariant
259 j3(i) = sxx(i)*syy(i)*szz(i) - szz(i)*sxy(i)**2
260 ! Drucker equivalent stress
261 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
262 ! Checking equivalent stress values
263 IF (fdr(i) > zero) THEN
264 sigdr(i) = kdr * exp((one/six)*log(fdr(i))) ! FDR(I)**(1/6)
265 ELSE
266 sigdr(i) = zero
267 ENDIF
268 ! Computation of the stress triaxiality and the etaT factor
269 IF (sigdr(i)>zero) THEN
270 triax(i) = (trsig(i)*third)/sigdr(i)
271 ELSE
272 triax(i) = zero
273 ENDIF
274 IF (trsig(i)<zero) THEN
275 etat(i) = zero
276 ELSE
277 etat(i) = one
278 ENDIF
279 ENDDO
280c
281 !========================================================================
282 ! - COMPUTATION OF YIELD FONCTION
283 !========================================================================
284 DO i=1,nel
285 fdam(i) = two*q1*fs(i)*cosh(q2*etat(i)*trsig(i)/yld(i)/two) - (q1*fs(i))**2
286 phi(i) = (sigdr(i) / yld(i))**2 - one + fdam(i)
287 ENDDO
288c
289 ! Checking plastic behavior for all elements
290 nindx = 0
291 DO i=1,nel
292 IF ((phi(i) >= zero).AND.(off(i) == one).AND.(ft(i)<fr)) THEN
293 nindx=nindx+1
294 index(nindx)=i
295 ENDIF
296 ENDDO
297c
298 !====================================================================
299 ! - PLASTIC CORRECTION WITH CUTTING PLANE METHOD (SEMI-IMPLICIT)
300 !====================================================================
301c
302 ! Number of iterations
303 niter = 3
304c
305 ! Loop over yielding elements
306#include "vectorize.inc"
307 DO ii=1,nindx
308c
309 ! Number of the element with plastic behaviour
310 i = index(ii)
311c
312 ! Initialization of the iterative Newton procedure
313 sigdr2(i) = sigdr(i)**2
314 yld2i(i) = one / yld(i)**2
315 dpxx(i) = zero
316 dpyy(i) = zero
317 dpzz(i) = zero
318 dpxy(i) = zero
319 ENDDO
320c
321 ! Loop over the iterations
322 DO iter = 1, niter
323#include "vectorize.inc"
324 DO ii=1,nindx
325 i = index(ii)
326c
327 ! Note: in this part, the purpose is to compute for each iteration
328 ! a plastic multiplier allowing to update internal variables to satisfy
329 ! the consistency condition using the cutting plane semi-implicit
330 ! iterative procedure.
331 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
332 ! -> PHI : current value of yield function (known)
333 ! -> dphi_dlambda : derivative of phi with respect to dlambda by taking
334 ! into account of internal variables kinetic :
335 ! plasticity, temperature and damage (to compute)
336c
337 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
338 !-------------------------------------------------------------
339c
340 ! Derivative with respect to the equivalent stress and trace
341 yld2i(i) = one/(yld(i)**2)
342 dphi_dsig = two*sigdr(i)*yld2i(i)
343 fsinh = sinh(q2*etat(i)*trsig(i)/(yld(i)*two))
344 dphi_dtrsig = q1*q2*etat(i)*fs(i)*fsinh/yld(i)
345c
346 ! Computation of the Eulerian norm of the stress tensor
347 normsig = sqrt(signxx(i)*signxx(i)
348 . + signyy(i)*signyy(i)
349 . + two*signxy(i)*signxy(i))
350 normsig = max(normsig,one)
351c
352 ! Derivative with respect to Fdr
353 fdr(i) = (j2(i)/(normsig**2))**3 - cdr*((j3(i)/(normsig**3))**2)
354 dphi_dfdr = dphi_dsig*kdr*(one/six)*exp(-(five/six)*log(fdr(i)))
355 dsdrdj2 = dphi_dfdr*three*(j2(i)/(normsig**2))**2
356 dsdrdj3 = -dphi_dfdr*two*cdr*(j3(i)/(normsig**3))
357 ! dJ3/dSig
358 dj3dsxx = two_third*(syy(i)*szz(i))/(normsig**2)
359 . - third*(sxx(i)*szz(i))/(normsig**2)
360 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
361 dj3dsyy = - third*(syy(i)*szz(i))/(normsig**2)
362 . + two_third*(sxx(i)*szz(i))/(normsig**2)
363 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
364 dj3dszz = - third*(syy(i)*szz(i))/(normsig**2)
365 . - third*(sxx(i)*szz(i))/(normsig**2)
366 . + two_third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
367 dj3dsxy = two*(sxx(i)*sxy(i) + sxy(i)*syy(i))/(normsig**2)
368 ! dphi/dsig
369 normxx = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx + dphi_dtrsig
370 normyy = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy + dphi_dtrsig
371 normzz = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz + dphi_dtrsig
372 normxy = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
373c
374 ! 2 - Computation of DPHI_DLAMBDA
375 !---------------------------------------------------------
376c
377 ! a) Derivative with respect stress increments tensor DSIG
378 ! --------------------------------------------------------
379 trdfds = normxx + normyy + normzz
380 dfdsig2 = normxx * (a11*normxx + a12*normyy)
381 . + normyy * (a11*normyy + a12*normxx)
382 . + normxy * normxy * g
383c
384 ! b) Derivatives with respect to plastic strain P
385 ! ------------------------------------------------
386c
387 ! i) Derivative of the yield stress with respect to plastic strain dYLD / dPLA
388 ! ----------------------------------------------------------------------------
389 hardp(i) = hard + qvoce*bvoce*exp(-bvoce*pla(i))
390 IF (igurson == 2) THEN
391 IF (pla_nl(i) - pla(i) < zero) THEN
392 hardp(i) = hardp(i) + hkhi
393 ENDIF
394 ENDIF
395 dyld_dpla = hardp(i)*frate(i)*ftherm(i)
396c
397 ! ii) Derivative of dPLA with respect to DLAM
398 ! -------------------------------------------
399 sdv_dfdsig = sxx(i) * normxx
400 . + syy(i) * normyy
401 . + szz(i) * normzz
402 . + sxy(i) * normxy
403 sig_dfdsig = signxx(i) * normxx
404 . + signyy(i) * normyy
405 . + signxy(i) * normxy
406 dpla_dlam(i) = sig_dfdsig / ((one - ft(i))*yld(i))
407c
408 ! c) Derivatives with respect to the temperature TEMP
409 ! ---------------------------------------------------
410 IF (jthe == 0 .and. cp > zero) THEN
411 ! i) Derivative of the yield stress with respect to temperature dYLD / dTEMP
412 ! ---------------------------------------------------------------------------
413 dyld_dtemp = -fhard(i)*frate(i)*mtemp
414 ! ii) Derivative of the temperature TEMP with respect to DLAM
415 ! -----------------------------------------------------------
416 dtemp_dlam = weitemp(i)*(eta/(rho(i)*cp))*sig_dfdsig
417 ELSE
418 dyld_dtemp = zero
419 dtemp_dlam = zero
420 ENDIF
421c
422 ! d) Derivative with respect to the yield stress
423 ! ----------------------------------------------
424 sigdr2(i) = sigdr(i)**2
425 dphi_dyld = -two*sigdr2(i)/yld(i)**3-dphi_dtrsig*trsig(i)/yld(i)
426c
427 ! e) Derivatives with respect to the damage variables: FG,FN,FSH
428 ! --------------------------------------------------------------
429c
430 ! i) Derivative of PHI with respect to the damage FS
431 fcosh = cosh(q2*etat(i)*trsig(i)/(yld(i)*two))
432 dphi_dfs = two*q1*fcosh - two*q1*q1*fs(i)
433c
434 ! ii) Derivative of FS with respect to FT
435 ! -----------------------------------------------------------
436 IF (ft(i) >= fc) THEN
437 dfs_dft = ((one/q1)-fc)*(fr-fc)
438 ELSE
439 dfs_dft = one
440 ENDIF
441c
442 ! iii) Derivative of FN with respect to LAM
443 ! -----------------------------------------------------------
444 IF ((pla(i)>=ed).AND.(ft(i)<fr)) THEN
445 ! Case for positive stress triaxiality
446 IF (triax(i)>=zero) THEN
447 dfn_dlam = an*dpla_dlam(i)
448 ! Case for negative stress triaxiality
449 ELSEIF ((triax(i)<zero).AND.(triax(i)>=-third)) THEN
450 dfn_dlam = an*max(one + three*triax(i),zero)*dpla_dlam(i)
451 ! Other cases
452 ELSE
453 dfn_dlam = zero
454 ENDIF
455 ELSE
456 dfn_dlam = zero
457 ENDIF
458c
459 ! iv) Derivative of FSH with respect to LAM
460 ! -----------------------------------------------------------
461 IF ((sigdr(i) > zero).AND.(ft(i)>zero).AND.(ft(i)<fr)) THEN
462 sigvm = sqrt(max(em20, three*(j2(i)/(normsig**2))))
463 omega = one - ((twenty7 *(j3(i)/(normsig**3)))/(two*sigvm**3))**2
464 omega = max(omega,zero)
465 omega = min(omega,one)
466 dfsh_dlam = kw*omega*ft(i)*sdv_dfdsig/sigdr(i)
467 ELSE
468 dfsh_dlam = zero
469 ENDIF
470c
471 ! v) Derivative of FG with respect to LAM
472 ! -----------------------------------------------------------
473 IF ((ft(i)>zero).AND.(ft(i)<fr).AND.(trdfds>zero)) THEN
474 dfg_dlam = (one-ft(i))*trdfds
475 ELSE
476 dfg_dlam = zero
477 ENDIF
478c
479 ! vi) Derivative of FT with respect to LAM
480 ! -----------------------------------------------------------
481 dft_dlam = dfn_dlam + dfg_dlam + dfsh_dlam
482c
483 ! e) Derivative of PHI with respect to DLAM
484 dphi_dlam(i) = - dfdsig2 + (dphi_dyld*dyld_dpla*dpla_dlam(i))
485 . + dphi_dfs*dfs_dft*dft_dlam
486 IF (jthe == 0 .and. cp > zero) THEN
487 dphi_dlam(i) = dphi_dlam(i) + dphi_dyld*dyld_dtemp*dtemp_dlam
488 ENDIF
489 dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20) ,dphi_dlam(i))
490c
491 ! 3 - Computation of plastic multiplier and variables update
492 !----------------------------------------------------------
493c
494 ! Computation of the plastic multiplier
495 dlam = -phi(i)/dphi_dlam(i)
496c
497 ! Plastic strains tensor update
498 dpxx(i) = dlam * normxx
499 dpyy(i) = dlam * normyy
500 dpzz(i) = dlam * normzz
501 dpxy(i) = dlam * normxy
502 trdep(i) = dpxx(i) + dpyy(i) + dpzz(i)
503c
504 ! Cumulated plastic strain and strain rate update
505 ddep = (dlam/((one - ft(i))*yld(i)))*sig_dfdsig
506 dpla(i) = max(zero, dpla(i) + ddep)
507 pla(i) = pla(i) + ddep
508c
509 ! Damage variables update
510 ! Growth damage
511 IF ((ft(i)>zero).AND.(ft(i)<fr).AND.(trdep(i)>zero)) THEN
512 fg(i) = fg(i) + (one-ft(i))*trdep(i)
513 ENDIF
514 fg(i) = max(fg(i),zero)
515 ! Nucleation damage
516 IF ((pla(i) >= ed).AND.(ft(i)<fr)) THEN
517 ! Case for positive stress triaxiality
518 IF (triax(i)>=zero) THEN
519 fn(i) = fn(i) + an*ddep
520 ! Case for negative stress triaxiality
521 ELSEIF ((triax(i)<zero).AND.(triax(i)>=-third)) THEN
522 fn(i) = fn(i) + an*max(one + three*triax(i),zero)*ddep
523 ENDIF
524 ENDIF
525 fn(i) = max(fn(i),zero)
526 ! Shear damage
527 IF ((sigdr(i) > zero).AND.(ft(i)>zero).AND.(ft(i)<fr)) THEN
528 sigvm = sqrt(max(em20, three*(j2(i)/(normsig**2))))
529 omega = one - ((twenty7 *(j3(i)/(normsig**3)))/(two*sigvm**3))**2
530 omega = max(zero,omega)
531 omega = min(one,omega)
532 sdpla = sxx(i)*dpxx(i) + syy(i)*dpyy(i) + szz(i)*dpzz(i)
533 . + sxy(i)*dpxy(i)
534 fsh(i) = fsh(i) + kw*omega*ft(i)*(sdpla/sigdr(i))
535 ENDIF
536 fsh(i) = max(fsh(i),zero)
537 ! Total damage
538 ft(i) = f0 + fg(i) + fn(i) + fsh(i)
539 ft(i) = min(ft(i),fr)
540 ! Effective damage
541 IF (ft(i) < fc)THEN
542 fs(i) = ft(i)
543 ELSE
544 fs(i) = fc + (one/q1 - fc) * (ft(i)-fc)/(fr-fc)
545 ENDIF
546 fs(i) = min(fs(i),one/q1)
547c
548 ! Temperature update
549 IF (jthe == 0 .AND. cp > zero) THEN
550 dtemp = weitemp(i)*yld(i)*(one-ft(i))*ddep*eta/(rho(i)*cp)
551 temp(i) = temp(i) + dtemp
552 ftherm(i) = one - mtemp*(temp(i) - tref)
553 ENDIF
554c
555 ! Hardening law update
556 fhard(i) = yld0 + hard * pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
557 IF (igurson == 2) THEN
558 IF (pla_nl(i) - pla(i) < zero) THEN
559 fhard(i) = fhard(i) - hkhi*(pla_nl(i) - pla(i))
560 ENDIF
561 ENDIF
562c
563 ! Yield stress update
564 yld(i) = fhard(i) * frate(i) * ftherm(i)
565 yld(i) = max(yld(i), em10)
566c
567 ! Elasto-plastic stresses update
568 signxx(i) = signxx(i) - (a11*dpxx(i) + a12*dpyy(i))
569 signyy(i) = signyy(i) - (a11*dpyy(i) + a12*dpxx(i))
570 signxy(i) = signxy(i) - dpxy(i)*g
571 trsig(i) = signxx(i) + signyy(i)
572 sigm(i) = -trsig(i) * third
573 sxx(i) = signxx(i) + sigm(i)
574 syy(i) = signyy(i) + sigm(i)
575 szz(i) = sigm(i)
576 sxy(i) = signxy(i)
577c
578 ! Drucker equivalent stress update
579 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) + sxy(i)**2
580 j3(i) = sxx(i) * syy(i) * szz(i) - szz(i) * sxy(i)**2
581 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
582 sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
583 ! Computation of the stress triaxiality and the etaT factor
584 triax(i) = (trsig(i)*third)/sigdr(i)
585 IF (trsig(i)<zero) THEN
586 etat(i) = zero
587 ELSE
588 etat(i) = one
589 ENDIF
590c
591 ! Yield function value update
592 sigdr2(i) = sigdr(i)**2
593 yld2i(i) = one / yld(i)**2
594 fcosh = cosh(q2*etat(i)*trsig(i)/(yld(i)*two))
595 fdam(i) = two*q1*fs(i)*fcosh - (q1*fs(i))**2
596 phi(i) = sigdr2(i) * yld2i(i) - one + fdam(i)
597c
598 ! Transverse strain update
599 dezz(i) = dezz(i) + nnu*dpxx(i) + nnu*dpyy(i) + dpzz(i)
600c
601 ENDDO
602 ! End of the loop over yielding elements
603 ENDDO
604 ! End of the loop over the iterations -
605 !===================================================================
606 ! - END OF PLASTIC CORRECTION WITH CUTTING PLANE ITERATIVE METHOD
607 !===================================================================
608c
609 ! Storing new values
610 DO i=1,nel
611 ! Standard outputs
612 dmg(i,1) = ft(i)/fr ! Normalized total damage
613 dmg(i,2) = fg(i) ! Growth damage
614 dmg(i,3) = fn(i) ! Nucleation damage
615 dmg(i,4) = fsh(i) ! Shear damage
616 dmg(i,5) = min(ft(i),fr) ! Total damage
617 dmg(i,6) = min(fs(i),one/q1) ! Effective damage
618 seq(i) = sigdr(i) ! Equivalent stress
619 ! If element is broken
620 IF (ft(i) >= fr) THEN
621 IF (off(i) == one) off(i) = four_over_5
622 dpla(i) = zero
623 signxx(i) = zero
624 signyy(i) = zero
625 signxy(i) = zero
626 signyz(i) = zero
627 signzx(i) = zero
628 seq(i) = zero
629 ENDIF
630 ! Plastic strain-rate (filtered)
631 dpdt = dpla(i) / max(em20,timestep)
632 epsd(i) = afiltr * dpdt + (one - afiltr) * epsd(i)
633 ! Coefficient for hourglass
634 IF (dpla(i) > zero) THEN
635 et(i) = hardp(i)*frate(i) / (hardp(i)*frate(i) + young)
636 ELSE
637 et(i) = one
638 ENDIF
639 ! Computation of the sound speed
640 soundsp(i) = sqrt((a11)/rho(i))
641 ! Storing the yield stress
642 sigy(i) = yld(i)
643 ! Computation of the thickness variation
644 thk(i) = thk(i) + dezz(i)*thkly(i)*off(i)
645 ENDDO
646c
647 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
mat104c_ldam_newton
subroutine mat104c_ldam_newton(nel, ngl, nuparam, nuvar, time, timestep, uparam, uvar, jthe, off, gs, rho, pla, dpla, epsd, soundsp, depsxx, depsyy, depsxy, depsyz, depszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, thkly, thk, sigy, et, tempel, dpla_nl, dmg, temp, seq, pla_nl, l_planl, plap_nl, l_epsdnl)
Definition mat104c_ldam_newton.F:37