OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_puck_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/.
23C modified Puck model ------
24!||====================================================================
25!|| fail_puck_s ../engine/source/materials/fail/puck/fail_puck_s.F
26!||--- called by ------------------------------------------------------
27!|| mmain ../engine/source/materials/mat_share/mmain.F90
28!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
29!|| usermat_solid ../engine/source/materials/mat_share/usermat_solid.F
30!||====================================================================
31 SUBROUTINE fail_puck_s(
32 1 NEL ,NUVAR ,ILAY ,NPT0 ,
33 2 TIME ,TIMESTEP,UPARAM ,
34 3 NGL ,OFF ,NOFF ,SIGNXX ,
35 4 SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
36 5 UVAR ,NUPARAM ,DFMAX ,LF_DAMMX,TDELE )
37C-----------------------------------------------
38C I m p l i c i t T y p e s
39C-----------------------------------------------
40#include "implicit_f.inc"
41C---------+---------+---+---+--------------------------------------------
42C VAR | SIZE |TYP| RW| DEFINITION
43C---------+---------+---+---+--------------------------------------------
44C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
45C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
46C---------+---------+---+---+--------------------------------------------
47C TIME | 1 | F | R | CURRENT TIME
48C TIMESTEP| 1 | F | R | CURRENT TIME STEP
49C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
50C EPSPXX | NEL | F | R | STRAIN RATE XX
51C EPSPYY | NEL | F | R | STRAIN RATE YY
52C ... | | | |
53C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
54C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
55C ... | | | |
56C EPSXX | NEL | F | R | STRAIN XX
57C EPSYY | NEL | F | R | STRAIN YY
58C ... | | | |
59C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
60C SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
61C ... | | | |
62C---------+---------+---+---+--------------------------------------------
63C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
64C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
65C ... | | | |
66C SIGVXX | NEL | F | W | VISCOUS STRESS XX
67C SIGVYY | NEL | F | W | VISCOUS STRESS YY
68C ... | | | |
69C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
70C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
71C---------+---------+---+---+--------------------------------------------
72C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
73C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
74C---------+---------+---+---+--------------------------------------------
75#include "mvsiz_p.inc"
76#include "com01_c.inc"
77#include "units_c.inc"
78#include "comlock.inc"
79#include "param_c.inc"
80C-----------------------------------------------
81C I N P U T A r g u m e n t s
82C-----------------------------------------------
83 INTEGER NEL,NUPARAM, NUVAR,ILAY,NPT0
84 INTEGER NGL(NEL)
85 INTEGER, INTENT(IN) :: LF_DAMMX
86 my_real
87 . TIME,TIMESTEP(NEL),UPARAM(NUPARAM),
88 . signxx(nel),signyy(nel),signzz(nel),
89 . signxy(nel),signyz(nel),signzx(nel)
90C-----------------------------------------------
91C I N P U T O U T P U T A r g u m e n t s
92C-----------------------------------------------
93 INTEGER NOFF(NEL)
94 my_real UVAR(NEL,NUVAR), OFF(NEL), TDELE(NEL)
95 my_real ,DIMENSION(NEL,LF_DAMMX),INTENT(INOUT) :: DFMAX
96C-----------------------------------------------
97C L o c a l V a r i a b l e s
98C-----------------------------------------------
99 INTEGER
100 . i,j,idel,idel_l,iflag,indx(mvsiz),iadbuf,nindx,
101 . nindex,index(mvsiz),ifail,jj,imatly,
102 . iunidir,ifabric,indx0(mvsiz),nindx0,fsmooth
103 my_real
104 . sigt1,sigt2, sigc1,sigc2,
105 . fsig12,f1,fa,fb,fc,pn12,pp12,pn22,
106 . fac,tmax,fa_2,fb_2,fc_2,fcut,asrate,
107 . sxx(nel),syy(nel),szz(nel),
108 . sxy(nel),syz(nel),szx(nel)
109C--------------------------------------------------------------
110 sigt1 = uparam(1 )
111 sigt2 = uparam(2)
112 sigc1 = uparam(3)
113 sigc2 = uparam(4)
114 fsig12 = uparam(5)
115 pp12 = uparam(6)
116 pn12 = uparam(7)
117 pn22 = uparam(8)
118 tmax = uparam(9)
119 iflag = int(uparam(11))
120 fcut = uparam(12)
121 IF (fcut > zero) THEN
122 fsmooth = 1
123 asrate = two*pi*fcut*timestep(1)
124 asrate = asrate/(one+asrate)
125 ELSE
126 fsmooth = 0
127 ENDIF
128 indx(1:mvsiz) = 0
129 index(1:mvsiz) = 0
130C-----------------------------------------------
131C USER VARIABLES INITIALIZATION
132C-----------------------------------------------
133 IF(isigi == 0)THEN
134 IF ((uvar(1,8) == zero).AND.(off(1) == one)) THEN
135 DO i=1,nel
136 uvar(i,8) = one
137 ENDDO
138 ENDIF
139 ENDIF
140C-----------------------------------------------
141 iunidir = 0
142 ifabric = 0
143 DO i=1,nel
144 IF(off(i)<em01) off(i)=zero
145 IF(off(i)<one) off(i)=off(i)*four_over_5
146 IF (fsmooth > 0) THEN
147 sxx(i) = asrate*signxx(i) + (one - asrate)*uvar(i,9)
148 syy(i) = asrate*signyy(i) + (one - asrate)*uvar(i,10)
149 szz(i) = asrate*signzz(i) + (one - asrate)*uvar(i,11)
150 sxy(i) = asrate*signxy(i) + (one - asrate)*uvar(i,12)
151 syz(i) = asrate*signyz(i) + (one - asrate)*uvar(i,13)
152 szx(i) = asrate*signzx(i) + (one - asrate)*uvar(i,14)
153 uvar(i,9) = sxx(i)
154 uvar(i,10) = syy(i)
155 uvar(i,11) = szz(i)
156 uvar(i,12) = sxy(i)
157 uvar(i,13) = syz(i)
158 uvar(i,14) = szx(i)
159 ELSE
160 sxx(i) = signxx(i)
161 syy(i) = signyy(i)
162 szz(i) = signzz(i)
163 sxy(i) = signxy(i)
164 syz(i) = signyz(i)
165 szx(i) = signzx(i)
166 ENDIF
167 END DO
168C-------------------------------
169C
170C OFF = 0. si la matrice ou la fibre a rompu
171C-------------------------------
172 nindx = 0
173 nindx0= 0
174 indx = 0
175 indx0 = 0
176 DO i=1,nel
177 f1 = zero
178 fa = zero
179 fb = zero
180 fc = zero
181 fa_2 = zero
182 fb_2 = zero
183 fc_2 = zero
184C
185 IF (off(i) == one )THEN
186C
187 IF (iflag == 1) THEN
188C-------------------------------
189C OFF = 0 when one layer fiber or matrix criteria is reached
190C-------------------------------
191 IF (uvar(i,8) < one) THEN
192 uvar(i,8) = exp(-(time - uvar(i,7))/tmax)
193 IF (uvar(i,8) < em02) uvar(i,8) = zero
194 signxx(i) = uvar(i,1)*uvar(i,8)
195 signyy(i) = uvar(i,2)*uvar(i,8)
196 signzz(i) = uvar(i,3)*uvar(i,8)
197 signxy(i) = uvar(i,4)*uvar(i,8)
198 signyz(i) = uvar(i,5)*uvar(i,8)
199 signzx(i) = uvar(i,6)*uvar(i,8)
200 IF (uvar(i,8) == zero )THEN
201 off(i)=four_over_5
202 nindx=nindx+1
203 indx(nindx)=i
204 tdele(i) = time
205 ENDIF
206 ELSE
207C
208C fiber criteria
209C
210 IF (sxx(i) >= zero) THEN
211 f1 = sxx(i)/sigt1
212 dfmax(i,2) = max(dfmax(i,2),f1)
213 dfmax(i,2) = min(dfmax(i,2),one)
214 ELSE
215 f1 = -sxx(i)/sigc1
216 dfmax(i,3) = max(dfmax(i,3),f1)
217 dfmax(i,3) = min(dfmax(i,3),one)
218 ENDIF
219C
220C matrice criteria direction 2
221C
222 IF (syy(i) >= zero) THEN
223 fac = one - pp12*sigt2/fsig12
224 fac = fac*syy(i)/sigt2
225 fa = sqrt((sxy(i)/fsig12)**2 + fac*fac)
226 . + pp12*syy(i)/fsig12
227 dfmax(i,4) = max(dfmax(i,4),fa)
228 dfmax(i,4) = min(dfmax(i,4),one)
229 ELSE
230 fac = half/(one + pn22)/fsig12
231 fc = (sxy(i)*fac)**2 + (syy(i)/sigc2)**2
232 fc =-fc*sigc2/syy(i)
233 dfmax(i,6) = max(dfmax(i,6),fc)
234 dfmax(i,6) = min(dfmax(i,6),one)
235C plane angle
236c FAC = SQRT(ONE + TWO*PN12(I)*SIGC2(I)/FSIG12(I)) - ONE
237c FAC = HALF*FAC*FSIG12(I)/PN12(I)
238c CC = (FAC/FSIG12(I))**2
239c CC = CC*(SIGNXY(I)/SIGNYY(I))**2 + ONE
240c FAC = HALF/ONE + PN22(I))
241c CC = SQRT(FAC*CC)
242c ANGLE(I) = ACOS(CC)
243 ENDIF
244
245 fb = sqrt(sxy(i)**2 + (pn12*syy(i))**2 )
246 . + pn12*syy(i)
247 fb = fb/fsig12
248 dfmax(i,5) = max(dfmax(i,5),fb)
249 dfmax(i,5) = min(dfmax(i,5),one)
250C
251C matrice criteria direction 3
252C
253 IF (szz(i) >= zero) THEN
254 fac = one - pp12*sigt2/fsig12
255 fac = fac*szz(i)/sigt2
256 fa_2 = sqrt((szx(i)/fsig12)**2 + fac*fac)
257 . + pp12*szz(i)/fsig12
258 dfmax(i,4) = max(dfmax(i,4),fa_2)
259 dfmax(i,4) = min(dfmax(i,4),one)
260 ELSE
261 fac = half/(one + pn22)/fsig12
262 fc_2 = (szx(i)*fac)**2 + (szz(i)/sigc2)**2
263 fc_2 =-fc_2*sigc2/szz(i)
264 dfmax(i,6) = max(dfmax(i,6),fc_2)
265 dfmax(i,6) = min(dfmax(i,6),one)
266C plane angle
267c FAC = SQRT(ONE + TWO*PN12(I)*SIGC2(I)/FSIG12(I)) - ONE
268c FAC = HALF*FAC*FSIG12(I)/PN12(I)
269c CC = (FAC/FSIG12(I))**2
270c CC = CC*(SIGNZX(I)/SIGNZZ(I))**2 + ONE
271c FAC = HALF/ONE + PN22(I))
272c CC = SQRT(FAC*CC)
273c ANGLE(I) = ACOS(CC)
274 ENDIF
275
276 fb_2 = sqrt(szx(i)**2 + (pn12*szz(i))**2 )
277 . + pn12*szz(i)
278 fb_2 = fb_2/fsig12
279 dfmax(i,5) = max(dfmax(i,5),fb_2)
280 dfmax(i,5) = min(dfmax(i,5),one)
281C
282 dfmax(i,1) = min(one,max(dfmax(i,1),f1,fa,fb,fc,fa_2,fb_2,fc_2))
283C
284 IF (f1 >= one .OR. fa >= one .OR. fb >= one .OR.
285 . fc >= one .OR. fa_2 >= one .OR. fb_2 >= one .OR.
286 . fc_2 >= one ) THEN
287 uvar(i,1) = signxx(i)
288 uvar(i,2) = signyy(i)
289 uvar(i,3) = signzz(i)
290 uvar(i,4) = signxy(i)
291 uvar(i,5) = signyz(i)
292 uvar(i,6) = signzx(i)
293 uvar(i,7) = time
294 uvar(i,8) = four_over_5
295 ENDIF
296 ENDIF
297C
298 ELSE ! iflag/= 1
299C-------------------------------
300C OFF = 0. all layer fiber or matrix criteria is reateched
301C-------------------------------
302 IF(uvar(i,8) == zero )THEN
303 signxx(i) = zero
304 signyy(i) = zero
305 signzz(i) = zero
306 signxy(i) = zero
307 signzx(i) = zero
308 signyz(i) = zero
309 ELSE IF(uvar(i,8) < one) THEN
310 uvar(i,8)= exp(-(time - uvar(i,7))/tmax)
311 IF(uvar(i,8) < em02)uvar(i,8) = zero
312 signxx(i) = uvar(i,1)*uvar(i,8)
313 signyy(i) = uvar(i,2)*uvar(i,8)
314 signzz(i) = uvar(i,3)*uvar(i,8)
315 signxy(i) = uvar(i,4)*uvar(i,8)
316 signyz(i) = uvar(i,5)*uvar(i,8)
317 signzx(i) = uvar(i,6)*uvar(i,8)
318 IF(uvar(i,8) == zero )THEN
319 noff(i) = noff(i) + 1
320 IF(npt0 == 1 .OR. noff(i) == npt0 ) THEN
321 off(i) = four_over_5
322 nindx=nindx+1
323 indx(nindx)=i
324 tdele(i) = time
325 ENDIF
326 ENDIF
327 ELSE
328C
329C fiber criteria
330C
331 IF(sxx(i) > 0 ) THEN
332 f1 = sxx(i)/sigt1
333 dfmax(i,2) = max(dfmax(i,2),f1)
334 dfmax(i,2) = min(dfmax(i,2),one)
335 ELSE
336 f1 = -sxx(i)/sigc1
337 dfmax(i,3) = max(dfmax(i,3),f1)
338 dfmax(i,3) = min(dfmax(i,3),one)
339 ENDIF
340C
341C matrice criteria --direction 2
342C
343 IF(syy(i) >= zero) THEN
344 fac = one - pp12*sigt2/fsig12
345 fac = fac*syy(i)/sigt2
346 fa = sqrt((sxy(i)/fsig12)**2 + fac*fac)
347 . + pp12*syy(i)/fsig12
348 dfmax(i,4) = max(dfmax(i,4),fa)
349 dfmax(i,4) = min(dfmax(i,4),one)
350 ELSE
351 fac = half/(one + pn22)/fsig12
352 fc = (sxy(i)*fac)**2 + (syy(i)/sigc2)**2
353 fc =-fc*sigc2/syy(i)
354 dfmax(i,6) = max(dfmax(i,6),fc)
355 dfmax(i,6) = min(dfmax(i,6),one)
356C plane angle
357c FAC = SQRT(ONE + TWO*PN12(I)*SIGC2(I)/FSIG12(I)) - ONE
358c FAC = HALF*FAC*FSIG12(I)/PN12(I)
359c CC = (FAC/FSIG12(I))**2
360c CC = CC*(SIGNXY(I)/SIGNYY(I))**2 + ONE
361c FAC = HALF/ONE + PN22(I))
362c CC = SQRT(FAC*CC)
363c ANGLE(I) = ACOS(CC)
364 ENDIF
365 fb = sqrt(sxy(i)**2 + (pn12*syy(i))**2 )
366 . + pn12*syy(i)
367 fb = fb/fsig12
368 dfmax(i,5) = max(dfmax(i,5),fb)
369 dfmax(i,5) = min(dfmax(i,5),one)
370C
371C matrice criteria direction 3
372C
373 IF(szz(i) >= zero) THEN
374 fac = one - pp12*sigt2/fsig12
375 fac = fac*szz(i)/sigt2
376 fa_2 = sqrt((szx(i)/fsig12)**2 + fac*fac)
377 . + pp12*szz(i)/fsig12
378 dfmax(i,4) = max(dfmax(i,4),fa_2)
379 dfmax(i,4) = min(dfmax(i,4),one)
380 ELSE
381 fac = half/(one + pn22)/fsig12
382 fc_2 = (szx(i)*fac)**2 + (szz(i)/sigc2)**2
383 fc_2 =-fc_2*sigc2/szz(i)
384 dfmax(i,6) = max(dfmax(i,6),fc_2)
385 dfmax(i,6) = min(dfmax(i,6),one)
386C plane angle
387c FAC = SQRT(ONE + TWO*PN12(I)*SIGC2(I)/FSIG12(I)) - ONE
388c FAC = HALF*FAC*FSIG12(I)/PN12(I)
389c CC = (FAC/FSIG12(I))**2
390c CC = CC*(SIGNZX(I)/SIGNZZ(I))**2 + ONE
391c FAC = HALF/ONE + PN22(I))
392c CC = SQRT(FAC*CC)
393c ANGLE(I) = ACOS(CC)
394 ENDIF
395
396 fb_2 = sqrt(szx(i)**2 + (pn12*szz(i))**2 )
397 . + pn12*szz(i)
398 fb_2 = fb_2/fsig12
399 dfmax(i,5) = max(dfmax(i,5),fb_2)
400 dfmax(i,5) = min(dfmax(i,5),one)
401C
402 dfmax(i,1) = min(one,max(dfmax(i,1),f1,fa,fb,fc,fa_2,fb_2,fc_2))
403C
404 IF(f1 >= one .OR. fa >= one .OR. fb >= one .OR.
405 . fc >= one .OR. fa_2 >= one .OR. fb_2 >= one .OR.
406 . fc_2 >= one )THEN
407 uvar(i,1) = signxx(i)
408 uvar(i,2) = signyy(i)
409 uvar(i,3) = signzz(i)
410 uvar(i,4) = signxy(i)
411 uvar(i,5) = signyz(i)
412 uvar(i,6) = signzx(i)
413 uvar(i,7) = time
414 uvar(i,8) = four_over_5
415 nindx0= nindx0+1
416 indx0(nindx0)=i
417 ENDIF
418 ENDIF
419 ENDIF ! iflag choice
420 ENDIF
421 ENDDO
422C--------------------------------------------
423C print
424C--------------------------------------------
425 IF(nindx > 0)THEN
426 DO j=1,nindx
427 i = indx(j)
428#include "lockon.inc"
429 WRITE(iout, 1200) ngl(i),time
430 WRITE(istdo,1200) ngl(i),time
431#include "lockoff.inc"
432 END DO
433 ENDIF
434C
435 IF(nindx0 > 0)THEN
436 DO j=1,nindx0
437 i = indx0(j)
438#include "lockon.inc"
439 WRITE(iout, 1100) ngl(i),ilay,time
440 WRITE(istdo,1100) ngl(i),ilay,time
441#include "lockoff.inc"
442 END DO
443 ENDIF
444C--------------------------------------------
445 1100 FORMAT(1x,'FAILURE ELEMENT #',i10,1x,
446 .'IP #',i10,1x, 'AT TIME #:',1pe20.13)
447
448 1200 FORMAT(1x,'DELETE SOLID ELEMENT (PUCK MODEL) #',i10,1x,
449 .'AT TIME # ',1pe20.13)
450C--------------------------------------------
451 RETURN
452 END
fail_puck_s
subroutine fail_puck_s(nel, nuvar, ilay, npt0, time, timestep, uparam, ngl, off, noff, signxx, signyy, signzz, signxy, signyz, signzx, uvar, nuparam, dfmax, lf_dammx, tdele)
Definition fail_puck_s.F:37
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21