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

Go to the source code of this file.

Functions/Subroutines

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)

Function/Subroutine Documentation

◆ fail_visual_s()

subroutine fail_visual_s ( integer nel,
integer nuparam,
integer nuvar,
time,
timestep,
uparam,
epsxx,
epsyy,
epszz,
epsxy,
epsyz,
epszx,
dimension(nel) signxx,
dimension(nel) signyy,
dimension(nel) signzz,
dimension(nel) signxy,
dimension(nel) signyz,
dimension(nel) signzx,
uvar,
dimension(nel) off,
integer, dimension(nel) ngl,
dfmax,
integer ismstr )

Definition at line 31 of file fail_visual_s.F.

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
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21