OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps72c.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "com01_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps72c (nel, nuparam, nuvar, time, timestep, uparam, rho0, thkly, depsxx, depsyy, depsxy, depsyz, depszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, soundsp, thk, pla, uvar, off, etse, gs, yld, hardm, seq, dpla, dmg, inloc, dplanl, loff)

Function/Subroutine Documentation

◆ sigeps72c()

subroutine sigeps72c ( integer nel,
integer nuparam,
integer nuvar,
time,
timestep,
uparam,
rho0,
thkly,
depsxx,
depsyy,
depsxy,
depsyz,
depszx,
sigoxx,
sigoyy,
sigoxy,
sigoyz,
sigozx,
signxx,
signyy,
signxy,
signyz,
signzx,
soundsp,
thk,
pla,
uvar,
off,
etse,
gs,
yld,
hardm,
seq,
dpla,
dmg,
integer inloc,
dplanl,
intent(in) loff )

Definition at line 28 of file sigeps72c.F.

37C-----------------------------------------------
38C I m p l i c i t T y p e s
39C-----------------------------------------------
40#include "implicit_f.inc"
41C-----------------------------------------------
42C G l o b a l P a r a m e t e r s
43C-----------------------------------------------
44#include "mvsiz_p.inc"
45C---------+---------+---+---+--------------------------------------------
46C C o m m o n B l o c k s
47C-----------------------------------------------
48#include "com01_c.inc"
49C-----------------------------------------------
50C I N P U T A r g u m e n t s
51C-----------------------------------------------
52C
53 INTEGER NEL,NUPARAM,NUVAR,INLOC
54 my_real
55 . time,timestep(nel),uparam(nuparam),
56 . rho0(nel),thkly(nel),pla(nel),
57 . depsxx(nel),depsyy(nel),
58 . depsxy(nel),depsyz(nel),depszx(nel),
59 . sigoxx(nel),sigoyy(nel),
60 . sigoxy(nel),sigoyz(nel),sigozx(nel),
61 . gs(*),hardm(nel),seq(nel),dmg(nel),
62 . dplanl(nel)
63 my_real, DIMENSION(NEL), INTENT(IN) :: loff
64C-----------------------------------------------
65C O U T P U T A r g u m e n t s
66C-----------------------------------------------
67 my_real
68 . signxx(nel),signyy(nel),
69 . signxy(nel),signyz(nel),signzx(nel),
70 . soundsp(nel),etse(nel),dpla(nel)
71C-----------------------------------------------
72C I N P U T O U T P U T A r g u m e n t s
73C-----------------------------------------------
74 my_real
75 . uvar(nel,nuvar), off(nel),thk(nel),yld(nel)
76C-----------------------------------------------
77C L o c a l V a r i a b l e s
78C-----------------------------------------------
79 INTEGER I,II,NINDX,INDEX(MVSIZ),NITER,ITER
80 my_real
81 . e,nu,g,a11,a12,
82 . sigy,eps0,nexp,
83 . ff,gg,hh,nn,
84 . c1,c2,c3,mexp,dc
85 my_real
86 . eta,cos3theta,theta,f1,f2,f3,cc(mvsiz),
87 . beta(mvsiz),dam(mvsiz),dezz(mvsiz),
88 . epsf,p,sighl(mvsiz),h(mvsiz),svm(mvsiz),
89 . normxx,normyy,normxy,dpxx(mvsiz),dpyy(mvsiz),
90 . dpzz(mvsiz),dpxy(mvsiz),deelzz(mvsiz),
91 . dpla_dlam(mvsiz),dlam,sig_dfdsig,dfdsig2,
92 . phi(mvsiz),dphi_dlam(mvsiz),ddep
93C-----------------------------------------------
94C USER VARIABLES INITIALIZATION
95C-----------------------------------------------
96C
97 ! Elastic parameters
98 e = uparam(1) ! Young modulus
99 nu = uparam(2) ! Poisson's ratio
100 g = uparam(3) ! Shear modulus
101 a11 = uparam(5) ! Diag. component of plane stress elastic matrix
102 a12 = uparam(6) ! Non diag component of plane stress elastic matrix
103 ! Hardening parameters
104 sigy = uparam(11) ! Initial yield stress
105 eps0 = uparam(12) ! Initial plastic strain (> 0)
106 nexp = uparam(13) ! Hardening exponent
107 ! Hill yield criterion parameters
108 ff = uparam(14) ! First Hill coefficient
109 gg = uparam(15) ! Second Hill coefficient
110 hh = uparam(16) ! Third Hill coefficient
111 nn = uparam(17) ! Fourth Hill coefficient
112 ! Modified Mohr-Coulomb yield criterion parameters
113 c1 = uparam(20) ! First failure parameter
114 c2 = uparam(21) ! Second failure parameter
115 c3 = uparam(22) ! Third failure parameter
116 mexp = uparam(23) ! Damage exponent
117 dc = uparam(24) ! Critical damage value
118c
119 ! Initial value
120 IF (isigi == 0) THEN
121 IF (time == zero) THEN
122 DO i=1,nel
123 uvar(i,1) = zero
124 ENDDO
125 ENDIF
126 ENDIF
127c
128 ! Recovering internal variables
129 DO i=1,nel
130 ! Checking deletion flag value
131 IF (off(i) < one) off(i) = four_over_5*off(i)
132 IF (off(i) < em01) off(i) = zero
133 ! Hourglass coefficient
134 etse(i) = one
135 h(i) = zero
136 ! Plastic strain increment
137 dpla(i) = zero
138 dpxx(i) = zero
139 dpyy(i) = zero
140 dpzz(i) = zero
141 dpxy(i) = zero
142 ! Damage variable
143 dam(i) = uvar(i,1)
144 ENDDO
145c
146 ! Computing yield stress and checking damage criteria
147 DO i = 1,nel
148 ! Computation of the BETA factor
149 beta(i) = one
150 IF (off(i) == one) THEN
151 IF ((dam(i) <= dc) .AND. (dam(i) >= one)) THEN
152 beta(i) = one/(max(dc - one,em20))
153 beta(i) = (dc - dam(i))*beta(i)
154 beta(i) = beta(i)**mexp
155 ENDIF
156 ENDIF
157 ! Computation of the yield stress
158 cc(i) = pla(i) + eps0
159 yld(i) = beta(i)*sigy*exp(nexp*log(cc(i)))
160 yld(i) = max(yld(i), em10)
161 ENDDO
162c
163 !========================================================================
164 ! - COMPUTATION OF TRIAL VALUES
165 !========================================================================
166 DO i = 1,nel
167c
168 ! Computation of the trial stress tensor
169 signxx(i) = sigoxx(i) + a11*depsxx(i)+a12*depsyy(i)
170 signyy(i) = sigoyy(i) + a12*depsxx(i)+a11*depsyy(i)
171 signxy(i) = sigoxy(i) + g*depsxy(i)
172 signyz(i) = sigoyz(i) + gs(i)*depsyz(i)
173 signzx(i) = sigozx(i) + gs(i)*depszx(i)
174C
175 ! Hill equivalent stress
176 sighl(i) = (ff+hh)*signyy(i)**2 + (gg+hh)*signxx(i)**2
177 . - two*hh*signxx(i)*signyy(i) + two*nn*signxy(i)**2
178 sighl(i) = sqrt(max(zero,sighl(i)))
179C
180 ENDDO
181c
182 !========================================================================
183 ! - COMPUTATION OF YIELD FONCTION
184 !========================================================================
185 phi(1:nel) = sighl(1:nel) - yld(1:nel)
186 ! Checking plastic behavior for all elements
187 nindx = 0
188 index = 0
189 DO i=1,nel
190 IF ((phi(i)>zero).AND.(off(i) == one)) THEN
191 nindx = nindx+1
192 index(nindx) = i
193 ENDIF
194 ENDDO
195c
196 !====================================================================
197 ! - PLASTIC CORRECTION WITH CUTTING PLANE ALGORITHM (NEWTON ITERATION)
198 !====================================================================
199c
200 ! Number of maximum Newton iterations
201 niter = 3
202c
203 ! Loop over the iterations
204 DO iter = 1, niter
205#include "vectorize.inc"
206 ! Loop over yielding elements
207 DO ii=1,nindx
208 i = index(ii)
209c
210 ! Note : in this part, the purpose is to compute for each iteration
211 ! a plastic multiplier allowing to update internal variables to satisfy
212 ! the consistency condition using the backward Euler implicit method
213 ! with a Newton-Raphson iterative procedure
214 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
215 ! -> PHI : current value of yield function (known)
216 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
217 ! into account of internal variables kinetic :
218 ! plasticity, temperature and damage (to compute)
219c
220 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
221 !-------------------------------------------------------------
222 normxx = (gg*signxx(i) + hh*(signxx(i)-signyy(i)))/(max(sighl(i),em20))
223 normyy = (ff*signyy(i) + hh*(signyy(i)-signxx(i)))/(max(sighl(i),em20))
224 normxy = two*nn*signxy(i)/(max(sighl(i),em20))
225c
226 ! 2 - Computation of DPHI_DLAMBDA
227 !---------------------------------------------------------
228c
229 ! a) Derivative with respect stress increments tensor DSIG
230 ! --------------------------------------------------------
231 dfdsig2 = normxx * (a11*normxx + a12*normyy)
232 . + normyy * (a11*normyy + a12*normxx)
233 . + normxy * normxy * g
234c
235 ! b) Derivatives with respect to plastic strain P
236 ! ------------------------------------------------
237c
238 ! i) Derivative of the yield stress with respect to plastic strain dYLD / dPLA
239 ! ----------------------------------------------------------------------------
240 h(i) = beta(i)*nexp*sigy*exp((nexp-1)*log(cc(i)))
241c
242 ! ii) Derivative of dPLA with respect to DLAM
243 ! -------------------------------------------
244 sig_dfdsig = signxx(i) * normxx
245 . + signyy(i) * normyy
246 . + signxy(i) * normxy
247 dpla_dlam(i) = sig_dfdsig / max(yld(i),em20)
248c
249 ! 3 - Computation of plastic multiplier and variables update
250 !----------------------------------------------------------
251c
252 ! Derivative of PHI with respect to DLAM
253 dphi_dlam(i) = - dfdsig2 - h(i)*dpla_dlam(i)
254 dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20) ,dphi_dlam(i))
255c
256 ! Computation of the plastic multiplier
257 dlam = -phi(i)/dphi_dlam(i)
258c
259 ! Plastic strains tensor update
260 dpxx(i) = dlam * normxx
261 dpyy(i) = dlam * normyy
262 dpxy(i) = dlam * normxy
263c
264 ! Elasto-plastic stresses update
265 signxx(i) = signxx(i) - (a11*dpxx(i) + a12*dpyy(i))
266 signyy(i) = signyy(i) - (a11*dpyy(i) + a12*dpxx(i))
267 signxy(i) = signxy(i) - dpxy(i)*g
268c
269 ! Cumulated plastic strain and strain rate update
270 ddep = dlam*dpla_dlam(i)
271 dpla(i) = max(zero, dpla(i) + ddep)
272 pla(i) = max(pla(i) + ddep,zero)
273c
274 ! Update the hardening yield stress
275 cc(i) = pla(i) + eps0
276 yld(i) = beta(i)*sigy*exp(nexp*log(cc(i)))
277 yld(i) = max(yld(i),em10)
278c
279 ! Update Hill equivalent stress
280 sighl(i) = (ff+hh)*signyy(i)**2 + (gg+hh)*signxx(i)**2
281 . - two*hh*signxx(i)*signyy(i) + two*nn*signxy(i)**2
282 sighl(i) = sqrt(max(zero,sighl(i)))
283c
284 ! Update yield function value
285 phi(i) = sighl(i) - yld(i)
286c
287 ! Transverse strain update
288 dpzz(i) = dpzz(i) - (dpxx(i)+dpyy(i))
289c
290 ENDDO
291 ! End of the loop over the yielding elements
292 ENDDO
293 ! End of the loop over the iterations
294 !===================================================================
295 ! - END OF PLASTIC CORRECTION WITH CUTTING PLANE ALGORITHM
296 !===================================================================
297C
298 ! Update and store new internal variables
299 DO i=1,nel
300 ! Modified Mohr Criteria Failure Model
301 IF (off(i) == one) THEN
302 ! Pressure
303 p = third*(signxx(i) + signyy(i))
304 ! Von Mises equivalent stress
305 svm(i) = signxx(i)**2 + signyy(i)**2 - signxx(i)*signyy(i) +
306 . three*signxy(i)**2
307 svm(i) = sqrt(svm(i))
308C
309 ! Computation of triaxiality
310 eta = p/max(svm(i),em20)
311 IF (eta < -two_third) eta = -two_third
312 IF (eta > two_third) eta = two_third
313 ! Computation of Lode angle
314 cos3theta = -half*twenty7*eta*(eta**2 - third)
315 IF (cos3theta < -one ) cos3theta = -one
316 IF (cos3theta > one ) cos3theta = one
317 theta = one - two*acos(cos3theta)/pi
318C
319 ! Computation of the failure factor
320 f1 = cos(theta*pi/six)
321 f2 = sin(theta*pi/six)
322 f3 = c3 + (sqrt(three)/(two - sqrt(three)))*(one - c3)*(one/max(f1,em20) - one)
323C
324 ! Computation of the failure plastic strain
325 epsf = (sigy/max(c2,em20))*f3*(f1*sqrt(third*(one + c1**2)) + c1*(eta + f2*third))
326 epsf = max(epsf,em20)**(-one/nexp)
327C
328 ! Computation of the damage variable
329 IF (inloc > 0) THEN
330 dam(i) = dam(i) + max(dplanl(i),zero)/max(em20,epsf)
331 ELSE
332 dam(i) = dam(i) + dpla(i)/max(em20,epsf)
333 ENDIF
334 IF (dam(i) >= dc) THEN
335 dam(i) = dc
336 off(i) = four_over_5
337 ENDIF
338C
339 ! Store the new value of damage
340 uvar(i,1) = dam(i)
341 ENDIF
342 ! Equivalent stress
343 seq(i) = sighl(i)
344 ! Normalized damage
345 dmg(i) = dam(i)/dc
346 ! Sound-speed
347 soundsp(i) = sqrt(a11/rho0(i))
348 ! Coefficient for hourglass
349 IF (dpla(i) > zero) THEN
350 etse(i) = h(i) / (h(i) + e)
351 hardm(i) = h(i)
352 ELSE
353 etse(i) = one
354 hardm(i) = zero
355 ENDIF
356 ! Computation of the thickness variation
357 deelzz(i) = -nu*(signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/e
358 IF (inloc > 0) THEN
359 IF (loff(i) == one) THEN
360 normxx = (gg*signxx(i) + hh*(signxx(i)-signyy(i)))/(max(sighl(i),em20))
361 normyy = (ff*signyy(i) + hh*(signyy(i)-signxx(i)))/(max(sighl(i),em20))
362 dezz(i) = deelzz(i) - max(dplanl(i),zero)*(normxx+normyy)
363 ENDIF
364 ELSE
365 dezz(i) = deelzz(i) + dpzz(i)
366 ENDIF
367 thk(i) = thk(i) + dezz(i)*thkly(i)*off(i)
368 ENDDO
369C
my_real
#define my_real
Definition cppsort.cpp:32
max
#define max(a, b)
Definition macros.h:21