OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps120_tab_vm.F File Reference
#include "implicit_f.inc"
#include "com04_c.inc"
#include "units_c.inc"
#include "comlock.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps120_tab_vm (nel, ngl, nuparam, nuvar, time, timestep, uparam, uvar, off, pla, epsd, soundsp, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, depsxx, depsyy, depszz, depsxy, depsyz, depszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, numtabl, itable, table, nvartmp, vartmp, temp, inloc, dplanl, dmg, dmg_scale)

Function/Subroutine Documentation

◆ sigeps120_tab_vm()

subroutine sigeps120_tab_vm ( integer nel,
integer, dimension(nel), intent(in) ngl,
integer nuparam,
integer nuvar,
time,
timestep,
intent(in) uparam,
intent(inout) uvar,
intent(inout) off,
intent(inout) pla,
intent(out) epsd,
intent(out) soundsp,
intent(in) epspxx,
intent(in) epspyy,
intent(in) epspzz,
intent(in) epspxy,
intent(in) epspyz,
intent(in) epspzx,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depszz,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
intent(in) sigoxx,
intent(in) sigoyy,
intent(in) sigozz,
intent(in) sigoxy,
intent(in) sigoyz,
intent(in) sigozx,
intent(out) signxx,
intent(out) signyy,
intent(out) signzz,
intent(out) signxy,
intent(out) signyz,
intent(out) signzx,
integer numtabl,
integer, dimension(ntable), intent(in) itable,
type(ttable), dimension(ntable) table,
integer nvartmp,
integer, dimension(nel,nvartmp), intent(inout) vartmp,
intent(in) temp,
integer inloc,
intent(inout) dplanl,
intent(inout) dmg,
intent(inout) dmg_scale )

Definition at line 33 of file sigeps120_tab_vm.F.

42C-----------------------------------------------
43C M o d u l e s
44C-----------------------------------------------
45 USE table_mod
46 USE interface_table_mod
47C-----------------------------------------------
48C I m p l i c i t T y p e s
49C-----------------------------------------------
50#include "implicit_f.inc"
51C-----------------------------------------------
52C C o m m o n B l o c k s
53C-----------------------------------------------
54#include "com04_c.inc"
55#include "units_c.inc"
56#include "comlock.inc"
57C----------------------------------------------------------------
58C D u m m y A R G U M E N T S
59C----------------------------------------------------------------
60 INTEGER :: NEL,NUPARAM,NUVAR,NVARTMP,NUMTABL,INLOC
61 my_real :: time,timestep
62 INTEGER ,DIMENSION(NEL) ,INTENT(IN) :: NGL
63 INTEGER ,DIMENSION(NTABLE) ,INTENT(IN) :: ITABLE
64 my_real ,DIMENSION(NUPARAM),INTENT(IN) :: uparam
65 my_real ,DIMENSION(NEL) ,INTENT(IN) :: temp,
66 . epspxx,epspyy,epspzz,epspxy,epspyz,epspzx,
67 . depsxx,depsyy,depszz,depsxy,depsyz,depszx,
68 . sigoxx,sigoyy,sigozz,sigoxy,sigoyz,sigozx
69c
70 my_real ,DIMENSION(NEL) ,INTENT(OUT) :: epsd,soundsp
71 my_real ,DIMENSION(NEL) ,INTENT(OUT) ::
72 . signxx,signyy,signzz,signxy,signyz,signzx
73 INTEGER ,DIMENSION(NEL,NVARTMP) ,INTENT(INOUT) :: VARTMP
74 my_real ,DIMENSION(NEL,NUVAR) ,INTENT(INOUT) :: uvar
75 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: off,pla,dplanl,dmg,dmg_scale
76 TYPE(TTABLE), DIMENSION(NTABLE) :: TABLE
77C----------------------------------------------------------------
78C L O C A L V A R I B L E S
79C----------------------------------------------------------------
80 INTEGER :: I,II,ITER,NITER,NINDX,NINDF,ITRX,IDAM
81 my_real :: xscale,yscale
82 INTEGER :: IPOS1(NEL,1), IPOS2(NEL,2), IPOS3(NEL,3)
83 INTEGER :: NDIM_YLD,FUNC_YLD
84 my_real :: xvec1(nel,1), xvec2(nel,2), xvec3(nel,3)
85 INTEGER ,DIMENSION(NEL) :: INDX,INDF
86c material parameters
87 my_real :: young,nu,ssp,g,g2,bulk,yld0,c11,c12,
88 . a1f,a2f,as,d1c,d2c,d1f,d2f,d_trx,d_jc,dm,expn
89c Johnson & Cook rate-dependency
90 my_real cc,epsdref,epsdmax,frate,fcut,alpha,alphai,dtinv
91c Local variables
92 my_real facr,triax,i1p,ip,jp,tp,fact,facd,dfac,
93 . dlam,ddep,dfdlam,dpla_dlam,dphi_dlam,
94 . gamc,gamf,np_xx,np_yy,np_zz,np_xy,np_yz,np_zx,
95 . nf_xx,nf_yy,nf_zz,nf_xy,nf_yz,nf_zx,dpxx,dpyy,dpzz,dpxy,dpyz,dpzx
96 my_real ,DIMENSION(NEL) :: i1,j2,yld,phi,dam,jcrate,
97 . sxx,syy,szz,sxy,syz,szx,dpla,jcr_log,gamma,dgamm,hardp,fvm,
98 . yld_i,hardp_i,pla_nl,dpdt_nl
99c--------------------------------
100c Material parameters
101c--------------------------------
102c UPARAM(1) = Young modulus E
103c UPARAM(2) = Poisson's ratio nu
104c UPARAM(3) = shear modulus G
105c UPARAM(4) = bulk modulus K
106c UPARAM(5) = Yld0: instant yield stress
107c UPARAM(6) = Q : expontial term coefficient in the hardening law
108c UPARAM(7) = BETA : exponent of the exponential term
109c UPARAM(8) = P : multiplier of the linear term in the hardening law
110c UPARAM(9) = A1F : parameter of the yield function
111c UPARAM(10) = A2F : parameter of the yield function
112c UPARAM(11) = A1H : distortionnal yield hardening coefficiant
113c UPARAM(12) = A2H : distortionnal yield hardening coefficiant
114c UPARAM(13) = AS : parameter of the potential function
115c UPARAM(14) = D1C: first damage strain parameter in initial damage threshold
116c UPARAM(15) = D2C: second damage strain parameter in initial damage threshold
117c UPARAM(16) = D1F: first damage strain parameter in final damage threshold
118c UPARAM(17) = D2F: second damage strain parameter in final damage threshold
119c UPARAM(18) = D_TRX : triaxiality factor in damage formula
120c UPARAM(19) = D_JC: JC strain rate coefficient in damage formula
121c UPARAM(20) = EXPN exponent in damage evolution
122c UPARAM(21) = CC : Johnson-Cook strain rate-dependency coefficient
123c UPARAM(22) = EPSDREF quasi-static reference strain rate
124c UPARAM(23) = EPSDMAX maximal reference strain rate
125c UPARAM(24) = FCUT : cut frequency for strain rate filtering
126c UPARAM(25) = IFORM = 1: Yield formulation flag => Drucker-Prager in tension
127c IFORM = 2: Yield formulation flag => Von Mises in tension
128c UPARAM(26) = ITRX = 1 : pressure dependent for all T values
129c ITRX = 2 : no pressure dependency when T < 0
130c UPARAM(27) = IDAM = 1 : damage model without turning point
131c IDAM = 2 : damage model with turning point
132c UPARAM(28) = SSP : sound speed
133c UPARAM(29) = Table Id
134c UPARAM(30) = Xscale for yld function
135c UPARAM(31) = Yscale for yld function
136c UPARAM(32) = Elastic modulus C11
137c UPARAM(33) = Elastic modulus C12
138c--------------------------------
139c UVAR(1) : ARCL = PLA / (1-DAM) or non local plastic strain
140c UVAR(2) : DAM
141c UVAR(3) : TOTAL STRAIN RATE
142c UVAR(4) : plastic function residual (for output/convergence check)
143c UVAR(5) : non local plastic strain
144C=======================================================================
145 young = uparam(1)
146 nu = uparam(2)
147 g = uparam(3)
148 bulk = uparam(4)
149c Parameters of yield function/hardening law and flow potential
150 yld0 = uparam(5)
151 a1f = uparam(9)
152 a2f = uparam(10)
153c A1H = UPARAM(11) ! not used
154c A2H = UPARAM(12) ! not used
155 as = uparam(13)
156c Johnson-Cook failure parameters
157 d1c = uparam(14)
158 d2c = uparam(15)
159 d1f = uparam(16)
160 d2f = uparam(17)
161 d_trx = uparam(18)
162 d_jc = uparam(19)
163 expn = uparam(20)
164c Johnson-Cook strain rate-dependency and distortional hardening
165 cc = uparam(21)
166 epsdref = uparam(22)
167 epsdmax = uparam(23)
168 fcut = uparam(24)
169c
170 itrx = nint(uparam(26))
171 idam = nint(uparam(27))
172 ssp = uparam(28)
173 xscale = uparam(30)
174 yscale = uparam(31)
175 c11 = uparam(32)
176 c12 = uparam(33)
177c
178 func_yld = itable(1)
179 ndim_yld = table(func_yld)%NDIM
180c
181 g2 = g * two
182 alpha = 0.005
183 alphai = one - alpha
184c
185 IF (inloc > 0) THEN
186 dtinv = one / max(timestep, em20)
187 DO i = 1,nel
188 pla_nl(i) = uvar(i,5) + max(dplanl(i),zero)
189 uvar(i,5) = pla_nl(i)
190 dpdt_nl(i) = min(max(dplanl(i),zero)*dtinv ,epsdmax)
191 ENDDO
192 ENDIF
193c------------------------------------------------------------------
194c Computation of the nominal trial stress tensor (non damaged)
195c------------------------------------------------------------------
196 dam(1:nel) = dmg(1:nel)
197 dpla(1:nel) = zero
198c
199 DO i = 1,nel
200 IF (off(i) < 0.001) off(i) = zero
201 IF (off(i) < one) off(i) = off(i)*four_over_5
202 IF (off(i) == one) THEN
203 signxx(i) = sigoxx(i) + c11*depsxx(i) + c12*(depsyy(i) + depszz(i))
204 signyy(i) = sigoyy(i) + c11*depsyy(i) + c12*(depsxx(i) + depszz(i))
205 signzz(i) = sigozz(i) + c11*depszz(i) + c12*(depsxx(i) + depsyy(i))
206 signxy(i) = sigoxy(i) + depsxy(i)*g
207 signyz(i) = sigoyz(i) + depsyz(i)*g
208 signzx(i) = sigozx(i) + depszx(i)*g
209 i1(i) = signxx(i) + signyy(i) + signzz(i)
210c
211 sxx(i) = signxx(i) - i1(i) * third
212 syy(i) = signyy(i) - i1(i) * third
213 szz(i) = signzz(i) - i1(i) * third
214 sxy(i) = signxy(i)
215 syz(i) = signyz(i)
216 szx(i) = signzx(i)
217 j2(i) = (sxx(i)**2 + syy(i)**2 + szz(i)**2)*half
218 . + sxy(i)**2 + syz(i)**2 + szx(i)**2
219 END IF
220 ENDDO
221c Johnson & Cook rate dependency (total strain rate)
222 jcr_log(1:nel) = zero
223 jcrate(1:nel) = one
224 IF (epsdref > zero) THEN
225 IF (inloc == 1) THEN ! non local plastic strain rate
226 jcr_log(1:nel) = log(dpdt_nl(1:nel) / epsdref)
227 jcrate(1:nel) = one + cc * jcr_log(1:nel)
228 ELSE ! total strain rate
229 DO i = 1,nel
230 IF (off(i) == one) THEN
231 epsd(i) = (epspxx(i)**2 + epspyy(i)**2 + epspzz(i)**2 )*two
232 . + epspxy(i)**2 + epspyz(i)**2 + epspzx(i)**2
233 epsd(i) = min(sqrt(epsd(i)) ,epsdmax)
234 epsd(i) = alpha*epsd(i) + alphai*uvar(i,3)
235 uvar(i,3) = epsd(i)
236 IF (epsd(i) > epsdref) THEN
237 jcr_log(i) = log(epsd(i) / epsdref)
238 jcrate(i) = one + cc * jcr_log(i)
239 END IF
240 END IF
241 ENDDO
242 END IF
243 END IF
244c----------------------------------------------------
245c Computation of the initial yield stress
246c----------------------------------------------------
247 IF (ndim_yld == 1) THEN
248 xvec1(1:nel,1) = pla(1:nel)
249 ipos1(1:nel,1) = vartmp(1:nel,1)
250c
251 CALL table_vinterp(table(func_yld),nel,nel,ipos1,xvec1,yld,hardp)
252c
253 yld(1:nel) = yld(1:nel) * yscale * jcrate(1:nel)
254 hardp(1:nel) = hardp(1:nel) * yscale * jcrate(1:nel)
255 vartmp(1:nel,1) = ipos1(1:nel,1)
256
257 ELSE IF (ndim_yld == 2) THEN
258 xvec2(1:nel,1) = pla(1:nel)
259 xvec2(1:nel,2) = epsd(1:nel)
260 ipos2(1:nel,1) = vartmp(1:nel,1)
261 ipos2(1:nel,2) = 1
262c
263 CALL table_vinterp(table(func_yld),nel,nel,ipos2,xvec2,yld,hardp)
264c
265 yld(1:nel) = yld(1:nel) * yscale
266 hardp(1:nel) = hardp(1:nel) * yscale
267 vartmp(1:nel,1) = ipos2(1:nel,1)
268
269 ELSE
270 xvec3(1:nel,1) = pla(1:nel)
271 xvec3(1:nel,2) = epsd(1:nel)
272 xvec3(1:nel,3) = temp(1:nel)
273 ipos3(1:nel,1) = vartmp(1:nel,1)
274 ipos3(1:nel,2) = 1
275 ipos3(1:nel,3) = 1
276c
277 CALL table_vinterp(table(func_yld),nel,nel,ipos3,xvec3,yld,hardp)
278c
279 yld(1:nel) = yld(1:nel) * yscale
280 hardp(1:nel) = hardp(1:nel) * yscale
281 vartmp(1:nel,1) = ipos3(1:nel,1)
282
283 END IF
284c
285 DO i = 1,nel
286 fvm(i) = max(zero, i1(i) + half*sqr3*a1f*yld0 / a2f )
287 phi(i) = j2(i) + third*a2f * fvm(i)**2
288 . - (yld(i)**2 + fourth*(a1f*yld0)**2/a2f)
289 ENDDO
290c------------------------------------------------------------------
291c Check plasticity
292c------------------------------------------------------------------
293 nindx = 0
294 nindf = 0
295 DO i = 1,nel
296 IF (phi(i) > zero .and. off(i) == one) THEN
297 nindx = nindx + 1
298 indx(nindx) = i
299 END IF
300 ENDDO
301c
302 IF (nindx > 0) THEN
303
304 niter = 4
305 dpla(1:nel) = zero
306
307 DO iter = 1,niter
308c
309 DO ii = 1,nindx
310 i = indx(ii)
311
312 ! normal to non associated plastic flow surface = d_psi/d_sig
313 i1p = max(zero, i1(i))
314 jp = third * as * i1p
315 tp = two * jp
316c
317 np_xx = sxx(i) + tp
318 np_yy = syy(i) + tp
319 np_zz = szz(i) + tp
320 np_xy = sxy(i)
321 np_yz = syz(i)
322 np_zx = szx(i)
323c
324 ! normal to plastic yield surface = d_phi/d_sig
325 ip = two_third * a2f * fvm(i)
326 nf_xx = sxx(i) + ip
327 nf_yy = syy(i) + ip
328 nf_zz = szz(i) + ip
329 nf_xy = sxy(i)
330 nf_yz = syz(i)
331 nf_zx = szx(i)
332c
333 ! derivative of yld criterion : dphi/dlam = d_phi/d_sig * dsig/dlam
334 dfdlam =-nf_xx * (c11*np_xx + c12 * (np_yy + np_zz))
335 . - nf_yy * (c11*np_yy + c12 * (np_xx + np_zz))
336 . - nf_zz * (c11*np_zz + c12 * (np_xx + np_yy))
337 . -(nf_xy*np_xy + nf_yz*np_yz + nf_zx*np_zx)*two*g2
338c
339 dpla_dlam = two*sqrt((j2(i) + tp*as*i1(i)))
340 dphi_dlam = dfdlam - two*yld(i)*hardp(i)*dpla_dlam
341c----------------------------------------------------------------
342 dlam = -phi(i) / dphi_dlam
343c----------------------------------------------------------------
344 ! Plastic strains tensor update
345 dpxx = dlam*np_xx
346 dpyy = dlam*np_yy
347 dpzz = dlam*np_zz
348 dpxy = dlam*np_xy
349 dpyz = dlam*np_yz
350 dpzx = dlam*np_zx
351 ! Elasto-plastic stresses update
352 signxx(i) = signxx(i) - (c11*dpxx + c12*(dpyy + dpzz))
353 signyy(i) = signyy(i) - (c11*dpyy + c12*(dpxx + dpzz))
354 signzz(i) = signzz(i) - (c11*dpzz + c12*(dpxx + dpyy))
355 signxy(i) = signxy(i) - g2*dpxy
356 signyz(i) = signyz(i) - g2*dpyz
357 signzx(i) = signzx(i) - g2*dpzx
358c
359 ! Plastic strain increments update
360 ddep = dlam * dpla_dlam
361 dpla(i)= max(zero, dpla(i) + ddep)
362 pla(i) = pla(i) + ddep
363 ENDDO ! II = 1,NINDX
364c-----------------------------------------
365c Update yield from tabulated data
366c-----------------------------------------
367c
368 IF (ndim_yld == 1) THEN
369 DO ii = 1, nindx
370 i = indx(ii)
371 xvec1(ii,1) = pla(i)
372 ipos1(ii,1) = vartmp(i,1)
373 ENDDO
374c
375 CALL table_vinterp(table(func_yld),nel,nel,ipos1,xvec1,yld_i,hardp_i)
376c
377 DO ii = 1, nindx
378 i = indx(ii)
379 yld(i) = yld_i(ii)*yscale
380 hardp(i) = hardp_i(ii)*yscale
381 vartmp(i,1) = ipos1(ii,1)
382 ENDDO
383
384 ELSE IF (ndim_yld == 2) THEN
385 DO ii = 1, nindx
386 i = indx(ii)
387 xvec2(ii,1) = pla(i)
388 xvec2(ii,2) = epsd(i)
389 ipos2(ii,1) = vartmp(i,1)
390 ENDDO
391c
392 CALL table_vinterp(table(func_yld),nel,nindx,ipos2,xvec2,yld_i,hardp_i)
393c
394 DO ii = 1, nindx
395 i = indx(ii)
396 yld(i) = yld_i(ii)*yscale
397 hardp(i) = hardp_i(ii)*yscale
398 vartmp(i,1) = ipos2(ii,1)
399 ENDDO
400
401 ELSE
402 DO ii = 1, nindx
403 i = indx(ii)
404 xvec3(ii,1) = pla(i)
405 xvec3(ii,2) = epsd(i)
406 xvec3(ii,3) = temp(i)
407 ipos3(ii,1) = vartmp(i,1)
408 ENDDO
409c
410 CALL table_vinterp(table(func_yld),nel,nel,ipos3,xvec3,yld_i,hardp_i)
411c
412 DO ii = 1, nindx
413 i = indx(ii)
414 yld(i) = yld_i(ii)*yscale
415 hardp(i) = hardp_i(ii)*yscale
416 vartmp(i,1) = ipos3(ii,1)
417 ENDDO
418 END IF
419c----------------------
420c Update stress invariants and criterion function value
421c----------------------
422 DO ii = 1,nindx
423 i = indx(ii)
424 i1(i) = signxx(i) + signyy(i) + signzz(i)
425 sxx(i) = signxx(i) - i1(i)*third
426 syy(i) = signyy(i) - i1(i)*third
427 szz(i) = signzz(i) - i1(i)*third
428 sxy(i) = signxy(i)
429 syz(i) = signyz(i)
430 szx(i) = signzx(i)
431 j2(i) =(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) * half
432 . + sxy(i)**2 + syz(i)**2 + szx(i)**2
433 fvm(i) = max(zero, i1(i) + half*sqr3*a1f*yld0 / a2f )
434 phi(i) = j2(i) + third*a2f * fvm(i)**2
435 . - (yld(i)**2 + fourth*(a1f*yld0)**2/a2f)
436 uvar(i,4) = phi(i) / yld(i)**2
437 END DO
438c
439 END DO ! Newton iterations
440c
441 END IF ! NINDX > 0
442c---------------------------------------------------------------------
443c Update damage
444c---------------------------------------------------------------------
445 IF (idam > 0) THEN
446 IF (idam == 1) THEN
447 IF (inloc == 1) THEN
448 dgamm(1:nel) = dplanl(1:nel)
449 gamma(1:nel) = pla_nl(1:nel)
450 ELSE
451 dgamm(1:nel) = dpla(1:nel)
452 gamma(1:nel) = pla(1:nel)
453 END IF
454 ELSE
455 IF (inloc == 1) THEN
456 dgamm(1:nel) = dplanl(1:nel) * dmg_scale(1:nel)
457 uvar(1:nel,1) = uvar(1:nel,1) + dgamm(1:nel)
458 gamma(1:nel) = uvar(1:nel,1)
459 ELSE
460 dgamm(1:nel) = dpla(1:nel) * dmg_scale(1:nel)
461 uvar(1:nel,1) = uvar(1:nel,1) + dgamm(1:nel)
462 gamma(1:nel) = uvar(1:nel,1)
463 END IF
464 END IF
465c
466 DO i = 1,nel
467c
468 triax = zero
469 fact = one
470 facr = one + d_jc * jcr_log(i) ! total strain rate factor on damage
471 IF ( j2(i)/= zero)triax = third * i1(i) / sqrt(three*j2(i))
472 IF (itrx == 1 ) THEN
473 fact = exp(-d_trx*triax)
474 ELSE
475 IF (triax > zero )fact = exp(-d_trx*triax)
476 ENDIF
477 gamc = (d1c + d2c * fact) * facr
478 gamf = (d1f + d2f * fact) * facr
479 IF (gamma(i) > gamc) THEN
480 IF (dam(i) == zero) THEN
481 nindf = nindf + 1
482 indf(nindf) = i
483 END IF
484 IF (expn == one) THEN
485 dam(i) = dam(i) + dgamm(i) / (gamf - gamc)
486 ELSE
487 dfac = (gamma(i) - gamc) / (gamf - gamc)
488 dam(i) = dam(i) + expn * dfac**(expn-one) * dgamm(i) / (gamf - gamc)
489 END IF
490 IF (dam(i) >= one) THEN
491 dam(i) = one
492 off(i) = four_over_5
493 END IF
494 dmg(i) = dam(i)
495 END IF ! GAMMA > GAMC
496 ENDDO
497c Update damaged stress
498 DO i = 1, nel
499 dmg_scale(i) = one - dam(i)
500 END DO
501 END IF
502c-------------------------
503 soundsp(1:nel) = ssp
504c-------------------------
505 IF (nindf > 0) THEN
506 DO ii=1,nindf
507#include "lockon.inc"
508 WRITE(iout, 1000) ngl(indf(ii))
509 WRITE(istdo,1100) ngl(indf(ii)),time
510#include "lockoff.inc"
511 ENDDO
512 END IF
513c-----------------------------------------------------------------
514 1000 FORMAT(1x,'START DAMAGE IN SOLID ELEMENT NUMBER ',i10)
515 1100 FORMAT(1x,'START DAMAGE IN SOLID ELEMENT NUMBER ',i10,1x,' AT TIME :',g11.4)
516c-----------------------------------------------------------------
517 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
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