OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat104_ldam_newton.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 mat104_ldam_newton (nel, ngl, nuparam, nuvar, volume, time, timestep, uparam, uvar, jthe, off, rho0, rho, pla, dpla, epsd, soundsp, depsxx, depsyy, depszz, depsxy, depsyz, depszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigy, et, dpla_nl, dmg, temp, fheat, seq, pla_nl, l_planl, plap_nl, l_epsdnl, jlag)

Function/Subroutine Documentation

◆ mat104_ldam_newton()

subroutine mat104_ldam_newton ( integer nel,
integer, dimension(nel), intent(in) ngl,
integer nuparam,
integer nuvar,
intent(in) volume,
time,
timestep,
intent(in) uparam,
intent(inout) uvar,
integer jthe,
intent(inout) off,
intent(in) rho0,
intent(in) rho,
intent(inout) pla,
intent(inout) dpla,
intent(inout) epsd,
intent(out) soundsp,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depszz,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
intent(in) sigoxx,
intent(in) sigoyy,
intent(in) sigozz,
intent(in) sigoxy,
intent(in) sigoyz,
intent(in) sigozx,
intent(out) signxx,
intent(out) signyy,
intent(out) signzz,
intent(out) signxy,
intent(out) signyz,
intent(out) signzx,
intent(out) sigy,
intent(out) et,
intent(in) dpla_nl,
intent(inout) dmg,
intent(inout) temp,
intent(inout) fheat,
intent(inout) seq,
intent(in) pla_nl,
integer, intent(in) l_planl,
intent(in) plap_nl,
integer, intent(in) l_epsdnl,
integer, intent(in) jlag )

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