OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat104c_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!|| mat104c_nldam_newton ../engine/source/materials/mat/mat104/mat104c_nldam_newton.F
25!||--- called by ------------------------------------------------------
26!|| sigeps104c ../engine/source/materials/mat/mat104/sigeps104c.F
27!||====================================================================
28 SUBROUTINE mat104c_nldam_newton(
29 1 NEL ,NGL ,NUPARAM ,NUVAR ,
30 2 TIME ,TIMESTEP,UPARAM ,UVAR ,JTHE ,OFF ,
31 3 GS ,RHO ,PLA ,DPLA ,EPSD ,SOUNDSP ,
32 4 DEPSXX ,DEPSYY ,DEPSXY ,DEPSYZ ,DEPSZX ,
33 5 SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
34 6 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,THKLY ,
35 7 THK ,SIGY ,ET ,TEMPEL ,DPLA_NL ,DMG ,
36 8 TEMP ,SEQ ,PLA_NL ,PLAP_NL )
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 my_real
51 . TIME,TIMESTEP
52 my_real,DIMENSION(NUPARAM), INTENT(IN) ::
53 . UPARAM
54 my_real,DIMENSION(NEL), INTENT(IN) ::
55 . RHO,TEMPEL,
56 . depsxx,depsyy,depsxy,depsyz,depszx,
57 . sigoxx,sigoyy,sigoxy,sigoyz,sigozx,
58 . gs,thkly,pla_nl,plap_nl,dpla_nl
59c
60 my_real ,DIMENSION(NEL), INTENT(OUT) ::
61 . soundsp,sigy,et,
62 . signxx,signyy,signxy,signyz,signzx
63c
64 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
65 . pla,dpla,epsd,off,thk,temp,seq
66 my_real ,DIMENSION(NEL,6), INTENT(INOUT) ::
67 . dmg
68 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) ::
69 . uvar
70 !=======================================================================
71 ! Local Variables
72 !=======================================================================
73 INTEGER I,II,IGURSON,ITER,NITER,NINDX,INDEX(NEL)
74c
75 my_real ::
76 . young,bulk,lam,g,g2,nu,cdr,kdr,hard,yld0,qvoce,bvoce,jcc,
77 . epsp0,mtemp,tini,tref,eta,cp,dpis,dpad,asrate,afiltr,hkhi,
78 . q1,q2,ed,an,epn,kw,fr,fc,f0,a11,a12,nnu
79 my_real ::
80 . dti,h,ldav,sigvm,omega,
81 . dtemp,fcosh,fsinh,dpdt,dlam,ddep
82 my_real ::
83 . dsdrdj2,dsdrdj3,
84 . dj3dsxx,dj3dsyy,dj3dsxy,dj3dszz,
85 . dj2dsxx,dj2dsyy,dj2dsxy,normsig,
86 . dfdsxx,dfdsyy,dfdsxy,
87 . sdpla,dphi_dtrsig,dfdsig2,sdv_dfdsig,
88 . dphi_dsig,dphi_dyld,dphi_dfdr,df_dfs,dfs_dft,dphi_dft,
89 . dphi_dfs,dfn_dlam,dfsh_dlam,dfg_dlam,dft_dlam,
90 . dfn,dfsh,dfg,dft,dyld_dpla,dyld_dtemp,dtemp_dlam
91c
92 my_real, DIMENSION(NEL) ::
93 . dsigxx,dsigyy,dsigxy,trsig,trdep,trdfds,
94 . sxx,syy,sxy,szz,sigm,j2,j3,sigdr,yld,weitemp,
95 . hardp,fhard,frate,ftherm,dtherm,fdr,phi0,triax,
96 . fdam,phi,ft,fs,fg,fn,fsh,dpla_dlam,dphi_dlam,dezz,etat,
97 . normxx,normyy,normxy,normzz,sig_dfdsig,
98 . dpxx,dpyy,dpxy,dpzz,sigdr2,yld2i,dlam_nl
99c
100 !=======================================================================
101 ! DRUCKER - VOCE - JOHNSON-COOK MATERIAL LAW WITH GURSON DAMAGE
102 ! USING NON LOCAL PEERLINGS METHOD
103 !=======================================================================
104 !UVAR(1) YLD YIELD STRESS
105 !=======================================================================
106c
107 !=======================================================================
108 ! - INITIALISATION OF COMPUTATION ON TIME STEP
109 !=======================================================================
110 ! Recovering model parameters
111 ! Elastic parameters
112 young = uparam(1)
113 bulk = uparam(2)
114 g = uparam(3)
115 g2 = uparam(4)
116 lam = uparam(5)
117 nu = uparam(6)
118 nnu = uparam(7)
119 a11 = uparam(9)
120 a12 = uparam(10)
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)
143c
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 ! User inputs
166 yld(i) = uvar(i,1) ! Previous yield stress
167 ! Damage variables
168 fg(i) = dmg(i,2) ! growth damage
169 fn(i) = dmg(i,3) ! Nucleation damage
170 fsh(i) = dmg(i,4) ! Shear damage
171 ft(i) = dmg(i,5) ! Total damage
172 fs(i) = dmg(i,6) ! Effective damage
173 ! Standard inputs
174 dpla(i) = zero ! Initialization of the plastic strain increment
175 et(i) = one ! Initialization of hourglass coefficient
176 hardp(i) = zero ! Initialization of hardening modulus
177 dezz(i) = zero ! Initialization of the transverse strain
178 dlam_nl(i) = zero ! Initialization of the non-local plastic multiplier
179 sig_dfdsig(i) = zero ! Initialization of the stress - normal scalar product
180 normxx(i) = zero ! Initialization of the x component of the normal tensor
181 normyy(i) = zero ! Initialization of the y component of the normal tensor
182 normzz(i) = zero ! Initialization of the z component of the normal tensor
183 ENDDO
184c
185 ! Initialization of damage, temperature and self-heating weight factor
186 IF (time == zero) THEN
187 temp(1:nel) = tini
188 IF (isigi == 0) THEN
189 dmg(1:nel,5) = f0
190 ft(1:nel) = f0
191 dmg(1:nel,1) = f0/fr
192 IF (f0<fc) THEN
193 dmg(1:nel,6) = f0
194 ELSE
195 dmg(1:nel,6) = fc + (one/q1-fc)*(f0-fc)/(fr-fc)
196 ENDIF
197 fs(1:nel) = dmg(1:nel,6)
198 ENDIF
199 ENDIF
200 IF (cp > zero) THEN
201 IF (jthe == 0) THEN
202 DO i=1,nel
203 ! Computation of the weight factor
204 IF (plap_nl(i) < dpis) THEN
205 weitemp(i) = zero
206 ELSEIF (plap_nl(i) > dpad) THEN
207 weitemp(i) = one
208 ELSE
209 weitemp(i) = ((plap_nl(i)-dpis)**2 )
210 . * (three*dpad - two*plap_nl(i) - dpis)
211 . / ((dpad-dpis)**3)
212 ENDIF
213 ENDDO
214 ELSE
215 temp(1:nel) = tempel(1:nel)
216 ENDIF
217 ENDIF
218c
219 !========================================================================
220 ! NON-LOCAL VARIABLES UPDATE
221 !========================================================================
222 DO i=1,nel
223c
224 ! Previous value of Drucker equivalent stress
225 trsig(i) = sigoxx(i) + sigoyy(i)
226 sigm(i) = -trsig(i) * third
227 sxx(i) = sigoxx(i) + sigm(i)
228 syy(i) = sigoyy(i) + sigm(i)
229 szz(i) = sigm(i)
230 sxy(i) = sigoxy(i)
231 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) + sxy(i)**2
232 j3(i) = sxx(i) * syy(i) * szz(i) - szz(i) * sxy(i)**2
233 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
234 IF (fdr(i) > zero) THEN
235 sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
236 ELSE
237 sigdr(i) = zero
238 ENDIF
239c
240 ! Computation of the stress triaxiality and the etaT factor
241 IF (sigdr(i)>zero) THEN
242 triax(i) = (trsig(i)*third)/sigdr(i)
243 ELSE
244 triax(i) = zero
245 ENDIF
246 IF (trsig(i)<zero) THEN
247 etat(i) = zero
248 ELSE
249 etat(i) = one
250 ENDIF
251c
252 ! Normal to the previous yield surface
253 IF (yld(i)>zero) THEN
254 yld2i(i) = one / yld(i)**2
255 dphi_dsig = two * sigdr(i) * yld2i(i)
256 fsinh = sinh(q2*etat(i)*trsig(i)/(yld(i)*two))
257 dphi_dtrsig = q1*q2*etat(i)*fs(i)*fsinh/yld(i)
258 ELSE
259 yld2i(i) = zero
260 dphi_dsig = zero
261 fsinh = zero
262 dphi_dtrsig = zero
263 ENDIF
264c
265 ! Computation of the Eulerian norm of the stress tensor
266 normsig = sqrt(sigoxx(i)*sigoxx(i)
267 . + sigoyy(i)*sigoyy(i)
268 . + two*sigoxy(i)*sigoxy(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))/(normsig**2)
282 . - third*(sxx(i)*szz(i))/(normsig**2)
283 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
284 dj3dsyy = - third*(syy(i)*szz(i))/(normsig**2)
285 . + two_third*(sxx(i)*szz(i))/(normsig**2)
286 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
287 dj3dszz = - third*(syy(i)*szz(i))/(normsig**2)
288 . - third*(sxx(i)*szz(i))/(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))/(normsig**2)
291c
292 ! dPhi/dSig
293 normxx(i) = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx + dphi_dtrsig
294 normyy(i) = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy + dphi_dtrsig
295 normzz(i) = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz + dphi_dtrsig
296 normxy(i) = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
297 trdfds(i) = normxx(i) + normyy(i) + normzz(i)
298 sig_dfdsig(i) = sigoxx(i)*normxx(i) + sigoyy(i)*normyy(i)
299 . + sigoxy(i)*normxy(i)
300c
301 ! Computation of the non-local plastic multiplier
302 IF (sig_dfdsig(i)>zero) THEN
303 dlam_nl(i) = ((one - ft(i))*yld(i)*dpla_nl(i))/max(sig_dfdsig(i),em20)
304 ELSE
305 dlam_nl(i) = zero
306 ENDIF
307c
308 ! Damage growth update
309 IF ((trdfds(i)>zero).AND.(ft(i)>zero).AND.(ft(i)<fr)) THEN
310 fg(i) = fg(i) + (one-ft(i))*dlam_nl(i)*trdfds(i)
311 ENDIF
312 fg(i) = max(fg(i),zero)
313c
314 ! Nucleation damage update
315 IF ((pla_nl(i) >= ed).AND.(ft(i)<fr)) THEN
316 ! Case for positive stress triaxiality
317 IF (triax(i)>=zero) THEN
318 fn(i) = fn(i) + an*dpla_nl(i)
319 ! Case for negative stress triaxiality
320 ELSEIF ((triax(i)<zero).AND.(triax(i)>=-third)) THEN
321 fn(i) = fn(i) + an*max(one + three*triax(i),zero)*dpla_nl(i)
322 ENDIF
323 ENDIF
324 fn(i) = max(fn(i),zero)
325c
326 ! Shear damage update
327 IF ((sigdr(i) > zero).AND.(ft(i)>zero).AND.(ft(i)<fr)) THEN
328 sigvm = sqrt(max(em20,three*(j2(i)/(normsig**2))))
329 omega = one - ((twenty7 *(j3(i)/(normsig**3)))/(two*(sigvm**3)))**2
330 omega = max(omega,zero)
331 omega = min(omega,one)
332 sdpla = (sxx(i)*normxx(i)+syy(i)*normyy(i)+ szz(i)*normzz(i)
333 . + sxy(i)*normxy(i))
334 . * dlam_nl(i)
335 fsh(i) = fsh(i) + kw*omega*ft(i)*(sdpla/sigdr(i))
336 ENDIF
337 fsh(i) = max(fsh(i),zero)
338c
339 ! Total damage update
340 ft(i) = f0 + fg(i) + fn(i) + fsh(i)
341 ft(i) = min(ft(i), fr)
342 IF (ft(i) >= fr) THEN
343 IF (off(i)==one) off(i) = four_over_5
344 ENDIF
345c
346 ! Effective update
347 IF (ft(i) < fc)THEN
348 fs(i) = ft(i)
349 ELSE
350 fs(i) = fc + (one/q1 - fc) * (ft(i)-fc)/(fr-fc)
351 ENDIF
352 fs(i) = min(fs(i),one/q1)
353c
354 ! Temperature update
355 IF (cp > zero ) THEN
356 IF (jthe == 0) THEN
357 dtemp = weitemp(i)*(one-ft(i))*yld(i)*dpla_nl(i)*eta/(rho(i)*cp)
358 temp(i) = temp(i) + dtemp
359 ENDIF
360 ENDIF
361 ENDDO
362c
363 ! Computation of the initial yield stress
364 DO i=1,nel
365 ! a) - Hardening law
366 fhard(i) = yld0 + hard*pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
367 ! b) - Correction factor for strain-rate dependence (Johnson-Cook)
368 frate(i) = one
369 IF (epsd(i) > epsp0) frate(i) = frate(i) + jcc*log(epsd(i)/epsp0)
370 ! c) - Correction factor for thermal effects
371 ftherm(i) = one
372 IF (cp > zero) 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 signxx(i) = sigoxx(i) + (a11*depsxx(i) + a12*depsyy(i))
386 signyy(i) = sigoyy(i) + (a11*depsyy(i) + a12*depsxx(i))
387 signxy(i) = sigoxy(i) + (depsxy(i)*g)
388 signyz(i) = sigoyz(i) + (depsyz(i)*gs(i))
389 signzx(i) = sigozx(i) + (depszx(i)*gs(i))
390 ! Computation of the trace of the trial stress tensor
391 trsig(i) = signxx(i) + signyy(i)
392 sigm(i) = -trsig(i) * third
393 ! Computation of the deviatoric trial stress tensor
394 sxx(i) = signxx(i) + sigm(i)
395 syy(i) = signyy(i) + sigm(i)
396 szz(i) = sigm(i)
397 sxy(i) = signxy(i)
398 ! second deviatoric invariant
399 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) + sxy(i)**2
400 ! Third deviatoric invariant
401 j3(i) = sxx(i)*syy(i)*szz(i) - szz(i)*sxy(i)**2
402 ! Drucker equivalent stress
403 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
404 ! Checking equivalent stress values
405 IF (fdr(i) > zero) THEN
406 sigdr(i) = kdr * exp((one/six)*log(fdr(i))) ! FDR(I)**(1/6)
407 ELSE
408 sigdr(i) = zero
409 ENDIF
410 ! Computation of the stress triaxiality and the etaT factor
411 IF (sigdr(i)>zero) THEN
412 triax(i) = (trsig(i)*third)/sigdr(i)
413 ELSE
414 triax(i) = zero
415 ENDIF
416 IF (trsig(i)<zero) THEN
417 etat(i) = zero
418 ELSE
419 etat(i) = one
420 ENDIF
421c
422 ! Update out of plane strain increment
423 IF (off(i) == one) THEN
424 dezz(i) = -nnu*depsxx(i)-nnu*depsyy(i)
425 IF (sig_dfdsig(i) > em01) THEN
426 dezz(i) = dezz(i) + nnu*(dlam_nl(i)*normxx(i))
427 . + nnu*(dlam_nl(i)*normyy(i))
428 . + dlam_nl(i)*normzz(i)
429 ENDIF
430 ENDIF
431 ENDDO
432c
433 !========================================================================
434 ! - COMPUTATION OF YIELD FONCTION
435 !========================================================================
436 DO i=1,nel
437 fdam(i) = two*q1*fs(i)*cosh(q2*etat(i)*trsig(i)/yld(i)/two) - (q1*fs(i))**2
438 phi(i) = (sigdr(i) / yld(i))**2 - one + fdam(i)
439 ENDDO
440c
441 ! Checking plastic behavior for all elements
442 nindx = 0
443 DO i=1,nel
444 IF ((phi(i) >= zero).AND.(off(i) == one).AND.(ft(i)<fr)) THEN
445 nindx=nindx+1
446 index(nindx)=i
447 ENDIF
448 ENDDO
449c
450 !====================================================================
451 ! - PLASTIC CORRECTION WITH CUTTING PLANE METHOD (SEMI-IMPLICIT)
452 !====================================================================
453c
454 ! Number of iterations
455 niter = 3
456c
457 ! Loop over yielding elements
458#include "vectorize.inc"
459 DO ii=1,nindx
460c
461 ! Number of the element with plastic behaviour
462 i = index(ii)
463c
464 ! Initialization of the iterative Newton procedure
465 sigdr2(i) = sigdr(i)**2
466 yld2i(i) = one / yld(i)**2
467 dpxx(i) = zero
468 dpyy(i) = zero
469 dpzz(i) = zero
470 dpxy(i) = zero
471 ENDDO
472c
473 ! Loop over the iterations
474 DO iter = 1, niter
475#include "vectorize.inc"
476 DO ii=1,nindx
477 i = index(ii)
478c
479 ! Note: in this part, the purpose is to compute for each iteration
480 ! a plastic multiplier allowing to update internal variables to satisfy
481 ! the consistency condition using the cutting plane semi-implicit
482 ! iterative procedure.
483 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
484 ! -> PHI : current value of yield function (known)
485 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
486 ! into account of internal variables kinetic :
487 ! plasticity, temperature and damage (to compute)
488c
489 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
490 !-------------------------------------------------------------
491c
492 ! Derivative with respect to the equivalent stress and trace
493 yld2i(i) = one/(yld(i)**2)
494 dphi_dsig = two*sigdr(i)*yld2i(i)
495 fsinh = sinh(q2*etat(i)*trsig(i)/(yld(i)*two))
496 dphi_dtrsig = q1*q2*etat(i)*fs(i)*fsinh/yld(i)
497c
498 ! Computation of the Eulerian norm of the stress tensor
499 normsig = sqrt(signxx(i)*signxx(i)
500 . + signyy(i)*signyy(i)
501 . + two*signxy(i)*signxy(i))
502 normsig = max(normsig,one)
503c
504 ! Derivative with respect to Fdr
505 fdr(i) = (j2(i)/(normsig**2))**3 - cdr*((j3(i)/(normsig**3))**2)
506 dphi_dfdr = dphi_dsig*kdr*(one/six)*exp(-(five/six)*log(fdr(i)))
507 dsdrdj2 = dphi_dfdr*three*(j2(i)/(normsig**2))**2
508 dsdrdj3 = -dphi_dfdr*two*cdr*(j3(i)/(normsig**3))
509 ! dJ3/dSig
510 dj3dsxx = two_third*(syy(i)*szz(i))/(normsig**2)
511 . - third*(sxx(i)*szz(i))/(normsig**2)
512 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
513 dj3dsyy = - third*(syy(i)*szz(i))/(normsig**2)
514 . + two_third*(sxx(i)*szz(i))/(normsig**2)
515 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
516 dj3dszz = - third*(syy(i)*szz(i))/(normsig**2)
517 . - third*(sxx(i)*szz(i))/(normsig**2)
518 . + two_third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
519 dj3dsxy = two*(sxx(i)*sxy(i) + sxy(i)*syy(i))/(normsig**2)
520 ! dPhi/dSig
521 normxx(i) = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx + dphi_dtrsig
522 normyy(i) = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy + dphi_dtrsig
523 normzz(i) = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz + dphi_dtrsig
524 normxy(i) = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
525c
526 ! 2 - Computation of DPHI_DLAMBDA
527 !---------------------------------------------------------
528c
529 ! a) Derivative with respect stress increments tensor DSIG
530 ! --------------------------------------------------------
531 trdfds(i) = normxx(i) + normyy(i) + normzz(i)
532 dfdsig2 = normxx(i) * (a11*normxx(i) + a12*normyy(i))
533 . + normyy(i) * (a11*normyy(i) + a12*normxx(i))
534 . + normxy(i) * normxy(i) * g
535c
536 ! b) Derivatives with respect to plastic strain P
537 ! ------------------------------------------------
538c
539 ! i) Derivative of the yield stress with respect to plastic strain dYLD / dPLA
540 ! ----------------------------------------------------------------------------
541 hardp(i) = hard + qvoce*bvoce*exp(-bvoce*pla(i))
542 dyld_dpla = hardp(i)*frate(i)*ftherm(i)
543c
544 ! ii) Derivative of dPLA with respect to DLAM
545 ! -------------------------------------------
546 sig_dfdsig(i) = signxx(i) * normxx(i)
547 . + signyy(i) * normyy(i)
548 . + signxy(i) * normxy(i)
549 dpla_dlam(i) = sig_dfdsig(i) / ((one - ft(i))*yld(i))
550c
551 ! c) Derivative with respect to the yield stress
552 ! ----------------------------------------------
553 sigdr2(i) = sigdr(i)**2
554 dphi_dyld = -two*sigdr2(i)/yld(i)**3-dphi_dtrsig*trsig(i)/yld(i)
555c
556 ! d) Derivative of PHI with respect to DLAM
557 dphi_dlam(i) = - dfdsig2 + (dphi_dyld*dyld_dpla*dpla_dlam(i))
558 dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20) ,dphi_dlam(i))
559c
560 ! Computation of the plastic multiplier
561 dlam = -phi(i)/dphi_dlam(i)
562c
563 ! Plastic strains tensor update
564 dpxx(i) = dlam * normxx(i)
565 dpyy(i) = dlam * normyy(i)
566 dpzz(i) = dlam * normzz(i)
567 dpxy(i) = dlam * normxy(i)
568 trdep(i) = dpxx(i) + dpyy(i) + dpzz(i)
569c
570 ! Elasto-plastic stresses update
571 signxx(i) = signxx(i) - (a11*dpxx(i) + a12*dpyy(i))
572 signyy(i) = signyy(i) - (a11*dpyy(i) + a12*dpxx(i))
573 signxy(i) = signxy(i) - dpxy(i)*g
574 trsig(i) = signxx(i) + signyy(i)
575 sigm(i) = -trsig(i) * third
576 sxx(i) = signxx(i) + sigm(i)
577 syy(i) = signyy(i) + sigm(i)
578 szz(i) = sigm(i)
579 sxy(i) = signxy(i)
580c
581 ! Cumulated plastic strain and strain rate update
582 ddep = (dlam/((one - ft(i))*yld(i)))*sig_dfdsig(i)
583 dpla(i) = max(zero, dpla(i) + ddep)
584 pla(i) = pla(i) + ddep
585c
586 ! Drucker equivalent stress update
587 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) + sxy(i)**2
588 j3(i) = sxx(i) * syy(i) * szz(i) - szz(i) * sxy(i)**2
589 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
590 sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
591 ! Computation of the stress triaxiality and the etaT factor
592 triax(i) = (trsig(i)*third)/sigdr(i)
593 IF (trsig(i)<zero) THEN
594 etat(i) = zero
595 ELSE
596 etat(i) = one
597 ENDIF
598c
599 ! Hardening law update
600 fhard(i) = yld0 + hard * pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
601c
602 ! Yield stress update
603 yld(i) = fhard(i) * frate(i) * ftherm(i)
604 yld(i) = max(yld(i), em10)
605c
606 ! Yield function value update
607 sigdr2(i) = sigdr(i)**2
608 yld2i(i) = one / yld(i)**2
609 fcosh = cosh(q2*etat(i)*trsig(i)/(yld(i)*two))
610 fdam(i) = two*q1*fs(i)*fcosh - (q1*fs(i))**2
611 phi(i) = sigdr2(i) * yld2i(i) - one + fdam(i)
612c
613 ENDDO
614 ENDDO
615 ! End of the loop over the iterations
616 !===================================================================
617 ! - END OF PLASTIC CORRECTION WITH CUTTING PLANE ITERATIVE METHOD
618 !===================================================================
619c
620 ! Storing new values
621 DO i=1,nel
622 ! USR Outputs
623 uvar(i,1) = yld(i) ! Yield stress
624 ! Standard outputs
625 dmg(i,1) = ft(i)/fr ! Normalized total damage
626 dmg(i,2) = fg(i) ! Growth damage
627 dmg(i,3) = fn(i) ! Nucleation damage
628 dmg(i,4) = fsh(i) ! Shear damage
629 dmg(i,5) = min(ft(i),fr) ! Total damage
630 dmg(i,6) = min(fs(i),one/q1) ! Effective damage
631 seq(i) = sigdr(i) ! Equivalent stress
632 ! If element is broken
633 IF (ft(i) >= fr) THEN
634 dpla(i) = zero
635 signxx(i) = zero
636 signyy(i) = zero
637 signxy(i) = zero
638 signyz(i) = zero
639 signzx(i) = zero
640 seq(i) = zero
641 ENDIF
642 ! Plastic strain-rate (filtered)
643 dpdt = dpla(i) / max(em20,timestep)
644 epsd(i) = afiltr * dpdt + (one - afiltr) * epsd(i)
645 ! Coefficient for hourglass
646 IF (dpla(i) > zero) THEN
647 et(i) = hardp(i)*frate(i) / (hardp(i)*frate(i) + young)
648 ELSE
649 et(i) = one
650 ENDIF
651 ! Computation of the sound speed
652 soundsp(i) = sqrt((a11)/rho(i))
653 ! Storing the yield stress
654 sigy(i) = yld(i)
655 ! Computation of the thickness variation
656 thk(i) = thk(i) + dezz(i)*thkly(i)*off(i)
657 ENDDO
658c
659 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
mat104c_nldam_newton
subroutine mat104c_nldam_newton(nel, ngl, nuparam, nuvar, time, timestep, uparam, uvar, jthe, off, gs, rho, pla, dpla, epsd, soundsp, depsxx, depsyy, depsxy, depsyz, depszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, thkly, thk, sigy, et, tempel, dpla_nl, dmg, temp, seq, pla_nl, plap_nl)
Definition mat104c_nldam_newton.F:37
second
real function second()
SECOND Using ETIME
Definition second_EXT_ETIME.f:35