OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_tensstrain_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_tensstrain_c ../engine/source/materials/fail/tensstrain/fail_tensstrain_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_tensstrain_c(
32 1 NEL ,NFUNC ,NUPARAM ,NUVAR ,IFUNC ,
33 2 UPARAM ,UVAR ,NPF ,TF ,TIME ,
34 3 NGL ,IPG ,ILAY ,IPT ,EPSP ,
35 4 EPSXX ,EPSYY ,EPSXY ,EPSYZ ,EPSZX ,
36 5 OFF ,FOFF ,DFMAX ,TDEL ,
37 6 DMG_FLAG ,DMG_SCALE,ALDT ,TSTAR ,ISMSTR )
38C-----------------------------------------------
39C tension strain failure model
40C-----------------------------------------------
41C I m p l i c i t T y p e s
42C-----------------------------------------------
43#include "implicit_f.inc"
44C-----------------------------------------------
45C G l o b a l P a r a m e t e r s
46C-----------------------------------------------
47#include "units_c.inc"
48#include "comlock.inc"
49C---------+---------+---+---+--------------------------------------------
50C VAR | SIZE |TYP| RW| DEFINITION
51C---------+---------+---+---+--------------------------------------------
52C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
53C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
54C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
55C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
56C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
57C---------+---------+---+---+--------------------------------------------
58C NFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW not used
59C IFUNC | NFUNC | I | R | FUNCTION INDEX not used
60C NPF | * | I | R | FUNCTION ARRAY
61C TF | * | F | R | FUNCTION ARRAY
62C---------+---------+---+---+--------------------------------------------
63C EPSP | NEL | F | R | EQUIVALENT STRAIN RATE
64C EPSXX | NEL | F | R | STRAIN XX
65C EPSYY | NEL | F | R | STRAIN YY
66C ... | | | |
67C---------+---------+---+---+--------------------------------------------
68C OFF | NEL | F | R | DELETED ELEMENT FLAG (=1. ON, =0. OFF)
69C FOFF | NEL | I |R/W| DELETED INTEGRATION POINT FLAG (=1 ON, =0 OFF)
70C DFMAX | NEL | F |R/W| MAX DAMAGE FACTOR
71C TDEL | NEL | F | W | FAILURE TIME
72C DMG_FLAG| 1 | I | W | STRESS REDUCTION FLAG DUE TO DAMAGE
73C DMG_SCALE| NEL | F | W | STRESS REDUCTION FACTOR
74C---------+---------+---+---+--------------------------------------------
75C NGL ELEMENT ID
76C IPG CURRENT GAUSS POINT (in plane)
77C ILAY CURRENT LAYER
78C IPT CURRENT INTEGRATION POINT IN THE LAYER
79C---------+---------+---+---+--------------------------------------------
80C I N P U T A r g u m e n t s
81C-----------------------------------------------
82 INTEGER ,INTENT(IN) :: NEL,NUPARAM,NUVAR,IPG,ILAY,IPT,ISMSTR
83 INTEGER ,DIMENSION(NEL) ,INTENT(IN) :: NGL
84 my_real ,INTENT(IN) :: TIME
85 my_real ,DIMENSION(NEL) ,INTENT(IN) :: OFF,EPSP,ALDT,
86 . EPSXX,EPSYY,EPSXY,EPSYZ,EPSZX,TSTAR
87 my_real,DIMENSION(NUPARAM) ,INTENT(IN) :: 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 ,INTENT(OUT) :: DMG_FLAG
92 INTEGER ,DIMENSION(NEL) ,INTENT(INOUT) :: FOFF
93 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: dfmax
94 my_real ,DIMENSION(NEL) ,INTENT(OUT) :: tdel,dmg_scale
95 my_real ,DIMENSION(NEL,NUVAR) ,INTENT(INOUT) :: uvar
96C-----------------------------------------------
97C VARIABLES FOR FUNCTION INTERPOLATION
98C-----------------------------------------------
99 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
100 my_real finter ,tf(*)
101 EXTERNAL finter
102C-----------------------------------------------
103C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
104C Y : y = f(x)
105C X : x
106C DYDX : f'(x) = dy/dx
107C IFUNC(J): FUNCTION INDEX
108C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
109C NPF,TF : FUNCTION PARAMETER
110C-----------------------------------------------
111C L o c a l V a r i a b l e s
112C-----------------------------------------------
113 INTEGER I,J,NINDXD,NINDXF,NINDX1,NINDX2,S_FLAG,SCALE_FLAG,STRDEF,STRFLAG
114 INTEGER ,DIMENSION(NEL) :: INDXD,INDXF,INDX1,INDX2
115 my_real R1,R2,DYDX,EPSP1,EPSP2,EPSF1,EPSF2,E22,LAMBDA,FAC,DF,
116 . EL_REF,SC_EL,SCALE_TEMP,RFAC,RFAC2,EPST_A,EPST_B,
117 . UNIT_T,EPSP_UNIT
118 my_real ,DIMENSION(NEL) :: EPS_MAX,EPST1,EPST2,DAMAGE,EPSRATE
119C=======================================================================
120 dmg_flag = 1
121 s_flag = 1
122 scale_flag = 1
123 epsf1 = uparam(1)
124 epsf2 = uparam(2)
125 epsp1 = uparam(3)
126 epsp2 = uparam(4)
127 scale_temp = one
128
129 damage(1:nel)=zero
130C-----------------------------------------------
131C! Initialization of initial el.length
132C-----------------------------------------------
133 sc_el = uparam(5)
134 el_ref = uparam(6)
135 scale_temp = uparam(7)
136 s_flag = int(uparam(8))
137 unit_t = uparam(9)
138 strdef= int(uparam(10))
139 strflag = 0
140c
141 IF (uvar(1,1) == zero) THEN
142 DO i=1,nel
143 uvar(i,1) = aldt(i)
144 SELECT CASE (s_flag)
145 CASE (2)
146 IF (ifunc(2) /= 0)THEN
147 lambda = uvar(i,1) / el_ref
148 fac = sc_el * finter(ifunc(2),lambda,npf,tf,df)
149 uvar(i,1) = fac
150 ELSE
151 uvar(i,1) = one
152 ENDIF
153 CASE (3)
154 IF (ifunc(2) /= 0)THEN
155 lambda = uvar(i,1) / el_ref
156 fac = sc_el * finter(ifunc(2),lambda,npf,tf,df)
157 uvar(i,1) = fac
158 ELSE
159 uvar(i,1) = one
160 ENDIF
161 CASE DEFAULT ! old formulation, backward compatibility only
162 uvar(i,1) = epsp1
163 END SELECT
164 ENDDO
165 ENDIF
166c
167 IF (strdef == 2) THEN ! failure defined as engineering strain
168 IF (ismstr == 10 .or. ismstr == 12) THEN
169 strflag = 1 ! Cauchy-Green to engineering
170 ELSE IF (ismstr == 0 .or. ismstr == 2 .or. ismstr == 4) THEN
171 strflag = 2 ! true strain to engineering
172 END IF
173 ELSE IF (strdef == 3) THEN ! failure defined as true strain
174 IF (ismstr == 1 .or. ismstr == 3 .or. ismstr == 11) THEN
175 strflag = 3 ! engineering to true strain
176 ELSE IF (ismstr == 10 .or. ismstr == 12) THEN
177 strflag = 4 ! Cauchy-Green to true strain
178 END IF
179 END IF
180c------------------------------------
181c Equivalent or principal strain
182c------------------------------------
183 SELECT CASE (s_flag)
184c
185 CASE (1,2) ! equivalent strain criterion
186 IF (strflag == 1 .or. strflag == 4) THEN
187 DO i=1,nel !COMPUTE THE 2 PRINCIPAL VALUES
188 epst_a = half*(epsxx(i)+epsyy(i))
189 epst_b = sqrt( (half*(epsxx(i)-epsyy(i)))**2 + (half*epsxy(i))**2)
190 epst1(i) = epst_a + epst_b
191 epst2(i) = epst_a - epst_b
192 ENDDO
193 ELSE
194 DO i=1,nel
195 eps_max(i) = half*(epsxx(i)+epsyy(i)
196 . + sqrt( (epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
197 . + epsxy(i)*epsxy(i) ) )
198 ENDDO
199 END IF
200c
201 CASE (3) ! max tensile principal strain criterion
202 DO i=1,nel
203 epst_a = half*(epsxx(i)+epsyy(i))
204 epst_b = sqrt((half*(epsxx(i)-epsyy(i)))**2 + (half*epsxy(i))**2)
205 epst1(i) = epst_a + epst_b
206 epst2(i) = epst_a - epst_b
207 eps_max(i) = epst1(i)
208 ENDDO
209 END SELECT
210c----------------------------------------------
211c Max strain transformation following input definition
212c--------------------------
213 SELECT CASE (strflag)
214 CASE (1) ! transform Cauchy-Green to engineering
215 DO i=1,nel
216 IF (off(i) == one ) THEN
217 epst1(i) = sqrt(epst1(i) + one) - one
218 epst2(i) = sqrt(epst2(i) + one) - one
219 eps_max(i) = sqrt(epst1(i)**2 + epst2(i)**2)
220 epsrate(i) = sqrt(epsp(i) + one) - one
221 END IF
222 ENDDO
223c
224 CASE (2) ! transform true strain to engineering
225 DO i=1,nel
226 IF (off(i) == one ) THEN
227 eps_max(i) = exp(eps_max(i)) - one
228 epsrate(i) = exp(epsp(i)) - one
229 END IF
230 ENDDO
231c
232 CASE (3) ! transform engineering to true strain
233 DO i=1,nel
234 IF (off(i) == one ) THEN
235 eps_max(i) = log(eps_max(i) + one)
236 epsrate(i) = log(epsp(i) + one)
237 END IF
238 ENDDO
239c
240 CASE (4) ! transform Cauchy-Green to true strain
241 DO i=1,nel
242 IF (off(i) == one ) THEN
243 epst1(i) = log(sqrt(epst1(i)+one))
244 epst2(i) = log(sqrt(epst2(i)+one))
245 eps_max(i) = sqrt(epst1(i)**2 + epst2(i)**2)
246 epsrate(i) = log(sqrt(epsp(i)+one))
247 END IF
248 ENDDO
249c
250 CASE DEFAULT
251 ! no transformation : failure strain measure is defined by Ismstr
252 epsrate(1:nel) = epsp(1:nel)
253 END SELECT
254c----------------------------------------------
255c Calculate failure criterion
256c-------------------
257c
258 nindxd = 0
259 nindxf = 0
260 nindx1 = 0
261 nindx2 = 0
262c
263 SELECT CASE (s_flag)
264 CASE (2)
265C-------------------
266c Equivalent strain. New format: Element length and temperature
267C-------------------
268 DO i=1,nel
269 IF (off(i) == one .and. foff(i) == 1) THEN
270 IF (ifunc(1) > 0) THEN ! strain rate
271 epsp_unit = epsrate(i) * unit_t
272 rfac = finter(ifunc(1),epsp_unit,npf,tf,dydx)
273 rfac = max(rfac,em20)
274 ELSE
275 rfac = one
276 ENDIF
277 IF (ifunc(3) > 0) THEN ! temperatuer
278 rfac2 = finter(ifunc(3),tstar(i),npf,tf,dydx)
279 rfac2 = max(rfac2,em20)
280 ELSE
281 rfac2 = one
282 ENDIF
283 r1 = epsf1*rfac*rfac2*uvar(i,1)
284 r2 = epsf2*rfac*rfac2*uvar(i,1)
285 IF (eps_max(i) > r1) THEN
286 IF (dfmax(i) == zero) THEN
287 nindxd = nindxd + 1
288 indxd(nindxd) = i
289 ENDIF
290 damage(i) = min(one, (eps_max(i)-r1)/(r2-r1))
291 ENDIF
292 IF (eps_max(i) >= r2) THEN
293 foff(i) = 0
294 tdel(i) = time
295 nindxf = nindxf + 1
296 indxf(nindxf) = i
297 ENDIF
298 ENDIF
299 ENDDO
300
301 CASE (3)
302C-------------------
303c Principal strain. New format: No failure in compression
304C-------------------
305 DO i=1,nel
306 IF (off(i) == one .and. foff(i) == 1) THEN
307 IF (ifunc(1) > 0) THEN ! strain rate
308 epsp_unit = epsrate(i) * unit_t
309 rfac = finter(ifunc(1),epsp_unit,npf,tf,dydx)
310 rfac = max(rfac,em20)
311 ELSE
312 rfac = one
313 ENDIF
314 IF (ifunc(3) > 0) THEN ! temperatuer
315 rfac2 = finter(ifunc(3),tstar(i),npf,tf,dydx)
316 rfac2 = max(rfac2,em20)
317 ELSE
318 rfac2 = one
319 ENDIF
320 r1 = epsf1*rfac*rfac2*uvar(i,1)
321 r2 = epsf2*rfac*rfac2*uvar(i,1)
322 IF ( eps_max(i) > r1) THEN
323 IF (dfmax(i) == zero) THEN
324 nindxd = nindxd + 1
325 indxd(nindxd) = i
326 ENDIF
327 damage(i) = min(one, (eps_max(i)-r1)/(r2-r1))
328 ENDIF
329 IF ( eps_max(i) >= r2) THEN
330 foff(i) = 0
331 tdel(i) = time
332 nindxf = nindxf + 1
333 indxf(nindxf) = i
334 ENDIF
335 ENDIF
336 ENDDO
337
338C-------------------
339c Equivalent strain. No change for downward compatibility.
340C-------------------
341 CASE DEFAULT
342
343 DO i=1,nel
344 IF (off(i) == one .and. foff(i) == 1) THEN
345 IF (ifunc(1) > 0) THEN
346 epsp_unit = epsrate(i) * unit_t
347 rfac = finter(ifunc(1),epsp_unit,npf,tf,dydx)
348 rfac = max(rfac,em20)
349 ELSE
350 rfac = one
351 ENDIF
352 r1 = epsf1*rfac
353 r2 = epsf2*rfac
354 IF (eps_max(i) > r1) THEN
355 IF (dfmax(i) == zero) THEN
356 nindxd = nindxd + 1
357 indxd(nindxd) = i
358 ENDIF
359 damage(i) = min(one, (eps_max(i)-r1)/(r2-r1))
360 ENDIF
361 IF (eps_max(i) >= r2) THEN
362 foff(i) = 0
363 tdel(i) = time
364 nindxf = nindxf + 1
365 indxf(nindxf) = i
366 ENDIF
367c
368 e22 = zero
369 IF (epsp1 > zero .OR. epsp2 > zero) THEN
370 IF (epsp2 > zero) THEN
371 e22 = half*( epsxx(i)+epsyy(i)
372 . - sqrt( (epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
373 . + epsxy(i)*epsxy(i) ) ) !PRINCIPAL STRAIN MIN
374 IF (strflag == 1) THEN
375 e22 = sqrt(e22 + one) - one
376 ELSE IF (strflag == 2) THEN
377 e22 = exp(e22) - one
378 ELSE IF (strflag == 3) THEN
379 e22 = log(e22 + one)
380 ELSE IF (strflag == 4) THEN
381 e22 = log(sqrt(e22+one))
382 END IF
383 ENDIF
384 IF (eps_max(i) > uvar(i,1) .AND. uvar(i,1) /= zero) THEN
385 foff(i) = 0
386 tdel(i) = time
387 nindx1 = nindx1 + 1
388 indx1(nindx1) = i
389 ELSEIF (e22 > epsp2 .AND. epsp2 /= zero) THEN
390 foff(i) = 0
391 tdel(i) = time
392 nindx2 = nindx2 + 1
393 indx2(nindx2) = i
394 eps_max(i) = e22
395 ENDIF
396 dfmax(i) = max(dfmax(i),min( one, (eps_max(i) / uvar(i,1))))
397 scale_flag = 0
398 ENDIF
399c
400 ENDIF
401 ENDDO
402 END SELECT
403c-----------------------------------------------------------------------
404c--- Maximum Damage for output : 0 < DFMAX < 1
405 DO i=1,nel
406 dfmax(i) = max(dfmax(i), damage(i))
407 IF (scale_flag == 1) dmg_scale(i) = one - dfmax(i)
408 ENDDO
409c------------------------
410c print
411c------------------------
412 IF (nindxd > 0) THEN ! start damage
413 DO j=1,nindxd
414 i = indxd(j)
415 IF (s_flag == 3) THEN
416#include "lockon.inc"
417 WRITE(iout, 2001) ngl(i),ipg,ilay,ipt,eps_max(i)
418 WRITE(istdo,2001) ngl(i),ipg,ilay,ipt,eps_max(i)
419#include "lockoff.inc"
420 ELSE
421#include "lockon.inc"
422 WRITE(iout, 2002) ngl(i),ipg,ilay,ipt,eps_max(i)
423 WRITE(istdo,2002) ngl(i),ipg,ilay,ipt,eps_max(i)
424#include "lockoff.inc"
425 ENDIF
426 END DO
427 END IF
428c
429 IF (nindxf > 0) THEN ! failure
430 DO j=1,nindxf
431 i = indxf(j)
432 IF (s_flag == 3) THEN
433#include "lockon.inc"
434 WRITE(iout, 3001) ngl(i),ipg,ilay,ipt,eps_max(i)
435 WRITE(istdo,3101) ngl(i),ipg,ilay,ipt,time,eps_max(i)
436#include "lockoff.inc"
437 ELSE
438 WRITE(iout, 3002) ngl(i),ipg,ilay,ipt,eps_max(i)
439 WRITE(istdo,3102) ngl(i),ipg,ilay,ipt,time,eps_max(i)
440 ENDIF
441 END DO
442 END IF
443c
444 IF (nindx1 > 0) THEN
445 DO j=1,nindx1
446 i = indx1(j)
447#include "lockon.inc"
448 WRITE(iout, 4000) ngl(i),ipg,ilay,ipt,eps_max(i)
449 WRITE(istdo,4100) ngl(i),ipg,ilay,ipt,time,eps_max(i)
450#include "lockoff.inc"
451 END DO
452 END IF
453c
454 IF (nindx2 > 0) THEN
455 DO j=1,nindx2
456 i = indx2(j)
457#include "lockon.inc"
458 WRITE(iout, 5000) ngl(i),ipg,ilay,ipt,eps_max(i)
459 WRITE(istdo,5100) ngl(i),ipg,ilay,ipt,time,eps_max(i)
460#include "lockoff.inc"
461 END DO
462 END IF
463c-----------------------------------------------
464 2001 FORMAT(1x,'START DAMAGE (TENS) OF SHELL ELEMENT ',i10,1x,',GAUSS PT',
465 . i2,1x,',LAYER',i3,1x,',INTEGRATION PT',i3,/
466 . 1x,'PRINCIPAL STRAIN=',g11.4)
467 2002 FORMAT(1x,'START DAMAGE (TENS) OF SHELL ELEMENT ',i10,1x,',GAUSS PT',
468 . i2,1x,',LAYER',i3,1x,',INTEGRATION PT',i3,/
469 . 1x,'EQUIVALENT STRAIN=',g11.4)
470c---
471 3001 FORMAT(1x,'FAILURE (TENS) OF SHELL ELEMENT ',i10,1x,',GAUSS PT',
472 . i2,1x,',LAYER',i3,1x,',INTEGRATION PT',i3,/
473 . 1x,'PRINCIPAL STRAIN=',g11.4)
474 3101 FORMAT(1x,'failure(tens) of shell element ',I10,1X,',gauss pt',
475 . I2,1X,',layer',I3,1X,',integration pt',I3,
476 . 1X,'at time :',1PE12.4,/1X,'principal strain=',G11.4)
477 3002 FORMAT(1X,'failure(tens) of shell element ',i10,1x,',GAUSS PT',
478 . i2,1x,',LAYER',i3,1x,',INTEGRATION PT',i3,/
479 . 1x,'EQUIVALENT STRAIN=',g11.4)
480 3102 FORMAT(1x,'FAILURE (TENS) OF SHELL ELEMENT ',i10,1x,',GAUSS PT',
481 . i2,1x,',LAYER',i3,1x,',INTEGRATION PT',i3,
482 . 1x,'AT TIME :',1pe12.4,/1x,'EQUIVALENT STRAIN=',g11.4)
483 4000 FORMAT(1x,'FAILURE (TENS) OF SHELL ELEMENT ',i10,1x,',GAUSS PT',
484 . i2,1x,',LAYER',i3,1x,',INTEGRATION PT',i3,/
485 . 1x,'1st PRINCIPAL STRAIN=',g11.4)
486 4100 FORMAT(1x,'FAILURE (TENS) OF SHELL ELEMENT ',i10,1x,',GAUSS PT',
487 . i2,1x,',LAYER',i3,1x,',INTEGRATION PT',i3,
488 . 1x,'AT TIME :',1pe12.4,/1x,'1st PRINCIPAL STRAIN=',g11.4)
489 5000 FORMAT(1x,'FAILURE (TENS) OF SHELL ELEMENT ',i10,1x,',GAUSS PT',
490 . i2,1x,',LAYER',i3,1x,',INTEGRATION PT',i3,/
491 . 1x,'2nd PRINCIPAL STRAIN=',g11.4)
492 5100 FORMAT(1x,'FAILURE (TENS) OF SHELL ELEMENT ',i10,1x,',GAUSS PT',
493 . i2,1x,',LAYER',i3,1x,',INTEGRATION PT',i3,
494 . 1x,'AT TIME :',1pe12.4,/1x, '2ndPRINCIPAL STRAIN=',g11.4)
495c-----------------------------------------------
496 RETURN
497 END
fail_tensstrain_c
subroutine fail_tensstrain_c(nel, nfunc, nuparam, nuvar, ifunc, uparam, uvar, npf, tf, time, ngl, ipg, ilay, ipt, epsp, epsxx, epsyy, epsxy, epsyz, epszx, off, foff, dfmax, tdel, dmg_flag, dmg_scale, aldt, tstar, ismstr)
Definition fail_tensstrain_c.F:38
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21