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