OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_gene1_s.F File Reference
#include "implicit_f.inc"
#include "scr17_c.inc"
#include "units_c.inc"
#include "comlock.inc"
#include "com04_c.inc"
#include "com01_c.inc"
#include "tabsiz_c.inc"
#include "lockon.inc"
#include "lockoff.inc"

Go to the source code of this file.

Functions/Subroutines

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)

Function/Subroutine Documentation

◆ fail_gene1_s()

subroutine fail_gene1_s ( integer, intent(in) nel,
integer, intent(in) nuparam,
integer, intent(in) nuvar,
integer, intent(in) nfunc,
integer, dimension(nfunc), intent(in) ifunc,
intent(inout) loff,
integer, dimension(snpc), intent(in) npf,
dimension(stf), intent(in) tf,
intent(in) time,
intent(in) timestep,
intent(in) uparam,
integer, intent(in) ipg,
integer, dimension(nel), intent(in) ngl,
intent(in) dt,
intent(in) epsp,
intent(inout) uvar,
intent(inout) off,
integer, intent(in) npg,
intent(in) epsxx,
intent(in) epsyy,
intent(in) epszz,
intent(in) epsxy,
intent(in) epsyz,
intent(in) epszx,
intent(inout) signxx,
intent(inout) signyy,
intent(inout) signzz,
intent(inout) signxy,
intent(inout) signyz,
intent(inout) signzx,
intent(in) temp,
intent(in) voln,
intent(inout) dfmax,
intent(inout) tdele,
intent(in) aldt,
type (ttable), dimension(ntable) table,
integer, intent(in) irupt,
type (elbuf_struct_), target elbuf_tab,
integer, intent(in) ilay1,
integer, intent(in) ntablf,
integer, dimension(ntablf), intent(in) itablf,
integer, intent(in) lf_dammx,
integer, intent(in) niparam,
integer, dimension(niparam), intent(in) iparam )

Table dependency type (used only if tab_IDfld is a table).

Number of cycles for the stress reduction, Default = 10 (Integer)

Number of conditions to reach before the element is deleted, Default = 1 (Integer)

Interpolation type (in case of tabulated yield function)

Engineering / True input strain flag

Minimum pressure (positive in compression)

Maximum pressure (positive in compression)

Maximum principal stress

Failure time, Default = 1E+20 (Real)

Minimum time step

Reference strain rate value for fct_IDsm, Default = 1 (Real)

Von Mises stress

Reference strain rate value for fct_IDps, Default = 1 (Real)

Ordinate scale factor for fct_IDps or maximum principal strain if fct_IDps is not defined, Default = 1

Maximum effective strain

Maximum volumetric strain

Minimum principal strain

Damaged volume fraction to reach before the element is deleted (fully-integrated and higher order solid elements only), Default = 0.5 (Real)

Maximum temperature.

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)

Reference element size for fct_IDel, tab_IDfld (Itab=2), fct_IDg12, fct_IDg23, fct_IDg13 and fct_IDe1c, Default = 1.0 (Real)

Function identifier of the maximum equivalent stress versus strain rate

maximum principal strain VS strain-rate

in-plane shear strain VS element size

in-plane shear strain VS element size

major in plane-strain VS element size

element size regularization

Definition at line 38 of file fail_gene1_s.F.

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)
my_real
#define my_real
Definition cppsort.cpp:32
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