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

Go to the source code of this file.

Functions/Subroutines

subroutine fail_wind_xfem (nel, nuparam, nuvar, uparam, uvar, time, timestep, ssp, aldt, crklen, elcrkini, tdel, off, offly, signxx, signyy, signxy, signyz, signzx, ngl, ixel, ilay, ipt, npt, ixfem, iorth, dir_a, dir1_crk, dir2_crk, crkdir)

Function/Subroutine Documentation

◆ fail_wind_xfem()

subroutine fail_wind_xfem ( integer nel,
integer nuparam,
integer nuvar,
uparam,
uvar,
time,
timestep,
dimension(nel) ssp,
dimension(nel) aldt,
dimension(nel) crklen,
integer, dimension(nxlaymax,nel) elcrkini,
dimension(nel) tdel,
off,
integer, dimension(nel) offly,
dimension(nel) signxx,
dimension(nel) signyy,
dimension(nel) signxy,
dimension(nel) signyz,
dimension(nel) signzx,
integer, dimension(nel) ngl,
integer ixel,
integer ilay,
integer ipt,
integer npt,
integer ixfem,
integer iorth,
dir_a,
target dir1_crk,
target dir2_crk,
crkdir )

Definition at line 29 of file fail_wind_xfem.F.

37C-----------------------------------------------
38c Windshield failure model for XFEM (ref : PhD thesis Christian Alter 2018)
39c now : only for monolayer xfem
40C-----------------------------------------------
41C I m p l i c i t T y p e s
42C-----------------------------------------------
43#include "implicit_f.inc"
44C-----------------------------------------------
45C G l o b a l P a r a m e t e r s
46C-----------------------------------------------
47#include "units_c.inc"
48#include "comlock.inc"
49#include "com_xfem1.inc"
50#include "mvsiz_p.inc"
51C---------+---------+---+---+--------------------------------------------
52C VAR | SIZE |TYP| RW| DEFINITION
53C---------+---------+---+---+--------------------------------------------
54C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
55C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
56C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
57C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
58C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
59C---------+---------+---+---+--------------------------------------------
60C TIME | 1 | F | R | CURRENT TIME
61C TIMESTEP| 1 | F | R | CURRENT TIME STEP
62C---------+---------+---+---+--------------------------------------------
63C SIGNXX | NEL | F | R | STRESS XX
64C SIGNYY | NEL | F | R | STRESS YY
65C ... | | | |
66C---------+---------+---+---+--------------------------------------------
67C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
68C---------+---------+---+---+--------------------------------------------
69C ILAY CURRENT LAYER
70C IPT CURRENT INTEGRATION POINT IN THE LAYER (FOR OUTPUT ONLY)
71C---------+---------+---+---+--------------------------------------------
72C I N P U T A r g u m e n t s
73C-----------------------------------------------
74 INTEGER NEL,NUPARAM,NUVAR,IXFEM,IORTH,IXEL,ILAY,IPT,NPT
75 INTEGER ELCRKINI(NXLAYMAX,NEL),NGL(NEL)
76 my_real time,timestep
77 my_real, DIMENSION(NEL) :: ssp,aldt,tdel,crklen,
78 . signxx,signyy,signzz,signxy,signyz,signzx
79 my_real uparam(nuparam),dir_a(nel,2),crkdir(nel,2),dir1_crk(nxlaymax,mvsiz),
80 . dir2_crk(nxlaymax,nel)
81 TARGET :: dir1_crk,dir2_crk
82C-----------------------------------------------
83C I N P U T O U T P U T A r g u m e n t s
84C-----------------------------------------------
85 INTEGER OFFLY(NEL)
86 my_real uvar(nel,nuvar),off(nel)
87C-----------------------------------------------
88C L o c a l V a r i a b l e s
89C-----------------------------------------------
90 INTEGER I,J,NINDX,NINDXD,IDEB,ISRATE
91 INTEGER , DIMENSION(NEL) :: INDX,INDXD,RFLAG
92 my_real sigmax,aa,bb,cc,s1,s2,s3,cr,k_ic,k_th,v0,vc,tfact,reflen,
93 . cosi,cosx,cosy,sinx,siny,cos2,sin2,dmg1,dmg3,
94 . sp1,sp2,sp3,geored,cr_foil,cr_air,cr_core,cr_edge,alpha,alphai,exp_n,exp_m
95 my_real, DIMENSION(NEL) :: dadv,tpropg,ai,formf,sigp_akt,dmg_scale,
96 . sigp1,sigp2,sigdt1,sigdt2,sig_dtf1,sig_dtf2,sigp_min,sigp_max
97 my_real, DIMENSION(:), POINTER :: dir1,dir2
98 CHARACTER (LEN=3) :: XCHAR
99c--------------------------------------------------------------------
100c! cm(6): initial crack depth in contact to foil
101c! cm(7): initial crack depth in contact to air
102c! cm(8): initial crack depth on edge
103c! cm(9): K_IC - fracture toughness from literature
104c! cm(10): K_TH - fatigue threshold from literature
105c--------------------------------------------------------------------
106c UVAR(1) = SIGP1 ! principal stress DIR1 ! uvar 17
107c UVAR(1) = SIGP2 ! principal stress DIR2 ! uvar 18
108c UVAR(3) = SIGDT1 ! principal stress rate DIR1
109c UVAR(4) = SIGDT2 ! principal stress rate DIR2
110c UVAR(5) = SIGP_MIN ! uvar 24
111c UVAR(6) = SIGP_MAX ! uvar 14
112c UVAR(7) = FORMF ! uvar 25
113c UVAR(8) = ai ! uvar 26
114c UVAR(9) = sigp_akt
115c UVAR(10) = edge element flag
116c UVAR(11) = reduction factor flag
117c UVAR(12) = DADV
118c UVAR(13) = SIG_DTF1
119c UVAR(14) = SIG_DTF2
120C=======================================================================
121 xchar = ' '
122 IF (ixel > 0) THEN ! testing phantom elements
123 IF (ixel == 1) THEN
124 xchar = '1st'
125 ELSEIF (ixel == 2) THEN
126 xchar = '2nd'
127 ELSEIF (ixel == 3) THEN
128 xchar = '3rd'
129 ENDIF
130 ELSE
131 xchar = 'standard'
132 ENDIF
133c
134 rflag(1:nel) = 0
135 indx(1:nel) = 0
136 nindx = 0
137 nindxd = 0
138c--------------------------------------------------------------------
139 exp_n = uparam(1)
140 cr_foil = uparam(2)
141 cr_air = uparam(3)
142 cr_core = uparam(4)
143 cr_edge = uparam(5)
144 k_ic = uparam(6)
145 k_th = uparam(7)
146 v0 = uparam(8)
147 vc = uparam(9)
148 alpha = uparam(10)
149 geored = uparam(11)
150 reflen = nint(uparam(14))
151 israte = nint(uparam(16))
152 ideb = nint(uparam(17))
153c
154 alphai = one - alpha
155 exp_m = one / (one + exp_n)
156c---
157c IF (IXEL == 0) THEN
158 DO i = 1,nel
159 dadv(i) = min(one, geored / sqrt(aldt(i)/reflen) ) ! reduction factor for advancement
160 tpropg(i) = crklen(i) / min(vc,ssp(i))
161 tpropg(i) = max(tpropg(i), timestep)
162 uvar(i,12)= dadv(i)
163 END DO
164c ENDIF
165c--------------------------------------------------------------------
166 IF (time == zero) then
167c Crack depth modification and definition of Y (formfactor)
168 DO i=1,nel
169 IF (uvar(i,10) == one) THEN ! edge-element
170 ai(i) = cr_edge
171 formf(i) = 1.12
172 ELSE
173 formf(i) = one
174 IF (ipt == 1) THEN
175 ai(i) = cr_foil ! crack depth at foil side in mm (unit system: length)
176 ELSEIF (ipt == npt) THEN
177 ai(i) = cr_air ! crack depth at top side in mm (unit system: length)
178 ELSE
179 ai(i) = cr_core ! crack depth inside glass
180 ENDIF
181 ENDIF
182 sigp_akt(i) = (two*(exp_n + one)*k_ic**exp_n) /
183 . ((exp_n - two)*v0*(formf(i)*sqrt(pi))**exp_n*
184 . ai(i)**((exp_n-two)/two))
185 sigp_akt(i) = sigp_akt(i)**exp_m
186c
187 sigp_min(i) = k_th / (formf(i)*sqrt(pi*ai(i)))
188 sigp_max(i) = k_ic / (formf(i)*sqrt(pi*ai(i)))
189 uvar(i,5) = sigp_min(i)
190 uvar(i,6) = sigp_max(i)
191 uvar(i,7) = formf(i)
192 uvar(i,8) = ai(i)
193 uvar(i,9) = sigp_akt(i)
194 ENDDO
195 ELSE
196 sigp_min(1:nel) = uvar(1:nel,5)
197 sigp_max(1:nel) = uvar(1:nel,6)
198 formf(1:nel) = uvar(1:nel,7)
199 ai(1:nel) = uvar(1:nel,8)
200 sigp_akt(1:nel) = uvar(1:nel,9)
201 ENDIF
202c--------------------------------------------------------------------
203c Principal stress
204c--------------------------------------------------------------------
205 DO i = 1,nel
206 aa = (signxx(i) + signyy(i))*half
207 bb = (signxx(i) - signyy(i))*half
208 cc = signxy(i)
209 cr = sqrt(bb**2 + cc**2)
210 s1 = aa + cr
211 s2 = aa - cr
212 sigp1(i) = s1
213 sigp2(i) = s2
214c
215c IF (S12 == ZERO .and. S2 >= ZERO) THEN
216c MAJANG = ZERO
217c ELSEIF (S12 == ZERO .and. S2 < ZERO) THEN
218c MAJANG = HALF*PI
219c ELSE
220c MAJANG = ATAN((SIGP1(I) - SIGNXX(I)) / S12)
221c ENDIF
222 END DO
223c--------------------------------------------------------------------
224c Calculate stress rate + filtering
225c--------------------------------------------------------------------
226 IF (israte == 0) THEN
227c exponential moving average with smoohting coefficient = alpha
228 DO i = 1,nel
229 sigdt1(i) = abs(sigp1(i) - uvar(i,1)) / max(timestep, em20)
230 sigdt2(i) = abs(sigp2(i) - uvar(i,2)) / max(timestep, em20)
231 sigdt1(i) = alpha * sigdt1(i) + alphai * uvar(i,3)
232 sigdt2(i) = alpha * sigdt2(i) + alphai * uvar(i,4)
233 END DO
234 ELSE
235c arythmetic moving average over 50 cycles
236 DO i = 1,nel
237 IF (off(i) == one) THEN
238 sigdt1(i) = abs(sigp1(i) - uvar(i,1)) / max(timestep, em20)
239 uvar(i,21) = sigdt1(i)
240 DO j = 1,49
241 sigdt1(i) = sigdt1(i) + uvar(i,21+j)
242 ENDDO
243 sigdt1(i) = sigdt1(i) / 50
244 DO j = 70,22,-1
245 uvar(i,j) = uvar(i,j-1)
246 ENDDO
247c
248 sigdt2(i) = abs(sigp2(i) - uvar(i,2)) / max(timestep, em20)
249 uvar(i,71) = sigdt2(i)
250 DO j = 1,49
251 sigdt2(i) = sigdt2(i) + uvar(i,71+j)
252 ENDDO
253 sigdt2(i) = sigdt2(i) / 50
254 DO j = 120,72,-1
255 uvar(i,j) = uvar(i,j-1)
256 ENDDO
257 ELSE
258 sigdt1(i) = zero
259 sigdt2(i) = zero
260 ENDIF
261 END DO
262 ENDIF
263c--------------------------------------------------------------------
264c Max strength calculation
265c--------------------------------------------------------------------
266 IF (ipt == npt) THEN ! Top layer integration point
267 DO i = 1,nel
268 sig_dtf1(i) = sigp_akt(i) *abs(sigdt1(i))**exp_m
269 sig_dtf2(i) = sigp_akt(i) *abs(sigdt2(i))**exp_m
270 sig_dtf1(i) = max(sig_dtf1(i),sigp_min(i))
271 sig_dtf1(i) = min(sig_dtf1(i),sigp_max(i))
272 sig_dtf2(i) = max(sig_dtf2(i),sigp_min(i))
273 sig_dtf2(i) = min(sig_dtf2(i),sigp_max(i))
274 END DO
275 ELSEIF (ipt < npt) THEN ! inner integration points
276 DO i = 1,nel
277 IF (uvar(i,10) == one) THEN ! edge element => stress rate dependent
278 sig_dtf1(i) = sigp_akt(i)*abs(sigdt1(i))**exp_m
279 sig_dtf2(i) = sigp_akt(i)*abs(sigdt2(i))**exp_m
280 sig_dtf1(i) = max(sig_dtf1(i),sigp_min(i))
281 sig_dtf1(i) = min(sig_dtf1(i),sigp_max(i))
282 sig_dtf2(i) = max(sig_dtf2(i),sigp_min(i))
283 sig_dtf2(i) = min(sig_dtf2(i),sigp_max(i))
284 ELSE
285 sig_dtf1(i) = sigp_max(i)
286 sig_dtf2(i) = sigp_max(i)
287 ENDIF
288 END DO
289 ENDIF
290c--------------------------
291 IF (ixel == 0) THEN ! original element
292 DO i = 1,nel
293 IF (off(i) == one) THEN
294 IF (elcrkini(ilay,i) == 0) THEN
295 ! no cracks in the neibourhood => initialization
296 IF (sigp1(i) > sig_dtf1(i) .and. tdel(i) == zero .and. offly(i) == 1) THEN
297 ! start damage
298 tdel(i) = time
299 offly(i) = 0
300 elcrkini(ilay,i) = -5 ! for crack length estimation
301 nindx = nindx+1
302 indx(nindx) = i
303 rflag(i) = 1
304 ENDIF
305 IF (tdel(i) > 0) THEN ! damage started already
306 tfact = (time - tdel(i)) / tpropg(i)
307 tfact = min(one, tfact)
308 dmg_scale(i) = one - tfact
309 nindxd = nindxd+1
310 indxd(nindxd) = i
311 IF (tfact >= one) THEN ! end of damage in normal direction
312 elcrkini(ilay,i) = -1 ! failure completed (by initiation)
313 off(i) = four_over_5
314 IF (iorth > 0) THEN
315 cosx = dir_a(i,1)
316 sinx = dir_a(i,2)
317 cosy = crkdir(i,1)
318 siny = crkdir(i,2)
319 dir1_crk(ilay,i) = cosx*cosy - sinx*siny
320 dir2_crk(ilay,i) = cosx*siny + sinx*cosy
321 ELSE
322 dir1_crk(ilay,i) = crkdir(i,1)
323 dir2_crk(ilay,i) = crkdir(i,2)
324 ENDIF
325 ELSE
326 ENDIF
327 ENDIF
328c
329 ELSEIF (elcrkini(ilay,i) == 2) THEN
330 ! neighbor element cut at common edge => advancement
331 uvar(i,11) = one ! reduction factor flag for output
332 sig_dtf1(i) = sig_dtf1(i)*dadv(i) ! apply reduction factor to criteria
333 sig_dtf2(i) = sig_dtf2(i)*dadv(i)
334 sig_dtf1(i) = max(sig_dtf1(i),sigp_min(i)*dadv(i))
335 sig_dtf1(i) = min(sig_dtf1(i),sigp_max(i)*dadv(i))
336 sig_dtf2(i) = max(sig_dtf2(i),sigp_min(i)*dadv(i))
337 sig_dtf2(i) = min(sig_dtf2(i),sigp_max(i)*dadv(i))
338 IF (sigp1(i) > sig_dtf1(i) .and. tdel(i) == zero .and. offly(i) == 1) THEN
339 tdel(i) = time
340 offly(i) = 0
341 elcrkini(ilay,i) = 5 ! tag for crack length estimation
342 nindx = nindx+1
343 indx(nindx) = i
344 rflag(i) =-1
345 ENDIF
346 IF (offly(i) == 0) THEN ! damage started already
347 tfact = (time - tdel(i)) / tpropg(i)
348 nindxd = nindxd+1
349 indxd(nindxd) = i
350 IF (tfact >= one) THEN
351 elcrkini(ilay,i) = 1 ! advance existing crack
352 dmg_scale(i) = zero
353 off(i) = four_over_5
354 IF (iorth > 0) THEN
355 cosx = dir_a(i,1)
356 sinx = dir_a(i,2)
357 cosy = crkdir(i,1)
358 siny = crkdir(i,2)
359 dir1_crk(ilay,i) = cosx*cosy - sinx*siny
360 dir2_crk(ilay,i) = cosx*siny + sinx*cosy
361 ELSE
362 dir1_crk(ilay,i) = crkdir(i,1)
363 dir2_crk(ilay,i) = crkdir(i,2)
364 ENDIF
365 ELSE
366 dmg_scale(i) = one - tfact
367 ENDIF
368 ENDIF
369 ENDIF
370 ENDIF ! OFF == 1
371 ENDDO
372c--------------------------
373c calculate and save principal stress direction when first failure occurs
374c--------------------------
375 IF (nindx > 0) THEN
376 DO j=1,nindx
377 i = indx(j)
378 s1 = signxy(i)
379 s2 = sigp1(i) - signxx(i)
380 cr = sqrt(s1**2 + s2**2)
381 IF (cr > zero) THEN
382 crkdir(i,1) = s1 / cr
383 crkdir(i,2) = s2 / cr
384 ELSEIF (s1 > s2) THEN
385 crkdir(i,1) = zero
386 crkdir(i,2) = one
387 ELSE
388 crkdir(i,1) = one
389 crkdir(i,2) = zero
390 ENDIF
391 ENDDO
392c
393 IF (iorth > 0) THEN ! stress is in orthotropic frame (DIR_A)
394 DO j=1,nindx
395 i = indx(j)
396 cosx = dir_a(i,1)
397 sinx = dir_a(i,2)
398 cosy = crkdir(i,1)
399 siny = crkdir(i,2)
400 dir1_crk(ilay,i) = cosx*cosy - sinx*siny
401 dir2_crk(ilay,i) = cosx*siny + sinx*cosy
402 ENDDO
403 ELSE ! stress is in local element frame
404 DO j=1,nindx
405 i = indx(j)
406 dir1_crk(ilay,i) = crkdir(i,1)
407 dir2_crk(ilay,i) = crkdir(i,2)
408 ENDDO
409 ENDIF
410 ENDIF
411c
412c apply progressive damage and turn the stress tensor back to the local system
413c
414 IF (nindxd > 0) THEN
415 DO j=1,nindxd
416 i = indxd(j)
417 s1 = signxx(i)
418 s2 = signyy(i)
419 s3 = signxy(i)
420 cosx = crkdir(i,1)
421 sinx = crkdir(i,2)
422 cos2 = cosx * cosx
423 sin2 = sinx * sinx
424 cosi = cosx * sinx
425c rotate stress to previously saved crack direction
426 sp1 = cos2*s1 + sin2*s2 + two*cosi*s3
427 sp2 = sin2*s1 + cos2*s2 - two*cosi*s3
428 sp3 = cosi*(s2 - s1) + (cos2 - sin2)*s3
429
430 dmg1 = dmg_scale(i)
431 dmg3 = min(one, 0.6 + 0.4 * dmg1)
432c stress reduction
433 sp1 = sp1 * dmg1
434 sp3 = sp3 * dmg3
435c rotate reduced stress back to current element coordinate system
436 signxx(i) = cos2*sp1 + sin2*sp2 - two*cosi*sp3
437 signyy(i) = sin2*sp1 + cos2*sp2 + two*cosi*sp3
438 signxy(i) = cosi*(sp1 - sp2) + (cos2 - sin2)*sp3
439 ENDDO
440c
441 if (ideb==1) then
442 DO j=1,nindxd
443 i = indxd(j)
444 write(iout,'(A,I10,3E16.9)') 'DAMAGE ELEMENT N, CRLEN,TPROPG,DMG_SCALE=',
445 . ngl(i),crklen(i),tpropg(i),dmg_scale(i)
446 ENDDO
447 endif
448 ENDIF
449c
450 ELSEIF (ixel > 0) THEN ! phantom element
451 DO i = 1,nel
452 IF (off(i) == one) THEN
453 IF (elcrkini(ilay,i) == 2) THEN ! second crack arrived => delete phantom element
454 nindx = nindx+1
455 indx(nindx) = i
456 rflag(i) = 3
457 off(i) = four_over_5
458 ELSE
459 sigp_min(i) = k_th / (sqrt(pi*cr_foil))
460 sigp_max(i) = k_ic / (sqrt(pi*cr_foil))
461 sig_dtf1(i) = sigp_max(i) * dadv(i)
462 IF (sigp1(i) > sig_dtf1(i)) THEN
463 ! phantom element reaches criteria
464c delete phantom but should transfer information to remaining neighbor element
465 nindx = nindx+1
466 indx(nindx) = i
467 rflag(i) = 3
468 off(i) = four_over_5
469 ENDIF
470 ENDIF
471 ENDIF
472 ENDDO
473 ENDIF
474c-----------------------------------------------
475 DO i=1,nel
476 uvar(i,1) = sigp1(i)
477 uvar(i,2) = sigp2(i)
478 uvar(i,3) = sigdt1(i)
479 uvar(i,4) = sigdt2(i)
480 END DO
481c-----------------------------------------------
482 IF (nindx > 0) THEN
483 DO j=1,nindx
484 i = indx(j)
485#include "lockon.inc"
486 WRITE(iout, 3000)ngl(i),ilay,ipt
487 WRITE(istdo,3100)ngl(i),ilay,ipt,time
488c initialization std element
489 IF (rflag(i) == 1) WRITE(iout ,4000)
490 IF (rflag(i) == 1) WRITE(istdo,4000)
491c advancement std element
492 IF (rflag(i) == -1) WRITE(iout, 4200)
493 IF (rflag(i) == -1) WRITE(istdo,4200)
494 if (ideb==1 .and. abs(rflag(i))==1) then
495 write(iout,'(A,5E16.9)') 'SIGP1,SIG_DTF1,SIGP_AKT,SIGP_MAX,SIGDT1=',
496 . sigp1(i),sig_dtf1(i),sigp_akt(i),sigp_max(i),sigdt1(i)
497 write(iout,'(A,5E16.9)') 'SIGP2,SIG_DTF2,SIGP_AKT,SIGP_MIN,SIGDT2=',
498 . sigp2(i),sig_dtf2(i),sigp_akt(i),sigp_min(i),sigdt2(i)
499 if (rflag(i)==-1) then
500 write(iout,'(A,E16.9)') 'advancement => crit reduction, DADV =',dadv(i)
501 endif
502 endif
503c delete phantom
504 IF (rflag(i) == 3) WRITE(iout, 4400) xchar,ngl(i)
505 IF (rflag(i) == 3) WRITE(istdo,4500) xchar,ngl(i),time
506#include "lockoff.inc"
507 ENDDO
508 ENDIF
509c-----------------------------------------------
510 3000 FORMAT(1x,'FAILURE IN SHELL',i10,1x,'LAYER',i2,1x,'INT POINT',i2)
511 3100 FORMAT(1x,'FAILURE IN SHELL',i10,1x,'LAYER',i2,1x,'INT POINT',i2,
512 . 1x,'AT TIME ',1pe12.4)
513 4000 FORMAT(10x,'CRACK INITIALIZATION')
514 4200 FORMAT(10x,'CRACK ADVANCEMENT')
515 4400 FORMAT(1x,'DELETE ',a4,' PHANTOM ELEMENT, SHELL ID=',i10)
516 4500 FORMAT(1x,'DELETE ',a4,' PHANTOM ELEMENT, SHELL ID=',i10,/
517 . 1x,'AT TIME :',1pe12.4)
518c-----------------------------------------------
519 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21