sigeps97.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/.
!||====================================================================
!||    sigeps97         ../engine/source/materials/mat/mat097/sigeps97.F
!||--- called by ------------------------------------------------------
!||    mulaw            ../engine/source/materials/mat_share/mulaw.F90
!||--- uses       -----------------------------------------------------
!||    elbufdef_mod     ../common_source/modules/mat_elem/elbufdef_mod.F90
!||    i22bufbric_mod   ../common_source/modules/interfaces/cut-cell-search_mod.F
!||    i22tri_mod       ../common_source/modules/interfaces/cut-cell-search_mod.F
!||====================================================================
      SUBROUTINE sigeps97 (
     1     NEL    ,NUPARAM ,NUVAR   ,TBURN    ,
     2     TIME   ,UPARAM  ,BFRAC   ,
     3     RHO0   ,RHO     ,EINT    ,DELTAX ,
     4     EPSPXX ,EPSPYY ,EPSPZZ  ,
     8     SIGNXX ,SIGNYY  ,SIGNZZ  ,SIGNXY  ,SIGNYZ   ,SIGNZX ,
     9     SIGVXX ,SIGVYY  ,SIGVZZ  ,SIGVXY  ,SIGVYZ   ,SIGVZX ,
     A     SOUNDSP,VISCMAX ,UVAR    ,OFF     ,
     C     GEO    ,PID     ,ILAY    ,NG      ,ELBUF_TAB,
     E     VOLN   ,
     F     QNEW   ,QOLD    ,DPDE)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE elbufdef_mod   
      USE i22bufbric_mod
      USE i22tri_mod       
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   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "param_c.inc"
#include      "com01_c.inc"
#include      "com08_c.inc"
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   I N P U T   A r g u m e n t s
C-----------------------------------------------
      INTEGER NEL, NUPARAM, NUVAR, PID(*), ILAY, NG
      my_real 
     .   TIME         ,UPARAM(NUPARAM),
     .   RHO(NEL)     ,RHO0(NEL)  ,
     .   EINT(NEL)    ,QNEW(NEL)  ,
     .   QOLD(NEL)    ,
     .   EPSPXX(NEL)  ,EPSPYY(NEL),EPSPZZ(NEL),
     .   GEO(NPROPG,*),SSP        ,
     .   VOLN(*)
      TYPE (ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
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),
     .    sigvxx(nel),sigvyy(nel),sigvzz(nel),
     .    sigvxy(nel),sigvyz(nel),sigvzx(nel),
     .    soundsp(nel),viscmax(nel)
      my_real,INTENT(INOUT) :: dpde(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), tburn(nel), bfrac(nel), deltax(nel)
C-----------------------------------------------
C   VARIABLES FOR FUNCTION INTERPOLATION 
C-----------------------------------------------
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
c      my_real
c     .     B1, B2, R1, R2, W1, D, PCJ, E0, C1, VCJ,
c     .     EADD, TBEGIN, TEND,BHE,
c     .     PSH,REACTION_RATE,REACTION_RATE2,A_MIL,M_MIL,N_MIL,ALPHA_UNIT,
c     .     DPDMU,
c     .     R1V, R2V, ER1M, ER2M, MUP1, MU, AA,BB, P0, P, VOLD, POLD,
c     .     WDR1V,WDR2V,DR1V,DR2V,ER1V,ER2V, PNEW,VOLO,FACM,ESPE,DVOL,
c     .     W1DF,EINC,QA,QB,QAL,QBL,DD     
c      INTEGER :: IBID, IBFRAC, QOPT, I      
c      my_real :: XL, MASS, DF, TB
      
      my_real
     .     d, pcj, e0, p0, vcj,c,psh,
     .     a(5),r(5),al(5),bl(5),rl(5),
     .     pnew,espe,dvol,
     .     qa,qb,qal,qbl,dd,bhe,
     .     lambda1,lambda2,lambda3,lambda4,lambda5,
     .     lambda,
     .     dldv1,dldv2,dldv3,dldv4,dldv5,
     .     dldv,
     .     erlv1,erlv2,erlv3,erlv4,erlv5,
     .     p1,p2,p3,p4,p5,
     .     rhoc2,
     .     rv1,rv2,rv3,rv4,rv5,
     .     rhoc2_1,rhoc2_2,rhoc2_3,rhoc2_4,rhoc2_5,
     .     dv, vv, vv_,
     .     ww, denom,  p, v0
     
      INTEGER :: IBFRAC, I
      my_real :: XL,  DF, TB
      
      TYPE(G_BUFEL_)  ,POINTER :: GBUF
      TYPE(L_BUFEL_)  ,POINTER :: LBUF      
      TYPE(BUF_LAY_)  ,POINTER :: BUFLY  
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------
      GBUF   => elbuf_tab(ng)%GBUF           
      lbuf   => elbuf_tab(ng)%BUFLY(ilay)%LBUF(1,1,1)
      bufly  => elbuf_tab(ng)%BUFLY(ilay)      
      
 
C-----------------------------------------------
C   S t a n d a r d   J W L   E O S 
C-----------------------------------------------  
    !IN CASE OF SWITCHING LAW5 LATER FROM LECMAT/MMAIN TO LECMUSER/MULAW
      
c      E0                 = UPARAM(1)   
c      P0                 = UPARAM(2)          
c      B1                 = UPARAM(4)   
c      B2                 = UPARAM(5)   
c      R1                 = UPARAM(6)   
c      R2                 = UPARAM(7)   
c      W1                 = UPARAM(8)   
c      D                  = UPARAM(9)   
c      PCJ                = UPARAM(10)  
c      BHE                = UPARAM(11)
c      IBFRAC             = UPARAM(12)  
c      QOPT               = UPARAM(13)         
c      EADD               = UPARAM(14)  
c      TBEGIN             = UPARAM(15)  
c      TEND               = UPARAM(16)  
c      REACTION_RATE      = UPARAM(17)  
c      A_MIL              = UPARAM(18)  
c      M_MIL              = UPARAM(19)  
c      N_MIL              = UPARAM(20)  
c      REACTION_RATE2     = UPARAM(21)  
c      ALPHA_UNIT         = UPARAM(22)  
c      PSH                = UPARAM(23) 
c                 
c      IF(DT1==ZERO)THEN
c        DO I=1,NEL
c          EINT(I) = E0*VOLN(I)
c        ENDDO
c        UVAR(1:NEL,4)    = VOLN(1:NEL)   !VOLD
c        UVAR(1:NEL,5)    = P0            !POLD
c      ENDIF
c      DO I=1,NEL     
c        !--------------------------------!
c        ! Calculation of BFRAC in [0,1]  !
c        !--------------------------------!
c        XL = DELTAX(I)  !VOLN(I)**THIRD
c        IF(BFRAC(I) < ONE) THEN
c         TB = - TBURN(I)
c         BFRAC(I) = ZERO
c         IF(IBFRAC/=1 .AND. TIME > TB) BFRAC(I) = D*(TIME-TB)*TWO_THIRD/XL 
c         IF(IBFRAC/=2) BFRAC(I)  = MAX( BFRAC(I) , BHE * (ONE - RHO0(I)/RHO(I)) )
c         IF(BFRAC(I) < EM04) THEN
c           BFRAC(I) = ZERO
c         ELSEIF(BFRAC(I) > ONE) THEN
c           BFRAC(I) = ONE
c         ENDIF
c        ENDIF
c        !--------------------------------!
c        ! EOS SOLVING                    !
c        !--------------------------------!
c        DF            = RHO0(I)/RHO(I)
c        R1V           = B1*W1/(R1*DF)
c        R2V           = B2*W1/(R2*DF)
c        WDR1V         = B1-R1V
c        WDR2V         = B2-R2V
c        DR1V          = W1*EINT(I)/MAX(EM20,VOLN(I))    !w*Eint/V = w*E/v  where v=V/V0
c        ER1V          = EXP(-R1*DF)
c        ER2V          = EXP(-R2*DF)
c        P             = P0 - PSH + WDR1V*ER1V+WDR2V*ER2V+DR1V
c        P             = MAX(ZERO - PSH, P)
c        SSP           = B1*ER1V*( (-W1/DF/R1) + R1*DF - W1)
c     .                + B2*ER2V*( (-W1/DF/R2) + R2*DF - W1)
c     .                + DR1V  +   (P + PSH)*W1
c        SSP           = SSP * DF 
c        SSP           = SQRT(ABS(SSP)/RHO0(I))
c        SSP           = MAX(SSP,D*(ONE-BFRAC(I)))
c        QA            = GEO(14,PID(I))
c        QB            = GEO(15,PID(I))
c        DD            = -EPSPXX(I)-EPSPYY(I)-EPSPZZ(I)
c        QAL           = QA*XL
c        QAL           = QAL*QAL
c        QBL           = QB*XL
c        VISCMAX(I)    = RHO(I)*(QAL*MAX(ZERO,DD) + QBL*SSP)
c        QNEW(I)       = VISCMAX(I)*MAX(ZERO,DD)
c        viscmax(I) = zero
c        !QNEW(I) = zero !QOLD(I)
c        VOLD          = UVAR(I,4)
c        POLD          = -UVAR(I,5)
c        DVOL          = HALF*(VOLN(I)-VOLD)        
c        EINC          = DVOL*(POLD-PSH-QOLD(I)-QNEW(I))
c        !EINT(I)       = EINT(I)+EINC      
c        QOLD(I)       = QNEW(I)
c        VOLO          = VOLN(I)/DF
c        ESPE          = (EINT(I)+EINC)/MAX(EM20,VOLO)
c        W1DF          = BFRAC(I)*W1/DF
c        FACM          = BFRAC(I)*(WDR1V*ER1V+WDR2V*ER2V) 
c        PNEW          = P0 - PSH + (FACM+ESPE*W1DF)/(ONE +W1DF*DVOL/MAX(EM20,VOLO))
c        PNEW          = MAX(ZERO - PSH, PNEW)*OFF(I)
c        EINC          = EINC-(PNEW + PSH)*DVOL
c       ! EINT(I)       = (EINT(I)-(PNEW + PSH)*DVOL)!/MAX(EM20,VOLN(I))
c        UVAR(I,4)     = VOLN(I)  
c        UVAR(I,5)     = PNEW  
c        !--------------------------------!
c        ! Returning values               !
c        !--------------------------------!     
c        SIGNXX(I)  = -PNEW
c        SIGNYY(I)  = -PNEW
c        SIGNZZ(I)  = -PNEW
c        SIGNXY(I)  = ZERO
c        SIGNYZ(I)  = ZERO
c        SIGNZX(I)  = ZERO
c        SIGVXX(I)  = ZERO
c        SIGVYY(I)  = ZERO
c        SIGVZZ(I)  = ZERO
c        SIGVXY(I)  = ZERO
c        SIGVYZ(I)  = ZERO
c        SIGVZX(I)  = ZERO
c       
c        SOUNDSP(I) = SSP
c      ENDDO! next I   
 
 
C-----------------------------------------------
C   J W L - B A K E R   E O S 
C-----------------------------------------------
 
      p0     = uparam(01)    
      psh    = uparam(02)    
      ibfrac = nint(uparam(03))
      d      = uparam(04)    
      pcj    = uparam(05)    
      e0     = uparam(06)    
      ww     = uparam(07)    
      c      = uparam(08)    
      a(1:5) = uparam(09:13) 
      r(1:5) = uparam(14:18) 
      al(1:5)= uparam(19:23) 
      bl(1:5)= uparam(24:28) 
      rl(1:5)= uparam(29:33) 
      bhe    = uparam(34)    
      vcj    = uparam(35)
 
      
      IF(dt1==zero)THEN
        DO i=1,nel
          eint(i) = e0*voln(i)
        ENDDO
        uvar(1:nel,4)    = voln(1:nel)   !VOLD
        uvar(1:nel,5)    = p0            !POLD
      ENDIF
 
      DO i=1,nel     
        !--------------------------------!
        ! Calculation of BFRAC in [0,1]  !
        !--------------------------------!
        xl = deltax(i)  !VOLN(I)**THIRD
        IF(bfrac(i) < one) THEN
         tb = - tburn(i)
         bfrac(i) = zero
         IF(ibfrac/=1 .AND. time > tb) bfrac(i) = d*(time-tb)*two_third/xl 
         IF(ibfrac/=2) bfrac(i)  = max( bfrac(i) , bhe * (one - rho0(i)/rho(i)) )
         IF(bfrac(i) < em04) THEN
           bfrac(i) = zero
         ELSEIF(bfrac(i) > one) THEN
           bfrac(i) = one
         ENDIF
        ENDIF
        !--------------------------------!
        ! EOS SOLVING                    !
        !--------------------------------!
        df            = rho0(i)/rho(i)
        v0            = rho(i)*voln(i) / rho0(i)
        espe          = eint(i)/max(em20,v0)
        
        ! REMINDER
        !UVAR(I,4)     = VOLN(I)  
        !UVAR(I,5)     = PNEW  
 
                
        vv            = voln(i)   / v0   ! current relative volume
        vv_           = uvar(i,4) / v0   ! previous relative volume
        dv            = vv-vv_
 
        erlv1         = exp(-rl(1)*vv)
        erlv2         = exp(-rl(2)*vv)
        erlv3         = exp(-rl(3)*vv)
        erlv4         = exp(-rl(4)*vv)
        erlv5         = exp(-rl(5)*vv)
 
        lambda1       = (al(1)*vv+bl(1))*erlv1
        lambda2       = (al(2)*vv+bl(2))*erlv2
        lambda3       = (al(3)*vv+bl(3))*erlv3
        lambda4       = (al(4)*vv+bl(4))*erlv4
        lambda5       = (al(5)*vv+bl(5))*erlv5
        
        lambda        = lambda1 + lambda2 + lambda3 + lambda4 + lambda5 + ww
 
        dpde(i)       = lambda/df
        
        dldv1         = al(1)*erlv1-(al(1)*vv+bl(1))*rl(1)*erlv1
        dldv2         = al(2)*erlv2-(al(2)*vv+bl(2))*rl(2)*erlv2
        dldv3         = al(3)*erlv3-(al(3)*vv+bl(3))*rl(3)*erlv3
        dldv4         = al(4)*erlv4-(al(4)*vv+bl(4))*rl(4)*erlv4
        dldv5         = al(5)*erlv5-(al(5)*vv+bl(5))*rl(5)*erlv5
        
        dldv          = dldv1 + dldv2 + dldv3 + dldv4 + dldv5
        
        rv1           = r(1)*vv
        rv2           = r(2)*vv
        rv3           = r(3)*vv
        rv4           = r(4)*vv
        rv5           = r(5)*vv
 
        p1            = a(1)*(one-lambda/rv1)*exp(-rv1)
        p2            = a(2)*(one-lambda/rv2)*exp(-rv2)
        p3            = a(3)*(one-lambda/rv3)*exp(-rv3)
        p4            = a(4)*(one-lambda/rv4)*exp(-rv4)
        p5            = a(5)*(one-lambda/rv5)*exp(-rv5)
        
        p             = p1+p2+p3+p4+p5 + lambda*espe/vv + c*(one-lambda/ww)*exp((-ww-one)*log(vv))
         
        rhoc2_1       = a(1)*( (vv*dldv-lambda)/r(1) + r(1)*vv*vv - lambda*vv )*exp(-rv1)
        rhoc2_2       = a(2)*( (vv*dldv-lambda)/r(2) + r(2)*vv*vv - lambda*vv )*exp(-rv2)
        rhoc2_3       = a(3)*( (vv*dldv-lambda)/r(3) + r(3)*vv*vv - lambda*vv )*exp(-rv3)
        rhoc2_4       = a(4)*( (vv*dldv-lambda)/r(4) + r(4)*vv*vv - lambda*vv )*exp(-rv4)
        rhoc2_5       = a(5)*( (vv*dldv-lambda)/r(5) + r(5)*vv*vv - lambda*vv )*exp(-rv5)
        
        !EINC          = -HALF*(UVAR(I,5)+P+PSH+PSH)*Dv
        
        rhoc2         = rhoc2_1 + rhoc2_2 + rhoc2_3 + rhoc2_4 + rhoc2_5 
        rhoc2         = rhoc2 + c*((ww+one)*(one-lambda/ww)+vv*dldv/ww)*exp(-ww*log(vv))
        rhoc2         = rhoc2 + (espe)*lambda + lambda*vv*(p+psh) - (espe)*vv*dldv         
 
        ssp           = sqrt(max(rhoc2/rho0(i),em20)) 
        qa            = geo(14,pid(i))
        qb            = geo(15,pid(i))
        dd            = -epspxx(i)-epspyy(i)-epspzz(i)
        qal           = qa*xl
        qal           = qal*qal
        qbl           = qb*xl
        viscmax(i)    = rho(i)*(qal*max(zero,dd) + qbl*ssp)
        qnew(i)       = viscmax(i)*max(zero,dd)
        !viscmax(I)    = zero
        
        denom         = (one+half*dv*lambda/vv) 
        denom         = one/denom
        
        !CURRENT PRESSURE
        pnew          = p - lambda/vv*half*(uvar(i,5)+psh)*dv
        pnew          = denom * pnew
        pnew          = (one-bfrac(i))*p0 + bfrac(i)*pnew      
        pnew          = max(-psh, pnew-psh)*off(i)
        
        !CURRENT ENERGY
        !pold          = uvar(i,5)
        dvol          = half*(voln(i)-uvar(i,4))        
        !EINC          = DVOL*(-PNEW-PSH-POLD-PSH-QOLD(I)-QNEW(I))
        eint(i)        = eint(i) - (psh+psh)*dvol
        
        !BACKUP FOR NEXT CYCLE
        qold(i)       = qnew(i)
        uvar(i,4)     = voln(i)  
        uvar(i,5)     = pnew  
        
        !--------------------------------!
        ! Returning values               !
        !--------------------------------!     
 
        signxx(i)  = -pnew
        signyy(i)  = -pnew
        signzz(i)  = -pnew
        signxy(i)  = zero
        signyz(i)  = zero
        signzx(i)  = zero
 
        sigvxx(i)  = zero
        sigvyy(i)  = zero
        sigvzz(i)  = zero
        sigvxy(i)  = zero
        sigvyz(i)  = zero
        sigvzx(i)  = zero
       
        soundsp(i) = ssp
 
      enddo! next I   
 
 
C-----------------------------------------------
      RETURN
C-----------------------------------------------
      END SUBROUTINE sigeps97
C-----------------------------------------------