OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat104c_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 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)

Function/Subroutine Documentation

◆ mat104c_ldam_newton()

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