OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_johnson_xfem.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| fail_johnson_xfem ../engine/source/materials/fail/johnson_cook/fail_johnson_xfem.f
25!||--- called by ------------------------------------------------------
26!|| mulawc ../engine/source/materials/mat_share/mulawc.F90
27!|| usermat_shell ../engine/source/materials/mat_share/usermat_shell.F
28!||--- uses -----------------------------------------------------
29!|| crackxfem_mod ../engine/share/modules/crackxfem_mod.F
30!||====================================================================
31 SUBROUTINE fail_johnson_xfem(
32 1 NEL ,NUPARAM ,UPARAM ,NUVAR ,UVAR ,
33 2 TIME ,TENS ,DPLA ,EPSP ,TSTAR ,
34 3 NGL ,IPT ,NPTOT ,NPTT ,NPTTF ,
35 4 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
36 5 OFF ,OFFL ,NOFF ,DFMAX ,TDEL ,
37 6 ELCRKINI ,IXFEM ,IXEL ,ILAY ,IPTT )
38C-----------------------------------------------
39c Johnson cook failure model
40C-----------------------------------------------
41C M o d u l e s
42C-----------------------------------------------
43 USE crackxfem_mod
44C-----------------------------------------------
45C I m p l i c i t T y p e s
46C-----------------------------------------------
47#include "implicit_f.inc"
48C-----------------------------------------------
49C G l o b a l P a r a m e t e r s
50C-----------------------------------------------
51#include "units_c.inc"
52#include "comlock.inc"
53#include "com_xfem1.inc"
54C---------+---------+---+---+--------------------------------------------
55C VAR | SIZE |TYP| RW| DEFINITION
56C---------+---------+---+---+--------------------------------------------
57C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
58C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
59C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
60C---------+---------+---+---+--------------------------------------------
61C NPTOT | 1 | I | R | NUMBER OF LAYERS OR INTEGRATION POINTS
62C IPT | 1 | I | R | LAYER OR INTEGRATION POINT NUMBER
63C IFLAG | * | I | R | GEOMETRICAL FLAGS
64C NGL | NEL | I | R | ELEMENT NUMBER
65C---------+---------+---+---+--------------------------------------------
66C TIME | 1 | F | R | CURRENT TIME
67C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
68C---------+---------+---+---+--------------------------------------------
69C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
70C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
71C---------+---------+---+---+--------------------------------------------
72C PLA | NEL | F |R/W| PLASTIC STRAIN
73C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
74C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
75C---------+---------+---+---+--------------------------------------------
76C ILAY CURRENT LAYER
77C IPT CURRENT INTEGRATION POINT IN ALL LAYERS
78C IPTT CURRENT INTEGRATION POINT IN THE LAYER (FOR OUTPUT ONLY)
79C NPTT NUMBER OF INTEGRATION POINTS IN THE LAYER THICKNESS
80C NPTTF NUMBER OF FAILED INTEGRATION POINTS IN THE LAYER
81C NPTOT NUMBER OF INTEGRATION POINTS IN ALL LAYERS (TOTAL)
82C NOFF NUMBER OF FAILED INTEGRATION POINTS (TOTAL)
83C---------+---------+---+---+--------------------------------------------
84C I N P U T A r g u m e n t s
85C-----------------------------------------------
86 INTEGER NEL, NUPARAM, NUVAR,IXFEM,IXEL,ILAY,NPTT,
87 . IPT,IPTT,NPTOT
88 INTEGER NGL(NEL),ELCRKINI(NXLAYMAX,*),NOFF(NEL)
89 my_real TIME,UPARAM(*),DPLA(NEL),NPTTF(NEL),
90 . EPSP(NEL),TSTAR(NEL),TENS(NEL,5)
91C-----------------------------------------------
92C I N P U T O U T P U T A r g u m e n t s
93C-----------------------------------------------
94 my_real uvar(nel,nuvar),off(nel),offl(nel),dfmax(nel),tdel(nel),
95 . signxx(nel),signyy(nel),signxy(nel),signyz(nel),signzx(nel)
96C-----------------------------------------------
97C L o c a l V a r i a b l e s
98C-----------------------------------------------
99 INTEGER I,J,K,L,JJ,ISHELL,NINDX,NINDXP,IADR,LAYXFEM
100 INTEGER INDX(NEL),INDXP(NEL),RFLAG(NEL),RFLAGP(NEL)
101 my_real D1,D2,D3,D4,D5,EPSP0,P,EPSF,SVM,THK_TMP,
102 . DADV,CC,BB,CR,ORM,SS1,SS2,EPSF_MIN
103 CHARACTER (LEN=3) :: XCHAR
104C=======================================================================
105 iadr = (ipt-1)*nel
106 xchar = ' '
107c
108 d1 = uparam(1)
109 d2 = uparam(2)
110 d3 = uparam(3)
111 d4 = uparam(4)
112 d5 = uparam(5)
113 epsp0 = uparam(6)
114 dadv = uparam(10)
115 epsf_min = uparam(12)
116 ishell = int(uparam(7))
117c-----------------------------
118 layxfem = ixfem
119 IF (layxfem == 1 .and. ishell == 1) ishell=2
120c-----------------------------
121 nindx = 0
122 nindxp = 0
123 DO i=1,nel
124 rflag(i) = 0
125 rflagp(i)= 0
126 indxp(i) = 0
127 indx(i) = 0
128 ENDDO
129C---------
130 DO i=1,nel
131 tens(i,1) = signxx(i)
132 tens(i,2) = signyy(i)
133 tens(i,3) = signxy(i)
134 tens(i,4) = signyz(i)
135 tens(i,5) = signzx(i)
136 END DO
137C
138 IF (ixel > 0) THEN ! testing phantom elements
139 IF (ixel == 1) THEN
140 xchar = '1st'
141 ELSEIF (ixel == 2) THEN
142 xchar = '2nd'
143 ELSEIF (ixel == 3) THEN
144 xchar = '3rd'
145 ENDIF
146 ENDIF
147C-------------------------------------------------------------------
148 SELECT CASE (layxfem)
149c---------------
150 CASE (1) ! multilayer Xfem (with Ishell = 2)
151C---------------
152 DO i=1,nel
153 IF (off(i) == one) THEN
154 IF (dfmax(i) < one) THEN
155 p = third*(signxx(i) + signyy(i))
156 svm = sqrt(signxx(i)*signxx(i)
157 . + signyy(i)*signyy(i)
158 . - signxx(i)*signyy(i)
159 . + three*signxy(i)*signxy(i))
160 epsf = d3*p/max(em20,svm)
161 epsf = (d1 + d2*exp(epsf))*(one
162 . + d4*log(max(one,epsp(i)/epsp0)))
163 . * (one + d5*tstar(i))
164 epsf = max(epsf,epsf_min)
165 IF (epsf > zero) dfmax(i) = dfmax(i) + dpla(i)/epsf
166c
167 IF (ixel == 0) THEN
168 IF (elcrkini(ilay,i)==0 .AND. dfmax(i) >= one) THEN
169 nindxp = nindxp+1
170 indxp(nindxp) = i
171 offl(i) = zero
172 noff(i) = noff(i) + 1
173 npttf(i) = npttf(i) + one
174 rflagp(i) = 1
175 IF (int(npttf(i)) == nptt) THEN
176 nindx = nindx+1
177 indx(nindx) = i
178 elcrkini(ilay,i) = -1 ! one layer failed (by initiation)
179 rflag(i) = 1
180 ENDIF
181 IF (noff(i) == nptot) THEN ! total element rupture
182 off(i) = four_over_5
183 tdel(i)= time
184 ENDIF
185 ELSEIF (elcrkini(ilay,i) == 2 .AND. dfmax(i) >= dadv) THEN
186 nindxp = nindxp + 1
187 indxp(nindxp) = i
188 offl(i) = zero
189 noff(i) = noff(i) + 1
190 npttf(i) = npttf(i)+ one
191 rflagp(i) = 1
192 IF (int(npttf(i)) == nptt) THEN
193 nindx = nindx+1
194 indx(nindx) = i
195 elcrkini(ilay,i) = 1 ! one layer failed(by advancing)
196 rflag(i) = -1
197 ENDIF
198 IF (noff(i) == nptot) THEN
199 off(i) = four_over_5
200 tdel(i)= time
201 ENDIF
202 ENDIF
203 ELSEIF (dfmax(i) >= one) THEN ! phantom element failure (IXEL > 0)
204 nindxp = nindxp + 1
205 indxp(nindxp) = i
206 offl(i) = zero
207 npttf(i) = npttf(i) + one
208 rflagp(i) = 1
209 IF (int(npttf(i)) == nptt) THEN
210 nindx = nindx+1
211 indx(nindx) = i
212 off(i) = four_over_5
213 tdel(i) = time
214 rflag(i) = 3
215 ENDIF
216 ENDIF ! ixel
217 ELSE ! DFMAX(I) >= ONE
218 signxx(i) = zero
219 signyy(i) = zero
220 signxy(i) = zero
221 signyz(i) = zero
222 signzx(i) = zero
223 ENDIF
224 ENDIF ! IF (OFF(I) == ONE)
225 ENDDO ! DO I=1,NEL
226c
227 IF (nindxp > 0) THEN
228 DO j=1,nindxp
229 i = indxp(j)
230 signxx(i) = zero
231 signyy(i) = zero
232 signxy(i) = zero
233 signyz(i) = zero
234 signzx(i) = zero
235#include "lockon.inc"
236 IF (rflagp(i) == 1) WRITE(iout, 3800)ngl(i),ilay,iptt
237 IF (rflagp(i) == 1) WRITE(istdo,3900)ngl(i),ilay,iptt,time
238#include "lockoff.inc"
239 ENDDO
240 ENDIF ! NINDXP > 0
241 IF (nindx > 0) THEN
242 DO j=1,nindx
243 i = indx(j)
244#include "lockon.inc"
245c initialization std element
246 IF (rflag(i) == 1) WRITE(iout ,3000) ngl(i),ilay
247 IF (rflag(i) == 1) WRITE(istdo,3100) ngl(i),ilay,time
248c advancement std element
249 IF (rflag(i) == -1) WRITE(iout ,3200) ngl(i),ilay
250 IF (rflag(i) == -1) WRITE(istdo,3300) ngl(i),ilay,time
251c delete phantom element
252 IF (rflag(i) == 3) WRITE(iout, 3400)xchar,ngl(i),ilay
253 IF (rflag(i) == 3) WRITE(istdo,3500)xchar,ngl(i),ilay,time
254#include "lockoff.inc"
255 ENDDO
256 ENDIF ! NINDX > 0
257c
258C---------------
259 CASE (2) ! monolayer Xfem
260C---------------
261 IF (ishell == 1) THEN
262c
263 DO i=1,nel
264 IF (off(i) == one) THEN
265 p = third*(signxx(i) + signyy(i))
266 svm = sqrt(signxx(i)*signxx(i)
267 . + signyy(i)*signyy(i)
268 . - signxx(i)*signyy(i)
269 . + three*signxy(i)*signxy(i))
270 epsf = d3*p/max(em20,svm)
271 epsf = (d1 + d2*exp(epsf))*(one
272 . + d4*log(max(one,epsp(i)/epsp0)))
273 . * (one + d5*tstar(i))
274 epsf = max(epsf,epsf_min)
275 IF (epsf > zero) dfmax(i) = dfmax(i) + dpla(i)/epsf
276c
277 IF (ixel == 0) THEN
278 IF (elcrkini(ilay,i)==0 .AND. dfmax(i)>=one) THEN
279 elcrkini(ilay,i) = -1 ! Initialisation
280 nindx = nindx + 1
281 indx(nindx) = i
282 rflag(i) = 1
283 tdel(i)= time
284 off(i) = four_over_5
285 ELSEIF (elcrkini(ilay,i) == 2 .AND.
286 . dfmax(i) >= dadv) THEN
287 elcrkini(ilay,i) = 1 ! Advancement
288 nindx = nindx + 1
289 indx(nindx) = i
290 rflag(i) =-1
291 tdel(i)= time
292 off(i) = four_over_5
293 ENDIF
294 ELSEIF (dfmax(i) >= one) THEN ! IXEL > 0
295 nindx = nindx+1
296 indx(nindx) = i
297 rflag(i) = 3
298 off(i) = four_over_5
299 ENDIF ! IXEL
300 ENDIF
301 ENDDO ! I=1,NEL
302
303 IF (nindx > 0) THEN
304 DO j=1,nindx
305 i = indx(j)
306#include "lockon.inc"
307 WRITE(iout, 4800)ngl(i),iptt
308 WRITE(istdo,4900)ngl(i),iptt,time
309c initialization std element
310 IF (rflag(i) == 1) WRITE(iout, 4000) ngl(i)
311 IF (rflag(i) == 1) WRITE(istdo,4100) ngl(i),time
312c advancement std element
313 IF (rflag(i) ==-1) WRITE(iout, 4200) ngl(i)
314 IF (rflag(i) ==-1) WRITE(istdo,4300) ngl(i),time
315c delete phantom
316 IF (rflag(i) == 3) WRITE(iout, 4400) xchar,ngl(i)
317 IF (rflag(i) == 3) WRITE(istdo,4500) xchar,ngl(i),time
318#include "lockoff.inc"
319 ENDDO
320 ENDIF
321c
322 ELSEIF (ishell == 2) THEN ! monolayer Xfem
323c
324 DO i=1,nel
325 IF (off(i) == one) THEN
326 IF (dfmax(i) < one) THEN
327 p = third*(signxx(i) + signyy(i))
328 svm = sqrt(signxx(i)*signxx(i)
329 . + signyy(i)*signyy(i)
330 . - signxx(i)*signyy(i)
331 . + three*signxy(i)*signxy(i))
332 epsf = d3*p/max(em20,svm)
333 epsf = (d1 + d2*exp(epsf))*(one
334 . + d4*log(max(one,epsp(i)/epsp0)))
335 . * (one + d5*tstar(i))
336 epsf = max(epsf,epsf_min)
337 IF (epsf > zero) dfmax(i) = dfmax(i) + dpla(i)/epsf
338c
339 IF (ixel == 0) THEN
340 IF (elcrkini(ilay,i)==0 .and. dfmax(i) >= one) THEN
341 nindxp = nindxp + 1
342 indxp(nindxp) = i
343 noff(i) = noff(i) + 1
344 offl(i) = zero
345 rflagp(i) = 1
346 IF (noff(i) == nptot) THEN ! all integration points failed
347 nindx = nindx + 1
348 indx(nindx) = i
349 elcrkini(ilay,i) = -1 ! initialize crack
350 rflag(i) = 1
351 tdel(i)= time
352 off(i) = four_over_5
353 ENDIF
354 ELSEIF (elcrkini(ilay,i) == 2 .AND. dfmax(i) >= dadv) THEN
355 nindxp = nindxp + 1
356 indxp(nindxp) = i
357 noff(i) = noff(i) + 1
358 offl(i) = zero
359 dfmax(i) = one
360 rflagp(i) = 1
361 IF (noff(i) == nptot) THEN ! all integration points failed
362 nindx = nindx + 1
363 indx(nindx) = i
364 elcrkini(ilay,i) = 1 ! advance crack
365 rflag(i) = -1
366 tdel(i)= time
367 off(i) = four_over_5
368 ENDIF
369 ENDIF
370 ELSEIF (dfmax(i) >= one) THEN ! phantom element failure (IXEL > 0)
371 nindxp = nindxp + 1
372 indxp(nindxp) = i
373 noff(i) = noff(i) + 1
374 rflagp(i) = 1
375 IF (noff(i) == nptot) THEN
376 nindx = nindx + 1
377 indx(nindx) = i
378 off(i) = four_over_5
379 rflag(i) = 3
380 ENDIF
381 ENDIF ! IXEL
382 ELSE ! DFMAX >= 1
383 signxx(i) = zero
384 signyy(i) = zero
385 signxy(i) = zero
386 signyz(i) = zero
387 signzx(i) = zero
388 ENDIF
389 ENDIF ! IF (OFF(I) == ONE)
390 ENDDO ! DO I=1,NEL
391C
392 IF (nindxp > 0) THEN
393 DO j=1,nindxp
394 i = indxp(j)
395 signxx(i) = zero
396 signyy(i) = zero
397 signxy(i) = zero
398 signyz(i) = zero
399 signzx(i) = zero
400#include "lockon.inc"
401 IF (rflagp(i) == 1) WRITE(iout, 4800)ngl(i),iptt
402 IF (rflagp(i) == 1) WRITE(istdo,4900)ngl(i),iptt,time
403#include "lockoff.inc"
404 ENDDO
405 ENDIF ! NINDXP > 0
406 IF (nindx > 0) THEN
407 DO j=1,nindx
408 i = indx(j)
409#include "lockon.inc"
410c initialization std element
411 IF (rflag(i) == 1) WRITE(iout ,4000) ngl(i)
412 IF (rflag(i) == 1) WRITE(istdo,4100) ngl(i),time
413c advancement std element
414 IF (rflag(i) == -1) WRITE(iout, 4200) ngl(i)
415 IF (rflag(i) == -1) WRITE(istdo,4300) ngl(i),time
416c delete phantom
417 IF (rflag(i) == 3) WRITE(iout, 4400) xchar,ngl(i)
418 IF (rflag(i) == 3) WRITE(istdo,4500) xchar,ngl(i),time
419#include "lockoff.inc"
420 ENDDO
421 ENDIF ! NINDX > 0
422
423 ENDIF ! ISHELL
424c-----------------
425 END SELECT ! LAYXFEM
426c
427C---------Damage for output 0 < DFMAX < 1 --------------------
428 DO i=1,nel
429 dfmax(i) = min(one,dfmax(i))
430 ENDDO
431C------------------
432 2000 FORMAT(1x,'FOR SHELL ELEMENT (JOHNS)',i10,1x,'LAYER',i3,':',/
433 . 1x, 'STRESS TENSOR SET TO ZERO')
434 2100 FORMAT(1x,'FOR SHELL ELEMENT (JOHNS)',i10,1x,'LAYER',i3,':',/,
435 . 1x,'STRESS TENSOR SET TO ZERO',1x,'AT TIME :',1pe12.4)
436 2400 FORMAT(1x,1pg20.13,' % OF THICKNESS OF SHELL BROKEN ')
437 2410 FORMAT(1x,1pg20.13,' % OF THICKNESS OF SHELL ',i10,' BROKEN ')
438 2500 FORMAT(1x,' LOWER SKIN -> UPPER SKIN ')
439 2600 FORMAT(1x,' UPPER SKIN -> LOWER SKIN ')
440c--- multilayer
441 3000 FORMAT(1x,'CRACK INITIALIZATION IN SHELL ELEMENT',i10,1x,'LAYER',i3)
442 3100 FORMAT(1x,'CRACK INITIALIZATION IN SHELL ELEMENT',i10,1x,'LAYER',i3,/
443 . 1x,'AT TIME :',1pe12.4)
444 3200 FORMAT(1x,'CRACK ADVANCEMENT IN SHELL ELEMENT',i10,' LAYER',i3)
445 3300 FORMAT(1x,'CRACK ADVANCEMENT IN SHELL ELEMENT',i10,' LAYER',i3/
446 . 1x,'AT TIME :',1pe12.4)
447 3400 FORMAT(1x,'DELETE ',a4,' PHANTOM ELEMENT, SHELL ID=',i10,' LAYER',i3)
448 3500 FORMAT(1x,'DELETE ',a4,' PHANTOM ELEMENT, SHELL ID=',i10,' layer',I3,/
449 . 1X,'at time :',1PE12.4)
450 3800 FORMAT(1X,'jc failure in shell',I10,1X,'layer',I3,1X,'int point',I2)
451 3900 FORMAT(1X,'jc failure in shell',I10,1X,'layer',I3,1X,'int point',I2,/
452 . 1X,'at time :',1PE12.4)
453c--- monolayer
454 4000 FORMAT(1X,'crack initialization in shell element',I10)
455 4100 FORMAT(1X,'crack initialization in shell element',I10,1X,'at time :',1PE12.4)
456 4200 FORMAT(1X,'crack advancement in shell element',I10)
457 4300 FORMAT(1X,'crack advancement in shell element',I10,1X,'at time :',1PE12.4)
458 4400 FORMAT(1X,'delete ',A4,' phantom element, shell id=',I10)
459 4500 FORMAT(1X,'delete ',A4,' phantom element, shell id=',I10,1X,'at time :',1PE12.4)
460 4800 FORMAT(1X,'jc failure in shell',I10,1X,'int point',I2)
461 4900 FORMAT(1X,'jc failure in shell',I10,1X,'int point',I2,1X,'at time :',1PE12.4)
462c------------------
463 RETURN
464 END
fail_johnson_xfem
subroutine fail_johnson_xfem(nel, nuparam, uparam, nuvar, uvar, time, tens, dpla, epsp, tstar, ngl, ipt, nptot, nptt, npttf, signxx, signyy, signxy, signyz, signzx, off, offl, noff, dfmax, tdel, elcrkini, ixfem, ixel, ilay, iptt)
Definition fail_johnson_xfem.F:38
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
crackxfem_mod
Definition crackxfem_mod.F:121
jc
subroutine jc(p, t, a, b, cm, cn, tref, tm, epsm, sigmam, jc_yield, tan_jc)
Definition sigeps106.F:339