OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat104_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 mat104_ldam_nice (nel, ngl, nuparam, nuvar, volume, fheat, 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, jlag, seq, pla_nl, l_planl, plap_nl, l_epsdnl)

Function/Subroutine Documentation

◆ mat104_ldam_nice()

subroutine mat104_ldam_nice ( integer nel,
integer, dimension(nel), intent(in) ngl,
integer nuparam,
integer nuvar,
intent(in) volume,
intent(inout) fheat,
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,
integer, intent(in) jlag,
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 mat104_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 ,INTENT(IN) :: JLAG
50 INTEGER ,DIMENSION(NEL), INTENT(IN) :: NGL
51 INTEGER ,INTENT(IN) :: L_PLANL,L_EPSDNL
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
65 my_real ,DIMENSION(NEL), INTENT(INOUT) :: fheat
66c
67 my_real ,DIMENSION(NEL), INTENT(OUT) ::
68 . soundsp,sigy,et,
69 . signxx,signyy,signzz,signxy,signyz,signzx
70c
71 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
72 . pla,dpla,epsd,off,temp,seq
73 my_real ,DIMENSION(NEL,6), INTENT(INOUT) ::
74 . dmg
75 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) ::
76 . uvar
77 !=======================================================================
78 ! Local Variables
79 !=======================================================================
80 INTEGER I,K,II,IGURSON,NINDX,INDEX(NEL),NICE
81c
82 my_real ::
83 . young,bulk,lam,g,g2,nu,cdr,kdr,hard,yld0,qvoce,bvoce,jcc,
84 . epsp0,mtemp,tini,tref,eta,cp,dpis,dpad,asrate,afiltr,hkhi,
85 . q1,q2,ed,an,epn,kw,fr,fc,f0,normsig,dti
86 my_real ::
87 . h,ldav,trdep,trdfds,sigvm,sigdr2,yld2i,omega,
88 . dtemp,fcosh,fsinh,dpdt,dlam,ddep
89 my_real ::
90 . dpxx,dpyy,dpzz,dpxy,dpyz,dpzx,dsdrdj2,dsdrdj3,
91 . dj3dsxx,dj3dsyy,dj3dszz,dj3dsxy,dj3dsyz,dj3dszx,
92 . dj2dsxx,dj2dsyy,dj2dszz,dj2dsxy,dj2dsyz,dj2dszx,
93 . dfdsxx,dfdsyy,dfdszz,dfdsxy,dfdsyz,dfdszx,
94 . normxx,normyy,normzz,normxy,normyz,normzx,
95 . denom,dphi,sdpla,dphi_dtrsig,sdv_dfdsig,sig_dfdsig,dfdsig2,
96 . dphi_dsig,dphi_dyld,dphi_dfdr,dphi_dfs,dfs_dft,
97 . dfn_dlam,dfsh_dlam,dfg_dlam,dft_dlam,dfdpla,dyld_dpla_nl,
98 . dfn,dfsh,dfg,dft,dyld_dpla,dyld_dtemp,dtemp_dlam
99c
100 my_real, DIMENSION(NEL) ::
101 . dsigxx,dsigyy,dsigzz,dsigxy,dsigyz,dsigzx,trsig,
102 . sxx,syy,szz,sxy,syz,szx,sigm,j2,j3,sigdr,yld,weitemp,
103 . hardp,fhard,frate,ftherm,fdr,depxx,depyy,depzz,depxy,depyz,depzx,
104 . fdam,phi0,phi,ft,fs,fg,fn,fsh,dpla_dlam,triax,etat
105c
106 !=======================================================================
107 ! DRUCKER - VOCE - JOHNSON-COOK MATERIAL LAW WITH GURSON DAMAGE
108 ! USING LOCAL DAMAGE OR NON-LOCAL MICROMORPHIC METHOD
109 !=======================================================================
110 !UVAR(1) PHI YIELD FUNCTION VALUE
111 !=======================================================================
112c
113 !=======================================================================
114 ! - INITIALISATION OF COMPUTATION ON TIME STEP
115 !=======================================================================
116 ! Recovering model parameters
117 ! Elastic parameters
118 young = uparam(1) ! Young modulus
119 bulk = uparam(2) ! Bulk modulus
120 g = uparam(3) ! Shear modulus
121 g2 = uparam(4) ! 2*Shear modulus
122 lam = uparam(5) ! Lambda Hooke parameter
123 nu = uparam(6) ! Poisson ration
124 ! Plastic criterion and hardening parameters [Drucker, 1948]
125 nice = nint(uparam(11)) ! Choice of the Nice method
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)
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) ! Void volume fracture at failure
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 ! User inputs
169 phi0(i) = uvar(i,1) ! Previous yield function value
170 ! Damage variables
171 fg(i) = dmg(i,2) ! Growth damage
172 fn(i) = dmg(i,3) ! Nucleation damage
173 fsh(i) = dmg(i,4) ! Shear damage
174 ft(i) = dmg(i,5) ! Total damage
175 fs(i) = dmg(i,6) ! Effective damage
176 ! Standard inputs
177 dpla(i) = zero ! Initialization of the plastic strain increment
178 et(i) = one ! Initialization of hourglass coefficient
179 hardp(i) = zero ! Initialization of hardening modulus
180 ENDDO
181c
182 ! Initialization of damage, temperature and self-heating weight factor
183 IF (time == zero) THEN
184 IF (jthe == 0) temp(1:nel) = tini
185 IF (isigi == 0) THEN
186 dmg(1:nel,5) = f0
187 ft(1:nel) = f0
188 dmg(1:nel,1) = f0/fr
189 IF (f0<fc) THEN
190 dmg(1:nel,6) = f0
191 ELSE
192 dmg(1:nel,6) = fc + (one/q1-fc)*(f0-fc)/(fr-fc)
193 ENDIF
194 fs(1:nel) = dmg(1:nel,6)
195 ENDIF
196 ENDIF
197 IF (cp > zero .OR. 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 ENDIF
211c
212 ! Computation of the initial yield stress
213 DO i=1,nel
214 ! a) - Hardening law
215 fhard(i) = yld0 + hard*pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
216 IF (igurson == 2) THEN
217 IF (pla_nl(i) - pla(i) < zero) THEN
218 fhard(i) = fhard(i) - hkhi*(pla_nl(i) - pla(i))
219 ENDIF
220 ENDIF
221 ! b) - Correction factor for strain-rate dependence (Johnson-Cook)
222 frate(i) = one
223 IF (epsd(i) > epsp0) frate(i) = frate(i) + jcc*log(epsd(i)/epsp0)
224 ! c) - Correction factor for thermal effects
225 ftherm(i) = one - mtemp*(temp(i) - tref)
226 ! d) - Computation of the yield function
227 yld(i) = fhard(i)*frate(i)*ftherm(i)
228 ! e) - Checking values
229 yld(i) = max(em10, yld(i))
230 ENDDO
231c
232 !========================================================================
233 ! - COMPUTATION OF TRIAL VALUES
234 !========================================================================
235 DO i=1,nel
236c
237 ! Computation of the trial stress tensor
238 ldav = (depsxx(i) + depsyy(i) + depszz(i)) * lam
239 signxx(i) = sigoxx(i) + (depsxx(i)*g2 + ldav)
240 signyy(i) = sigoyy(i) + (depsyy(i)*g2 + ldav)
241 signzz(i) = sigozz(i) + (depszz(i)*g2 + ldav)
242 signxy(i) = sigoxy(i) + depsxy(i)*g
243 signyz(i) = sigoyz(i) + depsyz(i)*g
244 signzx(i) = sigozx(i) + depszx(i)*g
245 ! Computation of the trace of the trial stress tensor
246 trsig(i) = signxx(i) + signyy(i) + signzz(i)
247 sigm(i) = -trsig(i) * third
248 ! Computation of the deviatoric trial stress tensor
249 sxx(i) = signxx(i) + sigm(i)
250 syy(i) = signyy(i) + sigm(i)
251 szz(i) = signzz(i) + sigm(i)
252 sxy(i) = signxy(i)
253 syz(i) = signyz(i)
254 szx(i) = signzx(i)
255 ! second deviatoric invariant
256 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 )
257 . + sxy(i)**2 + syz(i)**2 + szx(i)**2
258 ! Third deviatoric invariant
259 j3(i) = sxx(i) * syy(i) * szz(i)
260 . + sxy(i) * syz(i) * szx(i) * two
261 . - sxx(i) * syz(i)**2
262 . - syy(i) * szx(i)**2
263 . - szz(i) * sxy(i)**2
264 ! Drucker equivalent stress
265 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
266 ! Checking equivalent stress values
267 IF (fdr(i) > zero) THEN
268 sigdr(i) = kdr * exp((one/six)*log(fdr(i))) ! fdr(i)**(1/6)
269 ELSE
270 sigdr(i) = zero
271 ENDIF
272 ! Computation of the stress triaxiality and the etaT factor
273 IF (sigdr(i)>zero) THEN
274 triax(i) = (trsig(i)*third)/sigdr(i)
275 ELSE
276 triax(i) = zero
277 ENDIF
278 IF (trsig(i)<zero) THEN
279 etat(i) = zero
280 ELSE
281 etat(i) = one
282 ENDIF
283 ENDDO
284c
285 !========================================================================
286 ! - COMPUTATION OF YIELD FONCTION
287 !========================================================================
288 DO i=1,nel
289 fdam(i) = two*q1*fs(i)*cosh(q2*etat(i)*trsig(i)/yld(i)/two) - (q1*fs(i))**2
290 phi(i) = (sigdr(i) / yld(i))**2 - one + fdam(i)
291 ENDDO
292c
293 ! Checking plastic behavior for all elements
294 nindx = 0
295 DO i=1,nel
296 IF ((phi(i) >= zero).AND.(off(i) == one).AND.(ft(i) < fr)) THEN
297 nindx=nindx+1
298 index(nindx)=i
299 ENDIF
300 ENDDO
301c
302 !====================================================================
303 ! - PLASTIC CORRECTION WITH NICE - EXPLICIT METHOD
304 !====================================================================
305#include "vectorize.inc"
306 DO ii=1,nindx
307c
308 ! Number of the element with plastic behaviour
309 i = index(ii)
310c
311 ! Computation of the trial stress increment
312 ldav = (depsxx(i) + depsyy(i) + depszz(i)) * lam
313 dsigxx(i) = depsxx(i)*g2 + ldav
314 dsigyy(i) = depsyy(i)*g2 + ldav
315 dsigzz(i) = depszz(i)*g2 + ldav
316 dsigxy(i) = depsxy(i)*g
317 dsigyz(i) = depsyz(i)*g
318 dsigzx(i) = depszx(i)*g
319 dlam = zero
320c
321 ! Computation of Drucker equivalent stress of the previous stress tensor
322 trsig(i) = sigoxx(i) + sigoyy(i) + sigozz(i)
323 sigm(i) = -trsig(i) * third
324 sxx(i) = sigoxx(i) + sigm(i)
325 syy(i) = sigoyy(i) + sigm(i)
326 szz(i) = sigozz(i) + sigm(i)
327 sxy(i) = sigoxy(i)
328 syz(i) = sigoyz(i)
329 szx(i) = sigozx(i)
330 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 )
331 . + sxy(i)**2 + syz(i)**2 + szx(i)**2
332 j3(i) = sxx(i)*syy(i)*szz(i) + sxy(i)*syz(i)*szx(i) * two
333 . - sxx(i)*syz(i)**2 - syy(i)*szx(i)**2 - szz(i)*sxy(i)**2
334 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
335 sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
336 ! Computation of the stress triaxiality and the etaT factor
337 triax(i) = (trsig(i)*third)/sigdr(i)
338 IF (trsig(i)<zero) THEN
339 etat(i) = zero
340 ELSE
341 etat(i) = one
342 ENDIF
343c
344 ! Note : in this part, the purpose is to compute for each iteration
345 ! a plastic multiplier allowing to update internal variables to satisfy
346 ! the consistency condition.
347 ! Its expression is : DLAMBDA = - (PHI_OLD + DPHI) / DPHI_DLAMBDA
348 ! -> PHI_OLD : old value of yield function (known)
349 ! -> DPHI : yield function prediction (to compute)
350 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
351 ! into account of internal variables kinetic :
352 ! plasticity, temperature and damage (to compute)
353c
354 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
355 !-------------------------------------------------------------
356c
357 ! Derivative with respect to the equivalent stress and trace
358 yld2i = one/(yld(i)**2)
359 dphi_dsig = two*sigdr(i)*yld2i
360 fsinh = sinh(q2*etat(i)*trsig(i)/(yld(i)*two))
361 dphi_dtrsig = q1*q2*etat(i)*fs(i)*fsinh/yld(i)
362c
363 ! Computation of the Eulerian norm of the stress tensor
364 normsig = sqrt(sigoxx(i)*sigoxx(i)
365 . + sigoyy(i)*sigoyy(i)
366 . + sigozz(i)*sigozz(i)
367 . + two*sigoxy(i)*sigoxy(i)
368 . + two*sigoyz(i)*sigoyz(i)
369 . + two*sigozx(i)*sigozx(i))
370 normsig = max(normsig,one)
371c
372 ! Derivative with respect to Fdr
373 fdr(i) = (j2(i)/(normsig**2))**3 - cdr*((j3(i)/(normsig**3))**2)
374 dphi_dfdr = dphi_dsig*kdr*(one/six)*exp(-(five/six)*log(fdr(i)))
375 dsdrdj2 = dphi_dfdr*three*(j2(i)/(normsig**2))**2
376 dsdrdj3 = -dphi_dfdr*two*cdr*(j3(i)/(normsig**3))
377 ! dJ3/dSig
378 dj3dsxx = two_third*(syy(i)*szz(i)-syz(i)**2)/(normsig**2)
379 . - third*(sxx(i)*szz(i)-szx(i)**2)/(normsig**2)
380 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
381 dj3dsyy = - third*(syy(i)*szz(i)-syz(i)**2)/(normsig**2)
382 . + two_third*(sxx(i)*szz(i)-szx(i)**2)/(normsig**2)
383 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
384 dj3dszz = - third*(syy(i)*szz(i)-syz(i)**2)/(normsig**2)
385 . - third*(sxx(i)*szz(i)-szx(i)**2)/(normsig**2)
386 . + two_third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
387 dj3dsxy = two*(sxx(i)*sxy(i) + sxy(i)*syy(i) + szx(i)*syz(i))/(normsig**2)
388 dj3dsyz = two*(sxy(i)*szx(i) + syy(i)*syz(i) + syz(i)*szz(i))/(normsig**2)
389 dj3dszx = two*(sxx(i)*szx(i) + sxy(i)*syz(i) + szx(i)*szz(i))/(normsig**2)
390 ! dPhi/dSig
391 normxx = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx + dphi_dtrsig
392 normyy = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy + dphi_dtrsig
393 normzz = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz + dphi_dtrsig
394 normxy = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
395 normyz = two*dsdrdj2*syz(i)/normsig + dsdrdj3*dj3dsyz
396 normzx = two*dsdrdj2*szx(i)/normsig + dsdrdj3*dj3dszx
397c
398 ! Restoring previous value of the yield function
399 phi(i) = phi0(i)
400c
401 ! Computation of yield surface trial increment DPHI
402 dphi = normxx * dsigxx(i)
403 . + normyy * dsigyy(i)
404 . + normzz * dsigzz(i)
405 . + normxy * dsigxy(i)
406 . + normyz * dsigyz(i)
407 . + normzx * dsigzx(i)
408c
409 ! 2 - Computation of DPHI_DLAMBDA
410 !--------------------------------------------------------
411 ! a) Derivative with respect stress increments tensor DSIG
412 ! --------------------------------------------------------
413 trdfds = normxx + normyy + normzz
414 dfdsig2 = normxx *(normxx*g2 + lam*trdfds)
415 . + normyy *(normyy*g2 + lam*trdfds)
416 . + normzz *(normzz*g2 + lam*trdfds)
417 . + normxy *normxy*g
418 . + normyz *normyz*g
419 . + normzx *normzx*g
420c
421 ! b) Derivatives with respect to plastic strain P
422 ! ------------------------------------------------
423c
424 ! i) Derivative of the yield stress with respect to plastic strain dYLD / dPLA
425 ! ----------------------------------------------------------------------------
426 hardp(i) = hard + qvoce*bvoce*exp(-bvoce*pla(i))
427 IF (igurson == 2) THEN
428 IF (pla_nl(i) - pla(i) < zero) THEN
429 hardp(i) = hardp(i) + hkhi
430 ENDIF
431 ENDIF
432 dyld_dpla = hardp(i)*frate(i)*ftherm(i)
433c
434 IF (igurson == 2) THEN
435 IF (pla_nl(i) - pla(i) < zero) THEN
436 dyld_dpla_nl = -hkhi*frate(i)*ftherm(i)
437 ELSE
438 dyld_dpla_nl = zero
439 ENDIF
440 ENDIF
441c
442 ! ii) Derivative of dPLA with respect to DLAM
443 ! -------------------------------------------
444 sdv_dfdsig = sxx(i) * normxx
445 . + syy(i) * normyy
446 . + szz(i) * normzz
447 . + sxy(i) * normxy
448 . + syz(i) * normyz
449 . + szx(i) * normzx
450 sig_dfdsig = sigoxx(i) * normxx
451 . + sigoyy(i) * normyy
452 . + sigozz(i) * normzz
453 . + sigoxy(i) * normxy
454 . + sigoyz(i) * normyz
455 . + sigozx(i) * normzx
456 dpla_dlam(i) = sig_dfdsig / ((one - ft(i))*yld(i))
457c
458 ! c) Derivatives with respect to the temperature TEMP
459 ! ---------------------------------------------------
460 IF (jthe == 0 .and. cp > zero) THEN
461 ! i) Derivative of the yield stress with respect to temperature dYLD / dTEMP
462 ! ---------------------------------------------------------------------------
463 dyld_dtemp = -fhard(i)*frate(i)*mtemp
464 ! ii) Derivative of the temperature TEMP with respect to DLAM
465 ! -----------------------------------------------------------
466 dtemp_dlam = weitemp(i)*(eta/(rho0(i)*cp))*sig_dfdsig
467 ELSE
468 dyld_dtemp = zero
469 dtemp_dlam = zero
470 ENDIF
471c
472 ! d) Derivative with respect to the yield stress
473 ! ----------------------------------------------
474 sigdr2 = sigdr(i)**2
475 dphi_dyld = -two*sigdr2/yld(i)**3-dphi_dtrsig*trsig(i)/yld(i)
476c
477 ! e) Derivatives with respect to the damage variables: FG,FN,FSH
478 ! --------------------------------------------------------------
479c
480 ! i) Derivative of PHI with respect to the damage FS
481 fcosh = cosh(q2*etat(i)*trsig(i)/(yld(i)*two))
482 dphi_dfs = two*q1*fcosh - two*q1*q1*fs(i)
483c
484 ! ii) Derivative of FS with respect to FT
485 ! -----------------------------------------------------------
486 IF (ft(i) >= fc) THEN
487 dfs_dft = ((one/q1)-fc)*(fr-fc)
488 ELSE
489 dfs_dft = one
490 ENDIF
491c
492 ! iii) Derivative of FN with respect to LAM
493 ! -----------------------------------------------------------
494 IF ((pla(i)>=ed).AND.(ft(i)<fr)) THEN
495 ! Case for positive stress triaxiality
496 IF (triax(i)>=zero) THEN
497 dfn_dlam = an*dpla_dlam(i)
498 ! Case for negative stress triaxiality
499 ELSEIF ((triax(i)<zero).AND.(triax(i)>=-third)) THEN
500 dfn_dlam = an*max(one + three*triax(i),zero)*dpla_dlam(i)
501 ! Other cases
502 ELSE
503 dfn_dlam = zero
504 ENDIF
505 ELSE
506 dfn_dlam = zero
507 ENDIF
508c
509 ! iv) Derivative of FSH with respect to LAM
510 ! -----------------------------------------------------------
511 IF ((sigdr(i) > zero).AND.(ft(i)>zero).AND.(ft(i)<fr)) THEN
512 sigvm = sqrt(max(em20, three*(j2(i)/(normsig**2))))
513 omega = one - ((twenty7 *(j3(i)/(normsig**3)))/(two*sigvm**3))**2
514 omega = max(zero,omega)
515 omega = min(one,omega)
516 dfsh_dlam = kw*omega*ft(i)*sdv_dfdsig/sigdr(i)
517 ELSE
518 dfsh_dlam = zero
519 ENDIF
520c
521 ! v) Derivative of FG with respect to LAM
522 ! -----------------------------------------------------------
523 IF ((ft(i)>zero).AND.(ft(i)<fr).AND.(trdfds>zero)) THEN
524 dfg_dlam = (one-ft(i))*trdfds
525 ELSE
526 dfg_dlam = zero
527 ENDIF
528c
529 ! vi) Derivative of FT with respect to LAM
530 ! -----------------------------------------------------------
531 dft_dlam = dfn_dlam + dfg_dlam + dfsh_dlam
532c
533 ! e) Derivative of PHI with respect to DLAM ( = -DENOM)
534 denom = dfdsig2 - (dphi_dyld * dyld_dpla * dpla_dlam(i))
535 . - dphi_dfs * dfs_dft * dft_dlam
536 IF (jthe == 0 .and. cp > zero) THEN
537 denom = denom - dphi_dyld * dyld_dtemp * dtemp_dlam
538 ENDIF
539 denom = sign(max(abs(denom),em20) ,denom)
540c
541 ! 3 - Computation of plastic multiplier and variables update
542 !----------------------------------------------------------
543c
544 ! Computation of the plastic multiplier
545 IF (igurson == 2) THEN
546 dphi = dphi + dphi_dyld*dyld_dpla_nl*dpla_nl(i)
547 ENDIF
548 dlam = (dphi + phi(i)) / denom
549 dlam = max(dlam, zero)
550c
551 ! Plastic strains tensor update
552 dpxx = dlam * normxx
553 dpyy = dlam * normyy
554 dpzz = dlam * normzz
555 dpxy = dlam * normxy
556 dpyz = dlam * normyz
557 dpzx = dlam * normzx
558 trdep = dpxx + dpyy + dpzz
559c
560 ! Cumulated plastic strain and strain rate update
561 ddep = (dlam/((one - ft(i))*yld(i)))*sig_dfdsig
562 dpla(i) = dpla(i) + ddep
563 pla(i) = pla(i) + ddep
564c
565 ! Damage variables update
566 ! Growth damage
567 IF ((ft(i)>zero).AND.(ft(i)<fr).AND.(trdep>zero)) THEN
568 fg(i) = fg(i) + (one-ft(i))*trdep
569 ENDIF
570 fg(i) = max(fg(i),zero)
571 ! Nucleation damage
572 IF ((pla(i) >= ed).AND.(ft(i)<fr)) THEN
573 ! Case for positive stress triaxiality
574 IF (triax(i)>=zero) THEN
575 fn(i) = fn(i) + an*ddep
576 ! Case for negative stress triaxiality
577 ELSEIF ((triax(i)<zero).AND.(triax(i)>=-third)) THEN
578 fn(i) = fn(i) + an*max(one + three*triax(i),zero)*ddep
579 ENDIF
580 ENDIF
581 fn(i) = max(fn(i),zero)
582 ! Shear damage
583 IF ((sigdr(i) > zero).AND.(ft(i)>zero).AND.(ft(i)<fr)) THEN
584 sigvm = sqrt(max(em20, three*(j2(i)/(normsig**2))))
585 omega = one - ((twenty7 *(j3(i)/(normsig**3)))/(two*sigvm**3))**2
586 omega = max(zero,omega)
587 omega = min(one,omega)
588 sdpla = sxx(i)*dpxx + syy(i)*dpyy + szz(i)*dpzz
589 . + sxy(i)*dpxy + syz(i)*dpyz + szx(i)*dpzx
590 fsh(i) = fsh(i) + kw*omega*ft(i)*(sdpla/sigdr(i))
591 ENDIF
592 fsh(i) = max(fsh(i),zero)
593 ! Total damage
594 ft(i) = f0 + fg(i) + fn(i) + fsh(i)
595 ft(i) = min(ft(i),fr)
596 ! Effective damage
597 IF (ft(i) < fc)THEN
598 fs(i) = ft(i)
599 ELSE
600 fs(i) = fc + (one/q1 - fc) * (ft(i)-fc)/(fr-fc)
601 ENDIF
602 fs(i) = min(fs(i),one/q1)
603c
604 ! Temperature update
605 IF (jthe == 0 .AND. cp > zero ) THEN
606 dtemp = weitemp(i)*yld(i)*(one-ft(i))*ddep*eta/(rho0(i)*cp)
607 temp(i) = temp(i) + dtemp
608 ENDIF
609 ftherm(i) = one - mtemp*(temp(i) - tref)
610c
611 ! hardening law update
612 fhard(i) = yld0 + hard * pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
613 IF (igurson == 2) THEN
614 IF (pla_nl(i) - pla(i) < zero) THEN
615 fhard(i) = fhard(i) - hkhi*(pla_nl(i) - pla(i))
616 ENDIF
617 ENDIF
618c
619 ! Yield stress update
620 yld(i) = fhard(i) * frate(i) * ftherm(i)
621 yld(i) = max(yld(i), em10)
622c
623 ! Elasto-plastic stresses update
624 ldav = trdep * lam
625 signxx(i) = sigoxx(i) + dsigxx(i) - (dpxx*g2 + ldav)
626 signyy(i) = sigoyy(i) + dsigyy(i) - (dpyy*g2 + ldav)
627 signzz(i) = sigozz(i) + dsigzz(i) - (dpzz*g2 + ldav)
628 signxy(i) = sigoxy(i) + dsigxy(i) - dpxy*g
629 signyz(i) = sigoyz(i) + dsigyz(i) - dpyz*g
630 signzx(i) = sigozx(i) + dsigzx(i) - dpzx*g
631 trsig(i) = signxx(i) + signyy(i) + signzz(i)
632 sigm(i) = -trsig(i) * third
633 sxx(i) = signxx(i) + sigm(i)
634 syy(i) = signyy(i) + sigm(i)
635 szz(i) = signzz(i) + sigm(i)
636 sxy(i) = signxy(i)
637 syz(i) = signyz(i)
638 szx(i) = signzx(i)
639c
640 ! Drucker equivalent stress update
641 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 )
642 . + sxy(i)**2 + syz(i)**2 + szx(i)**2
643 j3(i) = sxx(i) * syy(i) * szz(i) + two * sxy(i) * syz(i) * szx(i)
644 . - sxx(i) * syz(i)**2 - syy(i) * szx(i)**2 - szz(i) * sxy(i)**2
645 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
646 sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
647 ! Computation of the stress triaxiality and the etaT factor
648 triax(i) = (trsig(i)*third)/sigdr(i)
649 IF (trsig(i)<zero) THEN
650 etat(i) = zero
651 ELSE
652 etat(i) = one
653 ENDIF
654c
655 ! Yield function value update
656 sigdr2 = sigdr(i)**2
657 yld2i = one / yld(i)**2
658 fcosh = cosh(q2*etat(i)*trsig(i)/(yld(i)*two))
659 fdam(i) = two*q1*fs(i)*fcosh - (q1*fs(i))**2
660 phi(i) = sigdr2 * yld2i - one + fdam(i)
661 ENDDO
662 !===================================================================
663 ! - END OF PLASTIC CORRECTION WITH NICE - EXPLICIT 1 METHOD
664 !===================================================================
665c
666 ! Storing new values
667 DO i=1,nel
668 ! USR Outputs & coefficient for hourglass
669 IF (dpla(i) > zero) THEN
670 uvar(i,1) = phi(i) ! Yield function value
671 et(i) = hardp(i)*frate(i) / (hardp(i)*frate(i) + young)
672 ELSE
673 uvar(i,1) = zero
674 et(i) = one
675 ENDIF
676 ! Standard outputs
677 dmg(i,1) = ft(i)/fr ! Normalized total damage
678 dmg(i,2) = fg(i) ! Growth damage
679 dmg(i,3) = fn(i) ! Nucleation damage
680 dmg(i,4) = fsh(i) ! Shear damage
681 dmg(i,5) = min(ft(i),fr) ! total damage
682 dmg(i,6) = min(fs(i),one/q1) ! Effective damage
683 seq(i) = sigdr(i) ! Equivalent stress
684 ! If element is broken
685 IF (ft(i) >= fr) THEN
686 IF (off(i) == one) off(i) = four_over_5
687 dpla(i) = zero
688 signxx(i) = zero
689 signyy(i) = zero
690 signzz(i) = zero
691 signxy(i) = zero
692 signyz(i) = zero
693 signzx(i) = zero
694 seq(i) = zero
695 ENDIF
696 ! Plastic strain-rate (filtered)
697 dpdt = dpla(i) / max(em20,timestep)
698 epsd(i) = afiltr * dpdt + (one - afiltr) * epsd(i)
699 ! Computation of the sound speed
700 soundsp(i) = sqrt((bulk + four_over_3*g)/rho0(i))
701 ! Storing the yield stress
702 sigy(i) = yld(i)
703 ENDDO
704C-----------------------------------------------
705C plastic work dissipated as heat - in case of /heat/mat
706C-----------------------------------------------
707 IF (jthe < 0 .AND. jlag /= 0) THEN
708 DO i=1,nel
709 fheat(i) = fheat(i) + eta*(one-ft(i))*weitemp(i)*yld(i)*dpla(i)*volume(i)
710 ENDDO
711 END IF
712!-----------
713 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
second
real function second()
SECOND Using ETIME
Definition second_EXT_ETIME.f:35