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