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