sigeps93.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 Orthotropic Hill
!||====================================================================
!||    sigeps93       ../engine/source/materials/mat/mat093/sigeps93.F
!||--- called by ------------------------------------------------------
!||    mulaw          ../engine/source/materials/mat_share/mulaw.F90
!||--- calls      -----------------------------------------------------
!||    mstrain_rate   ../engine/source/materials/mat_share/mstrain_rate.F
!||    vinter         ../engine/source/tools/curve/vinter.F
!||====================================================================
       SUBROUTINE sigeps93(
     1      NEL      ,NUPARAM  ,NUVAR    ,MFUNC    ,KFUNC    , 
     2      NPF      ,TF       ,TIME     ,TIMESTEP ,UPARAM   ,     
     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      SOUNDSP  ,PLA      ,UVAR     ,RHO0     ,OFF      ,
     8      ET       ,YLD      ,SEQ      ,EPSP     ,ASRATE   ,
     9      NVARTMP  ,VARTMP   ,DPLA     )
C-----------------------------------------------
C   I M P L I C I T   T Y P E S
C-----------------------------------------------
#include "implicit_f.inc"
C----------------------------------------------------------------
C  I N P U T   A R G U M E N T S
C----------------------------------------------------------------
      INTEGER NEL,NUPARAM,NUVAR,NVARTMP
      my_real
     .      TIME,TIMESTEP,UPARAM(NUPARAM),
     .      RHO0(NEL),PLA(NEL),
     .      EPSPXX(NEL),EPSPYY(NEL),EPSPZZ(NEL),
     .      EPSPXY(NEL),EPSPYZ(NEL),EPSPZX(NEL),
     .      DEPSXX(NEL),DEPSYY(NEL),DEPSZZ(NEL),
     .      DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
     .      SIGOXX(NEL),SIGOYY(NEL),SIGOZZ(NEL),
     .      sigoxy(nel),sigoyz(nel),sigozx(nel),
     .      seq(nel),epsp(nel),asrate
 
      INTEGER, INTENT(INOUT) :: VARTMP(NEL,NVARTMP)
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),ET(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),yld(nel)
C----------------------------------------------------------------
C  VARIABLES FOR FUNCTION INTERPOLATION 
C----------------------------------------------------------------
      INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
      my_real
     .        TF(*)
C----------------------------------------------------------------
C  L O C A L  V A R I B L E S
C----------------------------------------------------------------
       INTEGER  I,II,J,NITER,ITER,NINDX,INDEX(NEL),
     .          J1,J2,ITAB,JJ(NEL),VP,IDEV,ISRATE,
     .          IAD1(NEL),IPOS1(NEL),ILEN1(NEL),
     .          iad2(nel),ipos2(nel),ilen2(nel)
       my_real
     .  d11,d12,d13,d22,d23,d33,g12,g23,g13,e11,e22,e33,
     .  ff,gg,hh,ll,mm,nn,sigy,qr1,qr2,cr1,cr2,fac,
     .  normxx,normyy,normxy,normzz,normzx,normyz,
     .  dlam,ddep,sig_dfdsig,dfdsig2,dpdt
       my_real
     .  sighl(nel),h(nel),y1(nel),y2(nel),
     .  dpxx(nel),dpyy(nel),dpzz(nel),dpxy(nel),
     .  dpzx(nel),dpyz(nel),dpla_dlam(nel),phi(nel),
     .  dydx1(nel),dydx2(nel),dphi_dlam(nel)
      my_real, 
     .   DIMENSION(:), ALLOCATABLE :: rate, yfac
C----------------------------------------------------------------
C-----------------------------------------------
C     USER VARIABLES INITIALIZATION
C-----------------------------------------------      
C  
      ! Elastic parameters
      d11  = uparam(4)  ! 1st Diag. component of plane stress elastic matrix
      d12  = uparam(5)  ! 2nd Diag. component of plane stress elastic matrix
      d13  = uparam(6)  ! Non Diag. component of plane stress elastic matrix
      d22  = uparam(7)  ! Component 13 of the compliance matrix    
      d23  = uparam(8)  ! Component 23 of the compliance matrix      
      d33  = uparam(9)  ! Component 33 of the compliance matrix    
      g12  = uparam(10) ! Shear modulus in 12  
      g13  = uparam(11) ! Shear modulus in 13   
      g23  = uparam(12) ! Shear modulus in 23
      e11  = uparam(13) ! Young modulus in 11
      e22  = uparam(14) ! Young modulus in 22
      e33  = uparam(15) ! Young modulus in 33
      ! Hill yield criterion coefficient
      ff   = uparam(19) ! First  Hill coefficient
      gg   = uparam(20) ! Second Hill coefficient
      hh   = uparam(21) ! Third  Hill coefficient
      ll   = uparam(22) ! Fourth Hill coefficient
      mm   = uparam(23) ! Fifth  Hill coefficient
      nn   = uparam(24) ! Sixth  Hill coefficient
      ! Continuous hardening yield stress
      sigy = uparam(25) ! Initial yield stress
      qr1  = uparam(26) ! Voce 1 first parameter
      cr1  = uparam(27) ! Voce 1 second parameter
      qr2  = uparam(28) ! Voce 2 first parameter
      cr2  = uparam(29) ! Voce 2 second parameter
      ! Strain-rate computation flag
      vp   = nint(uparam(30))
      ! Tabulated hardening yield stress
      itab = 0  
      IF (mfunc > 0) THEN 
        itab = 1
        ALLOCATE(rate(mfunc),yfac(mfunc))
        DO i = 1,mfunc
          rate(i) = uparam(30+i)
          yfac(i) = uparam(30+mfunc+i)
        ENDDO
      ENDIF
C      
      ! Total or deviatoric strain-rate computation 
      IF ((vp == 2) .OR. (vp == 3)) THEN
        idev   = vp-2
        israte = 1
        CALL mstrain_rate(nel    ,israte ,asrate ,epsp   ,idev   ,
     .                    epspxx ,epspyy ,epspzz ,epspxy ,epspyz ,
     .                    epspzx )
      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
        et(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
        ! Computing strain-rate (only if more than 1 strain-rate function is indicated in the input deck)
        IF ((itab == 1).AND.(mfunc > 1)) THEN
          ! Plastic strain-rate
          IF (vp == 1) epsp(i)   = uvar(i,1)
          ! looking for corresponding rate in the tables
          jj(i) = 1
          DO j = 2,mfunc-1
            IF (epsp(i) >= rate(j)) jj(i) = j
          ENDDO
        ENDIF
      ENDDO 
C
      ! Computing yield stress
      !  -> Continuous yield stress 
      IF (itab == 0) THEN
        DO i = 1,nel
          yld(i) = sigy 
     .           + qr1*(one - exp(-cr1*pla(i)))
     .           + qr2*(one - exp(-cr2*pla(i))) 
          h(i) = qr1*cr1*exp(-cr1*pla(i)) + qr2*cr2*exp(-cr2*pla(i))
        ENDDO 
      ELSE
      !  -> Tabulated yield stress 
        ! ->  No strain-rate effect
        IF (mfunc == 1) THEN
          ! Recovering position tables
          ipos1(1:nel) = vartmp(1:nel,1)
          iad1(1:nel) = npf(kfunc(1)) / 2 + 1
          ilen1(1:nel) = npf(kfunc(1)+1) / 2 - iad1(1:nel) - ipos1(1:nel)
          ! Computing the plastic strain interpolation
          CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1) 
          ! Storing the position
          vartmp(1:nel,1) = ipos1(1:nel)
          ! Computing the yield stress and its derivative
          yld(1:nel) = yfac(1)*y1(1:nel)
          h(1:nel)   = yfac(1)*dydx1(1:nel)
        ! ->  Strain-rate effect
        ELSE
          DO i=1,nel
            j1 = jj(i)
            j2 = j1 + 1           
            ! Recovering position tables for the first rate
            ipos1(i) = vartmp(i,j1)
            iad1(i) = npf(kfunc(j1)) / 2 + 1
            ilen1(i) = npf(kfunc(j1)+1) / 2 - iad1(i) - ipos1(i)            
            ! Recovering position tables for the second rate
            ipos2(i) = vartmp(i,j2)
            iad2(i) = npf(kfunc(j2)) / 2 + 1
            ilen2(i) = npf(kfunc(j2)+1) / 2 - iad2(i) - ipos2(i)
          ENDDO
          ! Computing the plastic strain interpolation
          CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
          CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,y2)
          ! Storing new positions and computing yield stress
          DO i=1,nel
            j1 = jj(i)
            j2 = j1+1
            ! Hardening yield stress
            y1(i)  = y1(i)*yfac(j1)
            y2(i)  = y2(i)*yfac(j2)
            fac    = (epsp(i) - rate(j1))/(rate(j2) - rate(j1))
            yld(i) = y1(i) + fac*(y2(i)-y1(i))
            yld(i) = max(yld(i),em20)
            ! Derivatives of yield stress
            dydx1(i) = dydx1(i)*yfac(j1)
            dydx2(i) = dydx2(i)*yfac(j2)
            h(i)     = dydx1(i)+fac*(dydx2(i)-dydx1(i))
            ! New positions
            vartmp(i,j1) = ipos1(i)
            vartmp(i,j2) = ipos2(i)
          ENDDO
        ENDIF     
      ENDIF 
c      
      !========================================================================
      ! - COMPUTATION OF TRIAL VALUES
      !========================================================================       
      DO i=1,nel      
c    
        ! Computation of the trial stress tensor        
        signxx(i) = sigoxx(i) + d11*depsxx(i) + d12*depsyy(i) + d13*depszz(i)
        signyy(i) = sigoyy(i) + d12*depsxx(i) + d22*depsyy(i) + d23*depszz(i)
        signzz(i) = sigozz(i) + d13*depsxx(i) + d23*depsyy(i) + d33*depszz(i)
        signxy(i) = sigoxy(i) + g12*depsxy(i)
        signyz(i) = sigoyz(i) + g23*depsyz(i)
        signzx(i) = sigozx(i) + g13*depszx(i)
C
        ! 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(zero,sighl(i)))
      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)-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 * (d11 * normxx + d12 * normyy + d13 * normzz )
     .            + normyy * (d12 * normxx + d22 * normyy + d23 * normzz )
     .            + normzz * (d13 * normxx + d23 * normyy + d33 * normzz )
     .            + normxy * normxy * g12
     .            + normyz * normyz * g23
     .            + normzx * normzx * g13      
c
          !   b) Derivatives with respect to plastic strain P 
          !   ------------------------------------------------  
c          
          !     i) Derivative of the yield stress with respect to plastic strain dYLD / dPLA
          !     ----------------------------------------------------------------------------
          ! Already computed 
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) - (d11*dpxx(i) + d12*dpyy(i) + d13*dpzz(i))
          signyy(i) = signyy(i) - (d12*dpxx(i) + d22*dpyy(i) + d23*dpzz(i))
          signzz(i) = signzz(i) - (d13*dpxx(i) + d23*dpyy(i) + d33*dpzz(i))
          signxy(i) = signxy(i) - dpxy(i)*g12
          signyz(i) = signyz(i) - dpyz(i)*g23
          signzx(i) = signzx(i) - dpzx(i)*g13
c          
          ! Cumulated plastic strain and strain rate update           
          ddep    = dlam*dpla_dlam(i)
          dpla(i) = max(zero, dpla(i) + ddep)
          pla(i)  = pla(i) + ddep   
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
          ! If the continuous hardening yield stress is chosen
          IF (itab == 0) THEN
            ! Update the hardening yield stress
            yld(i) = sigy 
     .             + qr1*(one - exp(-cr1*pla(i)))
     .             + qr2*(one - exp(-cr2*pla(i))) 
            h(i)   = qr1*cr1*exp(-cr1*pla(i)) + qr2*cr2*exp(-cr2*pla(i))     
c
            ! Update yield function value
            phi(i)   = sighl(i) - yld(i)
          ENDIF
c
        ENDDO
        ! End of the loop over the yielding elements
c
        ! If the tabulated yield stress is chosen
        IF (itab == 1) THEN
          ! ->  No strain-rate effect
          IF (mfunc == 1) THEN
            ! Recovering position tables
            ipos1(1:nel) = vartmp(1:nel,1)
            iad1(1:nel) = npf(kfunc(1)) / 2 + 1
            ilen1(1:nel) = npf(kfunc(1)+1) / 2 - iad1(1:nel) - ipos1(1:nel)
            ! Computing the plastic strain interpolation
            CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1) 
            ! Storing the position
            vartmp(1:nel,1) = ipos1(1:nel)
            ! Computing the yield stress and its derivative
            yld(1:nel)   =  yfac(1)*y1(1:nel)
            h(1:nel)     =  yfac(1)*dydx1(1:nel)
            ! Update the yield function 
            phi(1:nel)   = sighl(1:nel) - yld(1:nel)
          ! ->  Strain-rate effect
          ELSE
            DO i=1,nel
              j1 = jj(i)
              j2 = j1 + 1           
              ! Recovering position tables for the first rate
              ipos1(i) = vartmp(i,j1)
              iad1(i) = npf(kfunc(j1)) / 2 + 1
              ilen1(i) = npf(kfunc(j1)+1) / 2 - iad1(i) - ipos1(i)            
              ! Recovering position tables for the second rate
              ipos2(i) = vartmp(i,j2)
              iad2(i) = npf(kfunc(j2)) / 2 + 1
              ilen2(i) = npf(kfunc(j2)+1) / 2 - iad2(i) - ipos2(i)
            ENDDO
            ! Computing the plastic strain interpolation
            CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
            CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,y2)
            ! Storing new positions and computing yield stress
            DO i=1,nel
              j1 = jj(i)
              j2 = j1+1
              ! Hardening yield stress
              y1(i)  = y1(i)*yfac(j1)
              y2(i)  = y2(i)*yfac(j2)
              fac    = (epsp(i) - rate(j1))/(rate(j2) - rate(j1))
              yld(i) = y1(i) + fac*(y2(i)-y1(i))
              yld(i) = max(yld(i),em20)
              ! Derivatives of yield stress
              dydx1(i) = dydx1(i)*yfac(j1)
              dydx2(i) = dydx2(i)*yfac(j2)
              h(i)     = dydx1(i)+fac*(dydx2(i)-dydx1(i))
              ! New positions
              vartmp(i,j1) = ipos1(i)
              vartmp(i,j2) = ipos2(i)
              ! Update the yield function 
              phi(i)   = sighl(i) - yld(i)
            ENDDO
          ENDIF            
        ENDIF
c
      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
        ! Equivalent stress
        seq(i)     = sighl(i)
        ! Sound-speed
        soundsp(i) = sqrt(max(d11,d22,d33)/rho0(i))
        ! Coefficient for hourglass
        IF (dpla(i)>zero) THEN 
          et(i)    = h(i) / (h(i) + max(e11,e22,e33))
        ELSE
          et(i)  = one
        ENDIF
        ! Plastic strain-rate filtering (if needed)
        IF ((itab == 1).AND.(mfunc > 1).AND.(vp == 1)) THEN 
          dpdt      = dpla(i)/max(em20,timestep)
          uvar(i,1) = asrate * dpdt + (one - asrate) * uvar(i,1)
          epsp(i)   = uvar(i,1)
        ENDIF
      ENDDO 
      IF (ALLOCATED(rate)) DEALLOCATE(rate)
      IF (ALLOCATED(yfac)) DEALLOCATE(yfac)      
C
      END