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

Go to the source code of this file.

Functions/Subroutines

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

Function/Subroutine Documentation

◆ mat104c_nldam_newton()

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

Definition at line 28 of file mat104c_nldam_newton.F.

37 !=======================================================================
38 ! Implicit types
39 !=======================================================================
40#include "implicit_f.inc"
41 !=======================================================================
42 ! Common
43 !=======================================================================
44#include "com01_c.inc"
45 !=======================================================================
46 ! Dummy arguments
47 !=======================================================================
48 INTEGER NEL,NUPARAM,NUVAR,JTHE
49 INTEGER ,DIMENSION(NEL), INTENT(IN) :: NGL
50 my_real
51 . time,timestep
52 my_real,DIMENSION(NUPARAM), INTENT(IN) ::
53 . uparam
54 my_real,DIMENSION(NEL), INTENT(IN) ::
55 . rho,tempel,
56 . depsxx,depsyy,depsxy,depsyz,depszx,
57 . sigoxx,sigoyy,sigoxy,sigoyz,sigozx,
58 . gs,thkly,pla_nl,plap_nl,dpla_nl
59c
60 my_real ,DIMENSION(NEL), INTENT(OUT) ::
61 . soundsp,sigy,et,
62 . signxx,signyy,signxy,signyz,signzx
63c
64 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
65 . pla,dpla,epsd,off,thk,temp,seq
66 my_real ,DIMENSION(NEL,6), INTENT(INOUT) ::
67 . dmg
68 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) ::
69 . uvar
70 !=======================================================================
71 ! Local Variables
72 !=======================================================================
73 INTEGER I,II,IGURSON,ITER,NITER,NINDX,INDEX(NEL)
74c
75 my_real ::
76 . young,bulk,lam,g,g2,nu,cdr,kdr,hard,yld0,qvoce,bvoce,jcc,
77 . epsp0,mtemp,tini,tref,eta,cp,dpis,dpad,asrate,afiltr,hkhi,
78 . q1,q2,ed,an,epn,kw,fr,fc,f0,a11,a12,nnu
79 my_real ::
80 . dti,h,ldav,sigvm,omega,
81 . dtemp,fcosh,fsinh,dpdt,dlam,ddep
82 my_real ::
83 . dsdrdj2,dsdrdj3,
84 . dj3dsxx,dj3dsyy,dj3dsxy,dj3dszz,
85 . dj2dsxx,dj2dsyy,dj2dsxy,normsig,
86 . dfdsxx,dfdsyy,dfdsxy,
87 . sdpla,dphi_dtrsig,dfdsig2,sdv_dfdsig,
88 . dphi_dsig,dphi_dyld,dphi_dfdr,df_dfs,dfs_dft,dphi_dft,
89 . dphi_dfs,dfn_dlam,dfsh_dlam,dfg_dlam,dft_dlam,
90 . dfn,dfsh,dfg,dft,dyld_dpla,dyld_dtemp,dtemp_dlam
91c
92 my_real, DIMENSION(NEL) ::
93 . dsigxx,dsigyy,dsigxy,trsig,trdep,trdfds,
94 . sxx,syy,sxy,szz,sigm,j2,j3,sigdr,yld,weitemp,
95 . hardp,fhard,frate,ftherm,dtherm,fdr,phi0,triax,
96 . fdam,phi,ft,fs,fg,fn,fsh,dpla_dlam,dphi_dlam,dezz,etat,
97 . normxx,normyy,normxy,normzz,sig_dfdsig,
98 . dpxx,dpyy,dpxy,dpzz,sigdr2,yld2i,dlam_nl
99c
100 !=======================================================================
101 ! DRUCKER - VOCE - JOHNSON-COOK MATERIAL LAW WITH GURSON DAMAGE
102 ! USING NON LOCAL PEERLINGS METHOD
103 !=======================================================================
104 !UVAR(1) YLD YIELD STRESS
105 !=======================================================================
106c
107 !=======================================================================
108 ! - INITIALISATION OF COMPUTATION ON TIME STEP
109 !=======================================================================
110 ! Recovering model parameters
111 ! Elastic parameters
112 young = uparam(1)
113 bulk = uparam(2)
114 g = uparam(3)
115 g2 = uparam(4)
116 lam = uparam(5)
117 nu = uparam(6)
118 nnu = uparam(7)
119 a11 = uparam(9)
120 a12 = uparam(10)
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)
143c
144 ! Gurson damage model parameters parameters
145 igurson = nint(uparam(30)) ! Gurson switch flag:
146 ! = 0 => no damage model
147 ! = 1 => local damage model
148 ! = 2 => non local (Forest - micromorphic) damage model
149 ! = 3 => non local (Peerlings) damage model
150 q1 = uparam(31) ! Gurson yield criterion 1st parameter
151 q2 = uparam(32) ! Gurson yield criterion 2nd parameter
152 ed = uparam(34) ! Plastic strain trigger for nucleation
153 an = uparam(35) ! Nucleation rate
154 kw = uparam(36) ! Shear coefficient (Nahshon & Hutchinson)
155 fr = uparam(37) ! Failure void volume fraction
156 fc = uparam(38) ! Critical void volume fraction
157 f0 = uparam(39) ! Initial void volume fraction
158 hkhi = uparam(40) ! Micromorphic penalty parameter
159c
160 ! Recovering internal variables
161 DO i=1,nel
162 ! If the element is failing
163 IF (off(i) < em03) off(i) = zero
164 IF (off(i) < one) off(i) = off(i)*four_over_5
165 ! User inputs
166 yld(i) = uvar(i,1) ! Previous yield stress
167 ! Damage variables
168 fg(i) = dmg(i,2) ! Growth damage
169 fn(i) = dmg(i,3) ! Nucleation damage
170 fsh(i) = dmg(i,4) ! Shear damage
171 ft(i) = dmg(i,5) ! Total damage
172 fs(i) = dmg(i,6) ! Effective damage
173 ! Standard inputs
174 dpla(i) = zero ! Initialization of the plastic strain increment
175 et(i) = one ! Initialization of hourglass coefficient
176 hardp(i) = zero ! Initialization of hardening modulus
177 dezz(i) = zero ! Initialization of the transverse strain
178 dlam_nl(i) = zero ! Initialization of the non-local plastic multiplier
179 sig_dfdsig(i) = zero ! Initialization of the stress - normal scalar product
180 normxx(i) = zero ! Initialization of the x component of the normal tensor
181 normyy(i) = zero ! Initialization of the y component of the normal tensor
182 normzz(i) = zero ! Initialization of the z component of the normal tensor
183 ENDDO
184c
185 ! Initialization of damage, temperature and self-heating weight factor
186 IF (time == zero) THEN
187 temp(1:nel) = tini
188 IF (isigi == 0) THEN
189 dmg(1:nel,5) = f0
190 ft(1:nel) = f0
191 dmg(1:nel,1) = f0/fr
192 IF (f0<fc) THEN
193 dmg(1:nel,6) = f0
194 ELSE
195 dmg(1:nel,6) = fc + (one/q1-fc)*(f0-fc)/(fr-fc)
196 ENDIF
197 fs(1:nel) = dmg(1:nel,6)
198 ENDIF
199 ENDIF
200 IF (cp > zero) THEN
201 IF (jthe == 0) THEN
202 DO i=1,nel
203 ! Computation of the weight factor
204 IF (plap_nl(i) < dpis) THEN
205 weitemp(i) = zero
206 ELSEIF (plap_nl(i) > dpad) THEN
207 weitemp(i) = one
208 ELSE
209 weitemp(i) = ((plap_nl(i)-dpis)**2 )
210 . * (three*dpad - two*plap_nl(i) - dpis)
211 . / ((dpad-dpis)**3)
212 ENDIF
213 ENDDO
214 ELSE
215 temp(1:nel) = tempel(1:nel)
216 ENDIF
217 ENDIF
218c
219 !========================================================================
220 ! NON-LOCAL VARIABLES UPDATE
221 !========================================================================
222 DO i=1,nel
223c
224 ! Previous value of Drucker equivalent stress
225 trsig(i) = sigoxx(i) + sigoyy(i)
226 sigm(i) = -trsig(i) * third
227 sxx(i) = sigoxx(i) + sigm(i)
228 syy(i) = sigoyy(i) + sigm(i)
229 szz(i) = sigm(i)
230 sxy(i) = sigoxy(i)
231 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) + sxy(i)**2
232 j3(i) = sxx(i) * syy(i) * szz(i) - szz(i) * sxy(i)**2
233 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
234 IF (fdr(i) > zero) THEN
235 sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
236 ELSE
237 sigdr(i) = zero
238 ENDIF
239c
240 ! Computation of the stress triaxiality and the etaT factor
241 IF (sigdr(i)>zero) THEN
242 triax(i) = (trsig(i)*third)/sigdr(i)
243 ELSE
244 triax(i) = zero
245 ENDIF
246 IF (trsig(i)<zero) THEN
247 etat(i) = zero
248 ELSE
249 etat(i) = one
250 ENDIF
251c
252 ! Normal to the previous yield surface
253 IF (yld(i)>zero) THEN
254 yld2i(i) = one / yld(i)**2
255 dphi_dsig = two * sigdr(i) * yld2i(i)
256 fsinh = sinh(q2*etat(i)*trsig(i)/(yld(i)*two))
257 dphi_dtrsig = q1*q2*etat(i)*fs(i)*fsinh/yld(i)
258 ELSE
259 yld2i(i) = zero
260 dphi_dsig = zero
261 fsinh = zero
262 dphi_dtrsig = zero
263 ENDIF
264c
265 ! Computation of the Eulerian norm of the stress tensor
266 normsig = sqrt(sigoxx(i)*sigoxx(i)
267 . + sigoyy(i)*sigoyy(i)
268 . + two*sigoxy(i)*sigoxy(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))/(normsig**2)
282 . - third*(sxx(i)*szz(i))/(normsig**2)
283 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
284 dj3dsyy = - third*(syy(i)*szz(i))/(normsig**2)
285 . + two_third*(sxx(i)*szz(i))/(normsig**2)
286 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
287 dj3dszz = - third*(syy(i)*szz(i))/(normsig**2)
288 . - third*(sxx(i)*szz(i))/(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))/(normsig**2)
291c
292 ! dPhi/dSig
293 normxx(i) = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx + dphi_dtrsig
294 normyy(i) = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy + dphi_dtrsig
295 normzz(i) = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz + dphi_dtrsig
296 normxy(i) = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
297 trdfds(i) = normxx(i) + normyy(i) + normzz(i)
298 sig_dfdsig(i) = sigoxx(i)*normxx(i) + sigoyy(i)*normyy(i)
299 . + sigoxy(i)*normxy(i)
300c
301 ! Computation of the non-local plastic multiplier
302 IF (sig_dfdsig(i)>zero) THEN
303 dlam_nl(i) = ((one - ft(i))*yld(i)*dpla_nl(i))/max(sig_dfdsig(i),em20)
304 ELSE
305 dlam_nl(i) = zero
306 ENDIF
307c
308 ! Damage growth update
309 IF ((trdfds(i)>zero).AND.(ft(i)>zero).AND.(ft(i)<fr)) THEN
310 fg(i) = fg(i) + (one-ft(i))*dlam_nl(i)*trdfds(i)
311 ENDIF
312 fg(i) = max(fg(i),zero)
313c
314 ! Nucleation damage update
315 IF ((pla_nl(i) >= ed).AND.(ft(i)<fr)) THEN
316 ! Case for positive stress triaxiality
317 IF (triax(i)>=zero) THEN
318 fn(i) = fn(i) + an*dpla_nl(i)
319 ! Case for negative stress triaxiality
320 ELSEIF ((triax(i)<zero).AND.(triax(i)>=-third)) THEN
321 fn(i) = fn(i) + an*max(one + three*triax(i),zero)*dpla_nl(i)
322 ENDIF
323 ENDIF
324 fn(i) = max(fn(i),zero)
325c
326 ! Shear damage update
327 IF ((sigdr(i) > zero).AND.(ft(i)>zero).AND.(ft(i)<fr)) THEN
328 sigvm = sqrt(max(em20,three*(j2(i)/(normsig**2))))
329 omega = one - ((twenty7 *(j3(i)/(normsig**3)))/(two*(sigvm**3)))**2
330 omega = max(omega,zero)
331 omega = min(omega,one)
332 sdpla = (sxx(i)*normxx(i)+syy(i)*normyy(i)+ szz(i)*normzz(i)
333 . + sxy(i)*normxy(i))
334 . * dlam_nl(i)
335 fsh(i) = fsh(i) + kw*omega*ft(i)*(sdpla/sigdr(i))
336 ENDIF
337 fsh(i) = max(fsh(i),zero)
338c
339 ! Total damage update
340 ft(i) = f0 + fg(i) + fn(i) + fsh(i)
341 ft(i) = min(ft(i), fr)
342 IF (ft(i) >= fr) THEN
343 IF (off(i)==one) off(i) = four_over_5
344 ENDIF
345c
346 ! Effective update
347 IF (ft(i) < fc)THEN
348 fs(i) = ft(i)
349 ELSE
350 fs(i) = fc + (one/q1 - fc) * (ft(i)-fc)/(fr-fc)
351 ENDIF
352 fs(i) = min(fs(i),one/q1)
353c
354 ! Temperature update
355 IF (cp > zero ) THEN
356 IF (jthe == 0) THEN
357 dtemp = weitemp(i)*(one-ft(i))*yld(i)*dpla_nl(i)*eta/(rho(i)*cp)
358 temp(i) = temp(i) + dtemp
359 ENDIF
360 ENDIF
361 ENDDO
362c
363 ! Computation of the initial yield stress
364 DO i=1,nel
365 ! a) - Hardening law
366 fhard(i) = yld0 + hard*pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
367 ! b) - correction factor for strain-rate dependence (johnson-cook)
368 frate(i) = one
369 IF (epsd(i) > epsp0) frate(i) = frate(i) + jcc*log(epsd(i)/epsp0)
370 ! c) - Correction factor for thermal effects
371 ftherm(i) = one
372 IF (cp > zero) 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 signxx(i) = sigoxx(i) + (a11*depsxx(i) + a12*depsyy(i))
386 signyy(i) = sigoyy(i) + (a11*depsyy(i) + a12*depsxx(i))
387 signxy(i) = sigoxy(i) + (depsxy(i)*g)
388 signyz(i) = sigoyz(i) + (depsyz(i)*gs(i))
389 signzx(i) = sigozx(i) + (depszx(i)*gs(i))
390 ! Computation of the trace of the trial stress tensor
391 trsig(i) = signxx(i) + signyy(i)
392 sigm(i) = -trsig(i) * third
393 ! Computation of the deviatoric trial stress tensor
394 sxx(i) = signxx(i) + sigm(i)
395 syy(i) = signyy(i) + sigm(i)
396 szz(i) = sigm(i)
397 sxy(i) = signxy(i)
398 ! Second deviatoric invariant
399 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) + sxy(i)**2
400 ! Third deviatoric invariant
401 j3(i) = sxx(i)*syy(i)*szz(i) - szz(i)*sxy(i)**2
402 ! Drucker equivalent stress
403 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
404 ! Checking equivalent stress values
405 IF (fdr(i) > zero) THEN
406 sigdr(i) = kdr * exp((one/six)*log(fdr(i))) ! FDR(I)**(1/6)
407 ELSE
408 sigdr(i) = zero
409 ENDIF
410 ! Computation of the stress triaxiality and the etaT factor
411 IF (sigdr(i)>zero) THEN
412 triax(i) = (trsig(i)*third)/sigdr(i)
413 ELSE
414 triax(i) = zero
415 ENDIF
416 IF (trsig(i)<zero) THEN
417 etat(i) = zero
418 ELSE
419 etat(i) = one
420 ENDIF
421c
422 ! Update out of plane strain increment
423 IF (off(i) == one) THEN
424 dezz(i) = -nnu*depsxx(i)-nnu*depsyy(i)
425 IF (sig_dfdsig(i) > em01) THEN
426 dezz(i) = dezz(i) + nnu*(dlam_nl(i)*normxx(i))
427 . + nnu*(dlam_nl(i)*normyy(i))
428 . + dlam_nl(i)*normzz(i)
429 ENDIF
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 CUTTING PLANE METHOD (SEMI-IMPLICIT)
452 !====================================================================
453c
454 ! Number of iterations
455 niter = 3
456c
457 ! Loop over yielding elements
458#include "vectorize.inc"
459 DO ii=1,nindx
460c
461 ! Number of the element with plastic behaviour
462 i = index(ii)
463c
464 ! Initialization of the iterative Newton procedure
465 sigdr2(i) = sigdr(i)**2
466 yld2i(i) = one / yld(i)**2
467 dpxx(i) = zero
468 dpyy(i) = zero
469 dpzz(i) = zero
470 dpxy(i) = zero
471 ENDDO
472c
473 ! Loop over the iterations
474 DO iter = 1, niter
475#include "vectorize.inc"
476 DO ii=1,nindx
477 i = index(ii)
478c
479 ! Note: in this part, the purpose is to compute for each iteration
480 ! a plastic multiplier allowing to update internal variables to satisfy
481 ! the consistency condition using the cutting plane semi-implicit
482 ! iterative procedure.
483 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
484 ! -> PHI : current value of yield function (known)
485 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
486 ! into account of internal variables kinetic :
487 ! plasticity, temperature and damage (to compute)
488c
489 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
490 !-------------------------------------------------------------
491c
492 ! Derivative with respect to the equivalent stress and trace
493 yld2i(i) = one/(yld(i)**2)
494 dphi_dsig = two*sigdr(i)*yld2i(i)
495 fsinh = sinh(q2*etat(i)*trsig(i)/(yld(i)*two))
496 dphi_dtrsig = q1*q2*etat(i)*fs(i)*fsinh/yld(i)
497c
498 ! Computation of the Eulerian norm of the stress tensor
499 normsig = sqrt(signxx(i)*signxx(i)
500 . + signyy(i)*signyy(i)
501 . + two*signxy(i)*signxy(i))
502 normsig = max(normsig,one)
503c
504 ! Derivative with respect to Fdr
505 fdr(i) = (j2(i)/(normsig**2))**3 - cdr*((j3(i)/(normsig**3))**2)
506 dphi_dfdr = dphi_dsig*kdr*(one/six)*exp(-(five/six)*log(fdr(i)))
507 dsdrdj2 = dphi_dfdr*three*(j2(i)/(normsig**2))**2
508 dsdrdj3 = -dphi_dfdr*two*cdr*(j3(i)/(normsig**3))
509 ! dJ3/dSig
510 dj3dsxx = two_third*(syy(i)*szz(i))/(normsig**2)
511 . - third*(sxx(i)*szz(i))/(normsig**2)
512 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
513 dj3dsyy = - third*(syy(i)*szz(i))/(normsig**2)
514 . + two_third*(sxx(i)*szz(i))/(normsig**2)
515 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
516 dj3dszz = - third*(syy(i)*szz(i))/(normsig**2)
517 . - third*(sxx(i)*szz(i))/(normsig**2)
518 . + two_third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
519 dj3dsxy = two*(sxx(i)*sxy(i) + sxy(i)*syy(i))/(normsig**2)
520 ! dPhi/dSig
521 normxx(i) = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx + dphi_dtrsig
522 normyy(i) = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy + dphi_dtrsig
523 normzz(i) = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz + dphi_dtrsig
524 normxy(i) = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
525c
526 ! 2 - Computation of DPHI_DLAMBDA
527 !---------------------------------------------------------
528c
529 ! a) Derivative with respect stress increments tensor DSIG
530 ! --------------------------------------------------------
531 trdfds(i) = normxx(i) + normyy(i) + normzz(i)
532 dfdsig2 = normxx(i) * (a11*normxx(i) + a12*normyy(i))
533 . + normyy(i) * (a11*normyy(i) + a12*normxx(i))
534 . + normxy(i) * normxy(i) * g
535c
536 ! b) Derivatives with respect to plastic strain P
537 ! ------------------------------------------------
538c
539 ! i) Derivative of the yield stress with respect to plastic strain dYLD / dPLA
540 ! ----------------------------------------------------------------------------
541 hardp(i) = hard + qvoce*bvoce*exp(-bvoce*pla(i))
542 dyld_dpla = hardp(i)*frate(i)*ftherm(i)
543c
544 ! ii) Derivative of dPLA with respect to DLAM
545 ! -------------------------------------------
546 sig_dfdsig(i) = signxx(i) * normxx(i)
547 . + signyy(i) * normyy(i)
548 . + signxy(i) * normxy(i)
549 dpla_dlam(i) = sig_dfdsig(i) / ((one - ft(i))*yld(i))
550c
551 ! c) Derivative with respect to the yield stress
552 ! ----------------------------------------------
553 sigdr2(i) = sigdr(i)**2
554 dphi_dyld = -two*sigdr2(i)/yld(i)**3-dphi_dtrsig*trsig(i)/yld(i)
555c
556 ! d) Derivative of PHI with respect to DLAM
557 dphi_dlam(i) = - dfdsig2 + (dphi_dyld*dyld_dpla*dpla_dlam(i))
558 dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20) ,dphi_dlam(i))
559c
560 ! Computation of the plastic multiplier
561 dlam = -phi(i)/dphi_dlam(i)
562c
563 ! Plastic strains tensor update
564 dpxx(i) = dlam * normxx(i)
565 dpyy(i) = dlam * normyy(i)
566 dpzz(i) = dlam * normzz(i)
567 dpxy(i) = dlam * normxy(i)
568 trdep(i) = dpxx(i) + dpyy(i) + dpzz(i)
569c
570 ! Elasto-plastic stresses update
571 signxx(i) = signxx(i) - (a11*dpxx(i) + a12*dpyy(i))
572 signyy(i) = signyy(i) - (a11*dpyy(i) + a12*dpxx(i))
573 signxy(i) = signxy(i) - dpxy(i)*g
574 trsig(i) = signxx(i) + signyy(i)
575 sigm(i) = -trsig(i) * third
576 sxx(i) = signxx(i) + sigm(i)
577 syy(i) = signyy(i) + sigm(i)
578 szz(i) = sigm(i)
579 sxy(i) = signxy(i)
580c
581 ! Cumulated plastic strain and strain rate update
582 ddep = (dlam/((one - ft(i))*yld(i)))*sig_dfdsig(i)
583 dpla(i) = max(zero, dpla(i) + ddep)
584 pla(i) = pla(i) + ddep
585c
586 ! Drucker equivalent stress update
587 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) + sxy(i)**2
588 j3(i) = sxx(i) * syy(i) * szz(i) - szz(i) * sxy(i)**2
589 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
590 sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
591 ! Computation of the stress triaxiality and the etaT factor
592 triax(i) = (trsig(i)*third)/sigdr(i)
593 IF (trsig(i)<zero) THEN
594 etat(i) = zero
595 ELSE
596 etat(i) = one
597 ENDIF
598c
599 ! Hardening law update
600 fhard(i) = yld0 + hard * pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
601c
602 ! Yield stress update
603 yld(i) = fhard(i) * frate(i) * ftherm(i)
604 yld(i) = max(yld(i), em10)
605c
606 ! Yield function value update
607 sigdr2(i) = sigdr(i)**2
608 yld2i(i) = one / yld(i)**2
609 fcosh = cosh(q2*etat(i)*trsig(i)/(yld(i)*two))
610 fdam(i) = two*q1*fs(i)*fcosh - (q1*fs(i))**2
611 phi(i) = sigdr2(i) * yld2i(i) - one + fdam(i)
612c
613 ENDDO
614 ENDDO
615 ! End of the loop over the iterations
616 !===================================================================
617 ! - END OF PLASTIC CORRECTION WITH CUTTING PLANE ITERATIVE METHOD
618 !===================================================================
619c
620 ! Storing new values
621 DO i=1,nel
622 ! USR Outputs
623 uvar(i,1) = yld(i) ! Yield stress
624 ! Standard outputs
625 dmg(i,1) = ft(i)/fr ! Normalized total damage
626 dmg(i,2) = fg(i) ! Growth damage
627 dmg(i,3) = fn(i) ! Nucleation damage
628 dmg(i,4) = fsh(i) ! Shear damage
629 dmg(i,5) = min(ft(i),fr) ! Total damage
630 dmg(i,6) = min(fs(i),one/q1) ! Effective damage
631 seq(i) = sigdr(i) ! Equivalent stress
632 ! If element is broken
633 IF (ft(i) >= fr) THEN
634 dpla(i) = zero
635 signxx(i) = zero
636 signyy(i) = zero
637 signxy(i) = zero
638 signyz(i) = zero
639 signzx(i) = zero
640 seq(i) = zero
641 ENDIF
642 ! Plastic strain-rate (filtered)
643 dpdt = dpla(i) / max(em20,timestep)
644 epsd(i) = afiltr * dpdt + (one - afiltr) * epsd(i)
645 ! Coefficient for hourglass
646 IF (dpla(i) > zero) THEN
647 et(i) = hardp(i)*frate(i) / (hardp(i)*frate(i) + young)
648 ELSE
649 et(i) = one
650 ENDIF
651 ! Computation of the sound speed
652 soundsp(i) = sqrt((a11)/rho(i))
653 ! Storing the yield stress
654 sigy(i) = yld(i)
655 ! Computation of the thickness variation
656 thk(i) = thk(i) + dezz(i)*thkly(i)*off(i)
657 ENDDO
658c
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
for
for(i8=*sizetab-1;i8 >=0;i8--)
Definition mumps_common.c:91