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

Go to the source code of this file.

Functions/Subroutines

subroutine mat121_newton (nel, ngl, nuparam, nuvar, nfunc, ifunc, npf, tf, timestep, time, uparam, uvar, rho, pla, dpla, soundsp, epsd, off, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigy, et, seq)

Function/Subroutine Documentation

◆ mat121_newton()

subroutine mat121_newton ( integer nel,
integer, dimension(nel), intent(in) ngl,
integer nuparam,
integer nuvar,
integer nfunc,
integer, dimension(nfunc) ifunc,
integer, dimension(*) npf,
tf,
timestep,
time,
intent(in) uparam,
intent(inout) uvar,
intent(in) rho,
intent(inout) pla,
intent(inout) dpla,
intent(out) soundsp,
intent(inout) epsd,
intent(inout) off,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depszz,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
intent(in) epspxx,
intent(in) epspyy,
intent(in) epspzz,
intent(in) epspxy,
intent(in) epspyz,
intent(in) epspzx,
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) seq )

Definition at line 30 of file mat121_newton.F.

39 !=======================================================================
40 ! Implicit types
41 !=======================================================================
42#include "implicit_f.inc"
43 !=======================================================================
44 ! Common
45 !=======================================================================
46 !=======================================================================
47 ! Dummy arguments
48 !=======================================================================
49 INTEGER NEL,NUPARAM,NUVAR,NPF(*),NFUNC,IFUNC(NFUNC)
50 INTEGER ,DIMENSION(NEL), INTENT(IN) :: NGL
51 my_real
52 . time,timestep,tf(*)
53 my_real,DIMENSION(NUPARAM), INTENT(IN) ::
54 . uparam
55 my_real,DIMENSION(NEL), INTENT(IN) ::
56 . rho,
57 . depsxx,depsyy,depszz,depsxy,depsyz,depszx,
58 . epspxx,epspyy,epspzz,epspxy,epspyz,epspzx,
59 . sigoxx,sigoyy,sigozz,sigoxy,sigoyz,sigozx
60 my_real ,DIMENSION(NEL), INTENT(OUT) ::
61 . soundsp,sigy,et,
62 . signxx,signyy,signzz,signxy,signyz,signzx
63 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
64 . pla,dpla,epsd,off,seq
65 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) ::
66 . uvar
67 !=======================================================================
68 ! Local Variables
69 !=======================================================================
70 INTEGER I,II,Ivisc,ITER,NITER,NINDX,INDEX(NEL),IPOS(NEL),
71 . IAD(NEL),ILEN(NEL)
72 my_real
73 . young(nel),bulk(nel),g(nel),nu,nnu,tang(nel),lam(nel),g2(nel),
74 . afiltr,xscale_sig0,yscale_sig0,xscale_youn,yscale_youn,
75 . xscale_tang,yscale_tang
76 my_real
77 . dpdt,dlam,ddep,depxx,depyy,devepspxx,devepspyy,devepspzz,trepsp,ldav,
78 . normxx,normyy,normzz,normxy,normyz,normzx,denom,dfdsig2,dpdt_nl,depsdt,
79 . dtinv
80 my_real, DIMENSION(NEL) ::
81 . sxx,syy,szz,sxy,syz,szx,sigvm,yld,hardp,phi,dezz,yld0,dyld0depsd,
82 . dyoundepsd,dtangdepsd,trsig,dphi_dlam,test,dpxx,dpyy,dpxy,
83 . dpzz,dpyz,dpzx
84c
85 !=======================================================================
86 ! - INITIALISATION OF COMPUTATION ON TIME STEP
87 !=======================================================================
88 ! Recovering model parameters
89 ! Elastic parameters
90 young(1:nel) = uparam(1) ! Young modulus
91 bulk(1:nel) = uparam(2) ! Bulk modulus
92 g(1:nel) = uparam(3) ! Shear modulus
93 g2(1:nel) = uparam(4) ! 2*Shear modulus
94 lam(1:nel) = uparam(5) ! Lame coefficient
95 nu = uparam(6) ! Poisson ration
96 nnu = uparam(7) ! NU/(1-NU)
97 ! Flags for computation
98 ivisc = nint(uparam(12)) ! Viscosity formulation
99 ! Strain-rate filtering (if Ivisc = 0)
100 afiltr = min(one, uparam(14)*timestep)
101 ! Initial yield stress vs strain-rate curve
102 IF (ifunc(1) > 0) THEN
103 xscale_sig0 = uparam(16) ! Strain-rate scale factor
104 yscale_sig0 = uparam(17) ! Initial yield stress scale factor
105 yld0(1:nel) = zero
106 ELSE
107 yld0(1:nel) = uparam(17) ! Constant yield stress
108 dyld0depsd(1:nel) = zero
109 ENDIF
110 ! Young modulus vs strain-rate curve
111 xscale_youn = uparam(18) ! Strain-rate scale factor
112 yscale_youn = uparam(19) ! Young modulus scale factor
113 ! Tangent modulus vs strain-rate curve
114 IF (ifunc(3) > 0) THEN
115 xscale_tang = uparam(20) ! Strain-rate scale factor
116 yscale_tang = uparam(21) ! Tangent modulus scale factor
117 tang(1:nel) = zero
118 ELSE
119 tang(1:nel) = uparam(21) ! Constant tangent modulus
120 dtangdepsd(1:nel) = zero
121 ENDIF
122 dtinv = one/max(timestep,em20) ! Inverse of timestep
123c
124 ! Recovering internal variables
125 DO i=1,nel
126 IF (off(i) < em01) off(i) = zero
127 IF (off(i) < one) off(i) = off(i)*four_over_5
128 ! Standard inputs
129 dpla(i) = zero ! Initialization of the plastic strain increment
130 et(i) = one ! Initialization of hourglass coefficient
131 hardp(i) = zero ! Initialization of hourglass coefficient
132 dpxx(i) = zero ! Initialization of the XX plastic strain increment
133 dpyy(i) = zero ! Initialization of the YY plastic strain increment
134 dpzz(i) = zero ! Initialization of the ZZ plastic strain increment
135 dpxy(i) = zero ! Initialization of the XY plastic strain increment
136 dpyz(i) = zero ! Initialization of the YZ plastic strain increment
137 dpzx(i) = zero ! Initialization of the ZX plastic strain increment
138 dyld0depsd(i) = zero ! Initialization of the derivative of SIG0
139 dyoundepsd(i) = zero ! Initialization of the derivative of YOUN
140 dtangdepsd(i) = zero ! Initialization of the derivative of TANG
141 ENDDO
142c
143 ! Filling the strain-rate vector
144 IF (ivisc == 0) THEN
145 ! Compute effective strain-rate
146 DO i = 1,nel
147 trepsp = third*(epspxx(i) + epspyy(i) + epspzz(i))
148 devepspxx = epspxx(i) - trepsp
149 devepspyy = epspyy(i) - trepsp
150 devepspzz = epspzz(i) - trepsp
151 depsdt = two_third*(devepspxx**2 + devepspyy**2 + devepspzz**2 +
152 . two*(epspxy(i)**2) + two*(epspyz(i)**2) +
153 . two*(epspzx(i)**2))
154 depsdt = sqrt(max(depsdt,zero))
155 epsd(i) = afiltr * depsdt + (one - afiltr) * epsd(i)
156 ENDDO
157 ELSE
158 ! Reset plastic strain-rate
159 epsd(1:nel) = zero
160 ENDIF
161c
162 ! Compute the initial yield stress
163 IF (ifunc(1) > 0) THEN
164 ipos(1:nel) = 1
165 iad(1:nel) = npf(ifunc(1)) / 2 + 1
166 ilen(1:nel) = npf(ifunc(1)+1) / 2 - iad(1:nel) - ipos(1:nel)
167 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_sig0,dyld0depsd,yld0)
168 yld0(1:nel) = yscale_sig0*yld0(1:nel)
169 dyld0depsd(1:nel) = yscale_sig0*dyld0depsd(1:nel)
170 ENDIF
171 ! Compute the Young modulus
172 IF (ifunc(2) > 0) THEN
173 ipos(1:nel) = 1
174 iad(1:nel) = npf(ifunc(2)) / 2 + 1
175 ilen(1:nel) = npf(ifunc(2)+1) / 2 - iad(1:nel) - ipos(1:nel)
176 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_youn,dyoundepsd,young)
177 young(1:nel) = yscale_youn*young(1:nel)
178 g(1:nel) = half * young(1:nel) / (one + nu)
179 g2(1:nel) = young(1:nel) / (one + nu)
180 bulk(1:nel) = third * young(1:nel) / (one - nu*two)
181 lam(1:nel) = g2(1:nel) * nu /(one - two*nu)
182 ENDIF
183 ! Compute the Tangent modulus
184 IF (ifunc(3) > 0) THEN
185 ipos(1:nel) = 1
186 iad(1:nel) = npf(ifunc(3)) / 2 + 1
187 ilen(1:nel) = npf(ifunc(3)+1) / 2 - iad(1:nel) - ipos(1:nel)
188 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_tang,dtangdepsd,tang)
189 tang(1:nel) = yscale_tang*tang(1:nel)
190 dtangdepsd(1:nel) = yscale_tang*dtangdepsd(1:nel)
191 ENDIF
192 ! Check tangent modulus value + Assembling the yield stress
193 DO i = 1,nel
194 IF (tang(i) >= 0.99d0*young(i)) THEN
195 tang(i) = 0.99d0*young(i)
196 dtangdepsd(i) = zero
197 ENDIF
198 yld(i) = yld0(i) + (young(i)*tang(i)/(young(i)-tang(i)))*pla(i)
199 ENDDO
200c
201 !========================================================================
202 ! - COMPUTATION OF TRIAL VALUES
203 !========================================================================
204 DO i=1,nel
205 ! Computation of the trial stress tensor
206 ldav = (depsxx(i) + depsyy(i) + depszz(i)) * lam(i)
207 signxx(i) = sigoxx(i) + depsxx(i)*g2(i) + ldav
208 signyy(i) = sigoyy(i) + depsyy(i)*g2(i) + ldav
209 signzz(i) = sigozz(i) + depszz(i)*g2(i) + ldav
210 signxy(i) = sigoxy(i) + depsxy(i)*g(i)
211 signyz(i) = sigoyz(i) + depsyz(i)*g(i)
212 signzx(i) = sigozx(i) + depszx(i)*g(i)
213 ! Computation of the trace of the trial stress tensor
214 trsig(i) = signxx(i) + signyy(i) + signzz(i)
215 ! Computation of the deviatoric trial stress tensor
216 sxx(i) = signxx(i) - trsig(i) * third
217 syy(i) = signyy(i) - trsig(i) * third
218 szz(i) = signzz(i) - trsig(i) * third
219 sxy(i) = signxy(i)
220 syz(i) = signyz(i)
221 szx(i) = signzx(i)
222 ! Von Mises equivalent stress
223 sigvm(i) = three_half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + three*sxy(i)**2
224 . + three*syz(i)**2 + three*szx(i)**2
225 sigvm(i) = sqrt(sigvm(i))
226 ENDDO
227c
228 !========================================================================
229 ! - COMPUTATION OF YIELD FONCTION
230 !========================================================================
231 phi(1:nel) = sigvm(1:nel) - yld(1:nel)
232c
233 ! Checking plastic behavior for all elements
234 nindx = 0
235 DO i=1,nel
236 IF (phi(i) > zero .AND. off(i) == one) THEN
237 nindx=nindx+1
238 index(nindx)=i
239 ENDIF
240 ENDDO
241c
242 !====================================================================
243 ! - PLASTIC CORRECTION WITH CUTTING PLANE (NEWTON-ITERATION) METHOD
244 !====================================================================
245 IF (nindx > 0) THEN
246c
247 ! Number of Newton iterations
248 IF (ivisc == 0) THEN
249 niter = 3
250 ELSE
251 niter = 5
252 ENDIF
253c
254 ! Loop over the iterations
255 DO iter = 1, niter
256#include "vectorize.inc"
257 ! Loop over yielding elements
258 DO ii=1,nindx
259 i = index(ii)
260c
261 ! Note : in this part, the purpose is to compute for each iteration
262 ! a plastic multiplier allowing to update internal variables to satisfy
263 ! the consistency condition using the backward Euler implicit method
264 ! with a Newton-Raphson iterative procedure
265 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
266 ! -> PHI : current value of yield function (known)
267 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
268 ! into account of internal variables kinetic :
269 ! plasticity, temperature and damage (to compute)
270c
271 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
272 !-------------------------------------------------------------
273 normxx = three_half*sxx(i)/sigvm(i)
274 normyy = three_half*syy(i)/sigvm(i)
275 normzz = three_half*szz(i)/sigvm(i)
276 normxy = three*sxy(i)/sigvm(i)
277 normyz = three*syz(i)/sigvm(i)
278 normzx = three*szx(i)/sigvm(i)
279c
280 ! 2 - Computation of DPHI_DLAMBDA
281 !---------------------------------------------------------
282c
283 ! a) Derivative with respect stress increments tensor DSIG
284 ! --------------------------------------------------------
285 dfdsig2 = normxx * normxx * g2(i)
286 . + normyy * normyy * g2(i)
287 . + normzz * normzz * g2(i)
288 . + normxy * normxy * g(i)
289 . + normyz * normyz * g(i)
290 . + normzx * normzx * g(i)
291c
292 ! b) Derivatives with respect to plastic strain P
293 ! ------------------------------------------------
294 hardp(i) = (young(i)*tang(i)/(young(i)-tang(i)))
295 IF (ivisc == 1) THEN
296 hardp(i) = hardp(i) + dyld0depsd(i)*dtinv +
297 . ((young(i)*dtangdepsd(i)*(young(i) - tang(i))
298 . + young(i)*tang(i)*dtangdepsd(i))/
299 . ((young(i) - tang(i))**2))*dtinv*pla(i)
300 ENDIF
301c
302 ! 3 - Computation of plastic multiplier and variables update
303 !----------------------------------------------------------
304c
305 ! Derivative of PHI with respect to DLAM
306 dphi_dlam(i) = - dfdsig2 - hardp(i)
307 dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20),dphi_dlam(i))
308c
309 ! Computation of the plastic multiplier
310 dlam = -phi(i)/dphi_dlam(i)
311c
312 ! Plastic strains tensor update
313 dpxx(i) = dlam * normxx
314 dpyy(i) = dlam * normyy
315 dpzz(i) = dlam * normzz
316 dpxy(i) = dlam * normxy
317 dpyz(i) = dlam * normyz
318 dpzx(i) = dlam * normzx
319c
320 ! Elasto-plastic stresses update
321 signxx(i) = signxx(i) - dpxx(i)*g2(i)
322 signyy(i) = signyy(i) - dpyy(i)*g2(i)
323 signzz(i) = signzz(i) - dpzz(i)*g2(i)
324 signxy(i) = signxy(i) - dpxy(i)*g(i)
325 signyz(i) = signyz(i) - dpyz(i)*g(i)
326 signzx(i) = signzx(i) - dpzx(i)*g(i)
327c
328 ! Computation of the trace of the new stress tensor
329 trsig(i) = signxx(i) + signyy(i) + signzz(i)
330 ! Computation of the new deviatoric stress tensor
331 sxx(i) = signxx(i) - trsig(i) * third
332 syy(i) = signyy(i) - trsig(i) * third
333 szz(i) = signzz(i) - trsig(i) * third
334 sxy(i) = signxy(i)
335 syz(i) = signyz(i)
336 szx(i) = signzx(i)
337c
338 ! Cumulated plastic strain and strain rate update
339 dpla(i) = max(zero,dpla(i) + dlam)
340 pla(i) = max(zero,pla(i) + dlam)
341 IF (ivisc == 1) THEN
342 epsd(i) = dpla(i)*dtinv
343 ENDIF
344c
345 ! Von Mises equivalent stress update
346 sigvm(i) = three_half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + three*sxy(i)**2
347 . + three*syz(i)**2 + three*szx(i)**2
348 sigvm(i) = sqrt(sigvm(i))
349c
350 IF (ivisc == 0) THEN
351 ! Yield stress update
352 yld(i) = yld0(i) + (young(i)*tang(i)/(young(i)-tang(i)))*pla(i)
353 ! Yield function value update
354 phi(i) = sigvm(i) - yld(i)
355 ENDIF
356c
357 ENDDO
358 ! End of the loop over yielding elements
359c
360 ! Update variable for full viscoplastic formulation
361 IF (ivisc == 1) THEN
362 ! Compute the initial yield stress
363 ipos(1:nel) = 1
364 iad(1:nel) = npf(ifunc(1)) / 2 + 1
365 ilen(1:nel) = npf(ifunc(1)+1) / 2 - iad(1:nel) - ipos(1:nel)
366 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_sig0,dyld0depsd,yld0)
367 yld0(1:nel) = yscale_sig0*yld0(1:nel)
368 dyld0depsd(1:nel) = yscale_sig0*dyld0depsd(1:nel)
369 ! Compute the Tangent modulus
370 IF (ifunc(3) > 0) THEN
371 ipos(1:nel) = 1
372 iad (1:nel) = npf(ifunc(3)) / 2 + 1
373 ilen(1:nel) = npf(ifunc(3)+1) / 2 - iad(1:nel) - ipos(1:nel)
374 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_tang,dtangdepsd,tang)
375 tang(1:nel) = yscale_tang*tang(1:nel)
376 dtangdepsd(1:nel) = yscale_tang*dtangdepsd(1:nel)
377 ENDIF
378 ! Updating values
379 DO ii=1,nindx
380 i = index(ii)
381 ! Check tangent modulus value
382 IF (tang(i) >= 0.99d0*young(i)) THEN
383 tang(i) = 0.99d0*young(i)
384 dtangdepsd(i) = zero
385 ENDIF
386 ! Yield stress update
387 yld(i) = yld0(i) + (young(i)*tang(i)/(young(i)-tang(i)))*pla(i)
388 ! Yield function value update
389 phi(i) = sigvm(i) - yld(i)
390 ENDDO
391 ENDIF
392 ENDDO
393 ENDIF
394 ! End of the loop over the iterations
395 !=========================================================================
396 ! - END OF PLASTIC CORRECTION WITH CUTTING PLANE (NEWTON-ITERATION) METHOD
397 !=========================================================================
398c
399 ! Storing new values
400 DO i=1,nel
401 ! USR Outputs
402 seq(i) = sigvm(i) ! SIGEQ
403 ! Coefficient for hourglass
404 IF (dpla(i) > zero) THEN
405 et(i) = hardp(i) / (hardp(i) + young(i))
406 ELSE
407 et(i) = one
408 ENDIF
409 ! Computation of the sound speed
410 soundsp(i) = sqrt((bulk(i) + four_over_3*g(i))/rho(i))
411 ! Storing the yield stress
412 sigy(i) = yld(i)
413 ENDDO
414c
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
vinter2
subroutine vinter2(tf, iad, ipos, ilen, nel0, x, dydx, y)
Definition vinter.F:143