sigeps120_vm.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    sigeps120_vm   ../engine/source/materials/mat/mat120/sigeps120_vm.F
!||--- called by ------------------------------------------------------
!||    sigeps120      ../engine/source/materials/mat/mat120/sigeps120.F
!||====================================================================
      SUBROUTINE sigeps120_vm(
     1      NEL     ,NGL     ,NUPARAM ,NUVAR   ,TIME    ,TIMESTEP,
     2      UPARAM  ,UVAR    ,OFF     ,PLA     ,EPSD    ,SOUNDSP ,
     3      EPSPXX  ,EPSPYY  ,EPSPZZ  ,EPSPXY  ,EPSPYZ  ,EPSPZX  ,
     4      DEPSXX  ,DEPSYY  ,DEPSZZ  ,DEPSXY  ,DEPSYZ  ,DEPSZX  ,
     5      SIGOXX  ,SIGOYY  ,SIGOZZ  ,SIGOXY  ,SIGOYZ  ,SIGOZX  ,
     6      SIGNXX  ,SIGNYY  ,SIGNZZ  ,SIGNXY  ,SIGNYZ  ,SIGNZX  ,
     7      INLOC   ,DPLANL  ,DMG     , DMG_SCALE)
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "units_c.inc"
#include      "comlock.inc"
C----------------------------------------------------------------
C  D u m m y   A R G U M E N T S
C----------------------------------------------------------------
      INTEGER :: NEL,NUPARAM,NUVAR,INLOC
      my_real :: TIME,TIMESTEP
      INTEGER ,DIMENSION(NEL)    ,INTENT(IN) :: NGL
      my_real ,DIMENSION(NUPARAM),INTENT(IN) :: UPARAM
      my_real ,DIMENSION(NEL)    ,INTENT(IN) ::
     .         EPSPXX,EPSPYY,EPSPZZ,EPSPXY,EPSPYZ,EPSPZX,
     .         DEPSXX,DEPSYY,DEPSZZ,DEPSXY,DEPSYZ,DEPSZX,
     .         sigoxx,sigoyy,sigozz,sigoxy,sigoyz,sigozx
c
      my_real ,DIMENSION(NEL) ,INTENT(OUT) :: epsd,soundsp
      my_real ,DIMENSION(NEL) ,INTENT(OUT) :: 
     .         signxx,signyy,signzz,signxy,signyz,signzx
      my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: off,pla,dplanl,dmg,dmg_scale
      my_real ,DIMENSION(NEL,NUVAR) ,INTENT(INOUT) :: uvar
C----------------------------------------------------------------
C  L O C A L  V A R I B L E S
C----------------------------------------------------------------
      INTEGER  :: I,II,ITER,NITER,NINDX,NINDF,ITRX,IDAM
      INTEGER  ,DIMENSION(NEL) :: INDX,INDF
c     material parameters
      my_real :: young,nu,ssp,g,g2,bulk,yld0,qq,beta,hh,c11,c12,
     .   a1f,a2f,as,d1c,d2c,d1f,d2f,d_trx,d_jc,dm,expn
c     Johnson & Cook rate-dependency
      my_real cc,epsdref,epsdmax,frate,fcut,alpha,alphai
c     Local variables
      my_real facr,triax,rvoce,fvm,i1p,ip,jp,tp,fact,facd,dfac,dtinv,
     .   dlam,ddep,dfdlam,dyld_dpla,dpla_dlam,dphi_dlam,
     .   gamc,gamf,np_xx,np_yy,np_zz,np_xy,np_yz,np_zx,
     .   nf_xx,nf_yy,nf_zz,nf_xy,nf_yz,nf_zx,dpxx,dpyy,dpzz,dpxy,dpyz,dpzx
      my_real ,DIMENSION(NEL) :: i1,j2,yld,phi,dam,jcrate,
     .   sxx,syy,szz,sxy,syz,szx,dpla,jcr_log,gamma,dgamm,pla_nl,dpdt_nl
c--------------------------------
c     Material parameters
c--------------------------------
c     UPARAM(1)  = Young modulus E
c     UPARAM(2)  = Poisson's ratio nu
c     UPARAM(3)  = shear modulus   G
c     UPARAM(4)  = bulk modulus    K
c     UPARAM(5)  = Yld0: instant yield stress 
c     UPARAM(6)  = Q : expontial term coefficient in the hardening law
c     UPARAM(7)  = BETA : exponent of the exponential term 
c     UPARAM(8)  = P : multiplier of the linear term in the hardening law
c     UPARAM(9)  = A1F : parameter of the yield function
c     UPARAM(10) = A2F : parameter of the yield function
c     UPARAM(11) = A1H : distortionnal yield hardening coefficiant
c     UPARAM(12) = A2H : distortionnal yield hardening coefficiant
c     UPARAM(13) = AS : parameter of the potential function
c     UPARAM(14) = D1C: first  damage strain parameter in initial damage threshold
c     UPARAM(15) = D2C: second damage strain parameter in initial damage threshold 
c     UPARAM(16) = D1F: first  damage strain parameter in final damage threshold
c     UPARAM(17) = D2F: second damage strain parameter in final damage threshold
c     UPARAM(18) = D_TRX : triaxiality factor in damage formula
c     UPARAM(19) = D_JC: JC strain rate coefficient in damage formula
c     UPARAM(20) = EXPN exponent in damage evolution
c     UPARAM(21) = CC : Johnson-Cook strain rate-dependency coefficient
c     UPARAM(22) = EPSDREF quasi-static reference strain rate
c     UPARAM(23) = EPSDMAX maximal reference strain rate
c     UPARAM(24) = FCUT : cut frequency for strain rate filtering
c     UPARAM(25) = IFORM = 1: Yield formulation flag  => Drucker-Prager in tension
c                  IFORM = 2: Yield formulation flag  => Von Mises in tension
c     UPARAM(26) = ITRX = 1 : pressure dependent for all T values
c                  ITRX = 2 : no pressure dependency when T < 0
c     UPARAM(27) = IDAM = 1 : damage model without turning point
c                  IDAM = 2 : damage model with turning point
c     UPARAM(28) = SSP  : sound speed
c     UPARAM(29) = Table Id
c     UPARAM(30) = Xscale for yld function
c     UPARAM(31) = Yscale for yld function
c     UPARAM(32) = Elastic modulus C11
c     UPARAM(33) = Elastic modulus C12
c--------------------------------     
c     UVAR(1)    : ARCL = PLA / (1-DAM) or non local plastic strain
c     UVAR(2)    : DAM
c     UVAR(3)    : TOTAL STRAIN RATE
c     UVAR(4)    : plastic function residual (for output/convergence check)
c     UVAR(5)    : non local plastic strain
C=======================================================================
      young = uparam(1)
      nu    = uparam(2)
      g     = uparam(3)
      bulk  = uparam(4)
c     Parameters of yield function/hardening law and flow potential
      yld0  = uparam(5)
      qq    = uparam(6)
      beta  = uparam(7)
      hh    = uparam(8)
      a1f   = uparam(9) 
      a2f   = uparam(10)
c      A1H   = UPARAM(11)
c      A2H   = UPARAM(12)
      as    = uparam(13)
c     Johnson-Cook failure parameters
      d1c   = uparam(14)
      d2c   = uparam(15)
      d1f   = uparam(16)
      d2f   = uparam(17)
      d_trx = uparam(18)
      d_jc  = uparam(19)
      expn  = uparam(20) 
c     Johnson-Cook strain rate-dependency and distortional hardening
      cc      = uparam(21) 
      epsdref = uparam(22) 
      epsdmax = uparam(23) 
      fcut    = uparam(24)
c
      itrx    = nint(uparam(26)) 
      idam    = nint(uparam(27))
      ssp     = uparam(28)
      c11     = uparam(32)
      c12     = uparam(33)      
c
      g2     = g * two
      alpha  = 0.005
      alphai = one - alpha
c
      IF (inloc > 0) THEN
        dtinv = one / max(timestep, em20)
        DO i = 1,nel
          pla_nl(i)  = uvar(i,5) + max(dplanl(i),zero)
          uvar(i,5)  = pla_nl(i)
          dpdt_nl(i) = min(max(dplanl(i),zero)*dtinv ,epsdmax)
        ENDDO
      ENDIF
c------------------------------------------------------------------
c     Computation of the nominal trial stress tensor (non damaged)
c------------------------------------------------------------------
      dam(1:nel) = dmg(1:nel)  
      dpla(1:nel)   = zero          
c
      DO i = 1,nel
        IF (off(i) < 0.001) off(i) = zero
        IF (off(i) < one)    off(i) = off(i)*four_over_5
        IF (off(i) == one) THEN
          signxx(i) = sigoxx(i)  + c11*depsxx(i) + c12*(depsyy(i) + depszz(i))
          signyy(i) = sigoyy(i)  + c11*depsyy(i) + c12*(depsxx(i) + depszz(i))
          signzz(i) = sigozz(i)  + c11*depszz(i) + c12*(depsxx(i) + depsyy(i))
          signxy(i) = sigoxy(i)  + depsxy(i)*g 
          signyz(i) = sigoyz(i)  + depsyz(i)*g
          signzx(i) = sigozx(i)  + depszx(i)*g
          i1(i)     = signxx(i) + signyy(i) + signzz(i)
c
          sxx(i) = signxx(i) - i1(i) * third
          syy(i) = signyy(i) - i1(i) * third
          szz(i) = signzz(i) - i1(i) * third
          sxy(i) = signxy(i)
          syz(i) = signyz(i)
          szx(i) = signzx(i)
          j2(i)   = (sxx(i)**2 + syy(i)**2 + szz(i)**2)*half
     .            +  sxy(i)**2 + syz(i)**2 + szx(i)**2
        END IF
      ENDDO
c     Johnson & Cook rate dependency  (total strain rate)
      jcr_log(1:nel) = zero
      jcrate(1:nel)  = one
      IF (epsdref > zero) THEN
        IF (inloc == 1) THEN   ! non local plastic strain rate
          jcr_log(1:nel) = log(dpdt_nl(1:nel) / epsdref)
          jcrate(1:nel)  = one + cc * jcr_log(1:nel)
        ELSE                   ! total strain rate
          DO i = 1,nel
            IF (off(i) == one) THEN
              epsd(i) = (epspxx(i)**2 + epspyy(i)**2 + epspzz(i)**2 )*two
     .                +  epspxy(i)**2 + epspyz(i)**2 + epspzx(i)**2
              epsd(i) = min(sqrt(epsd(i)) ,epsdmax)
              epsd(i) = alpha*epsd(i) + alphai*uvar(i,3)
              uvar(i,3) = epsd(i)
              IF (epsd(i) > epsdref) THEN
                jcr_log(i) = log(epsd(i) / epsdref)
                jcrate(i)  = one + cc * jcr_log(i)
              END IF
            END IF
          ENDDO
        END IF
      END IF
c------------------------------------------------------------------
c     Actual yield value, yield criterion and plastic flow function
c------------------------------------------------------------------
      DO i = 1,nel
        rvoce  = qq*(one - exp(-beta*pla(i))) + hh*pla(i)
        yld(i) = (yld0 + rvoce) * jcrate(i)
        fvm    = max(zero, i1(i) + half*sqr3*a1f*yld0 / a2f )
        phi(i) = j2(i) + third*a2f * fvm**2 
     .         - (yld(i)**2 + fourth*(a1f*yld0)**2/a2f)
      ENDDO
c------------------------------------------------------------------
c     Plasticity
c------------------------------------------------------------------
      nindx = 0
      nindf = 0
      DO i = 1,nel
        IF (phi(i) > zero .and. off(i) == one) THEN
          nindx = nindx + 1
          indx(nindx) = i
        END IF
      ENDDO
c
      DO ii = 1,nindx
        i = indx(ii)
        dpla(i) = zero
        i1p   = max(zero, i1(i))
        fvm   = max(zero, i1(i) + half*sqr3*a1f*yld0 / a2f)
        niter = 4
c
        DO iter = 1,niter
 
          ! normal to non associated plastic flow surface = d_psi/d_sig
          jp    = third * as * i1p
          tp    = two  * jp
c
          np_xx = sxx(i) + tp
          np_yy = syy(i) + tp
          np_zz = szz(i) + tp
          np_xy = sxy(i)
          np_yz = syz(i)
          np_zx = szx(i)
c
          ! normal to plastic yield surface = d_phi/d_sig
          ip = two_third * a2f * fvm
          nf_xx = sxx(i) + ip
          nf_yy = syy(i) + ip
          nf_zz = szz(i) + ip
          nf_xy = sxy(i)
          nf_yz = syz(i)
          nf_zx = szx(i)
c         
          ! derivative of yld criterion : dphi/dlam = d_phi/d_sig * dsig/dlam
          dfdlam =-nf_xx * (c11*np_xx + c12 * (np_yy + np_zz))
     .           - nf_yy * (c11*np_yy + c12 * (np_xx + np_zz))
     .           - nf_zz * (c11*np_zz + c12 * (np_xx + np_yy))
     .           -(nf_xy*np_xy + nf_yz*np_yz + nf_zx*np_zx)*two*g2
c
          dyld_dpla = (hh + beta*qq*exp(-beta*pla(i))) * jcrate(i)
          dpla_dlam = two*sqrt((j2(i) + tp*as*i1(i)))
          dphi_dlam = dfdlam - two*yld(i)*dyld_dpla*dpla_dlam
c----------------------------------------------------------------
          dlam = -phi(i) / dphi_dlam
c----------------------------------------------------------------
          ! Plastic strains tensor update
          dpxx = dlam*np_xx                          
          dpyy = dlam*np_yy                          
          dpzz = dlam*np_zz                          
          dpxy = dlam*np_xy                          
          dpyz = dlam*np_yz                          
          dpzx = dlam*np_zx                 
          ! Elasto-plastic stresses update   
          signxx(i) = signxx(i) - (c11*dpxx + c12*(dpyy + dpzz)) 
          signyy(i) = signyy(i) - (c11*dpyy + c12*(dpxx + dpzz))
          signzz(i) = signzz(i) - (c11*dpzz + c12*(dpxx + dpyy))
          signxy(i) = signxy(i) - g2*dpxy     
          signyz(i) = signyz(i) - g2*dpyz  
          signzx(i) = signzx(i) - g2*dpzx
c
          ! Plastic strain increments update
          ddep   = dlam * dpla_dlam
          dpla(i)= max(zero, dpla(i) + ddep)
          pla(i) = pla(i) + ddep
c----------------------        
c         criterion update  
c----------------------        
          i1(i)  = signxx(i) + signyy(i) + signzz(i)
          i1p    = max(zero, i1(i))
          sxx(i) = signxx(i) - i1(i)*third
          syy(i) = signyy(i) - i1(i)*third
          szz(i) = signzz(i) - i1(i)*third
          sxy(i) = signxy(i)
          syz(i) = signyz(i)
          szx(i) = signzx(i)
          j2(i)  =(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) * half
     .           + sxy(i)**2 + syz(i)**2 + szx(i)**2
          rvoce  = qq*(one - exp(-beta*pla(i))) + hh*pla(i)
          yld(i) = (yld0 + rvoce) * jcrate(i)
          fvm    = max(zero, i1(i) + half*sqr3*a1f*yld0 / a2f)
          phi(i) = j2(i) + third*a2f * fvm**2 
     .           - (yld(i)**2 + fourth*(a1f*yld0)**2/a2f)
        END DO  ! Newton iterations
c
        uvar(i,4) = phi(i) / yld0**2
      ENDDO   ! II = 1,NINDX     
c---------------------------------------------------------------------
c     Update damage
c---------------------------------------------------------------------
      IF (idam > 0) THEN
        IF (idam == 1) THEN
          IF (inloc == 1) THEN
            dgamm(1:nel) = dplanl(1:nel)
            gamma(1:nel) = pla_nl(1:nel)
          ELSE
            dgamm(1:nel) = dpla(1:nel)
            gamma(1:nel) = pla(1:nel)
          END IF
        ELSE
          IF (inloc == 1) THEN
            dgamm(1:nel)  = dplanl(1:nel) * dmg_scale(1:nel)  
            uvar(1:nel,1) = uvar(1:nel,1) + dgamm(1:nel)
            gamma (1:nel) = uvar(1:nel,1)
          ELSE
            dgamm(1:nel)  = dpla(1:nel) * dmg_scale(1:nel)  
            uvar(1:nel,1) = uvar(1:nel,1) + dgamm(1:nel)
            gamma(1:nel) = uvar(1:nel,1)
          END IF
        END IF
c
        DO i = 1,nel
c
          triax = zero
          fact = one
          facr = one + d_jc * jcr_log(i)  ! total strain rate factor on damage
          IF ( j2(i)/= zero)triax = third * i1(i) / sqrt(three*j2(i))
          IF (itrx == 1 ) THEN
            fact  = exp(-d_trx*triax)
          ELSE
           IF (triax > zero )fact  = exp(-d_trx*triax)
          ENDIF
          gamc = (d1c + d2c * fact) * facr
          gamf = (d1f + d2f * fact) * facr
          IF (gamma(i) > gamc) THEN
            IF (dam(i) == zero) THEN
              nindf = nindf + 1
              indf(nindf) = i
            END IF
            IF (expn == one) THEN
              dam(i) = dam(i) + dgamm(i) / (gamf - gamc)
            ELSE
              dfac   = (gamma(i) - gamc) / (gamf - gamc)
              dam(i) = dam(i) + expn * dfac**(expn-one) * dgamm(i) / (gamf - gamc)
            END IF
            IF (dam(i) >= one) THEN
              dam(i) = one
              off(i) = four_over_5
            END IF
            dmg(i) = dam(i)
          END IF ! GAMMA > GAMC
        ENDDO
c       Update damaged stress
        DO i = 1, nel
          dmg_scale(i) = one - dam(i)   
        END DO
      END IF
c-------------------------
      soundsp(1:nel) = ssp
c-------------------------
      IF (nindf > 0) THEN
        DO ii=1,nindf
#include "lockon.inc"
          WRITE(iout, 1000) ngl(indf(ii))
          WRITE(istdo,1100) ngl(indf(ii)),time
#include "lockoff.inc"
        ENDDO
      END IF
c-----------------------------------------------------------------
 1000 FORMAT(1x,'START DAMAGE IN SOLID ELEMENT NUMBER ',i10)
 1100 FORMAT(1x,'START DAMAGE IN SOLID ELEMENT NUMBER ',i10,1x,' AT TIME :',g11.4) 
c-----------------------------------------------------------------
      RETURN
      END