OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fail_tab2_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_tab2_s ../engine/source/materials/fail/tabulated/fail_tab2_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!|| table_vinterp ../engine/source/tools/curve/table_tools.F
32!|| vinter2 ../engine/source/tools/curve/vinter.F
33!||--- uses -----------------------------------------------------
34!|| interface_table_mod ../engine/share/modules/table_mod.F
35!|| table_mod ../engine/share/modules/table_mod.F
36!||====================================================================
37 SUBROUTINE fail_tab2_s (
38 1 NEL ,NUPARAM ,NUVAR ,NFUNC ,IFUNC ,
39 2 NPF ,TABLE ,TF ,TIME ,UPARAM ,
40 3 NGL ,ALDT ,DPLA ,EPSP ,UVAR ,
41 4 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
42 5 TEMP ,OFF ,DFMAX ,TDELE ,DMG_SCALE,
43 6 UELR ,IPG ,NPG ,LOFF ,NTABLF ,ITABLF ,
44 7 NOFF ,VOLN )
45C---------+---------+---+---+--------------------------------------------
46 USE table_mod
47 USE interface_table_mod
48C-----------------------------------------------
49C I m p l i c i t T y p e s
50C-----------------------------------------------
51#include "implicit_f.inc"
52C-----------------------------------------------
53C C o m m o n B l o c k s
54C-----------------------------------------------
55#include "scr17_c.inc"
56#include "units_c.inc"
57#include "comlock.inc"
58#include "com04_c.inc"
59#include "tabsiz_c.inc"
60C-----------------------------------------------
61C I N P U T A r g u m e n t s
62C-----------------------------------------------
63 INTEGER, INTENT(IN) :: NEL,NUPARAM,NUVAR,NGL(NEL),NPG,IPG,NTABLF
64 INTEGER, DIMENSION(NTABLF) ,INTENT(IN) :: ITABLF
65 INTEGER, INTENT(INOUT) :: NOFF(NEL)
66 my_real, INTENT(IN) ::
67 . TIME,UPARAM(NUPARAM),ALDT(NEL),
68 . DPLA(NEL),EPSP(NEL),TEMP(NEL),
69 . SIGNXX(NEL),SIGNYY(NEL),SIGNZZ(NEL),
70 . signxy(nel),signyz(nel),signzx(nel),
71 . voln(nel)
72 my_real, INTENT(INOUT) ::
73 . uvar(nel,nuvar),dfmax(nel),tdele(nel),
74 . dmg_scale(nel),off(nel),uelr(nel),
75 . loff(nel)
76 TYPE (TTABLE), INTENT(IN), DIMENSION(NTABLE) :: TABLE
77C!-----------------------------------------------
78C! VARIABLES FOR FUNCTION INTERPOLATION
79C!-----------------------------------------------
80 INTEGER, INTENT(IN) :: NPF(SNPC),NFUNC,IFUNC(NFUNC)
81 my_real, INTENT(IN) :: TF(STF)
82 my_real
83 . finter
84 EXTERNAL finter
85C-----------------------------------------------
86C L o c a l V a r i a b l e s
87C-----------------------------------------------
88 INTEGER I,J,INDX(NEL),NINDX,ITAB_EPSF,FAILIP,
89 . itab_inst,itab_size,ireg,ipos(nel,3),
90 . ndim,ipos2(nel),iad(nel),ilen(nel),
91 . log_scale1,log_scale2
92 my_real
93 . fcrit ,dn,dcrit,ecrit,exp_ref,expo,el_ref,
94 . sr_ref1,fscale_el,shrf,biaxf ,sr_ref2,
95 . fscale_sr,cjc,fscale_dlim,temp_ref, fscale_temp, volfrac
96 my_real
97 . lambda,fac,df,inst(nel) ,dc(nel) ,l0(nel) ,
98 . triax(nel) ,xi(nel) ,epsf(nel) ,epsl(nel) ,
99 . depsf(nel) ,depsl(nel) ,xvec(nel,3) ,dpl_def ,
100 . cos3theta ,det,p,svm ,sxx,syy,szz ,sizefac(nel),
101 . ratefac(nel),dsize(nel) ,softexp(nel),dlim(nel) ,
102 . tempfac(nel),tempfac2(nel),dft(nel) ,var(nel)
103C!--------------------------------------------------------------
104 !=======================================================================
105 ! - INITIALISATION OF COMPUTATION ON TIME STEP
106 !=======================================================================
107 ! Recovering failure criterion parameters
108 fcrit = uparam(1)
109 failip = min(nint(uparam(2)),npg)
110 dn = uparam(4)
111 dcrit = uparam(5)
112 ecrit = uparam(6)
113 exp_ref = uparam(7)
114 expo = uparam(8)
115 ireg = nint(uparam(9))
116 el_ref = uparam(10)
117 sr_ref1 = uparam(11)
118 fscale_el = uparam(12)
119 shrf = uparam(13)
120 biaxf = uparam(14)
121 sr_ref2 = uparam(15)
122 fscale_sr = uparam(16)
123 cjc = uparam(17)
124 fscale_dlim = uparam(18)
125 temp_ref = uparam(19)
126 fscale_temp = uparam(20)
127 log_scale1 = nint(uparam(21))
128 log_scale2 = nint(uparam(22))
129 volfrac = uparam(23)
130c
131 itab_epsf = itablf(1)
132 itab_inst = itablf(2)
133 itab_size = itablf(3)
134c
135 ! Checking element failure and recovering user variable
136 DO i=1,nel
137 ! If necking control is activated
138 IF ((itab_inst > 0).OR.(ecrit > zero)) THEN
139 ! Instability damage
140 inst(i) = uvar(i,1)
141 ! Necking critical damage
142 IF (uvar(i,2) == zero) uvar(i,2) = one
143 dc(i) = uvar(i,2)
144 ENDIF
145 END DO
146c
147 !====================================================================
148 ! - LOOP OVER THE ELEMENT TO COMPUTE THE STRESS STATE QUANTITIES
149 !====================================================================
150 DO i=1,nel
151c
152 ! Computation of hydrostatic stress, Von Mises stress, and stress triaxiality
153 p = third*(signxx(i) + signyy(i) + signzz(i))
154 sxx = signxx(i) - p
155 syy = signyy(i) - p
156 szz = signzz(i) - p
157 svm = half*(sxx**2 + syy**2 + szz**2)
158 . + signxy(i)**2 + signzx(i)**2 + signyz(i)**2
159 svm = sqrt(max(three*svm,zero))
160 triax(i) = p/max(em20,svm)
161 IF (triax(i) < -one) triax(i) = -one
162 IF (triax(i) > one) triax(i) = one
163c
164 ! Computation of Lode parameter
165 det = sxx*syy*szz + two*signxy(i)*signzx(i)*signyz(i)
166 . - sxx*signyz(i)**2-szz*signxy(i)**2 - syy*signzx(i)**2
167 cos3theta = half*twenty7*det/max(em20,svm**3)
168 IF (cos3theta < -one) cos3theta = -one
169 IF (cos3theta > one) cos3theta = one
170 xi(i) = one - two*acos(cos3theta)/pi
171c
172 ENDDO
173c
174 !====================================================================
175 ! - COMPUTE FACTORS FOR ELEMENT SIZE, STRAIN RATE AND TEMPERATURE
176 !====================================================================
177 ! At initial time, save the element size
178 IF (uvar(1,3) == zero) uvar(1:nel,3) = aldt(1:nel)
179 l0(1:nel) = uvar(1:nel,3)
180c
181 ! Compute the softening exponent
182 IF (ifunc(1) > 0) THEN
183 DO i=1,nel
184 lambda = l0(i)/exp_ref
185 softexp(i) = finter(ifunc(1),lambda,npf,tf,df)
186 softexp(i) = expo*softexp(i)
187 ENDDO
188 ELSE
189 softexp(1:nel) = expo
190 ENDIF
191c
192 ! Compute the temperature dependency factor
193 IF (ifunc(4) > 0) THEN
194 var(1:nel) = temp(1:nel)/temp_ref
195 ipos2(1:nel) = 1
196 iad(1:nel) = npf(ifunc(4)) / 2 + 1
197 ilen(1:nel) = npf(ifunc(4)+1) / 2 - iad(1:nel) - ipos2(1:nel)
198 CALL vinter2(tf,iad,ipos2,ilen,nel,var,dft,tempfac)
199 tempfac(1:nel) = fscale_temp*tempfac(1:nel)
200 tempfac2(1:nel) = tempfac(1:nel)
201 ELSE
202 tempfac(1:nel) = one
203 tempfac2(1:nel) = one
204 ENDIF
205c
206 ! Compute the element size regularization factor
207 IF (itab_size > 0) THEN
208 ! Element size scaling dependency
209 ndim = table(itab_size)%NDIM
210 IF (ireg == 1) THEN
211 SELECT CASE (ndim)
212 ! Scale factor vs element size
213 CASE(1)
214 xvec(1:nel,1) = l0(1:nel)/el_ref
215 xvec(1:nel,2:3) = zero
216 ipos(1:nel,1:3) = 1
217 ! Scale factor vs element size vs strain rate
218 CASE(2)
219 xvec(1:nel,1) = l0(1:nel)/el_ref
220 IF (log_scale1 > 0) THEN
221 DO i = 1,nel
222 xvec(i,2) = log(max(epsp(i),em20)/sr_ref1)
223 ENDDO
224 ELSE
225 xvec(1:nel,2) = epsp(1:nel)/sr_ref1
226 ENDIF
227 xvec(1:nel,3) = zero
228 ipos(1:nel,1:3) = 1
229 END SELECT
230 ELSEIF (ireg == 2) THEN
231 SELECT CASE (ndim)
232 ! Scale factor vs element size
233 CASE(1)
234 xvec(1:nel,1) = l0(1:nel)/el_ref
235 xvec(1:nel,2:3) = zero
236 ipos(1:nel,1:3) = 1
237 ! Scale factor vs element size vs triaxiality
238 CASE(2)
239 xvec(1:nel,1) = l0(1:nel)/el_ref
240 xvec(1:nel,2) = triax(1:nel)
241 xvec(1:nel,3) = zero
242 ipos(1:nel,1:3) = 1
243 ! Scale factor vs element size vs triaxiality vs Lode parameter
244 CASE(3)
245 xvec(1:nel,1) = l0(1:nel)/el_ref
246 xvec(1:nel,2) = triax(1:nel)
247 xvec(1:nel,3) = xi(1:nel)
248 ipos(1:nel,1:3) = 1
249 END SELECT
250 ENDIF
251 CALL table_vinterp(table(itab_size),nel,nel,ipos,xvec,sizefac,dsize)
252 sizefac(1:nel) = sizefac(1:nel)*fscale_el
253 IF (ireg == 1) THEN
254 DO i = 1,nel
255 IF (triax(i) < shrf) THEN
256 sizefac(i) = one
257 ELSEIF (triax(i) > biaxf) THEN
258 sizefac(i) = one
259 ENDIF
260 ENDDO
261 ENDIF
262 ELSE
263 sizefac(1:nel) = one
264 ENDIF
265c
266 ! Compute the strain rate dependency factor
267 IF (ifunc(2) > 0) THEN
268 IF (log_scale2 > 0) THEN
269 DO i = 1,nel
270 var(i) = log(max(epsp(i),em20)/sr_ref2)
271 ENDDO
272 ELSE
273 var(1:nel) = epsp(1:nel)/sr_ref2
274 ENDIF
275 ipos2(1:nel) = 1
276 iad(1:nel) = npf(ifunc(2)) / 2 + 1
277 ilen(1:nel) = npf(ifunc(2)+1) / 2 - iad(1:nel) - ipos2(1:nel)
278 CALL vinter2(tf,iad,ipos2,ilen,nel,var,dft,ratefac)
279 ratefac(1:nel) = fscale_sr*ratefac(1:nel)
280 ELSEIF (cjc > zero) THEN
281 DO i=1,nel
282 IF (epsp(i) > sr_ref2) THEN
283 ratefac(i) = one + cjc*log(epsp(i)/sr_ref2)
284 ELSE
285 ratefac(i) = one
286 ENDIF
287 ENDDO
288 ELSE
289 ratefac(1:nel) = one
290 ENDIF
291c
292 ! Compute the damage limit value
293 IF (ifunc(3) > 0) THEN
294 DO i = 1,nel
295 lambda = triax(i)
296 dlim(i) = finter(ifunc(3),lambda,npf,tf,df)
297 dlim(i) = fscale_dlim*dlim(i)
298 dlim(i) = min(dlim(i),one)
299 dlim(i) = max(dlim(i),zero)
300 ENDDO
301 ELSE
302 dlim(1:nel) = one
303 ENDIF
304c
305 !====================================================================
306 ! - COMPUTATION OF PLASTIC STRAIN AT FAILURE
307 !====================================================================
308 IF (itab_epsf > 0) THEN
309 ! Failure plastic strain map dependency
310 ndim = table(itab_epsf)%NDIM
311 SELECT CASE (ndim)
312 ! Failure plastic strain vs triaxiality
313 CASE (1)
314 xvec(1:nel,1) = triax(1:nel)
315 xvec(1:nel,2:3) = zero
316 ipos(1:nel,1:3) = 1
317 ! Failure plastic strain vs triaxiality vs Lode parameter
318 CASE (2)
319 xvec(1:nel,1) = triax(1:nel)
320 xvec(1:nel,2) = xi(1:nel)
321 xvec(1:nel,3) = zero
322 ipos(1:nel,1:3) = 1
323 ! Failure plastic strain vs triaxiality vs Lode parameter vs temperature
324 CASE (3)
325 xvec(1:nel,1) = triax(1:nel)
326 xvec(1:nel,2) = xi(1:nel)
327 xvec(1:nel,3) = temp(1:nel)/temp_ref
328 ipos(1:nel,1:3) = 1
329 tempfac(1:nel) = one
330 END SELECT
331 CALL table_vinterp(table(itab_epsf),nel,nel,ipos,xvec,epsf,depsf)
332 epsf(1:nel) = epsf(1:nel)*fcrit
333 ELSE
334 epsf(1:nel) = fcrit
335 ENDIF
336c
337 !====================================================================
338 ! - COMPUTATION OF PLASTIC STRAIN AT NECKING
339 !====================================================================
340 IF (itab_inst > 0) THEN
341 ! Instability plastic strain map dependency
342 ndim = table(itab_inst)%NDIM
343 SELECT CASE (ndim)
344 ! Instability plastic strain vs triaxiality
345 CASE(1)
346 xvec(1:nel,1) = triax(1:nel)
347 xvec(1:nel,2:3) = zero
348 ipos(1:nel,1:3) = 1
349 ! Instability plastic strain vs triaxiality vs Lode
350 CASE(2)
351 xvec(1:nel,1) = triax(1:nel)
352 xvec(1:nel,2) = xi(1:nel)
353 xvec(1:nel,3) = zero
354 ipos(1:nel,1:3) = 1
355 ! Instability plastic strain vs triaxiality vs Lode vs temperature
356 CASE(3)
357 xvec(1:nel,1) = triax(1:nel)
358 xvec(1:nel,2) = xi(1:nel)
359 xvec(1:nel,3) = temp(1:nel)/temp_ref
360 ipos(1:nel,1:3) = 1
361 tempfac2(1:nel) = one
362 END SELECT
363 CALL table_vinterp(table(itab_inst),nel,nel,ipos,xvec,epsl,depsl)
364 epsl(1:nel) = epsl(1:nel)*ecrit
365 ELSEIF (ecrit > zero) THEN
366 epsl(1:nel) = ecrit
367 ENDIF
368c
369 !====================================================================
370 ! - COMPUTATION OF THE DAMAGE VARIABLE EVOLUTION
371 !====================================================================
372 ! Initialization of element failure index
373 nindx = 0
374 indx(1:nel) = 0
375c
376 ! Loop over the elements
377 DO i=1,nel
378c
379 ! If the element is not broken
380 IF (loff(i) == one .AND. off(i) == one .AND. dpla(i) > zero) THEN
381c
382 ! Needs to initialize damage at a very small value the first time
383 IF (dfmax(i) == zero) dfmax(i) = em20
384 IF (inst(i) == zero) inst(i) = em20
385c
386 ! Compute failure strain damage variable
387 dpl_def = dpla(i)/max(epsf(i)*ratefac(i)*sizefac(i)*tempfac(i),em20)
388 dfmax(i) = dfmax(i) + dpl_def*dn*(dfmax(i)**(one-(one/dn)))
389 dfmax(i) = min(dfmax(i),dlim(i))
390 IF (dfmax(i) >= one) THEN
391 nindx = nindx + 1
392 indx(nindx) = i
393 noff(i) = noff(i) + 1
394 loff(i) = zero
395 uelr(i) = uelr(i) + voln(i)
396 IF ((npg > 1).AND.(volfrac > zero)) THEN
397 IF (uelr(i)/(npg*voln(i)) >= volfrac) THEN
398 off(i) = zero
399 tdele(i) = time
400 ENDIF
401 ELSE
402 IF (noff(i) >= failip) THEN
403 off(i) = zero
404 tdele(i) = time
405 ENDIF
406 ENDIF
407 ENDIF
408c
409 ! Compute the control necking instability damage
410 IF ((itab_inst > 0).OR.(ecrit > zero)) THEN
411 dpl_def = dpla(i)/max(epsl(i)*ratefac(i)*sizefac(i)*tempfac2(i),em20)
412 inst(i) = inst(i) + dpl_def*dn*(inst(i)**(one-(one/dn)))
413 inst(i) = min(inst(i),one)
414 IF ((inst(i) >= one).AND.(dc(i) == one)) THEN
415 dc(i) = dfmax(i)
416 ENDIF
417 ENDIF
418c
419 ENDIF
420 ENDDO
421c
422 !====================================================================
423 ! - UPDATE UVAR AND THE STRESS TENSOR
424 !====================================================================
425 DO i = 1,nel
426 IF ((itab_inst > 0).OR.(ecrit > zero)) THEN
427 uvar(i,1) = inst(i)
428 uvar(i,2) = dc(i)
429 IF (dfmax(i) >= dc(i)) THEN
430 IF (dc(i) < one) THEN
431 dmg_scale(i) = one - ((dfmax(i)-dc(i))/max(one-dc(i),em20))**softexp(i)
432 ELSE
433 dmg_scale(i) = zero
434 ENDIF
435 ELSE
436 dmg_scale(i) = one
437 ENDIF
438 ELSE
439 IF (dfmax(i) >= dcrit) THEN
440 IF (dcrit < one) THEN
441 dmg_scale(i) = one - ((dfmax(i)-dcrit)/max(one-dcrit,em20))**softexp(i)
442 ELSE
443 dmg_scale(i) = zero
444 ENDIF
445 ELSE
446 dmg_scale(i) = one
447 ENDIF
448 ENDIF
449 ENDDO
450c
451 !====================================================================
452 ! - PRINTOUT DATA ABOUT FAILED ELEMENTS
453 !====================================================================
454 IF (nindx > 0) THEN
455 DO j=1,nindx
456 i = indx(j)
457#include "lockon.inc"
458 WRITE(iout, 1000) ngl(i),ipg,time
459 WRITE(istdo,1000) ngl(i),ipg,time
460 IF (off(i) == zero) THEN
461 WRITE(iout, 2000) ngl(i),time
462 WRITE(istdo,2000) ngl(i),time
463 ENDIF
464#include "lockoff.inc"
465 END DO
466 END IF
467c-----------------------------------------------------------------------
468 1000 FORMAT(1x,'FOR SOLID ELEMENT NUMBER el#',i10,
469 . ' FAILURE (TAB2) AT GAUSS POINT ',i5,
470 . ' AT TIME :',1pe12.4)
471 2000 FORMAT(1x,'-- RUPTURE OF SOLID ELEMENT :',i10,
472 . ' AT TIME :',1pe12.4)
473c
474 END
ecrit
subroutine ecrit(timers, partsav, ms, v, in, r, dmas, weight, enintot, ekintot, a, ar, fxbipm, fxbrpm, monvol, xmom_sms, sensors, qfricint, ipari, weight_md, wfexth, iflag, ms_2d, multi_fvm, mas_nd, kend, h3d_data, dynain_data, usreint, output)
Definition ecrit.F:52
fail_tab2_s
subroutine fail_tab2_s(nel, nuparam, nuvar, nfunc, ifunc, npf, table, tf, time, uparam, ngl, aldt, dpla, epsp, uvar, signxx, signyy, signzz, signxy, signyz, signzx, temp, off, dfmax, tdele, dmg_scale, uelr, ipg, npg, loff, ntablf, itablf, noff, voln)
Definition fail_tab2_s.F:45
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
usermat_solid
subroutine usermat_solid(timers, lft, llt, nft, mtn, jcvt, pm, off, sig, eint, rho, qold, vol, strain, sigl, gama, uvar, bufmat, tf, npf, stifn, mat, ngl, nuvar, dt2t, neltst, ityptst, offg, geo, pid, epsd, el_temp, wxx, wyy, wzz, jsph, mumax, ssp, aire, voln, vd2, deltax, vis, d1, d2, d3, d4, d5, d6, pnew, psh, q, ssp_eq, dvol, sold1, sold2, sold3, sold4, sold5, sold6, rx, ry, rz, sx, sy, sz, tx, ty, tz, ipla, sigy, defp, ismstr, mfxx, mfxy, mfxz, mfyx, mfyy, mfyz, mfzx, mfzy, mfzz, ipm, isorth, fbuf, nfail, npg, sigdd, dxy, dyx, dyz, dzy, dzx, dxz, fr_wav, isrot, v, varnl, w, ix, x, jthe, et, mssa, dmels, iptr, ipts, iptt, table, fvd2, fdeltax, fssp, fqvis, iparg1, igeo, sigv, al_imp, signor, istrain, ng, elbuf_tab, vbuf, ilay, vk, iparg, bufvois, vdx, vdy, vdz, ihet, conde, itask, iexpan, vol_avg, amu, epsth3, epsth, svisc, nel, etotsh, iselect, tstar, muold, amu2, dpdm, rhoref, rhosp, nloc_dmg, ity, jtur, mat_elem, idel7nok, svis, dt, glob_therm, damp_buf, idamp_freq_range)
Definition usermat_solid.F:131
vinter2
subroutine vinter2(tf, iad, ipos, ilen, nel0, x, dydx, y)
Definition vinter.F:143