sigeps72c.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/.
!||====================================================================
!||    sigeps72c   ../engine/source/materials/mat/mat072/sigeps72c.F
!||--- called by ------------------------------------------------------
!||    mulawc      ../engine/source/materials/mat_share/mulawc.F90
!||====================================================================
      SUBROUTINE sigeps72c(
     1     NEL      ,NUPARAM  ,NUVAR    ,
     2     TIME     ,TIMESTEP ,UPARAM   ,RHO0     ,THKLY    ,
     3     DEPSXX   ,DEPSYY   ,DEPSXY   ,DEPSYZ   ,DEPSZX   ,
     4     SIGOXX   ,SIGOYY   ,SIGOXY   ,SIGOYZ   ,SIGOZX   ,
     5     SIGNXX   ,SIGNYY   ,SIGNXY   ,SIGNYZ   ,SIGNZX   ,
     6     SOUNDSP  ,THK      ,PLA      ,UVAR     ,OFF      ,
     7     ETSE     ,GS       ,YLD      ,HARDM    ,SEQ      ,
     8     DPLA     ,DMG      ,INLOC    ,DPLANL   ,LOFF     )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.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---------+---------+---+---+--------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com01_c.inc"
C-----------------------------------------------
C   I N P U T   A r g u m e n t s
C-----------------------------------------------
C
      INTEGER NEL,NUPARAM,NUVAR,INLOC
      my_real
     .   TIME,TIMESTEP(NEL),UPARAM(NUPARAM),
     .   RHO0(NEL),THKLY(NEL),PLA(NEL),
     .   DEPSXX(NEL),DEPSYY(NEL),
     .   DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
     .   SIGOXX(NEL),SIGOYY(NEL),
     .   SIGOXY(NEL),SIGOYZ(NEL),SIGOZX(NEL),
     .   gs(*),hardm(nel),seq(nel),dmg(nel),
     .   dplanl(nel)
      my_real, DIMENSION(NEL), INTENT(IN) :: loff
C-----------------------------------------------
C   O U T P U T   A r g u m e n t s
C-----------------------------------------------
      my_real
     .    signxx(nel),signyy(nel),
     .    signxy(nel),signyz(nel),signzx(nel),
     .    soundsp(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),thk(nel),yld(nel)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,II,J,NMAX,NINDX,INDEX(MVSIZ),NITER,ITER
      my_real
     . E,NU,G,A11,A12,
     . sigy,eps0,nexp,
     . ff,gg,hh,nn,
     . c1,c2,c3,mexp,dc
      my_real
     . eta,cos3theta,theta,f1,f2,f3,cc(mvsiz),
     . beta(mvsiz),dam(mvsiz),dezz(mvsiz),
     . epsf,p,sighl(mvsiz),h(mvsiz),svm(mvsiz),
     . normxx,normyy,normxy,dpxx(mvsiz),dpyy(mvsiz),
     . dpzz(mvsiz),dpxy(mvsiz),deelzz(mvsiz),
     . dpla_dlam(mvsiz),dlam,sig_dfdsig,dfdsig2,
     . phi(mvsiz),dphi_dlam(mvsiz),ddep
C-----------------------------------------------
C     USER VARIABLES INITIALIZATION
C-----------------------------------------------
C             
      ! Elastic parameters
      e    = uparam(1)   ! Young modulus
      nu   = uparam(2)   ! Poisson's ratio
      g    = uparam(3)   ! Shear modulus
      a11  = uparam(5)   ! Diag. component of plane stress elastic matrix
      a12  = uparam(6)   ! Non diag component of plane stress elastic matrix
      ! Hardening parameters
      sigy = uparam(11)  ! Initial yield stress 
      eps0 = uparam(12)  ! Initial plastic strain (> 0)
      nexp = 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
      ! 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
        ! 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(nexp*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
        signxx(i) = sigoxx(i) + a11*depsxx(i)+a12*depsyy(i)
        signyy(i) = sigoyy(i) + a12*depsxx(i)+a11*depsyy(i)
        signxy(i) = sigoxy(i) + g*depsxy(i)
        signyz(i) = sigoyz(i) + gs(i)*depsyz(i)
        signzx(i) = sigozx(i) + gs(i)*depszx(i)
C
        ! Hill equivalent stress
        sighl(i) = (ff+hh)*signyy(i)**2 + (gg+hh)*signxx(i)**2
     .           - two*hh*signxx(i)*signyy(i) + 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
      index = 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 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 = index(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) + hh*(signxx(i)-signyy(i)))/(max(sighl(i),em20))
          normyy = (ff*signyy(i) + hh*(signyy(i)-signxx(i)))/(max(sighl(i),em20))
          normxy = two*nn*signxy(i)/(max(sighl(i),em20))
c          
          ! 2 - Computation of DPHI_DLAMBDA
          !---------------------------------------------------------
c        
          !   a) Derivative with respect stress increments tensor DSIG
          !   --------------------------------------------------------
          dfdsig2 = normxx * (a11*normxx + a12*normyy)
     .            + normyy * (a11*normyy + a12*normxx)
     .            + normxy * normxy * 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)*nexp*sigy*exp((nexp-1)*log(cc(i)))
c
          !     ii) Derivative of dPLA with respect to DLAM
          !     -------------------------------------------   
          sig_dfdsig = signxx(i) * normxx
     .               + signyy(i) * normyy
     .               + signxy(i) * normxy          
          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
          dpxy(i) = dlam * normxy
c          
          ! Elasto-plastic stresses update   
          signxx(i) = signxx(i) - (a11*dpxx(i) + a12*dpyy(i))
          signyy(i) = signyy(i) - (a11*dpyy(i) + a12*dpxx(i))
          signxy(i) = signxy(i) - dpxy(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(nexp*log(cc(i)))
          yld(i) = max(yld(i),em10)
c
          ! Update Hill equivalent stress          
          sighl(i) = (ff+hh)*signyy(i)**2 + (gg+hh)*signxx(i)**2
     .             - two*hh*signxx(i)*signyy(i) + two*nn*signxy(i)**2    
          sighl(i) = sqrt(max(zero,sighl(i)))          
c
          ! Update yield function value
          phi(i)   = sighl(i) - yld(i)
c      
          ! Transverse strain update
          dpzz(i) = dpzz(i) - (dpxx(i)+dpyy(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
      DO i=1,nel
        ! Modified Mohr Criteria Failure Model
        IF (off(i) == one) THEN
          ! Pressure
          p = third*(signxx(i) + signyy(i))
          ! Von Mises equivalent stress
          svm(i) = signxx(i)**2 + signyy(i)**2 - signxx(i)*signyy(i) +
     .             three*signxy(i)**2
          svm(i) = sqrt(svm(i))
C
          ! Computation of triaxiality
          eta = p/max(svm(i),em20)
          IF (eta < -two_third) eta = -two_third
          IF (eta > two_third)  eta =  two_third
          ! Computation of Lode angle
          cos3theta = -half*twenty7*eta*(eta**2 - third)
          IF (cos3theta < -one ) cos3theta = -one
          IF (cos3theta > one )  cos3theta =  one
          theta = one - two*acos(cos3theta)/pi
C             
          ! Computation of the failure factor
          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/nexp)
C               
          ! Computation of the damage variable
          IF (inloc > 0) THEN 
            dam(i) = dam(i) + max(dplanl(i),zero)/max(em20,epsf)
          ELSE
            dam(i) = dam(i) + dpla(i)/max(em20,epsf)
          ENDIF
          IF (dam(i) >= dc) THEN 
            dam(i) = dc
            off(i) = four_over_5
          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(a11/rho0(i)) 
        ! Coefficient for hourglass
        IF (dpla(i) > zero) THEN 
          etse(i)  = h(i) / (h(i) + e)
          hardm(i) = h(i)
        ELSE
          etse(i)  = one
          hardm(i) = zero
        ENDIF
        ! Computation of the thickness variation 
        deelzz(i)  = -nu*(signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/e
        IF (inloc > 0) THEN 
          IF (loff(i) == one) THEN 
            normxx   = (gg*signxx(i) + hh*(signxx(i)-signyy(i)))/(max(sighl(i),em20))
            normyy   = (ff*signyy(i) + hh*(signyy(i)-signxx(i)))/(max(sighl(i),em20))
            dezz(i)  = deelzz(i) - max(dplanl(i),zero)*(normxx+normyy)
          ENDIF
        ELSE
          dezz(i)  = deelzz(i) + dpzz(i)
        ENDIF
        thk(i)     = thk(i) + dezz(i)*thkly(i)*off(i)  
      ENDDO     
C
      END
C