sigeps96.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/.
!||====================================================================
!||    sigeps96   ../engine/source/materials/mat/mat096/sigeps96.F
!||--- calls      -----------------------------------------------------
!||    vinter2    ../engine/source/tools/curve/vinter.F
!||====================================================================
      SUBROUTINE sigeps96 (
     1     NEL     ,NGL     ,NUPARAM ,NUVAR   ,NFUNC   ,IFUNC   , 
     2     NPF     ,TF      ,UPARAM  ,UVAR    ,JTHE    ,
     3     RHO     ,TEMP    ,DEFP    ,SOUNDSP ,OFF     ,EPSD    ,
     4     EPSPXX  ,EPSPYY  ,EPSPZZ  ,EPSPXY  ,EPSPYZ  ,EPSPZX  ,
     5     DEPSXX  ,DEPSYY  ,DEPSZZ  ,DEPSXY  ,DEPSYZ  ,DEPSZX  ,
     6     SIGOXX  ,SIGOYY  ,SIGOZZ  ,SIGOXY  ,SIGOYZ  ,SIGOZX  ,
     7     SIGNXX  ,SIGNYY  ,SIGNZZ  ,SIGNXY  ,SIGNYZ  ,SIGNZX  )
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
C-----------------------------------------------
C---------+---------+---+---+--------------------------------------------
C VAR     | SIZE    |TYP| RW| DEFINITION
C---------+---------+---+---+--------------------------------------------
C NEL     |  1      | I | R | SIZE OF THE ELEMENT GROUP NEL
C NUPARAM |  1      | I | R | SIZE OF THE USER PARAMETER ARRAY
C NUVAR   |  1      | I | R | NUMBER OF USER ELEMENT VARIABLES
C UPARAM  | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
C UVAR    |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
C RHO     | NEL     | F | R | DENSITY
C TEMP    | NEL     | F | R | ELEMENT TEMPERATURE (TSTAR)
C DEFP    | NEL     | F | R | PLASTIC STRAIN
C SOUNDSP | NEL     | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
C---------+---------+---+---+--------------------------------------------
C EPSPXX  | NEL     | F | R | STRAIN RATE XX
C EPSPYY  | NEL     | F | R | STRAIN RATE YY
C ...     |         |   |   |
C DEPSXX  | NEL     | F | R | STRAIN INCREMENT XX
C DEPSYY  | NEL     | F | R | STRAIN INCREMENT YY
C ...     |         |   |   |
C---------+---------+---+---+--------------------------------------------
C SIGOXX  | NEL     | F | R | OLD ELASTO PLASTIC STRESS XX
C SIGOYY  | NEL     | F | R | OLD ELASTO PLASTIC STRESS YY
C ...     |         |   |   |
C SIGNXX  | NEL     | F | W | NEW ELASTO PLASTIC STRESS XX
C SIGNYY  | NEL     | F | W | NEW ELASTO PLASTIC STRESS YY
C ...     |         |   |   |
C---------+---------+---+---+--------------------------------------------
C   I N P U T   A r g u m e n t s
C-----------------------------------------------
      INTEGER NEL,NUPARAM,NUVAR,NFUNC,JTHE
      INTEGER NPF(*),NGL(NEL),IFUNC(NFUNC)
      my_real TF(*),UPARAM(NUPARAM)
      my_real,DIMENSION(NEL), INTENT(IN) :: RHO,TEMP,OFF,
     .   EPSPXX,EPSPYY,EPSPZZ,EPSPXY,EPSPYZ,EPSPZX,
     .   DEPSXX,DEPSYY,DEPSZZ,DEPSXY,DEPSYZ,DEPSZX,
     .   SIGOXX,SIGOYY,SIGOZZ,SIGOXY,SIGOYZ,SIGOZX
C-----------------------------------------------
C   O U T P U T   A r g u m e n t s
C-----------------------------------------------
       my_real ,DIMENSION(NEL), INTENT(OUT) :: soundsp,
     .    signxx,signyy,signzz,signxy,signyz,signzx
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 ,DIMENSION(NEL,NUVAR), INTENT(INOUT), TARGET :: uvar
       my_real ,DIMENSION(NEL,2), INTENT(INOUT), TARGET :: defp
       my_real ,DIMENSION(NEL), INTENT(INOUT) :: epsd
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,ITER,NITER,IVARF,IFUNC_E,IFUNC_YLD
      INTEGER, DIMENSION(NEL) :: IPOS1,ILEN1,IAD1,IPOS2,ILEN2,IAD2
c
      my_real tanb,tanp,eini,kini,gini,nu,yxi,kep,h,ho,hm,res,dres,
     .  sigy,siga,sigb,sigr,ra1,ra2,rb,rc,na,fac,q2,pp1,pp2,pla2,pla,dpla,
     .  dfdp,dfdc,dgdp,denom,ldot,dwpx2,qx2,depslv,depspv,depspd,
     .  e1,e2,e3,e4,e5,e6,tr,sy1,sy2,sy3,sy4,sy5,sy6,
     .  dfds1,dfds2,dfds3,dfds4,dfds5,dfds6,depsp1,depsp2,depsp3,depsp4,
     .  depsp5,depsp6,depsl1,depsl2,depsl3,depsl4,depsl5,depsl6,epsp,
     .  qy,fy,py,f1,f2,x,x1,x2,dx,alpha
      my_real, DIMENSION(NEL) :: k,g,g2,ra,depsv,dav,d1,d2,d3,
     .  ds1,ds2,ds3,ds4,ds5,ds6,dp,de1,de2,de3,de4,de5,de6,dydx,
     .  fo,po,so1,so2,so3,so4,so5,so6,s1,s2,s3,s4,s5,s6,p,yld,
     .  s1n,s2n,s3n,s4n,s5n,s6n,pn,qn,q,f,fac_e,frate   
      my_real small,big,tol 
      my_real ,DIMENSION(:), POINTER :: epspd,epspv
C=======================================================================
c      UPARAM(1)   E    : Young modulus
c      UPARAM(2)   G    : Shear modulus
c      UPARAM(3)   K    : Bulk modulus
c      UPARAM(4)   NU   : Poisson coefficient
c      UPARAM(5)   TANB : Tangent beta (yield friction angle)
c      UPARAM(6)   TANP : Tangent Psi  (plastic flow angle)
c      UPARAM(7)   YXI  : sqrt(1 + TANB^2)
c      UPARAM(8)   SIGY : initial yield value in pure shear (p=0)
c      UPARAM(9)   SIGA : Max yield gain during early hardening peak
c      UPARAM(10)  SIGB : Max Yield drop during softening phase
c      UPARAM(11)  SIGR : Rubbery modulus (hardening modulus in rubbery state)
c      UPARAM(12)  RA1  : 1st exponent coefficient for SIGA evolution
c      UPARAM(13)  RA2  : 2nd exponent coefficient for SIGA evolution
c      UPARAM(14)  NA   : exponent coefficient for temperature dependency of SIGA
c      UPARAM(15)  RB   : exponent coefficient for SIGB evolution
c      UPARAM(16)  RC   : coefficient for SIGR evolution
c      UPARAM(20)  NUP  : Plastic deformation Poisson coefficient
c---------------------     
c      UVAR(1)     YIELD
c      UVAR(2)     EPSD
C=======================================================================
      niter = 3
      small = em10
      big   = ep20
      tol   = em20
c---------------------     
      epspd => defp(1:nel,1)   ! Von Mises Equivalent Plastic Strain  
      epspv => defp(1:nel,2)   ! Volumetric Plastic Strain
c---------------------     
      eini = uparam(1)
      gini = uparam(2)
      kini = uparam(3)
      nu   = uparam(4)
      tanb = uparam(5)
      tanp = uparam(6)
      yxi  = uparam(7)
      sigy = uparam(8) 
      siga = uparam(9) 
      sigb = uparam(10)
      sigr = uparam(11)
      ra1  = uparam(12)
      ra2  = uparam(13)
      na   = uparam(14)
      rb   = uparam(15)
      rc   = uparam(16)
      kep  = uparam(17)    ! 1 / (1 + nup^2)   : d_epsp = kep*ldot
      alpha= uparam(19)    ! strain rate smoothing coefficient
c
      ifunc_e   = ifunc(1)
      ifunc_yld = ifunc(2)
c-------------------------------------------
c     TEMPERATURE AND STRAIN RATE DEPENDENCY (YOUNG MODULUS, YIELD)
c-------------------------------------------
      ivarf = 2
      IF (jthe == 1. and. ifunc_e > 0) THEN
        DO i=1,nel
          ipos1(i) = nint(uvar(i, ivarf + 1))
          iad1(i)  = npf(ifunc_e) / 2 + 1
          ilen1(i) = npf(ifunc_e + 1) / 2 - iad1(i) - ipos1(i)
        ENDDO
c
        CALL vinter2(tf,iad1,ipos1,ilen1,nel,temp,dydx,fac_e)
        uvar(1:nel,ivarf + 1) = ipos1(1:nel)
c
        k(1:nel)  = kini*fac_e(1:nel)
        g(1:nel)  = gini*fac_e(1:nel)
        g2(1:nel) = g(1:nel)*two
        ivarf = ivarf + 1
      ELSE
        k(1:nel)  = kini
        g(1:nel)  = gini
        g2(1:nel) = gini*two
      ENDIF
c---
      IF (ifunc_yld > 0) THEN
        DO i=1,nel
          tr = (epspxx(i) + epspyy(i) + epspzz(i))*third
          e1 = epspxx(i) - tr
          e2 = epspyy(i) - tr
          e3 = epspzz(i) - tr
          e4 = half*epspxy(i)
          e5 = half*epspyz(i)
          e6 = half*epspzx(i)
          epsp = half*(e1**2 + e2**2 + e3**2) + e4**2 + e5**2 + e6**2
          epsp = sqrt(three*epsp)/three_half
          epsp = alpha*epsp + (one -alpha)*epsd(i)
          epsd(i) = epsp
        ENDDO
c
        DO i=1,nel
          ipos2(i) = nint(uvar(i, ivarf + 1))
          iad2(i)  = npf(ifunc_yld) / 2 + 1
          ilen2(i) = npf(ifunc_yld + 1) / 2 - iad2(i) - ipos2(i)
        ENDDO
c
        CALL vinter2(tf,iad2,ipos2,ilen2,nel,epsd,dydx,frate)
        uvar(1:nel,ivarf + 1) = ipos2(1:nel)
      ELSE
         frate(1:nel) = one
      ENDIF
      uvar(1:nel,2) = epsd(1:nel)
c---
      IF (jthe == 1) THEN
        ra(1:nel) = ra1*temp(1:nel)**na + ra2
      ELSE
        ra(1:nel) = ra1 + ra2
      ENDIF
c------------------------------------------
c     ELASTIC TRIAL INCREMENT
c------------------------------------------
      DO i=1,nel
        pla  = epspd(i)
        pla2 = pla**2
        pp1  = one - rc*pla2
        pp2  = pp1 + two
        yld(i) = sigy*frate(i) + siga*(one - exp(-ra(i)*pla)) 
     .         - sigb *(one - exp(-rb*pla))
     .         + sigr*pla*third*pp2 / pp1
c
        depsv(i) = (depsxx(i) + depsyy(i) + depszz(i))
        dav(i) = depsv(i)*third
        d1(i)  = depsxx(i)-dav(i)
        d2(i)  = depsyy(i)-dav(i)
        d3(i)  = depszz(i)-dav(i)
c       Increments elastiques pression et deviateur        
        dp(i) =-k(i)* depsv(i)
        ds1(i)= g2(i)* d1(i) 
        ds2(i)= g2(i)* d2(i)
        ds3(i)= g2(i)* d3(i)
        ds4(i)= g(i) * depsxy(i)
        ds5(i)= g(i) * depsyz(i)
        ds6(i)= g(i) * depszx(i)
c       Valeurs elastiques pression et deviateur   
        po(i)  =-(sigoxx(i)+sigoyy(i)+sigozz(i))*third
        so1(i) = sigoxx(i) + po(i)
        so2(i) = sigoyy(i) + po(i)
        so3(i) = sigozz(i) + po(i)
        so4(i) = sigoxy(i)
        so5(i) = sigoyz(i)
        so6(i) = sigozx(i)
 
        s1(i) = so1(i) + ds1(i)
        s2(i) = so2(i) + ds2(i)
        s3(i) = so3(i) + ds3(i)
        s4(i) = so4(i) + ds4(i)
        s5(i) = so5(i) + ds5(i)
        s6(i) = so6(i) + ds6(i)
        p(i)  = po(i)  + dp(i)
c       Contraintes elastiques
        signxx(i) = s1(i) - p(i)
        signyy(i) = s2(i) - p(i)
        signzz(i) = s3(i) - p(i)
        signxy(i) = s4(i)
        signyz(i) = s5(i)
        signzx(i) = s6(i)
      ENDDO
c------------------------------------------
c     EVALUATION FONCTION CRITERE DE PLASTICITE
c------------------------------------------
c     old criterion
      DO i=1,nel          
        q2   = three_half*(so1(i)**2+so2(i)**2+so3(i)**2)
     .       + three*(so4(i)**2+so5(i)**2+so6(i)**2)
        q(i) = sqrt(q2)
        fo(i)= q(i) - max(zero, po(i)*tanb + yld(i))
      ENDDO
c------------------------------------------
      DO i=1,nel          
        q2  = three_half*(s1(i)**2+s2(i)**2+s3(i)**2)
     .      + three*(s4(i)**2+s5(i)**2+s6(i)**2)
        q(i)= sqrt(q2)
        f(i)= q(i) - max(zero, p(i)*tanb + yld(i))
      ENDDO
c------------------------------------------
c     ITERATION POUR TROUVER LE PT D'INTERSECTION AVEC LE CRITERE
c------------------------------------------
      DO i=1,nel
        IF (f(i) > tol) THEN
c
          IF (fo(i) < zero) THEN
            x1 = one
            x2 = zero
            f1 = fo(i)
            f2 = f(i)
            DO iter=1,niter
              x   = (x1*f2 - x2*f1) / (f2 - f1)  ! regula falsi
              sy1 = s1(i) - x * ds1(i)
              sy2 = s2(i) - x * ds2(i)    
              sy3 = s3(i) - x * ds3(i)
              sy4 = s4(i) - x * ds4(i)
              sy5 = s5(i) - x * ds5(i)
              sy6 = s6(i) - x * ds6(i)
              py  = p(i)  - x * dp(i) 
              qy  = three_half*(sy1**2+sy2**2+sy3**2) + three*(sy4**2+sy5**2+sy6**2)
              qy  = sqrt(qy)
              fy  = qy - max(zero, py*tanb + yld(i))
c
              IF (fy*f1 > 0) THEN
                x1 = x
                f1 = fy
              ELSEIF (fy*f2 > 0) THEN
                x2 = x
                f2 = fy
              ELSE
                EXIT
              ENDIF
            ENDDO
c
c           constraints
            ds1(i) = x*ds1(i)
            ds2(i) = x*ds2(i)
            ds3(i) = x*ds3(i)
            ds4(i) = x*ds4(i)
            ds5(i) = x*ds5(i)
            ds6(i) = x*ds6(i)
            dp(i)  = x*dp(i)
c
            s1(i) = sy1
            s2(i) = sy2
            s3(i) = sy3
            s4(i) = sy4
            s5(i) = sy5
            s6(i) = sy6
            p(i)  = py
            q(i)  = qy
c           deformations
            de1(i) = x*d1(i)
            de2(i) = x*d2(i)
            de3(i) = x*d3(i)
            de4(i) = x*depsxy(i)
            de5(i) = x*depsyz(i)
            de6(i) = x*depszx(i)
            depsv(i)=x*depsv(i)
          ELSE
            de1(i) = d1(i)
            de2(i) = d2(i)
            de3(i) = d3(i)
            de4(i) = depsxy(i)
            de5(i) = depsyz(i)
            de6(i) = depszx(i)
          ENDIF
        ENDIF
      ENDDO
c---------------------------------
c     PLASTIC FLOW  : FULL NEWTON ITERATIONS
c---------------------------------
      DO i=1,nel
        IF (f(i) > zero) THEN
          dfdp =-tanb  
          dgdp =-tanp  
          dfdc =-one
c
          IF (q(i) > em20) THEN
            fac   = three_half/q(i)
            dfds1 = s1(i)*fac
            dfds2 = s2(i)*fac
            dfds3 = s3(i)*fac
            dfds4 = s4(i)*fac
            dfds5 = s5(i)*fac
            dfds6 = s6(i)*fac
          ELSE
            dfds1 =zero
            dfds2 =zero
            dfds3 =zero
            dfds4 =zero
            dfds5 =zero
            dfds6 =zero        
          ENDIF
c
          ho = siga*ra(i)*exp(-ra(i)*epspd(i)) - sigb*rb*exp(-rb*epspd(i))
     .       + sigr*(one + two_third*rc*pla2 * pp2 / pp1**2)
          denom = three*g(i) + k(i)*dfdp*dgdp - dfdc*ho
          IF (denom /= zero) THEN
            ldot = f(i) / denom
          ELSE
            ldot = zero
          ENDIF
c-----------------------------------------------------------------------
          niter = 4
          DO iter=1,niter
c          
            depspd = ldot
            depspv = ldot*dgdp*third
c
            depsp1 = ldot*dfds1 
            depsp2 = ldot*dfds2 
            depsp3 = ldot*dfds3 
            depsp4 = ldot*dfds4 * two
            depsp5 = ldot*dfds5 * two
            depsp6 = ldot*dfds6 * two
c
c           Deviatoric elastic deformations
c
            depsl1 = de1(i) - depsp1
            depsl2 = de2(i) - depsp2
            depsl3 = de3(i) - depsp3
            depsl4 = de4(i) - depsp4
            depsl5 = de5(i) - depsp5
            depsl6 = de6(i) - depsp6
            depslv = depsv(i) + depspv
c
c           Stress increment starting from yield surface
c
            ds1(i)= g2(i)*depsl1
            ds2(i)= g2(i)*depsl2
            ds3(i)= g2(i)*depsl3
            ds4(i)= g(i) *depsl4
            ds5(i)= g(i) *depsl5
            ds6(i)= g(i) *depsl6
            dp(i) =-k(i) *depslv
c
c           New stress
c
            s1n(i) = s1(i) + ds1(i)
            s2n(i) = s2(i) + ds2(i)
            s3n(i) = s3(i) + ds3(i)
            s4n(i) = s4(i) + ds4(i)
            s5n(i) = s5(i) + ds5(i)
            s6n(i) = s6(i) + ds6(i)
            pn(i)  = p(i)  + dp(i)
c
c           Update yield and plastic increment
c           
            q2   = three_half*(s1n(i)**2 + s2n(i)**2 + s3n(i)**2)
     .           + three*(s4n(i)**2 + s5n(i)**2 + s6n(i)**2)
            qn(i)= sqrt(q2)
c
            dpla = ldot*kep
            pla  = epspd(i) + dpla
            pla2 = pla**2
            pp1  = one - rc*pla2
            pp2  = pp1 + two
            yld(i) = sigy*frate(i) + siga*(one - exp(-ra(i)*pla))
     .             - sigb *(one - exp(-rb*pla))
     .             + sigr*pla*third*pp2 / pp1
            yld(i) = max(yld(i), zero)
 
c           H = SIGA*RA(I)*EXP(-RA(I)*EPSPD(I)) - SIGB*RB*EXP(-RB*EPSPD(I))
c    .        + SIGR*(ONE + TWO_THIRD*RC*PLA2 * PP2 / PP1**2)
c            HM = (HO + H) * HALF
            hm = ho
 
            res  = qn(i) - max(zero, pn(i)*tanb + yld(i))
            dres = -(three*g(i) + k(i)*dfdp*dgdp - dfdc*hm)
            
            IF (dres /= zero) THEN
              ldot = ldot - res / dres
            ENDIF
c
          ENDDO !  ITER=1,NITER
c-----------------------------------------------------------------------
c         CHECK REPROJECTION
c
          fy = qn(i) - max(zero, pn(i)*tanb + yld(i))
c
          IF (fy > small) THEN
            IF (tanb*pn(i) + yld(i) <= zero) THEN
              pn(i)  =-yld(i)/tanb
              s1n(i) = zero
              s2n(i) = zero
              s3n(i) = zero
              s4n(i) = zero
              s5n(i) = zero
              s6n(i) = zero
c              WRITE(6,*)'TRI-TRACTION FAILURE'
            ELSE 
              IF (qn(i) > small) THEN
                x = (pn(i)*tanb + yld(i)) / qn(i)
                IF (x < one-em02 .OR. x > one) THEN
                  WRITE(7,*)'REPROJ Q',x,fy,pn(i),qn(i)
                  WRITE(6,*)'REPROJ Q',x,fy,pn(i),qn(i)
                ENDIF
                s1n(i) = x*s1n(i)
                s2n(i) = x*s2n(i)
                s3n(i) = x*s3n(i)
                s4n(i) = x*s4n(i)
                s5n(i) = x*s5n(i)
                s6n(i) = x*s6n(i)
              ENDIF
            ENDIF
          ENDIF
          epspd(i) = epspd(i) + ldot*kep
          epspv(i) = epspv(i) + ldot*dgdp
c
c         Final stress update
c
          signxx(i) = s1n(i) - pn(i)
          signyy(i) = s2n(i) - pn(i)
          signzz(i) = s3n(i) - pn(i)
          signxy(i) = s4n(i)
          signyz(i) = s5n(i)
          signzx(i) = s6n(i)
c
        ENDIF   ! F > 0       
c
        uvar(i,1)= yld(i)
c
      ENDDO     ! I=1,NEL
c----------------------------------------------
c     END OF PLASTICITY PROJECTION LOOP
c----------------------------------------------
      DO i=1,nel
        soundsp(i) = sqrt((k(i) + four_over_3*g(i))/rho(i))
      ENDDO
c-----------
      RETURN
      END