OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat104c_nldam_nice.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_nice (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_nice()

subroutine mat104c_nldam_nice ( 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_nice.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
70c
71 !=======================================================================
72 ! Local Variables
73 !=======================================================================
74 INTEGER I,II,IGURSON,NINDX,INDEX(NEL)
75c
76 my_real ::
77 . young,bulk,lam,g,g2,nu,cdr,kdr,hard,yld0,qvoce,bvoce,jcc,
78 . epsp0,mtemp,tini,tref,eta,cp,dpis,dpad,asrate,afiltr,hkhi,
79 . q1,q2,ed,an,kw,fr,fc,f0,a11,a12,nnu
80 my_real ::
81 . dti,sigvm,sigdr2,yld2i,omega,
82 . dtemp,fcosh,fsinh,dpdt,dlam,ddep
83 my_real ::
84 . dpxx,dpyy,dpxy,dpzz,dsdrdj2,dsdrdj3,
85 . dj3dsxx,dj3dsyy,dj3dsxy,dj3dszz,
86 .
87 .
88 . denom,dphi,sdpla,dfdsig2,
89 . dphi_dsig,dphi_dyld,dphi_dfdr,
90 . normsig,
91 . dyld_dpla
92c
93 my_real, DIMENSION(NEL) ::
94 . dsigxx,dsigyy,dsigxy,trsig,trdep,
95 . sxx,syy,sxy,szz,sigm,j2,j3,sigdr,yld,weitemp,trdfds,
96 . hardp,fhard,frate,ftherm,fdr,triax,etat,sigdr_old,
97 . fdam,phi0,phi,ft,fs,fg,fn,fsh,dpla_dlam,dezz,dphi_dtrsig,
98 . normxx,normyy,normxy,normzz,sig_dfdsig,dlam_nl
99c
100 !=======================================================================
101 ! DRUCKER - VOCE - JOHNSON-COOK MATERIAL LAW WITH GURSON DAMAGE
102 ! USING NON LOCAL PEERLINGS METHOD
103 !=======================================================================
104 !UVAR(1) PHI YIELD FUNCTION VALUE
105 !UVAR(2) YLD YIELD STRESS
106 !UVAR(3) EPSD PLASTIC STRAIN RATE SAVED FOR FILTERING
107 !=======================================================================
108c
109 !=======================================================================
110 ! - INITIALISATION OF COMPUTATION ON TIME STEP
111 !=======================================================================
112 ! Recovering model parameters
113 ! Elastic parameters
114 young = uparam(1) ! young modulus
115 bulk = uparam(2) ! Bulk modulus
116 g = uparam(3) ! Shear modulus
117 g2 = uparam(4) ! 2*Shear modulus
118 lam = uparam(5) ! Lambda Hooke parameter
119 nu = uparam(6) ! Poisson ration
120 nnu = uparam(7) ! NU/(1-NU)
121 a11 = uparam(9) ! Diagonal term, elastic matrix in plane stress
122 a12 = uparam(10) ! Non-diagonal term, elastic matrix in plane stress
123 ! Plastic criterion and hardening parameters [Drucker, 1948]
124 cdr = uparam(12) ! Drucker coefficient
125 kdr = uparam(13) ! Drucker 1/K coefficient
126 tini = uparam(14) ! Initial temperature
127 hard = uparam(15) ! Linear hardening
128 yld0 = uparam(16) ! Initial yield stress
129 qvoce = uparam(17) ! 1st Voce parameter
130 bvoce = uparam(18) ! 2nd Voce parameter
131 ! Strain-rate dependence parameters
132 jcc = uparam(20) ! Johnson-Cook strain rate coefficient
133 epsp0 = uparam(21) ! Johnson-Cook reference strain rate
134 ! Thermal softening and self-heating parameters
135 mtemp = uparam(22) ! Thermal softening slope
136 tref = uparam(23) ! Reference temperature
137 eta = uparam(24) ! Taylor-Quinney coefficient
138 cp = uparam(25) ! Thermal mass capacity
139 dpis = uparam(26) ! Isothermal plastic strain rate
140 dpad = uparam(27) ! Adiabatic plastic strain rate
141 ! Plastic strain-rate filtering parameters
142 asrate = uparam(28) ! Plastic strain rate filtering frequency
143 afiltr = asrate*timestep/(asrate*timestep + one)
144 dti = one / max(timestep, em20)
145c
146 ! Gurson damage model parameters parameters
147 igurson = nint(uparam(30)) ! Gurson switch flag:
148 ! = 0 => no damage model
149 ! = 1 => local damage model
150 ! = 2 => non local (Forest - micromorphic) damage model
151 ! = 3 => non local (Peerlings) damage model
152 q1 = uparam(31) ! Gurson yield criterion 1st parameter
153 q2 = uparam(32) ! Gurson yield criterion 2nd parameter
154 ed = uparam(34) ! Plastic strain trigger for nucleation
155 an = uparam(35) ! Nucleation rate
156 kw = uparam(36) ! Shear coefficient (Nahshon & Hutchinson)
157 fr = uparam(37) ! Failure void volume fraction
158 fc = uparam(38) ! Critical void volume fraction
159 f0 = uparam(39) ! Initial void volume fraction
160 hkhi = uparam(40) ! Micromorphic penalty parameter
161c
162 ! Recovering internal variables
163 DO i=1,nel
164 ! If the element is failing
165 IF (off(i) < em03) off(i) = zero
166 IF (off(i) < one) off(i) = off(i)*four_over_5
167 ! User inputs
168 phi0(i) = uvar(i,1) ! Previous yield function value
169 yld(i) = uvar(i,2) ! Previous yield stress
170 ! Damage variables
171 fg(i) = dmg(i,2) ! Growth damage
172 fn(i) = dmg(i,3) ! Nucleation damage
173 fsh(i) = dmg(i,4) ! Shear damage
174 ft(i) = dmg(i,5) ! Total damage
175 fs(i) = dmg(i,6) ! Effective damage
176 ! Standard inputs
177 dpla(i) = zero ! Initialization of the plastic strain increment
178 et(i) = one ! Initialization of hourglass coefficient
179 hardp(i) = zero ! Initialization of hardening modulus
180 dezz(i) = zero ! Initialization of the transverse strain increment
181 dlam_nl(i) = zero ! Initialization of the non-local plastic multiplier
182 sig_dfdsig(i) = zero ! Initialization of the stress - normal scalar product
183 normxx(i) = zero ! Initialization of the x component of the normal tensor
184 normyy(i) = zero ! Initialization of the y component of the normal tensor
185 normzz(i) = zero ! Initialization of the z component of the normal tensor
186 ENDDO
187!
188 epsd(1:nel) = uvar(1:nel,3)
189!
190 ! Initialization of damage, temperature and self-heating weight factor
191 IF (time == zero) THEN
192 temp(1:nel) = tini
193 IF (isigi == 0) THEN
194 dmg(1:nel,5) = f0
195 ft(1:nel) = f0
196 dmg(1:nel,1) = f0/fr
197 IF (f0<fc) THEN
198 dmg(1:nel,6) = f0
199 ELSE
200 dmg(1:nel,6) = fc + (one/q1-fc)*(f0-fc)/(fr-fc)
201 ENDIF
202 fs(1:nel) = dmg(1:nel,6)
203 ENDIF
204 ENDIF
205 IF (cp > zero) THEN
206 IF (jthe == 0) THEN
207 DO i=1,nel
208 ! Computation of the weight factor
209 IF (plap_nl(i) < dpis) THEN
210 weitemp(i) = zero
211 ELSEIF (plap_nl(i) > dpad) THEN
212 weitemp(i) = one
213 ELSE
214 weitemp(i) = ((plap_nl(i)-dpis)**2 )
215 . * (three*dpad - two*plap_nl(i) - dpis)
216 . / ((dpad-dpis)**3)
217 ENDIF
218 ENDDO
219 ELSE
220 temp(1:nel) = tempel(1:nel)
221 ENDIF
222 ENDIF
223c
224 !========================================================================
225 ! NON-LOCAL VARIABLES UPDATE
226 !========================================================================
227 DO i=1,nel
228c
229 ! Previous value of Drucker equivalent stress
230 trsig(i) = sigoxx(i) + sigoyy(i)
231 sigm(i) = -trsig(i) * third
232 sxx(i) = sigoxx(i) + sigm(i)
233 syy(i) = sigoyy(i) + sigm(i)
234 szz(i) = sigm(i)
235 sxy(i) = sigoxy(i)
236 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + sxy(i)**2
237 j3(i) = sxx(i) * syy(i) * szz(i) - szz(i) * sxy(i)**2
238 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
239 IF (fdr(i) > zero) THEN
240 sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
241 ELSE
242 sigdr(i) = zero
243 ENDIF
244 sigdr_old(i) = sigdr(i)
245c
246 ! Computation of the stress triaxiality and the etaT factor
247 IF (sigdr(i)>zero) THEN
248 triax(i) = (trsig(i)*third)/sigdr(i)
249 ELSE
250 triax(i) = zero
251 ENDIF
252 IF (trsig(i)<zero) THEN
253 etat(i) = zero
254 ELSE
255 etat(i) = one
256 ENDIF
257c
258 ! Normal to the previous yield surface
259 IF (yld(i)>zero) THEN
260 yld2i = one / yld(i)**2
261 dphi_dsig = two * sigdr(i) * yld2i
262 fsinh = sinh(q2*etat(i)*trsig(i)/(yld(i)*two))
263 dphi_dtrsig(i) = q1*q2*etat(i)*fs(i)*fsinh/yld(i)
264 ELSE
265 yld2i = zero
266 dphi_dsig = zero
267 fsinh = zero
268 dphi_dtrsig(i) = zero
269 ENDIF
270c
271 ! Computation of the Eulerian norm of the stress tensor
272 normsig = sqrt(sigoxx(i)*sigoxx(i)
273 . + sigoyy(i)*sigoyy(i)
274 . + two*sigoxy(i)*sigoxy(i))
275 normsig = max(normsig,one)
276c
277 ! Computation of the normal to the yield surface
278 fdr(i) = (j2(i)/(normsig**2))**3 - cdr*((j3(i)/(normsig**3))**2)
279 IF (fdr(i) > zero) THEN
280 dphi_dfdr = dphi_dsig*kdr*(one/six)*exp(-(five/six)*log(fdr(i)))
281 ELSE
282 dphi_dfdr = zero
283 ENDIF
284 dsdrdj2 = dphi_dfdr*three*(j2(i)/(normsig**2))**2
285 dsdrdj3 = -dphi_dfdr*two*cdr*(j3(i)/(normsig**3))
286 ! dJ3/dS
287 dj3dsxx = two_third*(syy(i)*szz(i))/(normsig**2)
288 . - third*(sxx(i)*szz(i))/(normsig**2)
289 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
290 dj3dsyy = - third*(syy(i)*szz(i))/(normsig**2)
291 . + two_third*(sxx(i)*szz(i))/(normsig**2)
292 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
293 dj3dszz = - third*(syy(i)*szz(i))/(normsig**2)
294 . - third*(sxx(i)*szz(i))/(normsig**2)
295 . + two_third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
296 dj3dsxy = two*(sxx(i)*sxy(i) + sxy(i)*syy(i))/(normsig**2)
297 ! dPhi/dSig
298 normxx(i) = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx + dphi_dtrsig(i)
299 normyy(i) = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy + dphi_dtrsig(i)
300 normzz(i) = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz + dphi_dtrsig(i)
301 normxy(i) = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
302 trdfds(i) = normxx(i) + normyy(i) + normzz(i)
303 sig_dfdsig(i) = sigoxx(i)*normxx(i) + sigoyy(i)*normyy(i)
304 . + sigoxy(i)*normxy(i)
305c
306 ! Computation of the non-local plastic multiplier
307 IF (sig_dfdsig(i) > zero) THEN
308 dlam_nl(i) = ((one - ft(i))*yld(i)*dpla_nl(i))/max(sig_dfdsig(i),em20)
309 ELSE
310 dlam_nl(i) = zero
311 ENDIF
312c
313 ! Damage growth update
314 IF ((trdfds(i)>zero).AND.(ft(i)>zero).AND.(ft(i)<fr)) THEN
315 fg(i) = fg(i) + (one-ft(i))*dlam_nl(i)*trdfds(i)
316 ENDIF
317 fg(i) = max(fg(i),zero)
318c
319 ! Nucleation damage update
320 IF ((pla_nl(i) >= ed).AND.(ft(i)<fr)) THEN
321 ! Case for positive stress triaxiality
322 IF (triax(i)>=zero) THEN
323 fn(i) = fn(i) + an*dpla_nl(i)
324 ! Case for negative stress triaxiality
325 ELSEIF ((triax(i)<zero).AND.(triax(i)>=-third)) THEN
326 fn(i) = fn(i) + an*max(one + three*triax(i),zero)*dpla_nl(i)
327 ENDIF
328 ENDIF
329 fn(i) = max(fn(i),zero)
330c
331 ! Shear damage update
332 IF ((sigdr(i) > zero).AND.(ft(i)>zero).AND.(ft(i)<fr)) THEN
333 sigvm = sqrt(max(em20,three*(j2(i)/(normsig**2))))
334 omega = one - ((twenty7 *(j3(i)/(normsig**3)))/(two*(sigvm**3)))**2
335 omega = max(omega,zero)
336 omega = min(omega,one)
337 sdpla = (sxx(i)*normxx(i)+syy(i)*normyy(i)+ szz(i)*normzz(i)
338 . + sxy(i)*normxy(i))*dlam_nl(i)
339 fsh(i) = fsh(i) + kw*omega*ft(i)*(sdpla/sigdr(i))
340 ENDIF
341 fsh(i) = max(fsh(i),zero)
342c
343 ! Total damage update
344 ft(i) = f0 + fg(i) + fn(i) + fsh(i)
345 ft(i) = min(ft(i), fr)
346 IF (ft(i) >= fr) THEN
347 IF (off(i)==one) off(i) = four_over_5
348 ENDIF
349c
350 ! Effective update
351 IF (ft(i) < fc)THEN
352 fs(i) = ft(i)
353 ELSE
354 fs(i) = fc + (one/q1 - fc) * (ft(i)-fc)/(fr-fc)
355 ENDIF
356 fs(i) = min(fs(i),one/q1)
357c
358 ! Temperature update
359 IF (cp > zero) THEN
360 IF (jthe == 0) THEN
361 dtemp = weitemp(i)*(one-ft(i))*yld(i)*dpla_nl(i)*eta/(rho(i)*cp)
362 temp(i) = temp(i) + dtemp
363 ENDIF
364 ENDIF
365 ENDDO
366c
367 ! Computation of the initial yield stress
368 DO i=1,nel
369 ! a) - Hardening law
370 fhard(i) = yld0 + hard*pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
371 ! b) - Correction factor for strain-rate dependence (Johnson-Cook)
372 frate(i) = one
373 IF (epsd(i) > epsp0) frate(i) = frate(i) + jcc*log(epsd(i)/epsp0)
374 ! c) - Correction factor for thermal effects
375 ftherm(i) = one
376 IF (cp > zero) ftherm(i) = one - mtemp*(temp(i) - tref)
377 ! d) - Computation of the yield function
378 yld(i) = fhard(i)*frate(i)*ftherm(i)
379 ! e) - Checking values
380 yld(i) = max(em10, yld(i))
381 ENDDO
382c
383 !========================================================================
384 ! - COMPUTATION OF TRIAL VALUES
385 !========================================================================
386 DO i=1,nel
387c
388 ! Computation of the trial stress tensor
389 signxx(i) = sigoxx(i) + (a11*depsxx(i) + a12*depsyy(i))
390 signyy(i) = sigoyy(i) + (a11*depsyy(i) + a12*depsxx(i))
391 signxy(i) = sigoxy(i) + (depsxy(i)*g)
392 signyz(i) = sigoyz(i) + (depsyz(i)*gs(i))
393 signzx(i) = sigozx(i) + (depszx(i)*gs(i))
394 ! Computation of the trace of the trial stress tensor
395 trsig(i) = signxx(i) + signyy(i)
396 sigm(i) = -trsig(i) * third
397 ! Computation of the deviatoric trial stress tensor
398 sxx(i) = signxx(i) + sigm(i)
399 syy(i) = signyy(i) + sigm(i)
400 szz(i) = sigm(i)
401 sxy(i) = signxy(i)
402 ! Second deviatoric invariant
403 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + sxy(i)**2
404 ! Third deviatoric invariant
405 j3(i) = sxx(i)*syy(i)*szz(i) - szz(i)*sxy(i)**2
406 ! drucker equivalent stress
407 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
408 ! Checking equivalent stress values
409 IF (fdr(i) > zero) THEN
410 sigdr(i) = kdr * exp((one/six)*log(fdr(i))) ! FDR(I)**(1/6)
411 ELSE
412 sigdr(i) = zero
413 ENDIF
414 ! Computation of the stress triaxiality and the etaT factor
415 IF (sigdr(i)>zero) THEN
416 triax(i) = (trsig(i)*third)/sigdr(i)
417 ELSE
418 triax(i) = zero
419 ENDIF
420 IF (trsig(i)<zero) THEN
421 etat(i) = zero
422 ELSE
423 etat(i) = one
424 ENDIF
425c
426 ! Update out of plane strain increment
427 IF (off(i) == one) THEN
428 dezz(i) = -nnu*depsxx(i)-nnu*depsyy(i)
429 IF (sig_dfdsig(i) > em01) THEN
430 dezz(i) = dezz(i) + nnu*(dlam_nl(i)*normxx(i))
431 . + nnu*(dlam_nl(i)*normyy(i))
432 . + dlam_nl(i)*normzz(i)
433 ENDIF
434 ENDIF
435 ENDDO
436c
437 !========================================================================
438 ! - COMPUTATION OF YIELD FONCTION
439 !========================================================================
440 DO i=1,nel
441 fdam(i) = two*q1*fs(i)*cosh(q2*etat(i)*trsig(i)/yld(i)/two) - (q1*fs(i))**2
442 phi(i) = (sigdr(i) / yld(i))**2 - one + fdam(i)
443 ENDDO
444c
445 ! Checking plastic behavior for all elements
446 nindx = 0
447 DO i=1,nel
448 IF ((phi(i) >= zero).AND.(off(i) == one).AND.(ft(i)<fr)) THEN
449 nindx = nindx+1
450 index(nindx)=i
451 ENDIF
452 ENDDO
453c
454 !====================================================================
455 ! - PLASTIC CORRECTION WITH NICE - EXPLICIT METHOD
456 !====================================================================
457#include "vectorize.inc"
458 DO ii=1,nindx
459c
460 ! Number of the element with plastic behaviour
461 i = index(ii)
462c
463 ! Computation of the trial stress increment
464 dsigxx(i) = a11*depsxx(i) + a12*depsyy(i)
465 dsigyy(i) = a11*depsyy(i) + a12*depsxx(i)
466 dsigxy(i) = depsxy(i)*g
467 dlam = zero
468c
469 ! Note : in this part, the purpose is to compute in one iteration
470 ! a plastic multiplier allowing to update internal variables to satisfy
471 ! the consistency condition.
472 ! Its expression is : DLAMBDA = - (PHI_OLD + DPHI) / DPHI_DLAMBDA
473 ! -> PHI_OLD : old value of yield function (known)
474 ! -> DPHI : yield function prediction (to compute)
475 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
476 ! into account of internal variables kinetic :
477 ! plasticity, temperature and damage (to compute)
478c
479 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
480 !-------------------------------------------------------------
481c
482 ! DPHI_DSIG already computed above
483c
484 ! Restoring previous value of the yield function
485 phi(i) = phi0(i)
486c
487 ! Computation of yield surface trial increment DPHI
488 dphi = normxx(i) * dsigxx(i)
489 . + normyy(i) * dsigyy(i)
490 . + normxy(i) * dsigxy(i)
491c
492 ! 2 - Computation of DPHI_DLAMBDA
493 !---------------------------------------------------------
494c
495 ! a) Derivative with respect stress increments tensor DSIG
496 ! --------------------------------------------------------
497 dfdsig2 = normxx(i) * (a11*normxx(i) + a12*normyy(i))
498 . + normyy(i) * (a11*normyy(i) + a12*normxx(i))
499 . + normxy(i) * normxy(i) * g
500c
501 ! b) Derivatives with respect to plastic strain P
502 ! ------------------------------------------------
503c
504 ! i) Derivative of the yield stress with respect to plastic strain dYLD / dPLA
505 ! ----------------------------------------------------------------------------
506 hardp(i) = hard + qvoce*bvoce*exp(-bvoce*pla(i))
507 dyld_dpla = hardp(i)*frate(i)*ftherm(i)
508c
509 ! ii) Derivative of dPLA with respect to DLAM
510 ! -------------------------------------------
511 dpla_dlam(i) = sig_dfdsig(i) / ((one - ft(i))*yld(i))
512c
513 ! c) Derivative with respect to the yield stress
514 ! ----------------------------------------------
515 sigdr2 = sigdr_old(i)**2
516 trsig(i) = sigoxx(i) + sigoyy(i)
517 dphi_dyld = -two*sigdr2/yld(i)**3-dphi_dtrsig(i)*trsig(i)/yld(i)
518c
519 ! d) Derivative of PHI with respect to DLAM ( = -DENOM)
520 denom = dfdsig2 - (dphi_dyld * dyld_dpla * dpla_dlam(i))
521 denom = sign(max(abs(denom),em20) ,denom)
522c
523 ! 3 - Computation of plastic multiplier and variables update
524 !----------------------------------------------------------
525c
526 ! Computation of the plastic multiplier
527 dlam = (dphi + phi(i)) / denom
528 dlam = max(dlam, zero)
529c
530 ! Plastic strains tensor update
531 dpxx = dlam * normxx(i)
532 dpyy = dlam * normyy(i)
533 dpzz = dlam * normzz(i)
534 dpxy = dlam * normxy(i)
535 trdep(i) = dpxx + dpyy + dpzz
536c
537 ! Elasto-plastic stresses update
538 signxx(i) = sigoxx(i) + dsigxx(i) - (a11*dpxx + a12*dpyy)
539 signyy(i) = sigoyy(i) + dsigyy(i) - (a11*dpyy + a12*dpxx)
540 signxy(i) = sigoxy(i) + dsigxy(i) - dpxy*g
541 trsig(i) = signxx(i) + signyy(i)
542 sigm(i) = -trsig(i) * third
543 sxx(i) = signxx(i) + sigm(i)
544 syy(i) = signyy(i) + sigm(i)
545 szz(i) = sigm(i)
546 sxy(i) = signxy(i)
547c
548 ! Cumulated plastic strain and strain rate update
549 ddep = (dlam/((one - ft(i))*yld(i)))*sig_dfdsig(i)
550 dpla(i) = dpla(i) + ddep
551 pla(i) = pla(i) + ddep
552c
553 ! Drucker equivalent stress update
554 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + sxy(i)**2
555 j3(i) = sxx(i) * syy(i) * szz(i) - szz(i) * sxy(i)**2
556 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
557 sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
558 ! Computation of the stress triaxiality and the etaT factor
559 triax(i) = (trsig(i)*third)/sigdr(i)
560 IF (trsig(i)<zero) THEN
561 etat(i) = zero
562 ELSE
563 etat(i) = one
564 ENDIF
565c
566 ! Hardening law update
567 fhard(i) = yld0 + hard * pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
568c
569 ! Yield stress update
570 yld(i) = fhard(i) * frate(i) * ftherm(i)
571 yld(i) = max(yld(i), em10)
572c
573 ! Yield function value update
574 sigdr2 = sigdr(i)**2
575 yld2i = one / yld(i)**2
576 fcosh = cosh(q2*etat(i)*trsig(i)/(yld(i)*two))
577 fdam(i) = two*q1*fs(i)*fcosh - (q1*fs(i))**2
578 phi(i) = sigdr2 * yld2i - one + fdam(i)
579c
580 ENDDO
581 !===================================================================
582 ! - END OF PLASTIC CORRECTION WITH NICE - EXPLICIT METHOD
583 !===================================================================
584c
585 ! Storing new values
586 DO i=1,nel
587 ! USR Outputs & coefficient for hourglass
588 IF (dpla(i) > zero) THEN
589 uvar(i,1) = phi(i) ! Yield function value
590 et(i) = hardp(i)*frate(i) / (hardp(i)*frate(i) + young)
591 ELSE
592 uvar(i,1) = zero
593 et(i) = one
594 ENDIF
595 uvar(i,2) = yld(i) ! Yield stress
596 ! Standard outputs
597 dmg(i,1) = ft(i)/fr ! Normalized total damage
598 dmg(i,2) = fg(i) ! Growth damage
599 dmg(i,3) = fn(i) ! Nucleation damage
600 dmg(i,4) = fsh(i) ! Shear damage
601 dmg(i,5) = min(ft(i),fr) ! Total damage
602 dmg(i,6) = min(fs(i),one/q1) ! Effective damage
603 seq(i) = sigdr(i) ! Equivalent stress
604 ! If element is broken
605 IF (ft(i) >= fr) THEN
606 dpla(i) = zero
607 signxx(i) = zero
608 signyy(i) = zero
609 signxy(i) = zero
610 signyz(i) = zero
611 signzx(i) = zero
612 seq(i) = zero
613 ENDIF
614 ! Plastic strain-rate (filtered)
615 dpdt = dpla(i) / max(em20,timestep)
616 epsd(i) = afiltr * dpdt + (one - afiltr) * epsd(i)
617 uvar(i,3) = epsd(i)
618 ! Computation of the sound speed
619 soundsp(i) = sqrt((a11)/rho(i))
620 ! Storing the yield stress
621 sigy(i) = yld(i)
622 ! Computation of the thickness variation
623 thk(i) = thk(i) + dezz(i)*thkly(i)*off(i)
624 ENDDO
625c
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