OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat104c_nldam_newton.F File Reference
#include "implicit_f.inc"
#include "com01_c.inc"
#include "vectorize.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine 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)

Function/Subroutine Documentation

◆ mat104c_nldam_newton()

subroutine mat104c_nldam_newton ( integer nel,
integer, dimension(nel), intent(in) ngl,
integer nuparam,
integer nuvar,
time,
timestep,
intent(in) uparam,
intent(inout) uvar,
integer jthe,
intent(inout) off,
intent(in) gs,
intent(in) rho,
intent(inout) pla,
intent(inout) dpla,
intent(inout) epsd,
intent(out) soundsp,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
intent(in) sigoxx,
intent(in) sigoyy,
intent(in) sigoxy,
intent(in) sigoyz,
intent(in) sigozx,
intent(out) signxx,
intent(out) signyy,
intent(out) signxy,
intent(out) signyz,
intent(out) signzx,
intent(in) thkly,
intent(inout) thk,
intent(out) sigy,
intent(out) et,
intent(in) tempel,
intent(in) dpla_nl,
intent(inout) dmg,
intent(inout) temp,
intent(inout) seq,
intent(in) pla_nl,
intent(in) plap_nl )

Definition at line 28 of file mat104c_nldam_newton.F.

37 !=======================================================================
38 ! Implicit types
39 !=======================================================================
40#include "implicit_f.inc"
41 !=======================================================================
42 ! Common
43 !=======================================================================
44#include "com01_c.inc"
45 !=======================================================================
46 ! Dummy arguments
47 !=======================================================================
48 INTEGER NEL,NUPARAM,NUVAR,JTHE
49 INTEGER ,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,kw,fr,fc,f0,a11,a12,nnu
79 my_real ::
80 . dti,sigvm,omega,
81 . dtemp,fcosh,fsinh,dpdt,dlam,ddep
82 my_real ::
83 . dsdrdj2,dsdrdj3,
84 . dj3dsxx,dj3dsyy,dj3dsxy,dj3dszz,
85 . normsig,
86 .
87 . sdpla,dphi_dtrsig,dfdsig2,
88 . dphi_dsig,dphi_dyld,dphi_dfdr,
89 .
90 . dyld_dpla
91c
92 my_real, DIMENSION(NEL) ::
93 . trsig,trdep,trdfds,
94 . sxx,syy,sxy,szz,sigm,j2,j3,sigdr,yld,weitemp,
95 . hardp,fhard,frate,ftherm,fdr,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 !UVAR(2) EPSD PLASTIC STRAIN RATE SAVED FOR FILTERING
106 !=======================================================================
107c
108 !=======================================================================
109 ! - INITIALISATION OF COMPUTATION ON TIME STEP
110 !=======================================================================
111 ! Recovering model parameters
112 ! Elastic parameters
113 young = uparam(1)
114 bulk = uparam(2)
115 g = uparam(3)
116 g2 = uparam(4)
117 lam = uparam(5)
118 nu = uparam(6)
119 nnu = uparam(7)
120 a11 = uparam(9)
121 a12 = uparam(10)
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)
144c
145 ! Gurson damage model parameters parameters
146 igurson = nint(uparam(30)) ! Gurson switch flag:
147 ! = 0 => no damage model
148 ! = 1 => local damage model
149 ! = 2 => non local (Forest - micromorphic) damage model
150 ! = 3 => non local (Peerlings) damage model
151 q1 = uparam(31) ! Gurson yield criterion 1st parameter
152 q2 = uparam(32) ! Gurson yield criterion 2nd parameter
153 ed = uparam(34) ! Plastic strain trigger for nucleation
154 an = uparam(35) ! Nucleation rate
155 kw = uparam(36) ! Shear coefficient (Nahshon & Hutchinson)
156 fr = uparam(37) ! Failure void volume fraction
157 fc = uparam(38) ! Critical void volume fraction
158 f0 = uparam(39) ! Initial void volume fraction
159 hkhi = uparam(40) ! Micromorphic penalty parameter
160c
161 ! Recovering internal variables
162 DO i=1,nel
163 ! If the element is failing
164 IF (off(i) < em03) off(i) = zero
165 IF (off(i) < one) off(i) = off(i)*four_over_5
166 ! User inputs
167 yld(i) = uvar(i,1) ! Previous yield stress
168 ! Damage variables
169 fg(i) = dmg(i,2) ! Growth damage
170 fn(i) = dmg(i,3) ! Nucleation damage
171 fsh(i) = dmg(i,4) ! Shear damage
172 ft(i) = dmg(i,5) ! Total damage
173 fs(i) = dmg(i,6) ! Effective damage
174 ! Standard inputs
175 dpla(i) = zero ! Initialization of the plastic strain increment
176 et(i) = one ! Initialization of hourglass coefficient
177 hardp(i) = zero ! Initialization of hardening modulus
178 dezz(i) = zero ! Initialization of the transverse strain
179 dlam_nl(i) = zero ! Initialization of the non-local plastic multiplier
180 sig_dfdsig(i) = zero ! Initialization of the stress - normal scalar product
181 normxx(i) = zero ! Initialization of the x component of the normal tensor
182 normyy(i) = zero ! Initialization of the y component of the normal tensor
183 normzz(i) = zero ! Initialization of the z component of the normal tensor
184 ENDDO
185!
186 epsd(1:nel) = uvar(1:nel,2)
187!
188 ! Initialization of damage, temperature and self-heating weight factor
189 IF (time == zero) THEN
190 temp(1:nel) = tini
191 IF (isigi == 0) THEN
192 dmg(1:nel,5) = f0
193 ft(1:nel) = f0
194 dmg(1:nel,1) = f0/fr
195 IF (f0<fc) THEN
196 dmg(1:nel,6) = f0
197 ELSE
198 dmg(1:nel,6) = fc + (one/q1-fc)*(f0-fc)/(fr-fc)
199 ENDIF
200 fs(1:nel) = dmg(1:nel,6)
201 ENDIF
202 ENDIF
203 IF (cp > zero) THEN
204 IF (jthe == 0) THEN
205 DO i=1,nel
206 ! Computation of the weight factor
207 IF (plap_nl(i) < dpis) THEN
208 weitemp(i) = zero
209 ELSEIF (plap_nl(i) > dpad) THEN
210 weitemp(i) = one
211 ELSE
212 weitemp(i) = ((plap_nl(i)-dpis)**2 )
213 . * (three*dpad - two*plap_nl(i) - dpis)
214 . / ((dpad-dpis)**3)
215 ENDIF
216 ENDDO
217 ELSE
218 temp(1:nel) = tempel(1:nel)
219 ENDIF
220 ENDIF
221c
222 !========================================================================
223 ! NON-LOCAL VARIABLES UPDATE
224 !========================================================================
225 DO i=1,nel
226c
227 ! Previous value of Drucker equivalent stress
228 trsig(i) = sigoxx(i) + sigoyy(i)
229 sigm(i) = -trsig(i) * third
230 sxx(i) = sigoxx(i) + sigm(i)
231 syy(i) = sigoyy(i) + sigm(i)
232 szz(i) = sigm(i)
233 sxy(i) = sigoxy(i)
234 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) + sxy(i)**2
235 j3(i) = sxx(i) * syy(i) * szz(i) - szz(i) * sxy(i)**2
236 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
237 IF (fdr(i) > zero) THEN
238 sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
239 ELSE
240 sigdr(i) = zero
241 ENDIF
242c
243 ! Computation of the stress triaxiality and the etaT factor
244 IF (sigdr(i)>zero) THEN
245 triax(i) = (trsig(i)*third)/sigdr(i)
246 ELSE
247 triax(i) = zero
248 ENDIF
249 IF (trsig(i)<zero) THEN
250 etat(i) = zero
251 ELSE
252 etat(i) = one
253 ENDIF
254c
255 ! Normal to the previous yield surface
256 IF (yld(i)>zero) THEN
257 yld2i(i) = one / yld(i)**2
258 dphi_dsig = two * sigdr(i) * yld2i(i)
259 fsinh = sinh(q2*etat(i)*trsig(i)/(yld(i)*two))
260 dphi_dtrsig = q1*q2*etat(i)*fs(i)*fsinh/yld(i)
261 ELSE
262 yld2i(i) = zero
263 dphi_dsig = zero
264 fsinh = zero
265 dphi_dtrsig = zero
266 ENDIF
267c
268 ! Computation of the Eulerian norm of the stress tensor
269 normsig = sqrt(sigoxx(i)*sigoxx(i)
270 . + sigoyy(i)*sigoyy(i)
271 . + two*sigoxy(i)*sigoxy(i))
272 normsig = max(normsig,one)
273c
274 ! Computation of the normal to the yield surface
275 fdr(i) = (j2(i)/(normsig**2))**3 - cdr*((j3(i)/(normsig**3))**2)
276 IF (fdr(i) > zero) THEN
277 dphi_dfdr = dphi_dsig*kdr*(one/six)*exp(-(five/six)*log(fdr(i)))
278 ELSE
279 dphi_dfdr = zero
280 ENDIF
281 dsdrdj2 = dphi_dfdr*three*(j2(i)/(normsig**2))**2
282 dsdrdj3 = -dphi_dfdr*two*cdr*(j3(i)/(normsig**3))
283 !dj3/ds
284 dj3dsxx = two_third*(syy(i)*szz(i))/(normsig**2)
285 . - third*(sxx(i)*szz(i))/(normsig**2)
286 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
287 dj3dsyy = - third*(syy(i)*szz(i))/(normsig**2)
288 . + two_third*(sxx(i)*szz(i))/(normsig**2)
289 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
290 dj3dszz = - third*(syy(i)*szz(i))/(normsig**2)
291 . - third*(sxx(i)*szz(i))/(normsig**2)
292 . + two_third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
293 dj3dsxy = two*(sxx(i)*sxy(i) + sxy(i)*syy(i))/(normsig**2)
294c
295 ! dPhi/dSig
296 normxx(i) = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx + dphi_dtrsig
297 normyy(i) = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy + dphi_dtrsig
298 normzz(i) = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz + dphi_dtrsig
299 normxy(i) = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
300 trdfds(i) = normxx(i) + normyy(i) + normzz(i)
301 sig_dfdsig(i) = sigoxx(i)*normxx(i) + sigoyy(i)*normyy(i)
302 . + sigoxy(i)*normxy(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))/max(sig_dfdsig(i),em20)
307 ELSE
308 dlam_nl(i) = zero
309 ENDIF
310c
311 ! Damage growth update
312 IF ((trdfds(i)>zero).AND.(ft(i)>zero).AND.(ft(i)<fr)) 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))
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 (cp > zero ) THEN
359 IF (jthe == 0) THEN
360 dtemp = weitemp(i)*(one-ft(i))*yld(i)*dpla_nl(i)*eta/(rho(i)*cp)
361 temp(i) = temp(i) + dtemp
362 ENDIF
363 ENDIF
364 ENDDO
365c
366 ! Computation of the initial yield stress
367 DO i=1,nel
368 ! a) - Hardening law
369 fhard(i) = yld0 + hard*pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
370 ! b) - Correction factor for strain-rate dependence (Johnson-Cook)
371 frate(i) = one
372 IF (epsd(i) > epsp0) frate(i) = frate(i) + jcc*log(epsd(i)/epsp0)
373 ! c) - Correction factor for thermal effects
374 ftherm(i) = one
375 IF (cp > zero) ftherm(i) = one - mtemp*(temp(i) - tref)
376 ! d) - Computation of the yield function
377 yld(i) = fhard(i)*frate(i)*ftherm(i)
378 ! e) - Checking values
379 yld(i) = max(em10, yld(i))
380 ENDDO
381c
382 !========================================================================
383 ! - COMPUTATION OF TRIAL VALUES
384 !========================================================================
385 DO i=1,nel
386c
387 ! Computation of the trial stress tensor
388 signxx(i) = sigoxx(i) + (a11*depsxx(i) + a12*depsyy(i))
389 signyy(i) = sigoyy(i) + (a11*depsyy(i) + a12*depsxx(i))
390 signxy(i) = sigoxy(i) + (depsxy(i)*g)
391 signyz(i) = sigoyz(i) + (depsyz(i)*gs(i))
392 signzx(i) = sigozx(i) + (depszx(i)*gs(i))
393 ! Computation of the trace of the trial stress tensor
394 trsig(i) = signxx(i) + signyy(i)
395 sigm(i) = -trsig(i) * third
396 ! Computation of the deviatoric trial stress tensor
397 sxx(i) = signxx(i) + sigm(i)
398 syy(i) = signyy(i) + sigm(i)
399 szz(i) = sigm(i)
400 sxy(i) = signxy(i)
401 ! Second deviatoric invariant
402 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) + sxy(i)**2
403 ! Third deviatoric invariant
404 j3(i) = sxx(i)*syy(i)*szz(i) - szz(i)*sxy(i)**2
405 ! Drucker equivalent stress
406 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
407 ! Checking equivalent stress values
408 IF (fdr(i) > zero) THEN
409 sigdr(i) = kdr * exp((one/six)*log(fdr(i))) ! FDR(I)**(1/6)
410 ELSE
411 sigdr(i) = zero
412 ENDIF
413 ! Computation of the stress triaxiality and the etaT factor
414 IF (sigdr(i)>zero) THEN
415 triax(i) = (trsig(i)*third)/sigdr(i)
416 ELSE
417 triax(i) = zero
418 ENDIF
419 IF (trsig(i)<zero) THEN
420 etat(i) = zero
421 ELSE
422 etat(i) = one
423 ENDIF
424c
425 ! Update out of plane strain increment
426 IF (off(i) == one) THEN
427 dezz(i) = -nnu*depsxx(i)-nnu*depsyy(i)
428 IF (sig_dfdsig(i) > em01) THEN
429 dezz(i) = dezz(i) + nnu*(dlam_nl(i)*normxx(i))
430 . + nnu*(dlam_nl(i)*normyy(i))
431 . + dlam_nl(i)*normzz(i)
432 ENDIF
433 ENDIF
434 ENDDO
435c
436 !========================================================================
437 ! - COMPUTATION OF YIELD FONCTION
438 !========================================================================
439 DO i=1,nel
440 fdam(i) = two*q1*fs(i)*cosh(q2*etat(i)*trsig(i)/yld(i)/two) - (q1*fs(i))**2
441 phi(i) = (sigdr(i) / yld(i))**2 - one + fdam(i)
442 ENDDO
443c
444 ! Checking plastic behavior for all elements
445 nindx = 0
446 DO i=1,nel
447 IF ((phi(i) >= zero).AND.(off(i) == one).AND.(ft(i)<fr)) THEN
448 nindx=nindx+1
449 index(nindx)=i
450 ENDIF
451 ENDDO
452c
453 !====================================================================
454 ! - PLASTIC CORRECTION WITH CUTTING PLANE METHOD (SEMI-IMPLICIT)
455 !====================================================================
456c
457 ! Number of iterations
458 niter = 3
459c
460 ! Loop over yielding elements
461#include "vectorize.inc"
462 DO ii=1,nindx
463c
464 ! Number of the element with plastic behaviour
465 i = index(ii)
466c
467 ! Initialization of the iterative Newton procedure
468 sigdr2(i) = sigdr(i)**2
469 yld2i(i) = one / yld(i)**2
470 dpxx(i) = zero
471 dpyy(i) = zero
472 dpzz(i) = zero
473 dpxy(i) = zero
474 ENDDO
475c
476 ! Loop over the iterations
477 DO iter = 1, niter
478#include "vectorize.inc"
479 DO ii=1,nindx
480 i = index(ii)
481c
482 ! Note: in this part, the purpose is to compute for each iteration
483 ! a plastic multiplier allowing to update internal variables to satisfy
484 ! the consistency condition using the cutting plane semi-implicit
485 ! iterative procedure.
486 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
487 ! -> PHI : current value of yield function (known)
488 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
489 ! into account of internal variables kinetic :
490 ! plasticity, temperature and damage (to compute)
491c
492 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
493 !-------------------------------------------------------------
494c
495 ! Derivative with respect to the equivalent stress and trace
496 yld2i(i) = one/(yld(i)**2)
497 dphi_dsig = two*sigdr(i)*yld2i(i)
498 fsinh = sinh(q2*etat(i)*trsig(i)/(yld(i)*two))
499 dphi_dtrsig = q1*q2*etat(i)*fs(i)*fsinh/yld(i)
500c
501 ! Computation of the Eulerian norm of the stress tensor
502 normsig = sqrt(signxx(i)*signxx(i)
503 . + signyy(i)*signyy(i)
504 . + two*signxy(i)*signxy(i))
505 normsig = max(normsig,one)
506c
507 ! Derivative with respect to Fdr
508 fdr(i) = (j2(i)/(normsig**2))**3 - cdr*((j3(i)/(normsig**3))**2)
509 dphi_dfdr = dphi_dsig*kdr*(one/six)*exp(-(five/six)*log(fdr(i)))
510 dsdrdj2 = dphi_dfdr*three*(j2(i)/(normsig**2))**2
511 dsdrdj3 = -dphi_dfdr*two*cdr*(j3(i)/(normsig**3))
512 ! dJ3/dSig
513 dj3dsxx = two_third*(syy(i)*szz(i))/(normsig**2)
514 . - third*(sxx(i)*szz(i))/(normsig**2)
515 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
516 dj3dsyy = - third*(syy(i)*szz(i))/(normsig**2)
517 . + two_third*(sxx(i)*szz(i))/(normsig**2)
518 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
519 dj3dszz = - third*(syy(i)*szz(i))/(normsig**2)
520 . - third*(sxx(i)*szz(i))/(normsig**2)
521 . + two_third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
522 dj3dsxy = two*(sxx(i)*sxy(i) + sxy(i)*syy(i))/(normsig**2)
523 ! dPhi/dSig
524 normxx(i) = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx + dphi_dtrsig
525 normyy(i) = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy + dphi_dtrsig
526 normzz(i) = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz + dphi_dtrsig
527 normxy(i) = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
528c
529 ! 2 - Computation of DPHI_DLAMBDA
530 !---------------------------------------------------------
531c
532 ! a) Derivative with respect stress increments tensor DSIG
533 ! --------------------------------------------------------
534 trdfds(i) = normxx(i) + normyy(i) + normzz(i)
535 dfdsig2 = normxx(i) * (a11*normxx(i) + a12*normyy(i))
536 . + normyy(i) * (a11*normyy(i) + a12*normxx(i))
537 . + normxy(i) * normxy(i) * g
538c
539 ! b) Derivatives with respect to plastic strain P
540 ! ------------------------------------------------
541c
542 ! i) Derivative of the yield stress with respect to plastic strain dYLD / dPLA
543 ! ----------------------------------------------------------------------------
544 hardp(i) = hard + qvoce*bvoce*exp(-bvoce*pla(i))
545 dyld_dpla = hardp(i)*frate(i)*ftherm(i)
546c
547 ! ii) Derivative of dPLA with respect to DLAM
548 ! -------------------------------------------
549 sig_dfdsig(i) = signxx(i) * normxx(i)
550 . + signyy(i) * normyy(i)
551 . + signxy(i) * normxy(i)
552 dpla_dlam(i) = sig_dfdsig(i) / ((one - ft(i))*yld(i))
553c
554 ! c) Derivative with respect to the yield stress
555 ! ----------------------------------------------
556 sigdr2(i) = sigdr(i)**2
557 dphi_dyld = -two*sigdr2(i)/yld(i)**3-dphi_dtrsig*trsig(i)/yld(i)
558c
559 ! d) Derivative of PHI with respect to DLAM
560 dphi_dlam(i) = - dfdsig2 + (dphi_dyld*dyld_dpla*dpla_dlam(i))
561 dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20) ,dphi_dlam(i))
562c
563 ! Computation of the plastic multiplier
564 dlam = -phi(i)/dphi_dlam(i)
565c
566 ! Plastic strains tensor update
567 dpxx(i) = dlam * normxx(i)
568 dpyy(i) = dlam * normyy(i)
569 dpzz(i) = dlam * normzz(i)
570 dpxy(i) = dlam * normxy(i)
571 trdep(i) = dpxx(i) + dpyy(i) + dpzz(i)
572c
573 ! Elasto-plastic stresses update
574 signxx(i) = signxx(i) - (a11*dpxx(i) + a12*dpyy(i))
575 signyy(i) = signyy(i) - (a11*dpyy(i) + a12*dpxx(i))
576 signxy(i) = signxy(i) - dpxy(i)*g
577 trsig(i) = signxx(i) + signyy(i)
578 sigm(i) = -trsig(i) * third
579 sxx(i) = signxx(i) + sigm(i)
580 syy(i) = signyy(i) + sigm(i)
581 szz(i) = sigm(i)
582 sxy(i) = signxy(i)
583c
584 ! Cumulated plastic strain and strain rate update
585 ddep = (dlam/((one - ft(i))*yld(i)))*sig_dfdsig(i)
586 dpla(i) = max(zero, dpla(i) + ddep)
587 pla(i) = pla(i) + ddep
588c
589 ! Drucker equivalent stress update
590 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) + sxy(i)**2
591 j3(i) = sxx(i) * syy(i) * szz(i) - szz(i) * sxy(i)**2
592 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
593 sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
594 ! Computation of the stress triaxiality and the etaT factor
595 triax(i) = (trsig(i)*third)/sigdr(i)
596 IF (trsig(i)<zero) THEN
597 etat(i) = zero
598 ELSE
599 etat(i) = one
600 ENDIF
601c
602 ! Hardening law update
603 fhard(i) = yld0 + hard * pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
604c
605 ! Yield stress update
606 yld(i) = fhard(i) * frate(i) * ftherm(i)
607 yld(i) = max(yld(i), em10)
608c
609 ! Yield function value update
610 sigdr2(i) = sigdr(i)**2
611 yld2i(i) = one / yld(i)**2
612 fcosh = cosh(q2*etat(i)*trsig(i)/(yld(i)*two))
613 fdam(i) = two*q1*fs(i)*fcosh - (q1*fs(i))**2
614 phi(i) = sigdr2(i) * yld2i(i) - one + fdam(i)
615c
616 ENDDO
617 ENDDO
618 ! End of the loop over the iterations
619 !===================================================================
620 ! - END OF PLASTIC CORRECTION WITH CUTTING PLANE ITERATIVE METHOD
621 !===================================================================
622c
623 ! Storing new values
624 DO i=1,nel
625 ! USR Outputs
626 uvar(i,1) = yld(i) ! Yield stress
627 ! Standard outputs
628 dmg(i,1) = ft(i)/fr ! Normalized total damage
629 dmg(i,2) = fg(i) ! Growth damage
630 dmg(i,3) = fn(i) ! Nucleation damage
631 dmg(i,4) = fsh(i) ! Shear damage
632 dmg(i,5) = min(ft(i),fr) ! Total damage
633 dmg(i,6) = min(fs(i),one/q1) ! Effective damage
634 seq(i) = sigdr(i) ! Equivalent stress
635 ! If element is broken
636 IF (ft(i) >= fr) THEN
637 dpla(i) = zero
638 signxx(i) = zero
639 signyy(i) = zero
640 signxy(i) = zero
641 signyz(i) = zero
642 signzx(i) = zero
643 seq(i) = zero
644 ENDIF
645 ! Plastic strain-rate (filtered)
646 dpdt = dpla(i) / max(em20,timestep)
647 epsd(i) = afiltr * dpdt + (one - afiltr) * epsd(i)
648 uvar(i,2) = epsd(i)
649 ! Coefficient for hourglass
650 IF (dpla(i) > zero) THEN
651 et(i) = hardp(i)*frate(i) / (hardp(i)*frate(i) + young)
652 ELSE
653 et(i) = one
654 ENDIF
655 ! Computation of the sound speed
656 soundsp(i) = sqrt((a11)/rho(i))
657 ! Storing the yield stress
658 sigy(i) = yld(i)
659 ! Computation of the thickness variation
660 thk(i) = thk(i) + dezz(i)*thkly(i)*off(i)
661 ENDDO
662c
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21