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

Go to the source code of this file.

Functions/Subroutines

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

Function/Subroutine Documentation

◆ mat104_nodam_newton()

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

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