OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_biquad_ib.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_biquad_ib ../engine/source/materials/fail/biquad/fail_biquad_ib.F
25!||--- called by ------------------------------------------------------
26!|| fail_beam18 ../engine/source/elements/beam/fail_beam18.F
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!||====================================================================
30 SUBROUTINE fail_biquad_ib(
31 . NEL ,NGL ,NUPARAM ,UPARAM ,
32 . TIME ,OFF ,DFMAX ,TDEL ,
33 . IOUT ,ISTDO ,NFUNC ,IFUNC ,DAMSCL ,
34 . UVAR ,NUVAR ,SNPC ,NPF ,
35 . STF ,TF ,IPT ,FOFF ,
36 . SIGNXX ,SIGNXY ,SIGNXZ ,DPLA ,AL , UELR,NPG)
37C-----------------------------------------------C-----------------------------------------------
38c BIQUAD- failure model for integrated beams (TYPE 18) IRUP=30
39C-----------------------------------------------
40C I m p l i c i t T y p e s
41C-----------------------------------------------
42#include "implicit_f.inc"
43C-----------------------------------------------
44C G l o b a l P a r a m e t e r s
45C-----------------------------------------------
46#include "comlock.inc"
47C-----------------------------------------------
48C D u m m y A r g u m e n t s
49C-----------------------------------------------
50 INTEGER ,INTENT(IN) :: NEL ! size of element group
51 INTEGER ,INTENT(IN) :: NUPARAM ! size of parameter array
52 INTEGER ,INTENT(IN) :: IPT ! current integration point
53 INTEGER ,INTENT(IN) :: IOUT ! output file unit
54 INTEGER ,INTENT(IN) :: ISTDO ! output file unit
55 INTEGER ,INTENT(IN) :: NFUNC ! number of functions
56 INTEGER ,INTENT(IN) :: SNPC
57 INTEGER ,INTENT(IN) :: STF
58 INTEGER ,INTENT(IN) :: NUVAR
59 INTEGER ,DIMENSION(NEL) ,INTENT(IN) :: NGL ! table of element identifiers
60 INTEGER ,DIMENSION(100) ,INTENT(IN) :: IFUNC ! table of functions identifiers
61 INTEGER ,DIMENSION(SNPC) ,INTENT(IN) :: NPF
62 INTEGER ,DIMENSION(NEL) ,INTENT(INOUT) :: FOFF ! integration point deactivation flag
63 my_real ,INTENT(IN) :: time ! current time
64 my_real ,DIMENSION(NUPARAM) ,INTENT(IN) :: uparam ! failure model parameter array
65 my_real ,DIMENSION(NEL) ,INTENT(IN) :: dpla ! plastic strain
66 my_real ,DIMENSION(NEL) ,INTENT(IN) :: al
67 my_real ,DIMENSION(STF) ,INTENT(IN) :: tf
68 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: off ! element deactivation flag
69 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: dfmax ! maximum damage
70 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: damscl !1-damage
71 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: tdel ! deactivation time
72 my_real ,DIMENSION(NEL) ,INTENT(IN) :: signxx,signxy,signxz
73 my_real ,DIMENSION(NEL,NUVAR),INTENT(INOUT) :: uvar
74 my_real, DIMENSION(NEL) ,INTENT(INOUT) :: uelr
75 INTEGER ,INTENT(IN) :: NPG
76C-----------------------------------------------
77C L o c a l V a r i a b l e s
78C-----------------------------------------------
79 INTEGER :: I,J,NINDX,NINDXD,ICOUP,SEL,FAILIP
80 INTEGER ,DIMENSION(NEL) :: INDX,INDXD
81 my_real :: dcrit, exp,d1,df,triax,scale,sxx,syy,szz
82
83 my_real p1x,p1y,s1x,s1y,s2y, a1, b1, c1, ref_el_len, lambda,fac ,
84 . p , svm,triaxs
85 my_real x_1(3) , x_2(3)
86 EXTERNAL finter
87 my_real eps_fail,finter
88 my_real ,DIMENSION(NEL) :: damage
89C=======================================================================
90 eps_fail = zero
91 x_1(1) = uparam(1)
92 x_1(2) = uparam(2)
93 x_1(3) = uparam(3)
94 x_2(1) = uparam(4)
95 x_2(2) = uparam(5)
96 x_2(3) = uparam(6)
97
98 sel = int(uparam(11)+0.0001)
99 IF (sel == 3) sel = 2
100 ref_el_len = uparam(13)
101 icoup = nint(uparam(14))
102 dcrit = uparam(15)
103 exp = uparam(16)
104 failip =min(nint(uparam(17)),npg) ! number of failed integration point prior to solid element deletion. default = 1 (integer)
105
106 ! At initial time, compute the element size regularization factor
107 IF (nfunc > 0) THEN
108 IF (nuvar == 3) THEN
109 IF (uvar(1,3) == zero) THEN
110 DO i = 1,nel
111 uvar(i,3) = al(i)
112 lambda = al(i) / ref_el_len
113 fac = finter(ifunc(1),lambda,npf,tf,df)
114 uvar(i,3) = fac
115 ENDDO
116 ENDIF
117 ELSEIF (nuvar == 9) THEN
118 IF (uvar(1,9) == zero) THEN
119 DO i = 1,nel
120 uvar(i,9) = al(i)
121 lambda = uvar(i,9) / ref_el_len
122 fac = finter(ifunc(1),lambda,npf,tf,df)
123 uvar(i,9) = fac
124 ENDDO
125 ENDIF
126 ENDIF
127 ENDIF
128 !====================================================================
129 ! - LOOP OVER THE ELEMENT TO UPDATE DAMAGE VARIABLE
130 !====================================================================
131 ! Initialize deleted element counter
132 nindx = 0
133C
134 DO i = 1,nel
135 IF (off(i) == one .AND. foff(i) == 1 .AND. dpla(i) /= zero) THEN
136 ! Compute hydrostatic stress
137 p = third*signxx(i)
138 ! Compute Von Mises stress
139 svm = sqrt(three*(half*signxx(i)**2 + signxy(i)**2 + signxz(i)**2))
140 ! Compute stress triaxiality
141 triaxs = p/max(em20,svm)
142 IF (triaxs < -two_third) triaxs = -two_third
143 IF (triaxs > two_third) triaxs = two_third
144C
145 ! Compute the corresponding plastic strain at failure
146 ! -> Low stress triaxialities parabolic curve
147 IF (triaxs <= third) THEN
148 eps_fail = x_1(1) + x_1(2)*triaxs + x_1(3)*triaxs**2
149 IF ((nuvar == 3).AND.(nfunc > 0)) eps_fail = eps_fail* uvar(i,3)
150 ! -> High stress triaxiality parabolic curve
151 ELSE
152 ! Raw curve is used
153 IF (sel == 1) THEN
154 eps_fail = x_2(1) + x_2(2)*triaxs + x_2(3)*triaxs**2
155 IF ((nuvar == 3).AND.(nfunc > 0)) eps_fail = eps_fail*uvar(i,3)
156 ! plain strain is global minimum
157 ELSEIF (sel == 2) THEN
158 IF (triaxs <= one/sqr3) THEN ! triax < 0.57735
159 p1x = third
160 p1y = x_1(1) + x_1(2)*p1x + x_1(3)*p1x**2
161 s1x = one/sqr3
162 s1y = x_2(1) + x_2(2)/sqr3 + x_2(3)*(one/sqr3)**2
163 a1 = (p1y - s1y)/(p1x - s1x)**2
164 b1 = -two*a1*s1x
165 c1 = a1*s1x**2 + s1y
166 eps_fail = c1 + b1*triaxs + a1*triaxs**2
167 IF ((nuvar == 3).AND.(nfunc > 0)) eps_fail = eps_fail*uvar(i,3)
168 ELSE ! triax > 0.57735
169 p1x = two*third
170 p1y = x_2(1) + x_2(2)*p1x + x_2(3)*p1x**2
171 s1x = one/sqr3
172 s1y = x_2(1) + x_2(2)/sqr3 + x_2(3)*(one/sqr3)**2
173 a1 = (p1y - s1y)/(p1x - s1x)**2
174 b1 = -two*a1*s1x
175 c1 = a1*s1x**2 + s1y
176 eps_fail = c1 + b1*triaxs + a1*triaxs**2
177 IF ((nuvar == 3).AND.(nfunc> 0)) eps_fail = eps_fail*uvar(i,3)
178 ENDIF
179 ENDIF
180 ENDIF
181C
182 ! Update damage variable
183 dfmax(i) = dfmax(i) + dpla(i)/max(eps_fail,em6)
184 dfmax(i) = min(one,dfmax(i))
185 ! Check element failure
186 IF (dfmax(i) >= one ) THEN
187 nindx = nindx + 1
188 indx(nindx) = i
189 tdel(i) = time
190 dfmax(i) = one
191 damscl(i) = zero
192 foff(i) = 0
193 uelr(i) = uelr(i) + one
194 IF (nint(uelr(i)) >= failip) then
195 off(i) = four_over_5
196 ENDIF
197 ENDIF
198 ENDIF
199 ENDDO
200c------------------------
201c STRESS SOFTENING
202c------------------------
203 IF (icoup == 1) THEN
204 DO i = 1,nel
205 IF (dfmax(i) >= dcrit .AND.off(i) == one .AND. foff(i) == 1) THEN
206 IF (dcrit < one) THEN
207 damscl(i) = one - ((dfmax(i)-dcrit)/max(one-dcrit,em20))**exp
208 ELSE
209 damscl(i) = zero
210 ENDIF
211 ELSE
212 damscl(i) = one
213 ENDIF
214 ENDDO
215 ELSE
216 DO i=1,nel
217 damscl(i) = one - dfmax(i)
218 END DO
219
220 ENDIF
221
222c--------------------
223 IF (nindx > 0) THEN
224 DO j=1,nindx
225 i = indx(j)
226#include "lockon.inc"
227 WRITE(iout, 1000) ngl(i),ipt,time
228 WRITE(istdo,1000) ngl(i),ipt,time
229#include "lockoff.inc"
230 IF (off(i) == four_over_5) THEN
231 WRITE(iout, 1111)ngl(i),time
232 WRITE(istdo,1111)ngl(i),time
233 ENDIF
234 END DO
235 END IF ! NINDX
236c------------------
237 1000 FORMAT(5x,'FAILURE (BIQUAD) OF BEAM ELEMENT ',i10,1x,',INTEGRATION PT',i5
238 . ,2x,'AT TIME :',1pe12.4,1x )
239 2000 FORMAT(5x,'start damage(biquad) of beam element ',I10,1X,',integration pt',I5
240 . ,2X,'at time :',1PE12.4 )
241 1111 FORMAT(1X,'deleted beam element ',I10,1X,'at time :',1PE12.4)
242c------------------
243 RETURN
244 END
245
my_real
#define my_real
Definition cppsort.cpp:32
fail_biquad_ib
subroutine fail_biquad_ib(nel, ngl, nuparam, uparam, time, off, dfmax, tdel, iout, istdo, nfunc, ifunc, damscl, uvar, nuvar, snpc, npf, stf, tf, ipt, foff, signxx, signxy, signxz, dpla, al, uelr, npg)
Definition fail_biquad_ib.F:37
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21