OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_biquad_s.F File Reference
#include "implicit_f.inc"
#include "scr17_c.inc"
#include "tabsiz_c.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_biquad_s (nel, nuparam, nuvar, mfunc, kfunc, aldt, npf, tf, time, uparam, tdele, ngl, dpla, uvar, off, dfmax, dmgscl, signxx, signyy, signzz, signxy, signyz, signzx)

Function/Subroutine Documentation

◆ fail_biquad_s()

subroutine fail_biquad_s ( integer, intent(in) nel,
integer, intent(in) nuparam,
integer, intent(in) nuvar,
integer mfunc,
integer, dimension(mfunc) kfunc,
dimension(nel), intent(in) aldt,
integer, dimension(snpc) npf,
tf,
intent(in) time,
dimension(nuparam), intent(in) uparam,
dimension(nel), intent(inout) tdele,
integer, dimension(nel), intent(in) ngl,
dimension(nel), intent(in) dpla,
dimension(nel,nuvar), intent(inout) uvar,
dimension(nel), intent(inout) off,
dimension(nel), intent(inout) dfmax,
dimension(nel), intent(inout) dmgscl,
dimension(nel), intent(in) signxx,
dimension(nel), intent(in) signyy,
dimension(nel), intent(in) signzz,
dimension(nel), intent(in) signxy,
dimension(nel), intent(in) signyz,
dimension(nel), intent(in) signzx )

Definition at line 34 of file fail_biquad_s.F.

39C-----------------------------------------------
40C I m p l i c i t T y p e s
41C-----------------------------------------------
42#include "implicit_f.inc"
43C-----------------------------------------------
44C C o m m o n B l o c k s
45C-----------------------------------------------
46#include "scr17_c.inc"
47#include "tabsiz_c.inc"
48#include "units_c.inc"
49#include "comlock.inc"
50C-----------------------------------------------
51C I N P U T A r g u m e n t s
52C-----------------------------------------------
53 INTEGER, INTENT(IN) :: NEL,NUPARAM,NUVAR,NGL(NEL)
54 my_real, INTENT(IN) :: time,uparam(nuparam),
55 . signxx(nel),signyy(nel),signzz(nel),
56 . signxy(nel),signyz(nel),signzx(nel),
57 . dpla(nel),aldt(nel)
58 my_real, INTENT(INOUT) :: uvar(nel,nuvar),dfmax(nel),
59 . tdele(nel),off(nel),dmgscl(nel)
60C!-----------------------------------------------
61C! VARIABLES FOR FUNCTION INTERPOLATION
62C!-----------------------------------------------
63 INTEGER NPF(SNPC), MFUNC, KFUNC(MFUNC)
64 my_real finter ,tf(stf)
65 EXTERNAL finter
66C-----------------------------------------------
67C L o c a l V a r i a b l e s
68C-----------------------------------------------
69 INTEGER I,J,INDX(NEL),NINDX,SEL,FLAG_PERTURB,ICOUP
70 my_real
71 . df,p,triaxs,svm,sxx,syy,szz,eps_fail,p1x,p1y,s1x,s1y,
72 . s2y,a1,b1,c1,ref_el_len,lambda,fac,x_1(3),x_2(3),dcrit,exp
73C
74 !=======================================================================
75 ! - INITIALISATION OF COMPUTATION ON TIME STEP
76 !=======================================================================
77 ! Recovering failure criterion parameters
78 x_1(1) = uparam(1)
79 x_1(2) = uparam(2)
80 x_1(3) = uparam(3)
81 x_2(1) = uparam(4)
82 x_2(2) = uparam(5)
83 x_2(3) = uparam(6)
84 flag_perturb = uparam(8)
85 sel = int(uparam(11)+0.0001)
86 IF (sel == 3) sel = 2
87 ref_el_len = uparam(13)
88 icoup = nint(uparam(14))
89 dcrit = uparam(15)
90 exp = uparam(16)
91 eps_fail = zero
92
93C
94 ! At initial time, compute the element size regularization factor
95 IF (mfunc > 0) THEN
96 IF (nuvar == 3) THEN
97 IF (uvar(1,3) == zero) THEN
98 DO i = 1,nel
99 uvar(i,3) = aldt(i)
100 lambda = uvar(i,3) / ref_el_len
101 fac = finter(kfunc(1),lambda,npf,tf,df)
102 uvar(i,3) = fac
103 ENDDO
104 ENDIF
105 ELSEIF (nuvar == 9) THEN
106 IF (uvar(1,9) == zero) THEN
107 DO i = 1,nel
108 uvar(i,9) = aldt(i)
109 lambda = uvar(i,9) / ref_el_len
110 fac = finter(kfunc(1),lambda,npf,tf,df)
111 uvar(i,9) = fac
112 ENDDO
113 ENDIF
114 ENDIF
115 ENDIF
116C
117 ! Update element deletion flag
118 DO i = 1,nel
119 IF (off(i) < em01) off(i) = zero
120 IF (off(i) < one .AND. off(i) > zero) off(i) = off(i)*four_over_5
121 END DO
122C
123 !====================================================================
124 ! - LOOP OVER THE ELEMENT TO UPDATE DAMAGE VARIABLE
125 !====================================================================
126 ! Initialize deleted element counter
127 nindx = 0
128C
129 ! -> If perturbation option is not used, constant parameters
130 IF (flag_perturb == 0) THEN
131 DO i = 1,nel
132 IF (off(i) == one .AND. dpla(i) /= zero) THEN
133 ! Compute hydrostatic stress
134 p = third*(signxx(i) + signyy(i) + signzz(i))
135C
136 ! Compute Von Mises stress
137 sxx = signxx(i) - p
138 syy = signyy(i) - p
139 szz = signzz(i) - p
140 svm = half*(sxx**2 + syy**2 + szz**2)
141 . + signxy(i)**2 + signzx(i)**2 + signyz(i)**2
142 svm = sqrt(three*svm)
143C
144 ! Compute stress triaxiality
145 triaxs = p/max(em20,svm)
146 IF (triaxs < -two_third) triaxs = -two_third
147 IF (triaxs > two_third) triaxs = two_third
148C
149 ! Compute the corresponding plastic strain at failure
150 ! -> Low stress triaxialities parabolic curve
151 IF (triaxs <= third) THEN
152 eps_fail = x_1(1) + x_1(2)*triaxs + x_1(3)*triaxs**2
153 IF ((nuvar == 3).AND.(mfunc > 0)) eps_fail = eps_fail*uvar(i,3)
154 ! -> High stress triaxiality parabolic curve
155 ELSE
156 ! Raw curve is used
157 IF (sel == 1) THEN
158 eps_fail = x_2(1) + x_2(2)*triaxs + x_2(3)*triaxs**2
159 IF ((nuvar == 3).AND.(mfunc > 0)) eps_fail = eps_fail*uvar(i,3)
160 ! Plain strain is global minimum
161 ELSEIF (sel == 2) THEN
162 IF (triaxs <= one/sqr3) THEN ! triax < 0.57735
163 p1x = third
164 p1y = x_1(1) + x_1(2)*p1x + x_1(3)*p1x**2
165 s1x = one/sqr3
166 s1y = x_2(1) + x_2(2)/sqr3 + x_2(3)*(one/sqr3)**2
167 a1 = (p1y - s1y)/(p1x - s1x)**2
168 b1 = -two*a1*s1x
169 c1 = a1*s1x**2 + s1y
170 eps_fail = c1 + b1*triaxs + a1*triaxs**2
171 IF ((nuvar == 3).AND.(mfunc > 0)) eps_fail = eps_fail*uvar(i,3)
172 ELSE ! triax > 0.57735
173 p1x = two*third
174 p1y = x_2(1) + x_2(2)*p1x + x_2(3)*p1x**2
175 s1x = one/sqr3
176 s1y = x_2(1) + x_2(2)/sqr3 + x_2(3)*(one/sqr3)**2
177 a1 = (p1y - s1y)/(p1x - s1x)**2
178 b1 = -two*a1*s1x
179 c1 = a1*s1x**2 + s1y
180 eps_fail = c1 + b1*triaxs + a1*triaxs**2
181 IF ((nuvar == 3).AND.(mfunc > 0)) eps_fail = eps_fail*uvar(i,3)
182 ENDIF
183 ENDIF
184 ENDIF
185C
186 ! Update damage variable
187 dfmax(i) = dfmax(i) + dpla(i)/max(eps_fail,em6)
188 dfmax(i) = min(one,dfmax(i))
189 ! Check element failure
190 IF (dfmax(i) >= one .AND. off(i) == one) THEN
191 off(i) = four_over_5
192 nindx = nindx+1
193 indx(nindx) = i
194 tdele(i) = time
195 ENDIF
196 ENDIF
197 ENDDO
198 ! -> Perturbation flag is used, parameters vary for each Gauss points
199 ELSE
200 DO i = 1,nel
201 IF (off(i) == one .AND. dpla(i) /= zero) THEN
202 ! Compute hydrostatic stress
203 p = third*(signxx(i) + signyy(i) + signzz(i))
204C
205 ! Compute Von Mises stress
206 sxx = signxx(i) - p
207 syy = signyy(i) - p
208 szz = signzz(i) - p
209 svm = half*(sxx**2 + syy**2 + szz**2)
210 . + signxy(i)**2 + signzx(i)**2 + signyz(i)**2
211 svm = sqrt(three*svm)
212C
213 ! Compute stress triaxiality
214 triaxs = p/max(em20,svm)
215 IF (triaxs < -two_third) triaxs = -two_third
216 IF (triaxs > two_third) triaxs = two_third
217C
218 ! Compute the corresponding plastic strain at failure
219 ! -> Low stress triaxialities parabolic curve
220 IF (triaxs <= third) THEN
221 eps_fail = uvar(i,3) + uvar(i,4)*triaxs + uvar(i,5)*triaxs**2
222 IF ((nuvar == 9).AND.(mfunc > 0)) eps_fail = eps_fail*uvar(i,9)
223 ! -> High stress triaxiality parabolic curve
224 ELSE
225 ! Raw curve is used
226 IF (sel == 1) THEN
227 eps_fail = uvar(i,6) + uvar(i,7)*triaxs + uvar(i,8)*triaxs**2
228 IF ((nuvar == 9).AND.(mfunc > 0)) eps_fail = eps_fail*uvar(i,9)
229 ! Plain strain is global minimum
230 ELSEIF (sel == 2) THEN
231 IF (triaxs <= one/sqr3) THEN ! triax < 0.57735
232 p1x = third
233 p1y = uvar(i,3) + uvar(i,4)*p1x + uvar(i,5)*p1x**2
234 s1x = one/sqr3
235 s1y = uvar(i,6) + uvar(i,7)/sqr3 + uvar(i,8)*(one/sqr3)**2
236 a1 = (p1y - s1y)/(p1x - s1x)**2
237 b1 = -two*a1*s1x
238 c1 = a1*s1x**2 + s1y
239 eps_fail = c1 + b1*triaxs + a1*triaxs**2
240 IF ((nuvar == 9).AND.(mfunc > 0)) eps_fail = eps_fail*uvar(i,9)
241 ELSE ! triax > 0.57735
242 p1x = two*third
243 p1y = uvar(i,3) + uvar(i,4)*p1x + uvar(i,5)*p1x**2
244 s1x = one/sqr3
245 s1y = uvar(i,6) + uvar(i,7)/sqr3 + uvar(i,8)*(one/sqr3)**2
246 a1 = (p1y - s1y)/(p1x - s1x)**2
247 b1 = -two*a1*s1x
248 c1 = a1*s1x**2 + s1y
249 eps_fail = c1 + b1*triaxs + a1*triaxs**2
250 IF ((nuvar == 9).AND.(mfunc > 0)) eps_fail = eps_fail*uvar(i,9)
251 ENDIF
252 ENDIF
253 ENDIF
254C
255 ! Update damage variable
256 dfmax(i) = dfmax(i) + dpla(i)/max(eps_fail,em6)
257 dfmax(i) = min(one,dfmax(i))
258 ! Check element failure
259 IF (dfmax(i) >= one .AND. off(i) == one) THEN
260 off(i) = four_over_5
261 nindx = nindx+1
262 indx(nindx) = i
263 tdele(i) = time
264 ENDIF
265 ENDIF
266 ENDDO
267 ENDIF
268c
269c------------------------
270c STRESS SOFTENING
271c------------------------
272 IF (icoup == 1) THEN
273 DO i = 1,nel
274 IF (dfmax(i) >= dcrit) THEN
275 IF (dcrit < one) THEN
276 dmgscl(i) = one - ((dfmax(i)-dcrit)/max(one-dcrit,em20))**exp
277 ELSE
278 dmgscl(i) = zero
279 ENDIF
280 ELSE
281 dmgscl(i) = one
282 ENDIF
283 ENDDO
284 ENDIF
285c
286 !====================================================================
287 ! - PRINTOUT ELEMENT FAILURE MESSAGES
288 !====================================================================
289 IF (nindx > 0) THEN
290 DO j=1,nindx
291 i = indx(j)
292#include "lockon.inc"
293 WRITE(iout, 1000) ngl(i),time
294 WRITE(istdo,1100) ngl(i),time
295#include "lockoff.inc"
296 ENDDO
297 ENDIF
298C------------------
299 1000 FORMAT(1x,' -- RUPTURE OF SOLID ELEMENT (BIQUAD)',i10,
300 . ' AT TIME : ',1pe12.4)
301 1100 FORMAT(1x,' -- RUPTURE OF SOLID ELEMENT (BIQUAD)',i10,
302 . ' AT TIME : ',1pe12.4)
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