OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_fld_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_fld_c ../engine/source/materials/fail/fld/fail_fld_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!|| finterfld ../engine/source/materials/fail/fld/fail_fld_c.F
31!|| vinter2 ../engine/source/tools/curve/vinter.F
32!||====================================================================
33 SUBROUTINE fail_fld_c(
34 1 NEL ,NUPARAM ,NFUNC ,IFUNC ,
35 2 NPF ,TF ,TIME ,UPARAM ,
36 3 NGL ,IPG ,ILAY ,IPTT ,
37 4 EPSXX ,EPSYY ,EPSXY ,LF_DAMMX ,
38 5 DEPSXX ,DEPSYY ,DEPSXY ,PLA ,
39 6 ZT ,OFF ,FOFF ,TDEL ,
40 7 FLD_IDX ,DAM ,DFMAX ,DT1 ,
41 8 NIPARAM ,IPARAM ,NUVAR ,UVAR )
42C-----------------------------------------------
43c FLD failure model
44C-----------------------------------------------
45C I m p l i c i t T y p e s
46C-----------------------------------------------
47#include "implicit_f.inc"
48C-----------------------------------------------
49C G l o b a l P a r a m e t e r s
50C-----------------------------------------------
51#include "units_c.inc"
52#include "comlock.inc"
53C---------+---------+---+---+--------------------------------------------
54C VAR | SIZE |TYP| RW| DEFINITION
55C---------+---------+---+---+--------------------------------------------
56C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
57C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
58C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
59C---------+---------+---+---+--------------------------------------------
60C TIME | 1 | F | R | CURRENT TIME
61C TIMESTEP| 1 | F | R | CURRENT TIME STEP
62C---------+---------+---+---+--------------------------------------------
63C EPSXX | NEL | F | R | STRAIN XX
64C EPSYY | NEL | F | R | STRAIN YY
65C ... | | | |
66C---------+---------+---+---+--------------------------------------------
67C OFF | NEL | F | R | DELETED ELEMENT FLAG (=1. ON, =0. OFF)
68C FOFF | NEL | I |R/W| DELETED INTEGRATION POINT FLAG (=1 ON, =0 OFF)
69C DFMAX | NEL | F |R/W| MAX DAMAGE FACTOR
70C TDEL | NEL | F | W | FAILURE TIME
71C---------+---------+---+---+--------------------------------------------
72C NGL ELEMENT ID
73C IPG CURRENT GAUSS POINT (in plane)
74C ILAY CURRENT LAYER
75C IPT CURRENT INTEGRATION POINT IN THE LAYER (FOR OUTPUT ONLY)
76C---------+---------+---+---+--------------------------------------------
77C I N P U T A r g u m e n t s
78C-----------------------------------------------
79 INTEGER, INTENT(IN) :: NEL,NUPARAM,NFUNC,IPG,ILAY,IPTT,NIPARAM,NUVAR,
80 . LF_DAMMX
81 INTEGER, DIMENSION(NEL) :: NGL
82 INTEGER, DIMENSION(NFUNC) :: IFUNC
83 INTEGER, DIMENSION(NIPARAM) :: IPARAM
84 my_real, INTENT(IN) :: TIME,ZT,DT1
85 my_real, DIMENSION(NEL), INTENT(IN) :: OFF,
86 . EPSXX,EPSYY,EPSXY,DEPSXX,DEPSYY,DEPSXY,PLA
87 my_real, DIMENSION(NUPARAM) :: uparam
88C-----------------------------------------------
89C I N P U T O U T P U T A r g u m e n t s
90C-----------------------------------------------
91 INTEGER ,DIMENSION(NEL), INTENT(INOUT) :: FOFF,FLD_IDX
92 my_real ,DIMENSION(NEL,LF_DAMMX), INTENT(INOUT) :: dfmax
93 my_real ,DIMENSION(NEL), INTENT(INOUT) :: dam,tdel
94 my_real, DIMENSION(NEL,NUVAR), INTENT(INOUT) :: uvar
95C-----------------------------------------------
96C VARIABLES FOR FUNCTION INTERPOLATION
97C-----------------------------------------------
98 INTEGER NPF(*)
99 my_real finter , finterfld ,tf(*)
100 EXTERNAL finter
101C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
102C Y : y = f(x)
103C X : x
104C DYDX : f'(x) = dy/dx
105C IFUNC(J): FUNCTION INDEX
106C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
107C NPF,TF : FUNCTION PARAMETER
108C-----------------------------------------------
109C L o c a l V a r i a b l e s
110C-----------------------------------------------
111 INTEGER :: I,II,J,IENG,LENF,NINDX,NINDXF,IFAIL_SH,ISTRESS,IMARGIN
112 INTEGER ,DIMENSION(NEL) :: INDX,INDXF
113 my_real :: RANI,R1,R2,S1,S2,SS,Q,E12,FACT_MARGIN,FACT_LOOSEMETAL,
114 . FCUT,ALPHA
115 my_real ,ALLOCATABLE, DIMENSION(:) :: XF
116 my_real ,DIMENSION(NEL) :: EMAJ,EMIN,EM,DEMAJ,DEMIN,BETA
117 INTEGER, DIMENSION(NEL) :: IPOS,ILENP,IADP
118 my_real, DIMENSION(NEL) :: TAB_LOC,Y_LOC,DYDX_LOC
119C=======================================================================
120c
121 !=================================================================
122 ! - INITIALISATION OF COMPUTATION ON TIME STEP
123 !=================================================================
124 ! Recovering failure criterion parameters
125 ! -> Real parameters
126 fact_margin = uparam(1)
127 rani = uparam(3)
128 fact_loosemetal = uparam(4)
129 fcut = uparam(5)
130 IF (uparam(6) > zero) THEN
131 alpha = uparam(6)
132 ELSE
133 alpha = two*pi*fcut*dt1/(one + two*pi*fcut*dt1)
134 ENDIF
135 ! -> Integer parameters
136 ifail_sh = iparam(1)
137 imargin = iparam(2)
138 ieng = iparam(3)
139c
140 ! Shell element deletion
141 istress = 0
142 IF (ifail_sh == 1) THEN
143 istress = 0
144 ELSEIF (ifail_sh == 2) THEN
145 istress = 0
146 ELSEIF (ifail_sh == 3) THEN ! membrane criterion only
147 istress = 0
148 ELSEIF (ifail_sh == 4) THEN ! no element suppression
149 istress = -1
150 ENDIF
151c
152 ! Check if the element is already broken
153 nindx = 0
154 nindxf = 0
155 DO i = 1,nel
156 IF (off(i) == one .and. foff(i) == 1) THEN
157 nindx = nindx + 1
158 indx(nindx) = i
159 ENDIF
160 ENDDO
161c
162 !=================================================================
163 ! - minor and major(true) strain deformation
164 !=================================================================
165#include "vectorize.inc"
166 DO j = 1,nindx
167 i = indx(j)
168 ! For classic linear formulation
169 e12 = half*epsxy(i)
170 s1 = half*(epsxx(i) + epsyy(i))
171 s2 = half*(epsxx(i) - epsyy(i))
172 q = sqrt(s2**2 + e12**2)
173 emaj(i) = s1 + q
174 emin(i) = s1 - q
175 IF (emin(i) >= emaj(i)) THEN
176 ss = emin(i)
177 emin(i) = emaj(i)
178 emaj(i) = ss
179 ENDIF
180 ! For non-linear path formulation
181 e12 = half*depsxy(i)
182 s1 = half*(depsxx(i) + depsyy(i))
183 s2 = half*(depsxx(i) - depsyy(i))
184 q = sqrt(s2**2 + e12**2)
185 demaj(i) = s1 + q
186 demin(i) = s1 - q
187 ! Strain ratio filtering value for non-linear path formulation
188 demaj(i) = alpha*demaj(i) + (one - alpha)*uvar(i,2)
189 demin(i) = alpha*demin(i) + (one - alpha)*uvar(i,3)
190 beta(i) = demin(i)/sign(max(abs(demaj(i)),em20),demaj(i))
191 uvar(i,2) = demaj(i)
192 uvar(i,3) = demin(i)
193 IF (ieng == 2) THEN
194 dfmax(i,4) = beta(i)
195 ENDIF
196 ENDDO
197c
198 !=================================================================
199 ! FAILURE MAJOR STRAIN FROM INPUT CURVE AND DAMAGE RATIO
200 !=================================================================
201 ! -> Engineering strains input
202 IF (ieng == 1) THEN ! transform input fld curve to true strain
203 ii = npf(ifunc(1))
204 lenf = npf(ifunc(1)+ 1) - npf(ifunc(1))
205 ALLOCATE(xf(lenf))
206 DO i = 1,lenf
207 xf(i) = log(tf(ii + i-1) + one)
208 ENDDO
209c
210#include "vectorize.inc"
211 DO j = 1,nindx
212 i = indx(j)
213 em(i) = finterfld(emin(i),lenf,xf)
214 dam(i) = emaj(i) / em(i)
215 dfmax(i,2) = dam(i)
216 dfmax(i,1) = min(one, dam(i))
217 ENDDO
218 DEALLOCATE(xf)
219 ! -> True strains input
220 ELSE
221 ! -> Classical formulation
222 IF (ieng == 0) THEN
223#include "vectorize.inc"
224 DO j = 1,nindx
225 i = indx(j)
226 ipos(j) = 1
227 iadp(j) = npf(ifunc(1)) / 2 + 1
228 ilenp(j) = npf(ifunc(1)+1) / 2 -iadp(j) - ipos(j)
229 tab_loc(j) = emin(i)
230 ENDDO
231 CALL vinter2(tf,iadp,ipos,ilenp,nindx,tab_loc,dydx_loc,y_loc)
232#include "vectorize.inc"
233 DO j = 1,nindx
234 i = indx(j)
235 em(i) = y_loc(j)
236 dam(i) = emaj(i) / em(i)
237 dfmax(i,2) = dam(i)
238 dfmax(i,1) = min(one,dam(i))
239 ENDDO
240 ! -> Non-linear path formulation
241 ELSEIF (ieng == 2) THEN
242#include "vectorize.inc"
243 DO j = 1,nindx
244 i = indx(j)
245 ipos(j) = 1
246 iadp(j) = npf(ifunc(1)) / 2 + 1
247 ilenp(j) = npf(ifunc(1)+1) / 2 -iadp(j) - ipos(j)
248 tab_loc(j) = beta(i)
249 ENDDO
250 CALL vinter2(tf,iadp,ipos,ilenp,nindx,tab_loc,dydx_loc,y_loc)
251#include "vectorize.inc"
252 DO j = 1,nindx
253 i = indx(j)
254 em(i) = y_loc(j)
255 dam(i) = max(pla(i) / em(i),dam(i))
256 dfmax(i,2) = dam(i)
257 dfmax(i,1) = min(one,dam(i))
258 ENDDO
259 ENDIF
260 ENDIF
261c
262 !=================================================================
263 ! FLD ZONE INDEX CALCULATION FOR ANIM OUTPUT
264 !=================================================================
265 r1 = fact_loosemetal
266 r2 = rani/(rani+one)
267
268 IF (ieng < 2) THEN
269 IF (imargin == 3) THEN
270#include "vectorize.inc"
271 DO j = 1,nindx
272 i = indx(j)
273 IF (emaj(i) >= em(i)) THEN
274 fld_idx(i) = 6 ! zone 6 = failure
275 ELSEIF (emaj(i) >= em(i)*(one - fact_margin)) THEN
276 fld_idx(i) = 5 ! zone 5 = margin to fail
277 ELSEIF (emaj(i)**2 + emin(i)**2 < r1**2) THEN
278 fld_idx(i) = 1 ! zone 1 = radius 0.02
279 ELSEIF (emaj(i) >= abs(emin(i))) THEN
280 fld_idx(i) = 4 ! zone 4 = safe (45 deg line)
281 ELSEIF (emaj(i) >= r2*abs(emin(i))) THEN
282 fld_idx(i) = 3 ! zone 3 = angle atan(r/(1+r)) - compression
283 ELSE
284 fld_idx(i) = 2 ! zone 2 - high wrinkle tendency
285 ENDIF
286 dfmax(i,3) = fld_idx(i)
287 ENDDO
288 ELSE
289#include "vectorize.inc"
290 DO j = 1,nindx
291 i = indx(j)
292 IF (emaj(i) >= em(i)) THEN
293 fld_idx(i) = 6 ! zone 6 = failure
294 ELSEIF (emaj(i) >= em(i) - fact_margin) THEN
295 fld_idx(i) = 5 ! zone 5 = margin to fail
296 ELSEIF (emaj(i)**2 + emin(i)**2 < r1**2) THEN
297 fld_idx(i) = 1 ! zone 1 = radius 0.02
298 ELSEIF (emaj(i) >= abs(emin(i))) THEN
299 fld_idx(i) = 4 ! zone 4 = safe (45 deg line)
300 ELSEIF (emaj(i) >= r2*abs(emin(i))) THEN
301 fld_idx(i) = 3 ! zone 3 = angle atan(r/(1+r)) - compression
302 ELSE
303 fld_idx(i) = 2 ! zone 2 - high wrinkle tendency
304 ENDIF
305 dfmax(i,3) = fld_idx(i)
306 ENDDO
307 ENDIF
308 ELSE
309 IF (imargin == 3) THEN
310#include "vectorize.inc"
311 DO j = 1,nindx
312 i = indx(j)
313 IF (pla(i) >= em(i)) THEN
314 fld_idx(i) = max(6,fld_idx(i)) ! zone 6 = failure
315 ELSEIF (pla(i) >= em(i)*(one - fact_margin)) THEN
316 fld_idx(i) = max(5,fld_idx(i)) ! zone 5 = margin to fail
317 ELSEIF (pla(i)**2 + beta(i)**2 < r1**2) THEN
318 fld_idx(i) = max(1,fld_idx(i)) ! zone 1 = radius 0.02
319 ELSEIF (pla(i) >= abs(beta(i))) THEN
320 fld_idx(i) = max(4,fld_idx(i)) ! zone 4 = safe (45 deg line)
321 ELSEIF (pla(i) >= r2*abs(beta(i))) THEN
322 fld_idx(i) = max(3,fld_idx(i)) ! zone 3 = angle atan(r/(1+r)) - compression
323 ELSE
324 fld_idx(i) = max(2,fld_idx(i)) ! zone 2 - high wrinkle tendency
325 ENDIF
326 dfmax(i,3) = fld_idx(i)
327 ENDDO
328 ELSE
329#include "vectorize.inc"
330 DO j = 1,nindx
331 i = indx(j)
332 IF (pla(i) >= em(i)) THEN
333 fld_idx(i) = max(6,fld_idx(i)) ! zone 6 = failure
334 ELSEIF (pla(i) >= em(i) - fact_margin) THEN
335 fld_idx(i) = max(5,fld_idx(i)) ! zone 5 = margin to fail
336 ELSEIF (pla(i)**2 + beta(i)**2 < r1**2) THEN
337 fld_idx(i) = max(1,fld_idx(i)) ! zone 1 = radius 0.02
338 ELSEIF (pla(i) >= abs(beta(i))) THEN
339 fld_idx(i) = max(4,fld_idx(i)) ! zone 4 = safe (45 deg line)
340 ELSEIF (pla(i) >= r2*abs(beta(i))) THEN
341 fld_idx(i) = max(3,fld_idx(i)) ! zone 3 = angle atan(r/(1+r)) - compression
342 ELSE
343 fld_idx(i) = max(2,fld_idx(i)) ! zone 2 - high wrinkle tendency
344 ENDIF
345 dfmax(i,3) = fld_idx(i)
346 ENDDO
347 ENDIF
348 ENDIF
349c
350 !=================================================================
351 ! PRINTING OUT ELEMENT DELETION MESSAGES
352 !=================================================================
353 IF ((ifail_sh == 3 .and. zt == zero) .or. ifail_sh < 3) THEN
354 IF (ieng < 2) THEN
355#include "vectorize.inc"
356 DO j = 1,nindx
357 i = indx(j)
358 IF (emaj(i) >= em(i)) THEN
359 nindxf = nindxf + 1
360 indxf(nindxf) = i
361 foff(i) = istress
362 tdel(i) = time
363 ENDIF
364 ENDDO
365 ELSE
366#include "vectorize.inc"
367 DO j = 1,nindx
368 i = indx(j)
369 IF (pla(i) >= em(i)) THEN
370 nindxf = nindxf + 1
371 indxf(nindxf) = i
372 foff(i) = istress
373 tdel(i) = time
374 ENDIF
375 ENDDO
376 ENDIF
377 ENDIF
378 IF (nindxf > 0) THEN
379 DO j=1,nindxf
380 i = indxf(j)
381#include "lockon.inc"
382 WRITE(iout, 2000) ngl(i),ipg,ilay,iptt
383 WRITE(istdo,2100) ngl(i),ipg,ilay,iptt,time
384#include "lockoff.inc"
385 END DO
386 END IF
387c------------------------
388 2000 FORMAT(1x,'FAILURE (FLD) OF SHELL ELEMENT ',i10,1x,',GAUSS PT',i2,1x,',LAYER',i3,
389 . 1x,',INTEGRATION PT',i3)
390 2100 FORMAT(1x,'FAILURE (FLD) OF SHELL ELEMENT ',i10,1x,',GAUSS PT',i2,1x,',LAYER',i3,
391 . 1x,',INTEGRATION PT',i3,1x,'AT TIME :',1pe12.4)
392c------------------------
393 END
394c
395!||====================================================================
396!|| finterfld ../engine/source/materials/fail/fld/fail_fld_c.F
397!||--- called by ------------------------------------------------------
398!|| dam_fld_sol ../engine/source/output/h3d/h3d_results/h3d_sol_skin_scalar.F
399!|| fail_fld_c ../engine/source/materials/fail/fld/fail_fld_c.F
400!|| fail_fld_tsh ../engine/source/materials/fail/fld/fail_fld_tsh.F
401!|| fail_fld_xfem ../engine/source/materials/fail/fld/fail_fld_xfem.F
402!|| idx_fld_sol ../engine/source/output/h3d/h3d_results/h3d_sol_skin_scalar.F
403!||====================================================================
404 my_real FUNCTION finterfld(EPST,LENF,XF)
405C-----------------------------------------------
406C I m p l i c i t T y p e s
407C-----------------------------------------------
408#include "implicit_f.inc"
409C-----------------------------------------------
410C D u m m y A r g u m e n t s
411C-----------------------------------------------
412 INTEGER lenf
413 my_real xf(lenf),epst
414C-----------------------------------------------
415C L o c a l V a r i a b l e s
416C-----------------------------------------------
417 INTEGER i
418 my_real dx1,dx2,DERI
419C=======================================================================
420 finterfld = zero
421 dx2 = xf(1) - epst
422 DO i = 2,lenf-2,2
423 dx1 = -dx2
424 dx2 = xf(i+1) - epst
425 IF (dx2 >= zero .OR. i == lenf-2) THEN
426 deri = (xf(i+2) - xf(i)) / (xf(i+1) - xf(i-1))
427 IF (dx1 <= dx2) THEN
428 finterfld = xf(i) + dx1 * deri
429 ELSE
430 finterfld = xf(i+2) - dx2 * deri
431 ENDIF
432 EXIT
433 ENDIF
434 ENDDO
435c-----------
436 RETURN
437 END
my_real
#define my_real
Definition cppsort.cpp:32
fail_fld_c
subroutine fail_fld_c(nel, nuparam, nfunc, ifunc, npf, tf, time, uparam, ngl, ipg, ilay, iptt, epsxx, epsyy, epsxy, lf_dammx, depsxx, depsyy, depsxy, pla, zt, off, foff, tdel, fld_idx, dam, dfmax, dt1, niparam, iparam, nuvar, uvar)
Definition fail_fld_c.F:42
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
vinter2
subroutine vinter2(tf, iad, ipos, ilen, nel0, x, dydx, y)
Definition vinter.F:143