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