OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_visual_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_visual_s ../engine/source/materials/fail/visual/fail_visual_s.F
25!||--- called by ------------------------------------------------------
26!|| mmain ../engine/source/materials/mat_share/mmain.F90
27!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
28!|| mulaw8 ../engine/source/materials/mat_share/mulaw8.F90
29!|| usermat_solid ../engine/source/materials/mat_share/usermat_solid.F
30!||====================================================================
31 SUBROUTINE fail_visual_s(
32 . NEL ,NUPARAM ,NUVAR ,TIME ,TIMESTEP,UPARAM ,
33 . EPSXX ,EPSYY ,EPSZZ ,EPSXY ,EPSYZ ,EPSZX ,
34 . SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
35 . UVAR ,OFF ,NGL ,DFMAX ,ISMSTR )
36C--------------------------------------------------------------------
37C /FAIL/VISUAL - Visua failure criteria
38C--------------------------------------------------------------------
39C I m p l i c i t T y p e s
40C-----------------------------------------------
41#include "implicit_f.inc"
42C---------+---------+---+---+--------------------------------------------
43C VAR | SIZE |TYP| RW| DEFINITION
44C---------+--------+--+--+-------------------------------------------
45C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
46C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
47C NUVAR | 1 | I | R | NUMBER OF FAILURE ELEMENT VARIABLES
48C---------+--------+--+--+-------------------------------------------
49C TIME | 1 | F | R | CURRENT TIME
50C TIMESTEP| 1 | F | R | CURRENT TIME STEP
51C UPARAM | NUPARAM | F | R | USER FAILURE PARAMETER ARRAY
52C---------+--------+--+--+-------------------------------------------
53C EPSXX | NEL | F | R | STRAIN XX
54C EPSYY | NEL | F | R | STRAIN YY
55C ... | | | |
56C SIGNXX | NEL | F |R/W| NEW ELASTO PLASTIC STRESS XX
57C SIGNYY | NEL | F |R/W| NEW ELASTO PLASTIC STRESS YY
58C ... | | | |
59C---------+--------+--+--+-------------------------------------------
60C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
61C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
62C---------+--------+--+--+-------------------------------------------
63C I N P U T A r g u m e n t s
64C-----------------------------------------------
65#include "units_c.inc"
66#include "comlock.inc"
67C-----------------------------------------------
68 INTEGER NEL,NUPARAM,NUVAR,ISMSTR
69 INTEGER :: NGL(NEL)
70 my_real TIME,TIMESTEP,UPARAM(NUPARAM),
71 . EPSXX(NEL) ,EPSYY(NEL) ,EPSZZ(NEL) ,
72 . epsxy(nel) ,epsyz(nel) ,epszx(nel) ,dfmax(nel)
73C-----------------------------------------------
74C I N P U T O U T P U T A r g u m e n t s
75C-----------------------------------------------
76 my_real uvar(nel,nuvar)
77 my_real ,DIMENSION(NEL) :: off,
78 . signxx,signyy,signzz,signxy,signyz,signzx
79C-----------------------------------------------
80C L o c a l V a r i a b l e s
81C-----------------------------------------------
82 INTEGER I,J,NINDX,INDX(NEL),TYPE_MAX,F_FLAG,STRDEF,STRFLAG
83 my_real I1,I2,I3,E11,E22,E33,E_EQ,E_EQ1,E_EQ2,Q,R,R_INTER,PHI,
84 . C_MIN,C_MAX,EMA,DAMAGE
85 DOUBLE PRECISION :: A0(2),A1(2),A2(2),C0,C1,C2,C3,C4,C5,C6,C7,C8,C9,
86 . x1,x2,x3,y1,y2,y3,z1,z2,z3,f,ff,d,dd,d2,dp,e,g
87 DOUBLE PRECISION, PARAMETER :: PI8 = 0.3926990817d0
88 DOUBLE PRECISION, PARAMETER :: PI38 = 1.178097245d0
89 DOUBLE PRECISION, PARAMETER :: SPI8 = 0.3826834324d0
90 DOUBLE PRECISION, PARAMETER :: SPI38 = 0.9238795325d0
91C=======================================================================
92C USER VARIABLES
93
94c! User variable # 1, to store the previous damage value
95c! User variable # 2, to store the previous stress or strain value (for EMA filtering)
96c! User variable # 3-8, Storage values for the Butterworth filter
97C-----------------------------------------------
98C...
99 type_max = int(uparam(1))
100 c_min = uparam(2)
101 c_max = uparam(3)
102 ema = uparam(4)
103 ff = uparam(5)
104 f_flag = int(uparam(6))
105 strdef = nint(uparam(7))
106 nindx = 0
107 f = min(ff,zep4/max(em20,timestep))
108c----------------------------------------------
109c strain transformation flag following input definition
110c-------------------
111 strflag = 0
112 IF (strdef == 2) THEN ! failure defined as engineering strain
113 IF (ismstr == 10 .or. ismstr == 12) THEN
114 strflag = 1
115 ELSE IF (ismstr == 0 .or. ismstr == 2 .or. ismstr == 4) THEN
116 strflag = 2
117 END IF
118 ELSE IF (strdef == 3) THEN ! failure defined as true strain
119 IF (ismstr == 1 .or. ismstr == 3 .or. ismstr == 11) THEN
120 strflag = 3
121 ELSE IF (ismstr == 10 .or. ismstr == 12) THEN
122 strflag = 4
123 END IF
124 END IF
125
126c! ***********************
127c! *** NEW ... 26.06.2019
128c! ***********************
129 DO i=1,nel
130 IF (uvar(i,1) < one) THEN
131
132 IF (type_max == 1) THEN ! stress based
133c
134 i1 = signxx(i)+signyy(i)+signzz(i)
135 i2 = signxx(i)*signyy(i)+signyy(i)*signzz(i)+signzz(i)*signxx(i)-
136 . signxy(i)*signxy(i)-signzx(i)*signzx(i)-signyz(i)*signyz(i)
137 i3 = signxx(i)*signyy(i)*signzz(i)-signxx(i)*signyz(i)*signyz(i)-
138 . signyy(i)*signzx(i)*signzx(i)-signzz(i)*signxy(i)*signxy(i)+
139 . two*signxy(i)*signzx(i)*signyz(i)
140 q = (three*i2 - i1*i1)/9.0
141 q = min(q, zero)
142 r = (two*i1*i1*i1-9.0*i1*i2+27.0*i3)/54.0 ! (2*I3^3-9*I1*I2+27*I3)/54
143 r_inter = min(r/sqrt(max(em20,(-q**3))),one)
144 phi = acos(max(r_inter,-one))
145 e11 = two*sqrt(-q)*cos(phi/three)+third*i1
146 e22 = two*sqrt(-q)*cos((phi+two*pi)/three)+third*i1
147 e33 = two*sqrt(-q)*cos((phi+4.0*pi)/three)+third*i1
148
149 IF (e11 < e22) THEN
150 r_inter = e11
151 e11 = e22
152 e22 = r_inter
153 ENDIF
154 IF (e22 < e33)THEN
155 r_inter = e22
156 e22 = e33
157 e33 = r_inter
158 ENDIF
159 IF (e11 < e22)THEN
160 r_inter = e11
161 e11 = e22
162 e22 = r_inter
163 ENDIF
164c
165 ELSE ! TYPE_MAX - strain based criterion
166
167c! strains (deviatoric strain are in vector notation ==> gamma...==> division by 2 to get eps tensor components)
168
169 i1 = epsxx(i)+epsyy(i)+epszz(i)
170 i2 = epsxx(i)*epsyy(i)+epsyy(i)*epszz(i)+epszz(i)*epsxx(i)-
171 . half*epsxy(i)*half*epsxy(i)-half*epszx(i)*half*epszx(i)-
172 . half*epsyz(i)*half*epsyz(i)
173 i3 = epsxx(i)*epsyy(i)*epszz(i)-epsxx(i)*half*epsyz(i)*half*epsyz(i)-
174 . epsyy(i)*half*epszx(i)*half*epszx(i)-epszz(i)*half*epsxy(i)*half*epsxy(i)+
175 . two*half*epsxy(i)*half*epszx(i)*half*epsyz(i)
176 q = (three*i2 - i1*i1)/9.0
177 r = (two*i1*i1*i1-9.0*i1*i2+27.0*i3)/54.0 ! (2*I3^3-9*I1*I2+27*I3)/54
178
179 r_inter = min(r/sqrt(max(em20,(-q**3))),one)
180 phi = acos(max(r_inter,-one))
181
182 e11 = two*sqrt(-q)*cos(phi/three)+third*i1
183 e22 = two*sqrt(-q)*cos((phi+two*pi)/three)+third*i1
184 e33 = two*sqrt(-q)*cos((phi+4.0*pi)/three) +third*i1
185c
186 IF (strflag == 1) THEN
187 e11 = sqrt(e11 + one) - one
188 e22 = sqrt(e22 + one) - one
189 e33 = sqrt(e33 + one) - one
190 ELSE IF (strflag == 2) THEN
191 e11 = exp(e11) - one
192 e22 = exp(e22) - one
193 e33 = exp(e33) - one
194 ELSE IF (strflag == 3) THEN
195 e11 = log(e11 + one)
196 e22 = log(e22 + one)
197 e33 = log(e33 + one)
198 ELSE IF (strflag == 4) THEN
199 e11 = log(sqrt(e11 + one))
200 e22 = log(sqrt(e22 + one))
201 e33 = log(sqrt(e33 + one))
202 END IF
203c
204 IF (e11 < e22) THEN
205 r_inter = e11
206 e11 = e22
207 e22 = r_inter
208 ENDIF
209 IF (e22 < e33)THEN
210 r_inter = e22
211 e22 = e33
212 e33 = r_inter
213 ENDIF
214 IF (e11 < e22)THEN
215 r_inter = e11
216 e11 = e22
217 e22 = r_inter
218 ENDIF
219c
220 ENDIF ! TYPE_MAX
221c
222c! --- EMA or Butterworth filtering
223 IF (ema == one .and. ff /= zero .and. f_flag > 1) THEN
224C-----------------------------------------------
225C INITIALISATION OF THE FILTER-COEFFICIENTS
226C-----------------------------------------------
227C
228 d = tan(pi*f*timestep)
229 dd = d*d
230 d2 = two*d
231 dp = one + dd
232 e = d2*spi8
233 g = e + dp
234 g = one/g
235C
236 c0 = dd * g
237 c1 = two* c0
238 c2 = c0
239 c3 = two * g - c1
240 c4 = (e - dp) * g
241C
242 e = d2*spi38
243 g = e + dp
244 g = one/g
245C
246 c5 = dd * g
247 c6 = two * c5
248 c7 = c5
249 c8 = two * g - c6
250 c9 = (e - dp) * g
251
252C-----------------------------------------------
253C BUTTERWORTH FILTERING
254C-----------------------------------------------
255
256 a0(1) = uvar(i,3)*uvar(i,9)
257 a0(2) = uvar(i,4)*uvar(i,9)
258 a1(1) = uvar(i,5)*uvar(i,10)
259 a1(2) = uvar(i,6)*uvar(i,10)
260 a2(1) = uvar(i,7)*uvar(i,11)
261 a2(2) = uvar(i,8)*uvar(i,11)
262
263 x1 = a0(2)
264 x2 = a0(1)
265
266 x3 = e11
267 y1 = a1(2)
268 y2 = a1(1)
269 y3 = c0 * x3
270 y3 = y3 + c1 * x2
271 y3 = y3 + c2 * x1
272 y3 = y3 + c3 * y2
273 y3 = y3 + c4 * y1
274 z1 = a2(2)
275 z2 = a2(1)
276 z3 = c5 * y3
277 z3 = z3 + c6 * y2
278 z3 = z3 + c7 * y1
279 z3 = z3 + c8 * z2
280 z3 = z3 + c9 * z1
281C
282 a0(2) = x2
283 a0(1) = x3
284 a1(2) = y2
285 a1(1) = y3
286 a2(2) = z2
287 a2(1) = z3
288
289 IF ((x3 /= zero).AND.(x2 /= zero)) THEN
290 uvar(i,3) = a0(1)/x2
291 uvar(i,4) = a0(2)/x2
292 uvar(i,5) = a1(1)/(c0*x3)
293 uvar(i,6) = a1(2)/(c0*x3)
294 uvar(i,7) = a2(1)/(c0*y3)
295 uvar(i,8) = a2(2)/(c0*y3)
296 uvar(i,9) = x2
297 uvar(i,10) = c0*x3
298 uvar(i,11) = c0*y3
299 ELSE
300 uvar(i,3) = a0(1)
301 uvar(i,4) = a0(2)
302 uvar(i,5) = a1(1)
303 uvar(i,6) = a1(2)
304 uvar(i,7) = a2(1)
305 uvar(i,8) = a2(2)
306 uvar(i,9) = one
307 uvar(i,10) = one
308 uvar(i,11) = one
309 ENDIF
310
311 e11 = a2(1)
312
313 ELSE
314c!
315c! What it should be is like this:
316c!
317c! Value = USER_INPUT * 2 * Pi * DT1
318c! Alpha = Value / (Value + 1)
319c! Actual_filtered_stress = Alpha * actual_Stress + (1-Alpha) * previous_filtered_stress
320c!
321c! ALPHA = EMA / ( EMA + ONE)
322c! E11 = ALPHA * E11 + ( ONE - ALPHA ) * UVAR(I,2)
323c! UVAR(I,2) = E11
324
325c! EMA = 0 ==> 1 ==> no filtering ; EMA = 1e+20 extreme filtering
326c! E11 = E11*(TWO/(ONE+EMA)) + UVAR(I,2) * (ONE-(TWO/(ONE+EMA)))
327 e11 = ema * e11 + ( one - ema ) * uvar(i,2)
328 uvar(i,2) = e11
329 ENDIF
330c!
331 damage = max(zero , min(one ,(e11-c_min)/max(em6,(c_max-c_min)) ))
332 uvar(i,1) = max(uvar(i,1),damage)
333 dfmax(i) = uvar(i,1)
334
335 IF (uvar(i,1) >= one) THEN
336 nindx = nindx+1
337 indx(nindx) = i
338 ENDIF
339
340 ENDIF ! UVAR(I,1) < ONE
341
342 ENDDO ! NEL
343c---------------------------------------------
344 DO j=1,nindx
345 i = indx(j)
346#include "lockon.inc"
347 WRITE(iout, 1000) ngl(i),time
348 WRITE(istdo,1100) ngl(i),time
349#include "lockoff.inc"
350
351 ENDDO
352c---------------------------------------------
353 1000 FORMAT(1x,'SOLID ELEMENT NUMBER (VISUAL) el#',i10,
354 . ' LIMIT REACHED AT TIME :',1pe12.4)
355 1100 FORMAT(1x,'SOLID ELEMENT NUMBER (VISUAL) el#',i10,
356 . ' LIMIT REACHED AT TIME :',1pe12.4)
357c-----------
358 RETURN
359 END
fail_visual_s
subroutine fail_visual_s(nel, nuparam, nuvar, time, timestep, uparam, epsxx, epsyy, epszz, epsxy, epsyz, epszx, signxx, signyy, signzz, signxy, signyz, signzx, uvar, off, ngl, dfmax, ismstr)
Definition fail_visual_s.F:36
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21