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