OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_hashin_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_hashin_c ../engine/source/materials/fail/hashin/fail_hashin_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!||====================================================================
29 SUBROUTINE fail_hashin_c(
30 1 NEL ,NUPARAM ,NUVAR ,UPARAM ,UVAR ,
31 2 TIME ,TIMESTEP ,IPG ,ILAY ,IPT ,
32 3 NGL ,DMG_FLAG ,DMG_SCALE ,DFMAX ,TDEL ,
33 4 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
34 5 OFF ,FOFF ,PLY_ID ,
35 6 EPSP ,FWAVE_EL ,DADV ,LF_DAMMX )
36C-----------------------------------------------
37c Hashin failure model
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 "units_c.inc"
46#include "comlock.inc"
47C---------+---------+---+---+--------------------------------------------
48C VAR | SIZE |TYP| RW| DEFINITION
49C---------+---------+---+---+--------------------------------------------
50C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
51C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
52C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
53C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
54C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
55C---------+---------+---+---+--------------------------------------------
56C TIME | 1 | F | R | CURRENT TIME
57C---------+---------+---+---+--------------------------------------------
58C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
59C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
60C ... | | | |
61C---------+---------+---+---+--------------------------------------------
62C OFF | NEL | F | R | DELETED ELEMENT FLAG (=1. ON, =0. OFF)
63C FOFF | NEL | I |R/W| DELETED INTEGRATION POINT FLAG (=1 ON, =0 OFF)
64C DFMAX | NEL | F |R/W| MAX DAMAGE FACTOR
65C TDEL | NEL | F | W | FAILURE TIME
66C DMG_FLAG| 1 | I | W | STRESS REDUCTION FLAG DUE TO DAMAGE
67C DMG_SCALE| NEL | F | W | STRESS REDUCTION FACTOR
68C---------+---------+---+---+--------------------------------------------
69C NGL ELEMENT ID
70C IPG CURRENT GAUSS POINT (in plane)
71C ILAY CURRENT LAYER
72C IPT CURRENT INTEGRATION POINT IN THE LAYER
73C---------+---------+---+---+--------------------------------------------
74C I N P U T A r g u m e n t s
75C-----------------------------------------------
76 INTEGER ,INTENT(IN) :: NEL,NUPARAM,NUVAR,IPG,ILAY,IPT,PLY_ID,LF_DAMMX
77 INTEGER ,DIMENSION(NEL) ,INTENT(IN) :: NGL
78 my_real ,INTENT(IN) :: TIME
79 my_real ,DIMENSION(NEL) ,INTENT(IN) :: OFF,TIMESTEP
80 my_real,DIMENSION(NUPARAM) ,INTENT(IN) :: UPARAM
81C-----------------------------------------------
82C I N P U T O U T P U T A r g u m e n t s
83C-----------------------------------------------
84 INTEGER ,INTENT(OUT) ::DMG_FLAG
85 INTEGER ,DIMENSION(*) :: FWAVE_EL
86 INTEGER ,DIMENSION(NEL) ,INTENT(INOUT) :: FOFF
87 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: DADV,
88 . SIGNXX,SIGNYY,SIGNXY,SIGNYZ,SIGNZX,EPSP
89 my_real ,DIMENSION(NEL,LF_DAMMX) ,INTENT(INOUT) :: DFMAX
90 my_real ,DIMENSION(NEL) ,INTENT(OUT) :: tdel,dmg_scale
91 my_real ,DIMENSION(NEL,NUVAR) ,INTENT(INOUT) :: uvar
92C-----------------------------------------------
93C L o c a l V a r i a b l e s
94C-----------------------------------------------
95 INTEGER :: I,J,JJ,NINDX,IFAIL_SH,IMODEL,IMODE,TMOD,IFRWAVE,FAC
96 INTEGER ,DIMENSION(NEL) :: INDX,FS
97 INTEGER ,DIMENSION(7,NEL) :: MODE
98 my_real :: SIGT1,SIGT2,SIGC1,SIGC2,CSIG,FSIG12,MSIG12,MSIG23,
99 . MSIG13,ANGLE,SDEL,TMAX,RATIO,F1,F2,F3,F4,F5,SIG,P,F6,F7,
100 . TSIG12,XSIG12,XSIG23,DAMMX, TELEM,K2M,INSTF,K0,K1,K2,K3,
101 . zep669741,dtinv,km,epspref,filt,alpha
102c=======================================================================C
103 zep669741 = six*em01 + six*em02 + nine*em3 + seven*em04 +
104 . four*em05 + em06
105C
106 imodel = int(uparam(1))
107 sigt1 = uparam(2)
108 sigt2 = uparam(3)
109 sigc1 = uparam(5)
110 sigc2 = uparam(6)
111 csig = uparam(7)
112 fsig12 = uparam(8)
113 msig12 = uparam(9)
114 msig13 = uparam(10)
115 msig23 = uparam(11)
116 angle = uparam(12)
117 sdel = uparam(13)
118 tmax = uparam(14)
119 ifail_sh = int(uparam(15))
120 ratio = uparam(17)
121 dmg_flag = int(uparam(18))
122 tmod = int(uparam(19))
123 epspref = uparam(20)
124 filt = uparam(21)
125 km = uparam(22)
126 ifrwave = nint(uparam(23))
127c-------------------------------
128 mode(1,1:nel) = 0
129 mode(2,1:nel) = 0
130 mode(3,1:nel) = 0
131 mode(4,1:nel) = 0
132 mode(5,1:nel) = 0
133 mode(6,1:nel) = 0
134 mode(7,1:nel) = 0
135 nindx = 0
136c
137 IF (time == zero) then
138 dadv(:) = one
139 ENDIF
140c
141 IF (ifrwave > 1) THEN ! front wave propagation
142 DO i=1,nel
143 IF (fwave_el(i) /= 0) dadv(i) = uparam(24)
144 END DO
145 ENDIF
146c-------------------------------
147 SELECT CASE (imodel)
148c-------------------------------
149 CASE (1) ! Matrix or fiber failure mode
150c-------------------------------
151 DO i=1,nel
152 IF (off(i) == one .and. foff(i) == 1) THEN
153 f1 = zero
154 f2 = zero
155 f3 = zero
156 f4 = zero
157 f5 = zero
158C alpha = number of cycles
159 dtinv = timestep(i)/max(timestep(i)**2,em20)
160 alpha = int(max(one,filt*dtinv))
161c alpha coef de moyenne mobile exponentielle
162 alpha = two/(alpha + one)
163 fs(i)=uvar(i,10)
164c smoothed strain rate
165 epsp(i)=(one - alpha)*uvar(i,11) + alpha*epsp(i)
166 uvar(i,11)=epsp(i)
167 IF (uvar(i,6) == zero) THEN
168c
169 ! Fiber failure criteria
170 ! -> Tension
171 sig = half*(signxx(i) + abs(signxx(i)))
172 f1 = (sig/sigt1)**2
173 . + (signxy(i)**2 + signzx(i)**2)/fsig12**2
174 dfmax(i,2) = max(dfmax(i,2),f1)
175 dfmax(i,2) = min(dfmax(i,2),one)
176 ! -> Compression
177 sig = -half*signyy(i)
178 sig = -signxx(i) + half*(sig + abs(sig))
179 sig = half*(sig + abs(sig))
180 f2 = (sig /sigc1)**2
181 dfmax(i,3) = max(dfmax(i,3),f2)
182 dfmax(i,3) = min(dfmax(i,3),one)
183c
184 ! Crush mode criterion
185 p = -third*(signxx(i) + signyy(i))
186 IF (p > 0) f3 = (p/csig)**2
187 dfmax(i,4) = max(dfmax(i,4),f3)
188 dfmax(i,4) = min(dfmax(i,4),one)
189c
190 ! Matrix failure criteria
191 ! -> Delamination mode
192 f5 = (signyz(i)/msig23)**2 + (signzx(i)/msig13)**2
193 f5 = sdel*sdel*f5
194 dfmax(i,6) = max(dfmax(i,6),f5)
195 dfmax(i,6) = min(dfmax(i,6),one)
196c
197 ! -> Failure matrix mode
198 IF (signyy(i) < zero) THEN
199 xsig12 = msig12 - signyy(i)*(tan(angle))
200 xsig23 = msig23 - signyy(i)*(tan(angle))
201 ELSE
202 xsig12 = msig12
203 xsig23 = msig23
204 ENDIF
205 sig = half*(signyy(i) + abs(signyy(i)))
206 f4 = (sig/sigt2)**2
207 . + (signyz(i)/xsig23)**2 + (signxy(i)/xsig12)**2
208 dfmax(i,5) = max(dfmax(i,5),f4)
209 dfmax(i,5) = min(dfmax(i,5),one)
210c
211 dammx = min(one,max(f1,f2,f3,f4,f5))
212 dfmax(i,1) = min(one,dammx)
213 IF (dammx >= dadv(i)) THEN
214 nindx = nindx+1
215 indx(nindx) = i
216 uvar(i,6) = time
217C
218 fac = hundred*max(abs(signxx(i)),abs(signyy(i)),
219 . abs(signxy(i)),abs(signyz(i)),abs(signzx(i)))
220 k0 = fac*epspref**2*zep669741*tmax**2 + zep9*epspref*tmax
221 k1 = zep9*epsp(i)
222 k2 = zep669741*fac*epsp(i)**2
223 k2m = zep669741*fac*km*epspref**2
224 k2 = max(k2m,k2,em20)
225 telem = (sqrt(k1**2+4*k2*k0)-k1)/2/k2
226 uvar(i,12) = telem
227 IF (f1>=dadv(i)) mode(1,i) = 1
228 IF (f2>=dadv(i)) mode(2,i) = 1
229 IF (f3>=dadv(i)) mode(3,i) = 1
230 IF (f4>=dadv(i)) mode(4,i) = 1
231 IF (f5>=dadv(i)) mode(5,i) = 1
232 IF (dmg_flag == 0) THEN
233 uvar(i,1) = signxx(i)
234 uvar(i,2) = signyy(i)
235 uvar(i,3) = signxy(i)
236 uvar(i,4) = signyz(i)
237 uvar(i,5) = signzx(i)
238 ENDIF
239 ENDIF
240 ENDIF
241c
242 IF (uvar(i,6) > zero .AND. time > uvar(i,6) ) THEN
243 IF(tmod == 1) THEN
244 telem = max((timestep(i)*1.0),tmax*(epspref/epsp(i))**2)
245 dmg_scale(i)= dmg_scale(i)*exp(-timestep(i)/telem)
246 ELSEIF(tmod == 2) THEN
247 telem = uvar(i, 12)
248 dmg_scale(i) = exp(-(time - uvar(i,6))/telem)
249 ELSE
250 dmg_scale(i) = exp(-(time - uvar(i,6))/tmax)
251 ENDIF
252 IF (dmg_scale(i) < em02) THEN
253 foff(i) = 0
254 tdel(i)= time
255 dmg_scale(i) = zero
256 ENDIF
257 IF (dmg_flag == 0) THEN
258 signxx(i) = uvar(i,1)*dmg_scale(i)
259 signyy(i) = uvar(i,2)*dmg_scale(i)
260 signxy(i) = uvar(i,3)*dmg_scale(i)
261 signyz(i) = uvar(i,4)*dmg_scale(i)
262 signzx(i) = uvar(i,5)*dmg_scale(i)
263 ENDIF
264 ENDIF
265 ENDIF
266 ENDDO
267c ----------------------------------------------
268 CASE (2) ! Fabric failure mode
269c-------------------------------
270 DO i=1,nel
271 IF (off(i) == one .and. foff(i) == 1) THEN
272 f1 = zero
273 f2 = zero
274 f3 = zero
275 f4 = zero
276 f5 = zero
277 f6 = zero
278 f7 = zero
279c
280C alpha = number of cycles
281 dtinv = timestep(i)/max(timestep(i)**2,em20)
282 alpha = int(max(one,filt*dtinv))
283c alpha coef de moyenne mobile exponentielle
284 alpha = two/(alpha + one)
285 fs(i)=uvar(i,10)
286c smoothed strain rate
287 epsp(i)=(one - alpha)*uvar(i,11) + alpha*epsp(i)
288 uvar(i,11)=epsp(i)
289 IF (uvar(i,6) == zero) THEN
290c
291 ! Fiber failure criteria
292 ! -> Tension
293 sig = half*(signxx(i) + abs(signxx(i)))
294 f1 = (sig/sigt1)**2
295 . + (signxy(i)**2 + signzx(i)**2)/fsig12**2
296 dfmax(i,2) = max(dfmax(i,2),f1)
297 dfmax(i,2) = min(dfmax(i,2),one)
298c
299 ! -> Shear
300 sig = half*(signyy(i) + abs(signyy(i)))
301 tsig12 = fsig12*sigt2/sigt1
302 f2 = (sig/sigt2)**2
303 . + (signxy(i)**2 + signyz(i)**2)/tsig12**2
304 dfmax(i,3) = max(dfmax(i,3),f2)
305 dfmax(i,3) = min(dfmax(i,3),one)
306c
307 ! Matrix compression failure criteria
308 IF (signxx(i) < zero) f3 = (signxx(i)/sigc1)**2
309 dfmax(i,4) = max(dfmax(i,4),f3)
310 dfmax(i,4) = min(dfmax(i,4),one)
311 IF (signyy(i) < zero) f4 = (signyy(i)/sigc2)**2
312 dfmax(i,5) = max(dfmax(i,5),f4)
313 dfmax(i,5) = min(dfmax(i,5),one)
314c
315 ! -> Crush mode criterion
316 p = -third*(signxx(i) + signyy(i) )
317 IF (p > zero) f5 = (p/csig)**2
318 dfmax(i,6) = max(dfmax(i,6),f5)
319 dfmax(i,6) = min(dfmax(i,6),one)
320c
321 ! Shear Matrix failure
322 f6 = (signxy(i)/msig12)**2
323 dfmax(i,7) = max(dfmax(i,7),f6)
324 dfmax(i,7) = min(dfmax(i,7),one)
325c
326 ! Matrix failure
327 f7 = (signyz(i)/msig23)**2 + (signzx(i)/msig13)**2
328 f7 = sdel*sdel*f7
329 dfmax(i,8) = max(dfmax(i,8),f7)
330 dfmax(i,8) = min(dfmax(i,8),one)
331c
332 dammx = min(one,max(f1,f2,f3,f4,f5,f6,f7))
333 dfmax(i,1) = min(one,dammx)
334 IF (dammx >= dadv(i)) THEN
335 nindx = nindx+1
336 indx(nindx) = i
337 uvar(i,6) = time
338C
339 fac = hundred*max(abs(signxx(i)),abs(signyy(i)),
340 . abs(signxy(i)),abs(signyz(i)),abs(signzx(i)))
341 k0 = fac*epspref**2*zep669741*tmax**2 + zep9*epspref*tmax
342 k1 = zep9*epsp(i)
343 k2 = zep669741*fac*epsp(i)**2
344 k2m = zep669741*fac*km*epspref**2
345 k2 = max(k2m,k2,em20)
346 telem = (sqrt(k1**2+4*k2*k0)-k1)/2/k2
347 uvar(i,12) = telem
348C
349 IF (f1>=dadv(i)) mode(1,i) = 1
350 IF (f2>=dadv(i)) mode(2,i) = 1
351 IF (f3>=dadv(i)) mode(3,i) = 1
352 IF (f4>=dadv(i)) mode(4,i) = 1
353 IF (f5>=dadv(i)) mode(5,i) = 1
354 IF (f6>=dadv(i)) mode(6,i) = 1
355 IF (f7>=dadv(i)) mode(7,i) = 1
356 IF (dmg_flag == 0) THEN
357 uvar(i,1) = signxx(i)
358 uvar(i,2) = signyy(i)
359 uvar(i,3) = signxy(i)
360 uvar(i,4) = signyz(i)
361 uvar(i,5) = signzx(i)
362 ENDIF
363 ENDIF
364 ENDIF
365c
366 IF (uvar(i,6) > zero .AND. time > uvar(i,6) ) THEN
367 dfmax(i,1) = one
368 IF(tmod == 1) THEN
369 telem = max((timestep(i)*1.0),tmax*(epspref/epsp(i))**2)
370 dmg_scale(i)=dmg_scale(i)*exp(-timestep(i)/telem)
371 ELSEIF(tmod == 2) THEN
372 telem = uvar(i,12)
373 dmg_scale(i)=exp(-(time - uvar(i,6))/telem)
374 ELSE
375 dmg_scale(i)= exp(-(time - uvar(i,6))/tmax)
376 ENDIF
377c
378 IF (dmg_scale(i) < em02) THEN
379 foff(i) = 0
380 tdel(i) = time
381 dmg_scale(i) = zero
382 ENDIF
383 IF (dmg_flag == 0) THEN
384 signxx(i) = uvar(i,1)*dmg_scale(i)
385 signyy(i) = uvar(i,2)*dmg_scale(i)
386 signxy(i) = uvar(i,3)*dmg_scale(i)
387 signyz(i) = uvar(i,4)*dmg_scale(i)
388 signzx(i) = uvar(i,5)*dmg_scale(i)
389 ENDIF
390 ENDIF
391 ELSEIF (foff(i) == 0) THEN
392 dmg_scale(i) = zero
393 ENDIF ! FOFF
394 ENDDO ! I=1,NEL
395c-------------
396 END SELECT
397c--------------------------------------------
398c- print
399c--------------------------------------------
400 IF (nindx > 0) THEN
401 DO j=1,nindx
402 i = indx(j)
403 DO jj=1,6
404 imode = mode(jj,i)
405 IF (imode == 1) THEN
406 imode = jj
407#include "lockon.inc"
408 WRITE(iout ,1000) ngl(i),ipg,ilay,ipt
409 WRITE(istdo ,1000) ngl(i),ipg,ilay,ipt
410 WRITE(iout ,2000) imode,time
411 WRITE(istdo ,2000) imode,time
412 WRITE(iout ,3000) signxx(i),signyy(i),signxy(i),
413 . signyz(i),signzx(i)
414 IF (dadv(i) < one) WRITE(iout,4000) dadv(i)
415#include "lockoff.inc"
416 ENDIF
417 ENDDO
418 END DO
419 ENDIF
420c--------------------------------------------
421 1000 FORMAT(1x,'FAILURE (HASHIN) OF SHELL ELEMENT ',i10,
422 . 1x,',GAUSS PT',i2,1x,',LAYER',i3,1x,',INTEGRATION PT',i3)
423 2000 FORMAT(10x,'HASHIN MODE #',i10,1x, ',time #:',1PE12.4)
424 3000 FORMAT(10x,'RESPONSIBLE STRESS: SIG_11= ',1pe12.4,
425 . 1x,'SIG_22= ',1pe12.4,1x,'SIG_12= ',1pe12.4,/,
426 . 10x,'SIG_13= ',1pe12.4,1x,'SIG_23= ',1pe12.4,1x)
427 4000 FORMAT(1x,'CRITERIA REDUCTION FACTOR= ',1pe12.4)
428c--------------------------------------------
429 RETURN
430 END
alpha
#define alpha
Definition eval.h:35
fail_hashin_c
subroutine fail_hashin_c(nel, nuparam, nuvar, uparam, uvar, time, timestep, ipg, ilay, ipt, ngl, dmg_flag, dmg_scale, dfmax, tdel, signxx, signyy, signxy, signyz, signzx, off, foff, ply_id, epsp, fwave_el, dadv, lf_dammx)
Definition fail_hashin_c.F:36
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21