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 . dpdt,dlam,ddep,depxx,depyy,devepspxx,devepspyy,devepspzz,trepsp,
77 . normxx,normyy,normxy,denom,dfdsig2,dpdt_nl,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,test,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 ! Initial yield stress vs strain-rate curve
100 IF (ifunc(1) > 0) THEN
101 xscale_sig0 = uparam(16) ! Strain-rate scale factor
102 yscale_sig0 = uparam(17) ! Initial yield stress scale factor
103 yld0(1:nel) = zero
104 ELSE
105 yld0(1:nel) = uparam(17) ! Constant yield stress
106 dyld0depsd(1:nel) = zero
107 ENDIF
108 ! Young modulus vs strain-rate curve
109 xscale_youn = uparam(18) ! Strain-rate scale factor
110 yscale_youn = uparam(19) ! Young modulus scale factor
111 ! Tangent modulus vs strain-rate curve
112 IF (ifunc(3) > 0) THEN
113 xscale_tang = uparam(20) ! Strain-rate scale factor
114 yscale_tang = uparam(21) ! Tangent modulus scale factor
115 tang(1:nel) = zero
116 ELSE
117 tang(1:nel) = uparam(21) ! Constant tangent modulus
118 dtangdepsd(1:nel) = zero
119 ENDIF
120 dtinv = one/max(timestep,em20) ! Inverse of timestep
121c
122 ! Recovering internal variables
123 DO i=1,nel
124 IF (off(i) < em01) off(i) = zero
125 IF (off(i) < one) off(i) = off(i)*four_over_5
126 ! User inputs
127 phi0(i) = uvar(i,1) ! Previous yield function value
128 ! Standard inputs
129 dpla(i) = zero ! Initialization of the plastic strain increment
130 et(i) = one ! initialization of hourglass coefficient
131 hardp(i) = zero ! Initialization of hourglass coefficient
132 dezz(i) = zero ! Initialization of the strain increment in Z direction
133 dpxx(i) = zero ! Initialization of the XX plastic strain increment
134 dpyy(i) = zero ! Initialization of the YY plastic strain increment
135 dpxy(i) = zero ! Initialization of the ZZ plastic strain increment
136 dyld0depsd(i) = zero ! Initialization of the derivative of SIG0
137 dyoundepsd(i) = zero ! Initialization of the derivative of YOUN
138 dtangdepsd(i) = zero ! Initialization of the derivative of TANG
139 ENDDO
140c
141 ! Filling the strain-rate vector
142 IF (ivisc == 0) THEN
143 ! Compute effective strain-rate
144 DO i = 1,nel
145 trepsp = third*(epspxx(i) + epspyy(i))
146 devepspxx = epspxx(i) - trepsp
147 devepspyy = epspyy(i) - trepsp
148 devepspzz = -trepsp
149 depsdt = two_third*(devepspxx**2 + devepspyy**2 + devepspzz**2 +
150 . two*(epspxy(i)**2))
151 depsdt = sqrt(depsdt)
152 epsd(i) = afiltr * depsdt + (one - afiltr) * epsd(i)
153 ENDDO
154 ELSE
155 ! Reset plastic strain-rate
156 epsd(1:nel) = zero
157 ENDIF
158c
159 ! Compute the initial yield stress
160 IF (ifunc(1) > 0) THEN
161 ipos(1:nel) = 1
162 iad(1:nel) = npf(ifunc(1)) / 2 + 1
163 ilen(1:nel) = npf(ifunc(1)+1) / 2 - iad(1:nel) - ipos(1:nel)
164 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_sig0,dyld0depsd,yld0)
165 yld0(1:nel) = yscale_sig0*yld0(1:nel)
166 dyld0depsd(1:nel) = yscale_sig0*dyld0depsd(1:nel)
167 ENDIF
168 IF (ivisc == 1) THEN
169 IF (uvar(1,2) == zero) THEN
170 dyld0depsd(1:nel) = yscale_sig0*dyld0depsd(1:nel)
171 ELSE
172 dyld0depsd(1:nel) = uvar(1:nel,2)
173 ENDIF
174 ENDIF
175 ! Compute the Young modulus
176 IF (ifunc(2) > 0) THEN
177 ipos(1:nel) = 1
178 iad(1:nel) = npf(ifunc(2)) / 2 + 1
179 ilen(1:nel) = npf(ifunc(2)+1) / 2 - iad(1:nel) - ipos(1:nel)
180 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_youn,dyoundepsd,young)
181 young(1:nel) = yscale_youn*young(1:nel)
182 g(1:nel) = half * young(1:nel) / (one + nu)
183 bulk(1:nel) = third * young(1:nel) / (one - nu*two)
184 a11(1:nel) = young(1:nel) / (one - nu*nu)
185 a12(1:nel) = a11(1:nel) * nu
186 ENDIF
187 ! Compute the Tangent modulus
188 IF (ifunc(3) > 0) THEN
189 ipos(1:nel) = 1
190 iad(1:nel) = npf(ifunc(3)) / 2 + 1
191 ilen(1:nel) = npf(ifunc(3)+1) / 2 - iad(1:nel) - ipos(1:nel)
192 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_tang,dtangdepsd,tang)
193 tang(1:nel) = yscale_tang*tang(1:nel)
194 IF (ivisc == 1) THEN
195 IF (uvar(1,3) == zero) THEN
196 dtangdepsd(1:nel) = yscale_tang*dtangdepsd(1:nel)
197 ELSE
198 dtangdepsd(1:nel) = uvar(1:nel,3)
199 ENDIF
200 ENDIF
201 ENDIF
202 ! Check tangent modulus value + Assembling the yield stress
203 DO i = 1,nel
204 IF (tang(i) >= 0.99d0*young(i)) THEN
205 tang(i) = 0.99d0*young(i)
206 ENDIF
207 yld(i) = yld0(i) + (young(i)*tang(i)/(young(i)-tang(i)))*pla(i)
208 ENDDO
209c
210 !========================================================================
211 ! - COMPUTATION OF TRIAL VALUES
212 !========================================================================
213 DO i=1,nel
214 ! Computation of the trial stress tensor
215 signxx(i) = sigoxx(i) + a11(i)*depsxx(i) + a12(i)*depsyy(i)
216 signyy(i) = sigoyy(i) + a11(i)*depsyy(i) + a12(i)*depsxx(i)
217 signxy(i) = sigoxy(i) + depsxy(i)*g(i)
218 signyz(i) = sigoyz(i) + depsyz(i)*gs(i)
219 signzx(i) = sigozx(i) + depszx(i)*gs(i)
220 ! Computation of the trace of the trial stress tensor
221 trsig(i) = signxx(i) + signyy(i)
222 ! Computation of the deviatoric trial stress tensor
223 sxx(i) = signxx(i) - trsig(i) * third
224 syy(i) = signyy(i) - trsig(i) * third
225 szz(i) = -trsig(i) * third
226 sxy(i) = signxy(i)
227 ! Von Mises equivalent stress
228 sigvm(i) = three_half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + three*sxy(i)**2
229 sigvm(i) = sqrt(sigvm(i))
230 ENDDO
231c
232 !========================================================================
233 ! - COMPUTATION OF YIELD FONCTION
234 !========================================================================
235 phi(1:nel) = sigvm(1:nel) - yld(1:nel)
236c
237 ! Checking plastic behavior for all elements
238 nindx = 0
239 DO i=1,nel
240 IF (phi(i) > zero .AND. off(i) == one) THEN
241 nindx=nindx+1
242 index(nindx)=i
243 ENDIF
244 ENDDO
245c
246 !====================================================================
247 ! - PLASTIC CORRECTION WITH NICE - EXPLICIT METHOD
248 !====================================================================
249 IF (nindx > 0) THEN
250#include "vectorize.inc"
251 ! Loop over yielding elements
252 DO ii=1,nindx
253c
254 ! Number of the element with plastic behaviour
255 i = index(ii)
256c
257 ! Computation of the trial stress increment
258 dsigxx(i) = a11(i)*depsxx(i) + a12(i)*depsyy(i)
259 dsigyy(i) = a11(i)*depsyy(i) + a12(i)*depsxx(i)
260 dsigxy(i) = depsxy(i)*g(i)
261c
262 ! Computation of the old Von Mises stress
263 trsig(i) = sigoxx(i) + sigoyy(i)
264 sxx(i) = sigoxx(i) - trsig(i) * third
265 syy(i) = sigoyy(i) - trsig(i) * third
266 szz(i) = -trsig(i) * third
267 sxy(i) = sigoxy(i)
268 sigvm(i) = three_half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + three*sxy(i)**2
269 sigvm(i) = sqrt(sigvm(i))
270c
271 ! Note : in this part, the purpose is to compute for each iteration
272 ! a plastic multiplier allowing to update internal variables to satisfy
273 ! the consistency condition.
274 ! Its expression is : DLAMBDA = - (PHI + DPHI) / DPHI_DLAMBDA
275 ! -> PHI : current value of yield function (known)
276 ! -> DPHI : yield function prediction (to compute)
277 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
278 ! into account of internal variables kinetic :
279 ! plasticity, temperature and damage (to compute)
280c
281 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
282 !-------------------------------------------------------------
283 normxx = three_half*sxx(i)/sigvm(i)
284 normyy = three_half*syy(i)/sigvm(i)
285 normxy = three*sxy(i)/sigvm(i)
286c
287 ! Restoring previous value of the yield function
288 phi(i) = phi0(i)
289c
290 ! Computation of yield surface trial increment DPHI
291 dphi = normxx * dsigxx(i)
292 . + normyy * dsigyy(i)
293 . + normxy * dsigxy(i)
294c
295 ! 2 - Computation of DPHI_DLAMBDA
296 !---------------------------------------------------------
297c
298 ! a) Derivative with respect stress increments tensor DSIG
299 ! --------------------------------------------------------
300 dfdsig2 = normxx * (a11(i)*normxx + a12(i)*normyy)
301 . + normyy * (a11(i)*normyy + a12(i)*normxx)
302 . + normxy * normxy * g(i)
303c
304 ! b) Derivatives with respect to plastic strain P
305 ! ------------------------------------------------
306 hardp(i) = (young(i)*tang(i)/(young(i)-tang(i)))
307 IF (ivisc == 1) THEN
308 hardp(i) = hardp(i) + dyld0depsd(i)*dtinv +
309 . ((young(i)*dtangdepsd(i)*(young(i) - tang(i))
310 . + young(i)*tang(i)*dtangdepsd(i))/
311 . ((young(i) - tang(i))**2))*dtinv*pla(i)
312 ENDIF
313c
314 ! 3 - Computation of plastic multiplier and variables update
315 !----------------------------------------------------------
316c
317 ! Derivative of PHI with respect to DLAM
318 dphi_dlam(i) = - dfdsig2 - hardp(i)
319 dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20),dphi_dlam(i))
320c
321 ! Computation of the plastic multiplier
322 dlam = -(phi(i)+dphi)/dphi_dlam(i)
323c
324 ! Plastic strains tensor update
325 dpxx(i) = dlam * normxx
326 dpyy(i) = dlam * normyy
327 dpxy(i) = dlam * normxy
328c
329 ! Elasto-plastic stresses update
330 signxx(i) = signxx(i) - (a11(i)*dpxx(i) + a12(i)*dpyy(i))
331 signyy(i) = signyy(i) - (a11(i)*dpyy(i) + a12(i)*dpxx(i))
332 signxy(i) = signxy(i) - dpxy(i)*g(i)
333 trsig(i) = signxx(i) + signyy(i)
334 sxx(i) = signxx(i) - trsig(i) * third
335 syy(i) = signyy(i) - trsig(i) * third
336 szz(i) = - trsig(i) * third
337 sxy(i) = signxy(i)
338c
339 ! Cumulated plastic strain and strain rate update
340 dpla(i) = max(zero,dpla(i) + dlam)
341 pla(i) = max(zero,pla(i) + dlam)
342 IF (ivisc == 1) THEN
343 epsd(i) = dpla(i)*dtinv
344 ENDIF
345c
346 ! Von Mises equivalent stress update
347 sigvm(i) = three_half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + three*sxy(i)**2
348 sigvm(i) = sqrt(sigvm(i))
349c
350 IF (ivisc == 0) THEN
351 ! Yield stress update
352 yld(i) = yld0(i) + (young(i)*tang(i)/(young(i)-tang(i)))*pla(i)
353 ! Yield function value update
354 phi(i) = sigvm(i) - yld(i)
355 ENDIF
356c
357 ! Transverse strain update
358 IF (inloc == 0) THEN
359 dezz(i) = dezz(i) - dpxx(i) - dpyy(i)
360 ENDIF
361c
362 ENDDO
363 ! End of the loop over yielding elements
364c
365 ! Update variable for full viscoplastic formulation
366 IF (ivisc == 1) THEN
367 ! Compute the initial yield stress
368 ipos(1:nel) = 1
369 iad(1:nel) = npf(ifunc(1)) / 2 + 1
370 ilen(1:nel) = npf(ifunc(1)+1) / 2 - iad(1:nel) - ipos(1:nel)
371 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_sig0,dyld0depsd,yld0)
372 yld0(1:nel) = yscale_sig0*yld0(1:nel)
373 dyld0depsd(1:nel) = yscale_sig0*dyld0depsd(1:nel)
374 ! Compute the Tangent modulus
375 IF (ifunc(3) > 0) THEN
376 ipos(1:nel) = 1
377 iad(1:nel) = npf(ifunc(3)) / 2 + 1
378 ilen(1:nel) = npf(ifunc(3)+1) / 2 - iad(1:nel) - ipos(1:nel)
379 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_tang,dtangdepsd,tang)
380 tang(1:nel) = yscale_tang*tang(1:nel)
381 dtangdepsd(1:nel) = yscale_tang*dtangdepsd(1:nel)
382 ENDIF
383 ! Updating values
384 DO ii=1,nindx
385 i = index(ii)
386 ! Check tangent modulus value
387 IF (tang(i) >= 0.99d0*young(i)) THEN
388 tang(i) = 0.99d0*young(i)
389 dtangdepsd(i) = zero
390 ENDIF
391 ! Yield stress update
392 yld(i) = yld0(i) + (young(i)*tang(i)/(young(i)-tang(i)))*pla(i)
393 ! Yield function value update
394 phi(i) = sigvm(i) - yld(i)
395 ENDDO
396 ENDIF
397 ENDIF
398 !===================================================================
399 ! - END OF PLASTIC CORRECTION WITH NICE - EXPLICIT 1 METHOD
400 !===================================================================
401c
402 ! Storing new values
403 DO i=1,nel
404 ! Storage of yield function value & coefficient for hourglass
405 IF (dpla(i) > zero) THEN
406 uvar(i,1) = phi(i)
407 et(i) = hardp(i) / (hardp(i) + young(i))
408 ELSE
409 uvar(i,1) = zero
410 et(i) = one
411 ENDIF
412 ! Storage of derivatives in case of viscoplastic formulation
413 IF (ivisc == 1) THEN
414 uvar(i,2) = dyld0depsd(i)
415 uvar(i,3) = dtangdepsd(i)
416 ENDIF
417 ! USR Outputs
418 seq(i) = sigvm(i) ! SIGEQ
419 ! computation of the sound speed
420 soundsp(i) = sqrt(a11(i)/rho(i))
421 ! Storing the yield stress
422 sigy(i) = yld(i)
423 ! Thickness variation
424 IF (inloc > 0) THEN
425 IF (loff(i) == one) THEN
426 dezz(i) = -nu*(signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/young(i)
427 dezz(i) = dezz(i) - max(dplanl(i),zero)*half*(signxx(i)+signyy(i))/max(yld(i),em20)
428 ENDIF
429 ELSE
430 dezz(i) = -nu*(signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/young(i) + dezz(i)
431 ENDIF
432 ! Computation of the thickness variation
433 thk(i) = thk(i) + dezz(i)*thkly(i)*off(i)
434 ENDDO
435c
436 END
the
end diagonal values have been computed in the(sparse) matrix id.SOL
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
for
for(i8=*sizetab-1;i8 >=0;i8--)
Definition mumps_common.c:91
vinter2
subroutine vinter2(tf, iad, ipos, ilen, nel0, x, dydx, y)
Definition vinter.F:143