OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat104_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 mat104_nldam_nice (nel, ngl, nuparam, nuvar, volume, fheat, time, timestep, uparam, uvar, jthe, off, rho0, rho, pla, dpla, epsd, soundsp, depsxx, depsyy, depszz, depsxy, depsyz, depszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigy, et, dpla_nl, dmg, temp, seq, pla_nl, plap_nl, jlag)

Function/Subroutine Documentation

◆ mat104_nldam_nice()

subroutine mat104_nldam_nice ( integer nel,
integer, dimension(nel), intent(in) ngl,
integer nuparam,
integer nuvar,
intent(in) volume,
intent(inout) fheat,
time,
timestep,
intent(in) uparam,
intent(inout) uvar,
integer jthe,
intent(inout) off,
intent(in) rho0,
intent(in) rho,
intent(inout) pla,
intent(inout) dpla,
intent(inout) epsd,
intent(out) soundsp,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depszz,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
intent(in) sigoxx,
intent(in) sigoyy,
intent(in) sigozz,
intent(in) sigoxy,
intent(in) sigoyz,
intent(in) sigozx,
intent(out) signxx,
intent(out) signyy,
intent(out) signzz,
intent(out) signxy,
intent(out) signyz,
intent(out) signzx,
intent(out) sigy,
intent(out) et,
intent(in) dpla_nl,
intent(inout) dmg,
intent(inout) temp,
intent(inout) seq,
intent(in) pla_nl,
intent(in) plap_nl,
integer, intent(in) jlag )

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