OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_tensstrain_s.F File Reference
#include "implicit_f.inc"
#include "scr17_c.inc"
#include "units_c.inc"
#include "comlock.inc"
#include "param_c.inc"
#include "impl1_c.inc"

Go to the source code of this file.

Functions/Subroutines

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)

Function/Subroutine Documentation

◆ fail_tensstrain_s()

subroutine fail_tensstrain_s ( integer nel,
integer nuparam,
integer nuvar,
integer nfunc,
integer, dimension(nfunc) ifunc,
integer, dimension(*) npf,
tf,
time,
timestep,
dimension(nuparam) uparam,
integer, dimension(nel) ngl,
dimension(nel) aldt,
dimension(nel) tstar,
integer ismstr,
dimension(nel) epsxx,
dimension(nel) epsyy,
dimension(nel) epszz,
dimension(nel) epsxy,
dimension(nel) epsyz,
dimension(nel) epszx,
dimension(nel) signxx,
dimension(nel) signyy,
dimension(nel) signzz,
dimension(nel) signxy,
dimension(nel) signyz,
dimension(nel) signzx,
dimension(nel) epsp,
uvar,
off,
dfmax,
tdele,
dimension(nel) mfxx,
dimension(nel) mfxy,
dimension(nel) mfxz,
dimension(nel) mfyx,
dimension(nel) mfyy,
dimension(nel) mfyz,
dimension(nel) mfzx,
dimension(nel) mfzy,
dimension(nel) mfzz,
intent(inout) dmg_scale )

Definition at line 33 of file fail_tensstrain_s.F.

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
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21