SANDIA REPORT

SAND2019-6025

Printed May 2019

The Mie-Grüneisen Power Equation of State

Allen C. Robinson

Prepared by
Sandia National Laboratories
Albuquerque, New Mexico 87185
Livermore, California 94550
Issued by Sandia National Laboratories, operated for the United States Department of Energy by National
Technology & Engineering Solutions of Sandia, LLC.
NOTICE: This report was prepared as an account of work sponsored by an agency of the United States Government.
Neither the United States Government, nor any agency thereof, nor any of their employees, nor any of their
contractors, subcontractors, or their employees, make any warranty, express or implied, or assume any legal liability
or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process
disclosed, or represent that its use would not infringe privately owned rights. Reference herein to any specific
commercial product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily
constitute or imply its endorsement, recommendation, or favoring by the United States Government, any agency
thereof, or any of their contractors or subcontractors. The views and opinions expressed herein do not necessarily
state or reflect those of the United States Government, any agency thereof, or any of their contractors.
Printed in the United States of America. This report has been reproduced directly from the best available copy.
Available to DOE and DOE contractors from
U.S. Department of Energy
Office of Scientific and Technical Information
P.O. Box 62
Oak Ridge, TN 37831

Telephone: (865) 576-8401

Facsimile: (865) 576-5728
E-Mail: [email protected]
Online ordering: http://www.osti.gov/scitech
Available to the public from
U.S. Department of Commerce
National Technical Information Service
5301 Shawnee Road
Alexandria, VA 22312

Telephone: (800) 553-6847

Facsimile: (703) 605-6900
E-Mail: [email protected]
Online order: https://classic.ntis.gov/help/order-methods

NT OF E
ME N
RT
ER
A
DEP

GY

• •
IC A
U NIT

ER
ED

ST M
A TES OF A

2
ABSTRACT
We describe details of a general Mie-Grüneisen equation of state and its numerical
implementation. The equation of state contains a polynomial Hugoniot reference curve, an
isentropic expansion and a tension cutoff.

3
ACKNOWLEDGMENT

The author thanks Ann E. Mattsson of Los Alamos National Laboratory for her encouragement to
document this model. Steven W. Bova and John H. Carpenter provided detailed reviews that led
to improvements in this document and the model implementation. John H. Carpenter suggested
the backward recurrence approach for the lower incomplete gamma function.

4
CONTENTS

Terms and Definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

1. The Mie-Grüneisen Power Equation of State 7

1.1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2. Compression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3. Expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.4. Pressure Cutoff in Tension . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.5. Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

References 13

5
TERMS AND DEFINITIONS

ρ0 Reference state density

η Volumetric Strain
T0 Reference state temperature
Cv Heat capacity at constant volume (assumed constant)
Γ Grüneisen parameter
α Γρ0 /(1 − η) = Γρ = Γ0 ρ0 (assumed constant)
C0 Reference sound speed
K0 = ρ0C02 Reference adiabatic bulk modulus
Kn Hugoniot coefficient
Pmin Minimum pressure on expansion isentrope

Us Shock velocity in Us − u p shock relationship

u p Particle velocity in Us − u p shock relationship
s Linear coefficient in Us − u p shock relationship

ρ Density
P Pressure
E Specific internal energy
T Temperature
S Entropy
p
Cs Sound speed, (∂ P/∂ ρ)S

6
1. THE MIE-GRÜNEISEN POWER
EQUATION OF STATE

1.1. INTRODUCTION

The classical Mie-Grüneisen equation of state assumes that the pressure is non-linearly related to
the density but linear in the specific energy relative to a reference curve, R. Thus

P(η, E) = PR (η) + α [E − ER (η)] (1.1)

where the volumetric strain η is

ρ0 v
η = 1− = 1− . (1.2)
ρ v0
ρ is the density and v is the specific volume. The subscript 0 refers to a particular point on the
reference curve where the strain is zero and the pressure is zero. This point is called the reference
state. In addition, we assume that α = Γρ0 /(1 − η) = Γρ = Γ0 ρ0 is constant. Γ is the Grüneisen
parameter. The heat capacity Cv = (∂ E/∂ T )v is also assumed constant giving

E(η, T ) = ER (η) +Cv [T − TR (η)] (1.3)

Historically, the Mie-Grüneisen equation of state has been primarily used with application to
compression of metal solids [2] and the compressive reference curve is assumed to be a Hugoniot.
However, the reference curve can be a different curve. For example, a Mie-Grüneisen equation of
state using a reference isentrope was matched to a Mie-Grüneisen Us −U p shock Hugoniot with
two parameters for the purpose of demonstrating relevant two dimensional isentropic jet flows
[3].
The sound speed is computed by differentiating Equation 1.1 with respect to density at constant
entropy and utilizing the thermodynamic relationship

dE = T dS − Pdv = T dS + Pv2 dρ. (1.4)

Thus
    
∂P dη ∂E dη
Cs2 = = PR0 (η) +α − ER0 (η)
∂ρ S dρ S∂ρ dρ
= v ρ0 PR0 (η) + α P − ρ0 ER0 (η)
2


(1.5)

7
since dη/dρ = ρ0 v2 . In addition,
 
∂P
= v2 ρ0 PR0 (η) − αCv TR0 (η)
 
(1.6)
∂ρ T

We now derive a useful and well known thermodynamics relation for integrating temperature
along the references curves [4, Equation 23]. First,
1 P
dS = dE + dv (1.7)
T T
   
1 ∂E 1 ∂E
= dT + + P dv (1.8)
T ∂T v T ∂v T
   
∂S ∂S
= dT + dv (1.9)
∂T v ∂v T
The equality of mixed differentials leads to
   
∂E ∂P
+P = T (1.10)
∂v T ∂T v
and finally  
∂P
T dS = Cv dT + T dv (1.11)
∂T v

For the Mie-Grüneisen equation of state with constant Cv and constant α = ρΓ = ρ0 Γ0 , we then
derive
T dS = Cv dT − Γ0Cv T dη (1.12)
Equation 1.12 is the key thermodynamic relationship used below to compute the temperature
along the reference curve.
The general Mie-Grüneisen power series equation of state described here contains polynomial
flexibility in the shape of the reference curves. The subscript R refers to a reference state curve
which is either a Hugoniot (subscript H) or an isentrope (subscript I). This particular model
provides three regions: a compressive region, a tension region, and a tensile pressure cutoff region
as described below.

1.2. COMPRESSION

For η > 0 we define:

PR = PH = K0 η 1 + K1 η + K2 η 2 + K3 η 3 + · · · + KM η M


(1.13)

where M is an integer giving the maximum number of terms in the polynomial and the subscript
H indicates that the reference curve is a Hugoniot. The first term, K0 = ρ0C02 , defines the
adiabatic bulk modulus at the reference point.

8
The energy on the Hugoniot curve is given by the well-known shock jump conditions [1]
PH η
ER = EH = + E0 . (1.14)
2ρ0
In order to compute the temperature on the Hugoniot, we use the method of Walsh and Christian
valid for constant α = ρΓ and constant Cv [4]. This method uses Equationn 1.14 and Equation
1.12 to obtain the temperature on the Hugoniot curve,
Substituting Equation 1.12 on the left hand side of Equation 1.7 and the shock jump condition
Equation 1.14 on the right hand side, we obtain a differential equation for the temperature along
the Hugoniot curve,
v0 η 2 d PH
 
dTH dEH
Cv − Γ0Cv TH = − v0 PH = . (1.15)
dη dη 2 dη η
Integration on the Hugoniot curve gives
 
eΓ0 η d PH
Z η
Γ0 η −Γ0 z 2
TH = T0 e + e z dz (1.16)
2Cv ρ0 0 dz z
which leads to

Γ0 η eΓ0 η
TH = T0 e + K0 (K1 I2 + 2K2 I3 + 3K3 I4 + · · · + MKM IM+1 ) (1.17)
2Cv ρ0
where
Z η
In = e−Γ0 z zn dz (1.18)
0
 n+1 Z Γ η
1 0
= e−z zn dz (1.19)
Γ0 0
 n+1
1
= γ(n + 1, Γ0 η) (1.20)
Γ0
and Z x
γ(n, x) = e−z zn−1 dz (1.21)
0
is the lower incomplete gamma function. The η derivative of TH is
dTH K0
K1 η 2 + 2K2 η 3 + 3K3 η 4 + · · · + MKM η M+1


= Γ0 TH + (1.22)
dη 2Cv ρ0

We need to evaluate γ(n, x) for arguments, x ≥ 0, and integral n ≥ 3. As a first cut one might think
to use the well known formula
!
n−1 k
x
γ(n, x) = (n − 1)! 1 − e−x ∑ (1.23)
k=0 k!
= (n − 1)! 1 − e−x en−1 (x)


(1.24)

9
where en (x) is the exponential sum function [5]. This identity can be shown by induction using
integration by parts (integrate the exponential factor). This formula would seem to allow for
efficient computation of the required M values of the lower incomplete gamma function.
However, the factorial term multiplied by the difference of two terms that might be close in
magnitude suggest that floating point precision could easily be lost as n increases. Numerical
experiments confirmed this and the algorithm was deemed unsuitable.
However, one can use integration by parts on Equation 1.21 (integrating the polynomial factor), to
derive
N−1
xn+k
Z x
−x 1
γ(n, x) = e ∑ + e−z zn+N−1 dz (1.25)
k=0 n(n + 1) · · · (n + k) n(n + 1) · · · (n + N − 1) 0
N N
= S +E (1.26)

where the sum SN has N terms and E N is the final integral term in Equation 1.25. A simple
estimate gives
xn+N−1 (1 − e−x )
EN ≤ . (1.27)
n(n + 1) · · · (n + N − 1)
An alternate notation is now convenient. If we define

γ2 (n, x) = x−n ex γ(n, x) (1.28)

then
N−1
xk ex x−n
Z x
γ2 (n, x) = ∑ + e−z zn+N−1 dz (1.29)
k=0 n(n + 1) · · · (n + k) n(n + 1) · · · (n + N − 1) 0

= S2N + E2N (1.30)

Similarly, we have
xN−1 ex (1 − e−x )
E2N ≤ (1.31)
n(n + 1) · · · (n + N − 1)
and to achieve double precision relative accuracy for γ2 (n, x), we stop adding terms when

E2N ≤ 10−15 S2N . (1.32)

We have also
eΓ0 η In = γ2 (n + 1, Γ0 η)η n+1 (1.33)
which leads to the convenient form

K0 M
TH = T0 eΓ0 η + ∑ iKiη i+2γ2(i + 2, Γ0η)
2Cv ρ0 i=1
(1.34)

Since we potentially have a large number of terms to evaluate, we look for even more efficiency
by looking for a stable backward recurrence for γ2 that will allow for computing γ2 (n, x) at the
same time that Equation 1.34 is evaluated as a telescoping sum.

10
Integrating by parts on Equation 1.21 (integrating the polynomial factor once), we derive
1
γ(n, x) = γ(n + 1, x) + e−x xn


(1.35)
n
and the corresponding recursion relation for γ2 is
x 1
γ2 (n, x) = γ2 (n + 1, x) + . (1.36)
n n
which is a backward recurrence with a starting value γ2 (M + 2, x) evaluated using Equation 1.29.
Let ε2 (n, x) be difference between the numerical solution and the exact solution for γ2 . This error
satisfies the backward recurrence relation
x
ε2 (n, x) = ε2 (n + 1, x) (1.37)
n
with solution
xM+2−n
ε2 (n, x) = ε2 (M + 2, x), 3 ≤ n ≤ M + 1 (1.38)
n(n + 1)(n + 2) · · · (M + 1)
Utilizing Equation 1.24 to switch to a relative error ε2r , we obtain
1 − e−x eM+1 (x) r
 
r
ε2 (n, x) = ε2 (M + 2, x), 3 ≤ n ≤ M + 1 (1.39)
1 − e−x en−1 (x)
and it is seen that the backwards recurrence algorithm is stable and has a well controlled error. As
a check, numerical results from the recurrence algorithm were compared to a direct evaluation of
Equation 1.29. The numerical stability of recurrence relations has an extensive literature [6].
As an example of the power series form for the Hugoniot, we can compute the power series
expansion of the linear Us − u p Hugoniot in terms of η. The relevant shock jump equations are,
respectively, mass, momentum conservation and the assumed linear Us −U p relation:
ρ0Us = ρ(Us − u p ), (1.40)
PH = ρUs u p (1.41)
and
Us = C0 + su p (1.42)
where Us is the shock speed, C0 is the sound speed at the initial (reference) state, s is the linear
coefficient and u p is the particle velocity. Eliminating Us from Equation 1.40 results in
C0 η
up = (1.43)
1 − sη
and substitution of Equation 1.42 and then Equation 1.43 in Equation 1.41 gives the Hugoniot
entirely in terms of η
ρ0C02 η
PH = (1.44)
(1 − sη)2
Expanding Equation 1.44 in a power series in η results in
PH = ρ0C02 η(1 + sη + (sη)2 + (sη)3 + · · ·)2 (1.45)
= ρ0C02 η(1 + 2sη + 3(sη)2 + 4(sη)3 + · · ·) (1.46)
so that Kn = (n + 1)sn for n ≥ 1.

11
1.3. EXPANSION
Pmin
For ηmin = K0 ≤ η < 0, the reference curve is defined by an isentrope with a single fixed K0 .

PR = PI = K0 η (1.47)

Equation 1.7 in the form

dE
= v0 P (1.48)
dη
yields
K0 η 2
ER = EI = + E0 (1.49)
2ρ0
and Equation 1.12 in the form
dT
= Γ0 T (1.50)
dη
results in
TR = TI = T0 eΓ0 η (1.51)

1.4. PRESSURE CUTOFF IN TENSION

Similarly, for η < ηmin = PKmin

0
, we define a continuation of the isentropic reference curve which
does not sustain additional tensile pressure. Thus

PR = PI = Pmin (1.52)
2
K0 ηmin Pmin
ER = EI = + E0 + (η − ηmin ) (1.53)
2ρ0 ρ0
TR = TI = T0 eΓ0 η (1.54)

1.5. CONCLUSION

We have documented details of a Mie-Grüneisen equation of state with a polynomial power series
compressive Hugoniot reference curve and a linear isentropic expansion reference curve with a
pressure cutoff in tension. A constant ρΓ product and constant heat capacity is assumed. The
numerical approaches utilized in the evaluating the analytic solution are fully described.

12
REFERENCES

[1] R. Courant and K. O. Friedrichs. Supersonic flow and shock waves. Springer-Verlag, 1948.
[2] M.H. Rice, R.G. McQueen, and J.M. Walsh. Compression of solids by strong shock waves.
In Frederick Seitz and David Turnbull, editors, Advances in Research and Applications,
volume 6 of Solid State Physics, pages 1 – 63. Academic Press, 1958.
[3] Allen C. Robinson. Evaluation techniques and properties of an exact solution to a subsonic
free surface jet flow. Technical report SAND2002-1015, Sandia National Laboratories,
Albuquerque, NM, April 2002.
[4] J. M. Walsh and R. H. Christian. Equation of state of metals from shock wave measurements.
Physical Review, 97(6):1544–1556, 1955.
[5] Eric W. Weisstein. Incomplete gamma function. From MathWorld- A Wolfram Web
Resource. http://mathworld.wolfram.com/IncompleteGammaFunction.html.
[6] Jet Wimp. Computation with Recurrence Relations. Pitman, 1984.

13
DISTRIBUTION

Email—External (encrypt for OUO)

Company Email
Name Company Name
Address

Ann E. Mattsson [email protected] Los Alamos National Laboratory

Email—Internal (encrypt for OUO)

Name Org. Sandia Email Address

Steven W. Bova 1443 [email protected]

James B. Carleton 1443 [email protected]

Timothy J. Fuller 1443 [email protected]

Glen Hansen 1443 [email protected]

Daniel A. Ibanez 1443 [email protected]

Edward Love 1443 [email protected]

Christopher B. Luchini 1443 [email protected]

Duncan A. O. McGregor 1443 [email protected]

John H. Niederhaus 1443 [email protected]

Michael Powell 1443 [email protected]

Angel E. Rodriquez 1443 [email protected]

Jason J. Sanchez 1443 [email protected]

Alan K. Stagg 1443 [email protected]

Thomas E. Voth 1443 [email protected]

John H. Carpenter 1444 [email protected]

Technical Library 01177 [email protected]

14
15
Sandia National Laboratories is a
multimission laboratory managed
and operated by National
Technology & Engineering
Solutions of Sandia LLC, a wholly
owned subsidiary of Honeywell
International Inc., for the U.S.
Department of Energy’s National
Nuclear Security Administration
under contract DE-NA0003525.