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,K,II,IGURSON,NINDX,INDEX(NEL),NICE
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,epn,kw,fr,fc,f0,a11,a12,nnu
80 my_real ::
81 . dti,h,ldav,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 . dj2dsxx,dj2dsyy,dj2dsxy,
87 . dfdsxx,dfdsyy,dfdsxy,
88 . denom,dphi,sdpla,sdv_dfdsig,dfdsig2,
89 . dphi_dsig,dphi_dyld,dphi_dfdr,dphi_dfs,dfs_dft,
90 . dfn_dlam,dfsh_dlam,dfg_dlam,dft_dlam,dfdpla,normsig,
91 . dfn,dfsh,dfg,dft,dyld_dpla,dyld_dtemp,dtemp_dlam
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 !=======================================================================
107c
108 !=======================================================================
109 ! - INITIALISATION OF COMPUTATION ON TIME STEP
110 !=======================================================================
111 ! Recovering model parameters
112 ! Elastic parameters
113 young = uparam(1) ! Young modulus
114 bulk = uparam(2) ! Bulk modulus
115 g = uparam(3) ! Shear modulus
116 g2 = uparam(4) ! 2*Shear modulus
117 lam = uparam(5) ! Lambda Hooke parameter
118 nu = uparam(6) ! Poisson ration
119 nnu = uparam(7) ! NU/(1-NU)
120 a11 = uparam(9) ! Diagonal term, elastic matrix in plane stress
121 a12 = uparam(10) ! Non-diagonal term, elastic matrix in plane stress
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 phi0(i) = uvar(i,1) ! Previous yield function value
168 yld(i) = uvar(i,2) ! Previous yield stress
169 ! Damage variables
170 fg(i) = dmg(i,2) ! Growth damage
171 fn(i) = dmg(i,3) ! Nucleation damage
172 fsh(i) = dmg(i,4) ! Shear damage
173 ft(i) = dmg(i,5) ! Total damage
174 fs(i) = dmg(i,6) ! Effective damage
175 ! Standard inputs
176 dpla(i) = zero ! Initialization of the plastic strain increment
177 et(i) = one ! Initialization of hourglass coefficient
178 hardp(i) = zero ! Initialization of hardening modulus
179 dezz(i) = zero ! Initialization of the transverse strain increment
180 dlam_nl(i) = zero ! Initialization of the non-local plastic multiplier
181 sig_dfdsig(i) = zero ! Initialization of the stress - normal scalar product
182 normxx(i) = zero ! Initialization of the x component of the normal tensor
183 normyy(i) = zero ! Initialization of the y component of the normal tensor
184 normzz(i) = zero ! Initialization of the z component of the normal tensor
185 ENDDO
186c
187 ! Initialization of damage, temperature and self-heating weight factor
188 IF (time == zero) THEN
189 temp(1:nel) = tini
190 IF (isigi == 0) THEN
191 dmg(1:nel,5) = f0
192 ft(1:nel) = f0
193 dmg(1:nel,1) = f0/fr
194 IF (f0<fc) THEN
195 dmg(1:nel,6) = f0
196 ELSE
197 dmg(1:nel,6) = fc + (one/q1-fc)*(f0-fc)/(fr-fc)
198 ENDIF
199 fs(1:nel) = dmg(1:nel,6)
200 ENDIF
201 ENDIF
202 IF (cp > zero) THEN
203 IF (jthe == 0) THEN
204 DO i=1,nel
205 ! Computation of the weight factor
206 IF (plap_nl(i) < dpis) THEN
207 weitemp(i) = zero
208 ELSEIF (plap_nl(i) > dpad) THEN
209 weitemp(i) = one
210 ELSE
211 weitemp(i) = ((plap_nl(i)-dpis)**2 )
212 . * (three*dpad - two*plap_nl(i) - dpis)
213 . / ((dpad-dpis)**3)
214 ENDIF
215 ENDDO
216 ELSE
217 temp(1:nel) = tempel(1:nel)
218 ENDIF
219 ENDIF
220c
221 !========================================================================
222 ! NON-LOCAL VARIABLES UPDATE
223 !========================================================================
224 DO i=1,nel
225c
226 ! Previous value of Drucker equivalent stress
227 trsig(i) = sigoxx(i) + sigoyy(i)
228 sigm(i) = -trsig(i) * third
229 sxx(i) = sigoxx(i) + sigm(i)
230 syy(i) = sigoyy(i) + sigm(i)
231 szz(i) = sigm(i)
232 sxy(i) = sigoxy(i)
233 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + sxy(i)**2
234 j3(i) = sxx(i) * syy(i) * szz(i) - szz(i) * sxy(i)**2
235 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
236 IF (fdr(i) > zero) THEN
237 sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
238 ELSE
239 sigdr(i) = zero
240 ENDIF
241 sigdr_old(i) = sigdr(i)
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 = one / yld(i)**2
258 dphi_dsig = two * sigdr(i) * yld2i
259 fsinh = sinh(q2*etat(i)*trsig(i)/(yld(i)*two))
260 dphi_dtrsig(i) = q1*q2*etat(i)*fs(i)*fsinh/yld(i)
261 ELSE
262 yld2i = zero
263 dphi_dsig = zero
264 fsinh = zero
265 dphi_dtrsig(i) = 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)
294 ! dPhi/dSig
295 normxx(i) = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx + dphi_dtrsig(i)
296 normyy(i) = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy + dphi_dtrsig(i)
297 normzz(i) = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz + dphi_dtrsig(i)
298 normxy(i) = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
299 trdfds(i) = normxx(i) + normyy(i) + normzz(i)
300 sig_dfdsig(i) = sigoxx(i)*normxx(i) + sigoyy(i)*normyy(i)
301 . + sigoxy(i)*normxy(i)
302c
303 ! Computation of the non-local plastic multiplier
304 IF (sig_dfdsig(i) > zero) THEN
305 dlam_nl(i) = ((one - ft(i))*yld(i)*dpla_nl(i))/max(sig_dfdsig(i),em20)
306 ELSE
307 dlam_nl(i) = zero
308 ENDIF
309c
310 ! Damage growth update
311 IF ((trdfds(i)>zero).AND.(ft(i)>zero).AND.(ft(i)<fr)) THEN
312 fg(i) = fg(i) + (one-ft(i))*dlam_nl(i)*trdfds(i)
313 ENDIF
314 fg(i) = max(fg(i),zero)
315c
316 ! Nucleation damage update
317 IF ((pla_nl(i) >= ed).AND.(ft(i)<fr)) THEN
318 ! Case for positive stress triaxiality
319 IF (triax(i)>=zero) THEN
320 fn(i) = fn(i) + an*dpla_nl(i)
321 ! Case for negative stress triaxiality
322 ELSEIF ((triax(i)<zero).AND.(triax(i)>=-third)) THEN
323 fn(i) = fn(i) + an*max(one + three*triax(i),zero)*dpla_nl(i)
324 ENDIF
325 ENDIF
326 fn(i) = max(fn(i),zero)
327c
328 ! Shear damage update
329 IF ((sigdr(i) > zero).AND.(ft(i)>zero).AND.(ft(i)<fr)) THEN
330 sigvm = sqrt(max(em20,three*(j2(i)/(normsig**2))))
331 omega = one - ((twenty7 *(j3(i)/(normsig**3)))/(two*(sigvm**3)))**2
332 omega = max(omega,zero)
333 omega = min(omega,one)
334 sdpla = (sxx(i)*normxx(i)+syy(i)*normyy(i)+ szz(i)*normzz(i)
335 . + sxy(i)*normxy(i))*dlam_nl(i)
336 fsh(i) = fsh(i) + kw*omega*ft(i)*(sdpla/sigdr(i))
337 ENDIF
338 fsh(i) = max(fsh(i),zero)
339c
340 ! Total damage update
341 ft(i) = f0 + fg(i) + fn(i) + fsh(i)
342 ft(i) = min(ft(i), fr)
343 IF (ft(i) >= fr) THEN
344 IF (off(i)==one) off(i) = four_over_5
345 ENDIF
346c
347 ! Effective update
348 IF (ft(i) < fc)THEN
349 fs(i) = ft(i)
350 ELSE
351 fs(i) = fc + (one/q1 - fc) * (ft(i)-fc)/(fr-fc)
352 ENDIF
353 fs(i) = min(fs(i),one/q1)
354c
355 ! Temperature update
356 IF (cp > zero) THEN
357 IF (jthe == 0) THEN
358 dtemp = weitemp(i)*(one-ft(i))*yld(i)*dpla_nl(i)*eta/(rho(i)*cp)
359 temp(i) = temp(i) + dtemp
360 ENDIF
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
373 IF (cp > zero) ftherm(i) = one - mtemp*(temp(i) - tref)
374 ! d) - Computation of the yield function
375 yld(i) = fhard(i)*frate(i)*ftherm(i)
376 ! e) - Checking values
377 yld(i) = max(em10, yld(i))
378 ENDDO
379c
380 !========================================================================
381 ! - COMPUTATION OF TRIAL VALUES
382 !========================================================================
383 DO i=1,nel
384c
385 ! Computation of the trial stress tensor
386 signxx(i) = sigoxx(i) + (a11*depsxx(i) + a12*depsyy(i))
387 signyy(i) = sigoyy(i) + (a11*depsyy(i) + a12*depsxx(i))
388 signxy(i) = sigoxy(i) + (depsxy(i)*g)
389 signyz(i) = sigoyz(i) + (depsyz(i)*gs(i))
390 signzx(i) = sigozx(i) + (depszx(i)*gs(i))
391 ! Computation of the trace of the trial stress tensor
392 trsig(i) = signxx(i) + signyy(i)
393 sigm(i) = -trsig(i) * third
394 ! Computation of the deviatoric trial stress tensor
395 sxx(i) = signxx(i) + sigm(i)
396 syy(i) = signyy(i) + sigm(i)
397 szz(i) = sigm(i)
398 sxy(i) = signxy(i)
399 ! Second deviatoric invariant
400 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + sxy(i)**2
401 ! Third deviatoric invariant
402 j3(i) = sxx(i)*syy(i)*szz(i) - szz(i)*sxy(i)**2
403 ! Drucker equivalent stress
404 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
405 ! Checking equivalent stress values
406 IF (fdr(i) > zero) THEN
407 sigdr(i) = kdr * exp((one/six)*log(fdr(i))) ! FDR(I)**(1/6)
408 ELSE
409 sigdr(i) = zero
410 ENDIF
411 ! Computation of the stress triaxiality and the etaT factor
412 IF (sigdr(i)>zero) THEN
413 triax(i) = (trsig(i)*third)/sigdr(i)
414 ELSE
415 triax(i) = zero
416 ENDIF
417 IF (trsig(i)<zero) THEN
418 etat(i) = zero
419 ELSE
420 etat(i) = one
421 ENDIF
422c
423 ! Update out of plane strain increment
424 IF (off(i) == one) THEN
425 dezz(i) = -nnu*depsxx(i)-nnu*depsyy(i)
426 IF (sig_dfdsig(i) > em01) THEN
427 dezz(i) = dezz(i) + nnu*(dlam_nl(i)*normxx(i))
428 . + nnu*(dlam_nl(i)*normyy(i))
429 . + dlam_nl(i)*normzz(i)
430 ENDIF
431 ENDIF
432 ENDDO
433c
434 !========================================================================
435 ! - COMPUTATION OF YIELD FONCTION
436 !========================================================================
437 DO i=1,nel
438 fdam(i) = two*q1*fs(i)*cosh(q2*etat(i)*trsig(i)/yld(i)/two) - (q1*fs(i))**2
439 phi(i) = (sigdr(i) / yld(i))**2 - one + fdam(i)
440 ENDDO
441c
442 ! Checking plastic behavior for all elements
443 nindx = 0
444 DO i=1,nel
445 IF ((phi(i) >= zero).AND.(off(i) == one).AND.(ft(i)<fr)) THEN
446 nindx = nindx+1
447 index(nindx)=i
448 ENDIF
449 ENDDO
450c
451 !====================================================================
452 ! - PLASTIC CORRECTION WITH NICE - EXPLICIT METHOD
453 !====================================================================
454#include "vectorize.inc"
455 DO ii=1,nindx
456c
457 ! Number of the element with plastic behaviour
458 i = index(ii)
459c
460 ! Computation of the trial stress increment
461 dsigxx(i) = a11*depsxx(i) + a12*depsyy(i)
462 dsigyy(i) = a11*depsyy(i) + a12*depsxx(i)
463 dsigxy(i) = depsxy(i)*g
464 dlam = zero
465c
466 ! Note : in this part, the purpose is to compute in one iteration
467 ! a plastic multiplier allowing to update internal variables to satisfy
468 ! the consistency condition.
469 ! Its expression is : DLAMBDA = - (PHI_OLD + DPHI) / DPHI_DLAMBDA
470 ! -> PHI_OLD : old value of yield function (known)
471 ! -> DPHI : yield function prediction (to compute)
472 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
473 ! into account of internal variables kinetic :
474 ! plasticity, temperature and damage (to compute)
475c
476 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
477 !-------------------------------------------------------------
478c
479 ! DPHI_DSIG already computed above
480c
481 ! Restoring previous value of the yield function
482 phi(i) = phi0(i)
483c
484 ! Computation of yield surface trial increment DPHI
485 dphi = normxx(i) * dsigxx(i)
486 . + normyy(i) * dsigyy(i)
487 . + normxy(i) * dsigxy(i)
488c
489 ! 2 - Computation of DPHI_DLAMBDA
490 !---------------------------------------------------------
491c
492 ! a) Derivative with respect stress increments tensor DSIG
493 ! --------------------------------------------------------
494 dfdsig2 = normxx(i) * (a11*normxx(i) + a12*normyy(i))
495 . + normyy(i) * (a11*normyy(i) + a12*normxx(i))
496 . + normxy(i) * normxy(i) * g
497c
498 ! b) Derivatives with respect to plastic strain P
499 ! ------------------------------------------------
500c
501 ! i) Derivative of the yield stress with respect to plastic strain dYLD / dPLA
502 ! ----------------------------------------------------------------------------
503 hardp(i) = hard + qvoce*bvoce*exp(-bvoce*pla(i))
504 dyld_dpla = hardp(i)*frate(i)*ftherm(i)
505c
506 ! ii) Derivative of dPLA with respect to DLAM
507 ! -------------------------------------------
508 dpla_dlam(i) = sig_dfdsig(i) / ((one - ft(i))*yld(i))
509c
510 ! c) Derivative with respect to the yield stress
511 ! ----------------------------------------------
512 sigdr2 = sigdr_old(i)**2
513 trsig(i) = sigoxx(i) + sigoyy(i)
514 dphi_dyld = -two*sigdr2/yld(i)**3-dphi_dtrsig(i)*trsig(i)/yld(i)
515c
516 ! d) Derivative of PHI with respect to DLAM ( = -DENOM)
517 denom = dfdsig2 - (dphi_dyld * dyld_dpla * dpla_dlam(i))
518 denom = sign(max(abs(denom),em20) ,denom)
519c
520 ! 3 - Computation of plastic multiplier and variables update
521 !----------------------------------------------------------
522c
523 ! Computation of the plastic multiplier
524 dlam = (dphi + phi(i)) / denom
525 dlam = max(dlam, zero)
526c
527 ! Plastic strains tensor update
528 dpxx = dlam * normxx(i)
529 dpyy = dlam * normyy(i)
530 dpzz = dlam * normzz(i)
531 dpxy = dlam * normxy(i)
532 trdep(i) = dpxx + dpyy + dpzz
533c
534 ! Elasto-plastic stresses update
535 signxx(i) = sigoxx(i) + dsigxx(i) - (a11*dpxx + a12*dpyy)
536 signyy(i) = sigoyy(i) + dsigyy(i) - (a11*dpyy + a12*dpxx)
537 signxy(i) = sigoxy(i) + dsigxy(i) - dpxy*g
538 trsig(i) = signxx(i) + signyy(i)
539 sigm(i) = -trsig(i) * third
540 sxx(i) = signxx(i) + sigm(i)
541 syy(i) = signyy(i) + sigm(i)
542 szz(i) = sigm(i)
543 sxy(i) = signxy(i)
544c
545 ! Cumulated plastic strain and strain rate update
546 ddep = (dlam/((one - ft(i))*yld(i)))*sig_dfdsig(i)
547 dpla(i) = dpla(i) + ddep
548 pla(i) = pla(i) + ddep
549c
550 ! Drucker equivalent stress update
551 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + sxy(i)**2
552 j3(i) = sxx(i) * syy(i) * szz(i) - szz(i) * sxy(i)**2
553 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
554 sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
555 ! Computation of the stress triaxiality and the etaT factor
556 triax(i) = (trsig(i)*third)/sigdr(i)
557 IF (trsig(i)<zero) THEN
558 etat(i) = zero
559 ELSE
560 etat(i) = one
561 ENDIF
562c
563 ! Hardening law update
564 fhard(i) = yld0 + hard * pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
565c
566 ! Yield stress update
567 yld(i) = fhard(i) * frate(i) * ftherm(i)
568 yld(i) = max(yld(i), em10)
569c
570 ! Yield function value update
571 sigdr2 = sigdr(i)**2
572 yld2i = one / yld(i)**2
573 fcosh = cosh(q2*etat(i)*trsig(i)/(yld(i)*two))
574 fdam(i) = two*q1*fs(i)*fcosh - (q1*fs(i))**2
575 phi(i) = sigdr2 * yld2i - one + fdam(i)
576c
577 ENDDO
578 !===================================================================
579 ! - END OF PLASTIC CORRECTION WITH NICE - EXPLICIT METHOD
580 !===================================================================
581c
582 ! Storing new values
583 DO i=1,nel
584 ! USR Outputs & coefficient for hourglass
585 IF (dpla(i) > zero) THEN
586 uvar(i,1) = phi(i) ! Yield function value
587 et(i) = hardp(i)*frate(i) / (hardp(i)*frate(i) + young)
588 ELSE
589 uvar(i,1) = zero
590 et(i) = one
591 ENDIF
592 uvar(i,2) = yld(i) ! Yield stress
593 ! Standard outputs
594 dmg(i,1) = ft(i)/fr ! Normalized total damage
595 dmg(i,2) = fg(i) ! Growth damage
596 dmg(i,3) = fn(i) ! Nucleation damage
597 dmg(i,4) = fsh(i) ! Shear damage
598 dmg(i,5) = min(ft(i),fr) ! Total damage
599 dmg(i,6) = min(fs(i),one/q1) ! Effective damage
600 seq(i) = sigdr(i) ! Equivalent stress
601 ! If element is broken
602 IF (ft(i) >= fr) THEN
603 dpla(i) = zero
604 signxx(i) = zero
605 signyy(i) = zero
606 signxy(i) = zero
607 signyz(i) = zero
608 signzx(i) = zero
609 seq(i) = zero
610 ENDIF
611 ! Plastic strain-rate (filtered)
612 dpdt = dpla(i) / max(em20,timestep)
613 epsd(i) = afiltr * dpdt + (one - afiltr) * epsd(i)
614 ! Computation of the sound speed
615 soundsp(i) = sqrt((a11)/rho(i))
616 ! Storing the yield stress
617 sigy(i) = yld(i)
618 ! Computation of the thickness variation
619 thk(i) = thk(i) + dezz(i)*thkly(i)*off(i)
620 ENDDO
621c
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