OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat115_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!|| mat115_newton ../engine/source/materials/mat/mat115/mat115_newton.F
25!||--- called by ------------------------------------------------------
26!|| sigeps115 ../engine/source/materials/mat/mat115/sigeps115.F
27!||====================================================================
28 SUBROUTINE mat115_newton(
29 1 NEL ,NGL ,NUPARAM ,NUVAR ,GRHO ,
30 2 TIME ,TIMESTEP,UPARAM ,UVAR ,OFF ,SIGY ,
31 3 RHO0 ,PLA ,DPLA ,SOUNDSP ,ET ,SEQ ,
32 4 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
33 5 SIGOXX ,SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
34 6 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX )
35 !=======================================================================
36 ! Implicit types
37 !=======================================================================
38#include "implicit_f.inc"
39 !=======================================================================
40 ! Common
41 !=======================================================================
42 !=======================================================================
43 ! Dummy arguments
44 !=======================================================================
45 INTEGER NEL,NUPARAM,NUVAR
46 INTEGER ,DIMENSION(NEL), INTENT(IN) :: NGL
47 my_real
48 . TIME,TIMESTEP
49 my_real,DIMENSION(NUPARAM), INTENT(IN) ::
50 . UPARAM
51 my_real,DIMENSION(2*NEL), INTENT(IN) ::
52 . grho
53 my_real,DIMENSION(NEL), INTENT(IN) ::
54 . rho0,
55 . depsxx,depsyy,depszz,depsxy,depsyz,depszx,
56 . sigoxx,sigoyy,sigozz,sigoxy,sigoyz,sigozx
57c
58 my_real ,DIMENSION(NEL), INTENT(OUT) ::
59 . soundsp,sigy,et,
60 . signxx,signyy,signzz,signxy,signyz,signzx
61c
62 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
63 . pla,dpla,off,seq
64 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) ::
65 . uvar
66c
67 !=======================================================================
68 ! Local Variables
69 !=======================================================================
70 INTEGER I,II,ITER,NITER,NINDX,INDEX(NEL),Istat
71c
72 my_real
73 . YOUNG,BULK,LAMHOOK,G,G2,NU,ALPHA,CFAIL,PFAIL,RHOF0
74c
75 my_real, DIMENSION(NEL) ::
76 . gamma,epsd,alpha2,beta,sigp
77c
78 real(kind=8) :: ldav
79 my_real
80 . h,trdep,dlam,ddep,dphi_dsigvm,dphi_dsigm,
81 . dpxx,dpyy,dpzz,dpxy,dpyz,dpzx,trdfds,
82 . normxx,normyy,normzz,normxy,normyz,normzx,
83 . dphi,sig_dfdsig,dfdsig2,dphi_dpla,
84 . i1,i2,i3,q,r,r_inter,psi,s11,s22,s33
85c
86 my_real, DIMENSION(NEL) ::
87 . sigvm,sigeq,dsigxx,dsigyy,dsigzz,dsigxy,dsigyz,dsigzx,trsig,
88 . sxx,syy,szz,sxy,syz,szx,sigm,j2,yld,hardp,phi,dpla_dlam,dphi_dlam,
89 . defvp
90 !=======================================================================
91 ! DRUCKER MATERIAL LAW WITH NO DAMAGE
92 !=======================================================================
93 !-
94 !DEPIJ PLASTIC STRAIN TENSOR COMPONENT
95 !DEPSIJ TOTAL STRAIN TENSOR COMPONENT (EL+PL)
96 !=======================================================================
97c
98 !=======================================================================
99 ! - INITIALISATION OF COMPUTATION ON TIME STEP
100 !=======================================================================
101 ! Recovering model parameters
102 ! Elastic parameters
103 young = uparam(1) ! Young modulus
104 bulk = uparam(2) ! Bulk modulus
105 g = uparam(3) ! Shear modulus
106 g2 = uparam(4) ! 2*Shear modulus
107 lamhook = uparam(5) ! Lambda Hooke parameter
108 nu = uparam(6) ! Poisson ration
109 ! Statistic variation flag
110 istat = uparam(12)
111 ! Plastic criterion and hardening parameters
112 alpha = uparam(13) ! Yield function shape parameter
113 cfail = uparam(14) ! Tensile volumic strain at failure
114 pfail = uparam(15) ! Principal stress at failure
115 ! Density of base material
116 IF (istat == 1) rhof0 = uparam(16)
117c
118 ! Recovering internal variables
119 DO i=1,nel
120 IF (off(i) < 0.1) off(i) = zero
121 IF (off(i) < one) off(i) = off(i)*four_over_5
122 ! Standard inputs
123 defvp(i) = zero ! Initialization of the volumic plastic strain increment
124 dpla(i) = zero ! Initialization of the plastic strain increment
125 et(i) = one ! Initialization of hourglass coefficient
126 hardp(i) = zero ! Initialization of hourglass coefficient
127 IF (istat == 0) THEN
128 gamma(i) = uparam(16) ! Yield stress parameter
129 epsd(i) = uparam(17) ! Densification strain
130 alpha2(i) = uparam(18) ! Yield stress parameter
131 beta(i) = uparam(19) ! Yield stress parameter
132 sigp(i) = uparam(20) ! Initial yield stress
133 ELSE
134 sigp(i) = uparam(17) + uparam(18)*((grho(i+nel)/rhof0)**uparam(19))
135 alpha2(i) = uparam(20) + uparam(21)*((grho(i+nel)/rhof0)**uparam(22))
136 gamma(i) = uparam(23) + uparam(24)*((grho(i+nel)/rhof0)**uparam(25))
137 beta(i) = one/(uparam(26) + uparam(27)*((grho(i+nel)/rhof0)**uparam(28)))
138 epsd(i) = (-(nine + (alpha**2))/(three*(alpha**2)))*log(grho(i+nel)/rhof0)
139 ENDIF
140 ENDDO
141c
142 ! Computation of the initial yield stress
143 DO i = 1,nel
144 ! a) - hardening law
145 yld(i) = sigp(i) + gamma(i)*(pla(i)/epsd(i)) +
146 . alpha2(i)*log(one/max((one - (pla(i)/epsd(i))**beta(i)),em20))
147 ! b) - Checking values
148 yld(i) = max(em10,yld(i))
149 ENDDO
150c
151 !========================================================================
152 ! - COMPUTATION OF TRIAL VALUES
153 !========================================================================
154 DO i=1,nel
155 ! Computation of the trial stress tensor
156 ldav = (depsxx(i) + depsyy(i) + depszz(i)) * lamhook
157 signxx(i) = sigoxx(i) + depsxx(i)*g2 + ldav
158 signyy(i) = sigoyy(i) + depsyy(i)*g2 + ldav
159 signzz(i) = sigozz(i) + depszz(i)*g2 + ldav
160 signxy(i) = sigoxy(i) + depsxy(i)*g
161 signyz(i) = sigoyz(i) + depsyz(i)*g
162 signzx(i) = sigozx(i) + depszx(i)*g
163 ! Computation of the trace of the trial stress tensor
164 trsig(i) = signxx(i) + signyy(i) + signzz(i)
165 sigm(i) = trsig(i) * third
166 ! Computation of the deviatoric trial stress tensor
167 sxx(i) = signxx(i) - sigm(i)
168 syy(i) = signyy(i) - sigm(i)
169 szz(i) = signzz(i) - sigm(i)
170 sxy(i) = signxy(i)
171 syz(i) = signyz(i)
172 szx(i) = signzx(i)
173 ! Second deviatoric invariant and Mises equivalent stress
174 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 )
175 . + sxy(i)**2 + syz(i)**2 + szx(i)**2
176 sigvm(i) = sqrt(three*j2(i))
177 ! Deshpande - Fleck equivalent stress
178 sigeq(i) = sqrt((sigvm(i)**2 + (alpha*sigm(i))**2)/(one + (alpha/three)**2))
179 ENDDO
180c
181 !========================================================================
182 ! - COMPUTATION OF YIELD FONCTION
183 !========================================================================
184 phi(1:nel) = sigeq(1:nel) - yld(1:nel)
185c
186 ! Checking plastic behavior for all elements
187 nindx = 0
188 DO i=1,nel
189 IF ((phi(i)>zero).AND.(off(i) == one)) THEN
190 nindx=nindx+1
191 index(nindx)=i
192 ENDIF
193 ENDDO
194c
195 !====================================================================
196 ! - PLASTIC CORRECTION WITH BACKWARD EULER METHOD (NEWTON RESOLUTION)
197 !====================================================================
198c
199 ! Number of maximum Newton iterations
200 niter = 3
201c
202 ! Loop over the iterations
203 DO iter = 1, niter
204#include "vectorize.inc"
205 ! Loop over yielding elements
206 DO ii=1,nindx
207 i = index(ii)
208c
209 ! Note : in this part, the purpose is to compute for each iteration
210 ! a plastic multiplier allowing to update internal variables to satisfy
211 ! the consistency condition using the backward Euler implicit method
212 ! with a Newton-Raphson iterative procedure
213 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
214 ! -> PHI : current value of yield function (known)
215 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
216 ! into account of internal variables kinetic :
217 ! plasticity, temperature and damage (to compute)
218c
219 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
220 !-------------------------------------------------------------
221c
222 ! Derivative with respect to the Von Mises equivalent stress
223 dphi_dsigvm = sigvm(i)/max((sigeq(i)*(one + (alpha/three)**2)),em20)
224c
225 ! Derivative with respect to the pressure stress
226 dphi_dsigm = (alpha**2)*sigm(i)/max((sigeq(i)*(one + (alpha/three)**2)),em20)
227c
228 ! dPhi/dSig
229 normxx = dphi_dsigvm*three_half*sxx(i)/(max(sigvm(i),em20)) + dphi_dsigm*third
230 normyy = dphi_dsigvm*three_half*syy(i)/(max(sigvm(i),em20)) + dphi_dsigm*third
231 normzz = dphi_dsigvm*three_half*szz(i)/(max(sigvm(i),em20)) + dphi_dsigm*third
232 normxy = two*dphi_dsigvm*three_half*sxy(i)/(max(sigvm(i),em20))
233 normyz = two*dphi_dsigvm*three_half*syz(i)/(max(sigvm(i),em20))
234 normzx = two*dphi_dsigvm*three_half*szx(i)/(max(sigvm(i),em20))
235c
236 ! 2 - Computation of DPHI_DLAMBDA
237 !---------------------------------------------------------
238c
239 ! a) Derivative with respect stress increments tensor DSIG
240 ! --------------------------------------------------------
241 trdfds = normxx + normyy + normzz
242 dfdsig2 = normxx * (normxx*g2 + lamhook*trdfds)
243 . + normyy * (normyy*g2 + lamhook*trdfds)
244 . + normzz * (normzz*g2 + lamhook*trdfds)
245 . + normxy * normxy * g
246 . + normyz * normyz * g
247 . + normzx * normzx * g
248c
249 ! b) Derivatives with respect to plastic strain P
250 ! ------------------------------------------------
251c
252 ! i) Derivative of the yield function with respect to plastic strain dPHI / dPLA
253 ! ----------------------------------------------------------------------------
254 hardp(i) = (gamma(i)/epsd(i)) +
255 . alpha2(i)*((one-(beta(i)/epsd(i))*(pla(i)/epsd(i))**(beta(i)-1))/
256 . max((one - (pla(i)/epsd(i))**beta(i)),em20))
257 dphi_dpla = hardp(i)
258c
259 ! ii) Derivative of dPLA with respect to DLAM
260 ! -------------------------------------------
261 sig_dfdsig = signxx(i) * normxx
262 . + signyy(i) * normyy
263 . + signzz(i) * normzz
264 . + signxy(i) * normxy
265 . + signyz(i) * normyz
266 . + signzx(i) * normzx
267 dpla_dlam(i) = sig_dfdsig / yld(i)
268c
269 ! 3 - Computation of plastic multiplier and variables update
270 !----------------------------------------------------------
271c
272 ! Derivative of PHI with respect to DLAM
273 dphi_dlam(i) = - dfdsig2 - dphi_dpla*dpla_dlam(i)
274 dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20) ,dphi_dlam(i))
275c
276 ! Computation of the plastic multiplier
277 dlam = -phi(i)/dphi_dlam(i)
278c
279 ! Plastic strains tensor update
280 dpxx = dlam * normxx
281 dpyy = dlam * normyy
282 dpzz = dlam * normzz
283 dpxy = dlam * normxy
284 dpyz = dlam * normyz
285 dpzx = dlam * normzx
286 trdep = dpxx + dpyy + dpzz
287c
288 ! Cumulated plastic strain and strain rate update
289 ddep = (dlam/yld(i))*sig_dfdsig
290 dpla(i) = max(zero, dpla(i) + ddep)
291 pla(i) = pla(i) + ddep
292 defvp(i) = defvp(i) + trdep
293c
294 ! Elasto-plastic stresses update
295 ldav = trdep * lamhook
296 signxx(i) = signxx(i) - (dpxx*g2 + ldav)
297 signyy(i) = signyy(i) - (dpyy*g2 + ldav)
298 signzz(i) = signzz(i) - (dpzz*g2 + ldav)
299 signxy(i) = signxy(i) - dpxy*g
300 signyz(i) = signyz(i) - dpyz*g
301 signzx(i) = signzx(i) - dpzx*g
302 trsig(i) = signxx(i) + signyy(i) + signzz(i)
303 sigm(i) = trsig(i) * third
304 sxx(i) = signxx(i) - sigm(i)
305 syy(i) = signyy(i) - sigm(i)
306 szz(i) = signzz(i) - sigm(i)
307 sxy(i) = signxy(i)
308 syz(i) = signyz(i)
309 szx(i) = signzx(i)
310 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 )
311 . + sxy(i)**2 + syz(i)**2 + szx(i)**2
312 sigvm(i) = sqrt(three*j2(i))
313 ! Deshpande - Fleck equivalent stress
314 sigeq(i) = sqrt((sigvm(i)**2 + (alpha*sigm(i))**2)/(one + (alpha/three)**2))
315c
316 ! Yield stress update
317 yld(i) = sigp(i) + gamma(i)*(pla(i)/epsd(i)) +
318 . alpha2(i)*log(one/max((one - (pla(i)/epsd(i))**beta(i)),em20))
319 yld(i) = max(yld(i),em10)
320c
321 ! Yield function value update
322 phi(i) = sigeq(i) - yld(i)
323c
324 ENDDO
325 ! End of the loop over the iterations
326 ENDDO
327 ! End of the loop over the iterations
328 !===================================================================
329 ! - END OF PLASTIC CORRECTION WITH NEWTON IMPLICIT ITERATIVE METHOD
330 !===================================================================
331c
332 ! Storing new values
333 DO i=1,nel
334 ! USR Outputs
335 IF (cfail > zero) THEN
336 uvar(i,1) = uvar(i,1) + defvp(i)
337 IF (uvar(i,1) >= cfail) THEN
338 IF (off(i) == one) off(i) = four_over_5
339 ENDIF
340 ENDIF
341 IF (pfail > zero) THEN
342 ! Computing the principal stresses
343 i1 = signxx(i)+signyy(i)+signzz(i)
344 i2 = signxx(i)*signyy(i)+signyy(i)*signzz(i)+signzz(i)*signxx(i)-
345 . signxy(i)*signxy(i)-signzx(i)*signzx(i)-signyz(i)*signyz(i)
346 i3 = signxx(i)*signyy(i)*signzz(i)-signxx(i)*signyz(i)*signyz(i)-
347 . signyy(i)*signzx(i)*signzx(i)-signzz(i)*signxy(i)*signxy(i)+
348 . two*signxy(i)*signzx(i)*signyz(i)
349 q = (three*i2 - i1*i1)/nine
350 r = (two*i1*i1*i1-nine*i1*i2+twenty7*i3)/cinquante4 ! (2*I3^3-9*I1*I2+27*I3)/54
351 r_inter = min(r/sqrt(max(em20,(-q**3))),one)
352 psi = acos(max(r_inter,-one))
353 s11 = two*sqrt(-q)*cos(psi/three)+third*i1
354 s22 = two*sqrt(-q)*cos((psi+two*pi)/three)+third*i1
355 s33 = two*sqrt(-q)*cos((psi+four*pi)/three)+third*i1
356 s11 = max(s11,s22,s33)
357 ! Failure criterion
358 IF (s11 >= pfail) THEN
359 IF (off(i) == one) off(i) = four_over_5
360 ENDIF
361 ENDIF
362 seq(i) = sigeq(i) ! SIGEQ
363 ! Coefficient for hourglass
364 IF (dpla(i) > zero) THEN
365 et(i) = hardp(i) / (hardp(i) + young)
366 ELSE
367 et(i) = one
368 ENDIF
369 ! Computation of the sound speed
370 soundsp(i) = sqrt((bulk + four_over_3*g)/grho(nel+i))
371 ! Storing the yield stress
372 sigy(i) = yld(i)
373 ENDDO
374c
375 END
alpha2
#define alpha2
Definition eval.h:48
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
mat115_newton
subroutine mat115_newton(nel, ngl, nuparam, nuvar, grho, time, timestep, uparam, uvar, off, sigy, rho0, pla, dpla, soundsp, et, seq, depsxx, depsyy, depszz, depsxy, depsyz, depszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx)
Definition mat115_newton.F:35