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

Go to the source code of this file.

Functions/Subroutines

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)

Function/Subroutine Documentation

◆ mat121c_nice()

subroutine mat121c_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(in) gs,
intent(inout) thk,
intent(in) thkly,
intent(inout) off,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
intent(in) epspxx,
intent(in) epspyy,
intent(in) epspxy,
intent(in) epspyz,
intent(in) epspzx,
intent(in) sigoxx,
intent(in) sigoyy,
intent(in) sigoxy,
intent(in) sigoyz,
intent(in) sigozx,
intent(out) signxx,
intent(out) signyy,
intent(out) signxy,
intent(out) signyz,
intent(out) signzx,
intent(out) sigy,
intent(out) et,
intent(in) dplanl,
intent(inout) seq,
integer inloc,
intent(in) loff )

Definition at line 30 of file mat121c_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,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
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