OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_inievo_c.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_inievo_c ../engine/source/materials/fail/inievo/fail_inievo_c.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!||--- calls -----------------------------------------------------
29!|| table_vinterp ../engine/source/tools/curve/table_tools.F
30!||--- uses -----------------------------------------------------
31!|| interface_table_mod ../engine/share/modules/table_mod.F
32!|| message_mod ../engine/share/message_module/message_mod.F
33!|| table_mod ../engine/share/modules/table_mod.F
34!||====================================================================
35 SUBROUTINE fail_inievo_c (
36 1 NEL ,NUPARAM ,NUVAR ,
37 2 TABLE ,NTABLF ,ITABLF ,TIME ,UPARAM ,
38 3 NGL ,ALDT ,DPLA ,EPSP ,UVAR ,
39 4 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
40 5 PLA ,TEMP ,SIGY ,FOFF ,DFMAX ,
41 6 TDELE ,IPT ,IPG ,DMG_FLAG ,DMG_SCALE,
42 7 DAMINI ,AREA ,INLOC ,NPG )
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
481 END
alpha2
#define alpha2
Definition eval.h:48
alpha
#define alpha
Definition eval.h:35
fail_inievo_c
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)
Definition fail_inievo_c.F:43
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