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

Go to the source code of this file.

Functions/Subroutines

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)

Function/Subroutine Documentation

◆ mat115_newton()

subroutine mat115_newton ( integer nel,
integer, dimension(nel), intent(in) ngl,
integer nuparam,
integer nuvar,
intent(in) grho,
time,
timestep,
intent(in) uparam,
intent(inout) uvar,
intent(inout) off,
intent(out) sigy,
intent(in) rho0,
intent(inout) pla,
intent(inout) dpla,
intent(out) soundsp,
intent(out) et,
intent(inout) seq,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depszz,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
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 )

Definition at line 28 of file mat115_newton.F.

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
my_real
#define my_real
Definition cppsort.cpp:32
alpha2
#define alpha2
Definition eval.h:48
alpha
#define alpha
Definition eval.h:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21