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