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 . 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
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:144