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

Go to the source code of this file.

Functions/Subroutines

subroutine mat121_nice (nel, ngl, nuparam, nuvar, nfunc, ifunc, npf, tf, timestep, time, uparam, uvar, rho, pla, dpla, soundsp, epsd, off, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigy, et, seq)

Function/Subroutine Documentation

◆ mat121_nice()

subroutine mat121_nice ( integer nel,
integer, dimension(nel), intent(in) ngl,
integer nuparam,
integer nuvar,
integer nfunc,
integer, dimension(nfunc) ifunc,
integer, dimension(*) npf,
tf,
timestep,
time,
intent(in) uparam,
intent(inout) uvar,
intent(in) rho,
intent(inout) pla,
intent(inout) dpla,
intent(out) soundsp,
intent(inout) epsd,
intent(inout) off,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depszz,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
intent(in) epspxx,
intent(in) epspyy,
intent(in) epspzz,
intent(in) epspxy,
intent(in) epspyz,
intent(in) epspzx,
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,
intent(out) sigy,
intent(out) et,
intent(inout) seq )

Definition at line 30 of file mat121_nice.F.

39 !=======================================================================
40 ! Implicit types
41 !=======================================================================
42#include "implicit_f.inc"
43 !=======================================================================
44 ! Common
45 !=======================================================================
46 !=======================================================================
47 ! dummy arguments
48 !=======================================================================
49 INTEGER NEL,NUPARAM,NUVAR,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,
57 . depsxx,depsyy,depszz,depsxy,depsyz,depszx,
58 . epspxx,epspyy,epspzz,epspxy,epspyz,epspzx,
59 . sigoxx,sigoyy,sigozz,sigoxy,sigoyz,sigozx
60 my_real ,DIMENSION(NEL), INTENT(OUT) ::
61 . soundsp,sigy,et,
62 . signxx,signyy,signzz,signxy,signyz,signzx
63 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
64 . pla,dpla,epsd,off,seq
65 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) ::
66 . uvar
67 !=======================================================================
68 ! Local Variables
69 !=======================================================================
70 INTEGER I,II,Ivisc,NINDX,INDEX(NEL),IPOS(NEL),IAD(NEL),ILEN(NEL)
71 my_real
72 . young(nel),bulk(nel),g(nel),nu,nnu,tang(nel),lam(nel),g2(nel),
73 . afiltr,xscale_sig0,yscale_sig0,xscale_youn,yscale_youn,
74 . xscale_tang,yscale_tang
75 my_real
76 . dpdt,dlam,ddep,depxx,depyy,devepspxx,devepspyy,devepspzz,trepsp,ldav,
77 . normxx,normyy,normzz,normxy,normyz,normzx,denom,dfdsig2,dpdt_nl,depsdt,
78 . dphi,dtinv
79 my_real, DIMENSION(NEL) ::
80 . sxx,syy,szz,sxy,syz,szx,sigvm,yld,hardp,phi,dezz,yld0,dyld0depsd,
81 . dyoundepsd,dtangdepsd,trsig,dphi_dlam,test,dpxx,dpyy,dpxy,
82 . dpzz,dpyz,dpzx,dsigxx,dsigyy,dsigzz,dsigxy,dsigyz,dsigzx,phi0
83c
84 !=======================================================================
85 ! - INITIALISATION OF COMPUTATION ON TIME STEP
86 !=======================================================================
87 ! Recovering model parameters
88 ! Elastic parameters
89 young(1:nel) = uparam(1) ! Young modulus
90 bulk(1:nel) = uparam(2) ! Bulk modulus
91 g(1:nel) = uparam(3) ! Shear modulus
92 g2(1:nel) = uparam(4) ! 2*Shear modulus
93 lam(1:nel) = uparam(5) ! Lame coefficient
94 nu = uparam(6) ! Poisson ration
95 nnu = uparam(7) ! NU/(1-NU)
96 ! Flags for computation
97 ivisc = nint(uparam(12)) ! Viscosity formulation
98 ! Strain-rate filtering (if Ivisc = 0)
99 afiltr = min(one, uparam(14)*timestep)
100 ! Initial yield stress vs strain-rate curve
101 IF (ifunc(1) > 0) THEN
102 xscale_sig0 = uparam(16) ! Strain-rate scale factor
103 yscale_sig0 = uparam(17) ! Initial yield stress scale factor
104 yld0(1:nel) = zero
105 ELSE
106 yld0(1:nel) = uparam(17) ! Constant yield stress
107 dyld0depsd(1:nel) = zero
108 ENDIF
109 ! Young modulus vs strain-rate curve
110 xscale_youn = uparam(18) ! Strain-rate scale factor
111 yscale_youn = uparam(19) ! Young modulus scale factor
112 ! Tangent modulus vs strain-rate curve
113 IF (ifunc(3) > 0) THEN
114 xscale_tang = uparam(20) ! Strain-rate scale factor
115 yscale_tang = uparam(21) ! Tangent modulus scale factor
116 tang(1:nel) = zero
117 ELSE
118 tang(1:nel) = uparam(21) ! Constant tangent modulus
119 dtangdepsd(1:nel) = zero
120 ENDIF
121 dtinv = one/max(timestep,em20) ! Inverse of timestep
122c
123 ! Recovering internal variables
124 DO i=1,nel
125 IF (off(i) < em01) off(i) = zero
126 IF (off(i) < one) off(i) = off(i)*four_over_5
127 ! User inputs
128 phi0(i) = uvar(i,1) ! Previous yield function value
129 ! Standard inputs
130 dpla(i) = zero ! Initialization of the plastic strain increment
131 et(i) = one ! Initialization of hourglass coefficient
132 hardp(i) = zero ! Initialization of hourglass coefficient
133 dpxx(i) = zero ! Initialization of the XX plastic strain increment
134 dpyy(i) = zero ! Initialization of the YY plastic strain increment
135 dpzz(i) = zero ! Initialization of the ZZ plastic strain increment
136 dpxy(i) = zero ! Initialization of the XY plastic strain increment
137 dpyz(i) = zero ! Initialization of the YZ plastic strain increment
138 dpzx(i) = zero ! Initialization of the ZX plastic strain increment
139 dyld0depsd(i) = zero ! Initialization of the derivative of SIG0
140 dyoundepsd(i) = zero ! Initialization of the derivative of YOUN
141 dtangdepsd(i) = zero ! Initialization of the derivative of TANG
142 ENDDO
143c
144 ! Filling the strain-rate vector
145 IF (ivisc == 0) THEN
146 ! Compute effective strain-rate
147 DO i = 1,nel
148 trepsp = third*(epspxx(i) + epspyy(i) + epspzz(i))
149 devepspxx = epspxx(i) - trepsp
150 devepspyy = epspyy(i) - trepsp
151 devepspzz = epspzz(i) - trepsp
152 depsdt = two_third*(devepspxx**2 + devepspyy**2 + devepspzz**2 +
153 . two*(epspxy(i)**2) + two*(epspyz(i)**2) +
154 . two*(epspzx(i)**2))
155 depsdt = sqrt(max(depsdt,zero))
156 epsd(i) = afiltr * depsdt + (one - afiltr) * epsd(i)
157 ENDDO
158 ELSE
159 ! Reset plastic strain-rate
160 epsd(1:nel) = zero
161 ENDIF
162c
163 ! Compute the initial yield stress
164 IF (ifunc(1) > 0) THEN
165 ipos(1:nel) = 1
166 iad(1:nel) = npf(ifunc(1)) / 2 + 1
167 ilen(1:nel) = npf(ifunc(1)+1) / 2 - iad(1:nel) - ipos(1:nel)
168 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_sig0,dyld0depsd,yld0)
169 yld0(1:nel) = yscale_sig0*yld0(1:nel)
170 dyld0depsd(1:nel) = yscale_sig0*dyld0depsd(1:nel)
171 ENDIF
172 IF (ivisc == 1) THEN
173 IF (uvar(1,2) == zero) THEN
174 dyld0depsd(1:nel) = yscale_sig0*dyld0depsd(1:nel)
175 ELSE
176 dyld0depsd(1:nel) = uvar(1:nel,2)
177 ENDIF
178 ENDIF
179 ! Compute the Young modulus
180 IF (ifunc(2) > 0) THEN
181 ipos(1:nel) = 1
182 iad(1:nel) = npf(ifunc(2)) / 2 + 1
183 ilen(1:nel) = npf(ifunc(2)+1) / 2 - iad(1:nel) - ipos(1:nel)
184 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_youn,dyoundepsd,young)
185 young(1:nel) = yscale_youn*young(1:nel)
186 g(1:nel) = half * young(1:nel) / (one + nu)
187 g2(1:nel) = young(1:nel) / (one + nu)
188 bulk(1:nel) = third * young(1:nel) / (one - nu*two)
189 lam(1:nel) = g2(1:nel) * nu /(one - two*nu)
190 ENDIF
191 ! Compute the Tangent modulus
192 IF (ifunc(3) > 0) THEN
193 ipos(1:nel) = 1
194 iad(1:nel) = npf(ifunc(3)) / 2 + 1
195 ilen(1:nel) = npf(ifunc(3)+1) / 2 - iad(1:nel) - ipos(1:nel)
196 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_tang,dtangdepsd,tang)
197 tang(1:nel) = yscale_tang*tang(1:nel)
198 IF (ivisc == 1) THEN
199 IF (uvar(1,3) == zero) THEN
200 dtangdepsd(1:nel) = yscale_tang*dtangdepsd(1:nel)
201 ELSE
202 dtangdepsd(1:nel) = uvar(1:nel,3)
203 ENDIF
204 ENDIF
205 ENDIF
206 ! Check tangent modulus value + Assembling the yield stress
207 DO i = 1,nel
208 IF (tang(i) >= 0.99d0*young(i)) THEN
209 tang(i) = 0.99d0*young(i)
210 ENDIF
211 yld(i) = yld0(i) + (young(i)*tang(i)/(young(i)-tang(i)))*pla(i)
212 ENDDO
213c
214 !========================================================================
215 ! - COMPUTATION OF TRIAL VALUES
216 !========================================================================
217 DO i=1,nel
218 ! Computation of the trial stress tensor
219 ldav = (depsxx(i) + depsyy(i) + depszz(i)) * lam(i)
220 signxx(i) = sigoxx(i) + depsxx(i)*g2(i) + ldav
221 signyy(i) = sigoyy(i) + depsyy(i)*g2(i) + ldav
222 signzz(i) = sigozz(i) + depszz(i)*g2(i) + ldav
223 signxy(i) = sigoxy(i) + depsxy(i)*g(i)
224 signyz(i) = sigoyz(i) + depsyz(i)*g(i)
225 signzx(i) = sigozx(i) + depszx(i)*g(i)
226 ! Computation of the trace of the trial stress tensor
227 trsig(i) = signxx(i) + signyy(i) + signzz(i)
228 ! Computation of the deviatoric trial stress tensor
229 sxx(i) = signxx(i) - trsig(i) * third
230 syy(i) = signyy(i) - trsig(i) * third
231 szz(i) = signzz(i) - trsig(i) * third
232 sxy(i) = signxy(i)
233 syz(i) = signyz(i)
234 szx(i) = signzx(i)
235 ! Von Mises equivalent stress
236 sigvm(i) = three_half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + three*sxy(i)**2
237 . + three*syz(i)**2 + three*szx(i)**2
238 sigvm(i) = sqrt(sigvm(i))
239 ENDDO
240c
241 !========================================================================
242 ! - COMPUTATION OF YIELD FONCTION
243 !========================================================================
244 phi(1:nel) = sigvm(1:nel) - yld(1:nel)
245c
246 ! Checking plastic behavior for all elements
247 nindx = 0
248 DO i=1,nel
249 IF (phi(i) > zero .AND. off(i) == one) THEN
250 nindx=nindx+1
251 index(nindx)=i
252 ENDIF
253 ENDDO
254c
255 !====================================================================
256 ! - PLASTIC CORRECTION WITH NICE - EXPLICIT METHOD
257 !====================================================================
258 IF (nindx > 0) THEN
259#include "vectorize.inc"
260 ! Loop over yielding elements
261 DO ii=1,nindx
262c
263 ! Number of the element with plastic behaviour
264 i = index(ii)
265c
266 ! Computation of the trial stress increment
267 ldav = (depsxx(i) + depsyy(i) + depszz(i)) * lam(i)
268 dsigxx(i) = depsxx(i)*g2(i) + ldav
269 dsigyy(i) = depsyy(i)*g2(i) + ldav
270 dsigzz(i) = depszz(i)*g2(i) + ldav
271 dsigxy(i) = depsxy(i)*g(i)
272 dsigyz(i) = depsyz(i)*g(i)
273 dsigzx(i) = depszx(i)*g(i)
274c
275 ! Computation of the old Von Mises stress
276 trsig(i) = sigoxx(i) + sigoyy(i) + sigozz(i)
277 sxx(i) = sigoxx(i) - trsig(i) * third
278 syy(i) = sigoyy(i) - trsig(i) * third
279 szz(i) = sigozz(i) - trsig(i) * third
280 sxy(i) = sigoxy(i)
281 syz(i) = sigoyz(i)
282 szx(i) = sigozx(i)
283 sigvm(i) = three_half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + three*sxy(i)**2
284 . + three*syz(i)**2 + three*szx(i)**2
285 sigvm(i) = sqrt(sigvm(i))
286c
287 ! Note : in this part, the purpose is to compute for each iteration
288 ! a plastic multiplier allowing to update internal variables to satisfy
289 ! the consistency condition.
290 ! Its expression is : DLAMBDA = - (PHI + DPHI) / DPHI_DLAMBDA
291 ! -> PHI : current value of yield function (known)
292 ! -> DPHI : yield function prediction (to compute)
293 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
294 ! into account of internal variables kinetic :
295 ! plasticity, temperature and damage (to compute)
296c
297 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
298 !-------------------------------------------------------------
299 normxx = three_half*sxx(i)/sigvm(i)
300 normyy = three_half*syy(i)/sigvm(i)
301 normzz = three_half*szz(i)/sigvm(i)
302 normxy = three*sxy(i)/sigvm(i)
303 normyz = three*syz(i)/sigvm(i)
304 normzx = three*szx(i)/sigvm(i)
305c
306 ! Restoring previous value of the yield function
307 phi(i) = phi0(i)
308c
309 ! Computation of yield surface trial increment DPHI
310 dphi = normxx * dsigxx(i)
311 . + normyy * dsigyy(i)
312 . + normzz * dsigzz(i)
313 . + normxy * dsigxy(i)
314 . + normyz * dsigyz(i)
315 . + normzx * dsigzx(i)
316c
317 ! 2 - Computation of DPHI_DLAMBDA
318 !---------------------------------------------------------
319c
320 ! a) Derivative with respect stress increments tensor DSIG
321 ! --------------------------------------------------------
322 dfdsig2 = normxx * normxx * g2(i)
323 . + normyy * normyy * g2(i)
324 . + normzz * normzz * g2(i)
325 . + normxy * normxy * g(i)
326 . + normyz * normyz * g(i)
327 . + normzx * normzx * g(i)
328c
329 ! b) Derivatives with respect to plastic strain P
330 ! ------------------------------------------------
331 hardp(i) = (young(i)*tang(i)/(young(i)-tang(i)))
332 IF (ivisc == 1) THEN
333 hardp(i) = hardp(i) + dyld0depsd(i)*dtinv +
334 . ((young(i)*dtangdepsd(i)*(young(i) - tang(i))
335 . + young(i)*tang(i)*dtangdepsd(i))/
336 . ((young(i) - tang(i))**2))*dtinv*pla(i)
337 ENDIF
338c
339 ! 3 - Computation of plastic multiplier and variables update
340 !----------------------------------------------------------
341c
342 ! Derivative of PHI with respect to DLAM
343 dphi_dlam(i) = - dfdsig2 - hardp(i)
344 dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20),dphi_dlam(i))
345c
346 ! Computation of the plastic multiplier
347 dlam = -(phi(i)+dphi)/dphi_dlam(i)
348c
349 ! Plastic strains tensor update
350 dpxx(i) = dlam * normxx
351 dpyy(i) = dlam * normyy
352 dpzz(i) = dlam * normzz
353 dpxy(i) = dlam * normxy
354 dpyz(i) = dlam * normyz
355 dpzx(i) = dlam * normzx
356c
357 ! Elasto-plastic stresses update
358 signxx(i) = signxx(i) - dpxx(i)*g2(i)
359 signyy(i) = signyy(i) - dpyy(i)*g2(i)
360 signzz(i) = signzz(i) - dpzz(i)*g2(i)
361 signxy(i) = signxy(i) - dpxy(i)*g(i)
362 signyz(i) = signyz(i) - dpyz(i)*g(i)
363 signzx(i) = signzx(i) - dpzx(i)*g(i)
364c
365 ! Computation of the trace of the new stress tensor
366 trsig(i) = signxx(i) + signyy(i) + signzz(i)
367 ! Computation of the new deviatoric stress tensor
368 sxx(i) = signxx(i) - trsig(i) * third
369 syy(i) = signyy(i) - trsig(i) * third
370 szz(i) = signzz(i) - trsig(i) * third
371 sxy(i) = signxy(i)
372 syz(i) = signyz(i)
373 szx(i) = signzx(i)
374c
375 ! Cumulated plastic strain and strain rate update
376 dpla(i) = max(zero,dpla(i) + dlam)
377 pla(i) = max(zero,pla(i) + dlam)
378 IF (ivisc == 1) THEN
379 epsd(i) = dpla(i)*dtinv
380 ENDIF
381c
382 ! Von Mises equivalent stress update
383 sigvm(i) = three_half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + three*sxy(i)**2
384 . + three*syz(i)**2 + three*szx(i)**2
385 sigvm(i) = sqrt(sigvm(i))
386c
387 IF (ivisc == 0) THEN
388 ! Yield stress update
389 yld(i) = yld0(i) + (young(i)*tang(i)/(young(i)-tang(i)))*pla(i)
390 ! Yield function value update
391 phi(i) = sigvm(i) - yld(i)
392 ENDIF
393c
394 ENDDO
395 ! End of the loop over yielding elements
396c
397 ! Update variable for full viscoplastic formulation
398 IF (ivisc == 1) THEN
399 ! Compute the initial yield stress
400 ipos(1:nel) = 1
401 iad(1:nel) = npf(ifunc(1)) / 2 + 1
402 ilen(1:nel) = npf(ifunc(1)+1) / 2 - iad(1:nel) - ipos(1:nel)
403 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_sig0,dyld0depsd,yld0)
404 yld0(1:nel) = yscale_sig0*yld0(1:nel)
405 dyld0depsd(1:nel) = yscale_sig0*dyld0depsd(1:nel)
406 ! Compute the Tangent modulus
407 IF (ifunc(3) > 0) THEN
408 ipos(1:nel) = 1
409 iad(1:nel) = npf(ifunc(3)) / 2 + 1
410 ilen(1:nel) = npf(ifunc(3)+1) / 2 - iad(1:nel) - ipos(1:nel)
411 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_tang,dtangdepsd,tang)
412 tang(1:nel) = yscale_tang*tang(1:nel)
413 dtangdepsd(1:nel) = yscale_tang*dtangdepsd(1:nel)
414 ENDIF
415 ! Updating values
416 DO ii=1,nindx
417 i = index(ii)
418 ! Check tangent modulus value
419 IF (tang(i) >= 0.99d0*young(i)) THEN
420 tang(i) = 0.99d0*young(i)
421 dtangdepsd(i) = zero
422 ENDIF
423 ! Yield stress update
424 yld(i) = yld0(i) + (young(i)*tang(i)/(young(i)-tang(i)))*pla(i)
425 ! Yield function value update
426 phi(i) = sigvm(i) - yld(i)
427 ENDDO
428 ENDIF
429 ENDIF
430 !===================================================================
431 ! - END OF PLASTIC CORRECTION WITH NICE - EXPLICIT 1 METHOD
432 !===================================================================
433c
434 ! Storing new values
435 DO i=1,nel
436 ! Storage of yield function value & coefficient for hourglass
437 IF (dpla(i) > zero) THEN
438 uvar(i,1) = phi(i)
439 et(i) = hardp(i) / (hardp(i) + young(i))
440 ELSE
441 uvar(i,1) = zero
442 et(i) = one
443 ENDIF
444 ! Storage of derivatives in case of viscoplastic formulation
445 IF (ivisc == 1) THEN
446 uvar(i,2) = dyld0depsd(i)
447 uvar(i,3) = dtangdepsd(i)
448 ENDIF
449 ! USR Outputs
450 seq(i) = sigvm(i) ! SIGEQ
451 ! Computation of the sound speed
452 soundsp(i) = sqrt((bulk(i) + four_over_3*g(i))/rho(i))
453 ! Storing the yield stress
454 sigy(i) = yld(i)
455 ENDDO
456c
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
vinter2
subroutine vinter2(tf, iad, ipos, ilen, nel0, x, dydx, y)
Definition vinter.F:143