m5law.F File Reference

#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "param_c.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	m5law (pm, sig, eint, rho, psh, p0, tburn, bfrac, voln, deltax, mat, nel, ssp, df, er1v, er2v, wdr1v, wdr2v, w1, rho0, amu, nummat, tt, dpde)

Function/Subroutine Documentation

m5law()

subroutine m5law	(	dimension(npropm,nummat), intent(in)	pm,
		dimension(nel,6), intent(inout)	sig,
		dimension(nel), intent(in)	eint,
		dimension(nel), intent(in)	rho,
		dimension(*), intent(inout)	psh,
		dimension(*), intent(inout)	p0,
		dimension(mvsiz), intent(in)	tburn,
		dimension(mvsiz), intent(inout)	bfrac,
		dimension(mvsiz), intent(in)	voln,
		dimension(nel), intent(in)	deltax,
		integer, dimension(*), intent(in)	mat,
		integer, intent(in)	nel,
		dimension(nel), intent(inout)	ssp,
		dimension(*), intent(inout)	df,
		dimension(nel), intent(inout)	er1v,
		dimension(nel), intent(inout)	er2v,
		dimension(nel), intent(inout)	wdr1v,
		dimension(nel), intent(inout)	wdr2v,
		dimension(nel), intent(inout)	w1,
		dimension(nel), intent(inout)	rho0,
		dimension(*), intent(in)	amu,
		integer, intent(in)	nummat,
		intent(in)	tt,
		dimension(nel), intent(inout)	dpde )

Parameters

[in]	nummat	size for PM array
[in]	mat	application : [1:NEL] -> mat_id
[in]	nel	number of element in the current group

Definition at line 28 of file m5law.F.

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