mat115_nice.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/.
!||====================================================================
!||    mat115_nice   ../engine/source/materials/mat/mat115/mat115_nice.F
!||--- called by ------------------------------------------------------
!||    sigeps115     ../engine/source/materials/mat/mat115/sigeps115.F
!||====================================================================
      SUBROUTINE mat115_nice(
     1     NEL     ,NGL     ,NUPARAM ,NUVAR   ,GRHO    , 
     2     TIME    ,TIMESTEP,UPARAM  ,UVAR    ,OFF     ,SIGY    ,
     3     RHO0    ,PLA     ,DPLA    ,SOUNDSP ,ET      ,SEQ     ,
     4     DEPSXX  ,DEPSYY  ,DEPSZZ  ,DEPSXY  ,DEPSYZ  ,DEPSZX  ,
     5     SIGOXX  ,SIGOYY  ,SIGOZZ  ,SIGOXY  ,SIGOYZ  ,SIGOZX  ,
     6     SIGNXX  ,SIGNYY  ,SIGNZZ  ,SIGNXY  ,SIGNYZ  ,SIGNZX  )
      !=======================================================================
      !      Implicit types
      !=======================================================================
#include      "implicit_f.inc"
      !=======================================================================
      !      Common
      !=======================================================================
      !=======================================================================
      !      Dummy arguments
      !=======================================================================
      INTEGER NEL,NUPARAM,NUVAR
      INTEGER ,DIMENSION(NEL), INTENT(IN)    :: NGL
      my_real 
     .   TIME,TIMESTEP
      my_real,DIMENSION(NUPARAM), INTENT(IN) :: 
     .   UPARAM
      my_real,DIMENSION(2*NEL), INTENT(IN)     :: 
     .   grho
      my_real,DIMENSION(NEL), INTENT(IN)     :: 
     .   rho0,
     .   depsxx,depsyy,depszz,depsxy,depsyz,depszx,
     .   sigoxx,sigoyy,sigozz,sigoxy,sigoyz,sigozx
c
      my_real ,DIMENSION(NEL), INTENT(OUT)   :: 
     .   soundsp,sigy,et,
     .   signxx,signyy,signzz,signxy,signyz,signzx
c
      my_real ,DIMENSION(NEL), INTENT(INOUT)       :: 
     .   pla,dpla,off,seq
      my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) :: 
     .   uvar
c
      !=======================================================================
      !      Local Variables
      !=======================================================================
      INTEGER I,K,II,NINDX,INDEX(NEL),Istat
c
       my_real ::  
     .   YOUNG,BULK,LAMHOOK,G,G2,NU,ALPHA,CFAIL,PFAIL,RHOF0
c
      my_real, DIMENSION(NEL) ::
     .   gamma,epsd,alpha2,beta,sigp
c
      real(kind=8) :: ldav
      my_real :: trdep,dlam,ddep,dphi_dsigvm,dphi_dsigm,
     .   dpxx,dpyy,dpzz,dpxy,dpyz,dpzx,trdfds,
     .   normxx,normyy,normzz,normxy,normyz,normzx,
     .   dphi,sig_dfdsig,dfdsig2,dphi_dpla,
     .   i1,i2,i3,q,r,r_inter,psi,s11,s22,s33
c     
      my_real, DIMENSION(NEL) ::
     .   sigvm,sigeq,dsigxx,dsigyy,dsigzz,dsigxy,dsigyz,dsigzx,trsig,
     .   sxx,syy,szz,sxy,syz,szx,sigm,j2,yld,hardp,phi0,phi,dpla_dlam,
     .   defvp,dphi_dlam
      !=======================================================================
      !       DRUCKER MATERIAL LAW WITH NO DAMAGE
      !=======================================================================
      !UVAR(1)   PHI    YIELD FUNCTION VALUE
      !-  
      !DEPIJ     PLASTIC STRAIN TENSOR COMPONENT
      !DEPSIJ    TOTAL   STRAIN TENSOR COMPONENT (EL+PL)
      !=======================================================================
c
      !=======================================================================
      ! - INITIALISATION OF COMPUTATION ON TIME STEP
      !=======================================================================
      ! Recovering model parameters
      ! Elastic parameters      
      young   = uparam(1)  ! Young modulus
      bulk    = uparam(2)  ! Bulk modulus
      g       = uparam(3)  ! Shear modulus
      g2      = uparam(4)  ! 2*Shear modulus
      lamhook = uparam(5)  ! Lambda Hooke parameter
      nu      = uparam(6)  ! Poisson ration
      ! Statistic variation flag
      istat   = uparam(12)
      ! Plastic criterion and hardening parameters
      alpha   = uparam(13) ! Yield function shape parameter
      cfail   = uparam(14) ! Tensile volumic strain at failure
      pfail   = uparam(15) ! Principal stress at failure  
      ! Density of base material
      IF (istat .NE. 0) rhof0 = uparam(16) 
c
      ! Recovering internal variables
      DO i=1,nel
        IF (off(i) < 0.1) off(i) = zero
        IF (off(i) < one)  off(i) = off(i)*four_over_5
        ! User inputs
        phi0(i)  = uvar(i,1) ! Previous yield function value
        ! Standard inputs
        defvp(i) = zero      ! Initialization of the volumic plastic strain increment
        dpla(i)  = zero      ! Initialization of the plastic strain increment
        et(i)    = one        ! Initialization of hourglass coefficient
        hardp(i) = zero      ! Initialization of hourglass coefficient
        IF (istat == 0) THEN
          gamma(i)   = uparam(16) ! Yield stress parameter
          epsd(i)    = uparam(17) ! Densification strain
          alpha2(i)  = uparam(18) ! Yield stress parameter
          beta(i)    = uparam(19) ! Yield stress parameter
          sigp(i)    = uparam(20) ! Initial yield stress
        ELSE
          sigp(i)    = uparam(17) + uparam(18)*((grho(i+nel)/rhof0)**uparam(19))
          alpha2(i)  = uparam(20) + uparam(21)*((grho(i+nel)/rhof0)**uparam(22))
          gamma(i)   = uparam(23) + uparam(24)*((grho(i+nel)/rhof0)**uparam(25))
          beta(i)    = one/(uparam(26) + uparam(27)*((grho(i+nel)/rhof0)**uparam(28)))
          epsd(i)    = (-(nine + (alpha**2))/(three*(alpha**2)))*log(grho(i+nel)/rhof0)
        ENDIF
      ENDDO    
c      
      ! Computation of the initial yield stress
      DO i = 1,nel
        ! a) - Hardening law
        yld(i) = sigp(i) + gamma(i)*(pla(i)/epsd(i)) + 
     .                alpha2(i)*log(one/max((one - (pla(i)/epsd(i))**beta(i)),em20))
        ! b) - Checking values
        yld(i) = max(em10,yld(i))
      ENDDO     
c      
      !========================================================================
      ! - COMPUTATION OF TRIAL VALUES
      !========================================================================
      DO i=1,nel
        ! Computation of the trial stress tensor
        ldav = (depsxx(i) + depsyy(i) + depszz(i)) * lamhook
        signxx(i) = sigoxx(i) + depsxx(i)*g2 + ldav
        signyy(i) = sigoyy(i) + depsyy(i)*g2 + ldav
        signzz(i) = sigozz(i) + depszz(i)*g2 + ldav
        signxy(i) = sigoxy(i) + depsxy(i)*g
        signyz(i) = sigoyz(i) + depsyz(i)*g
        signzx(i) = sigozx(i) + depszx(i)*g
        ! Computation of the trace of the trial stress tensor
        trsig(i)  = signxx(i) + signyy(i) + signzz(i)
        sigm(i)   = trsig(i) * third
        ! Computation of the deviatoric trial stress tensor
        sxx(i) = signxx(i) - sigm(i)
        syy(i) = signyy(i) - sigm(i)
        szz(i) = signzz(i) - sigm(i)
        sxy(i) = signxy(i)
        syz(i) = signyz(i)
        szx(i) = signzx(i)
        ! Second deviatoric invariant and Mises equivalent stress
        j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 )
     .          +       sxy(i)**2 + syz(i)**2 + szx(i)**2
        sigvm(i) = sqrt(three*j2(i))
        ! Deshpande - Fleck equivalent stress
        sigeq(i) = sqrt((sigvm(i)**2 + (alpha*sigm(i))**2)/(one + (alpha/three)**2))
      ENDDO
c
      !========================================================================
      ! - COMPUTATION OF YIELD FONCTION
      !========================================================================
      phi(1:nel) = sigeq(1:nel) - yld(1:nel)
c
      ! Checking plastic behavior for all elements
      nindx = 0
      DO i=1,nel         
        IF ((phi(i)>zero).AND.(off(i) == one)) THEN
          nindx=nindx+1
          index(nindx)=i
        ENDIF
      ENDDO
c      
      !====================================================================
      ! - PLASTIC CORRECTION WITH NICE - EXPLICIT 1 METHOD
      !==================================================================== 
#include "vectorize.inc" 
      DO ii = 1, nindx  
c      
        ! Number of the element with plastic behaviour                                          
        i = index(ii)    
c        
        ! Computation of the trial stress increment
        ldav = (depsxx(i) + depsyy(i) + depszz(i)) * lamhook
        dsigxx(i) = depsxx(i)*g2 + ldav
        dsigyy(i) = depsyy(i)*g2 + ldav
        dsigzz(i) = depszz(i)*g2 + ldav
        dsigxy(i) = depsxy(i)*g   
        dsigyz(i) = depsyz(i)*g     
        dsigzx(i) = depszx(i)*g   
        dlam      = zero        
c
        ! Computation of Von Mises equivalent stress of the previous stress tensor
        trsig(i)  = sigoxx(i) + sigoyy(i) + sigozz(i)
        sigm(i)   = trsig(i) * third
        sxx(i)    = sigoxx(i) - sigm(i)
        syy(i)    = sigoyy(i) - sigm(i)
        szz(i)    = sigozz(i) - sigm(i)
        sxy(i)    = sigoxy(i)
        syz(i)    = sigoyz(i)
        szx(i)    = sigozx(i)
        j2(i)     = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 )
     .            +         sxy(i)**2 + syz(i)**2 + szx(i)**2
        sigvm(i)  = sqrt(three*j2(i))
        ! Deshpande - Fleck equivalent stress
        sigeq(i)  = sqrt((sigvm(i)**2 + (alpha*sigm(i))**2)/(one + (alpha/three)**2))
 
c        
        ! Note     : in this part, the purpose is to compute for each iteration
        ! a plastic multiplier allowing to update internal variables to satisfy
        ! the consistency condition. 
        ! Its expression is : DLAMBDA = - (PHI + DPHI) / DPHI_DLAMBDA
          ! -> PHI  : current value of yield function (known)
        ! -> DPHI : yield function prediction (to compute)
        ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
        !                   into account of internal variables kinetic : 
        !                   plasticity, temperature and damage (to compute)
c        
          ! 1 - Computation of DPHI_DSIG the normal to the yield surface
        !-------------------------------------------------------------     
c
        ! Derivative with respect to the Von Mises equivalent stress
        dphi_dsigvm = sigvm(i)/max((sigeq(i)*(one + (alpha/three)**2)),em20)
c          
        ! Derivative with respect to the pressure stress
        dphi_dsigm  = (alpha**2)*sigm(i)/max((sigeq(i)*(one + (alpha/three)**2)),em20)         
c                  
        ! dPhi/dSig
        normxx  = dphi_dsigvm*three_half*sxx(i)/(max(sigvm(i),em20)) + dphi_dsigm*third
        normyy  = dphi_dsigvm*three_half*syy(i)/(max(sigvm(i),em20)) + dphi_dsigm*third
        normzz  = dphi_dsigvm*three_half*szz(i)/(max(sigvm(i),em20)) + dphi_dsigm*third
        normxy  = two*dphi_dsigvm*three_half*sxy(i)/(max(sigvm(i),em20))
        normyz  = two*dphi_dsigvm*three_half*syz(i)/(max(sigvm(i),em20))
        normzx  = two*dphi_dsigvm*three_half*szx(i)/(max(sigvm(i),em20)) 
c    
        ! Restoring previous value of the yield function
        phi(i) = phi0(i)
c
        ! Computation of yield surface trial increment DPHI       
        dphi = normxx * dsigxx(i)  
     .       + normyy * dsigyy(i)  
     .       + normzz * dsigzz(i)  
     .       + normxy * dsigxy(i)
     .       + normyz * dsigyz(i)
     .       + normzx * dsigzx(i)
c        
        ! 2 - Computation of DPHI_DLAMBDA
        !---------------------------------------------------------
c        
        !   a) Derivative with respect stress increments tensor DSIG
        !   --------------------------------------------------------
        
        trdfds  = normxx + normyy + normzz    
        dfdsig2 = normxx * (normxx*g2 + lamhook*trdfds)
     .          + normyy * (normyy*g2 + lamhook*trdfds) 
     .          + normzz * (normzz*g2 + lamhook*trdfds)   
     .          + normxy * normxy * g
     .          + normyz * normyz * g
     .          + normzx * normzx * g
c
        !   b) Derivatives with respect to plastic strain P 
        !   ------------------------------------------------
c        
        !     i) Derivative of the yield stress with respect to plastic strain dYLD / dPLA
        !     ----------------------------------------------------------------------------
        hardp(i)   = (gamma(i)/epsd(i)) + 
     .      alpha2(i)*((one-(beta(i)/epsd(i))*(pla(i)/epsd(i))**(beta(i)-1))/
     .      max((one - (pla(i)/epsd(i))**beta(i)),em20))    
        dphi_dpla  = hardp(i)
c        
        !     ii) Derivative of dPLA with respect to DLAM
        !     -------------------------------------------
        sig_dfdsig = signxx(i) * normxx
     .             + signyy(i) * normyy
     .             + signzz(i) * normzz
     .             + signxy(i) * normxy
     .             + signyz(i) * normyz
     .             + signzx(i) * normzx         
        dpla_dlam(i) = sig_dfdsig / yld(i)   
c        
        !  d) Derivative with respect to the yield stress
        !  ----------------------------------------------
c        
        !  e) Derivative of PHI with respect to DLAM ( = -DENOM)
        dphi_dlam(i) = - dfdsig2 - dphi_dpla*dpla_dlam(i)
        dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20) ,dphi_dlam(i))
c        
        ! 3 - Computation of plastic multiplier and variables update
        !----------------------------------------------------------
c
        ! computation of the plastic multiplier
        dlam  = -(dphi + phi(i)) / dphi_dlam(i)
        dlam  = max(dlam, zero)
c        
        ! Plastic strains tensor update
        dpxx  = dlam * normxx
        dpyy  = dlam * normyy
        dpzz  = dlam * normzz
        dpxy  = dlam * normxy
        dpyz  = dlam * normyz
        dpzx  = dlam * normzx  
        trdep = dpxx + dpyy + dpzz
c        
        ! Cumulated plastic strain and strain rate update
        ddep     = (dlam/yld(i))*sig_dfdsig
        dpla(i)  = dpla(i) + ddep 
        pla(i)   = pla(i)  + ddep
        defvp(i) = defvp(i) + trdep
c
        ! Elasto-plastic stresses update   
        ldav      = trdep * lamhook
        signxx(i) = sigoxx(i) + dsigxx(i) - (dpxx*g2 + ldav)
        signyy(i) = sigoyy(i) + dsigyy(i) - (dpyy*g2 + ldav)
        signzz(i) = sigozz(i) + dsigzz(i) - (dpzz*g2 + ldav)
        signxy(i) = sigoxy(i) + dsigxy(i) - dpxy*g
        signyz(i) = sigoyz(i) + dsigyz(i) - dpyz*g
        signzx(i) = sigozx(i) + dsigzx(i) - dpzx*g
        trsig(i)  = signxx(i) + signyy(i) + signzz(i)
        sigm(i)   = trsig(i) * third
        sxx(i)    = signxx(i) - sigm(i)
        syy(i)    = signyy(i) - sigm(i)
        szz(i)    = signzz(i) - sigm(i)
        sxy(i)    = signxy(i)
        syz(i)    = signyz(i)
        szx(i)    = signzx(i)
        j2(i)     = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 )
     .            +         sxy(i)**2 + syz(i)**2 + szx(i)**2
        sigvm(i)  = sqrt(three*j2(i))
        ! Deshpande - Fleck equivalent stress
        sigeq(i)  = sqrt((sigvm(i)**2 + (alpha*sigm(i))**2)/(one + (alpha/three)**2))
c        
        ! Yield stress update
        yld(i) = sigp(i) + gamma(i)*(pla(i)/epsd(i)) + 
     .                alpha2(i)*log(one/max((one - (pla(i)/epsd(i))**beta(i)),em20))
        yld(i) = max(yld(i),em10)
c
        ! Yield function value update
        phi(i) = sigeq(i) - yld(i)       
      ENDDO 
      !===================================================================
      ! - END OF PLASTIC CORRECTION WITH NICE - EXPLICIT 1 METHOD
      !===================================================================
c
      ! Storing new values
      DO i=1,nel
        ! USR Outputs
        IF (dpla(i) > zero) THEN 
          uvar(i,1) = phi(i)   ! Yield function value
        ELSE
          uvar(i,1) = zero
        ENDIF
        IF (cfail > zero) THEN 
          uvar(i,2) = uvar(i,2) + defvp(i)
          IF (uvar(i,2) >= cfail) THEN
            IF (off(i) == one) off(i) = four_over_5
          ENDIF
        ENDIF
        IF (pfail > zero) THEN 
          ! Computing the principal stresses
          i1 = signxx(i)+signyy(i)+signzz(i)
          i2 = signxx(i)*signyy(i)+signyy(i)*signzz(i)+signzz(i)*signxx(i)-
     .         signxy(i)*signxy(i)-signzx(i)*signzx(i)-signyz(i)*signyz(i)
          i3 = signxx(i)*signyy(i)*signzz(i)-signxx(i)*signyz(i)*signyz(i)-
     .         signyy(i)*signzx(i)*signzx(i)-signzz(i)*signxy(i)*signxy(i)+
     .         two*signxy(i)*signzx(i)*signyz(i)
          q  = (three*i2 - i1*i1)/nine
          r  = (two*i1*i1*i1-nine*i1*i2+twenty7*i3)/cinquante4     ! (2*I3^3-9*I1*I2+27*I3)/54
          r_inter = min(r/sqrt(max(em20,(-q**3))),one)
          psi = acos(max(r_inter,-one))
 
          s11 = two*sqrt(-q)*cos(psi/three)+third*i1
          s22 = two*sqrt(-q)*cos((psi+two*pi)/three)+third*i1
          s33 = two*sqrt(-q)*cos((psi+four*pi)/three)+third*i1
          s11 = max(s11,s22,s33)
          ! Failure criterion
          IF (s11 >= pfail) THEN 
            IF (off(i) == one) off(i) = four_over_5
          ENDIF
        ENDIF
        seq(i)     = sigeq(i) ! SIGEQ
        ! Coefficient for hourglass
        IF (dpla(i) > zero) THEN 
          et(i)    = hardp(i) / (hardp(i) + young)
        ELSE
          et(i)    = one
        ENDIF
        ! Computation of the sound speed
        soundsp(i) = sqrt((bulk + four_over_3*g)/grho(nel+i))
        ! Storing the yield stress
        sigy(i)    = yld(i)
      ENDDO
c
      END