OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat104_ldam_nice.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| mat104_ldam_nice ../engine/source/materials/mat/mat104/mat104_ldam_nice.F
25!||--- called by ------------------------------------------------------
26!|| sigeps104 ../engine/source/materials/mat/mat104/sigeps104.F
27!||====================================================================
28 SUBROUTINE mat104_ldam_nice(
29 1 NEL ,NGL ,NUPARAM ,NUVAR ,VOLUME ,FHEAT ,
30 2 TIME ,TIMESTEP,UPARAM ,UVAR ,JTHE ,OFF ,
31 3 RHO0 ,RHO ,PLA ,DPLA ,EPSD ,SOUNDSP ,
32 4 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
33 5 SIGOXX ,SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
34 6 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
35 7 SIGY ,ET ,DPLA_NL ,DMG ,TEMP ,JLAG ,
36 8 SEQ ,PLA_NL ,L_PLANL ,PLAP_NL ,L_EPSDNL)
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
714 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
mat104_ldam_nice
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)
Definition mat104_ldam_nice.F:37