OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_orthenerg_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_orthenerg_c ../engine/source/materials/fail/orthenerg/fail_orthenerg_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!||====================================================================
29 SUBROUTINE fail_orthenerg_c(
30 1 NEL ,NUPARAM ,NUVAR ,UPARAM ,UVAR ,
31 2 NGL ,TIME ,IPG ,ILAY ,IPT ,
32 3 DEPSXX ,DEPSYY ,DEPSXY ,DMG_FLAG ,DMG_SCALE,
33 4 ALDT ,FOFF ,DFMAX ,TDEL ,
34 5 SIGNXX ,SIGNYY ,SIGNXY ,IGTYP ,PLY_ID )
35C-----------------------------------------------
36C I m p l i c i t T y p e s
37C-----------------------------------------------
38#include "implicit_f.inc"
39C-----------------------------------------------
40C G l o b a l P a r a m e t e r s
41C-----------------------------------------------
42#include "units_c.inc"
43#include "comlock.inc"
44C-----------------------------------------------
45C I N P U T A r g u m e n t s
46C-----------------------------------------------
47 INTEGER, INTENT(IN) :: NEL,NUPARAM,NUVAR,IPG,ILAY,IPT,
48 . IGTYP,PLY_ID
49 INTEGER, DIMENSION(NEL), INTENT(IN) :: NGL
50 my_real, INTENT(IN) :: TIME
51 my_real, DIMENSION(NEL), INTENT(IN) :: DEPSXX,DEPSYY,
52 . depsxy,signxx,signyy,signxy,aldt
53 my_real, DIMENSION(NUPARAM), INTENT(IN) :: uparam
54C-----------------------------------------------
55C I N P U T O U T P U T A r g u m e n t s
56C-----------------------------------------------
57 INTEGER, INTENT(OUT) :: DMG_FLAG
58 INTEGER, DIMENSION(NEL), INTENT(INOUT) :: FOFF
59 my_real, DIMENSION(NEL), INTENT(INOUT) :: DFMAX,TDEL
60 my_real, DIMENSION(NEL,5), INTENT(OUT) :: DMG_SCALE
61 my_real, DIMENSION(NEL,NUVAR), INTENT(INOUT) :: UVAR
62C-----------------------------------------------
63C L o c a l V a r i a b l e s
64C-----------------------------------------------
65 INTEGER I,J,NINDX,FAILMOD,MODE,ISHAP11T,ISHAP11C,ISHAP22T,
66 . ishap22c,ishap12t,ishap12c,nmod
67 INTEGER ,DIMENSION(NEL) :: INDX
68 INTEGER ,DIMENSION(NEL,6) :: FMODE
69 my_real DAM(NEL,6),ENER(NEL,6),LE(NEL)
70 my_real sigma_11t,sigma_11c,sigma_22t,sigma_22c,
71 . sigma_12t,sigma_12c,g_11t,g_11c,g_22t,g_22c,g_12t,g_12c
72C
73 !=======================================================================
74 ! - INITIALISATION OF COMPUTATION ON TIME STEP
75 !=======================================================================
76 ! Recovering failure criterion parameters
77 sigma_11t = uparam(3) ! -> Critical stress for tension in direction 11
78 sigma_11c = uparam(4) ! -> Critical stress for compression in direction 11
79 sigma_22t = uparam(5) ! -> Critical stress for tension in direction 22
80 sigma_22c = uparam(6) ! -> Critical stress for compression in direction 22
81 sigma_12t = uparam(9) ! -> Critical stress for positive shear in plane 12
82 sigma_12c = uparam(10) ! -> Critical stress for negative shear in plane 12
83 g_11t = uparam(15) ! -> Fracture energy for tension in direction 11
84 g_11c = uparam(16) ! -> Fracture energy for compression in direction 11
85 g_22t = uparam(17) ! -> Fracture energy for tension in direction 22
86 g_22c = uparam(18) ! -> Fracture energy for compression in direction 22
87 g_12t = uparam(21) ! -> fracture energy for positive shear in plane 12
88 g_12c = uparam(22) ! -> Fracture energy for negative shear in plane 12
89 ishap11t = nint(uparam(27)) ! -> Softening shape flag for tension in direction 11
90 ishap11c = nint(uparam(28)) ! -> Softening shape flag for compression in direction 11
91 ishap22t = nint(uparam(29)) ! -> Softening shape flag for tension in direction 22
92 ishap22c = nint(uparam(30)) ! -> Softening shape flag for compression in direction 22
93 ishap12t = nint(uparam(33)) ! -> softening shape flag for positive shear in plane 12
94 ishap12c = nint(uparam(34)) ! -> Softening shape flag for negative shear in plane 12
95 nmod = nint(uparam(39)) ! -> Number of failed modes prior to integration point failure
96 nmod = min(nmod,6)
97C
98 ! Set DMG_FLAG to Orthotropic
99 dmg_flag = 3
100C
101 ! Save initial element length
102 IF (uvar(1,13) == zero) THEN
103 uvar(1:nel,13) = aldt(1:nel)
104 ENDIF
105 le(1:nel) = uvar(1:nel,13)
106c
107 ! Recover user variable value
108 DO j = 1,6
109 DO i = 1,nel
110 ! Damage variable per mode
111 dam(i,j) = uvar(i,j)
112 ! Dissipated energy per mode
113 ener(i,j) = uvar(i,j+6)
114 ENDDO
115 ENDDO
116c
117 !====================================================================
118 ! - COMPUTATION OF THE DAMAGE VARIABLE EVOLUTION
119 !====================================================================
120 ! Initialization of element failure index
121 nindx = 0
122c
123 ! Loop over the elements
124 DO i=1, nel
125c
126 IF (foff(i) == 1) THEN
127c
128 ! Initialization of failed modes index
129 fmode(i,1:6) = 0
130 failmod = 0
131c
132 ! Mode 1 failure - Tension XX
133 mode = 1
134 IF (signxx(i) > sigma_11t .AND. signxx(i) > zero .AND. dam(i,mode) < one) THEN
135 ! -> Linear stress softening
136 IF (ishap11t == 1) THEN
137 dam(i,mode) = dam(i,mode) + le(i)*depsxx(i)*sigma_11t/(two*g_11t)
138 dam(i,mode) = max(dam(i,mode),uvar(i,mode))
139 ! -> Exponential stress softening
140 ELSEIF (ishap11t == 2) THEN
141 ener(i,mode) = ener(i,mode) + signxx(i)*le(i)*depsxx(i)
142 ener(i,mode) = max(ener(i,mode),uvar(i,mode+6))
143 dam(i,mode) = (one - exp(-ener(i,mode)/g_11t))
144 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
145 ENDIF
146 dam(i,mode) = min(dam(i,mode),one)
147 ! Output damage variable
148 dfmax(i) = max(dfmax(i),dam(i,mode))
149 ENDIF
150 ! -> Check mode 1 failure
151 IF (dam(i,mode) >= one) THEN
152 failmod = failmod + 1
153 fmode(i,mode) = 1
154 ENDIF
155c
156 ! Mode 2 failure - Tension YY
157 mode = 2
158 IF (signyy(i) > sigma_22t .AND. signyy(i) > zero .AND. dam(i,mode) < one) THEN
159 ! -> Linear stress softening
160 IF (ishap22t == 1) THEN
161 dam(i,mode) = dam(i,mode) + le(i)*depsyy(i)*sigma_22t/(two*g_22t)
162 dam(i,mode) = max(dam(i,mode),uvar(i,mode))
163 ! -> Exponential stress softening
164 ELSEIF (ishap22t == 2) THEN
165 ener(i,mode) = ener(i,mode) + signyy(i)*le(i)*depsyy(i)
166 ener(i,mode) = max(ener(i,mode),uvar(i,mode+6))
167 dam(i,mode) = (one - exp(-ener(i,mode)/g_22t))
168 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
169 ENDIF
170 dam(i,mode) = min(dam(i,mode),one)
171 ! Output damage variable
172 dfmax(i) = max(dfmax(i),dam(i,mode))
173 ENDIF
174 ! -> Check mode 2 failure
175 IF (dam(i,mode) >= one) THEN
176 failmod = failmod + 1
177 fmode(i,mode) = 1
178 ENDIF
179c
180 ! Mode 3 failure - Positive shear in plane XY
181 mode = 3
182 IF (signxy(i) > sigma_12t .AND. signxy(i) > zero .AND. dam(i,mode) < one) THEN
183 ! -> Linear stress softening
184 IF (ishap12t == 1) THEN
185 dam(i,mode) = dam(i,mode) + le(i)*depsxy(i)*sigma_12t/(four*g_12t)
186 dam(i,mode) = max(dam(i,mode),uvar(i,mode))
187 ! -> Exponential stress softening
188 ELSEIF (ishap12t == 2) THEN
189 ener(i,mode) = ener(i,mode) + signxy(i)*le(i)*half*depsxy(i)
190 ener(i,mode) = max(ener(i,mode),uvar(i,mode+6))
191 dam(i,mode) = (one - exp(-ener(i,mode)/g_12t))
192 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
193 ENDIF
194 dam(i,mode) = min(dam(i,mode),one)
195 ! Output damage variable
196 dfmax(i) = max(dfmax(i),dam(i,mode))
197 ENDIF
198 ! -> Check mode 3 failure
199 IF (dam(i,mode) >= one) THEN
200 failmod = failmod + 1
201 fmode(i,mode) = 1
202 ENDIF
203c
204 ! Mode 4 failure - Compression XX
205 mode = 4
206 IF (-signxx(i) > sigma_11c .AND. signxx(i) < zero .AND. dam(i,mode) < one) THEN
207 ! -> Linear stress softening
208 IF (ishap11c == 1) THEN
209 dam(i,mode) = dam(i,mode) + le(i)*abs(depsxx(i))*sigma_11c/(two*g_11c)
210 dam(i,mode) = max(dam(i,mode),uvar(i,mode))
211 ! -> Exponential stress softening
212 ELSEIF (ishap11c == 2) THEN
213 ener(i,mode) = ener(i,mode) + abs(signxx(i))*le(i)*abs(depsxx(i))
214 ener(i,mode) = max(ener(i,mode),uvar(i,mode+6))
215 dam(i,mode) = (one - exp(-ener(i,mode)/g_11c))
216 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
217 ENDIF
218 dam(i,mode) = min(dam(i,mode),one)
219 ! Output damage variable
220 dfmax(i) = max(dfmax(i),dam(i,mode))
221 ENDIF
222 ! -> Check mode 4 failure
223 IF (dam(i,mode) >= one) THEN
224 failmod = failmod + 1
225 fmode(i,mode) = 1
226 ENDIF
227c
228 ! Mode 5 failure - Compression YY
229 mode = 5
230 IF (-signyy(i) > sigma_22c .AND. signyy(i) < zero .AND. dam(i,mode) < one) THEN
231 ! -> Linear stress softening
232 IF (ishap22c == 1) THEN
233 dam(i,mode) = dam(i,mode) + le(i)*abs(depsyy(i))*sigma_22c/(two*g_22c)
234 dam(i,mode) = max(dam(i,mode),uvar(i,mode))
235 ! -> Exponential stress softening
236 ELSEIF (ishap22c == 2) THEN
237 ener(i,mode) = ener(i,mode) + abs(signyy(i))*le(i)*abs(depsyy(i))
238 ener(i,mode) = max(ener(i,mode),uvar(i,mode+6))
239 dam(i,mode) = (one - exp(-ener(i,mode)/g_22c))
240 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
241 ENDIF
242 dam(i,mode) = min(dam(i,mode),one)
243 ! Output damage variable
244 dfmax(i) = max(dfmax(i),dam(i,mode))
245 ENDIF
246 ! -> Check mode 5 failure
247 IF (dam(i,mode) >= one) THEN
248 failmod = failmod + 1
249 fmode(i,mode) = 1
250 ENDIF
251c
252 ! Mode 6 failure - Negative shear in plane XY
253 mode = 6
254 IF (-signxy(i) > sigma_12c .AND. signxy(i) < zero .AND. dam(i,mode) < one) THEN
255 ! -> Linear stress softening
256 IF (ishap12c == 1) THEN
257 dam(i,mode) = dam(i,mode) + le(i)*abs(depsxy(i))*sigma_12c/(four*g_12c)
258 dam(i,mode) = max(dam(i,mode),uvar(i,mode))
259 ! -> Exponential stress softening
260 ELSEIF (ishap12c == 2) THEN
261 ener(i,mode) = ener(i,mode) + abs(signxy(i))*le(i)*half*abs(depsxy(i))
262 ener(i,mode) = max(ener(i,mode),uvar(i,mode+6))
263 dam(i,mode) = (one - exp(-ener(i,mode)/g_12c))
264 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
265 ENDIF
266 dam(i,mode) = min(dam(i,mode),one)
267 ! Output damage variable
268 dfmax(i) = max(dfmax(i),dam(i,mode))
269 ENDIF
270 ! -> Check mode 6 failure
271 IF (dam(i,mode) >= one) THEN
272 failmod = failmod + 1
273 fmode(i,mode) = 1
274 ENDIF
275c
276 ! If at least one mode has failed
277 IF (failmod >= nmod) THEN
278 foff(i) = 0
279 tdel(i) = time
280 nindx = nindx + 1
281 indx(nindx) = i
282 ENDIF
283 ENDIF
284 ENDDO
285c
286 !====================================================================
287 ! - COMPUTE STRESS SOFTENING SCALE FACTORS
288 !====================================================================
289 DO i=1, nel
290 IF (foff(i) == 1) THEN
291 ! Direction 1
292 IF (signxx(i) >= zero) THEN
293 dmg_scale(i,1) = one - dam(i,1)
294 ELSE
295 dmg_scale(i,1) = one - dam(i,4)
296 ENDIF
297 ! Direction 2
298 IF (signyy(i) >= zero) THEN
299 dmg_scale(i,2) = one - dam(i,2)
300 ELSE
301 dmg_scale(i,2) = one - dam(i,5)
302 ENDIF
303 ! Plane 12
304 IF (signxy(i) >= zero) THEN
305 dmg_scale(i,3) = one - dam(i,3)
306 ELSE
307 dmg_scale(i,3) = one - dam(i,6)
308 ENDIF
309 ENDIF
310 ENDDO
311c
312 !====================================================================
313 ! - UPDATE USER VARIABLES
314 !====================================================================
315 DO j = 1,6
316 DO i = 1,nel
317 ! Damage variable per mode
318 uvar(i,j) = dam(i,j)
319 ! Dissipated energy per mode
320 uvar(i,j+6) = ener(i,j)
321 ENDDO
322 ENDDO
323c
324 !====================================================================
325 ! - PRINTOUT DATA ABOUT FAILED MODES
326 !====================================================================
327 IF (nindx > 0) THEN
328 DO j=1,nindx
329 i = indx(j)
330#include "lockon.inc"
331 IF (igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52) THEN
332 WRITE(iout, 1200) ngl(i),ipg,ply_id,ipt
333 WRITE(istdo,1200) ngl(i),ipg,ply_id,ipt
334 ELSEIF (igtyp == 1 .OR. igtyp == 9) THEN
335 WRITE(iout, 1000) ngl(i),ipg,ipt
336 WRITE(istdo,1000) ngl(i),ipg,ipt
337 ELSE
338 WRITE(iout, 1100) ngl(i),ipg,ilay,ipt
339 WRITE(istdo,1100) ngl(i),ipg,ilay,ipt
340 ENDIF
341 IF (fmode(i,1)==1) WRITE(iout, 2000) '1 - TRACTION XX'
342 IF (fmode(i,2)==1) WRITE(iout, 2000) '2 - TRACTION YY'
343 IF (fmode(i,3)==1) WRITE(iout, 2000) '3 - POSITIVE SHEAR XY'
344 IF (fmode(i,4)==1) WRITE(iout, 2000) '4 - COMPRESSION XX'
345 IF (fmode(i,5)==1) WRITE(iout, 2000) '5 - COMPRESSION YY'
346 IF (fmode(i,6)==1) WRITE(iout, 2000) '6 - NEGATIVE SHEAR XY'
347#include "lockoff.inc"
348 ENDDO
349 ENDIF
350c-----------------------------------------------------------------------
351 1000 FORMAT(1x,'FAILURE (ORTHENERG) OF SHELL ELEMENT ',i10,1x,',GAUSS PT',
352 . i2,1x,',INTEGRATION PT',i3)
353 1100 FORMAT(1x,'FAILURE (ORTHENERG) OF SHELL ELEMENT ',i10,1x,',GAUSS PT',
354 . i2,1x,',LAYER',i3,1x,',INTEGRATION PT',i3)
355 1200 FORMAT(1x,'FAILURE (ORTHENERG) OF SHELL ELEMENT ',i10,1x,',GAUSS PT',
356 . i2,1x,',PLY ID',i10,1x,',INTEGRATION PT',i3)
357 2000 FORMAT(1x,'MODE',1x,a)
358c-----------------------------------------------------------------------
359 END
fail_orthenerg_c
subroutine fail_orthenerg_c(nel, nuparam, nuvar, uparam, uvar, ngl, time, ipg, ilay, ipt, depsxx, depsyy, depsxy, dmg_flag, dmg_scale, aldt, foff, dfmax, tdel, signxx, signyy, signxy, igtyp, ply_id)
Definition fail_orthenerg_c.F:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
for
for(i8=*sizetab-1;i8 >=0;i8--)
Definition mumps_common.c:91
usermat_shell
subroutine usermat_shell(timers, elbuf_str, mat_elem, jft, jlt, nel, pm, for, mom, gstr, thk, eint, off, dir_a, dir_b, mat, area, exx, eyy, exy, exz, eyz, kxx, kyy, kxy, geo, thk_ly, pid, tf, npf, mtn, dt1c, dm, bufmat, ssp, rho, viscmx, ipla, iofc, indx, ngl, thkly, matly, zcfac, ng, shf, gs, sigy, thk0, epsd_pg, posly, igeo, ipm, failwave, fwave_el, ifailure, aldt, tempel, die, r11, r12, r13, r21, r22, r23, r31, r32, r33, table, ixfem, elcrkini, dir1_crk, dir2_crk, iparg, jhbe, ismstr, jthe, tensx, ir, is, nlay, npt, ixlay, ixel, ithk, f_def, ishplyxfem, itask, pm_stack, isubstack, stack, alpe, ply_exx, ply_eyy, ply_exy, ply_exz, ply_eyz, ply_f, varnl, nloc_dmg, nlay_max, laynpt_max, dt)
Definition usermat_shell.F:116