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