OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_gene1_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_gene1_s ../engine/source/materials/fail/gene1/fail_gene1_s.f
25!||--- called by ------------------------------------------------------
26!|| mmain ../engine/source/materials/mat_share/mmain.F90
27!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
28!|| usermat_solid ../engine/source/materials/mat_share/usermat_solid.F
29!||--- calls -----------------------------------------------------
30!|| finter ../engine/source/tools/curve/finter.f
31!|| table2d_vinterp_log ../engine/source/tools/curve/table2d_vinterp_log.f
32!|| table_vinterp ../engine/source/tools/curve/table_tools.F
33!||--- uses -----------------------------------------------------
34!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
35!|| interface_table_mod ../engine/share/modules/table_mod.F
36!|| table_mod ../engine/share/modules/table_mod.F
37!||====================================================================
38 SUBROUTINE fail_gene1_s (
39 1 NEL ,NUPARAM ,NUVAR ,NFUNC ,IFUNC ,LOFF ,
40 2 NPF ,TF ,TIME ,TIMESTEP ,UPARAM ,IPG ,
41 3 NGL ,DT ,EPSP ,UVAR ,OFF ,NPG ,
42 4 EPSXX ,EPSYY ,EPSZZ ,EPSXY ,EPSYZ ,EPSZX ,
43 5 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
44 6 TEMP ,VOLN ,DFMAX ,TDELE ,ALDT ,TABLE ,
45 7 IRUPT ,ELBUF_TAB,ILAY1 ,NTABLF ,ITABLF ,LF_DAMMX ,
46 8 NIPARAM ,IPARAM )
47C-----------------------------------------------
48C M o d u l e s
49C-----------------------------------------------
50 USE table_mod
51 USE interface_table_mod
52 USE elbufdef_mod
53C!-----------------------------------------------
54C! I m p l i c i t T y p e s
55C!-----------------------------------------------
56#include "implicit_f.inc"
57C-----------------------------------------------
58C C o m m o n B l o c k s
59C-----------------------------------------------
60#include "scr17_c.inc"
61#include "units_c.inc"
62#include "comlock.inc"
63#include "com04_c.inc"
64#include "com01_c.inc"
65#include "tabsiz_c.inc"
66C!-----------------------------------------------
67 INTEGER, INTENT(IN) :: NEL,NUPARAM,NUVAR,IPG,NPG,NFUNC,NTABLF,
68 . NIPARAM,LF_DAMMX,IRUPT,ILAY1
69 INTEGER, DIMENSION(NFUNC) ,INTENT(IN) :: IFUNC
70 INTEGER, DIMENSION(NTABLF) ,INTENT(IN) :: ITABLF
71 INTEGER, DIMENSION(NEL) ,INTENT(IN) :: NGL
72 INTEGER, DIMENSION(NIPARAM),INTENT(IN) :: IPARAM
73 my_real, INTENT(IN) :: TIME,TIMESTEP
74 my_real, DIMENSION(NUPARAM), INTENT(IN) :: UPARAM
75 my_real, DIMENSION(NEL), INTENT(IN) :: epsxx,epsyy,epszz,
76 . epsxy,epsyz,epszx,dt,epsp,aldt,temp,voln
77 my_real, DIMENSION(NEL), INTENT(INOUT) :: signxx,signyy,signzz,
78 . signxy,signyz,signzx,off,loff,tdele
79 my_real, DIMENSION(NEL,LF_DAMMX), INTENT(INOUT) :: dfmax
80 my_real, DIMENSION(NEL,NUVAR), INTENT(INOUT) :: uvar
81 TYPE (TTABLE), DIMENSION(NTABLE) :: TABLE
82 TYPE (ELBUF_STRUCT_), TARGET :: ELBUF_TAB
83C!-----------------------------------------------
84C! VARIABLES FOR FUNCTION INTERPOLATION
85C!-----------------------------------------------
86 INTEGER, INTENT(IN) :: NPF(SNPC)
87 my_real, INTENT(IN) :: tf(stf)
88 my_real finter
89 EXTERNAL finter
90C!-----------------------------------------------
91C! L o c a l V a r i a b l e s
92C!-----------------------------------------------
93 INTEGER I,K,J,INDX(NEL),NINDX,INDX2(NEL),NINDX2,NSTEP,CRIT,NMOD,
94 . fct_ism,fct_ips,fct_idg12,fct_idg13,fct_ide1c,fct_idel,
95 . ncrit(nel),ipos(nel,2),ismooth,istrain,ir,is,it,
96 . ilay,tab_idfld,itab,ncs,indx3(nel),nindx3,ipmax(nel),ipmin(nel),
97 . is1max(nel),itmax(nel),imindt(nel),isigmax(nel),isigth(nel),
98 . iepsmax(nel),ieffeps(nel),ivoleps(nel),imineps(nel),ishear(nel),
99 . imix12(nel),imix13(nel),imxe1c(nel),ifld(nel),ithin(nel),
100 . imaxtemp(nel)
101 my_real
102 . minpres,maxpres,sigp1,tmax,dtmin,epsdot_sm,sigvm,sigth,
103 . kf,epsdot_ps,maxeps,effeps,voleps,mineps,epssh,epsdot_fld,
104 . thin,volfrac,maxtemp,fscale_el,el_ref,lambda,fac,df
105 my_real
106 . e1,e2,e3,e4,e5,e6,e42,e52,e62,i1,i2,i3,p(nel),sxx,syy,szz,svm(nel),
107 . q,r,r_inter,phi,e11(nel),e22(nel),e33(nel),vol_strain(nel),dav,e1d,
108 . e2d,e3d,e4d,e5d,e6d,s11(nel),s22(nel),s33(nel),eff_strain(nel),psi,
109 . epsmax(nel),sigmax(nel),facl(nel),sh12(nel),sh13(nel),e1c(nel),
110 . xvec(nel,2),e1fld(nel),dfld(nel),hardr(nel),vtot,vdam,denom,vfail(nel),
111 . triax(nel)
112 TYPE(buf_fail_) ,POINTER :: FBUF
113C!--------------------------------------------------------------
114 !=======================================================================
115 ! - INITIALISATION OF COMPUTATION ON TIME STEP
116 !=======================================================================
117 ! Recovering failure criterion parameters
118 ! -> Integer parameter, activated criteria
119 crit = iparam(1)
120 itab = iparam(2) !> Table dependency type (used only if tab_IDfld is a table).
121 nstep = iparam(3) !> Number of cycles for the stress reduction, Default = 10 (Integer)
122 ncs = iparam(4) !> Number of conditions to reach before the element is deleted, Default = 1 (Integer)
123 ismooth = iparam(5) !> Interpolation type (in case of tabulated yield function)
124 istrain = iparam(6) !> Engineering / True input strain flag
125 ! -> Real parameters
126 minpres = uparam(1) !> Minimum pressure (positive in compression)
127 maxpres = uparam(2) !> Maximum pressure (positive in compression)
128 sigp1 = uparam(3) !> Maximum principal stress
129 tmax = uparam(4) !> Failure time, Default = 1E+20 (Real)
130 dtmin = uparam(5) !> Minimum time step
131 epsdot_sm = uparam(6) !> Reference strain rate value for fct_IDsm, Default = 1 (Real)
132 sigvm = uparam(7) !> Von Mises stress
133 sigth = uparam(8)
134 kf = uparam(9) ! -> Tuler-Butcher
135 epsdot_ps = uparam(10) !> Reference strain rate value for fct_IDps, Default = 1 (Real)
136 maxeps = uparam(11) !> Ordinate scale factor for fct_IDps or maximum principal strain if fct_IDps is not defined, Default = 1
137 effeps = uparam(12) !> Maximum effective strain
138 voleps = uparam(13) !> Maximum volumetric strain
139 mineps = uparam(14) !> Minimum principal strain
140 epssh = uparam(15)
141 epsdot_fld = uparam(16)
142 thin = uparam(17) !Thinning failure value (Real)
143 volfrac = uparam(18) !> Damaged volume fraction to reach before the element is deleted (fully-integrated and higher order solid elements only), Default = 0.5 (Real)
144 maxtemp = uparam(20) !> Maximum temperature.
145 fscale_el = uparam(21) !> Element size function scale factor for fct_IDel, tab_IDfld (Itab=2), fct_IDg12, fct_IDg23, fct_IDg13 and fct_IDe1c, Default = 1.0 (Real)
146 el_ref = uparam(22) !> Reference element size for fct_IDel, tab_IDfld (Itab=2), fct_IDg12, fct_IDg23, fct_IDg13 and fct_IDe1c, Default = 1.0 (Real)
147c
148c function & tables
149 fct_ism = ifunc(1) !> Function identifier of the maximum equivalent stress versus strain rate
150 fct_ips = ifunc(2) !> maximum principal strain VS strain-rate
151 fct_idg12 = ifunc(3) !> in-plane shear strain VS element size
152 fct_idg13 = ifunc(4) !> in-plane shear strain VS element size
153 fct_ide1c = ifunc(5) !> major in plane-strain VS element size
154 fct_idel = ifunc(6) !> element size regularization
155 IF (ntablf > 0) THEN
156 tab_idfld = itablf(1)
157 ELSE
158 tab_idfld = 0
159 END IF
160c
161 ! Initialization of variable
162 nindx = 0
163 nindx2 = 0
164 nindx3 = 0
165 indx(1:nel) = 0
166 indx2(1:nel) = 0
167 indx3(1:nel) = 0
168 vfail(1:nel) = zero
169 ipmax(1:nel) = 0
170 ipmin(1:nel) = 0
171 is1max(1:nel) = 0
172 itmax(1:nel) = 0
173 imindt(1:nel) = 0
174 isigmax(1:nel) = 0
175 isigth(1:nel) = 0
176 iepsmax(1:nel) = 0
177 ieffeps(1:nel) = 0
178 ivoleps(1:nel) = 0
179 imineps(1:nel) = 0
180 ishear(1:nel) = 0
181 imix12(1:nel) = 0
182 imix13(1:nel) = 0
183 imxe1c(1:nel) = 0
184 ifld(1:nel) = 0
185 ithin(1:nel) = 0
186 imaxtemp(1:nel) = 0
187 ncrit(1:nel) = 0
188c
189 ! At initial time, compute the element size regularization factor
190 IF (uvar(1,1)==zero) THEN
191 IF (fct_idel > 0) THEN
192 DO i=1,nel
193 lambda = aldt(i)/el_ref
194 fac = finter(fct_idel,lambda,npf,tf,df)
195 uvar(i,1) = fac*fscale_el
196 ENDDO
197 ELSE
198 uvar(1:nel,1) = one
199 ENDIF
200 ENDIF
201 IF (uvar(1,5) == zero.AND.(off(1) /= zero)) uvar(1:nel,5) = one
202 IF (uvar(1,8) == zero) uvar(1:nel,8) = aldt(1:nel)
203c
204 ! Checking element failure and recovering user variable
205 DO i=1,nel
206 ! Integration point failure
207 IF (uvar(i,5) < one .AND. uvar(i,5) >= em08) THEN
208 uvar(i,5) = uvar(i,5) - one/nstep
209 ENDIF
210 IF (uvar(i,5) <= em08) uvar(i,5) = zero
211 signxx(i) = signxx(i)*uvar(i,5)
212 signyy(i) = signyy(i)*uvar(i,5)
213 signzz(i) = signzz(i)*uvar(i,5)
214 signxy(i) = signxy(i)*uvar(i,5)
215 signyz(i) = signyz(i)*uvar(i,5)
216 signzx(i) = signzx(i)*uvar(i,5)
217 ! Regularization factors for length, surface and volume
218 facl(i) = uvar(i,1)
219 ! Storage of integration point volume
220 IF (npg > 1) THEN
221 uvar(i,3) = voln(i)
222 IF (uvar(i,5) == zero) uvar(i,4) = voln(i)
223 ENDIF
224 ENDDO
225c
226 !====================================================================
227 ! - LOOP OVER THE ELEMENT TO COMPUTE THE STRESSES AND STRAINS
228 !====================================================================
229 DO i=1,nel
230c
231 ! For active element or Gauss point
232 IF ((uvar(i,5) == one).AND.(off(i)==one)) THEN
233 ! ----------------------------------------------------------------------------------------
234 ! Computation of volumetric strain, effective strain, shear strain and principal strains
235 ! ----------------------------------------------------------------------------------------
236 ! -> computation of tensiorial strain
237 e1 = epsxx(i)
238 e2 = epsyy(i)
239 e3 = epszz(i)
240 e4 = half*epsxy(i)
241 e5 = half*epsyz(i)
242 e6 = half*epszx(i)
243 ! -> computation of strain tensor invariants
244 e42 = e4*e4
245 e52 = e5*e5
246 e62 = e6*e6
247 i1 = e1 + e2 + e3
248 i2 = e1*e2 + e2*e3 + e3*e1 - e4*e4 - e5*e5 - e6*e6
249 i3 = e1*e2*e3 - e1*e52 - e2*e62 - e3*e42 + two*e4*e5*e6
250 ! -> computation of principal strains
251 q = (three*i2 - i1*i1)/nine
252 r = (two*i1*i1*i1-nine*i1*i2+twenty7*i3)/cinquante4 ! (2*I3^3-9*I1*I2+27*I3)/54
253 r_inter = min(r/sqrt(max(em20,(-q**3))),one)
254 phi = acos(max(r_inter,-one))
255 e11(i) = two*sqrt(max(zero,-q))*cos(phi/three)+third*i1
256 e22(i) = two*sqrt(max(zero,-q))*cos((phi+two*pi)/three)+third*i1
257 e33(i) = two*sqrt(max(zero,-q))*cos((phi+four*pi)/three)+third*i1
258 IF (e11(i) < e22(i)) THEN
259 r_inter = e11(i)
260 e11(i) = e22(i)
261 e22(i) = r_inter
262 ENDIF
263 IF (e22(i) < e33(i))THEN
264 r_inter = e22(i)
265 e22(i) = e33(i)
266 e33(i) = r_inter
267 ENDIF
268 IF (e11(i) < e22(i))THEN
269 r_inter = e11(i)
270 e11(i) = e22(i)
271 e22(i) = r_inter
272 ENDIF
273 ! -> computation of volumetric strain
274 vol_strain(i) = e11(i) + e22(i) + e33(i)
275c
276 ! -> computation of effective strain
277 dav = (epsxx(i)+epsyy(i)+epszz(i))*third
278 e1d = epsxx(i) - dav
279 e2d = epsyy(i) - dav
280 e3d = epszz(i) - dav
281 e4d = half*epsxy(i)
282 e5d = half*epsyz(i)
283 e6d = half*epszx(i)
284 eff_strain(i) = e1d**2 + e2d**2 + e3d**3 + two*(e4d**2 + e5d**2 + e6d**2)
285 eff_strain(i) = sqrt(two_third*eff_strain(i))
286c
287 ! --------------------------------------------------------------------------
288 ! Computation of hydrostatic stress, Von Mises stress and principal stresses
289 ! --------------------------------------------------------------------------
290 ! -> pressure stress (positive in compression)
291 p(i) = -third*(signxx(i) + signyy(i) + signzz(i))
292 ! -> equivalent stress of Von Mises
293 sxx = signxx(i) + p(i)
294 syy = signyy(i) + p(i)
295 szz = signzz(i) + p(i)
296 svm(i) = half*(sxx**2 + syy**2 + szz**2)
297 . + signxy(i)**2 + signzx(i)**2 + signyz(i)**2
298 svm(i) = sqrt(three*svm(i))
299 triax(i) = -p(i)/max(svm(i),em20)
300 ! -> computing the principal stresses
301 i1 = signxx(i)+signyy(i)+signzz(i)
302 i2 = signxx(i)*signyy(i)+signyy(i)*signzz(i)+signzz(i)*signxx(i)-
303 . signxy(i)*signxy(i)-signzx(i)*signzx(i)-signyz(i)*signyz(i)
304 i3 = signxx(i)*signyy(i)*signzz(i)-signxx(i)*signyz(i)*signyz(i)-
305 . signyy(i)*signzx(i)*signzx(i)-signzz(i)*signxy(i)*signxy(i)+
306 . two*signxy(i)*signzx(i)*signyz(i)
307 q = (three*i2 - i1*i1)/nine
308 r = (two*i1*i1*i1-nine*i1*i2+twenty7*i3)/cinquante4 ! (2*I3^3-9*I1*I2+27*I3)/54
309 r_inter = min(r/sqrt(max(em20,(-q**3))),one)
310 psi = acos(max(r_inter,-one))
311 s11(i) = two*sqrt(max(zero,-q))*cos(psi/three)+third*i1
312 s22(i) = two*sqrt(max(zero,-q))*cos((psi+two*pi)/three)+third*i1
313 s33(i) = two*sqrt(max(zero,-q))*cos((psi+four*pi)/three)+third*i1
314 IF (s11(i) < s22(i)) THEN
315 r_inter = s11(i)
316 s11(i) = s22(i)
317 s22(i) = r_inter
318 ENDIF
319 IF (s22(i) < s33(i)) THEN
320 r_inter = s22(i)
321 s22(i) = s33(i)
322 s33(i) = r_inter
323 ENDIF
324 IF (s11(i) < s22(i)) THEN
325 r_inter = s11(i)
326 s11(i) = s22(i)
327 s22(i) = r_inter
328 ENDIF
329c
330 ! For broken element or Gauss point
331 ELSE
332 e11(i) = zero
333 e22(i) = zero
334 e33(i) = zero
335 vol_strain(i) = zero
336 eff_strain(i) = zero
337 p(i) = zero
338 svm(i) = zero
339 triax(i) = zero
340 s11(i) = zero
341 s22(i) = zero
342 s33(i) = zero
343 ENDIF
344c
345 ENDDO
346c
347 ! -> Forming limit diagram
348 IF (ntablf > 0) THEN
349 IF (itab == 1) THEN
350 ! Diagram using true strains
351 IF (istrain == 0) THEN
352 ! In-plane tabulation with strain-rate
353 xvec(1:nel,1) = e22(1:nel)
354 xvec(1:nel,2) = epsp(1:nel)/epsdot_fld
355 ! -> Tensile yield stress in direction 1 (MD)
356 ipos(1:nel,1:2) = 1
357 CALL table2d_vinterp_log(table(tab_idfld),ismooth,nel,nel,ipos,xvec,e1fld,dfld,hardr)
358 ! Diagram using engineering strain
359 ELSE
360 ! In-plane tabulation with strain-rate
361 xvec(1:nel,1) = exp(e22(1:nel))-one
362 xvec(1:nel,2) = epsp(1:nel)/epsdot_fld
363 ! -> Tensile yield stress in direction 1 (MD)
364 ipos(1:nel,1:2) = 1
365 CALL table2d_vinterp_log(table(tab_idfld),ismooth,nel,nel,ipos,xvec,e1fld,dfld,hardr)
366 e1fld = log(one + e1fld)
367 ENDIF
368 ELSE
369 ! Diagram using true strains
370 IF (istrain == 0) THEN
371 ! In-plane tabulation with strain-rate
372 xvec(1:nel,1) = e22(1:nel)
373 xvec(1:nel,2) = aldt(1:nel)/el_ref
374 ! -> Tensile yield stress in direction 1 (MD)
375 ipos(1:nel,1:2) = 1
376 CALL table_vinterp(table(tab_idfld),nel,nel,ipos,xvec,e1fld,dfld)
377 ! Diagram using engineering strains
378 ELSE
379 ! In-plane tabulation with strain-rate
380 xvec(1:nel,1) = exp(e22(1:nel))-one
381 xvec(1:nel,2) = aldt(1:nel)/el_ref
382 ! -> Tensile yield stress in direction 1 (MD)
383 ipos(1:nel,1:2) = 1
384 CALL table_vinterp(table(tab_idfld),nel,nel,ipos,xvec,e1fld,dfld)
385 e1fld = log(one + e1fld)
386 ENDIF
387 ENDIF
388 ENDIF
389c
390 !====================================================================
391 ! - LOOP OVER THE ELEMENT TO CHECK THE EROSION CRITERIA
392 !====================================================================
393 DO i = 1,nel
394 nmod = 0
395 IF ((uvar(i,5) == one).AND.(off(i)==one)) THEN
396 ! -> minimum pressure
397 IF (btest(crit,1)) THEN
398 nmod = nmod + 1
399 dfmax(i,1+nmod) = max(p(i)/(minpres*facl(i)),dfmax(i,1+nmod))
400 dfmax(i,1+nmod) = min(dfmax(i,1+nmod),one)
401 IF (p(i) <= minpres*facl(i)) THEN
402 ncrit(i) = ncrit(i) + 1
403 ipmin(i) = 1
404 ENDIF
405 ENDIF
406 ! -> maximum pressure
407 IF (btest(crit,2)) THEN
408 nmod = nmod + 1
409 dfmax(i,1+nmod) = max(p(i)/(maxpres*facl(i)),dfmax(i,1+nmod))
410 dfmax(i,1+nmod) = min(dfmax(i,1+nmod),one)
411 IF (p(i) >= maxpres*facl(i)) THEN
412 ncrit(i) = ncrit(i) + 1
413 ipmax(i) = 1
414 ENDIF
415 ENDIF
416 ! -> maximal principal stress
417 IF (btest(crit,3)) THEN
418 nmod = nmod + 1
419 ! (unrestricted)
420 IF (sigp1 > zero) THEN
421 dfmax(i,1+nmod) = max(s11(i)/(sigp1*facl(i)),dfmax(i,1+nmod))
422 dfmax(i,1+nmod) = min(dfmax(i,1+nmod),one)
423 IF (s11(i) >= sigp1*facl(i)) THEN
424 ncrit(i) = ncrit(i) + 1
425 is1max(i) = 1
426 ENDIF
427 ! (restricted to positive stress triaxialities)
428 ELSE
429 IF (triax(i)>em10) THEN
430 dfmax(i,1+nmod) = max(s11(i)/(abs(sigp1)*facl(i)),dfmax(i,1+nmod))
431 dfmax(i,1+nmod) = min(dfmax(i,1+nmod),one)
432 IF (s11(i) >= abs(sigp1)*facl(i)) THEN
433 ncrit(i) = ncrit(i) + 1
434 is1max(i) = 1
435 ENDIF
436 ENDIF
437 ENDIF
438 ENDIF
439 ! -> maximum time
440 IF (btest(crit,4)) THEN
441 nmod = nmod + 1
442 dfmax(i,1+nmod) = max(time/tmax,dfmax(i,1+nmod))
443 dfmax(i,1+nmod) = min(dfmax(i,1+nmod),one)
444 IF (time >= tmax) THEN
445 ncrit(i) = ncrit(i) + 1
446 itmax(i) = 1
447 ENDIF
448 ENDIF
449 ! -> minimum timestep
450 IF (btest(crit,5)) THEN
451 nmod = nmod + 1
452 IF (time > zero) THEN
453 dfmax(i,1+nmod) = max(dtmin/dt(i),dfmax(i,1+nmod))
454 dfmax(i,1+nmod) = min(dfmax(i,1+nmod),one)
455 IF ((dt(i) <= dtmin).AND.(time > zero)) THEN
456 ncrit(i) = ncrit(i) + 1
457 imindt(i) = 1
458 ENDIF
459 ENDIF
460 ENDIF
461 ! -> equivalent stress
462 IF (btest(crit,6)) THEN
463 nmod = nmod + 1
464 IF (epsdot_sm /= zero) THEN
465 lambda = epsp(i)/epsdot_sm
466 sigmax(i) = finter(fct_ism,lambda,npf,tf,df)
467 sigmax(i) = sigmax(i)*sigvm
468 ELSE
469 sigmax(i) = sigvm
470 ENDIF
471 dfmax(i,1+nmod) = max(svm(i)/(sigmax(i)*facl(i)),dfmax(i,1+nmod))
472 dfmax(i,1+nmod) = min(dfmax(i,1+nmod),one)
473 IF (svm(i) >= sigmax(i)*facl(i)) THEN
474 ncrit(i) = ncrit(i) + 1
475 isigmax(i) = 1
476 ENDIF
477 ENDIF
478 ! -> Tuler-Butcher
479 IF (btest(crit,7)) THEN
480 nmod = nmod + 1
481 dfmax(i,1+nmod) = max(uvar(i,2)/(kf*facl(i)),dfmax(i,1+nmod))
482 dfmax(i,1+nmod) = min(dfmax(i,1+nmod),one)
483 IF (s11(i) > sigth) THEN
484 uvar(i,2) = uvar(i,2) + timestep*(s11(i) - sigth)**2
485 IF (uvar(i,2) >= kf*facl(i)) THEN
486 ncrit(i) = ncrit(i) + 1
487 isigth(i) = 1
488 ENDIF
489 ENDIF
490 ENDIF
491 ! -> maximal principal strain
492 IF (btest(crit,8)) THEN
493 nmod = nmod + 1
494 IF (epsdot_ps /= zero) THEN
495 lambda = epsp(i)/epsdot_ps
496 epsmax(i) = finter(fct_ips,lambda,npf,tf,df)
497 epsmax(i) = epsmax(i)*maxeps
498 ELSE
499 epsmax(i) = maxeps
500 ENDIF
501 dfmax(i,1+nmod) = max(e11(i)/(epsmax(i)*facl(i)),dfmax(i,1+nmod))
502 dfmax(i,1+nmod) = min(dfmax(i,1+nmod),one)
503 IF (e11(i) >= epsmax(i)*facl(i)) THEN
504 ncrit(i) = ncrit(i) + 1
505 iepsmax(i) = 1
506 ENDIF
507 ENDIF
508 ! -> maximum effective strain
509 IF (btest(crit,9)) THEN
510 nmod = nmod + 1
511 dfmax(i,1+nmod) = max(eff_strain(i)/(effeps*facl(i)),dfmax(i,1+nmod))
512 dfmax(i,1+nmod) = min(dfmax(i,1+nmod),one)
513 IF (eff_strain(i) >= effeps*facl(i)) THEN
514 ncrit(i) = ncrit(i) + 1
515 ieffeps(i) = 1
516 ENDIF
517 ENDIF
518 ! -> maximum volumetric strain
519 IF (btest(crit,10)) THEN
520 nmod = nmod + 1
521 IF (voleps > zero) THEN
522 dfmax(i,1+nmod) = max(vol_strain(i)/(voleps*facl(i)),dfmax(i,1+nmod))
523 dfmax(i,1+nmod) = min(dfmax(i,1+nmod),one)
524 IF (vol_strain(i) >= voleps*facl(i)) THEN
525 ncrit(i) = ncrit(i) + 1
526 ivoleps(i) = 1
527 ENDIF
528 ELSE
529 dfmax(i,1+nmod) = max(vol_strain(i)/(voleps*facl(i)),dfmax(i,1+nmod))
530 dfmax(i,1+nmod) = min(dfmax(i,1+nmod),one)
531 IF (vol_strain(i) <= voleps*facl(i)) THEN
532 ncrit(i) = ncrit(i) + 1
533 ivoleps(i) = 1
534 ENDIF
535 ENDIF
536 ENDIF
537 ! -> minimum principal strain
538 IF (btest(crit,11)) THEN
539 nmod = nmod + 1
540 IF (e33(i) /= zero) THEN
541 dfmax(i,1+nmod) = max(mineps*facl(i)/(e33(i)),dfmax(i,1+nmod))
542 ENDIF
543 dfmax(i,1+nmod) = min(dfmax(i,1+nmod),one)
544 IF (e33(i) <= mineps*facl(i)) THEN
545 ncrit(i) = ncrit(i) + 1
546 imineps(i) = 1
547 ENDIF
548 ENDIF
549 ! -> maximum tensorial shear strain
550 IF (btest(crit,12)) THEN
551 nmod = nmod + 1
552 dfmax(i,1+nmod) = max(((e11(i) - e33(i))/two)/(epssh*facl(i)),dfmax(i,1+nmod))
553 dfmax(i,1+nmod) = min(dfmax(i,1+nmod),one)
554 IF ((e11(i) - e33(i))/two >= epssh*facl(i)) THEN
555 ncrit(i) = ncrit(i) + 1
556 ishear(i) = 1
557 ENDIF
558 ENDIF
559 ! -> mixed mode
560 IF (btest(crit,13)) THEN
561 lambda = uvar(i,8)/el_ref
562 sh12(i) = finter(fct_idg12,lambda,npf,tf,df)
563 denom = sign(max(abs(e11(i)),em20),e11(i))
564 nmod = nmod + 1
565 IF (((e22(i)/denom)<=-half).AND.((e22(i)/denom)>=-two)) THEN
566 dfmax(i,1+nmod) = max(((e11(i) - e22(i))/two)/(sh12(i)),dfmax(i,1+nmod))
567 dfmax(i,1+nmod) = min(dfmax(i,1+nmod),one)
568 IF ((e11(i) - e22(i))/two >= sh12(i)) THEN
569 ncrit(i) = ncrit(i) + 1
570 imix12(i) = 1
571 ENDIF
572 ENDIF
573 ENDIF
574 IF (btest(crit,14)) THEN
575 lambda = uvar(i,8)/el_ref
576 sh13(i) = finter(fct_idg13,lambda,npf,tf,df)
577 denom = sign(max(abs(e11(i)),em20),e11(i))
578 nmod = nmod + 1
579 IF (((e22(i)/denom)<=one).AND.((e22(i)/denom)>=-half)) THEN
580 dfmax(i,1+nmod) = max(((e11(i) - e33(i))/two)/(sh13(i)),dfmax(i,1+nmod))
581 dfmax(i,1+nmod) = min(dfmax(i,1+nmod),one)
582 IF ((e11(i) - e33(i))/two >= sh13(i)) THEN
583 ncrit(i) = ncrit(i) + 1
584 imix13(i) = 1
585 ENDIF
586 ENDIF
587 ENDIF
588 IF (btest(crit,15)) THEN
589 lambda = uvar(i,8)/el_ref
590 e1c(i) = finter(fct_ide1c,lambda,npf,tf,df)
591 denom = sign(max(abs(e11(i)),em20),e11(i))
592 nmod = nmod + 1
593 IF (((e22(i)/denom)<=one).AND.((e22(i)/denom)>=-half)) THEN
594 dfmax(i,1+nmod) = max(e11(i)/e1c(i),dfmax(i,1+nmod))
595 dfmax(i,1+nmod) = min(dfmax(i,1+nmod),one)
596 IF (e11(i) >= e1c(i)) THEN
597 ncrit(i) = ncrit(i) + 1
598 imxe1c(i) = 1
599 ENDIF
600 ENDIF
601 ENDIF
602 ! -> Forming limit diagram
603 IF (btest(crit,16)) THEN
604 nmod = nmod + 1
605 IF (itab == 1) THEN
606 dfmax(i,1+nmod) = max(e11(i)/(e1fld(i)*facl(i)),dfmax(i,1+nmod))
607 dfmax(i,1+nmod) = min(dfmax(i,1+nmod),one)
608 IF (e11(i) >= e1fld(i)*facl(i)) THEN
609 ncrit(i) = ncrit(i) + 1
610 ifld(i) = 1
611 ENDIF
612 ELSE
613 dfmax(i,1+nmod) = max(e11(i)/(e1fld(i)),dfmax(i,1+nmod))
614 dfmax(i,1+nmod) = min(dfmax(i,1+nmod),one)
615 IF (e11(i) >= e1fld(i)) THEN
616 ncrit(i) = ncrit(i) + 1
617 ifld(i) = 1
618 ENDIF
619 ENDIF
620 ENDIF
621 ! -> maximum shell thinning
622 IF (btest(crit,17)) THEN
623 nmod = nmod + 1
624 dfmax(i,1+nmod) = max((epszz(i))/(-abs(thin)*facl(i)),dfmax(i,1+nmod))
625 dfmax(i,1+nmod) = min(dfmax(i,1+nmod),one)
626 IF (epszz(i) <= -abs(thin)*facl(i)) THEN
627 ncrit(i) = ncrit(i) + 1
628 ithin(i) = 1
629 ENDIF
630 ENDIF
631 ! -> maximum temperature
632 IF (btest(crit,18)) THEN
633 nmod = nmod + 1
634 dfmax(i,1+nmod) = max(temp(i)/maxtemp,dfmax(i,1+nmod))
635 dfmax(i,1+nmod) = min(dfmax(i,1+nmod),one)
636 IF (temp(i) >= maxtemp) THEN
637 ncrit(i) = ncrit(i) + 1
638 imaxtemp(i) = 1
639 ENDIF
640 ENDIF
641c
642 ! -> Checking failure
643 DO j = 1,nmod
644 dfmax(i,1) = max(dfmax(i,1),dfmax(i,1+j))
645 ENDDO
646 dfmax(i,1) = min(dfmax(i,1),one)
647 IF (ncrit(i) >= ncs) THEN
648 uvar(i,5) = uvar(i,5) - one/nstep
649 signxx(i) = signxx(i)*uvar(i,5)
650 signyy(i) = signyy(i)*uvar(i,5)
651 signzz(i) = signzz(i)*uvar(i,5)
652 signxy(i) = signxy(i)*uvar(i,5)
653 signyz(i) = signyz(i)*uvar(i,5)
654 signzx(i) = signzx(i)*uvar(i,5)
655 dfmax(i,1) = one
656 nindx = nindx + 1
657 indx(nindx) = i
658 ENDIF
659 ENDIF
660c
661 ENDDO
662c
663 !====================================================================
664 ! - LOOKING FOR ELEMENT DELETION
665 !====================================================================
666 ! Case of under-integrated solids
667 IF (npg == 1) THEN
668 DO i = 1,nel
669 IF ((uvar(i,5) == zero).AND.(off(i) /= zero)) THEN
670 nindx2 = nindx2 + 1
671 indx2(nindx2) = i
672 dfmax(i,1) = one
673 off(i) = zero
674 tdele(i) = time
675 ENDIF
676 ENDDO
677 ! Case of fully integrated solids
678 ELSE
679 IF ((ipg == npg).AND.(ilay1 == elbuf_tab%NLAY)) THEN
680 DO i = 1,nel
681 IF (off(i) == one) THEN
682 ! Initialization of volumes computation
683 vtot = zero
684 vdam = zero
685 ! Loop over integration points
686 DO ilay = 1, elbuf_tab%NLAY
687 DO ir = 1, elbuf_tab%NPTR
688 DO is = 1, elbuf_tab%NPTS
689 DO it = 1, elbuf_tab%NPTT
690 fbuf => elbuf_tab%BUFLY(ilay)%FAIL(ir,is,it)
691 ! Compute total volume of the element
692 vtot = vtot + fbuf%FLOC(irupt)%VAR(2*nel+i)
693 ! Compute damaged volume of the element
694 vdam = vdam + fbuf%FLOC(irupt)%VAR(3*nel+i)
695 ENDDO
696 ENDDO
697 ENDDO
698 ENDDO
699 ! Checking volume fraction criterion
700 IF ((vdam/vtot) >= volfrac) THEN
701 dfmax(i,1) = one
702 nindx3 = nindx3 + 1
703 indx3(nindx3) = i
704 vfail(i) = vdam/vtot
705 tdele(i) = time
706 off(i) = zero
707 ENDIF
708 ENDIF
709 ENDDO
710 ENDIF
711 ENDIF
712c------------------------
713c------------------------
714 IF (nindx > 0) THEN
715 DO j=1,nindx
716 i = indx(j)
717#include "lockon.inc"
718 IF (ncrit(i) == 1) THEN
719 WRITE(iout, 1000) ngl(i),ipg,time,ncrit(i)
720 WRITE(istdo,1000) ngl(i),ipg,time,ncrit(i)
721 ELSE
722 WRITE(iout, 1001) ngl(i),ipg,time,ncrit(i)
723 WRITE(istdo,1001) ngl(i),ipg,time,ncrit(i)
724 ENDIF
725 IF (ipmax(i) == 1) THEN
726 WRITE(iout, 1002) p(i),maxpres*facl(i)
727 WRITE(istdo,1002) p(i),maxpres*facl(i)
728 ENDIF
729 IF (ipmin(i) == 1) THEN
730 WRITE(iout, 1003) p(i),minpres*facl(i)
731 WRITE(istdo,1003) p(i),minpres*facl(i)
732 ENDIF
733 IF (is1max(i) == 1) THEN
734 WRITE(iout, 1004) s11(i),abs(sigp1)*facl(i)
735 WRITE(istdo,1004) s11(i),abs(sigp1)*facl(i)
736 ENDIF
737 IF (itmax(i) == 1) THEN
738 WRITE(iout, 1005) time,tmax
739 WRITE(istdo,1005) time,tmax
740 ENDIF
741 IF (imindt(i) == 1) THEN
742 WRITE(iout, 1006) dt(i),dtmin
743 WRITE(istdo,1006) dt(i),dtmin
744 ENDIF
745 IF (isigmax(i) == 1) THEN
746 WRITE(iout, 1007) svm(i),sigmax(i)*facl(i)
747 WRITE(istdo,1007) svm(i),sigmax(i)*facl(i)
748 ENDIF
749 IF (isigth(i) == 1) THEN
750 WRITE(iout, 1008) uvar(i,2),kf*facl(i)
751 WRITE(istdo,1008) uvar(i,2),kf*facl(i)
752 ENDIF
753 IF (iepsmax(i) == 1) THEN
754 WRITE(iout, 1009) e11(i),epsmax(i)*facl(i)
755 WRITE(istdo,1009) e11(i),epsmax(i)*facl(i)
756 ENDIF
757 IF (ieffeps(i) == 1) THEN
758 WRITE(iout, 1010) eff_strain(i),effeps*facl(i)
759 WRITE(istdo,1010) eff_strain(i),effeps*facl(i)
760 ENDIF
761 IF (ivoleps(i) == 1) THEN
762 IF (voleps >= zero) THEN
763 WRITE(iout, 1011) vol_strain(i),voleps*facl(i)
764 WRITE(istdo,1011) vol_strain(i),voleps*facl(i)
765 ELSE
766 WRITE(iout, 1012) vol_strain(i),voleps*facl(i)
767 WRITE(istdo,1012) vol_strain(i),voleps*facl(i)
768 ENDIF
769 ENDIF
770 IF (imineps(i) == 1) THEN
771 WRITE(iout, 1013) e33(i),mineps*facl(i)
772 WRITE(istdo,1013) e33(i),mineps*facl(i)
773 ENDIF
774 IF (ishear(i) == 1) THEN
775 WRITE(iout, 1014) (e11(i) - e33(i))/two,epssh*facl(i)
776 WRITE(istdo,1014) (e11(i) - e33(i))/two,epssh*facl(i)
777 ENDIF
778 IF (imix12(i) == 1) THEN
779 WRITE(iout, 1015) (e11(i) - e22(i))/two,sh12(i)
780 WRITE(istdo,1015) (e11(i) - e22(i))/two,sh12(i)
781 ENDIF
782 IF (imix13(i) == 1) THEN
783 WRITE(iout, 1016) (e11(i) - e33(i))/two,sh13(i)
784 WRITE(istdo,1016) (e11(i) - e33(i))/two,sh13(i)
785 ENDIF
786 IF (imxe1c(i) == 1) THEN
787 WRITE(iout, 1017) e11(i),e1c(i)
788 WRITE(istdo,1017) e11(i),e1c(i)
789 ENDIF
790 IF (ifld(i) == 1) THEN
791 IF (itab == 1) THEN
792 WRITE(iout, 1018) e11(i),e1fld(i)*facl(i)
793 WRITE(istdo,1018) e11(i),e1fld(i)*facl(i)
794 ELSE
795 WRITE(iout, 1018) e11(i),e1fld(i)
796 WRITE(istdo,1018) e11(i),e1fld(i)
797 ENDIF
798 ENDIF
799 IF (ithin(i) == 1) THEN
800 WRITE(iout, 1019) epszz(i),-abs(thin)*facl(i)
801 WRITE(istdo,1019) epszz(i),-abs(thin)*facl(i)
802 ENDIF
803 IF (imaxtemp(i) == 1) THEN
804 WRITE(iout, 1020) temp(i),maxtemp
805 WRITE(istdo,1020) temp(i),maxtemp
806 ENDIF
807#include "lockoff.inc"
808 END DO
809 END IF
810 IF (nindx2 > 0) THEN
811 DO j=1,nindx2
812 i = indx2(j)
813#include "lockon.inc"
814 WRITE(iout, 1021) ngl(i),time
815 WRITE(istdo,1021) ngl(i),time
816#include "lockoff.inc"
817 END DO
818 END IF
819 IF (nindx3 > 0) THEN
820 DO j=1,nindx3
821 i = indx3(j)
822#include "lockon.inc"
823 WRITE(iout, 1021) ngl(i),time
824 WRITE(istdo,1021) ngl(i),time
825 WRITE(iout, 3000) vfail(i),volfrac
826 WRITE(istdo,3000) vfail(i),volfrac
827#include "lockoff.inc"
828 END DO
829 END IF
830c------------------------
831 1000 FORMAT(1x,'>> FOR SOLID ELEMENT NUMBER (GENE1) el#',i10,', GAUSS POINT # ',i5,
832 . ', FAILURE START AT TIME: ',1pe12.4,', ',i3,' CRITERION HAS BEEN REACHED:')
833 1001 FORMAT(1x,'>> FOR SOLID ELEMENT NUMBER (GENE1) el#',i10,', GAUSS POINT # ',i5,
834 . ', FAILURE START AT TIME: ',1pe12.4,', ',i3,' CRITERIA HAVE BEEN REACHED:')
835 1002 FORMAT(1x,'HYDROSTATIC PRESSURE VALUE: ',1pe12.4,' > CRITICAL VALUE: ',1pe12.4)
836 1003 FORMAT(1x,'HYDROSTATIC PRESSURE VALUE: ',1pe12.4,' < CRITICAL VALUE: ',1pe12.4)
837 1004 FORMAT(1x,'1ST PRINCIPAL STRESS VALUE: ',1pe12.4,' > CRITICAL VALUE: ',1pe12.4)
838 1005 FORMAT(1x,'TIME VALUE: ',1pe12.4,' > CRITICAL VALUE: ',1pe12.4)
839 1006 FORMAT(1x,'ELEMENT TIMESTEP VALUE: ',1pe12.4,' < CRITICAL VALUE: ',1pe12.4)
840 1007 FORMAT(1x,'EQUIVALENT STRESS VALUE: ',1pe12.4,' > CRITICAL VALUE: ',1pe12.4)
841 1008 FORMAT(1x,'T-BUTCHER INTG. VALUE: ',1pe12.4,' > CRITICAL VALUE: ',1pe12.4)
842 1009 FORMAT(1x,'1ST PRINCIPAL STRAIN VALUE: ',1pe12.4,' > CRITICAL VALUE: ',1pe12.4)
843 1010 FORMAT(1x,'EFFECTIVE STRAIN VALUE: ',1pe12.4,' > CRITICAL VALUE: ',1pe12.4)
844 1011 FORMAT(1x,'VOLUMETRIC STRAIN VALUE: ',1pe12.4,' > CRITICAL VALUE: ',1pe12.4)
845 1012 FORMAT(1x,'VOLUMETRIC STRAIN VALUE: ',1pe12.4,' < CRITICAL VALUE: ',1pe12.4)
846 1013 FORMAT(1x,'3RD PRINCIPAL STRAIN VALUE: ',1pe12.4,' < CRITICAL VALUE: ',1pe12.4)
847 1014 FORMAT(1x,'MAX. SHEAR STRAIN VALUE: ',1pe12.4,' > CRITICAL VALUE: ',1pe12.4)
848 1015 FORMAT(1x,'IN-PLANE SH.STRAIN 12 VALUE: ',1pe12.4,' > CRITICAL VALUE: ',1pe12.4)
849 1016 FORMAT(1x,'TRANSV. SH.STRAIN 13 VALUE: ',1pe12.4,' > CRITICAL VALUE: ',1pe12.4)
850 1017 FORMAT(1x,'IN-PLANE PRINC.STRAIN VALUE: ',1pe12.4,' > CRITICAL VALUE: ',1pe12.4)
851 1018 FORMAT(1x,'1ST PRINCIPAL STRESS VALUE: ',1pe12.4,' > FORMING LIMIT VALUE : ',1pe12.4)
852 1019 FORMAT(1x,'THINNING VALUE: ',1pe12.4,' < CRITICAL VALUE: ',1pe12.4)
853 1020 FORMAT(1x,'TEMPERATURE VALUE: ',1pe12.4,' > CRITICAL VALUE: ',1pe12.4)
854 1021 FORMAT(1x,'-- RUPTURE OF SOLID ELEMENT : ',i10,' AT TIME :',1pe12.4)
855 3000 FORMAT(1x,'DAMAGED VOLUME FRACTION : ',1pe12.4,' > CRITICAL VALUE: ',1pe12.4)
856 END
fail_gene1_s
subroutine fail_gene1_s(nel, nuparam, nuvar, nfunc, ifunc, loff, npf, tf, time, timestep, uparam, ipg, ngl, dt, epsp, uvar, off, npg, epsxx, epsyy, epszz, epsxy, epsyz, epszx, signxx, signyy, signzz, signxy, signyz, signzx, temp, voln, dfmax, tdele, aldt, table, irupt, elbuf_tab, ilay1, ntablf, itablf, lf_dammx, niparam, iparam)
Definition fail_gene1_s.F:47
interface_table_mod::table2d_vinterp_log
Definition table_mod.F:261
interface_table_mod::table_vinterp
Definition table_mod.F:251
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
interface_table_mod
Definition table_mod.F:238
table_mod
Definition table_mod.F:112
table2d_vinterp_log
subroutine table2d_vinterp_log(table, ismooth, dimx, nel, ipos, xx, yy, dydx1, dydx2)
Definition table2d_vinterp_log.F:52