OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_tensstrain_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/.
23!||====================================================================
24!|| fail_tensstrain_s ../engine/source/materials/fail/tensstrain/fail_tensstrain_s.f
25!||--- called by ------------------------------------------------------
26!|| mmain ../engine/source/materials/mat_share/mmain.F90
27!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
28!|| mulaw8 ../engine/source/materials/mat_share/mulaw8.F90
29!|| usermat_solid ../engine/source/materials/mat_share/usermat_solid.F
30!||--- calls -----------------------------------------------------
31!|| finter ../engine/source/tools/curve/finter.F
32!||====================================================================
33 SUBROUTINE fail_tensstrain_s(
34 1 NEL ,NUPARAM,NUVAR ,NFUNC ,IFUNC ,
35 2 NPF ,TF ,TIME ,TIMESTEP,UPARAM ,
36 3 NGL ,ALDT ,TSTAR ,ISMSTR ,
37 4 EPSXX ,EPSYY ,EPSZZ ,EPSXY ,EPSYZ ,EPSZX ,
38 5 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
39 6 EPSP ,UVAR ,OFF ,DFMAX ,TDELE ,
40 7 MFXX ,MFXY ,MFXZ ,MFYX ,MFYY ,MFYZ ,
41 8 MFZX ,MFZY ,MFZZ ,DMG_SCALE)
42C-----------------------------------------------
43C Tensile Strain failure criterion
44C-----------------------------------------------
45C I m p l i c i t T y p e s
46C-----------------------------------------------
47#include "implicit_f.inc"
48C---------+---------+---+---+--------------------------------------------
49C VAR | SIZE |TYP| RW| DEFINITION
50C---------+---------+---+---+--------------------------------------------
51C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
52C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
53C NUVAR | 1 | I | R | NUMBER OF FAILURE ELEMENT VARIABLES
54C---------+---------+---+---+--------------------------------------------
55C NFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW
56C IFUNC | NFUNC | I | R | FUNCTION INDEX
57C NPF | * | I | R | FUNCTION ARRAY
58C TF | * | F | R | FUNCTION ARRAY
59C---------+---------+---+---+--------------------------------------------
60C TIME | 1 | F | R | CURRENT TIME
61C TIMESTEP| 1 | F | R | CURRENT TIME STEP
62C UPARAM | NUPARAM | F | R | USER FAILURE PARAMETER ARRAY
63C---------+---------+---+---+--------------------------------------------
64C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
65C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
66C ... | | | |
67C ... | | | |
68C---------+---------+---+---+--------------------------------------------
69C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
70C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
71C---------+---------+---+---+--------------------------------------------
72#include "scr17_c.inc"
73#include "units_c.inc"
74#include "comlock.inc"
75#include "param_c.inc"
76#include "impl1_c.inc"
77C-----------------------------------------------
78C I N P U T A r g u m e n t s
79C-----------------------------------------------
80 INTEGER :: NEL,NUPARAM,NUVAR,ISMSTR
81 INTEGER ,DIMENSION(NEL) :: NGL
82 my_real :: TIME,TIMESTEP
83 my_real ,DIMENSION(NUPARAM) :: UPARAM
84 my_real ,DIMENSION(NEL) :: SIGNXX,SIGNYY,SIGNZZ,SIGNXY,SIGNYZ,SIGNZX,
85 . EPSXX,EPSYY,EPSZZ,EPSXY,EPSYZ,EPSZX,EPSP,ALDT,TSTAR,
86 . MFXX,MFXY,MFXZ,MFYX,MFYY,MFYZ,MFZX,MFZY,MFZZ
87 my_real, DIMENSION(NEL), INTENT(INOUT) :: DMG_SCALE
88C-----------------------------------------------
89C I N P U T O U T P U T A r g u m e n t s
90C-----------------------------------------------
91 my_real uvar(nel,nuvar), off(nel), dfmax(nel),tdele(nel)
92C-----------------------------------------------
93C VARIABLES FOR FUNCTION INTERPOLATION
94C-----------------------------------------------
95 INTEGER NPF(*), IFUNC(NFUNC),NFUNC
96 my_real FINTER ,TF(*),DF
97 EXTERNAL FINTER
98C Y = FINTER(IFUNC(J),X,NPF,TF,DF)
99C Y : y = f(x)
100C X : x
101C DF : f'(x) = dy/dx
102C IFUNC(J): FUNCTION INDEX
103C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
104C NPF,TF : FUNCTION PARAMETER
105C-----------------------------------------------
106C L o c a l V a r i a b l e s
107C-----------------------------------------------
108 INTEGER I,J,NINDXD,NINDXF,S_FLAG,STRDEF,STRFLAG
109 INTEGER ,DIMENSION(NEL) :: IFLAG,INDXD,INDXF
110 my_real :: E1,E2,E3,E4,E5,E6,RFAC,RFAC2,E42,E52,E62,C,D,EPST,EPST2,
111 . R1,R2,Y,YP,DAV,DYDX,IE_SP,P,SCALE,CC,Y1,I1,I2,I3,Q,R,PHI,
112 . R_INTER,EL_REF,SC_EL,EPSP1,EPSP2,FAC,SCALE_TEMP,
113 . E11,E22,E33,EPSF1,EPSF2,LAMBDA,UNIT_T,EPSP_UNIT
114 my_real,DIMENSION(NEL) :: EPS11,EPS22,EPS33,EPS12,EPS23,EPS31,
115 . EPS_MAX,DAMAGE,RIEF1,RIEF2
116C=======================================================================
117 epsf1 = uparam(1)
118 epsf2 = uparam(2)
119 epsp1 = uparam(3)
120 epsp2 = uparam(4)
121 sc_el = uparam(5)
122 el_ref = uparam(6)
123 scale_temp = uparam(7)
124 s_flag = int(uparam(8))
125 unit_t = uparam(9)
126 strdef = int(uparam(10))
127c
128 damage(:nel) = zero
129 eps_max(:nel) = zero
130 nindxd = 0
131 nindxf = 0
132 strflag = 0
133C-----------------------------------------------
134c Initialization
135C-----------------------------------------------
136 IF (uvar(1,1) == zero) THEN
137 DO i=1,nel
138 SELECT CASE (s_flag)
139 CASE (1) ! default - old formulation
140 uvar(i,1) = epsp1
141 CASE (2)
142 IF (ifunc(2) > 0)THEN
143 lambda = aldt(i) / el_ref
144 fac = sc_el * finter(ifunc(2),lambda,npf,tf,df)
145 uvar(i,1) = fac
146 ELSE
147 uvar(i,1) = one
148 ENDIF
149 CASE (3)
150 IF (ifunc(2) > 0)THEN
151 lambda = aldt(i) / el_ref
152 fac = sc_el * finter(ifunc(2),lambda,npf,tf,df)
153 uvar(i,1) = fac
154 ELSE
155 uvar(i,1) = one
156 ENDIF
157 END SELECT
158 ENDDO
159 ENDIF
160c-----------------------------------------------
161 DO i=1,nel
162 IF (off(i) < em01) off(i)=zero
163 IF (off(i) < one) off(i)=off(i)*four_over_5
164 END DO
165c----------------------------------------------
166c Max strain transformation following input definition
167c-------------------
168 IF (strdef == 2) THEN ! failure defined as engineering strain
169 IF (ismstr == 10 .or. ismstr == 12) THEN
170 strflag = 1
171 ELSE IF (ismstr == 0 .or. ismstr == 2 .or. ismstr == 4) THEN
172c transform true strain to engineering
173 strflag = 2
174 END IF
175 ELSE IF (strdef == 3) THEN ! failure defined as true strain
176 IF (ismstr == 1 .or. ismstr == 3 .or. ismstr == 11) THEN
177c transform engineering to true strain
178 strflag = 3
179 ELSE IF (ismstr == 10 .or. ismstr == 12) THEN
180 strflag = 4
181 END IF
182 END IF
183c
184c-------------------------
185 SELECT CASE (s_flag)
186c-------------------------
187 CASE (1) ! Equivalent strain criterion (old). Backward compatibility only
188c-------------------
189 IF (strflag == 1) THEN ! Radioss strain is Cauchy-Green
190c transform total grad def strain to engineering
191 DO i=1,nel
192 eps11(i) = mfxx(i)
193 eps22(i) = mfyy(i)
194 eps33(i) = mfzz(i)
195 eps12(i) = mfxy(i) + mfyx(i)
196 eps23(i) = mfzy(i) + mfyz(i)
197 eps31(i) = mfxz(i) + mfzx(i)
198 END DO
199 ELSE IF (strflag == 2) THEN
200c transform engineering strain to true
201 DO i=1,nel
202 eps11(i) = log(epsxx(i) + one)
203 eps22(i) = log(epsyy(i) + one)
204 eps33(i) = log(epszz(i) + one)
205 eps12(i) = log(epsxy(i) + one)
206 eps23(i) = log(epsyz(i) + one)
207 eps31(i) = log(epszx(i) + one)
208 END DO
209 ELSE IF (strflag == 3) THEN
210c transform engineering strain to true
211 DO i=1,nel
212 eps11(i) = log(epsxx(i) + one)
213 eps22(i) = log(epsyy(i) + one)
214 eps33(i) = log(epszz(i) + one)
215 eps12(i) = log(epsxy(i) + one)
216 eps23(i) = log(epsyz(i) + one)
217 eps31(i) = log(epszx(i) + one)
218 END DO
219 ELSE IF (strflag == 4) THEN ! Radioss strain is Cauchy-Green
220c transform total grad def strain to true
221 DO i=1,nel
222 eps11(i) = log(mfxx(i) + one)
223 eps22(i) = log(mfyy(i) + one)
224 eps33(i) = log(mfzz(i) + one)
225 eps12(i) = log(mfxy(i) + mfyx(i) + one)
226 eps23(i) = log(mfxz(i) + mfzx(i) + one)
227 eps31(i) = log(mfzy(i) + mfyz(i) + one)
228 END DO
229 ELSE
230 eps11(1:nel) = epsxx(1:nel)
231 eps22(1:nel) = epsyy(1:nel)
232 eps33(1:nel) = epszz(1:nel)
233 eps12(1:nel) = epsxy(1:nel)
234 eps23(1:nel) = epsyz(1:nel)
235 eps31(1:nel) = epszx(1:nel)
236 END IF
237c calculate max strain
238 DO i=1,nel
239 IF (off(i) == one ) THEN
240 dav = (eps11(i)+eps22(i)+eps33(i))*third
241 e1 = eps11(i) - dav
242 e2 = eps22(i) - dav
243 e3 = eps33(i) - dav
244 e4 = half*eps12(i)
245 e5 = half*eps23(i)
246 e6 = half*eps31(i)
247 e42 = e4*e4
248 e52 = e5*e5
249 e62 = e6*e6
250 c = (e1*e2 + e1*e3 + e2*e3) - e42 - e52 - e62
251 d = - e1*e2*e3 + e1*e52 + e2*e62 + e3*e42 - two*e4*e5*e6
252 cc = c*third
253 epst = sqrt(-cc)
254 epst2 = epst * epst
255 y = (epst2 + c)* epst + d
256 y1 = -(epst2 + c)* epst + d
257 IF(abs(y) > em8 .OR.(abs(y) < em8 .AND. abs(y1) < abs(y)))THEN
258 epst = 1.75 * epst
259 epst2 = epst * epst
260 y = (epst2 + c)* epst + d
261 yp = three*epst2 + c
262 IF(yp /= zero)epst = epst - y/yp
263 epst2 = epst * epst
264 y = (epst2 + c)* epst + d
265 yp = three*epst2 + c
266 IF(yp /= zero)epst = epst - y/yp
267 epst2 = epst * epst
268 y = (epst2 + c)* epst + d
269 yp = three*epst2 + c
270 IF(yp /= zero)epst = epst - y/yp
271
272 epst2 = epst * epst
273 y = (epst2 + c)* epst + d
274 yp = three*epst2 + c
275 IF(yp /= zero)epst = epst - y/yp
276C
277 ENDIF
278 epst = epst + dav
279 eps_max(i) = epst
280 END IF
281 ENDDO
282c
283c test failure crit
284c
285 DO i=1,nel
286 IF (off(i) == one)THEN
287 rfac = one
288 IF (ifunc(1) > 0)THEN
289 epsp_unit = epsp(i) * unit_t
290 rfac = finter(ifunc(1),epsp_unit,npf,tf,dydx)
291 rfac = max(rfac,em20)
292 ENDIF
293 r1 = epsf1*rfac
294 r2 = epsf2*rfac
295c
296 IF(eps_max(i) > r1.AND.r1 < r2) THEN
297 IF (dfmax(i) == zero) THEN
298 nindxd = nindxd + 1
299 indxd(nindxd) = i
300#include "lockon.inc"
301 WRITE(iout, 2002) ngl(i),eps_max(i)
302 WRITE(istdo,2002) ngl(i),eps_max(i)
303#include "lockoff.inc"
304 ENDIF
305 damage(i)= (eps_max(i)-r1)/(r2-r1)
306 damage(i)= min(one,damage(i))
307 ENDIF
308 IF (eps_max(i) >= r2) THEN
309 damage(i)= one
310 off(i)=four_over_5
311 nindxf=nindxf+1
312 indxf(nindxf)=i
313 tdele(i) = time
314#include "lockon.inc"
315 WRITE(iout, 3000) ngl(i),eps_max(i)
316 WRITE(istdo,3100) ngl(i),eps_max(i),time
317#include "lockoff.inc"
318 ENDIF
319
320 IF (epsp1 > zero .OR. epsp2 > zero)THEN
321c! from: http://www.continuummechanics.org/principalstrain.html
322 e1 = eps11(i)
323 e2 = eps22(i)
324 e3 = eps33(i)
325 e4 = half*eps12(i)
326 e5 = half*eps23(i)
327 e6 = half*eps31(i)
328 e42 = e4*e4
329 e52 = e5*e5
330 e62 = e6*e6
331
332 i1 = e1 + e2 + e3
333 i2 = e1*e2 + e2*e3 + e3*e1 - e4*e4 - e5*e5 - e6*e6
334 i3 = e1*e2*e3 - e1*e52 - e2*e62 - e3*e42 + two*e4*e5*e6
335
336 q = (three*i2 - i1*i1)/nine
337 r = (two*i1*i1*i1-nine*i1*i2+twenty7*i3)/cinquante4 ! (2*I3^3-9*I1*I2+27*I3)/54
338
339 r_inter = min(r/sqrt(max(em20,(-q**3))),one)
340 phi = acos(max(r_inter,-one))
341
342 e11 = two*sqrt(-q)*cos(phi/three)+third*i1
343 e22 = two*sqrt(-q)*cos((phi+two*pi)/three)+third*i1
344 e33 = two*sqrt(-q)*cos((phi+four*pi)/three)+third*i1
345c
346 IF (e11 < e22) THEN
347 r_inter = e11
348 e11 = e22
349 e22 = r_inter
350 ENDIF
351 IF (e22 < e33)THEN
352 r_inter = e22
353 e22 = e33
354 e33 = r_inter
355 ENDIF
356 IF (e11 < e22)THEN
357 r_inter = e11
358 e11 = e22
359 e22 = r_inter
360 ENDIF
361 dfmax(i) = min(one,(e11 / uvar(i,1)))
362c
363 IF(e11 >= uvar(i,1)) THEN
364 damage(i)= one
365 off(i)=four_over_5
366 nindxf=nindxf+1
367 indxf(nindxf)=i
368 tdele(i) = time
369#include "lockon.inc"
370 WRITE(iout, 6000) ngl(i),e11
371 WRITE(istdo,6100) ngl(i),e11,time
372#include "lockoff.inc"
373 ENDIF
374c
375 IF (uvar(i,1) < zero .AND. abs(e11) >= abs(uvar(i,1))) THEN
376 damage(i)= one
377 off(i)=four_over_5
378 nindxf=nindxf+1
379 indxf(nindxf)=i
380 tdele(i) = time
381#include "lockon.inc"
382 WRITE(iout, 6000) ngl(i),e11
383 WRITE(istdo,6100) ngl(i),e11,time
384#include "lockoff.inc"
385 ENDIF
386
387 IF(epsp2 > zero)THEN
388 IF(abs(e22) >= epsp2) THEN
389 damage(i)= one
390 off(i)=four_over_5
391 nindxf=nindxf+1
392 indxf(nindxf)=i
393 tdele(i) = time
394#include "lockon.inc"
395 WRITE(iout, 7000) ngl(i),e22
396 WRITE(istdo,7100) ngl(i),e22,time
397#include "lockoff.inc"
398 ENDIF
399 ENDIF
400 ENDIF ! EPSP1 > ZERO .OR. EPSP2 > ZERO
401c
402 ENDIF ! OFF(I) == ONE
403 ENDDO
404c
405c----------------------------------------------------
406 CASE (2) ! Equivalent strain criterion (new)
407c----------------------------------------------------
408 IF (strflag == 1) THEN ! Radioss strain is Cauchy-Green
409c transform total grad def strain to engineering
410 DO i=1,nel
411 eps11(i) = mfxx(i)
412 eps22(i) = mfyy(i)
413 eps33(i) = mfzz(i)
414 eps12(i) = mfxy(i) + mfyx(i)
415 eps23(i) = mfzy(i) + mfyz(i)
416 eps31(i) = mfxz(i) + mfzx(i)
417 END DO
418 ELSE IF (strflag == 2) THEN
419c transform engineering strain to true
420 DO i=1,nel
421 eps11(i) = log(epsxx(i) + one)
422 eps22(i) = log(epsyy(i) + one)
423 eps33(i) = log(epszz(i) + one)
424 eps12(i) = log(epsxy(i) + one)
425 eps23(i) = log(epsyz(i) + one)
426 eps31(i) = log(epszx(i) + one)
427 END DO
428 ELSE IF (strflag == 3) THEN
429c transform engineering strain to true
430 DO i=1,nel
431 eps11(i) = log(epsxx(i) + one)
432 eps22(i) = log(epsyy(i) + one)
433 eps33(i) = log(epszz(i) + one)
434 eps12(i) = log(epsxy(i) + one)
435 eps23(i) = log(epsyz(i) + one)
436 eps31(i) = log(epszx(i) + one)
437 END DO
438 ELSE IF (strflag == 4) THEN ! Radioss strain is Cauchy-Green
439c transform total grad def strain to true
440 DO i=1,nel
441 eps11(i) = log(mfxx(i) + one)
442 eps22(i) = log(mfyy(i) + one)
443 eps33(i) = log(mfzz(i) + one)
444 eps12(i) = log(mfxy(i) + mfyx(i) + one)
445 eps23(i) = log(mfxz(i) + mfzx(i) + one)
446 eps31(i) = log(mfzy(i) + mfyz(i) + one)
447 END DO
448 ELSE
449 eps11(1:nel) = epsxx(1:nel)
450 eps22(1:nel) = epsyy(1:nel)
451 eps33(1:nel) = epszz(1:nel)
452 eps12(1:nel) = epsxy(1:nel)
453 eps23(1:nel) = epsyz(1:nel)
454 eps31(1:nel) = epszx(1:nel)
455 END IF
456c calculate eps max
457 DO i=1,nel
458 IF (off(i) == one ) THEN
459 dav = (eps11(i)+eps22(i)+eps33(i))*third
460 e1 = eps11(i) - dav
461 e2 = eps22(i) - dav
462 e3 = eps33(i) - dav
463 e4 = half*eps12(i)
464 e5 = half*eps23(i)
465 e6 = half*eps31(i)
466 e42 = e4*e4
467 e52 = e5*e5
468 e62 = e6*e6
469 c = (e1*e2 + e1*e3 + e2*e3) - e42 - e52 - e62
470 d = - e1*e2*e3 + e1*e52 + e2*e62 + e3*e42 - two*e4*e5*e6
471 cc = c*third
472 epst = sqrt(-cc)
473 epst2 = epst * epst
474 y = (epst2 + c)* epst + d
475 y1 = -(epst2 + c)* epst + d
476 IF(abs(y) > em8 .OR.(abs(y) < em8 .AND. abs(y1) < abs(y)))THEN
477 epst = 1.75 * epst
478 epst2 = epst * epst
479 y = (epst2 + c)* epst + d
480 yp = three*epst2 + c
481 IF(yp /= zero)epst = epst - y/yp
482 epst2 = epst * epst
483 y = (epst2 + c)* epst + d
484 yp = three*epst2 + c
485 IF(yp /= zero)epst = epst - y/yp
486 epst2 = epst * epst
487 y = (epst2 + c)* epst + d
488 yp = three*epst2 + c
489 IF(yp /= zero)epst = epst - y/yp
490
491 epst2 = epst * epst
492 y = (epst2 + c)* epst + d
493 yp = three*epst2 + c
494 IF(yp /= zero)epst = epst - y/yp
495 ENDIF
496 epst = epst + dav
497 eps_max(i) = epst
498 END IF
499 ENDDO
500c
501c test failure crit
502c
503 DO i=1,nel
504 IF (off(i) == one)THEN
505 rfac = one
506 IF(ifunc(1) > 0)THEN
507 epsp_unit = epsp(i) * unit_t
508 rfac = finter(ifunc(1),epsp_unit,npf,tf,dydx)
509 rfac = max(rfac,em20)
510 ENDIF
511 IF (ifunc(3) > 0) THEN ! temperature
512 rfac2 = finter(ifunc(3),tstar(i),npf,tf,dydx)
513 rfac2 = max(rfac2,em20)
514 ELSE
515 rfac2 = one
516 ENDIF
517 r1 = epsf1*rfac*rfac2*uvar(i,1)
518 r2 = epsf2*rfac*rfac2*uvar(i,1)
519c
520 IF (eps_max(i) > r1 .AND. r1 < r2) THEN
521 damage(i)= (eps_max(i)-r1)/(r2-r1)
522 damage(i)= min(one,damage(i))
523 ENDIF
524
525 IF (eps_max(i) >= r2) THEN
526 damage(i)= one
527 off(i)=four_over_5
528 nindxf=nindxf+1
529 indxf(nindxf)=i
530 tdele(i) = time
531#include "lockon.inc"
532 WRITE(iout, 3000) ngl(i),eps_max(i)
533 WRITE(istdo,3100) ngl(i),eps_max(i),time
534#include "lockoff.inc"
535 ENDIF
536 ENDIF
537 ENDDO
538c-------------------
539 CASE (3) ! Max principal tensile strain. No failure in compression
540c-------------------
541 DO i=1,nel
542 IF (off(i) == one ) THEN
543
544 e1 = epsxx(i)
545 e2 = epsyy(i)
546 e3 = epszz(i)
547 e4 = half*epsxy(i)
548 e5 = half*epsyz(i)
549 e6 = half*epszx(i)
550 e42 = e4*e4
551 e52 = e5*e5
552 e62 = e6*e6
553
554 i1 = e1 + e2 + e3
555 i2 = e1*e2 + e2*e3 + e3*e1 - e4*e4 - e5*e5 - e6*e6
556 i3 = e1*e2*e3 - e1*e52 - e2*e62 - e3*e42 + two*e4*e5*e6
557 q = (three*i2 - i1*i1)/nine
558 r = (two*i1*i1*i1-nine*i1*i2+twenty7*i3)/cinquante4 ! (2*I3^3-9*I1*I2+27*I3)/54
559
560 r_inter = min(r/sqrt(max(em20,(-q**3))),one)
561 phi = acos(max(r_inter,-one))
562
563 e11 = two*sqrt(-q)*cos(phi/three)+third*i1
564 e22 = two*sqrt(-q)*cos((phi+two*pi)/three)+third*i1
565 e33 = two*sqrt(-q)*cos((phi+four*pi)/three)+third*i1
566c
567 IF (strflag == 1) THEN
568 e11 = sqrt(e11 + one) - one
569 e22 = sqrt(e22 + one) - one
570 e33 = sqrt(e33 + one) - one
571 ELSE IF (strflag == 2) THEN
572 e11 = exp(e11) - one
573 e22 = exp(e22) - one
574 e33 = exp(e33) - one
575 ELSE IF (strflag == 3) THEN
576 e11 = log(e11 + one)
577 e22 = log(e22 + one)
578 e33 = log(e33 + one)
579 ELSE IF (strflag == 4) THEN
580 e11 = log(sqrt(e11+one))
581 e22 = log(sqrt(e22+one))
582 e33 = log(sqrt(e33+one))
583 END IF
584c
585 IF (e11 < e22) THEN
586 r_inter = e11
587 e11 = e22
588 e22 = r_inter
589 ENDIF
590 IF (e22 < e33)THEN
591 r_inter = e22
592 e22 = e33
593 e33 = r_inter
594 ENDIF
595 IF (e11 < e22)THEN
596 r_inter = e11
597 e11 = e22
598 e22 = r_inter
599 ENDIF
600 eps_max(i) = e11
601 END IF
602 ENDDO
603c test failure crit
604 DO i=1,nel
605 IF (off(i) == one ) THEN
606 r1 = epsf1*uvar(i,1)
607 r2 = epsf2*uvar(i,1)
608 IF (ifunc(1) > 0) THEN ! strain rate factor
609 epsp_unit = epsp(i) * unit_t
610 rfac = finter(ifunc(1),epsp_unit,npf,tf,dydx)
611 rfac = max(rfac,em20)
612 r1 = r1*rfac
613 r2 = r2*rfac
614 ENDIF
615 IF (ifunc(3) > 0) THEN ! temperature factor
616 rfac2 = finter(ifunc(3),tstar(i),npf,tf,dydx)
617 rfac2 = max(rfac2,em20)
618 r1 = r1*rfac2
619 r2 = r2*rfac2
620 ENDIF
621c
622 IF (eps_max(i) > r1 .AND. r1 < r2 .AND. eps_max(i) > zero) THEN
623 IF (dfmax(i) == zero) THEN
624 nindxd = nindxd + 1
625 indxd(nindxd) = i
626#include "lockon.inc"
627 WRITE(iout, 2001) ngl(i),eps_max(i)
628 WRITE(istdo,2001) ngl(i),eps_max(i)
629#include "lockoff.inc"
630 ENDIF
631 damage(i) = (eps_max(i)-r1)/(r2-r1)
632 damage(i) = min(one,damage(i))
633 ENDIF
634 IF (eps_max(i) >= r2 .AND. eps_max(i) > zero) THEN
635 damage(i)= one
636 off(i) = four_over_5
637 nindxf=nindxf+1
638 indxf(nindxf)=i
639 tdele(i) = time
640#include "lockon.inc"
641 WRITE(iout, 6000) ngl(i),eps_max(i)
642 WRITE(istdo,6100) ngl(i),eps_max(i),time
643#include "lockoff.inc"
644 ENDIF
645 ENDIF
646c
647 ENDDO
648c---------------
649 END SELECT
650c--------------------------
651c
652 IF (nindxf > 0 .AND. imconv == 1) THEN
653 DO j=1,nindxf
654#include "lockon.inc"
655 WRITE(iout, 1000) ngl(indxf(j))
656 WRITE(istdo,1100) ngl(indxf(j)),time
657#include "lockoff.inc"
658 END DO
659 END IF
660C ---
661 DO i=1,nel
662 r1 = epsf1
663 r2 = epsf2
664 IF (r1 < r2) THEN
665 dmg_scale(i) = one - damage(i)
666 END IF
667 dfmax(i) = max(dfmax(i),damage(i)) ! Maximum Damage storing for output : 0 < DFMAX < 1
668 ENDDO
669c-----------------------------------------------
670 1000 FORMAT(1x,'DELETE SOLID ELEMENT NUMBER ',i10)
671 1100 FORMAT(1x,'DELETE SOLID ELEMENT NUMBER ',i10,' AT TIME :',1pe20.13)
672 2001 FORMAT(1x,'START DAMAGE (TENS) OF ELEMENT ',i10,', 1st PRINCIPAL STRAIN = ',g11.4)
673 2002 FORMAT(1x,'START DAMAGE (TENS) OF ELEMENT ',i10,', EQUIVALENT STRAIN = ',g11.4)
674 3000 FORMAT(1x,'FAILURE (TENS) OF ELEMENT ',i10,', MAX EQUIVALENT STRAIN =',g11.4)
675 3100 FORMAT(1x,'FAILURE (TENS) OF ELEMENT ',i10,', MAX EQUIVALENT STRAIN =',g11.4,
676 . ' AT TIME :',1pe12.4)
677 6000 FORMAT(1x,'FAILURE (TENS) OF ELEMENT ',i10,', 1st PRINCIPAL STRAIN = ',g11.4)
678 6100 FORMAT(1x,'FAILURE (TENS) OF ELEMENT ',i10,', 1st PRINCIPAL STRAIN = ',g11.4,
679 . ' AT TIME :',1pe12.4)
680 7000 FORMAT(1x,'FAILURE (TENS) OF ELEMENT ',i10,', 2nd PRINCIPAL STRAIN = ',g11.4)
681 7100 FORMAT(1x,'FAILURE (TENS) OF ELEMENT ',i10,', 2nd PRINCIPAL STRAIN = ',g11.4,
682 . ' AT TIME :',1pe12.4)
683c-----------------------------------------------
684 RETURN
685 END
fail_tensstrain_s
subroutine fail_tensstrain_s(nel, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, ngl, aldt, tstar, ismstr, epsxx, epsyy, epszz, epsxy, epsyz, epszx, signxx, signyy, signzz, signxy, signyz, signzx, epsp, uvar, off, dfmax, tdele, mfxx, mfxy, mfxz, mfyx, mfyy, mfyz, mfzx, mfzy, mfzz, dmg_scale)
Definition fail_tensstrain_s.F:42
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21