OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_syazwan_c.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_syazwan_c ../engine/source/materials/fail/syazwan/fail_syazwan_c.F
25!||--- called by ------------------------------------------------------
26!|| mulawc ../engine/source/materials/mat_share/mulawc.F90
27!|| usermat_shell ../engine/source/materials/mat_share/usermat_shell.F
28!||--- calls -----------------------------------------------------
29!|| finter ../engine/source/tools/curve/finter.F
30!||====================================================================
31 SUBROUTINE fail_syazwan_c(
32 1 NEL ,UPARAM ,NUPARAM ,UVAR ,NUVAR ,
33 2 TIME ,NGL ,IPT ,DPLA ,PLA ,
34 3 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
35 4 EPSXX ,EPSYY ,EPSXY ,EPSYZ ,EPSZX ,
36 5 DFMAX ,NFUNC ,IFUNC ,ALDT ,FOFF ,
37 6 IPG ,DMG_FLAG ,DMG_SCALE,NPF ,TF )
38C-----------------------------------------------
39C I m p l i c i t T y p e s
40C-----------------------------------------------
41#include "implicit_f.inc"
42C-----------------------------------------------
43C G l o b a l P a r a m e t e r s
44C-----------------------------------------------
45#include "mvsiz_p.inc"
46#include "units_c.inc"
47#include "comlock.inc"
48#include "com01_c.inc"
49#include "tabsiz_c.inc"
50C-----------------------------------------------
51C I N P U T A r g u m e n t s
52C-----------------------------------------------
53 INTEGER, INTENT(IN) ::
54 . NEL ,NUPARAM,NUVAR ,IPT ,NFUNC ,
55 . IPG ,NGL(NEL),IFUNC(NFUNC)
56 my_real, INTENT(IN) ::
57 . TIME ,UPARAM(NUPARAM),DPLA(NEL) ,PLA(NEL),
58 . ALDT(NEL),EPSXX(NEL) ,EPSYY(NEL),
59 . epsxy(nel),epsyz(nel) ,epszx(nel)
60C-----------------------------------------------
61C I N P U T O U T P U T A r g u m e n t s
62C-----------------------------------------------
63 INTEGER, INTENT(INOUT) ::
64 . DMG_FLAG,FOFF(NEL)
65 my_real, INTENT(INOUT) ::
66 . UVAR(NEL,NUVAR),DFMAX(NEL),DMG_SCALE(NEL) ,
67 . SIGNXX(NEL),SIGNYY(NEL),SIGNXY(NEL),SIGNYZ(NEL),
68 . SIGNZX(NEL)
69C-----------------------------------------------
70C VARIABLES FOR FUNCTION INTERPOLATION
71C-----------------------------------------------
72 INTEGER, INTENT(IN) :: NPF(SNPC)
73 my_real, INTENT(IN) :: TF(STF)
74 my_real
75 . FINTER
76 EXTERNAL FINTER
77C-----------------------------------------------
78C L o c a l V a r i a b l e s
79C-----------------------------------------------
80 INTEGER I,J,NINDX,INST,DINIT,IFORM
81 INTEGER, DIMENSION(MVSIZ) :: INDX
82 my_real
83 . c1 ,c2 ,c3 ,c4 ,c5 ,c6 ,
84 . n_val ,softexp,ref_len,reg_scale,dam_sf,max_dam
85 my_real
86 . lambda ,dydx ,fac(nel),dc(nel),p ,svm ,
87 . triax ,cos3theta,lodep ,epsfail,e12 ,s1 ,
88 . s2 ,q ,emaj ,emin ,beta ,alpha ,
89 . equa1 ,equa2 ,e1 ,eps_c ,dinst(nel),epfmin,
90 . delta ,denom
91 !=======================================================================
92 ! - INITIALISATION OF COMPUTATION ON TIME STEP
93 !=======================================================================
94 ! Recovering failure criterion parameters
95 c1 = uparam(1)
96 c2 = uparam(2)
97 c3 = uparam(3)
98 c4 = uparam(4)
99 c5 = uparam(5)
100 c6 = uparam(6)
101 iform = nint(uparam(7))
102 dinit = nint(uparam(8))
103 dam_sf = uparam(9)
104 max_dam = uparam(10)
105 inst = nint(uparam(11))
106 n_val = uparam(12)
107 softexp = uparam(13)
108 ref_len = uparam(14)
109 reg_scale = uparam(15)
110 epfmin = uparam(16)
111c
112 ! Set flag for stress softening
113 IF (inst > 0) dmg_flag = 1
114c
115 ! At first timestep, initialization of the critical damage and
116 ! the element size scaling factor
117 IF (uvar(1,1) == zero) uvar(1:nel,1) = one
118 IF (uvar(1,3) == zero) THEN
119 IF (ifunc(1) > 0) THEN
120 DO i=1,nel
121 lambda = aldt(i)/ref_len
122 uvar(i,3) = finter(ifunc(1),lambda,npf,tf,dydx)
123 uvar(i,3) = uvar(i,3)*reg_scale
124 ENDDO
125 ELSE
126 uvar(1:nel,3) = one
127 ENDIF
128 ENDIF
129c
130 ! Checking element failure and recovering user variable
131 DO i=1,nel
132 ! Critical damage for necking instability
133 dc(i) = uvar(i,1)
134 ! Instability criterion
135 dinst(i) = uvar(i,2)
136 ! Recover element size scaling
137 fac(i) = uvar(i,3)
138 ENDDO
139c
140 !====================================================================
141 ! - DAMAGE INITIALIZATION
142 !====================================================================
143 IF (time == zero .AND. isigi /= 0 .AND. dinit > 0) THEN
144 DO i = 1,nel
145 IF (pla(i) > zero) THEN
146 ! Principal strains computation
147 e12 = half*epsxy(i)
148 s1 = half*(epsxx(i) + epsyy(i))
149 s2 = half*(epsxx(i) - epsyy(i))
150 q = sqrt(s2**2 + e12**2)
151 emaj = s1 + q
152 emin = s1 - q
153 ! Computation of the principal strains ratio BETA
154 beta = emin/max(emaj,em20)
155 beta = max(-one,beta)
156 beta = min( one,beta)
157 ! Estimation of stress triaxiality value
158 triax = (one/sqrt(three))*((one + beta)/sqrt(one + beta + beta**2))
159 IF (triax < -two_third) triax = -two_third
160 IF (triax > two_third) triax = two_third
161 ! Computation of Lode parameter
162 cos3theta = -half*twenty7*triax*(triax**2 - third)
163 IF (cos3theta < -one ) cos3theta = -one
164 IF (cos3theta > one ) cos3theta = one
165 lodep = one - two*acos(cos3theta)/pi
166 ! Computation of the plastic strain at failure
167 epsfail = c1 + c2*triax + c3*lodep + c4*(triax**2) +
168 . c5*(lodep**2) + c6*triax*lodep
169 epsfail = epsfail*fac(i)
170 epsfail = max(epfmin,epsfail)
171 ! Estimation of the damage variable value
172 dfmax(i) = pla(i)/max(epsfail,em6)
173 dfmax(i) = dfmax(i)*dam_sf
174 dfmax(i) = min(max_dam,dfmax(i))
175 ! Computation of instability curve and stress softening
176 IF (inst > 0 .AND. triax > em10) THEN
177 ! Computation of the principal strain ratio
178 alpha = (two*beta + one)/(two + beta)
179 ! Computation of the effective plastic strain
180 equa1 = one - alpha + alpha**2
181 equa2 = four - three*alpha - three*alpha**2 + four*alpha
182 e1 = (two*(two-alpha)*equa1/max(equa2,em20))*n_val
183 eps_c = e1*sqrt(four*third*(one + beta + beta**2))
184 eps_c = eps_c*fac(i)
185 ! Update instability criterion
186 dinst(i) = pla(i)/max(eps_c,em20)
187 dinst(i) = min(dinst(i),one)
188 uvar(i,2) = dinst(i)
189 ! Test if the instability curve is reached
190 IF (dinst(i) >= one) THEN
191 dc(i) = dfmax(i)
192 uvar(i,1) = dc(i)
193 ENDIF
194 ENDIF
195 ENDIF
196 ENDDO
197 ENDIF
198c
199 !====================================================================
200 ! - COMPUTATION OF THE DAMAGE VARIABLE EVOLUTION
201 !====================================================================
202 ! Initialization of element failure index
203 nindx = 0
204 indx(1:nel) = 0
205c
206 ! Loop over the elements
207 DO i=1,nel
208c
209 IF (foff(i) /= 0 .AND. dpla(i) > zero) THEN
210c
211 ! Computation of hydrostatic stress, Von Mises stress, and stress triaxiality
212 p = third*(signxx(i) + signyy(i))
213 ! Von Mises equivalent stress
214 svm = signxx(i)**2 + signyy(i)**2 - signxx(i)*signyy(i) +
215 . three*signxy(i)**2
216 svm = sqrt(max(svm,zero))
217 triax = p/max(em20,svm)
218 IF (triax < -two_third) triax = -two_third
219 IF (triax > two_third) triax = two_third
220c
221 ! Computation of Lode parameter
222 cos3theta = -half*twenty7*triax*(triax**2 - third)
223 IF (cos3theta < -one ) cos3theta = -one
224 IF (cos3theta > one ) cos3theta = one
225 lodep = one - two*acos(cos3theta)/pi
226c
227 ! Computation of the plastic strain at failure
228 epsfail = c1 + c2*triax + c3*lodep + c4*(triax**2) +
229 . c5*(lodep**2) + c6*triax*lodep
230 epsfail = epsfail*fac(i)
231 epsfail = max(epfmin,epsfail)
232c
233 ! Computation of the damage variable update
234 dfmax(i) = dfmax(i) + dpla(i)/max(epsfail,em20)
235 dfmax(i) = min(dfmax(i),one)
236c
237 ! Computation of instability curve and stress softening
238 IF (inst > 0) THEN
239 ! If the instability has not been reached
240 IF (dinst(i) < one .AND. triax > em10) THEN
241 ! Computation of the principal strains ratio BETA from stress triaxiality
242 delta = three*(triax**2)*(four - nine*(triax**2))
243 denom = two*(three*(triax**2)-one)
244 denom = sign(max(abs(denom),em20),denom)
245 beta = ((two - three*(triax**2))-sqrt(delta))/denom
246 beta = max(-one,beta)
247 beta = min( one,beta)
248 ! Computation of the principal stress ratio
249 alpha = (two*beta + one)/(two + beta)
250 ! Computation of the effective plastic strain
251 equa1 = one - alpha + alpha**2
252 equa2 = four - three*alpha - three*alpha**2 + four*alpha
253 e1 = (two*(two-alpha)*equa1/max(equa2,em20))*n_val
254 eps_c = e1*sqrt(four*third*(one + beta + beta**2))
255 eps_c = eps_c*fac(i)
256 ! Update instability criterion
257 IF (iform == 2) THEN
258 dinst(i) = max(dinst(i),pla(i)/max(eps_c,em20))
259 ELSE
260 dinst(i) = dinst(i) + dpla(i)/max(eps_c,em20)
261 ENDIF
262 dinst(i) = min(dinst(i),one)
263 uvar(i,2) = dinst(i)
264 ! Test if the instability curve is reached
265 IF (dinst(i) >= one) THEN
266 dc(i) = dfmax(i)
267 uvar(i,1) = dc(i)
268 ENDIF
269 ENDIF
270 ENDIF
271c
272 ! Check element failure
273 IF (dfmax(i) >= one) THEN
274 foff(i) = 0
275 nindx = nindx + 1
276 indx(nindx) = i
277 ENDIF
278 ENDIF
279 ENDDO
280c
281 !====================================================================
282 ! - UPDATE DAMAGE SCALE FACTOR FOR NECKING INSTABILITY
283 !====================================================================
284 IF (inst > 0) THEN
285 DO i = 1,nel
286 ! Computation of the damage scale factor
287 IF (dfmax(i) >= dc(i)) THEN
288 IF (dc(i) < one) THEN
289 dmg_scale(i) = one - ((dfmax(i)-dc(i))/max(one-dc(i),em20))**softexp
290 ELSE
291 dmg_scale(i) = zero
292 ENDIF
293 ELSE
294 dmg_scale(i) = one
295 ENDIF
296 ENDDO
297 ENDIF
298c
299 !====================================================================
300 ! - PRINTOUT DATA ABOUT FAILED ELEMENTS
301 !====================================================================
302 IF (nindx > 0) THEN
303 DO j = 1,nindx
304 i = indx(j)
305#include "lockon.inc"
306 WRITE(iout, 1000) ngl(i),ipg,ipt,time
307 WRITE(istdo,1000) ngl(i),ipg,ipt,time
308#include "lockoff.inc"
309 ENDDO
310 ENDIF
311C---------
312 1000 FORMAT(1x,'FOR SHELL ELEMENT NUMBER el#',i10,
313 . ' FAILURE (SYAZWAN) AT GAUSS POINT ',i3,' LAYER ',i3,
314 . ' AT TIME :',1pe12.4)
315c---------
316 END
alpha
#define alpha
Definition eval.h:35
fail_syazwan_c
subroutine fail_syazwan_c(nel, uparam, nuparam, uvar, nuvar, time, ngl, ipt, dpla, pla, signxx, signyy, signxy, signyz, signzx, epsxx, epsyy, epsxy, epsyz, epszx, dfmax, nfunc, ifunc, aldt, foff, ipg, dmg_flag, dmg_scale, npf, tf)
Definition fail_syazwan_c.F:38
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21