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