sigeps72.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/.
C  --New-- Hill model + MMC failure model
!||====================================================================
!||    sigeps72   ../engine/source/materials/mat/mat072/sigeps72.F
!||--- called by ------------------------------------------------------
!||    mulaw      ../engine/source/materials/mat_share/mulaw.F90
!||====================================================================
      SUBROUTINE sigeps72(
     1     NEL      ,NUPARAM  ,NUVAR    ,
     2     TIME     ,TIMESTEP ,UPARAM   ,RHO0     ,RHO      ,
     3     DEPSXX   ,DEPSYY   ,DEPSZZ   ,DEPSXY   ,DEPSYZ   ,DEPSZX   ,
     4     SIG0XX   ,SIG0YY   ,SIG0ZZ   ,SIG0XY   ,SIG0YZ   ,SIG0ZX   ,
     5     SIGNXX   ,SIGNYY   ,SIGNZZ   ,SIGNXY   ,SIGNYZ   ,SIGNZX   ,
     6     SOUNDSP  ,UVAR     ,OFF      ,NGL      ,YLD      ,PLA      ,
     7     DPLA     ,ETSE     ,SEQ      ,DMG      ,INLOC    ,DPLANL   )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
#include      "comlock.inc"
C-----------------------------------------------
C   G l o b a l   P a r a m e t e r s
C-----------------------------------------------
#include      "mvsiz_p.inc"
C---------+---------+---+---+--------------------------------------------
#include      "scr17_c.inc"
#include      "com01_c.inc"
#include      "com08_c.inc"
#include      "units_c.inc"
C
      INTEGER NEL,NUPARAM,NUVAR,IPT,NGL(NEL),INLOC
      my_real
     .   TIME,TIMESTEP,UPARAM(NUPARAM),
     .   RHO(NEL),RHO0(NEL),
     .   DEPSXX(NEL),DEPSYY(NEL),DEPSZZ(NEL),
     .   DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
     .   SIG0XX(NEL),SIG0YY(NEL),SIG0ZZ(NEL),
     .   sig0xy(nel),sig0yz(nel),sig0zx(nel)
C-----------------------------------------------
C   O U T P U T   A r g u m e n t s
C-----------------------------------------------
      my_real
     .    signxx(nel),signyy(nel),signzz(nel),
     .    signxy(nel),signyz(nel),signzx(nel),
     .    soundsp(nel),yld(nel),etse(nel),dpla(nel)
C-----------------------------------------------
C   I N P U T   O U T P U T   A r g u m e n t s 
C-----------------------------------------------
      my_real 
     .    uvar(nel,nuvar),off(nel),pla(nel),seq(nel),
     .    dmg(nel),dplanl(nel)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,II,J,NINDX,INDX(MVSIZ),NITER,ITER
      my_real
     .    e,nu,g,g2,a1,a2,bulk,lamhook,
     .    sigy,eps0,nexpo,
     .    ff,gg,hh,nn,ll,mm,
     .    c1,c2,c3,dc,mexp
      my_real
     .    cc(mvsiz),dam(mvsiz),beta(mvsiz),
     .    p,sd11,sd22,sd33,
     .    f1,f2,f3,epsf,det,svm,theta,cos3theta,eta,
     .    sighl(mvsiz),h(mvsiz),dav,
     .    phi(mvsiz),normxx,normyy,normzz,
     .    normxy,normzx,normyz,dpxx(mvsiz),dpyy(mvsiz),
     .    dpzz(mvsiz),dpxy(mvsiz),dpzx(mvsiz),dpyz(mvsiz),
     .    sig_dfdsig,dfdsig2,dphi_dlam(mvsiz),dpla_dlam(mvsiz),
     .    dlam,ddep
C-----------------------------------------------
C     USER VARIABLES INITIALIZATION
C-----------------------------------------------
      ! Elastic parameters
      e       = uparam(1)   ! Young modulus
      nu      = uparam(2)   ! Poisson's ratio
      g       = uparam(3)   ! Shear modulus
      g2      = uparam(4)   ! 2x Shear modulus
      a1      = uparam(7)   ! Diag. component of 3D elastic matrix
      a2      = uparam(8)   ! Non-diag. component of 3D elastic matrix
      bulk    = uparam(9)   ! Bulk modulus
      lamhook = uparam(10)  ! Hooke's lambda coefficient
      ! Hardening parameter
      sigy    = uparam(11)  ! Initial yield stress
      eps0    = uparam(12)  ! Initial plastic strain (> 0)
      nexpo   = uparam(13)  ! Hardening exponent
      ! Hill yield criterion parameters
      ff      = uparam(14)  ! First  Hill coefficient
      gg      = uparam(15)  ! Second Hill coefficient
      hh      = uparam(16)  ! Third  Hill coefficient
      nn      = uparam(17)  ! Fourth Hill coefficient
      ll      = uparam(18)  ! Fifth  Hill coefficient
      mm      = uparam(19)  ! Sixth  Hill coefficient
      ! Modified Mohr-Coulomb yield criterion parameters
      c1      = uparam(20)  ! First failure parameter
      c2      = uparam(21)  ! Second failure parameter
      c3      = uparam(22)  ! Third failure parameter
      mexp    = uparam(23)  ! Damage exponent
      dc      = uparam(24)  ! Critical damage value
c 
      ! Initial value
      IF (isigi == 0) THEN
        IF (time == zero) THEN
          DO i=1,nel
            uvar(i,1) = zero                              
          ENDDO
        ENDIF
      ENDIF
c
      ! Recovering internal variables
      DO i=1,nel
        ! Checking deletion flag value
        IF (off(i) < one)  off(i) = four_over_5*off(i)
        IF (off(i) < em01) off(i) = zero
        ! Hourglass coefficient
        etse(i) = one
        h(i)    = zero
        ! Plastic strain increment
        dpla(i) = zero
        dpxx(i) = zero
        dpyy(i) = zero
        dpzz(i) = zero
        dpxy(i) = zero
        dpyz(i) = zero
        dpzx(i) = zero
        ! Damage variable
        dam(i)  = uvar(i,1)
      ENDDO
c
      ! Computing yield stress and checking damage criteria
      DO i = 1,nel
        ! Computation of the BETA factor
        beta(i) = one
        IF (off(i) == one) THEN
          IF ((dam(i) <= dc) .AND. (dam(i) >= one)) THEN
            beta(i) = one/(max(dc - one,em20))
            beta(i) = (dc - dam(i))*beta(i)
            beta(i) = beta(i)**mexp
          ENDIF
        ENDIF       
        ! Computation of the yield stress
        cc(i)  = pla(i) + eps0
        yld(i) = beta(i)*sigy*exp(nexpo*log(cc(i)))
        yld(i) = max(yld(i), em10)
      ENDDO
c      
      !========================================================================
      ! - computation of trial values
      !======================================================================== 
      DO i = 1,nel
c    
        ! Computation of the trial stress tensor
        dav = (depsxx(i) + depsyy(i) + depszz(i)) * lamhook
        signxx(i) = sig0xx(i) + depsxx(i)*g2 + dav
        signyy(i) = sig0yy(i) + depsyy(i)*g2 + dav
        signzz(i) = sig0zz(i) + depszz(i)*g2 + dav
        signxy(i) = sig0xy(i) + depsxy(i)*g
        signyz(i) = sig0yz(i) + depsyz(i)*g
        signzx(i) = sig0zx(i) + depszx(i)*g
C
        ! Hill yield stress
        sighl(i) = ff*(signyy(i) - signzz(i))**2 + gg*(signzz(i) - signxx(i))**2 +
     .             hh*(signxx(i) - signyy(i))**2 + two*ll*signyz(i)**2 +
     .             two*mm*signzx(i)**2 + two*nn*signxy(i)**2
        sighl(i) = sqrt(max(zero,sighl(i)))
C        
      ENDDO  
c
      !========================================================================
      ! - COMPUTATION OF YIELD FONCTION
      !========================================================================
      phi(1:nel) = sighl(1:nel) - yld(1:nel)
      ! Checking plastic behavior for all elements
      nindx = 0
      indx  = 0
      DO i=1,nel         
        IF ((phi(i)>zero).AND.(off(i) == one)) THEN
          nindx = nindx+1
          indx(nindx) = i
        ENDIF
      ENDDO
c      
      !====================================================================
      ! - PLASTIC CORRECTION WITH CUTTING PLANE ALGORITHM (NEWTON ITERATION)
      !====================================================================       
c      
      ! Number of maximum Newton iterations
      niter = 3
c
      ! Loop over the iterations     
      DO iter = 1, niter
#include "vectorize.inc" 
        ! Loop over yielding elements
        DO ii=1,nindx 
          i = indx(ii)
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 using the backward Euler implicit method
          ! with a Newton-Raphson iterative procedure
          ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
            ! -> PHI          : current value of yield function (known)
          ! -> 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
          !------------------------------------------------------------- 
          
          normxx = (gg*(signxx(i)-signzz(i)) + hh*(signxx(i)-signyy(i)))/(max(sighl(i),em20))
          normyy = (ff*(signyy(i)-signzz(i)) + hh*(signyy(i)-signxx(i)))/(max(sighl(i),em20))
          normzz = (ff*(signzz(i)-signyy(i)) + gg*(signzz(i)-signxx(i)))/(max(sighl(i),em20))
          normxy = two*nn*signxy(i)/(max(sighl(i),em20))
          normyz = two*ll*signyz(i)/(max(sighl(i),em20))
          normzx = two*mm*signzx(i)/(max(sighl(i),em20))
c          
          ! 2 - Computation of DPHI_DLAMBDA
          !---------------------------------------------------------
c        
          !   a) Derivative with respect stress increments tensor DSIG
          !   --------------------------------------------------------
          dfdsig2 = normxx * normxx * g2
     .            + normyy * normyy * g2
     .            + normzz * normzz * g2
     .            + 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
          !     ----------------------------------------------------------------------------
          h(i) = beta(i)*nexpo*sigy*exp((nexpo-1)*log(cc(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 / max(yld(i),em20)
c
          ! 3 - Computation of plastic multiplier and variables update
          !----------------------------------------------------------
c          
          ! Derivative of PHI with respect to DLAM
          dphi_dlam(i) = - dfdsig2 - h(i)*dpla_dlam(i)
          dphi_dlam(i) = sign(max(abs(dphi_dlam(i)),em20) ,dphi_dlam(i))
c          
          ! Computation of the plastic multiplier
          dlam = -phi(i)/dphi_dlam(i)
c         
          ! Plastic strains tensor update
          dpxx(i) = dlam * normxx
          dpyy(i) = dlam * normyy
          dpzz(i) = dlam * normzz
          dpxy(i) = dlam * normxy
          dpyz(i) = dlam * normyz
          dpzx(i) = dlam * normzx  
c          
          ! Elasto-plastic stresses update
          signxx(i) = signxx(i) - dpxx(i)*g2
          signyy(i) = signyy(i) - dpyy(i)*g2
          signzz(i) = signzz(i) - dpzz(i)*g2
          signxy(i) = signxy(i) - dpxy(i)*g
          signyz(i) = signyz(i) - dpyz(i)*g
          signzx(i) = signzx(i) - dpzx(i)*g
c          
          ! Cumulated plastic strain and strain rate update           
          ddep    = dlam*dpla_dlam(i)
          dpla(i) = max(zero, dpla(i) + ddep)
          pla(i)  = max(pla(i) + ddep,zero)
c
          ! Update the hardening yield stress
          cc(i)  = pla(i) + eps0
          yld(i) = beta(i)*sigy*exp(nexpo*log(cc(i)))
          yld(i) = max(yld(i),em10)
c
          ! Update Hill equivalent stress
          sighl(i) = ff*(signyy(i) - signzz(i))**2 + gg*(signzz(i) - signxx(i))**2 +
     .               hh*(signxx(i) - signyy(i))**2 + two*ll*signyz(i)**2 +
     .               two*mm*signzx(i)**2 + two*nn*signxy(i)**2
          sighl(i) = sqrt(max(sighl(i),zero))
c
          ! Update yield function value
          phi(i)   = sighl(i) - yld(i)
c
        ENDDO
        ! End of the loop over the yielding elements
      ENDDO
      ! End of the loop over the iterations 
      !===================================================================
      ! - END OF PLASTIC CORRECTION WITH CUTTING PLANE ALGORITHM
      !===================================================================
C    
      ! Update and store new internal variables
      nindx = 0
      indx  = 0
      DO i=1,nel
        ! Modified Mohr Criteria Failure Model
        IF (off(i) == one) THEN
          ! Pressure
          p    = third*(signxx(i) + signyy(i) + signzz(i))
          ! Von Mises equivalent stress
          sd11 = signxx(i) - p
          sd22 = signyy(i) - p
          sd33 = signzz(i) - p
          svm  = half*(sd11**2 + sd22**2 + sd33**2) + 
     +                   signxy(i)**2 + signzx(i)**2 + signyz(i)**2
          svm  = sqrt(max(three*svm,zero))
          ! Determinant of the deviatoric stress tensor
          det  = sd11*sd22*sd33 + two*signxy(i)*signzx(i)*signyz(i)
     .         - sd11*signyz(i)**2-sd33*signxy(i)**2 - sd22*signzx(i)**2
          ! Computation of triaxiality
          eta = p/max(svm,em20)
          IF (eta < -one) eta = -one
          IF (eta > one ) eta = one
          ! Computation of Lode angle and Lode parameter
          cos3theta = half*twenty7*det/max(em20,svm**3)
          IF (cos3theta < -one) cos3theta = -one
          IF (cos3theta > one)  cos3theta = one
          theta = one - two*acos(cos3theta)/pi
C 
          ! Computation of failure coefficient
          f1 = cos(theta*pi/six)
          f2 = sin(theta*pi/six) 
          f3 = c3 + (sqrt(three)/(two - sqrt(three)))*(one - c3)*(one/max(f1,em20) - one)
C
          ! Computation of the failure plastic strain
          epsf = (sigy/max(c2,em20))*f3*(f1*sqrt(third*(one + c1**2)) + c1*(eta + f2*third))
          epsf = max(epsf,em20)**(-one/nexpo)
C           
          ! Computation of the damage variable
          IF (inloc > 0) THEN 
            dam(i) = dam(i) + max(dplanl(i),zero)/max(epsf,em20)
          ELSE
            dam(i) = dam(i) + dpla(i)/max(epsf,em20)
          ENDIF
          IF (dam(i) >= dc) THEN
            dam(i) = dc
            off(i) = four_over_5
            nindx = nindx+1
            indx(nindx)=i
          ENDIF  
C
          ! Store the new value of damage
          uvar(i,1) = dam(i)
        ENDIF
        ! Equivalent stress
        seq(i)     = sighl(i)
        ! Normalized damage
        dmg(i)     = dam(i)/dc
        ! Sound-speed
        soundsp(i) = sqrt((bulk + four_over_3*g)/rho0(i)) 
        ! Coefficient for hourglass
        IF (dpla(i) > zero) THEN 
          etse(i)  = h(i) / (h(i) + e)
        ELSE
          etse(i)  = one
        ENDIF
      ENDDO 
C
      ! Printout element deletion
      IF(nindx>0)THEN
        DO j=1,nindx
#include "lockon.inc"
          WRITE(iout, 1000) ngl(indx(j))
          WRITE(istdo,1100) ngl(indx(j)),tt
#include "lockoff.inc"
        ENDDO
      ENDIF
!---
 1000 FORMAT(1x,'rupture of solid element number ',I10)
 1100 FORMAT(1X,'rupture of solid element number ',I10,
     .          ' at time :',G11.4)                          
!---
      END
C