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

Go to the source code of this file.

Functions/Subroutines

subroutine fail_inievo_c (nel, nuparam, nuvar, table, ntablf, itablf, time, uparam, ngl, aldt, dpla, epsp, uvar, signxx, signyy, signxy, signyz, signzx, pla, temp, sigy, foff, dfmax, tdele, ipt, ipg, dmg_flag, dmg_scale, damini, area, inloc, npg)

Function/Subroutine Documentation

◆ fail_inievo_c()

subroutine fail_inievo_c ( integer, intent(in) nel,
integer, intent(in) nuparam,
integer, intent(in) nuvar,
type (ttable), dimension(ntable), intent(in) table,
integer, intent(in) ntablf,
integer, dimension(ntablf), intent(in) itablf,
intent(in) time,
dimension(nuparam), intent(in) uparam,
integer, dimension(nel), intent(in) ngl,
dimension(nel), intent(in) aldt,
dimension(nel), intent(in) dpla,
dimension(nel), intent(in) epsp,
dimension(nel,nuvar), intent(inout) uvar,
dimension(nel), intent(in) signxx,
dimension(nel), intent(in) signyy,
dimension(nel), intent(in) signxy,
dimension(nel), intent(in) signyz,
dimension(nel), intent(in) signzx,
dimension(nel), intent(in) pla,
dimension(nel), intent(in) temp,
dimension(nel), intent(in) sigy,
integer, dimension(nel), intent(inout) foff,
dimension(nel), intent(inout) dfmax,
dimension(nel), intent(inout) tdele,
integer, intent(in) ipt,
integer, intent(in) ipg,
integer, intent(inout) dmg_flag,
dimension(nel), intent(inout) dmg_scale,
dimension(nel), intent(inout) damini,
dimension(nel), intent(in) area,
integer, intent(in) inloc,
integer, intent(in) npg )

Definition at line 35 of file fail_inievo_c.F.

43C---------+---------+---+---+--------------------------------------------
44 USE table_mod
45 USE interface_table_mod
46 USE message_mod
47C-----------------------------------------------
48C I m p l i c i t T y p e s
49C-----------------------------------------------
50#include "implicit_f.inc"
51C-----------------------------------------------
52C C o m m o n B l o c k s
53C-----------------------------------------------
54#include "units_c.inc"
55#include "comlock.inc"
56#include "com04_c.inc"
57#include "tabsiz_c.inc"
58C-----------------------------------------------
59C I N P U T A r g u m e n t s
60C-----------------------------------------------
61 INTEGER, INTENT(IN) ::
62 . NEL,NUPARAM,NUVAR,NGL(NEL),IPT,IPG,NTABLF,
63 . INLOC,NPG
64 INTEGER, DIMENSION(NTABLF) ,INTENT(IN) :: ITABLF
65 INTEGER, INTENT(INOUT) ::
66 . DMG_FLAG,FOFF(NEL)
67 my_real, INTENT(IN) ::
68 . time,uparam(nuparam),aldt(nel),
69 . dpla(nel),epsp(nel),temp(nel),
70 . signxx(nel),signyy(nel),signxy(nel),
71 . signyz(nel),signzx(nel),pla(nel),
72 . sigy(nel),area(nel)
73 my_real, INTENT(INOUT) ::
74 . uvar(nel,nuvar),dfmax(nel),tdele(nel),
75 . dmg_scale(nel),damini(nel)
76 TYPE (TTABLE), INTENT(IN), DIMENSION(NTABLE) :: TABLE
77C-----------------------------------------------
78C L o c a l V a r i a b l e s
79C-----------------------------------------------
80 INTEGER I,J,INDX(NEL),NINDX,IPOS(NEL,2),
81 . NINIEVO,ISHEAR,ILEN
82 INTEGER, DIMENSION(:), ALLOCATABLE ::
83 . INITYPE,EVOTYPE,EVOSHAP,COMPTYP,TAB_ID,
84 . TAB_EL,FCRIT
85 my_real, DIMENSION(:), ALLOCATABLE ::
86 . sr_ref,fscale,ini_p1,el_ref,elscal,
87 . disp,ener,alpha2
88 my_real, DIMENSION(:,:), ALLOCATABLE ::
89 . dmgini,dmgevo
90 my_real
91 . l0(nel) ,triax(nel) ,epsf(nel) ,depsf(nel) ,
92 . xvec(nel,2) ,sizefac(nel),maxshear(nel),epsmod(nel) ,
93 . p(nel) ,svm(nel) ,dmgmax(nel) ,dmgmul(nel) ,
94 . center ,radius ,sigp1 ,sigp2 ,
95 . devsp1 ,devsp2 ,alpha(nel) ,plas_disp ,
96 . yld0 ,sigpmaj(nel),dsize(nel) ,fac ,
97 . df ,sxx ,syy ,szz
98C!--------------------------------------------------------------
99 !=======================================================================
100 ! - INITIALISATION OF COMPUTATION ON TIME STEP
101 !=======================================================================
102 ! Recovering failure criterion parameters
103 ninievo = uparam(1)
104 ishear = int(uparam(2))
105 ilen = int(uparam(3))
106 IF (inloc > 0) ilen = 1
107 ALLOCATE(initype(ninievo))
108 ALLOCATE(evotype(ninievo))
109 ALLOCATE(evoshap(ninievo))
110 ALLOCATE(comptyp(ninievo))
111 ALLOCATE(tab_id(ninievo))
112 ALLOCATE(sr_ref(ninievo))
113 ALLOCATE(fscale(ninievo))
114 ALLOCATE(ini_p1(ninievo))
115 ALLOCATE(tab_el(ninievo))
116 ALLOCATE(el_ref(ninievo))
117 ALLOCATE(elscal(ninievo))
118 ALLOCATE(disp(ninievo))
119 ALLOCATE(ener(ninievo))
120 ALLOCATE(alpha2(ninievo))
121 tab_id(1:ninievo) = itablf(1:ninievo)
122 tab_el(1:ninievo) = itablf(ninievo+1:ninievo*2)
123c
124 DO j = 1,ninievo
125 initype(j) = uparam(6 + 14*(j-1))
126 evotype(j) = uparam(7 + 14*(j-1))
127 evoshap(j) = uparam(8 + 14*(j-1))
128 comptyp(j) = uparam(9 + 14*(j-1))
129c TAB_ID(J) = NINT(UPARAM(10 + 14*(J-1)))
130 sr_ref(j) = uparam(11 + 14*(j-1))
131 fscale(j) = uparam(12 + 14*(j-1))
132 ini_p1(j) = uparam(13 + 14*(j-1))
133c TAB_EL(J) = NINT(UPARAM(14 + 14*(J-1)))
134 el_ref(j) = uparam(15 + 14*(j-1))
135 elscal(j) = uparam(16 + 14*(j-1))
136 disp(j) = uparam(17 + 14*(j-1))
137 alpha2(j) = uparam(18 + 14*(j-1))
138 ener(j) = uparam(19 + 14*(j-1))
139 ENDDO
140c
141 ! Set flag for stress softening
142 dmg_flag = 1
143c
144 ! Element characteristic length computation
145 ! -> Initial values
146 IF (uvar(1,1) == zero) THEN
147 ! -> Critical timestep formulation
148 IF (ilen == 1) THEN
149 uvar(1:nel,1) = aldt(1:nel)
150 ! -> Geometric formulation
151 ELSE
152 DO i = 1,nel
153 uvar(i,1) = sqrt(npg*area(i))
154 ENDDO
155 ENDIF
156 ENDIF
157 ! -> Current geometric formulation
158 IF (ilen == 2) THEN
159 DO i = 1,nel
160 uvar(i,1) = sqrt(npg*area(i))
161 ENDDO
162 ENDIF
163 l0(1:nel) = uvar(1:nel,1)
164 ! Positive stress triaxiality bounded plastic strain
165 epsmod(1:nel) = uvar(1:nel,2)
166 ! Damage initiation and evolution variable
167 ALLOCATE(dmgini(nel,ninievo))
168 ALLOCATE(dmgevo(nel,ninievo))
169 ALLOCATE(fcrit(nel))
170 DO j = 1,ninievo
171 DO i=1,nel
172 ! Initiation damage
173 dmgini(i,j) = uvar(i,3+(j-1)*3)
174 ! Evolution damage
175 dmgevo(i,j) = uvar(i,4+(j-1)*3)
176 ENDDO
177 ENDDO
178 ! Criterion number leading to element deletion
179 fcrit(1:nel) = 0
180c
181 !====================================================================
182 ! - LOOP OVER THE ELEMENT TO COMPUTE THE STRESS STATE QUANTITIES
183 !====================================================================
184 DO i=1,nel
185c
186 ! Computation of hydrostatic stress, Von Mises stress, and stress triaxiality
187 p(i) = -third*(signxx(i) + signyy(i))
188 ! Von Mises equivalent stress
189 sxx = signxx(i) + p(i)
190 syy = signyy(i) + p(i)
191 szz = p(i)
192 svm(i) = half*(sxx**2 + syy**2 + szz**2)
193 . + signxy(i)**2
194 IF (ishear > 0) THEN
195 svm(i) = svm(i) + signzx(i)**2 + signyz(i)**2
196 ENDIF
197 svm(i) = sqrt(three*svm(i))
198 triax(i) = -p(i)/max(em20,svm(i))
199 IF (triax(i) < -two_third) triax(i) = -two_third
200 IF (triax(i) > two_third) triax(i) = two_third
201c
202 ! Increase the modified plastic strain
203 IF (triax(i) > zero) epsmod(i) = epsmod(i) + dpla(i)
204c
205 ! Compute principal stresses
206 center = half*(signxx(i)+signyy(i))
207 radius = sqrt((half*(signxx(i)-signyy(i)))**2 + signxy(i)**2)
208 sigp1 = center + radius
209 sigp2 = center - radius
210 sigpmaj(i) = sigp1
211 maxshear(i) = (sigp1-sigp2)*half
212c
213 ! Compute the alpha parameter for FLD/MSFLD
214 devsp1 = sigp1 - third*(sigp1+sigp2)
215 devsp2 = sigp2 - third*(sigp1+sigp2)
216 alpha(i) = devsp2/sign(max(abs(devsp1),em20),devsp1)
217 ENDDO
218c
219 !====================================================================
220 ! - COMPUTE DAMAGE INITIATION AND EVOLUTION
221 !====================================================================
222 DO j = 1,ninievo
223 ! damage initiation type selection
224 SELECT CASE(initype(j))
225 ! Plastic strain vs triaxiality
226 CASE(1)
227 xvec(1:nel,1) = triax(1:nel)
228 ! Plastic strain vs shear influence (theta)
229 CASE(2)
230 DO i = 1,nel
231 xvec(i,1) = (svm(i) + ini_p1(j)*p(i))/max(maxshear(i),em08)
232 ENDDO
233 ! MSFLD / FLD
234 CASE(3,4)
235 xvec(1:nel,1) = alpha(1:nel)
236 ! Normalized principal stress
237 CASE(5)
238 DO i = 1,nel
239 xvec(i,1) = (svm(i) + ini_p1(j)*p(i))/max(sigpmaj(i),em08)
240 ENDDO
241 END SELECT
242 xvec(1:nel,2) = epsp(1:nel)/sr_ref(j)
243 ipos(1:nel,1:2) = 1
244 CALL table_vinterp(table(tab_id(j)),nel,nel,ipos,xvec,epsf,depsf)
245 epsf(1:nel) = epsf(1:nel)*fscale(j)
246c
247 ! Compute the element size regularization factor
248 IF (tab_el(j) > 0) THEN
249 xvec(1:nel,1) = l0(1:nel)/el_ref(j)
250 SELECT CASE (initype(j))
251 CASE(1)
252 xvec(1:nel,2) = triax(1:nel)
253 CASE(2)
254 DO i = 1,nel
255 xvec(i,2) = (svm(i) + ini_p1(j)*p(i))/max(maxshear(i),em08)
256 ENDDO
257 CASE(3,4)
258 xvec(1:nel,2) = alpha(1:nel)
259 CASE(5)
260 DO i = 1,nel
261 xvec(i,2) = (svm(i) + ini_p1(j)*p(i))/max(sigpmaj(i),em08)
262 ENDDO
263 END SELECT
264 ipos(1:nel,1:2) = 1
265 CALL table_vinterp(table(tab_el(j)),nel,nel,ipos,xvec,sizefac,dsize)
266 sizefac(1:nel) = sizefac(1:nel)*elscal(j)
267 epsf(1:nel) = epsf(1:nel)*sizefac(1:nel)
268 ENDIF
269c
270 ! Update damage initiation
271 SELECT CASE (initype(j))
272 CASE(1,2,5)
273 DO i = 1,nel
274 IF ((dpla(i) > zero).AND.(dmgini(i,j)<one).AND.(foff(i) /= 0)) THEN
275 dmgini(i,j) = dmgini(i,j) + dpla(i)/max(epsf(i),em20)
276 dmgini(i,j) = min(dmgini(i,j),one)
277 ENDIF
278 ENDDO
279 CASE(3)
280 IF (nint(ini_p1(j))>0) THEN
281 DO i = 1,nel
282 IF (((epsmod(i)-uvar(i,2)) > zero).AND.(dmgini(i,j)<one).AND.(foff(i) /= 0)) THEN
283 dmgini(i,j) = dmgini(i,j) + (epsmod(i)-uvar(i,2))/max(epsf(i),em20)
284 dmgini(i,j) = min(dmgini(i,j),one)
285 ENDIF
286 ENDDO
287 ELSE
288 DO i = 1,nel
289 IF (((epsmod(i)-uvar(i,2)) > zero).AND.(dmgini(i,j)<one).AND.(foff(i) /= 0)) THEN
290 dmgini(i,j) = max(dmgini(i,j),epsmod(i)/max(epsf(i),em20))
291 dmgini(i,j) = min(dmgini(i,j),one)
292 ENDIF
293 ENDDO
294 ENDIF
295 CASE(4)
296 IF (nint(ini_p1(j))>0) THEN
297 DO i = 1,nel
298 IF ((dpla(i) > zero).AND.(dmgini(i,j)<one).AND.(foff(i) /= 0)) THEN
299 dmgini(i,j) = dmgini(i,j) + dpla(i)/max(epsf(i),em20)
300 dmgini(i,j) = min(dmgini(i,j),one)
301 ENDIF
302 ENDDO
303 ELSE
304 DO i = 1,nel
305 IF ((dpla(i) > zero).AND.(dmgini(i,j)<one).AND.(foff(i) /= 0)) THEN
306 dmgini(i,j) = max(dmgini(i,j),pla(i)/max(epsf(i),em20))
307 dmgini(i,j) = min(dmgini(i,j),one)
308 ENDIF
309 ENDDO
310 ENDIF
311 END SELECT
312c
313 ! update damage evolution
314 SELECT CASE (evotype(j))
315 ! Plastic displacement at failure
316 CASE(1)
317 SELECT CASE (evoshap(j))
318 ! Linear shape
319 CASE(1)
320 DO i = 1,nel
321 IF ((dmgini(i,j) >= one).AND.(dpla(i)>zero).AND.
322 . (foff(i) /= 0).AND.(dmgevo(i,j)<one)) THEN
323 dmgevo(i,j) = dmgevo(i,j) + l0(i)*dpla(i)/disp(j)
324 dmgevo(i,j) = min(one,dmgevo(i,j))
325 IF (dmgevo(i,j) >= one) fcrit(i) = j
326 ENDIF
327 ENDDO
328 ! Exponential shape
329 CASE(2)
330 DO i = 1,nel
331 IF ((dmgini(i,j) >= one).AND.(dpla(i)>zero).AND.
332 . (foff(i) /= 0).AND.(dmgevo(i,j)<one)) THEN
333 IF (dmgevo(i,j) == zero) uvar(i,5+(j-1)*3) = pla(i)
334 plas_disp = (pla(i) - uvar(i,5+(j-1)*3))*l0(i)/disp(j)
335 dmgevo(i,j) = dmgevo(i,j) + (alpha2(j)/(one - exp(-alpha2(j))))*
336 . exp(-alpha2(j)*plas_disp)*
337 . dpla(i)*l0(i)/disp(j)
338 IF (dmgevo(i,j) > 0.999d0) dmgevo(i,j) = one
339 dmgevo(i,j) = min(one,dmgevo(i,j))
340 IF (dmgevo(i,j) >= one) fcrit(i) = j
341 ENDIF
342 ENDDO
343 END SELECT
344 ! Fracture energy failure
345 CASE(2)
346 SELECT CASE (evoshap(j))
347 ! Linear shape
348 CASE(1)
349 DO i = 1,nel
350 IF ((dmgini(i,j) >= one).AND.(dpla(i)>zero).AND.
351 . (foff(i) /= 0).AND.(dmgevo(i,j)<one)) THEN
352 IF (dmgevo(i,j) == zero) uvar(i,5+(j-1)*3) = sigy(i)
353 yld0 = uvar(i,5+(j-1)*3)
354 dmgevo(i,j) = dmgevo(i,j) + dpla(i)*l0(i)*yld0/(two*ener(j))
355 dmgevo(i,j) = min(one,dmgevo(i,j))
356 IF (dmgevo(i,j) >= one) fcrit(i) = j
357 ENDIF
358 ENDDO
359 ! Exponential shape
360 CASE(2)
361 DO i = 1,nel
362 IF ((dmgini(i,j) >= one).AND.(dpla(i)>zero).AND.
363 . (foff(i) /= 0).AND.(dmgevo(i,j)<one)) THEN
364 uvar(i,5+(j-1)*3) = uvar(i,5+(j-1)*3) + sigy(i)*l0(i)*dpla(i)
365 dmgevo(i,j) = one - exp(-(uvar(i,5+(j-1)*3))/ener(j))
366 IF (dmgevo(i,j) > 0.999d0) dmgevo(i,j) = one
367 dmgevo(i,j) = min(one,dmgevo(i,j))
368 IF (dmgevo(i,j) >= one) fcrit(i) = j
369 ENDIF
370 ENDDO
371 END SELECT
372 ! Failure criterion approach
373 CASE DEFAULT
374 DO i = 1,nel
375 IF ((dmgini(i,j) >= one).AND.(dpla(i)>zero).AND.
376 . (foff(i) /= 0).AND.(dmgevo(i,j)<one)) THEN
377 dmgevo(i,j) = dmgini(i,j)
378 dmgevo(i,j) = min(one,dmgevo(i,j))
379 IF (dmgevo(i,j) >= one) fcrit(i) = j
380 ENDIF
381 ENDDO
382 END SELECT
383 ENDDO
384c
385 !====================================================================
386 ! - COMPUTE GLOBAL DAMAGE VARIABLE AND DAMAGE SCALING
387 !====================================================================
388 dfmax(1:nel) = zero
389 dmgmax(1:nel) = zero
390 dmgmul(1:nel) = one
391 DO j = 1,ninievo
392 SELECT CASE (comptyp(j))
393 ! Maximum damage
394 CASE(1)
395 DO i = 1,nel
396 dmgmax(i) = max(dmgmax(i),dmgevo(i,j))
397 ENDDO
398 ! Multiplicative damage
399 CASE(2)
400 DO i = 1,nel
401 dmgmul(i) = dmgmul(i)*(one-dmgevo(i,j))
402 ENDDO
403 END SELECT
404 ENDDO
405 dmgmul(1:nel) = one - dmgmul(1:nel)
406 nindx = 0
407 indx(1:nel) = 0
408 DO i = 1,nel
409 IF (foff(i) /= 0) THEN
410 dfmax(i) = max(dmgmax(i),dmgmul(i))
411 IF (dfmax(i) >= one) THEN
412 nindx = nindx + 1
413 indx(nindx) = i
414 foff(i) = 0
415 tdele(i) = time
416 ENDIF
417 ENDIF
418 ENDDO
419c
420 !====================================================================
421 ! - UPDATE THE DAMAGE SCALING FACTOR
422 !====================================================================
423 dmg_scale(1:nel) = one - dfmax(1:nel)
424c
425 !====================================================================
426 ! - UPDATE THE USER VARIABLE
427 !====================================================================
428 ! Positive stress triaxiality bounded plastic strain
429 uvar(1:nel,2) = epsmod(1:nel)
430 damini(1:nel) = zero
431 DO j = 1,ninievo
432 ! Checking element failure and recovering user variable
433 DO i=1,nel
434 ! Damage initiation output
435 damini(i) = max(dmgini(i,j),damini(i))
436 ! Initiation damage
437 uvar(i,3+(j-1)*3) = dmgini(i,j)
438 ! Evolution damage
439 uvar(i,4+(j-1)*3) = dmgevo(i,j)
440 ENDDO
441 ENDDO
442c
443 !====================================================================
444 ! - PRINTOUT DATA ABOUT FAILED ELEMENTS
445 !====================================================================
446 IF (nindx > 0) THEN
447 DO j=1,nindx
448 i = indx(j)
449#include "lockon.inc"
450 WRITE(iout, 1000) ngl(i),fcrit(i),ipg,ipt,time
451 WRITE(istdo,1000) ngl(i),fcrit(i),ipg,ipt,time
452#include "lockoff.inc"
453 END DO
454 END IF
455c
456 !====================================================================
457 ! - TABLES DEALLOCATION
458 !====================================================================
459 IF (ALLOCATED(initype)) DEALLOCATE(initype)
460 IF (ALLOCATED(evotype)) DEALLOCATE(evotype)
461 IF (ALLOCATED(evoshap)) DEALLOCATE(evoshap)
462 IF (ALLOCATED(comptyp)) DEALLOCATE(comptyp)
463 IF (ALLOCATED(tab_id)) DEALLOCATE(tab_id)
464 IF (ALLOCATED(sr_ref)) DEALLOCATE(sr_ref)
465 IF (ALLOCATED(fscale)) DEALLOCATE(fscale)
466 IF (ALLOCATED(ini_p1)) DEALLOCATE(ini_p1)
467 IF (ALLOCATED(tab_el)) DEALLOCATE(tab_el)
468 IF (ALLOCATED(el_ref)) DEALLOCATE(el_ref)
469 IF (ALLOCATED(elscal)) DEALLOCATE(elscal)
470 IF (ALLOCATED(disp)) DEALLOCATE(disp)
471 IF (ALLOCATED(ener)) DEALLOCATE(ener)
472 IF (ALLOCATED(alpha2)) DEALLOCATE(alpha2)
473 IF (ALLOCATED(dmgini)) DEALLOCATE(dmgini)
474 IF (ALLOCATED(dmgevo)) DEALLOCATE(dmgevo)
475 IF (ALLOCATED(fcrit)) DEALLOCATE(fcrit)
476c------------------------
477 1000 FORMAT(1x,'FOR SHELL ELEMENT NUMBER ',i10,
478 . ' FAILURE (INIEVO) WITH CRITERION NUMBER ',i3,
479 . ' AT GAUSS POINT ',i3,' LAYER ',i3,' AT TIME :',1pe12.4)
480c
my_real
#define my_real
Definition cppsort.cpp:32
alpha2
#define alpha2
Definition eval.h:48
alpha
#define alpha
Definition eval.h:35
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
interface_table_mod::table_vinterp
Definition table_mod.F:251
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
interface_table_mod
Definition table_mod.F:238
message_mod
Definition message_mod.F:1249
table_mod
Definition table_mod.F:112