OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat121c_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!|| mat121c_nice ../engine/source/materials/mat/mat121/mat121c_nice.F
25!||--- called by ------------------------------------------------------
26!|| sigeps121c ../engine/source/materials/mat/mat121/sigeps121c.F
27!||--- calls -----------------------------------------------------
28!|| vinter2 ../engine/source/tools/curve/vinter.F
29!||====================================================================
30 SUBROUTINE mat121c_nice(
31 1 NEL ,NGL ,NUPARAM ,NUVAR ,NFUNC ,IFUNC ,NPF ,
32 2 TF ,TIMESTEP,TIME ,UPARAM ,UVAR ,RHO ,PLA ,
33 3 DPLA ,SOUNDSP ,EPSD ,GS ,THK ,THKLY ,OFF ,
34 4 DEPSXX ,DEPSYY ,DEPSXY ,DEPSYZ ,DEPSZX ,
35 5 EPSPXX ,EPSPYY ,EPSPXY ,EPSPYZ ,EPSPZX ,
36 6 SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
37 7 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
38 8 SIGY ,ET ,DPLANL ,SEQ ,INLOC ,LOFF )
39 !=======================================================================
40 ! Implicit types
41 !=======================================================================
42#include "implicit_f.inc"
43 !=======================================================================
44 ! Common
45 !=======================================================================
46 !=======================================================================
47 ! Dummy arguments
48 !=======================================================================
49 INTEGER NEL,NUPARAM,NUVAR,INLOC,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,DPLANL,GS,THKLY,LOFF,
57 . depsxx,depsyy,depsxy,depsyz,depszx,
58 . epspxx,epspyy,epspxy,epspyz,epspzx,
59 . sigoxx,sigoyy,sigoxy,sigoyz,sigozx
60 my_real ,DIMENSION(NEL), INTENT(OUT) ::
61 . soundsp,sigy,et,
62 . signxx,signyy,signxy,signyz,signzx
63 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
64 . pla,dpla,epsd,off,thk,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,a11(nel),a12(nel),nnu,tang(nel),
73 . afiltr,xscale_sig0,yscale_sig0,xscale_youn,yscale_youn,
74 . xscale_tang,yscale_tang
75 my_real
76 . dlam,devepspxx,devepspyy,devepspzz,trepsp,
77 . normxx,normyy,normxy,dfdsig2,depsdt,dphi,dtinv
78 my_real, DIMENSION(NEL) ::
79 . sxx,syy,szz,sxy,sigvm,yld,hardp,phi,dezz,yld0,dyld0depsd,phi0,
80 . dyoundepsd,dtangdepsd,trsig,dphi_dlam,dpxx,dpyy,dpxy,
81 . dsigxx,dsigyy,dsigxy
82c
83 !=======================================================================
84 ! - INITIALISATION OF COMPUTATION ON TIME STEP
85 !=======================================================================
86 ! Recovering model parameters
87 ! Elastic parameters
88 young(1:nel) = uparam(1) ! Young modulus
89 bulk(1:nel) = uparam(2) ! Bulk modulus
90 g(1:nel) = uparam(3) ! Shear modulus
91 nu = uparam(6) ! Poisson ration
92 nnu = uparam(7) ! NU/(1-NU)
93 a11(1:nel) = uparam(9) ! Diagonal term, elastic matrix in plane stress
94 a12(1:nel) = uparam(10) ! Non-diagonal term, elastic matrix in plane stress
95 ! Flags for computation
96 ivisc = nint(uparam(12)) ! Viscosity formulation
97 ! Strain-rate filtering (if Ivisc = 0)
98 afiltr = min(one, uparam(14)*timestep)
99!
100 if (ivisc == 0) epsd(1:nel) = uvar(1:nel,2)
101!
102 ! Initial yield stress vs strain-rate curve
103 IF (ifunc(1) > 0) THEN
104 xscale_sig0 = uparam(16) ! Strain-rate scale factor
105 yscale_sig0 = uparam(17) ! Initial yield stress scale factor
106 yld0(1:nel) = zero
107 ELSE
108 yld0(1:nel) = uparam(17) ! Constant yield stress
109 dyld0depsd(1:nel) = zero
110 ENDIF
111 ! Young modulus vs strain-rate curve
112 xscale_youn = uparam(18) ! Strain-rate scale factor
113 yscale_youn = uparam(19) ! Young modulus scale factor
114 ! Tangent modulus vs strain-rate curve
115 IF (ifunc(3) > 0) THEN
116 xscale_tang = uparam(20) ! Strain-rate scale factor
117 yscale_tang = uparam(21) ! Tangent modulus scale factor
118 tang(1:nel) = zero
119 ELSE
120 tang(1:nel) = uparam(21) ! Constant Tangent Modulus
121 dtangdepsd(1:nel) = zero
122 ENDIF
123 dtinv = one/max(timestep,em20) ! Inverse of timestep
124c
125 ! Recovering internal variables
126 DO i=1,nel
127 IF (off(i) < em01) off(i) = zero
128 IF (off(i) < one) off(i) = off(i)*four_over_5
129 ! User inputs
130 phi0(i) = uvar(i,1) ! Previous yield function value
131 ! Standard inputs
132 dpla(i) = zero ! Initialization of the plastic strain increment
133 et(i) = one ! Initialization of hourglass coefficient
134 hardp(i) = zero ! Initialization of hourglass coefficient
135 dezz(i) = zero ! Initialization of the strain increment in Z direction
136 dpxx(i) = zero ! Initialization of the XX plastic strain increment
137 dpyy(i) = zero ! Initialization of the YY plastic strain increment
138 dpxy(i) = zero ! Initialization of the ZZ 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))
149 devepspxx = epspxx(i) - trepsp
150 devepspyy = epspyy(i) - trepsp
151 devepspzz = -trepsp
152 depsdt = two_third*(devepspxx**2 + devepspyy**2 + devepspzz**2 +
153 . two*(epspxy(i)**2))
154 depsdt = sqrt(depsdt)
155 epsd(i) = afiltr * depsdt + (one - afiltr) * epsd(i)
156 uvar(i,2) = 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 bulk(1:nel) = third * young(1:nel) / (one - nu*two)
188 a11(1:nel) = young(1:nel) / (one - nu*nu)
189 a12(1:nel) = a11(1:nel) * 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 signxx(i) = sigoxx(i) + a11(i)*depsxx(i) + a12(i)*depsyy(i)
220 signyy(i) = sigoyy(i) + a11(i)*depsyy(i) + a12(i)*depsxx(i)
221 signxy(i) = sigoxy(i) + depsxy(i)*g(i)
222 signyz(i) = sigoyz(i) + depsyz(i)*gs(i)
223 signzx(i) = sigozx(i) + depszx(i)*gs(i)
224 ! Computation of the trace of the trial stress tensor
225 trsig(i) = signxx(i) + signyy(i)
226 ! Computation of the deviatoric trial stress tensor
227 sxx(i) = signxx(i) - trsig(i) * third
228 syy(i) = signyy(i) - trsig(i) * third
229 szz(i) = -trsig(i) * third
230 sxy(i) = signxy(i)
231 ! Von Mises equivalent stress
232 sigvm(i) = three_half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + three*sxy(i)**2
233 sigvm(i) = sqrt(sigvm(i))
234 ENDDO
235c
236 !========================================================================
237 ! - COMPUTATION OF YIELD FONCTION
238 !========================================================================
239 phi(1:nel) = sigvm(1:nel) - yld(1:nel)
240c
241 ! Checking plastic behavior for all elements
242 nindx = 0
243 DO i=1,nel
244 IF (phi(i) > zero .AND. off(i) == one) THEN
245 nindx=nindx+1
246 index(nindx)=i
247 ENDIF
248 ENDDO
249c
250 !====================================================================
251 ! - PLASTIC CORRECTION WITH NICE - EXPLICIT METHOD
252 !====================================================================
253 IF (nindx > 0) THEN
254#include "vectorize.inc"
255 ! Loop over yielding elements
256 DO ii=1,nindx
257c
258 ! Number of the element with plastic behaviour
259 i = index(ii)
260c
261 ! Computation of the trial stress increment
262 dsigxx(i) = a11(i)*depsxx(i) + a12(i)*depsyy(i)
263 dsigyy(i) = a11(i)*depsyy(i) + a12(i)*depsxx(i)
264 dsigxy(i) = depsxy(i)*g(i)
265c
266 ! Computation of the old Von Mises stress
267 trsig(i) = sigoxx(i) + sigoyy(i)
268 sxx(i) = sigoxx(i) - trsig(i) * third
269 syy(i) = sigoyy(i) - trsig(i) * third
270 szz(i) = -trsig(i) * third
271 sxy(i) = sigoxy(i)
272 sigvm(i) = three_half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + three*sxy(i)**2
273 sigvm(i) = sqrt(sigvm(i))
274c
275 ! Note : in this part, the purpose is to compute for each iteration
276 ! a plastic multiplier allowing to update internal variables to satisfy
277 ! the consistency condition.
278 ! Its expression is : DLAMBDA = - (PHI + DPHI) / DPHI_DLAMBDA
279 ! -> PHI : current value of yield function (known)
280 ! -> DPHI : yield function prediction (to compute)
281 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
282 ! into account of internal variables kinetic :
283 ! plasticity, temperature and damage (to compute)
284c
285 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
286 !-------------------------------------------------------------
287 normxx = three_half*sxx(i)/sigvm(i)
288 normyy = three_half*syy(i)/sigvm(i)
289 normxy = three*sxy(i)/sigvm(i)
290c
291 ! Restoring previous value of the yield function
292 phi(i) = phi0(i)
293c
294 ! Computation of yield surface trial increment DPHI
295 dphi = normxx * dsigxx(i)
296 . + normyy * dsigyy(i)
297 . + normxy * dsigxy(i)
298c
299 ! 2 - Computation of DPHI_DLAMBDA
300 !---------------------------------------------------------
301c
302 ! a) Derivative with respect stress increments tensor DSIG
303 ! --------------------------------------------------------
304 dfdsig2 = normxx * (a11(i)*normxx + a12(i)*normyy)
305 . + normyy * (a11(i)*normyy + a12(i)*normxx)
306 . + normxy * normxy * g(i)
307c
308 ! b) Derivatives with respect to plastic strain P
309 ! ------------------------------------------------
310 hardp(i) = (young(i)*tang(i)/(young(i)-tang(i)))
311 IF (ivisc == 1) THEN
312 hardp(i) = hardp(i) + dyld0depsd(i)*dtinv +
313 . ((young(i)*dtangdepsd(i)*(young(i) - tang(i))
314 . + young(i)*tang(i)*dtangdepsd(i))/
315 . ((young(i) - tang(i))**2))*dtinv*pla(i)
316 ENDIF
317c
318 ! 3 - Computation of plastic multiplier and variables update
319 !----------------------------------------------------------
320c
321 ! Derivative of PHI with respect to DLAM
322 dphi_dlam(i) = - dfdsig2 - hardp(i)
323 dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20),dphi_dlam(i))
324c
325 ! Computation of the plastic multiplier
326 dlam = -(phi(i)+dphi)/dphi_dlam(i)
327c
328 ! Plastic strains tensor update
329 dpxx(i) = dlam * normxx
330 dpyy(i) = dlam * normyy
331 dpxy(i) = dlam * normxy
332c
333 ! Elasto-plastic stresses update
334 signxx(i) = signxx(i) - (a11(i)*dpxx(i) + a12(i)*dpyy(i))
335 signyy(i) = signyy(i) - (a11(i)*dpyy(i) + a12(i)*dpxx(i))
336 signxy(i) = signxy(i) - dpxy(i)*g(i)
337 trsig(i) = signxx(i) + signyy(i)
338 sxx(i) = signxx(i) - trsig(i) * third
339 syy(i) = signyy(i) - trsig(i) * third
340 szz(i) = - trsig(i) * third
341 sxy(i) = signxy(i)
342c
343 ! Cumulated plastic strain and strain rate update
344 dpla(i) = max(zero,dpla(i) + dlam)
345 pla(i) = max(zero,pla(i) + dlam)
346 IF (ivisc == 1) THEN
347 epsd(i) = dpla(i)*dtinv
348 ENDIF
349c
350 ! Von Mises equivalent stress update
351 sigvm(i) = three_half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + three*sxy(i)**2
352 sigvm(i) = sqrt(sigvm(i))
353c
354 IF (ivisc == 0) THEN
355 ! Yield stress update
356 yld(i) = yld0(i) + (young(i)*tang(i)/(young(i)-tang(i)))*pla(i)
357 ! Yield function value update
358 phi(i) = sigvm(i) - yld(i)
359 ENDIF
360c
361 ! Transverse strain update
362 IF (inloc == 0) THEN
363 dezz(i) = dezz(i) - dpxx(i) - dpyy(i)
364 ENDIF
365c
366 ENDDO
367 ! End of the loop over yielding elements
368c
369 ! Update variable for full viscoplastic formulation
370 IF (ivisc == 1) THEN
371 ! Compute the initial yield stress
372 ipos(1:nel) = 1
373 iad(1:nel) = npf(ifunc(1)) / 2 + 1
374 ilen(1:nel) = npf(ifunc(1)+1) / 2 - iad(1:nel) - ipos(1:nel)
375 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_sig0,dyld0depsd,yld0)
376 yld0(1:nel) = yscale_sig0*yld0(1:nel)
377 dyld0depsd(1:nel) = yscale_sig0*dyld0depsd(1:nel)
378 ! Compute the Tangent modulus
379 IF (ifunc(3) > 0) THEN
380 ipos(1:nel) = 1
381 iad(1:nel) = npf(ifunc(3)) / 2 + 1
382 ilen(1:nel) = npf(ifunc(3)+1) / 2 - iad(1:nel) - ipos(1:nel)
383 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_tang,dtangdepsd,tang)
384 tang(1:nel) = yscale_tang*tang(1:nel)
385 dtangdepsd(1:nel) = yscale_tang*dtangdepsd(1:nel)
386 ENDIF
387 ! Updating values
388 DO ii=1,nindx
389 i = index(ii)
390 ! Check tangent modulus value
391 IF (tang(i) >= 0.99d0*young(i)) THEN
392 tang(i) = 0.99d0*young(i)
393 dtangdepsd(i) = zero
394 ENDIF
395 ! yield stress update
396 yld(i) = yld0(i) + (young(i)*tang(i)/(young(i)-tang(i)))*pla(i)
397 ! Yield function value update
398 phi(i) = sigvm(i) - yld(i)
399 ENDDO
400 ENDIF
401 ENDIF
402 !===================================================================
403 ! - END OF PLASTIC CORRECTION WITH NICE - EXPLICIT 1 METHOD
404 !===================================================================
405c
406 ! Storing new values
407 DO i=1,nel
408 ! Storage of yield function value & coefficient for hourglass
409 IF (dpla(i) > zero) THEN
410 uvar(i,1) = phi(i)
411 et(i) = hardp(i) / (hardp(i) + young(i))
412 ELSE
413 uvar(i,1) = zero
414 et(i) = one
415 ENDIF
416 ! Storage of derivatives in case of viscoplastic formulation
417 IF (ivisc == 1) THEN
418 uvar(i,2) = dyld0depsd(i)
419 uvar(i,3) = dtangdepsd(i)
420 ENDIF
421 ! USR Outputs
422 seq(i) = sigvm(i) ! SIGEQ
423 ! Computation of the sound speed
424 soundsp(i) = sqrt(a11(i)/rho(i))
425 ! Storing the yield stress
426 sigy(i) = yld(i)
427 ! Thickness variation
428 IF (inloc > 0) THEN
429 IF (loff(i) == one) THEN
430 dezz(i) = -nu*(signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/young(i)
431 dezz(i) = dezz(i) - max(dplanl(i),zero)*half*(signxx(i)+signyy(i))/max(yld(i),em20)
432 ENDIF
433 ELSE
434 dezz(i) = -nu*(signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/young(i) + dezz(i)
435 ENDIF
436 ! Computation of the thickness variation
437 thk(i) = thk(i) + dezz(i)*thkly(i)*off(i)
438 ENDDO
439c
440 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
mat121c_nice
subroutine mat121c_nice(nel, ngl, nuparam, nuvar, nfunc, ifunc, npf, tf, timestep, time, uparam, uvar, rho, pla, dpla, soundsp, epsd, gs, thk, thkly, off, depsxx, depsyy, depsxy, depsyz, depszx, epspxx, epspyy, epspxy, epspyz, epspzx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, sigy, et, dplanl, seq, inloc, loff)
Definition mat121c_nice.F:39
vinter2
subroutine vinter2(tf, iad, ipos, ilen, nel0, x, dydx, y)
Definition vinter.F:144