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

Go to the source code of this file.

Functions/Subroutines

subroutine mat115_nice (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_nice()

subroutine mat115_nice ( 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_nice.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,K,II,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 :: trdep,dlam,ddep,dphi_dsigvm,dphi_dsigm,
80 . dpxx,dpyy,dpzz,dpxy,dpyz,dpzx,trdfds,
81 . normxx,normyy,normzz,normxy,normyz,normzx,
82 . dphi,sig_dfdsig,dfdsig2,dphi_dpla,
83 . i1,i2,i3,q,r,r_inter,psi,s11,s22,s33
84c
85 my_real, DIMENSION(NEL) ::
86 . sigvm,sigeq,dsigxx,dsigyy,dsigzz,dsigxy,dsigyz,dsigzx,trsig,
87 . sxx,syy,szz,sxy,syz,szx,sigm,j2,yld,hardp,phi0,phi,dpla_dlam,
88 . defvp,dphi_dlam
89 !=======================================================================
90 ! DRUCKER MATERIAL LAW WITH NO DAMAGE
91 !=======================================================================
92 !UVAR(1) PHI YIELD FUNCTION VALUE
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 .NE. 0) 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 ! User inputs
123 phi0(i) = uvar(i,1) ! Previous yield function value
124 ! Standard inputs
125 defvp(i) = zero ! Initialization of the volumic plastic strain increment
126 dpla(i) = zero ! Initialization of the plastic strain increment
127 et(i) = one ! Initialization of hourglass coefficient
128 hardp(i) = zero ! Initialization of hourglass coefficient
129 IF (istat == 0) THEN
130 gamma(i) = uparam(16) ! Yield stress parameter
131 epsd(i) = uparam(17) ! Densification strain
132 alpha2(i) = uparam(18) ! Yield stress parameter
133 beta(i) = uparam(19) ! Yield stress parameter
134 sigp(i) = uparam(20) ! Initial yield stress
135 ELSE
136 sigp(i) = uparam(17) + uparam(18)*((grho(i+nel)/rhof0)**uparam(19))
137 alpha2(i) = uparam(20) + uparam(21)*((grho(i+nel)/rhof0)**uparam(22))
138 gamma(i) = uparam(23) + uparam(24)*((grho(i+nel)/rhof0)**uparam(25))
139 beta(i) = one/(uparam(26) + uparam(27)*((grho(i+nel)/rhof0)**uparam(28)))
140 epsd(i) = (-(nine + (alpha**2))/(three*(alpha**2)))*log(grho(i+nel)/rhof0)
141 ENDIF
142 ENDDO
143c
144 ! Computation of the initial yield stress
145 DO i = 1,nel
146 ! a) - Hardening law
147 yld(i) = sigp(i) + gamma(i)*(pla(i)/epsd(i)) +
148 . alpha2(i)*log(one/max((one - (pla(i)/epsd(i))**beta(i)),em20))
149 ! b) - Checking values
150 yld(i) = max(em10,yld(i))
151 ENDDO
152c
153 !========================================================================
154 ! - COMPUTATION OF TRIAL VALUES
155 !========================================================================
156 DO i=1,nel
157 ! Computation of the trial stress tensor
158 ldav = (depsxx(i) + depsyy(i) + depszz(i)) * lamhook
159 signxx(i) = sigoxx(i) + depsxx(i)*g2 + ldav
160 signyy(i) = sigoyy(i) + depsyy(i)*g2 + ldav
161 signzz(i) = sigozz(i) + depszz(i)*g2 + ldav
162 signxy(i) = sigoxy(i) + depsxy(i)*g
163 signyz(i) = sigoyz(i) + depsyz(i)*g
164 signzx(i) = sigozx(i) + depszx(i)*g
165 ! Computation of the trace of the trial stress tensor
166 trsig(i) = signxx(i) + signyy(i) + signzz(i)
167 sigm(i) = trsig(i) * third
168 ! Computation of the deviatoric trial stress tensor
169 sxx(i) = signxx(i) - sigm(i)
170 syy(i) = signyy(i) - sigm(i)
171 szz(i) = signzz(i) - sigm(i)
172 sxy(i) = signxy(i)
173 syz(i) = signyz(i)
174 szx(i) = signzx(i)
175 ! Second deviatoric invariant and Mises equivalent stress
176 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 )
177 . + sxy(i)**2 + syz(i)**2 + szx(i)**2
178 sigvm(i) = sqrt(three*j2(i))
179 ! Deshpande - Fleck equivalent stress
180 sigeq(i) = sqrt((sigvm(i)**2 + (alpha*sigm(i))**2)/(one + (alpha/three)**2))
181 ENDDO
182c
183 !========================================================================
184 ! - COMPUTATION OF YIELD FONCTION
185 !========================================================================
186 phi(1:nel) = sigeq(1:nel) - yld(1:nel)
187c
188 ! Checking plastic behavior for all elements
189 nindx = 0
190 DO i=1,nel
191 IF ((phi(i)>zero).AND.(off(i) == one)) THEN
192 nindx=nindx+1
193 index(nindx)=i
194 ENDIF
195 ENDDO
196c
197 !====================================================================
198 ! - PLASTIC CORRECTION WITH NICE - EXPLICIT 1 METHOD
199 !====================================================================
200#include "vectorize.inc"
201 DO ii = 1, nindx
202c
203 ! Number of the element with plastic behaviour
204 i = index(ii)
205c
206 ! Computation of the trial stress increment
207 ldav = (depsxx(i) + depsyy(i) + depszz(i)) * lamhook
208 dsigxx(i) = depsxx(i)*g2 + ldav
209 dsigyy(i) = depsyy(i)*g2 + ldav
210 dsigzz(i) = depszz(i)*g2 + ldav
211 dsigxy(i) = depsxy(i)*g
212 dsigyz(i) = depsyz(i)*g
213 dsigzx(i) = depszx(i)*g
214 dlam = zero
215c
216 ! Computation of Von Mises equivalent stress of the previous stress tensor
217 trsig(i) = sigoxx(i) + sigoyy(i) + sigozz(i)
218 sigm(i) = trsig(i) * third
219 sxx(i) = sigoxx(i) - sigm(i)
220 syy(i) = sigoyy(i) - sigm(i)
221 szz(i) = sigozz(i) - sigm(i)
222 sxy(i) = sigoxy(i)
223 syz(i) = sigoyz(i)
224 szx(i) = sigozx(i)
225 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 )
226 . + sxy(i)**2 + syz(i)**2 + szx(i)**2
227 sigvm(i) = sqrt(three*j2(i))
228 ! Deshpande - Fleck equivalent stress
229 sigeq(i) = sqrt((sigvm(i)**2 + (alpha*sigm(i))**2)/(one + (alpha/three)**2))
230
231c
232 ! Note : in this part, the purpose is to compute for each iteration
233 ! a plastic multiplier allowing to update internal variables to satisfy
234 ! the consistency condition.
235 ! Its expression is : DLAMBDA = - (PHI + DPHI) / DPHI_DLAMBDA
236 ! -> PHI : current value of yield function (known)
237 ! -> DPHI : yield function prediction (to compute)
238 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
239 ! into account of internal variables kinetic :
240 ! plasticity, temperature and damage (to compute)
241c
242 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
243 !-------------------------------------------------------------
244c
245 ! Derivative with respect to the Von Mises equivalent stress
246 dphi_dsigvm = sigvm(i)/max((sigeq(i)*(one + (alpha/three)**2)),em20)
247c
248 ! Derivative with respect to the pressure stress
249 dphi_dsigm = (alpha**2)*sigm(i)/max((sigeq(i)*(one + (alpha/three)**2)),em20)
250c
251 ! dPhi/dSig
252 normxx = dphi_dsigvm*three_half*sxx(i)/(max(sigvm(i),em20)) + dphi_dsigm*third
253 normyy = dphi_dsigvm*three_half*syy(i)/(max(sigvm(i),em20)) + dphi_dsigm*third
254 normzz = dphi_dsigvm*three_half*szz(i)/(max(sigvm(i),em20)) + dphi_dsigm*third
255 normxy = two*dphi_dsigvm*three_half*sxy(i)/(max(sigvm(i),em20))
256 normyz = two*dphi_dsigvm*three_half*syz(i)/(max(sigvm(i),em20))
257 normzx = two*dphi_dsigvm*three_half*szx(i)/(max(sigvm(i),em20))
258c
259 ! Restoring previous value of the yield function
260 phi(i) = phi0(i)
261c
262 ! Computation of yield surface trial increment DPHI
263 dphi = normxx * dsigxx(i)
264 . + normyy * dsigyy(i)
265 . + normzz * dsigzz(i)
266 . + normxy * dsigxy(i)
267 . + normyz * dsigyz(i)
268 . + normzx * dsigzx(i)
269c
270 ! 2 - Computation of DPHI_DLAMBDA
271 !---------------------------------------------------------
272c
273 ! a) Derivative with respect stress increments tensor DSIG
274 ! --------------------------------------------------------
275
276 trdfds = normxx + normyy + normzz
277 dfdsig2 = normxx * (normxx*g2 + lamhook*trdfds)
278 . + normyy * (normyy*g2 + lamhook*trdfds)
279 . + normzz * (normzz*g2 + lamhook*trdfds)
280 . + normxy * normxy * g
281 . + normyz * normyz * g
282 . + normzx * normzx * g
283c
284 ! b) Derivatives with respect to plastic strain P
285 ! ------------------------------------------------
286c
287 ! i) Derivative of the yield stress with respect to plastic strain dYLD / dPLA
288 ! ----------------------------------------------------------------------------
289 hardp(i) = (gamma(i)/epsd(i)) +
290 . alpha2(i)*((one-(beta(i)/epsd(i))*(pla(i)/epsd(i))**(beta(i)-1))/
291 . max((one - (pla(i)/epsd(i))**beta(i)),em20))
292 dphi_dpla = hardp(i)
293c
294 ! ii) Derivative of dPLA with respect to DLAM
295 ! -------------------------------------------
296 sig_dfdsig = signxx(i) * normxx
297 . + signyy(i) * normyy
298 . + signzz(i) * normzz
299 . + signxy(i) * normxy
300 . + signyz(i) * normyz
301 . + signzx(i) * normzx
302 dpla_dlam(i) = sig_dfdsig / yld(i)
303c
304 ! d) Derivative with respect to the yield stress
305 ! ----------------------------------------------
306c
307 ! e) Derivative of PHI with respect to DLAM ( = -DENOM)
308 dphi_dlam(i) = - dfdsig2 - dphi_dpla*dpla_dlam(i)
309 dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20) ,dphi_dlam(i))
310c
311 ! 3 - Computation of plastic multiplier and variables update
312 !----------------------------------------------------------
313c
314 ! Computation of the plastic multiplier
315 dlam = -(dphi + phi(i)) / dphi_dlam(i)
316 dlam = max(dlam, zero)
317c
318 ! Plastic strains tensor update
319 dpxx = dlam * normxx
320 dpyy = dlam * normyy
321 dpzz = dlam * normzz
322 dpxy = dlam * normxy
323 dpyz = dlam * normyz
324 dpzx = dlam * normzx
325 trdep = dpxx + dpyy + dpzz
326c
327 ! Cumulated plastic strain and strain rate update
328 ddep = (dlam/yld(i))*sig_dfdsig
329 dpla(i) = dpla(i) + ddep
330 pla(i) = pla(i) + ddep
331 defvp(i) = defvp(i) + trdep
332c
333 ! Elasto-plastic stresses update
334 ldav = trdep * lamhook
335 signxx(i) = sigoxx(i) + dsigxx(i) - (dpxx*g2 + ldav)
336 signyy(i) = sigoyy(i) + dsigyy(i) - (dpyy*g2 + ldav)
337 signzz(i) = sigozz(i) + dsigzz(i) - (dpzz*g2 + ldav)
338 signxy(i) = sigoxy(i) + dsigxy(i) - dpxy*g
339 signyz(i) = sigoyz(i) + dsigyz(i) - dpyz*g
340 signzx(i) = sigozx(i) + dsigzx(i) - dpzx*g
341 trsig(i) = signxx(i) + signyy(i) + signzz(i)
342 sigm(i) = trsig(i) * third
343 sxx(i) = signxx(i) - sigm(i)
344 syy(i) = signyy(i) - sigm(i)
345 szz(i) = signzz(i) - sigm(i)
346 sxy(i) = signxy(i)
347 syz(i) = signyz(i)
348 szx(i) = signzx(i)
349 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 )
350 . + sxy(i)**2 + syz(i)**2 + szx(i)**2
351 sigvm(i) = sqrt(three*j2(i))
352 ! Deshpande - Fleck equivalent stress
353 sigeq(i) = sqrt((sigvm(i)**2 + (alpha*sigm(i))**2)/(one + (alpha/three)**2))
354c
355 ! Yield stress update
356 yld(i) = sigp(i) + gamma(i)*(pla(i)/epsd(i)) +
357 . alpha2(i)*log(one/max((one - (pla(i)/epsd(i))**beta(i)),em20))
358 yld(i) = max(yld(i),em10)
359c
360 ! Yield function value update
361 phi(i) = sigeq(i) - yld(i)
362 ENDDO
363 !===================================================================
364 ! - END OF PLASTIC CORRECTION WITH NICE - EXPLICIT 1 METHOD
365 !===================================================================
366c
367 ! Storing new values
368 DO i=1,nel
369 ! USR Outputs
370 IF (dpla(i) > zero) THEN
371 uvar(i,1) = phi(i) ! Yield function value
372 ELSE
373 uvar(i,1) = zero
374 ENDIF
375 IF (cfail > zero) THEN
376 uvar(i,2) = uvar(i,2) + defvp(i)
377 IF (uvar(i,2) >= cfail) THEN
378 IF (off(i) == one) off(i) = four_over_5
379 ENDIF
380 ENDIF
381 IF (pfail > zero) THEN
382 ! Computing the principal stresses
383 i1 = signxx(i)+signyy(i)+signzz(i)
384 i2 = signxx(i)*signyy(i)+signyy(i)*signzz(i)+signzz(i)*signxx(i)-
385 . signxy(i)*signxy(i)-signzx(i)*signzx(i)-signyz(i)*signyz(i)
386 i3 = signxx(i)*signyy(i)*signzz(i)-signxx(i)*signyz(i)*signyz(i)-
387 . signyy(i)*signzx(i)*signzx(i)-signzz(i)*signxy(i)*signxy(i)+
388 . two*signxy(i)*signzx(i)*signyz(i)
389 q = (three*i2 - i1*i1)/nine
390 r = (two*i1*i1*i1-nine*i1*i2+twenty7*i3)/cinquante4 ! (2*I3^3-9*I1*I2+27*I3)/54
391 r_inter = min(r/sqrt(max(em20,(-q**3))),one)
392 psi = acos(max(r_inter,-one))
393
394 s11 = two*sqrt(-q)*cos(psi/three)+third*i1
395 s22 = two*sqrt(-q)*cos((psi+two*pi)/three)+third*i1
396 s33 = two*sqrt(-q)*cos((psi+four*pi)/three)+third*i1
397 s11 = max(s11,s22,s33)
398 ! Failure criterion
399 IF (s11 >= pfail) THEN
400 IF (off(i) == one) off(i) = four_over_5
401 ENDIF
402 ENDIF
403 seq(i) = sigeq(i) ! SIGEQ
404 ! Coefficient for hourglass
405 IF (dpla(i) > zero) THEN
406 et(i) = hardp(i) / (hardp(i) + young)
407 ELSE
408 et(i) = one
409 ENDIF
410 ! Computation of the sound speed
411 soundsp(i) = sqrt((bulk + four_over_3*g)/grho(nel+i))
412 ! Storing the yield stress
413 sigy(i) = yld(i)
414 ENDDO
415c
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