m5law.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/.
!||====================================================================
!||    m5law   ../engine/source/materials/mat/mat005/m5law.F
!||--- called by ------------------------------------------------------
!||    mmain   ../engine/source/materials/mat_share/mmain.F90
!||====================================================================
      SUBROUTINE m5law(PM   ,SIG   ,EINT  ,RHO    ,PSH   ,
     1                 P0   ,TBURN ,BFRAC ,VOLN   ,DELTAX,
     2                 MAT  ,NEL   ,SSP   ,DF     ,
     3                 ER1V ,ER2V  ,WDR1V ,WDR2V  ,W1    ,
     4                 RHO0 ,AMU   ,NUMMAT,TT, DPDE)
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"
#include      "param_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER,INTENT(IN)    :: NUMMAT                !< size for PM array
      INTEGER,INTENT(IN)    :: MAT(*)                !< application : [1:NEL] -> mat_id
      INTEGER,INTENT(IN)    :: NEL                   !< number of element in the current group
      my_real,INTENT(IN)    :: TT                    !< current time
      my_real,INTENT(IN)    :: pm(npropm,nummat)     !< material buffer (real)
      my_real,INTENT(INOUT) :: sig(nel,6)            !< stress tensor
      my_real,INTENT(IN)    :: eint(nel)             !< internal energy
      my_real,INTENT(IN)    :: rho(nel)              ! mass density
      my_real,INTENT(INOUT) :: psh(*)                !< pressure shift
      my_real,INTENT(INOUT) :: p0(*)                 !< initial pressure
      my_real,INTENT(IN)    :: tburn(mvsiz)          !<time of burn (detonation time)
      my_real,INTENT(INOUT) :: bfrac(mvsiz)          !< burn fraction
      my_real,INTENT(IN)    :: voln(mvsiz)           !< element volume
      my_real,INTENT(IN)    :: deltax(nel)           !< mesh size
      my_real,INTENT(INOUT) :: ssp(nel)              !< sound speed
      my_real,INTENT(INOUT) :: df(*)                 !< V/V0 (or rho0/rho or 1/(1+mu))
      my_real,INTENT(IN)    :: amu(*)                !< volumetric strain
      my_real,INTENT(INOUT) :: er1v(nel), er2v(nel), wdr1v(nel), wdr2v(nel), w1(nel), rho0(nel)   !< working arrays to save computation time
      my_real,INTENT(INOUT) :: dpde(nel)              !< partial derivative
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,MX,IBFRAC
      my_real A,B, R1, R2, W, VDET, PSH_PARAM, BULK, P0_PARAM, RHO0_PARAM !< eos parameters
      my_real TB                                                                   !< time of burn
      my_real BFRAC1, BFRAC2, BHE                                                  !< working arrays for burn fraction calculation
      my_real r1v(mvsiz), r2v(mvsiz), dr1v(mvsiz)                                  !< temporary arrays
      my_real p(mvsiz)                                                             !< pressure
C-----------------------------------------------
C   B o d y
C-----------------------------------------------
 
      ! User defined material properties
      mx         = mat(1)
      rho0_param = pm( 1,mx)
      a          = pm(33,mx)
      b          = pm(34,mx)
      r1         = pm(35,mx)
      r2         = pm(36,mx)
      w          = pm(45,mx)
      vdet       = pm(38,mx)
      bhe        = pm(40,mx)
      psh_param  = pm(88,mx)
      p0_param   = pm(31,mx)
      bulk       = pm(44,mx)
      ibfrac     = nint(pm(41,mx))
 
      ! Initialize the material properties used later from mjwl.F
      DO i=1,nel
        rho0(i) = rho0_param
        psh(i) = psh_param
        p0(i) = p0_param
        w1(i) = w
      ENDDO
 
      ! Relative Volume Calculation
      DO i=1,nel
        df(i) = rho0(i)/rho(i) ! DF = v = V/V0 = RHO0/RHO = 1/(MU+1)
      ENDDO
 
      ! Burn Fraction Calculation
      !   BFRAC1 : time control
      !   BFRAC2 : volumetric control
      DO i=1,nel
        IF(bfrac(i) < one) THEN
          tb=-tburn(i)
          bfrac1 = zero
          bfrac2 = zero
          IF(ibfrac /= 1 .AND. tt > tb) bfrac1=vdet*(tt-tb)/three_half/deltax(i)   !time control
          IF(ibfrac /= 2)bfrac2=bhe*(one-df(i))                                    !volumetric control
          bfrac(i) = max(bfrac1,bfrac2)
          IF(bfrac(i)<em04) bfrac(i)=zero
          IF(bfrac(i)>one) THEN
            bfrac(i) = one
          ENDIF
        ENDIF
      ENDDO
 
      ! working arrays to save CPU operations. ER1V ER2V, WDR1V, WDR2V used later from mjwl.F
      DO i=1,nel
        r1v(i)   = a*w/(r1*df(i))
        r2v(i)   = b*w/(r2*df(i))
        wdr1v(i) = a-r1v(i)
        wdr2v(i) = b-r2v(i)
        dr1v(i)  = w*eint(i)/max(em20,voln(i))    !w*Eint/V = w*E/v  where v=V/V0
        er1v(i)  = exp(-r1*df(i))
        er2v(i)  = exp(-r2*df(i))
      ENDDO
 
      ! Pressure Calculation
      IF (bulk == zero) THEN
         ! Behavior of unreacted explosive is not provided
         DO i=1,nel
            p(i) =  p0(i) + (wdr1v(i)*er1v(i)+wdr2v(i)*er2v(i)+dr1v(i))
         ENDDO
      ELSE
         ! Behavior of unreacted explosive is provided
         DO i=1,nel
            p(i) = (one - bfrac(i))*(p0(i)+bulk*amu(i)) + bfrac(i)*(wdr1v(i)*er1v(i)+wdr2v(i)*er2v(i)+dr1v(i))
         ENDDO
      ENDIF
      DO i=1,nel
            p(i) = max(zero, p(i)) - psh(i)  !PMIN = 0.0 (fluid)
      ENDDO
 
      !partial derivative at constant volume
      DO i=1,nel
        dpde(i) = w/df(i)
      ENDDO
 
      ! sound speed calculation
      DO i=1,nel
        ssp(i) = a*er1v(i)*( (-w/df(i)/r1) + r1*df(i) - w)
     .         + b*er2v(i)*( (-w/df(i)/r2) + r2*df(i) - w)
     .         + dr1v(i)  +   (p(i) + psh(i))*w
        ssp(i) = ssp(i) * df(i)
      ENDDO
      IF (bulk == zero) THEN
         DO i=1,nel
            ssp(i) = sqrt(abs(ssp(i))/rho0(i))
            ssp(i) = max(ssp(i),vdet*(one-bfrac(i)))
         ENDDO
      ELSE
         DO i=1,nel
            ssp(i) =  bfrac(i) * (ssp(i) / rho0(i)) + (one - bfrac(i)) * (bulk / rho0(i))
            ssp(i) = sqrt(abs(ssp(i)))
         ENDDO
      ENDIF
 
      ! Return the updated stress tensor. FLUID => NO DEVIATOR STRESS
      DO i=1,nel
        sig(i,1) = zero
        sig(i,2) = zero
        sig(i,3) = zero
        sig(i,4) = zero
        sig(i,5) = zero
        sig(i,6) = zero
      ENDDO
 
      RETURN
      END