OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat121_newton.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!|| mat121_newton ../engine/source/materials/mat/mat121/mat121_newton.F
25!||--- called by ------------------------------------------------------
26!|| sigeps121 ../engine/source/materials/mat/mat121/sigeps121.F
27!||--- calls -----------------------------------------------------
28!|| vinter2 ../engine/source/tools/curve/vinter.F
29!||====================================================================
30 SUBROUTINE mat121_newton(
31 1 NEL ,NGL ,NUPARAM ,NUVAR ,NFUNC ,IFUNC ,NPF ,
32 2 TF ,TIMESTEP,TIME ,UPARAM ,UVAR ,RHO ,PLA ,
33 3 DPLA ,SOUNDSP ,EPSD ,OFF ,
34 4 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
35 5 EPSPXX ,EPSPYY ,EPSPZZ ,EPSPXY ,EPSPYZ ,EPSPZX ,
36 6 SIGOXX ,SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
37 7 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
38 8 SIGY ,ET ,SEQ )
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
415 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
mat121_newton
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)
Definition mat121_newton.F:39
vinter2
subroutine vinter2(tf, iad, ipos, ilen, nel0, x, dydx, y)
Definition vinter.F:143