OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_sahraei_s.F File Reference
#include "implicit_f.inc"
#include "units_c.inc"
#include "param_c.inc"
#include "scr17_c.inc"
#include "comlock.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine fail_sahraei_s (nel, nfunc, ifunc, npf, tf, time, ngl, uparam, epsxx, epsyy, epszz, epsxy, epsyz, epszx, off, dfmax, tdele, deltax, nuvar, uvar)

Function/Subroutine Documentation

◆ fail_sahraei_s()

subroutine fail_sahraei_s ( integer nel,
integer, intent(in) nfunc,
integer, dimension(nfunc), intent(in) ifunc,
integer, dimension(*) npf,
tf,
time,
integer, dimension(nel) ngl,
uparam,
epsxx,
epsyy,
epszz,
epsxy,
epsyz,
epszx,
off,
dfmax,
tdele,
deltax,
integer nuvar,
uvar )

Definition at line 33 of file fail_sahraei_s.F.

39C-----------------------------------------------
40C MIT Wierzbicki Sahraei electric battery failure
41C /FAIL/SAHRAEI1 - tabulated rupture criteria for solids
42C-----------------------------------------------
43C I m p l i c i t T y p e s
44C-----------------------------------------------
45#include "implicit_f.inc"
46C---------+---------+---+---+--------------------------------------------
47C VAR | SIZE |TYP| RW| DEFINITION
48C---------+---------+---+---+--------------------------------------------
49C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
50C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
51C NUVAR | 1 | I | R | NUMBER OF FAILURE ELEMENT VARIABLES
52C---------+---------+---+---+--------------------------------------------
53C NFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW not used
54C IFUNC | NFUNC | I | R | FUNCTION INDEX not used
55C NPF | * | I | R | FUNCTION ARRAY
56C TF | * | F | R | FUNCTION ARRAY
57C---------+---------+---+---+--------------------------------------------
58C TIME | 1 | F | R | CURRENT TIME
59C TIMESTEP| 1 | F | R | CURRENT TIME STEP
60C UPARAM | NUPARAM | F | R | USER FAILURE PARAMETER ARRAY
61C---------+---------+---+---+--------------------------------------------
62C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
63C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
64C ... | | | |
65C ... | | | |
66C---------+---------+---+---+--------------------------------------------
67C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
68C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
69C---------+---------+---+---+--------------------------------------------
70#include "units_c.inc"
71#include "param_c.inc"
72#include "scr17_c.inc"
73#include "comlock.inc"
74C-----------------------------------------------
75C I N P U T A r g u m e n t s
76C-----------------------------------------------
77C
78 INTEGER NEL, NUPARAM,NUVAR
79 INTEGER NGL(NEL)
80 INTEGER, INTENT(IN) :: NFUNC
81 INTEGER, DIMENSION(NFUNC), INTENT(IN) :: IFUNC
82 my_real time,uparam(*),
83 . epsxx(nel) ,epsyy(nel) ,epszz(nel) ,uvar(nel,nuvar),
84 . epsxy(nel) ,epsyz(nel) ,epszx(nel) ,deltax(nel)
85C-----------------------------------------------
86C I N P U T O U T P U T A r g u m e n t s
87C-----------------------------------------------
88 my_real off(nel), dfmax(nel),tdele(nel)
89C-----------------------------------------------
90C VARIABLES FOR FUNCTION INTERPOLATION
91C-----------------------------------------------
92 INTEGER NPF(*)
93 my_real finter ,tf(*),dydx
94 EXTERNAL finter
95C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
96C Y : y = f(x)
97C X : x
98C DYDX : f'(x) = dy/dx
99C IFUNC(J): FUNCTION INDEX
100C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
101C NPF,TF : FUNCTION PARAMETER
102C-----------------------------------------------
103C L o c a l V a r i a b l e s
104C-----------------------------------------------
105 INTEGER I,II,J,IFAIL,NINDX1,NINDX2,NUM,DEN,
106 . FAIL_ORT,ORDI,COMP_DIR,NINDXOFF,IDEL
107 my_real
108 . vol_strain_limit,max_comp_strain,ratio_2,el_ref,damage_sp(nel)
109 DOUBLE PRECISION
110 . I1,I2,I3,E11,E22,E33,EXX,EYY,EZZ,EXY,EZX,EYZ,
111 . Q,R,R_INTER,PHI,RATIO(NEL),DAMAGE(NEL),E00,FAC,
112 . VOL_STRAIN,NUMERATOR,DENOMINATOR,ORDINATE,
113 . LAMBDA,PLAMAX(NEL)
114 INTEGER,DIMENSION(NEL) :: INDX1,INDX2,INDX3,INDXOFF
115C-----------------------------------------------
116c!#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
117c!#
118c!# NUMERATOR = 1 ==> Eps_xx
119c!# NUMERATOR = 2 ==> Eps_yy
120c!# NUMERATOR = 3 ==> Eps_zz
121c!# NUMERATOR = 4 ==> Eps_1 1st. 3d-principal
122c!# NUMERATOR = 5 ==> Eps_2 2nd. 3d-principal
123c!# NUMERATOR = 6 ==> Eps_3 3rd. 3d-principal
124c!#
125c!# Denominator = 1 ==> 2D e1 = x-z - plane
126c!# Denominator = 2 ==> 2D e1 = x-y - plane
127c!# Denominator = 3 ==> 2D e1 = y-z - plane
128c!# Denominator = 4 ==> Eps_1 1st. 3d-principal
129c!# Denominator = 5 ==> Eps_2 2nd. 3d-principal
130c!# Denominator = 6 ==> Eps_3 3rd. 3d-principal
131c!#
132c!# Ordinate = 1 ==> MAX(Eps_xx,Eps_yy,Eps_zz)
133c!# Ordinate = 2 ==> Eps_xx
134c!# Ordinate = 3 ==> Eps_yy
135c!# Ordinate = 4 ==> Eps_zz
136c!# Ordinate = 5 ==> Eps_1 1st. 3d-principal
137c!# Ordinate = 6 ==> 2D - e1 = x-z - plane
138c!# Ordinate = 7 ==> 2D - e1 = x-y - plane
139c!# Ordinate = 8 ==> 2D - e1 = y-z - plane
140c!#
141c!# COMP_DIR = 1 ==> Pressure stress in xx
142c!# COMP_DIR = 2 ==> Pressure stress in yy
143c!# COMP_DIR = 3 ==> Pressure stress in zz
144c!#
145c!# MAX_COMP_STRAIN = In-plane failure compression strain (not deleting, only damage=1)
146c!# if value is negative, the element will NOT be deleted, only DAMAGE set to 1
147c!# RATIO = RATIO of the two other failure strains (default = 1.0)
148c!#
149c!
150C-----------------------------------------------
151 !=======================================================================
152 ! - INITIALISATION OF COMPUTATION ON TIME STEP
153 !=======================================================================
154 ! Recovering model parameters
155 vol_strain_limit = uparam(1)
156 num = int(uparam(2))
157 den = int(uparam(3))
158 ordi = int(uparam(4))
159 comp_dir = int(uparam(5))
160 max_comp_strain = uparam(6)
161 ratio_2 = uparam(7)
162 idel = int(uparam(8))
163 el_ref = uparam(9)
164c
165 ! Initialization
166 IF ((uvar(1,1) == zero).AND.(nfunc == 2)) THEN
167 DO i=1,nel
168 lambda = deltax(i)/el_ref
169 uvar(i,1) = finter(ifunc(2),lambda,npf,tf,dydx)
170 ENDDO
171 ENDIF
172c
173 ! Recovering internal variables
174 DO i=1,nel
175 ! Checking failure flag
176 IF (off(i) < one) off(i) = off(i)*four_over_5
177 IF (off(i) < em01) off(i) = zero
178 ! Resetting damage variable
179 damage(i) = zero
180 damage_sp(i) = zero
181 ! Recovering element size factor if needed
182 IF (nfunc == 2) THEN
183 fac = uvar(i,1)
184 ELSE
185 fac = one
186 ENDIF
187 END DO
188c
189 ! Initialization of counter
190 nindx1 = 0
191 nindx2 = 0
192 nindxoff = 0
193 DO i = 1,nel
194 IF (off(i) == one) THEN
195 nindxoff = nindxoff + 1
196 indxoff(nindxoff) = i
197 ENDIF
198 ENDDO
199c
200 !=======================================================================
201 ! - COMPUTATION OF FAILURE PLASTIC STRAIN
202 !=======================================================================
203 ! Loop over elements
204#include "vectorize.inc"
205 DO ii = 1,nindxoff
206c
207 ! Number of the element
208 i = indxoff(ii)
209C
210 ! Initialization of FAIL_ORT
211 fail_ort = ordi
212C
213 ! Copying strain value
214 exx = epsxx(i)
215 eyy = epsyy(i)
216 ezz = epszz(i)
217 exy = half*epsxy(i)
218 eyz = half*epsyz(i)
219 ezx = half*epszx(i)
220C
221 ! Computation of strain tensor invariants (attention au formalisme : EPSXY = GAMMAXY = 2*DEFXY)
222 i1 = exx + eyy + ezz
223 i2 = exx*eyy + eyy*ezz + ezz*exx - exy*exy - ezx*ezx - eyz*eyz
224 i3 = exx*eyy*ezz - exx*eyz*eyz - eyy*ezx*ezx - ezz*exy*exy + two*exy*ezx*eyz
225 q = (three*i2 - i1*i1)/nine
226 r = (two*i1*i1*i1-nine*i1*i2+twenty7*i3)/cinquante4
227 r_inter = min(r/sqrt(max(em20,(-q**3))),one)
228 phi = acos(max(r_inter,-one))
229 ! Computation of principal strain
230 e11 = two*sqrt(max(-q,zero))*cos( phi/three ) + third*i1
231 e22 = two*sqrt(max(-q,zero))*cos((phi+two*pi)/three ) + third*i1
232 e33 = two*sqrt(max(-q,zero))*cos((phi+four*pi)/three) + third*i1
233 ! Sorting principal strains
234 IF (e11 < e22) THEN
235 r_inter = e11
236 e11 = e22
237 e22 = r_inter
238 ENDIF
239 IF (e22 < e33)THEN
240 r_inter = e22
241 e22 = e33
242 e33 = r_inter
243 ENDIF
244 IF (e11 < e22)THEN
245 r_inter = e11
246 e11 = e22
247 e22 = r_inter
248 ENDIF
249c
250 ! Computation of volumetric strain
251 vol_strain = e11 + e22 + e33 + e11*e22 + e11*e33 + e22*e33 + e11*e22*e33
252c
253 !-----------------------------------------------------------------------------------
254 ! Volumetric strain failure
255 !-----------------------------------------------------------------------------------
256 IF (abs(vol_strain) >= vol_strain_limit) THEN
257C
258 ! Numerator flag for strain-ratio
259 IF (num == 1) THEN
260 numerator = exx
261 ELSEIF (num == 2) THEN
262 numerator = eyy
263 ELSEIF (num == 3) THEN
264 numerator = ezz
265 ELSEIF (num == 4) THEN
266 numerator = e11
267 ELSEIF (num == 5) THEN
268 numerator = e22
269 ELSE !IF (NUM == 6) THEN
270 numerator = e33
271 ENDIF
272C
273 ! Denominator flag for strain-ratio
274 IF (den == 1) THEN
275 denominator = ((exx + ezz)/two) + sqrt(((exx - ezz)/two)**2 + ezx**2)
276 ELSEIF (den == 2) THEN
277 denominator = ((exx + eyy)/two) + sqrt(((exx - eyy)/two)**2 + exy**2)
278 ELSEIF (den == 3) THEN
279 denominator = ((eyy + ezz)/two) + sqrt(((eyy - ezz)/two)**2 + eyz**2)
280 ELSEIF (den == 4) THEN
281 denominator = e11
282 ELSEIF (den == 5) THEN
283 denominator = e22
284 ELSE !IF (DEN == 6) THEN
285 denominator = e33
286 ENDIF
287C
288 ! Checking values
289 ratio(i) = abs(numerator/max(denominator,em20))
290 IF (numerator > zero .AND. denominator > zero) ratio(i) = zero
291 IF (numerator < zero .AND. denominator < zero) fail_ort = 1000
292 plamax(i) = finter(ifunc(1),ratio(i),npf,tf,dydx)
293 plamax(i) = plamax(i) * fac
294C
295 ! Case according to the ordinate
296 IF (fail_ort == 1) THEN
297 ! Computation of DAMAGE
298 damage(i)= max(exx/(max(plamax(i),em20)),
299 . eyy/(max(plamax(i),em20)),
300 . ezz/(max(plamax(i),em20)))
301 ! If the maximum strain is reached by one of the strain tensor diagonal components
302 IF (exx >= plamax(i) .OR. eyy >= plamax(i) .OR. ezz >= plamax(i)) THEN
303 tdele(i) = time
304 nindx1 = nindx1 + 1
305 indx1(nindx1) = i
306 damage(i) = one
307 off(i) = four_over_5
308 ENDIF
309 ELSEIF (fail_ort == 2) THEN
310 damage(i) = exx/(max(plamax(i),em20))
311 IF (exx >= plamax(i) .AND. exx > zero) THEN
312 tdele(i) = time
313 nindx1 = nindx1 + 1
314 indx1(nindx1) = i
315 damage(i) = one
316 off(i) = four_over_5
317 ENDIF
318 ELSEIF (fail_ort == 3) THEN
319 damage(i) = eyy/(max(plamax(i),em20))
320 IF (eyy >= plamax(i) .AND. eyy > zero) THEN
321 tdele(i) = time
322 nindx1 = nindx1 + 1
323 indx1(nindx1) = i
324 damage(i) = one
325 off(i) = four_over_5
326 ENDIF
327 ELSEIF (fail_ort == 4) THEN
328 damage(i) = ezz/(max(plamax(i),em20))
329 IF (ezz >= plamax(i) .AND. ezz > zero) THEN
330 tdele(i) = time
331 nindx1 = nindx1 + 1
332 indx1(nindx1) = i
333 damage(i) = one
334 off(i) = four_over_5
335 ENDIF
336 ELSEIF (fail_ort == 5) THEN
337 damage(i) = e11/(max(plamax(i),em20))
338 IF (e11 >= plamax(i) .AND. e11 > zero) THEN
339 tdele(i) = time
340 nindx1 = nindx1 + 1
341 indx1(nindx1) = i
342 damage(i) = one
343 off(i) = four_over_5
344 ENDIF
345 ELSEIF (fail_ort == 6) THEN
346 e00 = ((exx + ezz)/two) + sqrt(((exx - ezz)/two)**2 + ezx**2)
347 damage(i) = e00/(max(plamax(i),em20))
348 IF (e00 >= plamax(i) .AND. e00 > zero) THEN
349 tdele(i) = time
350 nindx1 = nindx1 + 1
351 indx1(nindx1) = i
352 damage(i) = one
353 off(i) = four_over_5
354 ENDIF
355 ELSEIF (fail_ort == 7) THEN
356 e00 = ((exx + eyy)/two) + sqrt(((exx - eyy)/two)**2 + exy**2)
357 damage(i) = e00/(max(plamax(i),em20))
358 IF (e00 >= plamax(i) .AND. e00 > zero) THEN
359 tdele(i) = time
360 nindx1 = nindx1 + 1
361 indx1(nindx1) = i
362 damage(i) = one
363 off(i) = four_over_5
364 ENDIF
365 ELSEIF (fail_ort == 8) THEN
366 e00 = ((eyy + ezz)/two) + sqrt(((eyy - ezz)/two)**2 + eyz**2)
367 damage(i) = e00/(max(plamax(i),em20))
368 IF (e00 >= plamax(i) .AND. e00 > zero) THEN
369 tdele(i) = time
370 nindx1 = nindx1 + 1
371 indx1(nindx1) = i
372 damage(i) = one
373 off(i) = four_over_5
374 ENDIF
375 ENDIF
376 ENDIF
377C
378 !-----------------------------------------------------------------------------------
379 ! In-plane compression failure
380 !-----------------------------------------------------------------------------------
381 IF (comp_dir > 0 .AND. dfmax(i) < one) THEN
382 IF (comp_dir == 1) THEN
383 IF (eyy < (-abs(max_comp_strain))) THEN
384 IF (idel > 0) THEN
385 off(i) = four_over_5
386 tdele(i) = time
387 ENDIF
388 damage(i) = one
389 nindx2 = nindx2 + 1
390 indx2(nindx2) = i
391 indx3(nindx2) = 2
392 ELSEIF (ezz < (-abs(max_comp_strain)*ratio_2)) THEN
393 IF (idel > 0) THEN
394 off(i) = four_over_5
395 tdele(i) = time
396 ENDIF
397 damage(i) = one
398 nindx2 = nindx2 + 1
399 indx2(nindx2) = i
400 indx3(nindx2) = 3
401 ENDIF
402 ELSEIF (comp_dir == 2) THEN
403 IF (ezz < (-abs(max_comp_strain))) THEN
404 IF (idel > 0) THEN
405 off(i) = four_over_5
406 tdele(i) = time
407 ENDIF
408 damage(i) = one
409 nindx2 = nindx2 + 1
410 indx2(nindx2) = i
411 indx3(nindx2) = 3
412 ELSEIF (exx < (-abs(max_comp_strain)*ratio_2)) THEN
413 IF (idel > 0) THEN
414 off(i) = four_over_5
415 tdele(i) = time
416 ENDIF
417 damage(i) = one
418 nindx2 = nindx2 + 1
419 indx2(nindx2) = i
420 indx3(nindx2) = 1
421 ENDIF
422 ELSEIF (comp_dir == 3) THEN
423 IF (exx < (-abs(max_comp_strain))) THEN
424 IF (idel > 0) THEN
425 off(i) = four_over_5
426 tdele(i) = time
427 ENDIF
428 damage(i) = one
429 nindx2 = nindx2 + 1
430 indx2(nindx2) = i
431 indx3(nindx2) = 1
432 ELSEIF (eyy < (-abs(max_comp_strain)*ratio_2)) THEN
433 IF (idel > 0) THEN
434 off(i) = four_over_5
435 tdele(i) = time
436 ENDIF
437 damage(i) = one
438 nindx2 = nindx2 + 1
439 indx2(nindx2) = i
440 indx3(nindx2) = 2
441 ENDIF
442 ENDIF
443 ENDIF
444c
445 ! Storing damage maximal value
446 damage_sp(i) = damage(i)
447 IF (dfmax(i) <= one) dfmax(i) = max(dfmax(i),damage_sp(i))
448 dfmax(i) = min(dfmax(i),one)
449c
450 ENDDO
451c
452 !=======================================================================
453 ! - PRINTOUT ELEMENT FAILURE
454 !=======================================================================
455 ! Volumetric strain failure
456 IF (nindx1 > 0)THEN
457 DO j=1,nindx1
458 i = indx1(j)
459#include "lockon.inc"
460 WRITE(iout, 1000) ngl(i),plamax(i),ratio(i),time
461 WRITE(istdo,1100) ngl(i),plamax(i),ratio(i),time
462#include "lockoff.inc"
463 ENDDO
464 ENDIF
465 ! Compression limit failure
466 IF (nindx2 > 0)THEN
467 DO j=1,nindx2
468 i = indx2(j)
469 IF (max_comp_strain > zero) THEN
470#include "lockon.inc"
471 WRITE(iout, 2000) ngl(i),time
472 WRITE(istdo,2100) ngl(i),time
473#include "lockoff.inc"
474 ELSE
475#include "lockon.inc"
476 WRITE(iout, 2002) ngl(i),time
477 WRITE(istdo,2102) ngl(i),time
478#include "lockoff.inc"
479 ENDIF
480 IF (indx3(i) == 1) THEN
481#include "lockon.inc"
482 WRITE(iout, 2010) epsxx(i)
483 WRITE(istdo,2110) epsxx(i)
484#include "lockoff.inc"
485 ENDIF
486 IF (indx3(i) == 2) THEN
487#include "lockon.inc"
488 WRITE(iout, 2020) epsyy(i)
489 WRITE(istdo,2120) epsyy(i)
490#include "lockoff.inc"
491 ENDIF
492 IF (indx3(i) == 3) THEN
493#include "lockon.inc"
494 WRITE(iout, 2030) epszz(i)
495 WRITE(istdo,2130) epszz(i)
496#include "lockoff.inc"
497 ENDIF
498 ENDDO
499 ENDIF
500C-----------------------------------------------
501 1000 FORMAT(1x,'DELETE SOLID ELEMENT NUMBER (SAHRAEI) EL :',i10,
502 . ' STRAIN LIMIT OF ',1pe10.3,'REACHED. AT STRAIN RATIO : ',1pe10.3,
503 . ' AT TIME :',1pe20.13)
504 1100 FORMAT(1x,'DELETE SOLID ELEMENT NUMBER (SAHRAEI) EL :',i10,
505 . ' STRAIN LIMIT OF ',1pe10.3,' REACHED. AT STRAIN RATIO : ',1pe10.3,
506 . ' AT TIME :',1pe20.13)
507 2000 FORMAT(1x,'SOLID ELEMENT NUMBER (SAHRAEI) EL :',i10,
508 . ' MAX. IN-PLANE-COMPRESSION STRAIN REACHED AT TIME :',1pe20.13)
509 2100 FORMAT(1x,'SOLID ELEMENT NUMBER (SAHRAEI) EL :',i10,
510 . ' MAX. IN-PLANE-COMPRESSION STRAIN REACHED AT TIME :',1pe20.13)
511 2002 FORMAT(1x,'SOLID ELEMENT NUMBER (SAHRAEI) EL :',i10,
512 . ' MAX. IN-PLANE-COMPRESSION STRAIN REACHED AT TIME :',1pe20.13)
513 2102 FORMAT(1x,'SOLID ELEMENT NUMBER (SAHRAEI) EL :',i10,
514 . ' MAX. IN-PLANE-COMPRESSION STRAIN REACHED AT TIME :',1pe20.13)
515c
516 2010 FORMAT(1x,'MAX PRESSURE STRAIN IN LOCAL X-DIRECTION REACHED : EPSXX= ',1pe20.13)
517 2110 FORMAT(1x,'MAX PRESSURE STRAIN IN LOCAL X-DIRECTION REACHED : EPSXX= ',1pe20.13)
518c
519 2020 FORMAT(1x,'MAX PRESSURE STRAIN IN LOCAL Y-DIRECTION REACHED : EPSYY= ',1pe20.13)
520 2120 FORMAT(1x,'MAX PRESSURE STRAIN IN LOCAL Y-DIRECTION REACHED : EPSYY= ',1pe20.13)
521c
522 2030 FORMAT(1x,'MAX PRESSURE STRAIN IN LOCAL Z-DIRECTION REACHED : EPSZZ= ',1pe20.13)
523 2130 FORMAT(1x,'MAX PRESSURE STRAIN IN LOCAL Z-DIRECTION REACHED : EPSZZ= ',1pe20.13)
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