OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat121_nice.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!|| mat121_nice ../engine/source/materials/mat/mat121/mat121_nice.F
25!||--- called by ------------------------------------------------------
26!|| sigeps121 ../engine/source/materials/mat/mat121/sigeps121.F
27!||--- calls -----------------------------------------------------
28!|| vinter2 ../engine/source/tools/curve/vinter.F
29!||====================================================================
30 SUBROUTINE mat121_nice(
31 1 NEL ,NGL ,NUPARAM ,NUVAR ,NFUNC ,IFUNC ,NPF ,
32 2 TF ,TIMESTEP,TIME ,UPARAM ,UVAR ,RHO ,PLA ,
33 3 DPLA ,SOUNDSP ,EPSD ,OFF ,
34 4 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
35 5 EPSPXX ,EPSPYY ,EPSPZZ ,EPSPXY ,EPSPYZ ,EPSPZX ,
36 6 SIGOXX ,SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
37 7 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
38 8 SIGY ,ET ,SEQ )
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
457 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
mat121_nice
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)
Definition mat121_nice.F:39
vinter2
subroutine vinter2(tf, iad, ipos, ilen, nel0, x, dydx, y)
Definition vinter.F:143