OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps72.F File Reference
#include "implicit_f.inc"
#include "comlock.inc"
#include "mvsiz_p.inc"
#include "scr17_c.inc"
#include "com01_c.inc"
#include "com08_c.inc"
#include "units_c.inc"
#include "vectorize.inc"
#include "lockon.inc"
#include "lockoff.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps72 (nel, nuparam, nuvar, time, timestep, uparam, rho0, rho, depsxx, depsyy, depszz, depsxy, depsyz, depszx, sig0xx, sig0yy, sig0zz, sig0xy, sig0yz, sig0zx, signxx, signyy, signzz, signxy, signyz, signzx, soundsp, uvar, off, ngl, yld, pla, dpla, etse, seq, dmg, inloc, dplanl)

Function/Subroutine Documentation

◆ sigeps72()

subroutine sigeps72 ( integer nel,
integer nuparam,
integer nuvar,
time,
timestep,
uparam,
rho0,
rho,
depsxx,
depsyy,
depszz,
depsxy,
depsyz,
depszx,
sig0xx,
sig0yy,
sig0zz,
sig0xy,
sig0yz,
sig0zx,
signxx,
signyy,
signzz,
signxy,
signyz,
signzx,
soundsp,
uvar,
off,
integer, dimension(nel) ngl,
yld,
pla,
dpla,
etse,
seq,
dmg,
integer inloc,
dplanl )

Definition at line 29 of file sigeps72.F.

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