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 . dpdt,dlam,ddep,depxx,depyy,devepspxx,devepspyy,devepspzz,trepsp,
78 . normxx,normyy,normxy,denom,dfdsig2,dpdt_nl,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,test,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
138c
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) * epsd(i)
151 ENDDO
152 ELSE
153 ! Reset plastic strain-rate
154 epsd(1:nel) = zero
155 ENDIF
156c
157 ! Compute the initial yield stress
158 IF (ifunc(1) > 0) THEN
159 ipos(1:nel) = 1
160 iad(1:nel) = npf(ifunc(1)) / 2 + 1
161 ilen(1:nel) = npf(ifunc(1)+1) / 2 - iad(1:nel) - ipos(1:nel)
162 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_sig0,dyld0depsd,yld0)
163 yld0(1:nel) = yscale_sig0*yld0(1:nel)
164 dyld0depsd(1:nel) = yscale_sig0*dyld0depsd(1:nel)
165 ENDIF
166 ! Compute the Young modulus
167 IF (ifunc(2) > 0) THEN
168 ipos(1:nel) = 1
169 iad(1:nel) = npf(ifunc(2)) / 2 + 1
170 ilen(1:nel) = npf(ifunc(2)+1) / 2 - iad(1:nel) - ipos(1:nel)
171 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_youn,dyoundepsd,young)
172 young(1:nel) = yscale_youn*young(1:nel)
173 g(1:nel) = half * young(1:nel) / (one + nu)
174 bulk(1:nel) = third * young(1:nel) / (one - nu*two)
175 a11(1:nel) = young(1:nel) / (one - nu*nu)
176 a12(1:nel) = a11(1:nel) * nu
177 ENDIF
178 ! Compute the Tangent modulus
179 IF (ifunc(3) > 0) THEN
180 ipos(1:nel) = 1
181 iad(1:nel) = npf(ifunc(3)) / 2 + 1
182 ilen(1:nel) = npf(ifunc(3)+1) / 2 - iad(1:nel) - ipos(1:nel)
183 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_tang,dtangdepsd,tang)
184 tang(1:nel) = yscale_tang*tang(1:nel)
185 dtangdepsd(1:nel) = yscale_tang*dtangdepsd(1:nel)
186 ENDIF
187 ! Check tangent modulus value + Assembling the yield stress
188 DO i = 1,nel
189 IF (tang(i) >= 0.99d0*young(i)) THEN
190 tang(i) = 0.99d0*young(i)
191 dtangdepsd(i) = zero
192 ENDIF
193 yld(i) = yld0(i) + (young(i)*tang(i)/(young(i)-tang(i)))*pla(i)
194 ENDDO
195c
196 !========================================================================
197 ! - COMPUTATION OF TRIAL VALUES
198 !========================================================================
199 DO i=1,nel
200 ! Computation of the trial stress tensor
201 signxx(i) = sigoxx(i) + a11(i)*depsxx(i) + a12(i)*depsyy(i)
202 signyy(i) = sigoyy(i) + a11(i)*depsyy(i) + a12(i)*depsxx(i)
203 signxy(i) = sigoxy(i) + depsxy(i)*g(i)
204 signyz(i) = sigoyz(i) + depsyz(i)*gs(i)
205 signzx(i) = sigozx(i) + depszx(i)*gs(i)
206 ! Computation of the trace of the trial stress tensor
207 trsig(i) = signxx(i) + signyy(i)
208 ! Computation of the deviatoric trial stress tensor
209 sxx(i) = signxx(i) - trsig(i) * third
210 syy(i) = signyy(i) - trsig(i) * third
211 szz(i) = -trsig(i) * third
212 sxy(i) = signxy(i)
213 ! Von Mises equivalent stress
214 sigvm(i) = three_half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + three*sxy(i)**2
215 sigvm(i) = sqrt(sigvm(i))
216 ENDDO
217c
218 !========================================================================
219 ! - COMPUTATION OF YIELD FONCTION
220 !========================================================================
221 phi(1:nel) = sigvm(1:nel) - yld(1:nel)
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 (NEWTON-ITERATION) METHOD
234 !====================================================================
235 IF (nindx > 0) THEN
236c
237 ! Number of Newton iterations
238 IF (ivisc == 0) THEN
239 niter = 3
240 ELSE
241 niter = 5
242 ENDIF
243c
244 ! Loop over the iterations
245 DO iter = 1, niter
246#include "vectorize.inc"
247 ! Loop over yielding elements
248 DO ii=1,nindx
249c
250 ! Number of the element with plastic behaviour
251 i = index(ii)
252c
253 ! Note : in this part, the purpose is to compute for each iteration
254 ! a plastic multiplier allowing to update internal variables to satisfy
255 ! the consistency condition using the backward Euler implicit method
256 ! with a cutting plane iterative procedure
257 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
258 ! -> PHI : current value of yield function (known)
259 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
260 ! into account of internal variables kinetic :
261 ! plasticity, temperature and damage (to compute)
262c
263 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
264 !-------------------------------------------------------------
265 normxx = three_half*sxx(i)/sigvm(i)
266 normyy = three_half*syy(i)/sigvm(i)
267 normxy = three*sxy(i)/sigvm(i)
268c
269 ! 2 - Computation of DPHI_DLAMBDA
270 !---------------------------------------------------------
271c
272 ! a) Derivative with respect stress increments tensor DSIG
273 ! --------------------------------------------------------
274 dfdsig2 = normxx * (a11(i)*normxx + a12(i)*normyy)
275 . + normyy * (a11(i)*normyy + a12(i)*normxx)
276 . + normxy * normxy * g(i)
277c
278 ! b) Derivatives with respect to plastic strain P
279 ! ------------------------------------------------
280 hardp(i) = (young(i)*tang(i)/(young(i)-tang(i)))
281 IF (ivisc == 1) THEN
282 hardp(i) = hardp(i) + dyld0depsd(i)*dtinv +
283 . ((young(i)*dtangdepsd(i)*(young(i) - tang(i))
284 . + young(i)*tang(i)*dtangdepsd(i))/
285 . ((young(i) - tang(i))**2))*dtinv*pla(i)
286 ENDIF
287c
288 ! 3 - Computation of plastic multiplier and variables update
289 !----------------------------------------------------------
290c
291 ! Derivative of PHI with respect to DLAM
292 dphi_dlam(i) = - dfdsig2 - hardp(i)
293 dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20),dphi_dlam(i))
294c
295 ! Computation of the plastic multiplier
296 dlam = -phi(i)/dphi_dlam(i)
297c
298 ! Plastic strains tensor update
299 dpxx(i) = dlam * normxx
300 dpyy(i) = dlam * normyy
301 dpxy(i) = dlam * normxy
302c
303 ! Elasto-plastic stresses update
304 signxx(i) = signxx(i) - (a11(i)*dpxx(i) + a12(i)*dpyy(i))
305 signyy(i) = signyy(i) - (a11(i)*dpyy(i) + a12(i)*dpxx(i))
306 signxy(i) = signxy(i) - dpxy(i)*g(i)
307 trsig(i) = signxx(i) + signyy(i)
308 sxx(i) = signxx(i) - trsig(i) * third
309 syy(i) = signyy(i) - trsig(i) * third
310 szz(i) = - trsig(i) * third
311 sxy(i) = signxy(i)
312c
313 ! Cumulated plastic strain and strain rate update
314 dpla(i) = max(zero,dpla(i) + dlam)
315 pla(i) = max(zero,pla(i) + dlam)
316 IF (ivisc == 1) THEN
317 epsd(i) = dpla(i)*dtinv
318 ENDIF
319c
320 ! Von Mises equivalent stress update
321 sigvm(i) = three_half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + three*sxy(i)**2
322 sigvm(i) = sqrt(sigvm(i))
323c
324 IF (ivisc == 0) THEN
325 ! Yield stress update
326 yld(i) = yld0(i) + (young(i)*tang(i)/(young(i)-tang(i)))*pla(i)
327 ! Yield function value update
328 phi(i) = sigvm(i) - yld(i)
329 ENDIF
330c
331 ! Transverse strain update
332 IF (inloc == 0) THEN
333 dezz(i) = dezz(i) - dpxx(i) - dpyy(i)
334 ENDIF
335c
336 ENDDO
337 ! End of the loop over yielding elements
338c
339 ! Update variable for full viscoplastic formulation
340 IF (ivisc == 1) THEN
341 ! Compute the initial yield stress
342 ipos(1:nel) = 1
343 iad(1:nel) = npf(ifunc(1)) / 2 + 1
344 ilen(1:nel) = npf(ifunc(1)+1) / 2 - iad(1:nel) - ipos(1:nel)
345 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_sig0,dyld0depsd,yld0)
346 yld0(1:nel) = yscale_sig0*yld0(1:nel)
347 dyld0depsd(1:nel) = yscale_sig0*dyld0depsd(1:nel)
348 ! Compute the Tangent modulus
349 IF (ifunc(3) > 0) THEN
350 ipos(1:nel) = 1
351 iad(1:nel) = npf(ifunc(3)) / 2 + 1
352 ilen(1:nel) = npf(ifunc(3)+1) / 2 - iad(1:nel) - ipos(1:nel)
353 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_tang,dtangdepsd,tang)
354 tang(1:nel) = yscale_tang*tang(1:nel)
355 dtangdepsd(1:nel) = yscale_tang*dtangdepsd(1:nel)
356 ENDIF
357 ! Updating values
358 DO ii=1,nindx
359 i = index(ii)
360 ! Check tangent modulus value
361 IF (tang(i) >= 0.99d0*young(i)) THEN
362 tang(i) = 0.99d0*young(i)
363 dtangdepsd(i) = zero
364 ENDIF
365 ! Yield stress update
366 yld(i) = yld0(i) + (young(i)*tang(i)/(young(i)-tang(i)))*pla(i)
367 ! Yield function value update
368 phi(i) = sigvm(i) - yld(i)
369 ENDDO
370 ENDIF
371 ENDDO
372 ! End of the loop over the iterations
373 ENDIF
374 !=========================================================================
375 ! - END OF PLASTIC CORRECTION WITH CUTTING PLANE (NEWTON-ITERATION) METHOD
376 !=========================================================================
377c
378 ! Storing new values
379 DO i=1,nel
380 ! USR Outputs
381 seq(i) = sigvm(i) ! SIGEQ
382 ! Coefficient for hourglass
383 IF (dpla(i) > zero) THEN
384 et(i) = hardp(i) / (hardp(i) + young(i))
385 ELSE
386 et(i) = one
387 ENDIF
388 ! Computation of the sound speed
389 soundsp(i) = sqrt(a11(i)/rho(i))
390 ! Storing the yield stress
391 sigy(i) = yld(i)
392 ! Thickness variation
393 IF (inloc > 0) THEN
394 IF (loff(i) == one) THEN
395 dezz(i) = -nu*(signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/young(i)
396 dezz(i) = dezz(i) - max(dplanl(i),zero)*half*(signxx(i)+signyy(i))/max(yld(i),em20)
397 ENDIF
398 ELSE
399 dezz(i) = -nu*(signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/young(i) + dezz(i)
400 ENDIF
401 ! Computation of the thickness variation
402 thk(i) = thk(i) + dezz(i)*thkly(i)*off(i)
403 ENDDO
404c
405 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:143