OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat104c_ldam_nice.F File Reference
#include "implicit_f.inc"
#include "com01_c.inc"
#include "vectorize.inc"

Go to the source code of this file.

Functions/Subroutines

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)

Function/Subroutine Documentation

◆ mat104c_ldam_nice()

subroutine mat104c_ldam_nice ( integer nel,
integer, dimension(nel), intent(in) ngl,
integer nuparam,
integer nuvar,
time,
timestep,
intent(in) uparam,
intent(inout) uvar,
integer jthe,
intent(inout) off,
intent(in) gs,
intent(in) rho,
intent(inout) pla,
intent(inout) dpla,
intent(inout) epsd,
intent(out) soundsp,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
intent(in) sigoxx,
intent(in) sigoyy,
intent(in) sigoxy,
intent(in) sigoyz,
intent(in) sigozx,
intent(out) signxx,
intent(out) signyy,
intent(out) signxy,
intent(out) signyz,
intent(out) signzx,
intent(in) thkly,
intent(inout) thk,
intent(out) sigy,
intent(out) et,
intent(in) tempel,
intent(in) dpla_nl,
intent(inout) dmg,
intent(inout) temp,
intent(inout) seq,
intent(in) pla_nl,
integer, intent(in) l_planl,
intent(in) plap_nl,
integer, intent(in) l_epsdnl )

Definition at line 28 of file mat104c_ldam_nice.F.

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