OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_tbutcher_xfem.F File Reference
#include "implicit_f.inc"
#include "units_c.inc"
#include "comlock.inc"
#include "com_xfem1.inc"
#include "lockon.inc"
#include "lockoff.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine fail_tbutcher_xfem (nel, nuparam, uparam, nuvar, uvar, time, timestep, tens, dfmax, tdel, signxx, signyy, signxy, signyz, signzx, ngl, ipt, nptot, noff, off, offl, elcrkini, ixfem, ixel, ilay, iptt, nptt, npttf)

Function/Subroutine Documentation

◆ fail_tbutcher_xfem()

subroutine fail_tbutcher_xfem ( integer nel,
integer nuparam,
uparam,
integer nuvar,
uvar,
time,
timestep,
tens,
dfmax,
tdel,
signxx,
signyy,
signxy,
signyz,
signzx,
integer, dimension(nel) ngl,
integer ipt,
integer nptot,
integer, dimension(nel) noff,
off,
offl,
integer, dimension(nxlaymax,*) elcrkini,
integer ixfem,
integer ixel,
integer ilay,
integer iptt,
integer nptt,
npttf )

Definition at line 29 of file fail_tbutcher_xfem.F.

36C---------+---------+---+---+--------------------------------------------
37c Tuler Butcher Failure model
38C-----------------------------------------------
39C I m p l i c i t T y p e s
40C-----------------------------------------------
41#include "implicit_f.inc"
42C---------+---------+---+---+--------------------------------------------
43C G l o b a l P a r a m e t e r s
44C-----------------------------------------------
45#include "units_c.inc"
46#include "comlock.inc"
47#include "com_xfem1.inc"
48C-----------------------------------------------
49C VAR | SIZE |TYP| RW| DEFINITION
50C---------+---------+---+---+--------------------------------------------
51C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
52C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
53C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
54C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
55C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
56C---------+---------+---+---+--------------------------------------------
57C TIME | 1 | F | R | CURRENT TIME
58C TIMESTEP| 1 | F | R | CURRENT TIME STEP
59C---------+---------+---+---+--------------------------------------------
60C SIGNXX | NEL | F | R | NEW ELASTO PLASTIC STRESS XX
61C SIGNYY | NEL | F | R | NEW ELASTO PLASTIC STRESS YY
62C ... | | | |
63C---------+---------+---+---+--------------------------------------------
64C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
65C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
66C---------+---------+---+---+--------------------------------------------
67C NGL ELEMENT ID
68C IPG CURRENT GAUSS POINT (in plane)
69C ILAY CURRENT LAYER
70C IPT CURRENT INTEGRATION POINT IN ALL LAYERS
71C IPTT CURRENT INTEGRATION POINT IN THE LAYER (FOR OUTPUT ONLY)
72C NPTT NUMBER OF INTEGRATION POINTS IN THE LAYER THICKNESS
73C NPTTF NUMBER OF FAILED INTEGRATION POINTS IN THE LAYER
74C NPTOT NUMBER OF INTEGRATION POINTS IN ALL LAYERS (TOTAL)
75C NOFF NUMBER OF FAILED INTEGRATION POINTS (TOTAL)
76C-----------------------------------------------
77C I N P U T A r g u m e n t s
78C-----------------------------------------------
79 INTEGER NEL,NUPARAM,IXFEM,IXEL,ILAY,IPTT,NPTT,NUVAR,IPT,NPTOT
80 INTEGER NGL(NEL),ELCRKINI(NXLAYMAX,*),NOFF(NEL)
81 my_real
82 . time,timestep(nel),uparam(nuparam),npttf(nel),
83 . signxx(nel),signyy(nel),signzz(nel),signxy(nel),signyz(nel),
84 . signzx(nel),tens(nel,5),dfmax(nel),tdel(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 uvar(nel,nuvar), off(nel),offl(nel)
89 my_real, DIMENSION(:) ,POINTER :: sigsav
90C-----------------------------------------------
91C L o c a l V a r i a b l e s
92C-----------------------------------------------
93 INTEGER I,J,ISHELL,NINDX,NINDXP,LAYXFEM,IBRIT
94 INTEGER INDX(NEL),INDXP(NEL),RFLAG(NEL),RFLAGP(NEL)
95 my_real
96 . tba,tbk,kk,sigmax,dadv,sigr_ini,sigr_adv,
97 . sig1,sig2,s1,s2,brit_b,brit_c,damd(nel)
98 CHARACTER (LEN=3) :: XCHAR
99C=======================================================================
100 tba = uparam(1)
101 tbk = uparam(2)
102 sigr_ini = uparam(3)
103 ishell = int(uparam(4))
104 ibrit = int(uparam(6))
105 brit_b = uparam(8)
106 brit_c = uparam(9)
107 dadv = uparam(10)
108c-----------------------------
109 layxfem = ixfem
110 IF (layxfem == 1 .and. ishell == 1) ishell = 2
111c-----------------------------
112 sigr_adv = sigr_ini*dadv
113C
114 DO i=1,nel
115 damd(i) = zero
116 rflag(i) = 0
117 rflagp(i)= 0
118 indxp(i) = 0
119 indx(i) = 0
120 ENDDO
121 nindx = 0
122 nindxp = 0
123C-----------------------------------------------
124 DO i=1,nel
125 tens(i,1) = signxx(i)
126 tens(i,2) = signyy(i)
127 tens(i,3) = signxy(i)
128 tens(i,4) = signyz(i)
129 tens(i,5) = signzx(i)
130 END DO
131C
132 IF (ixel > 0) THEN ! testing phantom elements
133 IF (ixel == 1) THEN
134 xchar = '1st'
135 ELSEIF (ixel == 2) THEN
136 xchar = '2nd'
137 ELSEIF (ixel == 3) THEN
138 xchar = '3rd'
139 ENDIF
140 ELSE
141 xchar = ' '
142 ENDIF
143c--------------------------
144 SELECT CASE (layxfem)
145c---------------
146 CASE (1) ! multilayer XFEM
147c---------------
148 DO i=1,nel
149 IF (off(i) == one) THEN
150 IF (uvar(i,1) < tbk) THEN
151 s1 = half*(signxx(i) + signyy(i))
152 s2 = half*(signxx(i) - signyy(i))
153 sig1 = s1 + sqrt(s2**2 + signxy(i)**2)
154 sig2 = s1 - sqrt(s2**2 + signxy(i)**2)
155 sigmax = max(sig1,sig2)
156c if(ngl(i)==7847)print*,'ilay,ippt,elcrkini=',ilay,ippt,elcrkini(ILAY,I)
157 IF (ixel == 0) THEN
158 IF (elcrkini(ilay,i)==0 .and. sigmax > sigr_ini) THEN
159 uvar(i,1) = uvar(i,1) + timestep(i)*(sigmax - sigr_ini)**tba
160c if(ngl(i)==7847)then
161c print*,'ilay,iptt=',ilay,iptt
162c print*,'SIGMAX,uvar=',SIGMAX,UVAR(I,1)
163c endif
164 IF (uvar(i,1) >= tbk) THEN
165 nindxp = nindxp + 1
166 indxp(nindxp) = i
167 offl(i) = zero
168 noff(i) = noff(i) + 1
169 npttf(i) = npttf(i) + one
170 rflagp(i) = 1
171 IF (int(npttf(i)) == nptt) THEN
172 nindx = nindx + 1
173 indx(nindx) = i
174 elcrkini(ilay,i) = -1 ! one layer failed (by initiation)
175 rflag(i) = 1
176 ENDIF
177 IF (noff(i) == nptot) THEN
178 off(i) = four_over_5
179 tdel(i)= time
180 ENDIF
181 ENDIF
182 ELSEIF (elcrkini(ilay,i) == 2 .and. sigmax > sigr_adv) THEN
183 uvar(i,1) = uvar(i,1) + timestep(i)*(sigmax - sigr_adv)**tba
184c print*,'ilay,iptt=',ilay,iptt
185c print*,'SIGMAX,uvar=',SIGMAX,UVAR(I,1)
186 IF (uvar(i,1) >= tbk) THEN
187 nindxp = nindxp + 1
188 indxp(nindxp) = i
189 offl(i) = zero
190 noff(i) = noff(i) + 1
191 npttf(i) = npttf(i)+ one
192 rflagp(i) = 1
193 IF (int(npttf(i)) == nptt) THEN
194 nindx = nindx+1
195 indx(nindx) = i
196 elcrkini(ilay,i) = 1 ! one layer failed (by advancing)
197 rflag(i) = -1
198 ENDIF
199 IF (noff(i) == nptot) THEN
200 off(i) = four_over_5
201 tdel(i)= time
202 ENDIF
203 ENDIF
204 ENDIF
205 ELSEIF (sigmax > sigr_ini) THEN ! IXEL > 0
206 uvar(i,1) = uvar(i,1) + timestep(i)*(sigmax - sigr_ini)**tba
207 IF (uvar(i,1) >= tbk) THEN
208 nindxp = nindxp + 1
209 indxp(nindxp) = i
210 offl(i) = zero
211 npttf(i) = npttf(i) + one
212 rflagp(i) = 1
213 IF (int(npttf(i)) == nptt) THEN
214 nindx = nindx+1
215 indx(nindx) = i
216 off(i) = four_over_5
217 tdel(i) = time
218 rflag(i) = 3
219 ENDIF
220 ENDIF
221 ENDIF ! IXEL
222 ELSE ! UVAR(I,1) > TBK
223 signxx(i) = zero
224 signyy(i) = zero
225 signxy(i) = zero
226 signyz(i) = zero
227 signzx(i) = zero
228 ENDIF !
229 ENDIF ! OFF(I) == ONE
230 ENDDO
231c-----------------------------------------------------------------------
232 IF (nindxp > 0) THEN
233 DO j=1,nindxp
234 i = indxp(j)
235 signxx(i) = zero
236 signyy(i) = zero
237 signxy(i) = zero
238 signyz(i) = zero
239 signzx(i) = zero
240#include "lockon.inc"
241 IF (rflagp(i) == 1) WRITE(iout, 3800)ngl(i),ilay,iptt
242 IF (rflagp(i) == 1) WRITE(istdo,3900)ngl(i),ilay,iptt,time
243#include "lockoff.inc"
244 ENDDO
245 ENDIF ! NINDXP > 0
246c
247 IF (nindx > 0) THEN
248 DO j=1,nindx
249#include "lockon.inc"
250 i = indx(j)
251c initialization std element
252 IF (rflag(i) == 1) WRITE(iout ,3000) ngl(i),ilay
253 IF (rflag(i) == 1) WRITE(istdo,3100) ngl(i),ilay,time
254c advancement std element
255 IF (rflag(i) == -1) WRITE(iout ,3200) ngl(i),ilay
256 IF (rflag(i) == -1) WRITE(istdo,3300) ngl(i),ilay,time
257c delete phantom element
258 IF (rflag(i) == 3) WRITE(iout, 3400)xchar,ngl(i),ilay
259 IF (rflag(i) == 3) WRITE(istdo,3500)xchar,ngl(i),ilay,time
260#include "lockoff.inc"
261 ENDDO
262 ENDIF ! NINDX > 0
263c
264c---------------
265 CASE (2) ! monolayer XFEM
266c---------------
267 IF (ishell == 1) THEN
268c---
269 DO i=1,nel
270 IF (off(i) == one) THEN
271 s1 = half*(signxx(i) + signyy(i))
272 s2 = half*(signxx(i) - signyy(i))
273 sig1 = s1 + sqrt(s2**2 + signxy(i)**2)
274 sig2 = s1 - sqrt(s2**2 + signxy(i)**2)
275 sigmax = max(sig1,sig2)
276C
277 IF (ibrit == 1) THEN ! ductile rupture
278 IF (ixel == 0) THEN
279 IF (elcrkini(ilay,i) == 0 .and. sigmax > sigr_ini) THEN
280 uvar(i,1) = uvar(i,1) + timestep(i)*(sigmax - sigr_ini)**tba
281 IF (uvar(i,1) >= tbk) THEN
282 elcrkini(ilay,i) = -1 ! ready for initiation
283 nindxp = nindxp + 1
284 indxp(nindxp) = i
285 nindx = nindx + 1
286 indx(nindx) = i
287 rflagp(i) = 1
288 rflag(i) = 1
289 tdel(i)= time
290 off(i) = four_over_5
291 ENDIF
292 ELSEIF (elcrkini(ilay,i) == 2 .and. sigmax > sigr_adv) THEN
293 uvar(i,1) = uvar(i,1) + timestep(i)*(sigmax - sigr_adv)**tba
294 IF (uvar(i,1) >= tbk) THEN
295 elcrkini(ilay,i) = 1 ! ready to advance
296 nindxp = nindxp + 1
297 indxp(nindxp) = i
298 nindx = nindx + 1
299 indx(nindx) = i
300 rflagp(i) = 1
301 rflag(i) = -1
302 tdel(i)= time
303 off(i) = four_over_5
304 ENDIF
305 ENDIF
306 ELSEIF (sigmax > sigr_ini) THEN ! IXEL > 0
307 uvar(i,1) = uvar(i,1) + timestep(i)*(sigmax - sigr_ini)**tba
308 IF (uvar(i,1) >= tbk) THEN
309 off(i) = four_over_5
310 nindxp = nindxp + 1
311 indxp(nindxp) = i
312 nindx = nindx + 1
313 indx(nindx) = i
314 rflagp(i) = 1
315 rflag(i) = 3
316 ENDIF
317 ENDIF ! IXEL
318c
319 ELSEIF (ibrit == 2) THEN ! brittle rupture
320 kk = one/max(em10,tbk)
321 IF (sigmax > uvar(i,1)) THEN
322 damd(i) = kk*(sigmax-uvar(i,1))**brit_b
323 ENDIF
324 noff(i) = noff(i) + damd(i)*timestep(i)
325 IF (noff(i)<one) uvar(i,1) = sigr_ini*(one-noff(i))**brit_c
326 IF (ixel == 0) THEN
327 IF (elcrkini(ilay,i) == 0) THEN ! ready for initiation
328 IF (noff(i) >= one) THEN
329 elcrkini(ilay,i) = -1
330 nindxp = nindxp + 1
331 indxp(nindxp) = i
332 nindx=nindx+1
333 indx(nindx)=i
334 rflagp(i) = 1
335 rflag(i) = 1
336 tdel(i)= time
337 off(i) = four_over_5
338 END IF
339 ELSEIF (elcrkini(ilay,i) == 2) THEN ! ready for advancing
340 IF (noff(i) >= one) THEN
341 elcrkini(ilay,i) = 1
342 nindxp = nindxp + 1
343 indxp(nindxp) = i
344 nindx=nindx+1
345 indx(nindx)=i
346 rflagp(i) = 1
347 rflag(i) = -1
348 tdel(i)= time
349 off(i) = four_over_5
350 ENDIF
351 ENDIF
352 ELSEIF (noff(i) >= one) THEN ! IXEL > 0
353 nindxp = nindxp + 1
354 indxp(nindxp) = i
355 nindx=nindx+1
356 indx(nindx)=i
357 rflagp(i) = 1
358 rflag(i) = 3
359 off(i) = four_over_5
360 END IF ! IXEL
361c
362 ENDIF ! IF(IBRIT==1
363 ENDIF ! IF(ISHELL==1
364 ENDDO
365c
366 ELSEIF (ishell == 2)THEN
367c
368 DO i=1,nel
369 IF (off(i) == one) THEN
370 IF (uvar(i,1) < tbk) THEN
371 s1 = half*(signxx(i) + signyy(i))
372 s2 = half*(signxx(i) - signyy(i))
373 sig1 = s1 + sqrt(s2**2 + signxy(i)**2)
374 sig2 = s1 - sqrt(s2**2 + signxy(i)**2)
375 sigmax = max(sig1,sig2)
376c
377 IF (ixel == 0) THEN
378 IF (elcrkini(ilay,i)==0 .and. sigmax > sigr_ini) THEN
379 uvar(i,1) = uvar(i,1) + timestep(i)*(sigmax - sigr_ini)**tba
380 IF (uvar(i,1) >= tbk) THEN ! ready for initiation
381 noff(i) = noff(i) + 1
382 nindxp = nindxp + 1
383 indxp(nindxp) = i
384 nindx = nindx+1
385 indx(nindx) = i
386 rflagp(i) = 1
387 IF (noff(i) == nptot) THEN
388 elcrkini(ilay,i) = -1 ! one layer failed (by initiation)
389 rflag(i) = 1
390 tdel(i)= time
391 off(i) = four_over_5
392 ENDIF
393 ENDIF
394 ELSEIF (elcrkini(ilay,i) == 2 .and. sigmax > sigr_adv) THEN
395 uvar(i,1) = uvar(i,1) + timestep(i)*(sigmax - sigr_adv)**tba
396 IF (uvar(i,1) >= tbk) THEN ! ready for advancing
397 noff(i) = noff(i) + 1
398 nindxp = nindxp + 1
399 indxp(nindxp) = i
400 nindx = nindx+1
401 indx(nindx) = i
402 rflagp(i) = 1
403 IF (noff(i) == nptot) THEN
404 elcrkini(ilay,i) = 1 ! one layer failed (by advancing)
405 rflag(i) = -1
406 tdel(i)= time
407 off(i) = four_over_5
408 ENDIF
409 ENDIF
410 ENDIF
411 ELSEIF (sigmax > sigr_ini) THEN ! IXEL > 0
412 uvar(i,1) = uvar(i,1) + timestep(i)*(sigmax - sigr_ini)**tba
413 IF (uvar(i,1) >= tbk) THEN ! IXEL > 0
414 nindxp = nindxp + 1
415 indxp(nindxp) = i
416 nindx = nindx+1
417 indx(nindx) = i
418 noff(i) = noff(i) + 1
419 IF (noff(i) == nptot) THEN
420 off(i) = four_over_5
421 rflag(i) = 3
422 ENDIF
423 ENDIF
424 ENDIF ! IXEL
425 ELSE ! UVAR(I,1) >= TBK
426 signxx(i) = zero
427 signyy(i) = zero
428 signxy(i) = zero
429 signyz(i) = zero
430 signzx(i) = zero
431 ENDIF ! IF(UVAR(I,1) < TBK)
432 ENDIF ! IF(ISHELL==2.AND.OFF(I)==ONE)
433 ENDDO
434c----
435 IF (nindxp > 0) THEN
436 DO j=1,nindxp
437 i = indxp(j)
438 signxx(i) = zero
439 signyy(i) = zero
440 signxy(i) = zero
441 signyz(i) = zero
442 signzx(i) = zero
443#include "lockon.inc"
444 IF (rflagp(i) == 1) WRITE(iout, 4800)ngl(i),iptt
445 IF (rflagp(i) == 1) WRITE(istdo,4900)ngl(i),iptt,time
446#include "lockoff.inc"
447 ENDDO
448 ENDIF ! NINDXP > 0
449c
450 IF (nindx > 0) THEN
451 DO j=1,nindx
452 i = indx(j)
453#include "lockon.inc"
454c initialization std element
455 IF (rflag(i) == 1) WRITE(iout, 4000) ngl(i)
456 IF (rflag(i) == 1) WRITE(istdo,4100) ngl(i),time
457c advancement std element
458 IF (rflag(i) == -1) WRITE(iout, 4200) ngl(i)
459 IF (rflag(i) == -1) WRITE(istdo,4300) ngl(i),time
460c delete phantom
461 IF (rflag(i) == 3) WRITE(iout, 4400)xchar,ngl(i)
462 IF (rflag(i) == 3) WRITE(istdo,4500)xchar,ngl(i),time
463#include "lockoff.inc"
464 ENDDO
465 ENDIF
466c----
467 ENDIF ! ishell
468c
469c-----------------
470 END SELECT ! LAYXFEM
471c-----------------
472c
473c--- Maximum Damage storing for output : 0 < DFMAX < 1
474 DO i=1,nel
475 dfmax(i) = min(one, max(dfmax(i),uvar(i,1)/tbk))
476 ENDDO
477C-----------------------------------------------
478 2000 FORMAT(1x,'FOR SHELL ELEMENT',i10,1x,'LAYER',i3,':',/
479 . 1x, 'STRESS TENSOR SET TO ZERO (TBUTCHER)')
480 2100 FORMAT(1x,'FOR SHELL ELEMENT',i10,1x,'LAYER',i3,':',/,
481 . 1x, 'STRESS TENSOR SET TO ZERO (TBUTCHER)',1x,'AT TIME :',1pe20.13)
482c--- multilayer xfem
483 3000 FORMAT(1x,'CRACK INITIALIZATION IN SHELL ELEMENT',i10,1x,'LAYER',i3)
484 3100 FORMAT(1x,'CRACK INITIALIZATION IN SHELL ELEMENT',i10,1x,'LAYER',i3,/
485 . 1x,'AT TIME :',1pe12.4)
486 3200 FORMAT(1x,'crack advancement in shell element',I10,' layer',I3)
487 3300 FORMAT(1X,'crack advancement in shell element',I10,' layer',I3/
488 . 1X,'at time :',1PE12.4)
489 3400 FORMAT(1X,'delete ',A4,' phantom element, shell id=',i10,' LAYER',i3)
490 3500 FORMAT(1x,'DELETE ',a4,' PHANTOM ELEMENT, SHELL ID=',i10,' LAYER',i3,/
491 . 1x,'AT TIME :',1pe12.4)
492 3800 FORMAT(1x,'TBUTCHER FAILURE IN SHELL',i10,1x,'LAYER',i3,1x,'INT POINT',i2)
493 3900 FORMAT(1x,'TBUTCHER FAILURE IN SHELL',i10,1x,'LAYER',i3,1x,'INT POINT',i2,/
494 . 1x,'AT TIME :',1pe12.4)
495c--- monolayer xfem
496 4000 FORMAT(1x,'CRACK INITIALIZATION IN SHELL ELEMENT',i10)
497 4100 FORMAT(1x,'CRACK INITIALIZATION IN SHELL ELEMENT',i10,/
498 . 1x,'AT TIME :',1pe12.4)
499 4200 FORMAT(1x,'CRACK ADVANCEMENT IN SHELL ELEMENT',i10)
500 4300 FORMAT(1x,'CRACK ADVANCEMENT IN SHELL ELEMENT',i10,/
501 . 1x,'AT TIME :',1pe12.4)
502 4400 FORMAT(1x,'DELETE ',a4,' PHANTOM ELEMENT, SHELL ID=',i10)
503 4500 FORMAT(1x,'DELETE ',a4,' PHANTOM ELEMENT, SHELL ID=',i10,/
504 . 1x,'AT TIME :',1pe12.4)
505 4800 FORMAT(1x,'TBUTCHER FAILURE IN SHELL',i10,1x,'INT POINT',i2)
506 4900 FORMAT(1x,'TBUTCHER FAILURE IN SHELL',i10,1x,'INT POINT',i2,1x,'AT TIME :',1pe12.4)
507c-----------
508 RETURN
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
id
initmumps id
Definition multiplerhs_example.m:4