OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_biquad_c.F File Reference
#include "implicit_f.inc"
#include "units_c.inc"
#include "comlock.inc"
#include "mvsiz_p.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine fail_biquad_c (nel, nuvar, time, uparam, ngl, ipt, nptot, signxx, signyy, signxy, signyz, signzx, dpla, uvar, uel1, off, offl, dfmax, tdel, nfunc, ifunc, npf, tf, aldt, foff, ipg, dmg_flag, dmgscl)

Function/Subroutine Documentation

◆ fail_biquad_c()

subroutine fail_biquad_c ( integer nel,
integer nuvar,
time,
uparam,
integer, dimension(nel) ngl,
integer ipt,
integer nptot,
signxx,
signyy,
signxy,
signyz,
signzx,
dpla,
uvar,
uel1,
off,
offl,
dfmax,
tdel,
integer nfunc,
integer, dimension(nfunc) ifunc,
integer, dimension(*) npf,
tf,
aldt,
integer, dimension(nel), intent(inout) foff,
integer ipg,
integer, intent(inout) dmg_flag,
dimension(nel), intent(inout) dmgscl )

Definition at line 31 of file fail_biquad_c.F.

39C-----------------------------------------------
40C Biquadratic failure model
41C-----------------------------------------------
42C I m p l i c i t T y p e s
43C-----------------------------------------------
44#include "implicit_f.inc"
45C-----------------------------------------------
46C G l o b a l P a r a m e t e r s
47C-----------------------------------------------
48#include "units_c.inc"
49#include "comlock.inc"
50#include "mvsiz_p.inc"
51C---------+---------+---+---+--------------------------------------------
52C VAR | SIZE |TYP| RW| DEFINITION
53C---------+---------+---+---+--------------------------------------------
54C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
55C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
56C---------+---------+---+---+--------------------------------------------
57C NPT | 1 | I | R | NUMBER OF LAYERS OR INTEGRATION POINTS
58C IPT | 1 | I | R | LAYER OR INTEGRATION POINT NUMBER
59C ISHELL | * | I | R | GEOMETRICAL FLAGS
60C NGL | NEL | I | R | ELEMENT NUMBER
61C SHF | NEL | F | R | SHEAR FACTOR
62C---------+---------+---+---+--------------------------------------------
63C TIME | 1 | F | R | CURRENT TIME
64C---------+---------+---+---+--------------------------------------------
65C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
66C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
67C ... | | | |
68C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
69C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
70C---------+---------+---+---+--------------------------------------------
71C PLA | NEL | F |R/W| PLASTIC STRAIN
72C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
73C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
74C FOFF | NEL | I |R/W| DELETED INTEGRATION POINT FLAG (=1 ON, =0 OFF)
75C---------+---------+---+---+--------------------------------------------
76C IPT CURRENT INTEGRATION POINT IN ALL LAYERS
77C IPTT CURRENT INTEGRATION POINT IN THE LAYER (FOR OUTPUT ONLY)
78C NPTOT NUMBER OF INTEGRATION POINTS IN ALL LAYERS (TOTAL)
79C NOFF NUMBER OF FAILED INTEGRATION POINTS (TOTAL)
80C IGTYP PROPERTY TYPE
81C-----------------------------------------------
82C I N P U T A r g u m e n t s
83C-----------------------------------------------
84 INTEGER NEL,NUVAR,IPT,NPTOT,NFUNC,IPG
85 INTEGER NGL(NEL),IFUNC(NFUNC)
86 my_real time,uparam(*),dpla(nel),uel1(nel),aldt(nel)
87 INTEGER, INTENT(INOUT) :: DMG_FLAG
88C-----------------------------------------------
89C I N P U T O U T P U T A r g u m e n t s
90C-----------------------------------------------
91 INTEGER ,DIMENSION(NEL), INTENT(INOUT) :: FOFF
92 my_real
93 . uvar(nel,nuvar),off(nel),offl(nel),
94 . signxx(nel),signyy(nel),signxy(nel),signyz(nel),signzx(nel),
95 . dfmax(nel),tdel(nel)
96 my_real, INTENT(INOUT) :: dmgscl(nel)
97C-----------------------------------------------
98C VARIABLES FOR FUNCTION INTERPOLATION
99C-----------------------------------------------
100 INTEGER NPF(*)
101 my_real finter ,tf(*)
102 EXTERNAL finter
103C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
104C Y : y = f(x)
105C X : x
106C DYDX : f'(x) = dy/dx
107C IFUNC(J): FUNCTION INDEX
108C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
109C NPF,TF : FUNCTION PARAMETER
110C-----------------------------------------------
111C L o c a l V a r i a b l e s
112C-----------------------------------------------
113 my_real hydros, vmises, triaxs,ref_el_len
114 INTEGER I,J,L,NINDX1,NINDX2,FLAG_PERTURB,SEL
115 INTEGER FAIL_COUNT,IT,ICOUP
116 INTEGER, DIMENSION(MVSIZ) :: INDX1,INDX2
117
118 my_real x_1(3) , x_2(3),c3,dydx,dc(nel)
119 my_real eps_fail, eps_fail2, damage, inst, dcrit, exp
120 my_real p1x,p1y,s1x,s1y,s2y, a1, a2, b1, b2, c1, c2,fac,lambda
121C-----------------------------------------------
122c! UVAR1 = damage due to instability (triax between 1/3 and 2/3)
123c! UVAR2 = actual integration point ON or OFF
124c! UVAR3-8 = perturbated parameter
125c! UVAR3 (if perturbation is not used) or UVAR9 (if used) = initial element length and then scale factor
126C-----------------------------------------------
127 x_1(1) = uparam(1)
128 x_1(2) = uparam(2)
129 x_1(3) = uparam(3)
130 x_2(1) = uparam(4)
131 x_2(2) = uparam(5)
132 x_2(3) = uparam(6)
133 flag_perturb = uparam(8)
134 c3 = uparam(9)
135 l = int(uparam(10)+0.0001)
136 sel = int(uparam(11)+0.0001)
137 inst = uparam(12)
138 ref_el_len = uparam(13)
139 nindx1 = 0
140 nindx2 = 0
141 eps_fail = zero
142 eps_fail2 = zero
143 icoup = nint(uparam(14))
144 dcrit = uparam(15)
145 exp = uparam(16)
146C
147C! Initialization
148 IF (nfunc > 0) THEN
149 IF (nuvar == 3) THEN
150 IF (uvar(1,3)==zero) THEN
151 DO i=1,nel
152 uvar(i,3) = aldt(i)
153 lambda = uvar(i,3) / ref_el_len
154 fac = finter(ifunc(1),lambda,npf,tf,dydx)
155 uvar(i,3) = fac
156 ENDDO
157 ENDIF
158 ELSEIF (nuvar == 9) THEN
159 IF (uvar(1,9)==zero) THEN
160 DO i=1,nel
161 uvar(i,9) = aldt(i)
162 lambda = uvar(i,9) / ref_el_len
163 fac = finter(ifunc(1),lambda,npf,tf,dydx)
164 uvar(i,9) = fac
165 ENDDO
166 ENDIF
167 ENDIF
168 ENDIF
169C-----------------------------------------------
170 IF(flag_perturb == 0) THEN
171 DO i =1,nel
172 IF(off(i) == one .AND. dpla(i) /= zero .AND. foff(i) == 1 ) THEN
173 hydros= (signxx(i)+ signyy(i))/3.0
174 vmises= sqrt((signxx(i)**2)+(signyy(i)**2)-(signxx(i)*signyy(i))+(3.0*signxy(i)**2))
175 triaxs = hydros / max(em20,vmises)
176 damage = uvar(i,1)
177 IF (triaxs <= third) THEN ! triax < 1/3
178 eps_fail = x_1(1) + x_1(2) * triaxs + x_1(3) * triaxs**2
179 IF ((nuvar == 3).AND.(nfunc > 0)) eps_fail = eps_fail * uvar(i,3)
180 ELSE ! triax > 1/3
181 SELECT CASE (sel)
182 CASE(1)
183 eps_fail = x_2(1) + x_2(2) * triaxs + x_2(3) * triaxs**2
184 IF ((nuvar == 3).AND.(nfunc > 0)) eps_fail = eps_fail * uvar(i,3)
185 CASE(2)
186 IF (triaxs <= one/sqr3) THEN ! triax < 0.57735
187 p1x = third
188 p1y = x_1(1) + x_1(2) * p1x + x_1(3) * p1x**2
189 s1x = one/sqr3
190 s1y = x_2(1) + x_2(2) / sqr3 + x_2(3) * (one/sqr3)**2
191 a1 = (p1y - s1y) / (p1x - s1x)**2
192 b1 = -two * a1 * s1x
193 c1 = a1 * s1x**2 + s1y
194 eps_fail = c1 + b1 * triaxs + a1 * triaxs**2
195 IF ((nuvar == 3).AND.(nfunc > 0)) eps_fail = eps_fail * uvar(i,3)
196 ELSE ! triax > 0.57735
197 p1x = two * third
198 p1y = x_2(1) + x_2(2) * p1x + x_2(3) * p1x**2
199 s1x = one/sqr3
200 s1y = x_2(1) + x_2(2) / sqr3 + x_2(3) * (one/sqr3)**2
201 a1 = (p1y - s1y) / (p1x - s1x)**2
202 b1 = -two * a1 * s1x
203 c1 = a1 * s1x**2 + s1y
204 eps_fail = c1 + b1 * triaxs + a1 * triaxs**2
205 IF ((nuvar == 3).AND.(nfunc > 0)) eps_fail = eps_fail * uvar(i,3)
206 ENDIF
207 CASE(3)
208 IF (triaxs <= one/sqr3) THEN ! triax < 0.57735
209 p1x = third
210 p1y = x_1(1) + x_1(2) * p1x + x_1(3) * p1x**2
211 s1x = one/sqr3
212 s1y = x_2(1) + x_2(2) / sqr3 + x_2(3) * (one/sqr3)**2
213 s2y = inst
214 a1 = (p1y - s1y) / (p1x - s1x)**2
215 a2 = (p1y - s2y) / (p1x - s1x)**2
216 b1 = -two * a1 * s1x
217 b2 = -two * a2 * s1x
218 c1 = a1 * s1x**2 + s1y
219 c2 = a2 * s1x**2 + s2y
220 eps_fail = c1 + b1 * triaxs + a1 * triaxs**2
221 IF ((nuvar == 3).AND.(nfunc > 0)) eps_fail = eps_fail * uvar(i,3)
222 eps_fail2 = c2 + b2 * triaxs + a2 * triaxs**2 ! INSTABILITY
223 IF ((nuvar == 3).AND.(nfunc > 0)) eps_fail2 = eps_fail2 * uvar(i,3)
224 ELSE ! triax > 0.57735
225 p1x = two * third
226 p1y = x_2(1) + x_2(2) * p1x + x_2(3) * p1x**2
227 s1x = one/sqr3
228 s1y = x_2(1) + x_2(2) / sqr3 + x_2(3) * (one/sqr3)**2
229 s2y = inst
230 a1 = (p1y - s1y) / (p1x - s1x)**2
231 a2 = (p1y - s2y) / (p1x - s1x)**2
232 b1 = -two * a1 * s1x
233 b2 = -two * a2 * s1x
234 c1 = a1 * s1x**2 + s1y
235 c2 = a2 * s1x**2 + s2y
236 eps_fail = c1 + b1 * triaxs + a1 * triaxs**2
237 IF ((nuvar == 3).AND.(nfunc > 0)) eps_fail = eps_fail * uvar(i,3)
238 eps_fail2 = c2 + b2 * triaxs + a2 * triaxs**2 ! INSTABILITY
239 IF ((nuvar == 3).AND.(nfunc > 0)) eps_fail2 = eps_fail2 * uvar(i,3)
240 ENDIF
241 IF (eps_fail2> zero) THEN
242 damage = damage + dpla(i)/eps_fail2
243 damage = min(one,damage)
244 uvar(i,1) = damage
245 IF ((damage >= one).AND.(uvar(i,2) == zero)) THEN
246 uvar(i,2) = dfmax(i)
247 ENDIF
248 ENDIF
249 END SELECT
250 ENDIF
251
252 dfmax(i) = dfmax(i) + dpla(i)/max(eps_fail,em6)
253 dfmax(i) = min(one,dfmax(i))
254
255 IF (dfmax(i) >= one) THEN
256 foff(i) = 0
257 nindx1 = nindx1 + 1
258 indx1(nindx1) = i
259 ENDIF
260
261 IF (damage >= one .AND. offl(i) == one .AND. icoup == 0) THEN
262 uel1(i) = uel1(i) + one
263 offl(i) = four_over_5
264 IF (int(uel1(i)) == nptot) THEN
265 nindx2 = nindx2 + 1
266 indx2(nindx2) = i
267 tdel(i) = time
268 off(i) = zero
269 signxx(i) = zero
270 signyy(i) = zero
271 signxy(i) = zero
272 signyz(i) = zero
273 signzx(i) = zero
274 ENDIF
275 ENDIF
276 ENDIF
277 ENDDO
278 ELSEIF(flag_perturb == 1) THEN
279 DO i =1,nel
280 IF(off(i) == one .AND. dpla(i) /= zero .AND. foff(i) == 1 ) THEN
281 hydros= (signxx(i)+ signyy(i))/3.0
282 vmises= sqrt((signxx(i)**2)+(signyy(i)**2)-(signxx(i)*signyy(i))+(3.0*signxy(i)**2))
283 triaxs = hydros / max(em20,vmises)
284 damage = uvar(i,1)
285 IF (triaxs <= third) THEN
286 eps_fail = uvar(i,3) + uvar(i,4) * triaxs + uvar(i,5) * triaxs**2
287 IF ((nuvar == 9).AND.(nfunc > 0)) eps_fail = eps_fail * uvar(i,9)
288 ELSE ! triax > 1/3
289 SELECT CASE (sel)
290 CASE(1)
291 eps_fail = uvar(i,6) + uvar(i,7) * triaxs + uvar(i,8) * triaxs**2
292 IF ((nuvar == 9).AND.(nfunc > 0)) eps_fail = eps_fail * uvar(i,9)
293 CASE(2)
294 IF (triaxs <= one/sqr3) THEN
295 p1x = third
296 p1y = uvar(i,3) + uvar(i,4) * p1x + uvar(i,5) * p1x**2
297 s1x = one/sqr3
298 s1y = uvar(i,6) + uvar(i,7) / sqr3 + uvar(i,8) * (one/sqr3)**2
299 a1 = (p1y - s1y) / (p1x - s1x)**2
300 b1 = -two * a1 * s1x
301 c1 = a1 * s1x**2 + s1y
302 eps_fail = c1 + b1 * triaxs + a1 * triaxs**2
303 IF ((nuvar == 9).AND.(nfunc > 0)) eps_fail = eps_fail * uvar(i,9)
304 ELSE
305 p1x = two * third
306 p1y = uvar(i,3) + uvar(i,4) * p1x + uvar(i,5) * p1x**2
307 s1x = one/sqr3
308 s1y = uvar(i,6) + uvar(i,7) / sqr3 + uvar(i,8) * (one/sqr3)**2
309 a1 = (p1y - s1y) / (p1x - s1x)**2
310 b1 = -two * a1 * s1x
311 c1 = a1 * s1x**2 + s1y
312 eps_fail = c1 + b1 * triaxs + a1 * triaxs**2
313 IF ((nuvar == 9).AND.(nfunc > 0)) eps_fail = eps_fail * uvar(i,9)
314 ENDIF
315 CASE(3)
316 IF (triaxs <= one/sqr3) THEN
317 p1x = third
318 p1y = uvar(i,3) + uvar(i,4) * p1x + uvar(i,5) * p1x**2
319 s1x = one/sqr3
320 s1y = uvar(i,6) + uvar(i,7) / sqr3 + uvar(i,8) * (one/sqr3)**2
321 s2y = inst
322 a1 = (p1y - s1y) / (p1x - s1x)**2
323 a2 = (p1y - s2y) / (p1x - s1x)**2
324 b1 = -two * a1 * s1x
325 b2 = -two * a2 * s1x
326 c1 = a1 * s1x**2 + s1y
327 c2 = a2 * s1x**2 + s2y
328 eps_fail = c1 + b1 * triaxs + a1 * triaxs**2
329 IF ((nuvar == 9).AND.(nfunc > 0)) eps_fail = eps_fail * uvar(i,9)
330 eps_fail2 = c2 + b2 * triaxs + a2 * triaxs**2
331 IF ((nuvar == 9).AND.(nfunc > 0)) eps_fail2 = eps_fail2 * uvar(i,9)
332 ELSE
333 p1x = two * third
334 p1y = uvar(i,3) + uvar(i,4) * p1x + uvar(i,5) * p1x**2
335 s1x = one/sqr3
336 s1y = uvar(i,6) + uvar(i,7) / sqr3 + uvar(i,8) * (one/sqr3)**2
337 s2y = inst
338 a1 = (p1y - s1y) / (p1x - s1x)**2
339 a2 = (p1y - s2y) / (p1x - s1x)**2
340 b1 = -two * a1 * s1x
341 b2 = -two * a2 * s1x
342 c1 = a1 * s1x**2 + s1y
343 c2 = a2 * s1x**2 + s2y
344 eps_fail = c1 + b1 * triaxs + a1 * triaxs**2
345 IF ((nuvar == 9).AND.(nfunc > 0)) eps_fail = eps_fail * uvar(i,9)
346 eps_fail2 = c2 + b2 * triaxs + a2 * triaxs**2
347 IF ((nuvar == 9).AND.(nfunc > 0)) eps_fail2 = eps_fail2 * uvar(i,9)
348 ENDIF
349 IF (eps_fail2> zero) THEN
350 damage = damage + dpla(i)/eps_fail2
351 damage = min(one,damage)
352 uvar(i,1) = damage
353 IF ((damage >= one).AND.(uvar(i,2) == zero)) THEN
354 uvar(i,2) = dfmax(i)
355 ENDIF
356 ENDIF
357 END SELECT
358 ENDIF
359
360 dfmax(i) = dfmax(i) + dpla(i)/max(eps_fail,em6)
361 dfmax(i) = min(one,dfmax(i))
362
363 IF (dfmax(i) >= one) THEN
364 foff(i) = 0
365 nindx1 = nindx1 + 1
366 indx1(nindx1) = i
367 ENDIF
368
369 IF (damage >= one .AND. offl(i) == one .AND. icoup == 0) THEN
370 uel1(i) = uel1(i) + one
371 offl(i) = four_over_5
372 IF (int(uel1(i)) == nptot) THEN
373 nindx2 = nindx2 + 1
374 indx2(nindx2) = i
375 tdel(i) = time
376 off(i) = zero
377 signxx(i) = zero
378 signyy(i) = zero
379 signxy(i) = zero
380 signyz(i) = zero
381 signzx(i) = zero
382 ENDIF
383 ENDIF
384 ENDIF
385 ENDDO
386
387 ENDIF
388c
389c------------------------
390c STRESS SOFTENING
391c------------------------
392 IF (icoup > 0) THEN
393 dmg_flag = 1
394 IF (icoup == 1) THEN
395 DO i = 1,nel
396 IF (dfmax(i) >= dcrit) THEN
397 IF (dcrit < one) THEN
398 dmgscl(i) = one - ((dfmax(i)-dcrit)/max(one-dcrit,em20))**exp
399 ELSE
400 dmgscl(i) = zero
401 ENDIF
402 ELSE
403 dmgscl(i) = one
404 ENDIF
405 ENDDO
406 ELSE
407 DO i = 1,nel
408 dc(i) = uvar(i,2)
409 IF (dc(i) == zero) dc(i) = one
410 IF (dfmax(i) >= dc(i)) THEN
411 IF (dc(i) < one) THEN
412 dmgscl(i) = one - ((dfmax(i)-dc(i))/max(one-dc(i),em20))**exp
413 ELSE
414 dmgscl(i) = zero
415 ENDIF
416 ELSE
417 dmgscl(i) = one
418 ENDIF
419 ENDDO
420 ENDIF
421 ENDIF
422c------------------------
423c------------------------
424c------------------------
425c------------------------
426 IF (nindx1 > 0) THEN
427 DO j=1,nindx1
428 i = indx1(j)
429#include "lockon.inc"
430 WRITE(iout, 2000) ngl(i),ipg,ipt,time
431 WRITE(istdo,2000) ngl(i),ipg,ipt,time
432#include "lockoff.inc"
433 ENDDO
434 ENDIF
435 IF (nindx2 > 0) THEN
436 DO j=1,nindx2
437 i = indx2(j)
438#include "lockon.inc"
439 WRITE(iout, 3000) ngl(i)
440 WRITE(iout, 2200) ngl(i), time
441 WRITE(istdo,3000) ngl(i)
442 WRITE(istdo,2200) ngl(i), time
443#include "lockoff.inc"
444 END DO
445 END IF ! NINDX
446c------------------------
447C---------Damage for output 0 < DFMAX < 1 --------------------
448 2000 FORMAT(1x,'FOR SHELL ELEMENT (BIQUAD)',i10,1x,'GAUSS POINT',i3,
449 . 1x,'LAYER',i3,':',/,
450 . 1x,'STRESS TENSOR SET TO ZERO',1x,'AT TIME :',1pe12.4)
451 2200 FORMAT(1x,' *** RUPTURE OF SHELL ELEMENT (BIQUAD)',i10,1x,
452 . ' AT TIME :',1pe12.4)
453 3000 FORMAT(1x,'FOR SHELL ELEMENT (BIQUAD)',i10,
454 . 1x,'INSTABILITY REACHED.')
455c------------------------
456 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