OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_biquad_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_biquad_s ../engine/source/materials/fail/biquad/fail_biquad_s.F
25!||--- called by ------------------------------------------------------
26!|| mmain ../engine/source/materials/mat_share/mmain.F90
27!|| mmain8 ../engine/source/materials/mat_share/mmain8.F
28!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
29!|| mulaw8 ../engine/source/materials/mat_share/mulaw8.F90
30!|| usermat_solid ../engine/source/materials/mat_share/usermat_solid.f
31!||--- calls -----------------------------------------------------
32!|| finter ../engine/source/tools/curve/finter.F
33!||====================================================================
34 SUBROUTINE fail_biquad_s(
35 1 NEL ,NUPARAM ,NUVAR ,MFUNC ,KFUNC ,ALDT ,
36 2 NPF ,TF ,TIME ,UPARAM ,TDELE ,
37 3 NGL ,DPLA ,UVAR ,OFF ,DFMAX ,DMGSCL ,
38 4 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX )
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)
303 END
fail_biquad_s
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)
Definition fail_biquad_s.F:39
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
usermat_solid
subroutine usermat_solid(timers, lft, llt, nft, mtn, jcvt, pm, off, sig, eint, rho, qold, vol, strain, sigl, gama, uvar, bufmat, tf, npf, stifn, mat, ngl, nuvar, dt2t, neltst, ityptst, offg, geo, pid, epsd, el_temp, wxx, wyy, wzz, jsph, mumax, ssp, aire, voln, vd2, deltax, vis, d1, d2, d3, d4, d5, d6, pnew, psh, q, ssp_eq, dvol, sold1, sold2, sold3, sold4, sold5, sold6, rx, ry, rz, sx, sy, sz, tx, ty, tz, ipla, sigy, defp, ismstr, mfxx, mfxy, mfxz, mfyx, mfyy, mfyz, mfzx, mfzy, mfzz, ipm, isorth, fbuf, nfail, npg, sigdd, dxy, dyx, dyz, dzy, dzx, dxz, fr_wav, isrot, v, varnl, w, ix, x, jthe, et, mssa, dmels, iptr, ipts, iptt, table, fvd2, fdeltax, fssp, fqvis, iparg1, igeo, sigv, al_imp, signor, istrain, ng, elbuf_tab, vbuf, ilay, vk, iparg, bufvois, vdx, vdy, vdz, ihet, conde, itask, iexpan, vol_avg, amu, epsth3, epsth, svisc, nel, etotsh, iselect, tstar, muold, amu2, dpdm, rhoref, rhosp, nloc_dmg, ity, jtur, mat_elem, idel7nok, svis, dt, glob_therm, damp_buf, idamp_freq_range)
Definition usermat_solid.F:131