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