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

Go to the source code of this file.

Functions/Subroutines

subroutine mat121c_newton (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_newton()

subroutine mat121c_newton ( 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_newton.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,ITER,NITER,NINDX,INDEX(NEL),IPOS(NEL),
71 . IAD(NEL),ILEN(NEL)
72 my_real
73 . young(nel),bulk(nel),g(nel),nu,a11(nel),a12(nel),nnu,tang(nel),
74 . afiltr,xscale_sig0,yscale_sig0,xscale_youn,yscale_youn,
75 . xscale_tang,yscale_tang
76 my_real
77 . dlam,devepspxx,devepspyy,devepspzz,trepsp,
78 . normxx,normyy,normxy,dfdsig2,depsdt,dtinv
79 my_real, DIMENSION(NEL) ::
80 . sxx,syy,szz,sxy,sigvm,yld,hardp,phi,dezz,yld0,dyld0depsd,
81 . dyoundepsd,dtangdepsd,trsig,dphi_dlam,dpxx,dpyy,dpxy
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 ! Standard inputs
127 dpla(i) = zero ! Initialization of the plastic strain increment
128 et(i) = one ! Initialization of hourglass coefficient
129 hardp(i) = zero ! Initialization of hourglass coefficient
130 dezz(i) = zero ! Initialization of the strain increment in Z direction
131 dpxx(i) = zero ! Initialization of the XX plastic strain increment
132 dpyy(i) = zero ! Initialization of the YY plastic strain increment
133 dpxy(i) = zero ! Initialization of the ZZ plastic strain increment
134 dyld0depsd(i) = zero ! Initialization of the derivative of SIG0
135 dyoundepsd(i) = zero ! Initialization of the derivative of YOUN
136 dtangdepsd(i) = zero ! Initialization of the derivative of TANG
137 ENDDO
138!
139 ! Filling the strain rate vector
140 IF (ivisc == 0) THEN
141 ! Compute effective strain-rate
142 DO i = 1,nel
143 trepsp = third*(epspxx(i) + epspyy(i))
144 devepspxx = epspxx(i) - trepsp
145 devepspyy = epspyy(i) - trepsp
146 devepspzz = -trepsp
147 depsdt = two_third*(devepspxx**2 + devepspyy**2 + devepspzz**2 +
148 . two*(epspxy(i)**2))
149 depsdt = sqrt(depsdt)
150 epsd(i) = afiltr * depsdt + (one - afiltr) * uvar(i,1)
151 uvar(i,1) = epsd(i)
152 ENDDO
153 ELSE
154 ! Reset plastic strain-rate
155 epsd(1:nel) = zero
156 ENDIF
157c
158 ! Compute the initial yield stress
159 IF (ifunc(1) > 0) THEN
160 ipos(1:nel) = 1
161 iad(1:nel) = npf(ifunc(1)) / 2 + 1
162 ilen(1:nel) = npf(ifunc(1)+1) / 2 - iad(1:nel) - ipos(1:nel)
163 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_sig0,dyld0depsd,yld0)
164 yld0(1:nel) = yscale_sig0*yld0(1:nel)
165 dyld0depsd(1:nel) = yscale_sig0*dyld0depsd(1:nel)
166 ENDIF
167 ! Compute the Young modulus
168 IF (ifunc(2) > 0) THEN
169 ipos(1:nel) = 1
170 iad(1:nel) = npf(ifunc(2)) / 2 + 1
171 ilen(1:nel) = npf(ifunc(2)+1) / 2 - iad(1:nel) - ipos(1:nel)
172 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_youn,dyoundepsd,young)
173 young(1:nel) = yscale_youn*young(1:nel)
174 g(1:nel) = half * young(1:nel) / (one + nu)
175 bulk(1:nel) = third * young(1:nel) / (one - nu*two)
176 a11(1:nel) = young(1:nel) / (one - nu*nu)
177 a12(1:nel) = a11(1:nel) * nu
178 ENDIF
179 ! Compute the Tangent modulus
180 IF (ifunc(3) > 0) THEN
181 ipos(1:nel) = 1
182 iad(1:nel) = npf(ifunc(3)) / 2 + 1
183 ilen(1:nel) = npf(ifunc(3)+1) / 2 - iad(1:nel) - ipos(1:nel)
184 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_tang,dtangdepsd,tang)
185 tang(1:nel) = yscale_tang*tang(1:nel)
186 dtangdepsd(1:nel) = yscale_tang*dtangdepsd(1:nel)
187 ENDIF
188 ! Check tangent modulus value + Assembling the yield stress
189 DO i = 1,nel
190 IF (tang(i) >= 0.99d0*young(i)) THEN
191 tang(i) = 0.99d0*young(i)
192 dtangdepsd(i) = zero
193 ENDIF
194 yld(i) = yld0(i) + (young(i)*tang(i)/(young(i)-tang(i)))*pla(i)
195 ENDDO
196c
197 !========================================================================
198 ! - COMPUTATION OF TRIAL VALUES
199 !========================================================================
200 DO i=1,nel
201 ! Computation of the trial stress tensor
202 signxx(i) = sigoxx(i) + a11(i)*depsxx(i) + a12(i)*depsyy(i)
203 signyy(i) = sigoyy(i) + a11(i)*depsyy(i) + a12(i)*depsxx(i)
204 signxy(i) = sigoxy(i) + depsxy(i)*g(i)
205 signyz(i) = sigoyz(i) + depsyz(i)*gs(i)
206 signzx(i) = sigozx(i) + depszx(i)*gs(i)
207 ! Computation of the trace of the trial stress tensor
208 trsig(i) = signxx(i) + signyy(i)
209 ! Computation of the deviatoric trial stress tensor
210 sxx(i) = signxx(i) - trsig(i) * third
211 syy(i) = signyy(i) - trsig(i) * third
212 szz(i) = -trsig(i) * third
213 sxy(i) = signxy(i)
214 ! Von Mises equivalent stress
215 sigvm(i) = three_half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + three*sxy(i)**2
216 sigvm(i) = sqrt(sigvm(i))
217 ENDDO
218c
219 !========================================================================
220 ! - COMPUTATION OF YIELD FONCTION
221 !========================================================================
222 phi(1:nel) = sigvm(1:nel) - yld(1:nel)
223c
224 ! Checking plastic behavior for all elements
225 nindx = 0
226 DO i=1,nel
227 IF (phi(i) > zero .AND. off(i) == one) THEN
228 nindx=nindx+1
229 index(nindx)=i
230 ENDIF
231 ENDDO
232c
233 !====================================================================
234 ! - PLASTIC CORRECTION WITH CUTTING PLANE (NEWTON-ITERATION) METHOD
235 !====================================================================
236 IF (nindx > 0) THEN
237c
238 ! Number of Newton iterations
239 IF (ivisc == 0) THEN
240 niter = 3
241 ELSE
242 niter = 5
243 ENDIF
244c
245 ! Loop over the iterations
246 DO iter = 1, niter
247#include "vectorize.inc"
248 ! Loop over yielding elements
249 DO ii=1,nindx
250c
251 ! Number of the element with plastic behaviour
252 i = index(ii)
253c
254 ! Note : in this part, the purpose is to compute for each iteration
255 ! a plastic multiplier allowing to update internal variables to satisfy
256 ! the consistency condition using the backward Euler implicit method
257 ! with a cutting plane iterative procedure
258 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
259 ! -> PHI : current value of yield function (known)
260 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
261 ! into account of internal variables kinetic :
262 ! plasticity, temperature and damage (to compute)
263c
264 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
265 !-------------------------------------------------------------
266 normxx = three_half*sxx(i)/sigvm(i)
267 normyy = three_half*syy(i)/sigvm(i)
268 normxy = three*sxy(i)/sigvm(i)
269c
270 ! 2 - Computation of DPHI_DLAMBDA
271 !---------------------------------------------------------
272c
273 ! a) Derivative with respect stress increments tensor DSIG
274 ! --------------------------------------------------------
275 dfdsig2 = normxx * (a11(i)*normxx + a12(i)*normyy)
276 . + normyy * (a11(i)*normyy + a12(i)*normxx)
277 . + normxy * normxy * g(i)
278c
279 ! b) Derivatives with respect to plastic strain P
280 ! ------------------------------------------------
281 hardp(i) = (young(i)*tang(i)/(young(i)-tang(i)))
282 IF (ivisc == 1) THEN
283 hardp(i) = hardp(i) + dyld0depsd(i)*dtinv +
284 . ((young(i)*dtangdepsd(i)*(young(i) - tang(i))
285 . + young(i)*tang(i)*dtangdepsd(i))/
286 . ((young(i) - tang(i))**2))*dtinv*pla(i)
287 ENDIF
288c
289 ! 3 - Computation of plastic multiplier and variables update
290 !----------------------------------------------------------
291c
292 ! Derivative of PHI with respect to DLAM
293 dphi_dlam(i) = - dfdsig2 - hardp(i)
294 dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20),dphi_dlam(i))
295c
296 ! Computation of the plastic multiplier
297 dlam = -phi(i)/dphi_dlam(i)
298c
299 ! Plastic strains tensor update
300 dpxx(i) = dlam * normxx
301 dpyy(i) = dlam * normyy
302 dpxy(i) = dlam * normxy
303c
304 ! Elasto-plastic stresses update
305 signxx(i) = signxx(i) - (a11(i)*dpxx(i) + a12(i)*dpyy(i))
306 signyy(i) = signyy(i) - (a11(i)*dpyy(i) + a12(i)*dpxx(i))
307 signxy(i) = signxy(i) - dpxy(i)*g(i)
308 trsig(i) = signxx(i) + signyy(i)
309 sxx(i) = signxx(i) - trsig(i) * third
310 syy(i) = signyy(i) - trsig(i) * third
311 szz(i) = - trsig(i) * third
312 sxy(i) = signxy(i)
313c
314 ! Cumulated plastic strain and strain rate update
315 dpla(i) = max(zero,dpla(i) + dlam)
316 pla(i) = max(zero,pla(i) + dlam)
317 IF (ivisc == 1) THEN
318 epsd(i) = dpla(i)*dtinv
319 ENDIF
320c
321 ! Von Mises equivalent stress update
322 sigvm(i) = three_half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + three*sxy(i)**2
323 sigvm(i) = sqrt(sigvm(i))
324c
325 IF (ivisc == 0) THEN
326 ! Yield stress update
327 yld(i) = yld0(i) + (young(i)*tang(i)/(young(i)-tang(i)))*pla(i)
328 ! Yield function value update
329 phi(i) = sigvm(i) - yld(i)
330 ENDIF
331c
332 ! Transverse strain update
333 IF (inloc == 0) THEN
334 dezz(i) = dezz(i) - dpxx(i) - dpyy(i)
335 ENDIF
336c
337 ENDDO
338 ! End of the loop over yielding elements
339c
340 ! Update variable for full viscoplastic formulation
341 IF (ivisc == 1) THEN
342 ! Compute the initial yield stress
343 ipos(1:nel) = 1
344 iad(1:nel) = npf(ifunc(1)) / 2 + 1
345 ilen(1:nel) = npf(ifunc(1)+1) / 2 - iad(1:nel) - ipos(1:nel)
346 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_sig0,dyld0depsd,yld0)
347 yld0(1:nel) = yscale_sig0*yld0(1:nel)
348 dyld0depsd(1:nel) = yscale_sig0*dyld0depsd(1:nel)
349 ! Compute the Tangent modulus
350 IF (ifunc(3) > 0) THEN
351 ipos(1:nel) = 1
352 iad(1:nel) = npf(ifunc(3)) / 2 + 1
353 ilen(1:nel) = npf(ifunc(3)+1) / 2 - iad(1:nel) - ipos(1:nel)
354 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_tang,dtangdepsd,tang)
355 tang(1:nel) = yscale_tang*tang(1:nel)
356 dtangdepsd(1:nel) = yscale_tang*dtangdepsd(1:nel)
357 ENDIF
358 ! Updating values
359 DO ii=1,nindx
360 i = index(ii)
361 ! Check tangent modulus value
362 IF (tang(i) >= 0.99d0*young(i)) THEN
363 tang(i) = 0.99d0*young(i)
364 dtangdepsd(i) = zero
365 ENDIF
366 ! Yield stress update
367 yld(i) = yld0(i) + (young(i)*tang(i)/(young(i)-tang(i)))*pla(i)
368 ! Yield function value update
369 phi(i) = sigvm(i) - yld(i)
370 ENDDO
371 ENDIF
372 ENDDO
373 ! End of the loop over the iterations
374 ENDIF
375 !=========================================================================
376 ! - END OF PLASTIC CORRECTION WITH CUTTING PLANE (NEWTON-ITERATION) METHOD
377 !=========================================================================
378c
379 ! Storing new values
380 DO i=1,nel
381 ! USR Outputs
382 seq(i) = sigvm(i) ! SIGEQ
383 ! Coefficient for hourglass
384 IF (dpla(i) > zero) THEN
385 et(i) = hardp(i) / (hardp(i) + young(i))
386 ELSE
387 et(i) = one
388 ENDIF
389 ! Computation of the sound speed
390 soundsp(i) = sqrt(a11(i)/rho(i))
391 ! Storing the yield stress
392 sigy(i) = yld(i)
393 ! Thickness variation
394 IF (inloc > 0) THEN
395 IF (loff(i) == one) THEN
396 dezz(i) = -nu*(signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/young(i)
397 dezz(i) = dezz(i) - max(dplanl(i),zero)*half*(signxx(i)+signyy(i))/max(yld(i),em20)
398 ENDIF
399 ELSE
400 dezz(i) = -nu*(signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/young(i) + dezz(i)
401 ENDIF
402 ! Computation of the thickness variation
403 thk(i) = thk(i) + dezz(i)*thkly(i)*off(i)
404 ENDDO
405c
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