OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_orthenerg_s.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_s ../engine/source/materials/fail/orthenerg/fail_orthenerg_s.F
25!||--- called by ------------------------------------------------------
26!|| mmain ../engine/source/materials/mat_share/mmain.F90
27!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
28!|| usermat_solid ../engine/source/materials/mat_share/usermat_solid.F
29!||====================================================================
30 SUBROUTINE fail_orthenerg_s(
31 1 NEL ,NUPARAM ,NUVAR ,UPARAM ,UVAR ,NGL ,
32 2 NPG ,IPG ,ILAY ,OFF ,LOFF ,NOFF ,
33 3 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
34 4 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
35 5 TIME ,TDEL ,DFMAX ,ALDT ,DMG_SCALE)
36C-----------------------------------------------
37C I m p l i c i t T y p e s
38C-----------------------------------------------
39#include "implicit_f.inc"
40C-----------------------------------------------
41C G l o b a l P a r a m e t e r s
42C-----------------------------------------------
43#include "units_c.inc"
44#include "comlock.inc"
45C-----------------------------------------------
46C I N P U T A r g u m e n t s
47C-----------------------------------------------
48 INTEGER, INTENT(IN) :: NEL,NUPARAM,NUVAR,IPG,NPG,ILAY
49 INTEGER, DIMENSION(NEL), INTENT(IN) :: NGL
50 my_real, INTENT(IN) :: TIME
51 my_real, DIMENSION(NEL), INTENT(IN) :: DEPSXX,DEPSYY,
52 . DEPSZZ,DEPSXY,DEPSYZ,DEPSZX,SIGNXX,SIGNYY,SIGNZZ,
53 . signxy,signyz,signzx,aldt
54 my_real, DIMENSION(NUPARAM), INTENT(IN) :: uparam
55C-----------------------------------------------
56C I N P U T O U T P U T A r g u m e n t s
57C-----------------------------------------------
58 INTEGER, DIMENSION(NEL), INTENT(INOUT) :: NOFF
59 my_real, DIMENSION(NEL), INTENT(INOUT) :: dfmax,tdel,off,loff
60 my_real, DIMENSION(NEL,6), INTENT(INOUT) :: 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,NINDX2,FAILMOD,MODE,ISHAP11T,ISHAP11C,
66 . ishap22t,ishap22c,ishap33t,ishap33c,ishap12t,ishap12c,
67 . ishap23t,ishap23c,ishap31t,ishap31c,failip,nmod
68 INTEGER ,DIMENSION(NEL) :: INDX, INDX2
69 INTEGER ,DIMENSION(NEL,12) :: FMODE
70 my_real dam(nel,12),ener(nel,12),le(nel)
71 my_real sigma_11t,sigma_11c,sigma_22t,sigma_22c,sigma_33t,
72 . sigma_33c,sigma_12t,sigma_12c,sigma_23t,sigma_23c,
73 . sigma_31t,sigma_31c,g_11t,g_11c,g_22t,g_22c,g_33t,g_33c,
74 . g_12t,g_12c,g_23t,g_23c,g_31t,g_31c
75C
76 !=======================================================================
77 ! - INITIALISATION OF COMPUTATION ON TIME STEP
78 !=======================================================================
79 ! Recovering failure criterion parameters
80 failip = nint(uparam(2)) ! -> Percentage of failed thickness/ratio
81 failip = min(failip,npg)
82 sigma_11t = uparam(3) ! -> Critical stress for tension in direction 11
83 sigma_11c = uparam(4) ! -> Critical stress for compression in direction 11
84 sigma_22t = uparam(5) ! -> Critical stress for tension in direction 22
85 sigma_22c = uparam(6) ! -> Critical stress for compression in direction 22
86 sigma_33t = uparam(7) ! -> Critical stress for tension in direction 33
87 sigma_33c = uparam(8) ! -> Critical stress for compression in direction 33
88 sigma_12t = uparam(9) ! -> Critical stress for positive shear in plane 12
89 sigma_12c = uparam(10) ! -> Critical stress for negative shear in plane 12
90 sigma_23t = uparam(11) ! -> Critical stress for positive shear in plane 23
91 sigma_23c = uparam(12) ! -> Critical stress for negative shear in plane 23
92 sigma_31t = uparam(13) ! -> Critical stress for positive shear in plane 31
93 sigma_31c = uparam(14) ! -> Critical stress for negative shear in plane 31
94 g_11t = uparam(15) ! -> Fracture energy for tension in direction 11
95 g_11c = uparam(16) ! -> Fracture energy for compression in direction 11
96 g_22t = uparam(17) ! -> Fracture energy for tension in direction 22
97 g_22c = uparam(18) ! -> Fracture energy for compression in direction 22
98 g_33t = uparam(19) ! -> Fracture energy for tension in direction 33
99 g_33c = uparam(20) ! -> Fracture energy for compression in direction 33
100 g_12t = uparam(21) ! -> Fracture energy for positive shear in plane 12
101 g_12c = uparam(22) ! -> Fracture energy for negative shear in plane 12
102 g_23t = uparam(23) ! -> Fracture energy for positive shear in plane 23
103 g_23c = uparam(24) ! -> Fracture energy for negative shear in plane 23
104 g_31t = uparam(25) ! -> Fracture energy for positive shear in plane 31
105 g_31c = uparam(26) ! -> Fracture energy for negative shear in plane 31
106 ishap11t = nint(uparam(27)) ! -> Softening shape flag for tension in direction 11
107 ishap11c = nint(uparam(28)) ! -> Softening shape flag for compression in direction 11
108 ishap22t = nint(uparam(29)) ! -> Softening shape flag for tension in direction 22
109 ishap22c = nint(uparam(30)) ! -> Softening shape flag for compression in direction 22
110 ishap33t = nint(uparam(31)) ! -> Softening shape flag for tension in direction 33
111 ishap33c = nint(uparam(32)) ! -> Softening shape flag for compression in direction 33
112 ishap12t = nint(uparam(33)) ! -> Softening shape flag for positive shear in plane 12
113 ishap12c = nint(uparam(34)) ! -> Softening shape flag for negative shear in plane 12
114 ishap23t = nint(uparam(35)) ! -> Softening shape flag for positive shear in plane 23
115 ishap23c = nint(uparam(36)) ! -> Softening shape flag for negative shear in plane 23
116 ishap31t = nint(uparam(37)) ! -> Softening shape flag for positive shear in plane 31
117 ishap31c = nint(uparam(38)) ! -> Softening shape flag for negative shear in plane 31
118 nmod = nint(uparam(39)) ! -> Number of failed modes prior to integration point failure
119C
120 ! Save initial element length
121 IF (uvar(1,25) == zero) THEN
122 uvar(1:nel,25) = aldt(1:nel)
123 ENDIF
124 le(1:nel) = uvar(1:nel,25)
125c
126 ! Recover user variable value
127 DO j = 1,12
128 DO i = 1,nel
129 ! Damage variable per mode
130 dam(i,j) = uvar(i,j)
131 ! Dissipated energy per mode
132 ener(i,j) = uvar(i,j+12)
133 ENDDO
134 ENDDO
135c
136 !====================================================================
137 ! - COMPUTATION OF THE DAMAGE VARIABLE EVOLUTION
138 !====================================================================
139 ! Initialization of element failure index
140 nindx = 0
141 nindx2 = 0
142 indx(1:nel) = 0
143 indx2(1:nel) = 0
144c
145 ! Loop over the elements
146 DO i = 1, nel
147c
148 IF ((loff(i) == one).AND.(off(i) == one)) THEN
149c
150 ! Initialization of failed modes index
151 fmode(i,1:12) = 0
152 failmod = 0
153c
154 ! Mode 1 failure - Tension XX
155 mode = 1
156 IF (signxx(i) > sigma_11t .AND. signxx(i) > zero .AND. dam(i,mode) < one) THEN
157 ! -> Linear stress softening
158 IF (ishap11t == 1) THEN
159 dam(i,mode) = dam(i,mode) + le(i)*depsxx(i)*sigma_11t/(two*g_11t)
160 dam(i,mode) = max(dam(i,mode),uvar(i,mode))
161 ! -> Exponential stress softening
162 ELSEIF (ishap11t == 2) THEN
163 ener(i,mode) = ener(i,mode) + signxx(i)*le(i)*depsxx(i)
164 ener(i,mode) = max(ener(i,mode),uvar(i,mode+12))
165 dam(i,mode) = (one - exp(-ener(i,mode)/g_11t))
166 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
167 ENDIF
168 dam(i,mode) = min(dam(i,mode),one)
169 ! Output damage variable
170 dfmax(i) = max(dfmax(i),dam(i,mode))
171 ENDIF
172 ! -> Check mode 1 failure
173 IF (dam(i,mode) >= one) THEN
174 failmod = failmod + 1
175 fmode(i,mode) = 1
176 ENDIF
177c
178 ! Mode 2 failure - Tension YY
179 mode = 2
180 IF (signyy(i) > sigma_22t .AND. signyy(i) > zero .AND. dam(i,mode) < one) THEN
181 ! -> Linear stress softening
182 IF (ishap22t == 1) THEN
183 dam(i,mode) = dam(i,mode) + le(i)*depsyy(i)*sigma_22t/(two*g_22t)
184 dam(i,mode) = max(dam(i,mode),uvar(i,mode))
185 ! -> Exponential stress softening
186 ELSEIF (ishap22t == 2) THEN
187 ener(i,mode) = ener(i,mode) + signyy(i)*le(i)*depsyy(i)
188 ener(i,mode) = max(ener(i,mode),uvar(i,mode+12))
189 dam(i,mode) = (one - exp(-ener(i,mode)/g_22t))
190 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
191 ENDIF
192 dam(i,mode) = min(dam(i,mode),one)
193 ! Output damage variable
194 dfmax(i) = max(dfmax(i),dam(i,mode))
195 ENDIF
196 ! -> Check mode 2 failure
197 IF (dam(i,mode) >= one) THEN
198 failmod = failmod + 1
199 fmode(i,mode) = 1
200 ENDIF
201c
202 ! Mode 3 failure - Positive shear in plane XY
203 mode = 3
204 IF (signxy(i) > sigma_12t .AND. signxy(i) > zero .AND. dam(i,mode) < one) THEN
205 ! -> Linear stress softening
206 IF (ishap12t == 1) THEN
207 dam(i,mode) = dam(i,mode) + le(i)*depsxy(i)*sigma_12t/(four*g_12t)
208 dam(i,mode) = max(dam(i,mode),uvar(i,mode))
209 ! -> Exponential stress softening
210 ELSEIF (ishap12t == 2) THEN
211 ener(i,mode) = ener(i,mode) + signxy(i)*le(i)*half*depsxy(i)
212 ener(i,mode) = max(ener(i,mode),uvar(i,mode+12))
213 dam(i,mode) = (one - exp(-ener(i,mode)/g_12t))
214 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
215 ENDIF
216 dam(i,mode) = min(dam(i,mode),one)
217 ! Output damage variable
218 dfmax(i) = max(dfmax(i),dam(i,mode))
219 ENDIF
220 ! -> Check mode 3 failure
221 IF (dam(i,mode) >= one) THEN
222 failmod = failmod + 1
223 fmode(i,mode) = 1
224 ENDIF
225c
226 ! Mode 4 failure - Compression XX
227 mode = 4
228 IF (-signxx(i) > sigma_11c .AND. signxx(i) < zero .AND. dam(i,mode) < one) THEN
229 ! -> Linear stress softening
230 IF (ishap11c == 1) THEN
231 dam(i,mode) = dam(i,mode) + le(i)*abs(depsxx(i))*sigma_11c/(two*g_11c)
232 dam(i,mode) = max(dam(i,mode),uvar(i,mode))
233 ! -> Exponential stress softening
234 ELSEIF (ishap11c == 2) THEN
235 ener(i,mode) = ener(i,mode) + abs(signxx(i))*le(i)*abs(depsxx(i))
236 ener(i,mode) = max(ener(i,mode),uvar(i,mode+12))
237 dam(i,mode) = (one - exp(-ener(i,mode)/g_11c))
238 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
239 ENDIF
240 dam(i,mode) = min(dam(i,mode),one)
241 ! Output damage variable
242 dfmax(i) = max(dfmax(i),dam(i,mode))
243 ENDIF
244 ! -> Check mode 4 failure
245 IF (dam(i,mode) >= one) THEN
246 failmod = failmod + 1
247 fmode(i,mode) = 1
248 ENDIF
249c
250 ! Mode 5 failure - Compression YY
251 mode = 5
252 IF (-signyy(i) > sigma_22c .AND. signyy(i) < zero .AND. dam(i,mode) < one) THEN
253 ! -> Linear stress softening
254 IF (ishap22c == 1) THEN
255 dam(i,mode) = dam(i,mode) + le(i)*abs(depsyy(i))*sigma_22c/(two*g_22c)
256 dam(i,mode) = max(dam(i,mode),uvar(i,mode))
257 ! -> Exponential stress softening
258 ELSEIF (ishap22c == 2) THEN
259 ener(i,mode) = ener(i,mode) + abs(signyy(i))*le(i)*abs(depsyy(i))
260 ener(i,mode) = max(ener(i,mode),uvar(i,mode+12))
261 dam(i,mode) = (one - exp(-ener(i,mode)/g_22c))
262 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
263 ENDIF
264 dam(i,mode) = min(dam(i,mode),one)
265 ! Output damage variable
266 dfmax(i) = max(dfmax(i),dam(i,mode))
267 ENDIF
268 ! -> Check mode 5 failure
269 IF (dam(i,mode) >= one) THEN
270 failmod = failmod + 1
271 fmode(i,mode) = 1
272 ENDIF
273c
274 ! Mode 6 failure - Negative shear in plane XY
275 mode = 6
276 IF (-signxy(i) > sigma_12c .AND. signxy(i) < zero .AND. dam(i,mode) < one) THEN
277 ! -> Linear stress softening
278 IF (ishap12c == 1) THEN
279 dam(i,mode) = dam(i,mode) + le(i)*abs(depsxy(i))*sigma_12c/(four*g_12c)
280 dam(i,mode) = max(dam(i,mode),uvar(i,mode))
281 ! -> Exponential stress softening
282 ELSEIF (ishap12c == 2) THEN
283 ener(i,mode) = ener(i,mode) + abs(signxy(i))*le(i)*half*abs(depsxy(i))
284 ener(i,mode) = max(ener(i,mode),uvar(i,mode+12))
285 dam(i,mode) = (one - exp(-ener(i,mode)/g_12c))
286 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
287 ENDIF
288 dam(i,mode) = min(dam(i,mode),one)
289 ! Output damage variable
290 dfmax(i) = max(dfmax(i),dam(i,mode))
291 ENDIF
292 ! -> Check mode 6 failure
293 IF (dam(i,mode) >= one) THEN
294 failmod = failmod + 1
295 fmode(i,mode) = 1
296 ENDIF
297c
298 ! Mode 7 failure - Tension ZZ
299 mode = 7
300 IF (signzz(i) > sigma_33t .AND. signzz(i) > zero .AND. dam(i,mode) < one) THEN
301 ! -> Linear stress softening
302 IF (ishap33t == 1) THEN
303 dam(i,mode) = dam(i,mode) + le(i)*depszz(i)*sigma_33t/(two*g_33t)
304 dam(i,mode) = max(dam(i,mode),uvar(i,mode))
305 ! -> Exponential stress softening
306 ELSEIF (ishap33t == 2) THEN
307 ener(i,mode) = ener(i,mode) + signzz(i)*le(i)*depszz(i)
308 ener(i,mode) = max(ener(i,mode),uvar(i,mode+12))
309 dam(i,mode) = (one - exp(-ener(i,mode)/g_33t))
310 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
311 ENDIF
312 dam(i,mode) = min(dam(i,mode),one)
313 ! Output damage variable
314 dfmax(i) = max(dfmax(i),dam(i,mode))
315 ENDIF
316 ! -> Check mode 7 failure
317 IF (dam(i,mode) >= one) THEN
318 failmod = failmod + 1
319 fmode(i,mode) = 1
320 ENDIF
321c
322 ! Mode 8 failure - Compression ZZ
323 mode = 8
324 IF (-signzz(i) > sigma_33c .AND. signzz(i) < zero .AND. dam(i,mode) < one) THEN
325 ! -> Linear stress softening
326 IF (ishap33c == 1) THEN
327 dam(i,mode) = dam(i,mode) + le(i)*abs(depszz(i))*sigma_33c/(two*g_33c)
328 dam(i,mode) = max(dam(i,mode),uvar(i,mode))
329 ! -> Exponential stress softening
330 ELSEIF (ishap33c == 2) THEN
331 ener(i,mode) = ener(i,mode) + abs(signzz(i))*le(i)*abs(depszz(i))
332 ener(i,mode) = max(ener(i,mode),uvar(i,mode+12))
333 dam(i,mode) = (one - exp(-ener(i,mode)/g_33c))
334 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
335 ENDIF
336 dam(i,mode) = min(dam(i,mode),one)
337 ! Output damage variable
338 dfmax(i) = max(dfmax(i),dam(i,mode))
339 ENDIF
340 ! -> Check mode 8 failure
341 IF (dam(i,mode) >= one) THEN
342 failmod = failmod + 1
343 fmode(i,mode) = 1
344 ENDIF
345c
346 ! Mode 9 failure - Positive shear in plane YZ
347 mode = 9
348 IF (signyz(i) > sigma_23t .AND. signyz(i) > zero .AND. dam(i,mode) < one) THEN
349 ! -> Linear stress softening
350 IF (ishap23t == 1) THEN
351 dam(i,mode) = dam(i,mode) + le(i)*depsyz(i)*sigma_23t/(four*g_23t)
352 dam(i,mode) = max(dam(i,mode),uvar(i,mode))
353 ! -> Exponential stress softening
354 ELSEIF (ishap23t == 2) THEN
355 ener(i,mode) = ener(i,mode) + signyz(i)*le(i)*half*depsyz(i)
356 ener(i,mode) = max(ener(i,mode),uvar(i,mode+12))
357 dam(i,mode) = (one - exp(-ener(i,mode)/g_23t))
358 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
359 ENDIF
360 dam(i,mode) = min(dam(i,mode),one)
361 ! output damage variable
362 dfmax(i) = max(dfmax(i),dam(i,mode))
363 ENDIF
364 ! -> Check mode 9 failure
365 IF (dam(i,mode) >= one) THEN
366 failmod = failmod + 1
367 fmode(i,mode) = 1
368 ENDIF
369c
370 ! Mode 10 failure - Negative shear in plane YZ
371 mode = 10
372 IF (-signyz(i) > sigma_23c .AND. signyz(i) < zero .AND. dam(i,mode) < one) THEN
373 ! -> Linear stress softening
374 IF (ishap23c == 1) THEN
375 dam(i,mode) = dam(i,mode) + le(i)*abs(depsyz(i))*sigma_23c/(four*g_23c)
376 dam(i,mode) = max(dam(i,mode),uvar(i,mode))
377 ! -> Exponential stress softening
378 ELSEIF (ishap23c == 2) THEN
379 ener(i,mode) = ener(i,mode) + abs(signyz(i))*le(i)*half*abs(depsyz(i))
380 ener(i,mode) = max(ener(i,mode),uvar(i,mode+12))
381 dam(i,mode) = (one - exp(-ener(i,mode)/g_23c))
382 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
383 ENDIF
384 dam(i,mode) = min(dam(i,mode),one)
385 ! Output damage variable
386 dfmax(i) = max(dfmax(i),dam(i,mode))
387 ENDIF
388 ! -> Check mode 10 failure
389 IF (dam(i,mode) >= one) THEN
390 failmod = failmod + 1
391 fmode(i,mode) = 1
392 ENDIF
393c
394 ! Mode 11 failure - Positive shear in plane ZX
395 mode = 11
396 IF (signzx(i) > sigma_31t .AND. signzx(i) > zero .AND. dam(i,mode) < one) THEN
397 ! -> Linear stress softening
398 IF (ishap31t == 1) THEN
399 dam(i,mode) = dam(i,mode) + le(i)*depszx(i)*sigma_31t/(four*g_31t)
400 dam(i,mode) = max(dam(i,mode),uvar(i,mode))
401 ! -> Exponential stress softening
402 ELSEIF (ishap31t == 2) THEN
403 ener(i,mode) = ener(i,mode) + signzx(i)*le(i)*half*depszx(i)
404 ener(i,mode) = max(ener(i,mode),uvar(i,mode+12))
405 dam(i,mode) = (one - exp(-ener(i,mode)/g_31t))
406 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
407 ENDIF
408 dam(i,mode) = min(dam(i,mode),one)
409 ! Output damage variable
410 dfmax(i) = max(dfmax(i),dam(i,mode))
411 ENDIF
412 ! -> Check mode 11 failure
413 IF (dam(i,mode) >= one) THEN
414 failmod = failmod + 1
415 fmode(i,mode) = 1
416 ENDIF
417c
418 ! Mode 12 failure - Negative shear in plane YZ
419 mode = 12
420 IF (-signzx(i) > sigma_31c .AND. signzx(i) < zero .AND. dam(i,mode) < one) THEN
421 ! -> Linear stress softening
422 IF (ishap31c == 1) THEN
423 dam(i,mode) = dam(i,mode) + le(i)*abs(depszx(i))*sigma_31c/(four*g_31c)
424 dam(i,mode) = max(dam(i,mode),uvar(i,mode))
425 ! -> Exponential stress softening
426 ELSEIF (ishap31c == 2) THEN
427 ener(i,mode) = ener(i,mode) + abs(signzx(i))*le(i)*half*abs(depszx(i))
428 ener(i,mode) = max(ener(i,mode),uvar(i,mode+12))
429 dam(i,mode) = (one - exp(-ener(i,mode)/g_31c))
430 IF (dam(i,mode) > 0.999d0) dam(i,mode) = one
431 ENDIF
432 dam(i,mode) = min(dam(i,mode),one)
433 ! Output damage variable
434 dfmax(i) = max(dfmax(i),dam(i,mode))
435 ENDIF
436 ! -> Check mode 12 failure
437 IF (dam(i,mode) >= one) THEN
438 failmod = failmod + 1
439 fmode(i,mode) = 1
440 ENDIF
441c
442 ! If at least one mode has failed
443 IF (failmod >= nmod) THEN
444 loff(i) = zero
445 nindx = nindx + 1
446 indx(nindx) = i
447 noff(i) = noff(i) + 1
448 IF (noff(i) >= failip) THEN
449 off(i) = zero
450 tdel(i) = time
451 nindx2 = nindx2 + 1
452 indx2(nindx2) = i
453 ENDIF
454 ENDIF
455 ENDIF
456 ENDDO
457c
458 !====================================================================
459 ! - COMPUTE STRESS SOFTENING SCALE FACTORS
460 !====================================================================
461 DO i = 1, nel
462 IF ((loff(i) == one) .AND. (off(i) == one)) THEN
463 ! Direction 1
464 IF (signxx(i) >= zero) THEN
465 dmg_scale(i,1) = one - dam(i,1)
466 ELSE
467 dmg_scale(i,1) = one - dam(i,4)
468 ENDIF
469 ! Direction 2
470 IF (signyy(i) >= zero) THEN
471 dmg_scale(i,2) = one - dam(i,2)
472 ELSE
473 dmg_scale(i,2) = one - dam(i,5)
474 ENDIF
475 ! Direction 3
476 IF (signzz(i) >= zero) THEN
477 dmg_scale(i,3) = one - dam(i,7)
478 ELSE
479 dmg_scale(i,3) = one - dam(i,8)
480 ENDIF
481 ! Plane 12
482 IF (signxy(i) >= zero) THEN
483 dmg_scale(i,4) = one - dam(i,3)
484 ELSE
485 dmg_scale(i,4) = one - dam(i,6)
486 ENDIF
487 ! Plane 23
488 IF (signyz(i) >= zero) THEN
489 dmg_scale(i,5) = one - dam(i,9)
490 ELSE
491 dmg_scale(i,5) = one - dam(i,10)
492 ENDIF
493 ! Plane 31
494 IF (signzx(i) >= zero) THEN
495 dmg_scale(i,6) = one - dam(i,11)
496 ELSE
497 dmg_scale(i,6) = one - dam(i,12)
498 ENDIF
499 ENDIF
500 ENDDO
501c
502 !====================================================================
503 ! - UPDATE USER VARIABLES
504 !====================================================================
505 DO j = 1,12
506 DO i = 1,nel
507 ! Damage variable per mode
508 uvar(i,j) = dam(i,j)
509 ! Dissipated energy per mode
510 uvar(i,j+12) = ener(i,j)
511 ENDDO
512 ENDDO
513c
514 !====================================================================
515 ! - PRINTOUT DATA ABOUT FAILED MODES
516 !====================================================================
517 IF (nindx > 0) THEN
518 DO j=1,nindx
519 i = indx(j)
520#include "lockon.inc"
521 WRITE(iout, 1000) ngl(i),ipg,ilay
522 WRITE(istdo,1000) ngl(i),ipg,ilay
523 IF (fmode(i,1) == 1) WRITE(iout,2000) '1 - TRACTION XX'
524 IF (fmode(i,2) == 1) WRITE(iout,2000) '2 - TRACTION YY'
525 IF (fmode(i,3) == 1) WRITE(iout,2000) '3 - POSITIVE SHEAR XY'
526 IF (fmode(i,4) == 1) WRITE(iout,2000) '4 - COMPRESSION XX'
527 IF (fmode(i,5) == 1) WRITE(iout,2000) '5 - COMPRESSION YY'
528 IF (fmode(i,6) == 1) WRITE(iout,2000) '6 - NEGATIVE SHEAR XY'
529 IF (fmode(i,7) == 1) WRITE(iout,2000) '7 - TRACTION ZZ'
530 IF (fmode(i,8) == 1) WRITE(iout,2000) '8 - COMPRESSION ZZ'
531 IF (fmode(i,9) == 1) WRITE(iout,2000) '9 - POSITIVE SHEAR YZ'
532 IF (fmode(i,10) == 1) WRITE(iout,2000) '10 - NEGATIVE SHEAR YZ'
533 IF (fmode(i,11) == 1) WRITE(iout,2000) '11 - POSITIVE SHEAR ZX'
534 IF (fmode(i,12) == 1) WRITE(iout,2000) '12 - NEGATIVE SHEAR ZX'
535#include "lockoff.inc"
536 ENDDO
537 ENDIF
538c
539 IF (nindx2 > 0) THEN
540 DO j=1,nindx2
541 i = indx2(j)
542#include "lockon.inc"
543 WRITE(iout, 3000) ngl(i),time
544 WRITE(istdo,3000) ngl(i),time
545#include "lockoff.inc"
546 ENDDO
547 ENDIF
548c-----------------------------------------------------------------------
549 1000 FORMAT(1x,'FAILURE (ORTHENERG) OF SOLID ELEMENT ',i10,1x,',GAUSS PT',
550 . i5,1x,',LAYER',i3)
551 2000 FORMAT(1x,'MODE',1x,a)
552 3000 FORMAT(1x,'-- RUPTURE OF SOLID ELEMENT : ',i10,1x,'AT TIME : ',1pe12.4)
553c-----------------------------------------------------------------------
554 END
fail_orthenerg_s
subroutine fail_orthenerg_s(nel, nuparam, nuvar, uparam, uvar, ngl, npg, ipg, ilay, off, loff, noff, depsxx, depsyy, depszz, depsxy, depsyz, depszx, signxx, signyy, signzz, signxy, signyz, signzx, time, tdel, dfmax, aldt, dmg_scale)
Definition fail_orthenerg_s.F:36
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21