OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat104_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 mat104_nldam_newton (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_newton()

subroutine mat104_nldam_newton ( 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_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 ,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,II,IGURSON,ITER,NITER,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,hkhi,
81 . q1,q2,ed,an,epn,kw,fr,fc,f0
82 my_real ::
83 . dti,h,ldav,trdep,sigvm,omega,
84 . dtemp,fcosh,fsinh,dpdt,dlam,ddep
85 my_real ::
86 . dsdrdj2,dsdrdj3,
87 . dj3dsxx,dj3dsyy,dj3dszz,dj3dsxy,dj3dsyz,dj3dszx,
88 . dj2dsxx,dj2dsyy,dj2dszz,dj2dsxy,dj2dsyz,dj2dszx,
89 . dfdsxx,dfdsyy,dfdszz,dfdsxy,dfdsyz,dfdszx,
90 . normsig,sdpla,dphi_dtrsig,dfdsig2,sdv_dfdsig,
91 . dphi_dsig,dphi_dyld,dphi_dfdr,df_dfs,dfs_dft,dphi_dft,
92 . dphi_dfs,dfn_dlam,dfsh_dlam,dfg_dlam,dft_dlam,
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,
97 . sxx,syy,szz,sxy,syz,szx,sigm,j2,j3,sigdr,yld,weitemp,
98 . hardp,fhard,frate,ftherm,dtherm,fdr,phi0,trdfds,triax,
99 . fdam,phi,ft,fs,fg,fn,fsh,dpla_dlam,dphi_dlam,etat,
100 . normxx,normyy,normxy,normzz,normyz,normzx,sig_dfdsig,
101 . dpxx,dpyy,dpxy,dpyz,dpzx,dpzz,sigdr2,yld2i,dlam_nl
102c
103 !=======================================================================
104 ! DRUCKER - VOCE - JOHNSON-COOK MATERIAL LAW WITH GURSON DAMAGE
105 ! USING NON LOCAL PEERLINGS METHOD
106 !=======================================================================
107 !UVAR(1) 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
157 hkhi = uparam(40) ! Micromorphic penalty parameter
158c
159 ! Recovering internal variables
160 DO i=1,nel
161 ! If the element is failing
162 IF (off(i) < em03) off(i) = zero
163 IF (off(i) < one) off(i) = off(i)*four_over_5
164 ! User inputs
165 yld(i) = uvar(i,1) ! 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 ((jthe == 0 .AND. 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
236c
237 ! Computation of the stress triaxiality and the etaT factor
238 IF (sigdr(i)>zero) THEN
239 triax(i) = (trsig(i)*third)/sigdr(i)
240 ELSE
241 triax(i) = zero
242 ENDIF
243 IF (trsig(i)<zero) THEN
244 etat(i) = zero
245 ELSE
246 etat(i) = one
247 ENDIF
248c
249 ! Normal to the previous yield surface
250 IF (yld(i)>zero) THEN
251 yld2i(i) = one / yld(i)**2
252 dphi_dsig = two * sigdr(i) * yld2i(i)
253 fsinh = sinh(q2*etat(i)*trsig(i)/(yld(i)*two))
254 dphi_dtrsig = q1*q2*etat(i)*fs(i)*fsinh/yld(i)
255 ELSE
256 yld2i(i) = zero
257 dphi_dsig = zero
258 fsinh = zero
259 dphi_dtrsig = zero
260 ENDIF
261c
262 ! Computation of the Eulerian norm of the stress tensor
263 normsig = sqrt(sigoxx(i)*sigoxx(i)
264 . + sigoyy(i)*sigoyy(i)
265 . + sigozz(i)*sigozz(i)
266 . + two*sigoxy(i)*sigoxy(i)
267 . + two*sigoyz(i)*sigoyz(i)
268 . + two*sigozx(i)*sigozx(i))
269 normsig = max(normsig,one)
270c
271 ! Computation of the normal to the yield surface
272 fdr(i) = (j2(i)/(normsig**2))**3 - cdr*((j3(i)/(normsig**3))**2)
273 IF (fdr(i) > zero) THEN
274 dphi_dfdr = dphi_dsig*kdr*(one/six)*exp(-(five/six)*log(fdr(i)))
275 ELSE
276 dphi_dfdr = zero
277 ENDIF
278 dsdrdj2 = dphi_dfdr*three*(j2(i)/(normsig**2))**2
279 dsdrdj3 = -dphi_dfdr*two*cdr*(j3(i)/(normsig**3))
280 !dJ3/dS
281 dj3dsxx = two_third*(syy(i)*szz(i)-syz(i)**2)/(normsig**2)
282 . - third*(sxx(i)*szz(i)-szx(i)**2)/(normsig**2)
283 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
284 dj3dsyy = - third*(syy(i)*szz(i)-syz(i)**2)/(normsig**2)
285 . + two_third*(sxx(i)*szz(i)-szx(i)**2)/(normsig**2)
286 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
287 dj3dszz = - third*(syy(i)*szz(i)-syz(i)**2)/(normsig**2)
288 . - third*(sxx(i)*szz(i)-szx(i)**2)/(normsig**2)
289 . + two_third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
290 dj3dsxy = two*(sxx(i)*sxy(i) + sxy(i)*syy(i) + szx(i)*syz(i))/(normsig**2)
291 dj3dsyz = two*(sxy(i)*szx(i) + syy(i)*syz(i) + syz(i)*szz(i))/(normsig**2)
292 dj3dszx = two*(sxx(i)*szx(i) + sxy(i)*syz(i) + szx(i)*szz(i))/(normsig**2)
293 ! dPhi/dSig
294 normxx(i) = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx + dphi_dtrsig
295 normyy(i) = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy + dphi_dtrsig
296 normzz(i) = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz + dphi_dtrsig
297 normxy(i) = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
298 normyz(i) = two*dsdrdj2*syz(i)/normsig + dsdrdj3*dj3dsyz
299 normzx(i) = two*dsdrdj2*szx(i)/normsig + dsdrdj3*dj3dszx
300 trdfds(i) = normxx(i) + normyy(i) + normzz(i)
301 sig_dfdsig(i) = sigoxx(i)*normxx(i) + sigoyy(i)*normyy(i) + sigozz(i)*normzz(i)
302 . + sigoxy(i)*normxy(i) + sigoyz(i)*normyz(i) + sigozx(i)*normzx(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))/sig_dfdsig(i)
307 ELSE
308 dlam_nl(i) = zero
309 ENDIF
310c
311 ! Damage growth update
312 IF ((ft(i)>zero).AND.(ft(i)<fr).AND.(trdfds(i)>zero)) 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)+ syz(i)*normyz(i)+ szx(i)*normzx(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 (jthe == 0 .AND. cp > zero ) THEN
359 dtemp = weitemp(i)*(one-ft(i))*yld(i)*dpla_nl(i)*eta/(rho0(i)*cp)
360 temp(i) = temp(i) + dtemp
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 - mtemp*(temp(i) - tref)
373 ! d) - Computation of the yield function
374 yld(i) = fhard(i)*frate(i)*ftherm(i)
375 ! e) - Checking values
376 yld(i) = max(em10, yld(i))
377 ENDDO
378c
379 !========================================================================
380 ! - COMPUTATION OF TRIAL VALUES
381 !========================================================================
382 DO i=1,nel
383c
384 ! Computation of the trial stress tensor
385 ldav = (depsxx(i) + depsyy(i) + depszz(i)) * lam
386 signxx(i) = sigoxx(i) + (depsxx(i)*g2 + ldav)
387 signyy(i) = sigoyy(i) + (depsyy(i)*g2 + ldav)
388 signzz(i) = sigozz(i) + (depszz(i)*g2 + ldav)
389 signxy(i) = sigoxy(i) + depsxy(i)*g
390 signyz(i) = sigoyz(i) + depsyz(i)*g
391 signzx(i) = sigozx(i) + depszx(i)*g
392 ! Computation of the trace of the trial stress tensor
393 trsig(i) = signxx(i) + signyy(i) + signzz(i)
394 sigm(i) = -trsig(i) * third
395 ! Computation of the deviatoric trial stress tensor
396 sxx(i) = signxx(i) + sigm(i)
397 syy(i) = signyy(i) + sigm(i)
398 szz(i) = signzz(i) + sigm(i)
399 sxy(i) = signxy(i)
400 syz(i) = signyz(i)
401 szx(i) = signzx(i)
402 ! Second deviatoric invariant
403 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 )
404 . + sxy(i)**2 + syz(i)**2 + szx(i)**2
405 ! Third deviatoric invariant
406 j3(i) = sxx(i) * syy(i) * szz(i)
407 . + sxy(i) * syz(i) * szx(i) * two
408 . - sxx(i) * syz(i)**2
409 . - syy(i) * szx(i)**2
410 . - szz(i) * sxy(i)**2
411 ! Drucker equivalent stress
412 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
413 ! Checking equivalent stress values
414 IF (fdr(i) > zero) THEN
415 sigdr(i) = kdr * exp((one/six)*log(fdr(i))) ! FDR(I)**(1/6)
416 ELSE
417 sigdr(i) = zero
418 ENDIF
419 ! computation of the stress triaxiality and the etat factor
420 IF (sigdr(i)>zero) THEN
421 triax(i) = (trsig(i)*third)/sigdr(i)
422 ELSE
423 triax(i) = zero
424 ENDIF
425 IF (trsig(i)<zero) THEN
426 etat(i) = zero
427 ELSE
428 etat(i) = one
429 ENDIF
430 ENDDO
431c
432 !========================================================================
433 ! - COMPUTATION OF YIELD FONCTION
434 !=======================================================================
435 DO i=1,nel
436 fdam(i) = two*q1*fs(i)*cosh(q2*etat(i)*trsig(i)/yld(i)/two) - (q1*fs(i))**2
437 phi(i) = (sigdr(i) / yld(i))**2 - one + fdam(i)
438 ENDDO
439c
440 ! Checking plastic behavior for all elements
441 nindx = 0
442 DO i=1,nel
443 IF ((phi(i) >= zero).AND.(off(i) == one).AND.(ft(i) < fr)) THEN
444 nindx=nindx+1
445 index(nindx)=i
446 ENDIF
447 ENDDO
448c
449 !====================================================================
450 ! - PLASTIC CORRECTION WITH CUTTING PLANE METHOD (SEMI-IMPLICIT)
451 !====================================================================
452c
453 ! Number of iterations
454 niter = 3
455c
456 ! Loop over plastifying elements
457#include "vectorize.inc"
458 DO ii=1,nindx
459c
460 ! Number of the element with plastic behaviour
461 i = index(ii)
462c
463 ! Initialization of the iterative Newton procedure
464 sigdr2(i) = sigdr(i)**2
465 yld2i(i) = one / yld(i)**2
466 dpxx(i) = zero
467 dpyy(i) = zero
468 dpzz(i) = zero
469 dpxy(i) = zero
470 dpyz(i) = zero
471 dpzx(i) = zero
472 ENDDO
473c
474 ! Loop over the iterations
475 DO iter = 1, niter
476#include "vectorize.inc"
477 DO ii=1,nindx
478 i = index(ii)
479c
480 ! Note: in this part, the purpose is to compute for each iteration
481 ! a plastic multiplier allowing to update internal variables to satisfy
482 ! the consistency condition using the cutting plane semi-implicit
483 ! iterative procedure.
484 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
485 ! -> PHI : current value of yield function (known)
486 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
487 ! into account of internal variables kinetic :
488 ! plasticity, temperature and damage (to compute)
489c
490 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
491 !-------------------------------------------------------------
492c
493 ! Derivative with respect to the equivalent stress and trace
494 yld2i(i) = one/(yld(i)**2)
495 dphi_dsig = two*sigdr(i)*yld2i(i)
496 fsinh = sinh(q2*etat(i)*trsig(i)/(yld(i)*two))
497 dphi_dtrsig = q1*q2*etat(i)*fs(i)*fsinh/yld(i)
498c
499 ! Computation of the Eulerian norm of the stress tensor
500 normsig = sqrt(signxx(i)*signxx(i)
501 . + signyy(i)*signyy(i)
502 . + signzz(i)*signzz(i)
503 . + two*signxy(i)*signxy(i)
504 . + two*signyz(i)*signyz(i)
505 . + two*signzx(i)*signzx(i))
506 normsig = max(normsig,one)
507c
508 ! Derivative with respect to Fdr
509 fdr(i) = (j2(i)/(normsig**2))**3 - cdr*((j3(i)/(normsig**3))**2)
510 dphi_dfdr = dphi_dsig*kdr*(one/six)*exp(-(five/six)*log(fdr(i)))
511 dsdrdj2 = dphi_dfdr*three*(j2(i)/(normsig**2))**2
512 dsdrdj3 = -dphi_dfdr*two*cdr*(j3(i)/(normsig**3))
513 ! dJ3/dSig
514 dj3dsxx = two_third*(syy(i)*szz(i)-syz(i)**2)/(normsig**2)
515 . - third*(sxx(i)*szz(i)-szx(i)**2)/(normsig**2)
516 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
517 dj3dsyy = - third*(syy(i)*szz(i)-syz(i)**2)/(normsig**2)
518 . + two_third*(sxx(i)*szz(i)-szx(i)**2)/(normsig**2)
519 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
520 dj3dszz = - third*(syy(i)*szz(i)-syz(i)**2)/(normsig**2)
521 . - third*(sxx(i)*szz(i)-szx(i)**2)/(normsig**2)
522 . + two_third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
523 dj3dsxy = two*(sxx(i)*sxy(i) + sxy(i)*syy(i) + szx(i)*syz(i))/(normsig**2)
524 dj3dsyz = two*(sxy(i)*szx(i) + syy(i)*syz(i) + syz(i)*szz(i))/(normsig**2)
525 dj3dszx = two*(sxx(i)*szx(i) + sxy(i)*syz(i) + szx(i)*szz(i))/(normsig**2)
526 ! dPhi/dSig
527 normxx(i) = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx + dphi_dtrsig
528 normyy(i) = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy + dphi_dtrsig
529 normzz(i) = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz + dphi_dtrsig
530 normxy(i) = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
531 normyz(i) = two*dsdrdj2*syz(i)/normsig + dsdrdj3*dj3dsyz
532 normzx(i) = two*dsdrdj2*szx(i)/normsig + dsdrdj3*dj3dszx
533c
534 ! 2 - Computation of DPHI_DLAMBDA
535 !---------------------------------------------------------
536c
537 ! a) Derivative with respect stress increments tensor DSIG
538 ! --------------------------------------------------------
539 trdfds(i) = normxx(i) + normyy(i) + normzz(i)
540 dfdsig2 = normxx(i) * (normxx(i)*g2 + lam*trdfds(i))
541 . + normyy(i) * (normyy(i)*g2 + lam*trdfds(i))
542 . + normzz(i) * (normzz(i)*g2 + lam*trdfds(i))
543 . + normxy(i) * normxy(i) * g
544 . + normyz(i) * normyz(i) * g
545 . + normzx(i) * normzx(i) * g
546c
547 ! b) Derivatives with respect to plastic strain P
548 ! ------------------------------------------------
549c
550 ! i) Derivative of the yield stress with respect to plastic strain dYLD / dPLA
551 ! ----------------------------------------------------------------------------
552 hardp(i) = hard + qvoce*bvoce*exp(-bvoce*pla(i))
553 dyld_dpla = hardp(i)*frate(i)*ftherm(i)
554c
555 ! ii) Derivative of dPLA with respect to DLAM
556 ! -------------------------------------------
557 sig_dfdsig(i) = signxx(i) * normxx(i)
558 . + signyy(i) * normyy(i)
559 . + signzz(i) * normzz(i)
560 . + signxy(i) * normxy(i)
561 . + signyz(i) * normyz(i)
562 . + signzx(i) * normzx(i)
563 dpla_dlam(i) = sig_dfdsig(i) / ((one - ft(i))*yld(i))
564c
565 ! c) Derivative with respect to the yield stress
566 ! ----------------------------------------------
567 sigdr2(i) = sigdr(i)**2
568 dphi_dyld = -two*sigdr2(i)/yld(i)**3-dphi_dtrsig*trsig(i)/yld(i)
569c
570 ! d) Derivative of PHI with respect to DLAM
571 dphi_dlam(i) = - dfdsig2 + dphi_dyld*dyld_dpla*dpla_dlam(i)
572 dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20) ,dphi_dlam(i))
573c
574 ! 3 - Computation of plastic multiplier and variables update
575 !----------------------------------------------------------
576c
577 ! Computation of the plastic multiplier
578 dlam = -phi(i)/dphi_dlam(i)
579c
580 ! Plastic strains tensor update
581 dpxx(i) = dlam * normxx(i)
582 dpyy(i) = dlam * normyy(i)
583 dpzz(i) = dlam * normzz(i)
584 dpxy(i) = dlam * normxy(i)
585 dpyz(i) = dlam * normyz(i)
586 dpzx(i) = dlam * normzx(i)
587 trdep = dpxx(i) + dpyy(i) + dpzz(i)
588c
589 ! Elasto-plastic stresses update
590 ldav = trdep * lam
591 signxx(i) = signxx(i) - (dpxx(i)*g2 + ldav)
592 signyy(i) = signyy(i) - (dpyy(i)*g2 + ldav)
593 signzz(i) = signzz(i) - (dpzz(i)*g2 + ldav)
594 signxy(i) = signxy(i) - dpxy(i)*g
595 signyz(i) = signyz(i) - dpyz(i)*g
596 signzx(i) = signzx(i) - dpzx(i)*g
597 trsig(i) = signxx(i) + signyy(i) + signzz(i)
598 sigm(i) = -trsig(i) * third
599 sxx(i) = signxx(i) + sigm(i)
600 syy(i) = signyy(i) + sigm(i)
601 szz(i) = signzz(i) + sigm(i)
602 sxy(i) = signxy(i)
603 syz(i) = signyz(i)
604 szx(i) = signzx(i)
605c
606 ! Cumulated plastic strain and strain rate update
607 ddep = (dlam/((one - ft(i))*yld(i)))*sig_dfdsig(i)
608 dpla(i) = dpla(i) + ddep
609 pla(i) = pla(i) + ddep
610c
611 ! Drucker equivalent stress update
612 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 )
613 . + sxy(i)**2 + syz(i)**2 + szx(i)**2
614 j3(i) = sxx(i) * syy(i) * szz(i) + two * sxy(i) * syz(i) * szx(i)
615 . - sxx(i) * syz(i)**2 - syy(i) * szx(i)**2 - szz(i) * sxy(i)**2
616 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
617 sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
618 ! Computation of the stress triaxiality and the etaT factor
619 triax(i) = (trsig(i)*third)/sigdr(i)
620 IF (trsig(i)<zero) THEN
621 etat(i) = zero
622 ELSE
623 etat(i) = one
624 ENDIF
625c
626 ! Hardening law update
627 fhard(i) = yld0 + hard * pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
628c
629 ! Yield stress update
630 yld(i) = fhard(i) * frate(i) * ftherm(i)
631 yld(i) = max(yld(i), em10)
632c
633 ! Yield function value update
634 sigdr2(i) = sigdr(i)**2
635 yld2i(i) = one / yld(i)**2
636 fcosh = cosh(q2*etat(i)*trsig(i)/(yld(i)*two))
637 fdam(i) = two*q1*fs(i)*fcosh - (q1*fs(i))**2
638 phi(i) = sigdr2(i) * yld2i(i) - one + fdam(i)
639c
640 ENDDO
641 ! End of the loop over the iterations
642 ENDDO
643 !===================================================================
644 ! - END OF PLASTIC CORRECTION WITH CUTTING PLANE ITERATIVE METHOD
645 !===================================================================
646c
647 ! Storing new values
648 DO i=1,nel
649 ! USR Outputs
650 uvar(i,1) = yld(i) ! Yield stress
651 ! Standard outputs
652 dmg(i,1) = ft(i)/fr ! Normalized total damage
653 dmg(i,2) = fg(i) ! Growth damage
654 dmg(i,3) = fn(i) ! Nucleation damage
655 dmg(i,4) = fsh(i) ! Shear damage
656 dmg(i,5) = min(ft(i),fr) ! Total damage
657 dmg(i,6) = min(fs(i),one/q1) ! Effective damage
658 seq(i) = sigdr(i) ! Equivalent stress
659 ! If element is broken
660 IF (ft(i) >= fr) THEN
661 dpla(i) = zero
662 signxx(i) = zero
663 signyy(i) = zero
664 signzz(i) = zero
665 signxy(i) = zero
666 signyz(i) = zero
667 signzx(i) = zero
668 seq(i) = zero
669 ENDIF
670 ! Plastic strain-rate (filtered)
671 dpdt = dpla(i) / max(em20,timestep)
672 epsd(i) = afiltr * dpdt + (one - afiltr) * epsd(i)
673 ! Coefficient for hourglass
674 IF (dpla(i) > zero) THEN
675 et(i) = hardp(i)*frate(i) / (hardp(i)*frate(i) + young)
676 ELSE
677 et(i) = one
678 ENDIF
679 ! Computation of the sound speed
680 soundsp(i) = sqrt((bulk + four_over_3*g)/rho0(i))
681 ! Storing the yield stress
682 sigy(i) = yld(i)
683 ENDDO
684C-----------------------------------------------
685C plastic work dissipated as heat - in case of /heat/mat
686C-----------------------------------------------
687 IF (jthe < 0 .AND. jlag /= 0) THEN
688 DO i=1,nel
689 fheat(i) = fheat(i) + eta*(one-ft(i))*weitemp(i)*yld(i)*dpla(i)*volume(i)
690 ENDDO
691 END IF
692!-----------
my_real
#define my_real
Definition cppsort.cpp:32
the
end diagonal values have been computed in the(sparse) matrix id.SOL
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21