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

Go to the source code of this file.

Functions/Subroutines

subroutine mat104c_nodam_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, dplanl, temp, seq, inloc)

Function/Subroutine Documentation

◆ mat104c_nodam_newton()

subroutine mat104c_nodam_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) dplanl,
intent(inout) temp,
intent(inout) seq,
integer inloc )

Definition at line 28 of file mat104c_nodam_newton.F.

37 !=======================================================================
38 ! Implicit types
39 !=======================================================================
40#include "implicit_f.inc"
41 !=======================================================================
42 ! Common
43 !=======================================================================
44 !=======================================================================
45 ! Dummy arguments
46 !=======================================================================
47 INTEGER NEL,NUPARAM,NUVAR,JTHE,INLOC
48 INTEGER ,DIMENSION(NEL), INTENT(IN) :: NGL
49 my_real
50 . time,timestep
51 my_real,DIMENSION(NUPARAM), INTENT(IN) ::
52 . uparam
53 my_real,DIMENSION(NEL), INTENT(IN) ::
54 . rho,tempel,dplanl,
55 . depsxx,depsyy,depsxy,depsyz,depszx,
56 . sigoxx,sigoyy,sigoxy,sigoyz,sigozx,
57 . gs , thkly
58c
59 my_real ,DIMENSION(NEL), INTENT(OUT) ::
60 . soundsp,sigy,et,
61 . signxx,signyy,signxy,signyz,signzx
62c
63 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
64 . pla,dpla,epsd,off,thk,temp,seq
65 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) ::
66 . uvar
67 !=======================================================================
68 ! Local Variables
69 !=======================================================================
70 INTEGER I,II,ITER,NITER,NINDX,INDEX(NEL)
71c
72 my_real
73 . young,bulk,lam,g,g2,nu,cdr,kdr,hard,yld0,qvoce,bvoce,jcc,
74 . epsp0,mtemp,tini,tref,eta,cp,dpis,dpad,asrate,afiltr,a11,a12,
75 . nnu, normsig
76 my_real
77 .
78 . dtemp,dpdt,dlam,ddep
79 my_real
80 . dsdrdj2,dsdrdj3,
81 . dj3dsxx,dj3dsyy,dj3dsxy,dj3dszz,
82 .
83 .
84 . normxx,normyy,normxy,normzz,
85 . sig_dfdsig,dfdsig2,
86 . dphi_dsig,dphi_dyld,dphi_dfdr,dpdt_nl,
87 . dyld_dpla,dyld_dtemp,dtemp_dlam,dlam_nl
88c
89 my_real, DIMENSION(NEL) ::
90 . trsig,
91 . sxx,syy,sxy,szz,sigm,j2,j3,sigdr,yld,weitemp,
92 . hardp,fhard,frate,ftherm,fdr,phi,dpla_dlam,
93 . dezz,dphi_dlam,dpxx,dpyy,dpxy,dpzz,sigdr2,yld2i
94 !=======================================================================
95 ! DRUCKER - VOCE - JOHNSON-COOK MATERIAL LAW
96 !=======================================================================
97 !UVAR(1) EPSD PLASTIC STRAIN RATE SAVED FOR FILTERING
98c
99 !=======================================================================
100 ! - INITIALISATION OF COMPUTATION ON TIME STEP
101 !=======================================================================
102 ! Recovering model parameters
103 ! Elastic parameters
104 young = uparam(1) ! Young modulus
105 bulk = uparam(2) ! Bulk modulus
106 g = uparam(3) ! Shear modulus
107 g2 = uparam(4) ! 2*Shear modulus
108 lam = uparam(5) ! Lambda Hooke parameter
109 nu = uparam(6) ! Poisson ration
110 nnu = uparam(7) ! NU/(1-NU)
111 a11 = uparam(9) ! Diagonal term, elastic matrix in plane stress
112 a12 = uparam(10) ! Non-diagonal term, elastic matrix in plane stress
113 ! Plastic criterion and hardening parameters [Drucker, 1948]
114 cdr = uparam(12) ! Drucker coefficient
115 kdr = uparam(13) ! Drucker 1/K coefficient
116 tini = uparam(14) ! Initial temperature
117 hard = uparam(15) ! Linear hardening
118 yld0 = uparam(16) ! Initial yield stress
119 qvoce = uparam(17) ! 1st Voce parameter
120 bvoce = uparam(18) ! 2nd Voce parameter
121 ! Strain-rate dependence parameters
122 jcc = uparam(20) ! Johnson-Cook strain rate coefficient
123 epsp0 = uparam(21) ! Johnson-Cook reference strain rate
124 ! Thermal softening and self-heating parameters
125 mtemp = uparam(22) ! Thermal softening slope
126 tref = uparam(23) ! Reference temperature
127 eta = uparam(24) ! Taylor-Quinney coefficient
128 cp = uparam(25) ! Thermal mass capacity
129 dpis = uparam(26) ! Isothermal plastic strain rate
130 dpad = uparam(27) ! Adiabatic plastic strain rate
131 ! Plastic strain-rate filtering parameters
132 asrate = uparam(28) ! Plastic strain rate filtering frequency
133 afiltr = asrate*timestep/(asrate*timestep + one)
134c
135 ! Recovering internal variables
136 DO i=1,nel
137 IF (off(i) < em01) off(i) = zero
138 IF (off(i) < one) off(i) = off(i)*four_over_5
139 ! Standard inputs
140 dpla(i) = zero ! Initialization of the plastic strain increment
141 et(i) = one ! Initialization of hourglass coefficient
142 hardp(i) = zero ! Initialization of hardening modulus
143 ENDDO
144!
145 epsd(1:nel) = uvar(1:nel,1)
146!
147 ! Initialization of temperature and self-heating weight factor
148 IF (time == zero) THEN
149 temp(1:nel) = tini
150 ENDIF
151 IF (cp > zero) THEN
152 IF (jthe == 0) THEN
153 IF (inloc == 0) THEN
154 DO i=1,nel
155 IF (epsd(i) < dpis) THEN
156 weitemp(i) = zero
157 ELSEIF (epsd(i) > dpad) THEN
158 weitemp(i) = one
159 ELSE
160 weitemp(i) = ((epsd(i)-dpis)**2 )
161 . * (three*dpad - two*epsd(i) - dpis)
162 . / ((dpad-dpis)**3)
163 ENDIF
164 ENDDO
165 ENDIF
166 ELSE
167 temp(1:nel) = tempel(1:nel)
168 ENDIF
169 ENDIF
170c
171 ! Computation of the initial yield stress
172 DO i = 1,nel
173 ! a) - Hardening law
174 fhard(i) = yld0 + hard*pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
175 ! b) - Correction factor for strain-rate dependence (Johnson-Cook)
176 frate(i) = one
177 IF (epsd(i) > epsp0) frate(i) = frate(i) + jcc*log(epsd(i)/epsp0)
178 ! c) - Correction factor for thermal softening
179 ftherm(i) = one
180 IF (cp > zero) ftherm(i) = one - mtemp * (temp(i) - tref)
181 ! d) - Computation of the yield function and value check
182 yld(i) = fhard(i)*frate(i)*ftherm(i)
183 ! e) - Checking values
184 yld(i) = max(em10, yld(i))
185 ENDDO
186c
187 !========================================================================
188 ! - COMPUTATION OF TRIAL VALUES
189 !========================================================================
190 DO i=1,nel
191c
192 ! Computation of the trial stress tensor
193 signxx(i) = sigoxx(i) + a11*depsxx(i) + a12*depsyy(i)
194 signyy(i) = sigoyy(i) + a11*depsyy(i) + a12*depsxx(i)
195 signxy(i) = sigoxy(i) + depsxy(i)*g
196 signyz(i) = sigoyz(i) + depsyz(i)*gs(i)
197 signzx(i) = sigozx(i) + depszx(i)*gs(i)
198 ! Computation of the trace of the trial stress tensor
199 trsig(i) = signxx(i) + signyy(i)
200 sigm(i) = -trsig(i) * third
201 ! Computation of the deviatoric trial stress tensor
202 sxx(i) = signxx(i) + sigm(i)
203 syy(i) = signyy(i) + sigm(i)
204 szz(i) = sigm(i)
205 sxy(i) = signxy(i)
206 dezz(i) = -nnu*depsxx(i) - nnu*depsyy(i)
207 ! Second deviatoric invariant
208 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) + sxy(i)**2
209 ! Third deviatoric invariant
210 j3(i) = sxx(i)*syy(i)*szz(i) - szz(i)*sxy(i)**2
211 ! Drucker equivalent stress
212 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
213 ! Checking equivalent stress values
214 IF (fdr(i) > zero) THEN
215 sigdr(i) = kdr * exp((one/six)*log(fdr(i))) ! FDR(I)**(1/6)
216 ELSE
217 sigdr(i) = zero
218 ENDIF
219 ENDDO
220c
221 !========================================================================
222 ! - COMPUTATION OF YIELD FONCTION
223 !========================================================================
224 phi(1:nel) = (sigdr(1:nel) / yld(1:nel))**2 - one
225c
226 ! Checking plastic behavior for all elements
227 nindx = 0
228 DO i=1,nel
229 IF (phi(i) >= zero .AND. off(i) == one) THEN
230 nindx=nindx+1
231 index(nindx)=i
232 ENDIF
233 ENDDO
234c
235 !====================================================================
236 ! - PLASTIC CORRECTION WITH CUTTING PLANE METHOD (SEMI-IMPLICIT)
237 !====================================================================
238c
239 ! Number of iterations
240 niter = 3
241c
242 ! Loop over yielding elements
243#include "vectorize.inc"
244 DO ii=1,nindx
245c
246 ! Number of the element with plastic behaviour
247 i = index(ii)
248c
249 ! Initialization of the iterative Newton procedure
250 sigdr2(i) = sigdr(i)**2
251 yld2i(i) = one / yld(i)**2
252 dpxx(i) = zero
253 dpyy(i) = zero
254 dpzz(i) = zero
255 dpxy(i) = zero
256 ENDDO
257c
258 ! Loop over the iterations
259 DO iter = 1, niter
260#include "vectorize.inc"
261 ! Loop over yielding elements
262 DO ii=1,nindx
263 i = index(ii)
264c
265 ! Note: in this part, the purpose is to compute for each iteration
266 ! a plastic multiplier allowing to update internal variables to satisfy
267 ! the consistency condition using the cutting plane semi-implicit
268 ! iterative procedure.
269 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
270 ! -> PHI : current value of yield function (known)
271 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
272 ! into account of internal variables kinetic :
273 ! plasticity, temperature and damage (to compute)
274c
275 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
276 !-------------------------------------------------------------
277c
278 ! Derivative with respect to the equivalent stress
279 yld2i(i) = one/(yld(i)**2)
280 dphi_dsig = two*sigdr(i)*yld2i(i)
281c
282 ! Computation of the Eulerian norm of the stress tensor
283 normsig = sqrt(signxx(i)*signxx(i)
284 . + signyy(i)*signyy(i)
285 . + two*signxy(i)*signxy(i))
286 normsig = max(normsig,one)
287c
288 ! Derivative with respect to Fdr
289 fdr(i) = (j2(i)/(normsig**2))**3 - cdr*((j3(i)/(normsig**3))**2)
290 dphi_dfdr = dphi_dsig*kdr*(one/six)*exp(-(five/six)*log(fdr(i)))
291 dsdrdj2 = dphi_dfdr*three*(j2(i)/(normsig**2))**2
292 dsdrdj3 = -dphi_dfdr*two*cdr*(j3(i)/(normsig**3))
293 ! dJ3/dSig
294 dj3dsxx = two_third*(syy(i)*szz(i))/(normsig**2)
295 . - third*(sxx(i)*szz(i))/(normsig**2)
296 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
297 dj3dsyy = - third*(syy(i)*szz(i))/(normsig**2)
298 . + two_third*(sxx(i)*szz(i))/(normsig**2)
299 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
300 dj3dszz = - third*(syy(i)*szz(i))/(normsig**2)
301 . - third*(sxx(i)*szz(i))/(normsig**2)
302 . + two_third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
303 dj3dsxy = two*(sxx(i)*sxy(i) + sxy(i)*syy(i))/(normsig**2)
304 ! dPhi/dSig
305 normxx = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx
306 normyy = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy
307 normzz = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz
308 normxy = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
309c
310 ! 2 - Computation of DPHI_DLAMBDA
311 !---------------------------------------------------------
312c
313 ! a) Derivative with respect stress increments tensor DSIG
314 ! --------------------------------------------------------
315 dfdsig2 = normxx * (a11*normxx + a12*normyy)
316 . + normyy * (a11*normyy + a12*normxx)
317 . + normxy * normxy * g
318c
319 ! b) Derivatives with respect to plastic strain P
320 ! ------------------------------------------------
321c
322 ! i) Derivative of the yield stress with respect to plastic strain dYLD / dPLA
323 ! ----------------------------------------------------------------------------
324 hardp(i) = hard + qvoce*bvoce*exp(-bvoce*pla(i))
325 dyld_dpla = hardp(i)*frate(i)*ftherm(i)
326c
327 ! ii) Derivative of dPLA with respect to DLAM
328 ! -------------------------------------------
329 sig_dfdsig = signxx(i) * normxx
330 . + signyy(i) * normyy
331 . + signxy(i) * normxy
332 dpla_dlam(i) = sig_dfdsig / yld(i)
333c
334 ! c) Derivatives with respect to the temperature TEMP
335 ! ---------------------------------------------------
336 IF (jthe == 0 .AND. cp > zero .AND. inloc == 0) THEN
337 ! i) Derivative of the yield stress with respect to temperature dYLD / dTEMP
338 ! ---------------------------------------------------------------------------
339 dyld_dtemp = -fhard(i)*frate(i)*mtemp
340 ! ii) Derivative of the temperature TEMP with respect to DLAM
341 ! -----------------------------------------------------------
342 dtemp_dlam = weitemp(i)*eta/(rho(i)*cp)*sig_dfdsig
343 ELSE
344 dyld_dtemp = zero
345 dtemp_dlam = zero
346 ENDIF
347c
348 ! d) Derivative with respect to the yield stress
349 ! ----------------------------------------------
350 dphi_dyld = -two*sigdr2(i)/(yld(i)**3)
351c
352 ! 3 - Computation of plastic multiplier and variables update
353 !----------------------------------------------------------
354c
355 ! Derivative of PHI with respect to DLAM
356 dphi_dlam(i) = - dfdsig2 + (dphi_dyld*dyld_dpla*dpla_dlam(i))
357 IF (jthe == 0 .AND. cp > zero .AND. inloc == 0) THEN
358 dphi_dlam(i) = dphi_dlam(i) + dphi_dyld*dyld_dtemp*dtemp_dlam
359 ENDIF
360 dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20) ,dphi_dlam(i))
361c
362 ! Computation of the plastic multiplier
363 dlam = -phi(i)/dphi_dlam(i)
364c
365 ! Plastic strains tensor update
366 dpxx(i) = dlam * normxx
367 dpyy(i) = dlam * normyy
368 dpzz(i) = dlam * normzz
369 dpxy(i) = dlam * normxy
370c
371 ! Elasto-plastic stresses update
372 signxx(i) = signxx(i) - (a11*dpxx(i) + a12*dpyy(i))
373 signyy(i) = signyy(i) - (a11*dpyy(i) + a12*dpxx(i))
374 signxy(i) = signxy(i) - dpxy(i)*g
375 trsig(i) = signxx(i) + signyy(i)
376 sigm(i) = -trsig(i) * third
377 sxx(i) = signxx(i) + sigm(i)
378 syy(i) = signyy(i) + sigm(i)
379 szz(i) = sigm(i)
380 sxy(i) = signxy(i)
381c
382 ! Cumulated plastic strain and strain rate update
383 ddep = (dlam/yld(i))*sig_dfdsig
384 dpla(i) = max(zero, dpla(i) + ddep)
385 pla(i) = pla(i) + ddep
386c
387 ! Drucker equivalent stress update
388 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) + sxy(i)**2
389 j3(i) = sxx(i) * syy(i) * szz(i) - szz(i) * sxy(i)**2
390 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
391 sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
392c
393 ! Temperature update
394 IF (jthe == 0 .AND. cp > zero .AND. inloc == 0) THEN
395 dtemp = weitemp(i)*yld(i)*ddep*eta/(rho(i)*cp)
396 temp(i) = temp(i) + dtemp
397 ftherm(i) = one - mtemp*(temp(i) - tref)
398 ENDIF
399c
400 ! Hardening law update
401 fhard(i) = yld0 + hard*pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
402c
403 ! Yield stress update
404 yld(i) = fhard(i)*frate(i)*ftherm(i)
405 yld(i) = max(yld(i), em10)
406c
407 ! Yield function value update
408 sigdr2(i) = sigdr(i)**2
409 yld2i(i) = one / yld(i)**2
410 phi(i) = sigdr2(i) * yld2i(i) - one
411c
412 ! Transverse strain update
413 IF (inloc == 0) THEN
414 dezz(i) = dezz(i) + nnu*dpxx(i) + nnu*dpyy(i) + dpzz(i)
415 ENDIF
416c
417 ENDDO
418 ! End of the loop over yielding elements
419 ENDDO
420 ! End of the loop over the iterations
421 !===================================================================
422 ! - END OF PLASTIC CORRECTION WITH CUTTING PLANE ITERATIVE METHOD
423 !===================================================================
424c
425 ! Storing new values
426 DO i=1,nel
427 ! USR Outputs
428 seq(i) = sigdr(i)
429 ! Coefficient for hourglass
430 IF (dpla(i) > zero) THEN
431 et(i) = hardp(i)*frate(i) / (hardp(i)*frate(i) + young)
432 ELSE
433 et(i) = one
434 ENDIF
435 ! Plastic strain-rate (filtered)
436 dpdt = dpla(i) / max(em20,timestep)
437 epsd(i) = afiltr * dpdt + (one - afiltr) * epsd(i)
438 uvar(i,1) = epsd(i)
439 ! Computation of the sound speed
440 soundsp(i) = sqrt(a11/rho(i))
441 ! Storing the yield stress
442 sigy(i) = yld(i)
443 ! Non-local thickness variation + temperature
444 IF ((inloc > 0).AND.(off(i) == one).AND.(dplanl(i)>zero)) THEN
445 yld2i(i) = one/(yld(i)**2)
446 dphi_dsig = two*sigdr(i)*yld2i(i)
447 normsig = sqrt(signxx(i)*signxx(i)
448 . + signyy(i)*signyy(i)
449 . + two*signxy(i)*signxy(i))
450 normsig = max(normsig,one)
451 fdr(i) = (j2(i)/(normsig**2))**3 - cdr*((j3(i)/(normsig**3))**2)
452 IF (fdr(i) > zero) THEN
453 dphi_dfdr = dphi_dsig*kdr*(one/six)*exp(-(five/six)*log(fdr(i)))
454 ELSE
455 dphi_dfdr = zero
456 ENDIF
457 dsdrdj2 = dphi_dfdr*three*(j2(i)/(normsig**2))**2
458 dsdrdj3 = -dphi_dfdr*two*cdr*(j3(i)/(normsig**3))
459 dj3dsxx = two_third*(syy(i)*szz(i))/(normsig**2)
460 . - third*(sxx(i)*szz(i))/(normsig**2)
461 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
462 dj3dsyy = - third*(syy(i)*szz(i))/(normsig**2)
463 . + two_third*(sxx(i)*szz(i))/(normsig**2)
464 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
465 dj3dszz = - third*(syy(i)*szz(i))/(normsig**2)
466 . - third*(sxx(i)*szz(i))/(normsig**2)
467 . + two_third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
468 dj3dsxy = two*(sxx(i)*sxy(i) + sxy(i)*syy(i))/(normsig**2)
469 normxx = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx
470 normyy = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy
471 normzz = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz
472 normxy = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
473 sig_dfdsig = signxx(i) * normxx
474 . + signyy(i) * normyy
475 . + signxy(i) * normxy
476 IF (sig_dfdsig > em01) THEN
477 dlam_nl = (yld(i)*dplanl(i))/sig_dfdsig
478 dezz(i) = dezz(i) + nnu*(dlam_nl*normxx)
479 . + nnu*(dlam_nl*normyy)
480 . + dlam_nl*normzz
481 ENDIF
482 IF (cp > zero .AND. jthe == 0) THEN
483 ! Computation of the weight factor
484 dpdt_nl = dplanl(i)/max(timestep,em20)
485 IF (dpdt_nl < dpis) THEN
486 weitemp(i) = zero
487 ELSEIF (dpdt_nl > dpad) THEN
488 weitemp(i) = one
489 ELSE
490 weitemp(i) = ((dpdt_nl-dpis)**2 )
491 . * (three*dpad - two*dpdt_nl - dpis)
492 . / ((dpad-dpis)**3)
493 ENDIF
494 dtemp = weitemp(i)*dplanl(i)*yld(i)*eta/(rho(i)*cp)
495 temp(i) = temp(i) + dtemp
496 ENDIF
497 ENDIF
498 ! Computation of the thickness variation
499 thk(i) = thk(i) + dezz(i)*thkly(i)*off(i)
500 ENDDO
501c
my_real
#define my_real
Definition cppsort.cpp:32
max
#define max(a, b)
Definition macros.h:21