EOS Models

The mathematical descriptions of these models are presented below while the details of using them is presented in the description of the EOS API.

EOS Primer

An equation of state (EOS) is a constituitive model that generally relates thermodynamic quantities such as pressure, temperature, density, internal energy, free energy, and entropy consistent with the constraints of equillibrium thermodynamics.

singularity-eos contains a number of equations of state that are most useful in hydrodynamics codes. All of the equations of state presented here are considered “complete” equations of state in the sense that they contain all the information needed to derive a given thermodynamic quantity. However, not all variables are exposed since the emphasis is on only those that are needed for hydrodynamics codes. An incomplete equation of state is often sufficient for pure hydrodynamics since it can relate pressure to a given energy-density state, but it is missing information that relates the energy to temperature (which in turn provides a means to calculate entropy).

An excellent resource on equations of state in the context of hydrodynamics is the seminal work from Menikoff and Plohr. In particular, Appendix A contains a number of thermodynamic relationships that can be useful for computing additional quantities from those output in singularity-eos.

The Mie-Gruneisen form

Many of the following equations of state are considered to be of the “Mie-Gruneisen form”, which has many implications but for our purposes means that the general form of the EOS is

\[P - P_\mathrm{ref} = \rho \Gamma(\rho) (e - e_\mathrm{ref})\]

where ‘ref’ denotes quantities along some reference curve, \(P\) is the pressure, \(\rho\) is the density, \(\Gamma\) is the Gruneisen parameter, and \(e\) is the specific internal energy. In this sense, an EOS of this form uses the Gruneisen parameter to describe the pressure behavior of the EOS away from the reference curve. Coupled with a relationship between energy and temperature (sometimes as simple as a constant heat capacity), the complete equation of state can be constructed.

To some degree it is the complexity of the reference state and the heat capacity that will determine an EOS’s ability to capture the complex behavior of a material. At the simplest level, the ideal gas EOS uses a reference state at zero pressure and energy, while more complex equations of state such as the Davis EOS use the material’s isentrope. In ths way, the reference curve indicates the conditions under which you can expect the EOS to represent the intended behavior.

Nomenclature

The EOS models in singularity-eos are defined for the following sets of dependent and independent variables through various member functions described in the EOS API.

Function	Dependent Variable	Independent Variables
\(T(\rho, e)\)	Temperature	Density, Internal Energy
\(P(\rho, e)\)	Pressure	Density, Internal Energy
\(e(\rho, T)\)	Internal Energy	Density, Temperature
\(P(\rho, T)\)	Pressure	Density, Temperature
\(\rho(P, T)\)	Density	Pressure, Temperature
\(e(P, T)\)	Internal Energy	Pressure, Temperature
\(C_V(\rho, T)\)	Constant Volume Specific Heat Capacity	Density, Temperature
\(C_V(\rho, e)\)	Constant Volume Specific Heat Capacity	Density, Internal Energy
\(B_S(\rho, T)\)	Isentropic Bulk Modulus	Density, Temperature
\(B_S(\rho, e)\)	Isentropic Bulk Modulus	Density, Internal Energy
\(\Gamma(\rho, T)\)	Gruneisen Parameter	Density, Temperature
\(\Gamma(\rho, e)\)	Gruneisen Parameter	Density, Internal Energy

A point of note is that “specific” implies that the quantity is intensive on a per unit mass basis. It should be assumed that the internal energy is always specific since we are working in terms of density (the inverse of specific volume).

Disambiguation

Gruneisen Parameter

In this description of the EOS models, we use \(\Gamma\) to represent the Gruneisen coeficient since this is the most commonly-used symbol in the context of Mie-Gruneisen equations of state. The definition of the Gruneisen parameter is

\[\Gamma := \frac{1}{\rho} \left( \frac{\partial P}{\partial e} \right)_\rho\]

This should be differentiated from

\[\gamma := \frac{V}{P} \left( \frac{\partial P}{\partial V} \right)_S = \frac{B_S}{P}\]

though, which is the adiabatic exponent.

For an ideal gas, the adiabatic exponent is simply the ratio of the heat capacities,

\[\gamma_\mathrm{id} = \frac{C_P}{C_V} = \frac{B_S}{B_T}.\]

Here \(C_P\) is the specific heat capacity at constant pressure and \(B_T\) is the isothermal bulk modulus.

Units and conversions

The default units for singularity-eos are cgs which results in the following units for thermodynamic quantities:

Quantity	Default Units	cgµ conversion	mm-mg-µs conversion
\(P\)	µbar	10^-12 Mbar	10^-10 GPa
\(\rho\)	g/cm³	1	1 mg/mm³
\(e\)	erg/g	10^-12 Terg/g	10^-10 J/mg
\(T\)	K	1	1
\(C_V\)	erg/g/K	10^-12 Terg/g/K	10^-10 J/mg/K
\(B_S\)	µbar	10^-12 Mbar	10^-10 GPa
\(\Gamma\)	unitless	–	–

Note: some temperatures are measured in eV for which the conversion is 8.617333262e-05 eV/K.

Sesame units are equivalent to the mm-mg-µs unit system.

Implemented EOS models

Ideal Gas

The ideal gas (aka perfect or gamma-law gas) model in singularity-eos takes the form

\[P = \Gamma \rho e\]

\[e = C_V T,\]

where quantities are defined in the nomenclature section.

The settable parameters are the Gruneisen parameter and specific heat capacity.

Although this differs from the traditional representation of the ideal gas law as \(P\tilde{V} = RT\), the forms are equivalent by recognizing that \(\Gamma = \frac{R}{\tilde{C_V}}\) where \(R\) is the ideal gas constant in units of energy per mole per Kelvin and \(\tilde{C_\mathrm{V}}\) is the molar heat capacity, relatable to the specific heat capacity through the molecular weight of the gas. Since \(\tilde{C_\mathrm{V}} = \frac{5}{2} R\) for a diatomic ideal gas, the corresponding Gruneisen parameter should be 0.4.

A common point of confusion is \(\Gamma\) versus \(\gamma\) with the latter being the adiabatic exponent. For an ideal gas, they are related through

\[\Gamma = \gamma - 1\]

The IdealGas EOS constructor has two arguments, gm1, which is the Gruneisen parameter \(\Gamma\), and Cv, which is the specific heat \(C_V\):

IdealGas(Real gm1, Real Cv)

Gruneisen EOS

One of the most commonly-used EOS to represent solids is the Steinberg variation of the Mie-Gruneisen EOS, often just shortened to “Gruneisen” EOS. This EOS uses the Hugoniot as the reference curve and thus is extremly powerful because the basic shock response of a material can be modeled using minimal parameters.

The pressure follows the traditional Mie-Gruneisen form,

\[P(\rho, e) = P_H(\rho) + \rho\Gamma(\rho) \left(e - e_H(\rho) \right),\]

Here the subscript \(H\) is a reminder that the reference curve is a Hugoniot. Other quantities are defined in the nomenclature section.

The above is an incomplete equation of state because it only relates the pressure to the density and energy, the minimum required in a solution to the Euler equations. To complete the EOS and determine the temperature, a constant heat capacity is assumed so that

\[T(\rho, e) = \frac{e}{C_V} + T_0\]

The user should note that this implies that \(e=0\) at the reference temperature, \(T_0\). Given this simple relationship, the user should treat the temperature from this EOS as only a rough estimate.

Given a reference density, \(\rho_0\), we first parameterize the EOS using \(\eta\) as a measure of compression given by

\[\eta = 1 - \frac{\rho_0}{\rho}.\]

This is convenient because \(eta = 0\) when \(\rho = \rho_0\), \(\eta = 1\) at the infinite density limit, and \(\eta = -\infty\) at the zero density limit. The Gruneisen parameter, \(\Gamma\) can be expressed in terms of \(\eta\) as

\[\begin{split}\Gamma(\rho) = \begin{cases} \Gamma_0 & \eta \leq 0 \\ \Gamma_0 (1 - \eta) + b\eta & 0 \leq \eta < 1 \end{cases}\end{split}\]

When the unitless user parameter \(b=0\), the Gruneisen parameter is of a form where \(\rho\Gamma =\) constant in compression, i.e. when \(\eta > 0\).

The reference pressure along the Hugoniot is determined by

\[\begin{split}P_H(\rho) = P_0 + c_0^2 \eta \begin{cases} \rho & \rho < \rho_0 \\ \frac{\rho_0}{\left( 1 - s_1 \eta - s_2 \eta^2 - s_3 \eta^3 \right)^2} & \rho \geq \rho_0 \end{cases}\end{split}\]

where \(P_0\) is the reference pressure and \(c_0\), \(s_1\), \(s_2\), and \(s_3\) are fitting paramters to the \(U_s\)-\(u_p\) curve such that

\[U_s = c_0 + u_p \left( s_1 + s_2 \frac{u_p}{U_s} + s_3\left(\frac{u_p}{U_s}\right)^2 \right).\]

Here \(U_s\) is the shock velocity and \(u_p\) is the particle velocity. For many materials, this relationship is roughly linear so only the \(s_1\) parameter is needed. The units for \(c_0\) are velocity while the rest are unitless.

Finally the energy along the Hugoniot is given by

\[\begin{split}E_H(\rho) = \begin{cases} 0 & \rho < \rho_0 \\ \frac{\eta (P_H + P_0)}{2 \rho_0} & \rho \geq \rho_0 \end{cases}.\end{split}\]

One should note that in this form neither the expansion region nor the overall temperature are thermodynamically consistent with the rest of the EOS. Since the EOS is a fit to the principal Hugoniot, the EOS will obviously reproduce single shocks quite well, but it may not be as appropriate when there are multiple shocks or for modeling the release behavior of a material.

The constructor for the Gruneisen EOS has the signature

Gruneisen(const Real C0, const Real s1, const Real s2, const Real s3, const Real G0,
          const Real b, const Real rho0, const Real T0, const Real P0, const Real Cv,
          const Real rho_max)

where C0 is \(C_0\), s1 is \(s_1\), s2 is \(s_2\), s3 is \(s_3\), G0 is \(\Gamma_0\), b is \(b\), rho0 is \(\rho_0\), T0 is \(T_0\), P0 is \(P_0\), and Cv is \(C_v\). rho_max is the maximum value of density for which the reference pressure curve is valid. Input densities above rho_max are pinned to rho_max.

There is an overload of the Gruneisen class which computes rho_max automatically without the user needing to specify:

Gruneisen(const Real C0, const Real s1, const Real s2, const Real s3, const Real G0,
          const Real b, const Real rho0, const Real T0, const Real P0, const Real Cv)

JWL EOS

The Jones-Wilkins-Lee (JWL) EOS is used mainly for detonation products of high explosives. Similar to the other EOS here, the JWL EOS can be written in a Mie-Gruneisen form as

\[P(\rho, e) = P_S(\rho) + \rho w (e - e_S(\rho))\]

where the reference curve is an isentrope of the form

\[P_S(\rho) = A \mathrm{e}^{-R_1 \eta} + B \mathrm{e}^{-R_2 \eta}\]

\[e_S(\rho) = \frac{A}{\rho_0 R_1} \mathrm{e}^{-R_1 \eta} + \frac{B}{\rho_0 R_2} \mathrm{e}^{-R_2 \eta}.\]

Here \(\eta = \frac{\rho_0}{\rho}\) and \(R_1\), \(R_2\), \(A\), \(B\), and \(w\) are constants particular to the material. Note that the parameter \(w\) is simply the Gruneisen parameter and is assumed constant for the EOS (which is fairly reasonable since the detonation products are gasses).

Finally, to complete the EOS the energy is related to the temperature by

\[e = e_S(\rho) + C_V T\]

where \(C_V\) is the constant volume specific heat capacity.

The constructor for the JWL EOS is

JWL(const Real A, const Real B, const Real R1, const Real R2,
    const Real w, const Real rho0, const Real Cv)

where A is \(A\), B is \(B\), R1 is \(R_1\), R2 is \(R_2\), w is \(w\), rho0 is \(\rho_0\), and Cv is \(C_V\).

Davis EOS

The Davis reactants and products EOS are both of Mie-Gruneisen forms that use isentropes for the reference curves. The equations of state are typically used to represent high explosives and their detonation products and the reference curves are calibrated to several sets of experimental data.

For both the reactants and products EOS, the pressure and energy take the forms

\[P(\rho, e) = P_S(\rho) + \rho\Gamma(\rho) \left(e - e_S(\rho) \right)\]

\[e(\rho, P) = e_S(\rho) + \frac{1}{\rho \Gamma(\rho)} \left(P - P_S(\rho) \right),\]

where the subscript \(S\) denotes quantities along the reference isentrope and other quantities are defined in the nomenclature section.

Davis Reactants EOS

The Davis reactants EOS uses an isentrope passing through a reference state and assumes that the heat capacity varies linearly with entropy such that

\[C_V = C_{V,0} + \alpha(S - S_0),\]

where subscript \(0\) refers to the reference state and \(\alpha\) is a dimensionless constant specified by the user.

The \(e(\rho, P)\) lookup is quite awkward, so the energy is more-conveniently cast in terms of termperature such that

\[e(\rho, T) = e_S(\rho) + \frac{C_{V,0} T_S(\rho)}{1 + \alpha} \left( \left(\frac{T}{T_S(\rho)} \right)^{1 + \alpha} - 1 \right),\]

which can easily be inverted to find \(T(\rho, e)\).

The Gruneisen parameter takes on a linear form such that

\[\begin{split}\Gamma(\rho) = \Gamma_0 + \begin{cases} 0 & \rho < \rho_0 \\ Zy & \rho >= \rho_0 \end{cases}\end{split}\]

where \(Z\) and \(y\) are dimensionless parameters.

Finally, the pressure, energy, and temperature along the isentrope are given by

\[\begin{split}P_S(\rho) = P_0 + \frac{\rho_0 A^2}{4B} \begin{cases} \mathrm{e} \left( 4By \right) -1 & \rho < \rho_0 \\ \sum\limits_{j=1}^3 \frac{(4By)^j}{j!} + C\frac{(4By)^4}{4!} + \frac{y^2}{(1-y)^4} & \rho >= \rho0 \end{cases}\end{split}\]

\[e_S(\rho) = e_0 + \int\limits_{\rho_0}^{\rho} \frac{P_S(\bar{\rho})}{\bar{\rho^2}}~\mathrm{d}\bar{\rho}\]

\[\begin{split}T_S(\rho) = T_0 \begin{cases} \left(\frac{\rho}{\rho_0} \right)^{\Gamma_0} & \rho < \rho_0 \\ \mathrm{e} \left( -Zy \right) \left(\frac{\rho}{\rho_0} \right)^{\Gamma_0 + Z} & \rho >= \rho_0 \end{cases}\end{split}\]

where \(A\), \(B\), \(C\), \(y\), and \(Z\) are all user-settable parameters and again quantities with a subcript of \(0\) refer to the reference state. The variable \(\bar{\rho}\) is simply an integration variable. The parameter \(C\) is especially useful for ensuring that the high-pressure portion of the shock Hugoniot does not cross that of the products.

The settable parameters are the dimensionless parameters listed above as well as the pressure, density, temperature, energy, Gruneisen parameter, and constant volume specific heat capacity at the reference state.

The constructor for the Davis Reactants EOS is

DavisReactants(const Real rho0, const Real e0, const Real P0, const Real T0,
               const Real A, const Real B, const Real C, const Real G0, const Real Z,
               const Real alpha, const Real Cv0)

where rho0 is \(\rho_0\), e0 is \(e_0\), P0 is \(P_0\), T0 is \(T_0\), A is \(A\), B is \(B\), C is \(C\), G0 is \(\Gamma_0\), Z is \(Z\), alpha is \(\alpha\), and Cv0 is the specific heat capacity at the reference state.

Davis Products EOS

The Davis products EOS is created from the reference isentrope passing through the CJ state of the high explosive along with a constant heat capacity. The constant heat capacity leads to the energy being a simple funciton of the temperature deviation from the reference isentrope such that

\[e(\rho, T) = e_S(\rho) + C_{V,0} (T - T_S(\rho)).\]

The Gruneisen parameter is given by

\[\Gamma(\rho) = k - 1 + (1-b) F(\rho)\]

where \(b\) is a user-settable dimensionless parameter and \(F(\rho)\) is given by

\[F(\rho) = \frac{2a (\rho V_{\mathrm{C}})^n}{(\rho V_{\mathrm{C}})^{-n} + (\rho V_{\mathrm{C}})^n}.\]

Here the calibration parameters \(a\) and \(n\) are dimensionless while \(V_{\mathrm{C}}\) is given in units of specific volume.

Finally, the pressure, energy, and temperature along the isentrope are given by

\[P_S(\rho) = P_{\mathrm{C}} G(\rho) \frac{k - 1 + F(\rho)}{k - 1 + a}\]

\[e_S(\rho) = e_{\mathrm{C}} G(\rho) \frac{1}{\rho V_{\mathrm{C}}}\]

\[T_S(\rho) = T_{\mathrm{C}} G(\rho) \frac{1}{(\rho V_{\mathrm{C}})^{ba + 1}}\]

where

\[G(\rho) = \frac{ \left( \frac{1}{2}(\rho V_{\mathrm{C}})^{-n} + \frac{1}{2}(\rho V_{\mathrm{C}})^n \right)^{a/n}} {(\rho V_{\mathrm{C}})^{-(k+a)}}\]

and

\[e_{\mathrm{C}} = \frac{P_{\mathrm{C}} V_{\mathrm{C}}}{k - 1 + a}.\]

Here, there are four dimensionless parameters that are settable by the user, \(a\), \(b\), \(k\), and \(n\), while \(P_\mathrm{C}\), \(e_\mathrm{C}\), \(V_\mathrm{C}\) and \(T_\mathrm{C}\) are tuning parameters with units related to their non-subscripted counterparts.

The constructor for the Davis Products EOS is

DavisProducts(const Real a, const Real b, const Real k, const Real n, const Real vc,
              const Real pc, const Real Cv, const Real E0)

where a is \(a\), b is \(b\), k is \(k\), n is \(n\), vc is \(V_\mathrm{C}\), pc is \(P_\mathrm{C}\), Cv is \(C_{V,0}\), and E0 is \(e_\mathrm{C}\).

Spiner EOS

Spiner EOS is a tabulated reader for the Sesame database of material equations of state. Materials include things like water, dry air, iron, or steel. This model comes in two flavors: SpinerEOSDependsRhoT and SpinerEOSDependsRhoSie. The former tabulates all quantities of interest in terms of density and temperature. The latter also includes tables in terms of density and specific internal energy.

Tabulating in terms of density and pressure means that computing, e.g., pressure in terms of density and internal energy requires solving the equation:

\[e_0 = e(\rho, T)\]

for temperature \(T\) given density \(\rho\) and specific internal energy \(e_0\). This is in general not closed algebraically and must be solved using a root-find. SpinerEOSDependsRhoT performs this root find in-line, and the result is performant, thanks to library’s ability to take and cache initial guesses. SpinerEOSDependsRhoSie circumvents this issue by tabulating in terms of both specific internal energy and temperature.

Both models use (approximately) log-linear interpolation on a grid that is (approximately) uniformly spaced on a log scale. Thermodynamic derivatives are tabulated and interpolated, rather than computed from the interpolating function. This approach allows for significantly higher fidelity approximations of these derivatives.

Both SpinerEOS classes benefit from a lambda parameter, as described in the EOS API section`. In particular, if an array of size 2 is passed in to the scalar call (or one per point for the vector call), the model will leverage this scratch space to cache initial guesses for root finds.

To avoid race conditions, at least one array should be allocated per thread. Depending on the call pattern, one per point may be best. In the vector case, one per point is necessary.

The constructor for SpinerEOSDependsRhoT is given by two overloads:

SpinerEOSDependsRhoT(const std::string &filename, int matid,
                     bool reproduciblity_mode = false);
SpinerEOSDependsRhoT(const std::string &filename, const std::string &materialName,
                     bool reproducibility_mode = false);

where here filename is the input file, matid is the unique material ID in the database in the file, materialName is the name of the material in the file, and reproducability_mode is a boolean which slightly changes how initial guesses for root finds are computed. The constructor for SpinerEOSDependsRhoSie is identical.

`sp5` files and `sesame2spiner`

The SpinerEOS models use their own file format built on hdf5, which we call sp5. These files can be generated by hand, or they can be generated from the sesame database (assuming eospac is installed) via the tool sesame2spiner, which is packaged with singularity-eos. Buld sesame2spiner by specifying

-DSINGULARITY_USE_HDF5=ON -DSPINGULARITY_USE_EOSPAC=ON -DSINGULARITY_BUILD_SESAME2SPINER=ON

at configure time. The call to sesame2spiner is of the form

sesame2spiner -s output_file_name.sp5 input1.dat input2.dat ...

for any number of input files. Verbosity flags -p and -v are also available. Use -h for a help message. The -s flag is optional and the output file name defaults to materials.sp5.

Each input file corresponds to a material and consists of simple key-value pairs. For exampe the following input deck is for air:

matid = 5030
# These set the number of grid points per decade
# for each variable. The default is 50 points
# per decade.
numrho/decade = 40
numT/decade = 40
numSie/decade = 40
# Defaults pulled from the sesame file if possible
name = air
rhomin = 1e-2
rhomax = 10
Tmin = 252
Tmax = 1e4
siemin = 1e12
siemax = 1e16
# These shrink the logarithm of the bounds by a fraction of the
# total inteval <= 1.
# Note that these may be deprecated in the near future.
shrinklRhoBounds = 0.15
shrinklTBounds = 0.15
shrinkleBounds = 0.5

The only required value in an input file is the matid, in this case 5030. All other values will be inferred from the original sesame database if possible and if no value in the input file is provided. Comments are prefixed with #.

eospac uses environment variables and files to locate files in the sesame database, and sesame2spiner uses eospac. So the location of the sesame database need not be provided by the command line. For how to specify sesame file locations, see the eospac manual.

Stellar Collapse EOS

This model provides finite temperature nuclear equations of state suitable for core collapse supernova and compact object (such as neutron star) simulations. These models assume nuclear statistical equilibrium (NSE). It reads tabulated data in the Stellar Collapse format, as first presented by OConnor and Ott.

Like SpinerEOSDependsRhoT, StellarCollapse tabulateds all quantities in terms of density and temperature on a logarithmically spaced grid. And similarly, it requires an in-line root-find to compute quantities in terms of density and specific internal energy. Unlike most of the other models in singularity-eos, StellarCollapse also depends on a third quantity, the electron fraction,

\[Y_e = \frac{n_e}{n_p + n_n}\]

which measures the number fraction of electrons to baryons. Symmetric matter has a \(Y_e\) of 0.5, while cold neutron stars, have a \(Y_e\) approximately less than 0.1.

As with SpinerEOSDependsRhoT, the Stellar Collapse tables tabulate thermodynamic derivatives separately, rather than reconstruct them from interpolants. However, the tabulated values can contain artifacts, such as unphysical spikes. To mitigate this issue, the thermodynamic derivatives are cleaned via a median filter. The bulk modulus is then recomputed from these thermodynamic derivatives via:

\[B_S(\rho, T) = \rho \left(\frac{\partial P}{\partial\rho}\right)_e + \frac{P}{\rho} \left(\frac{\partial P}{\partial e}\right)_\rho\]

Note that StellarCollapse is a relativistic model, and thus the sound speed is given by

\[c_s^2 = \frac{B_S}{w}\]

where \(w = \rho h\) for specific entalpy \(h\) is the enthalpy by volume, rather than the density \(rho\). This ensures the sound speed is bounded from above by the speed of light.

The StellarCollapse model requires a lambda parameter of size 2, as described in the EOS API section`. The zeroth element of the lambda array contains the electron fraction. The first element is reserved for caching. It currently contains the natural log of the temperature, but this should not be assumed.

To avoid race conditions, at least one array should be allocated per thread. Depending on the call pattern, one per point may be best. In the vector case, one per point is necessary.

The StellarCollpase model can read files in either the original format found on the Stellar Collapse website, or in the sp5 format described above.

Warning

Note that the data contained in an sp5 file for the StellarCollapse EOS and the SpinerEOS models is not identical and the files are not interchangeable.

The constructor for the StellarCollapse EOS class looks like

StellarCollapse(const std::string &filename, bool use_sp5 = false,
                bool filter_bmod = true)

where filename is the file containing the tabulated model, use_sp5 specifies whether to read an sp5 file or a file in the original Stellar Collapse format, and filter_bmod specifies whether or not to apply the above-described median filter.

StellarCollapse also provides

void Save(const std::string &filename)

which saves the current EOS data in sp5 format.

EOSPAC EOS

This is a striaghtforward wrapper of the EOSPAC library for the Sesame database. The constructor for the EOSPAC model looks like

EOSPAC(int matid, bool invert_at_setup = false)

where matid is the unique material number in the database and invert_at_setup specifies whether or not pre-compute tables of temperature as a function of density and energy.