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 Theory

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 thermal information (see the section on complete equations of state).

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.

Some Notes on Thermodynamic Consistency

For the pure purpose of solving the Euler equations, an incomplete equation of state of the form \(P(\rho, e)\) is sufficient. In essence, this provides the mechanical response of a material subjected to different types of compression and expansion. However, this EOS is lacking thermal information without an appropriate heat capacity relationship.

As discussed by Mattsson, an equation of state can be considered “complete” if it is constructed from one of the thermodynamic potentials using their natural variables, i.e.

\[\begin{split}e =& e(\rho, S) \\ h =& h(P, S) \\ F =& F(\rho, T) \\ G =& G(P, T) \\\end{split}\]

where all potentials are specific. Here \(e\) is again the internal energy, \(h\) is the enthalpy, \(F\) is the Helmholtz free energy, and \(G\) is the Gibbs free energy. While equations of state formulated using the Helmholtz free energy can be particularly attractive (such as the sesame tables), finding a convenient form can be difficult. As such, it becomes imperitive to extend the Mie-Gruneisen form so that it can form a complete EOS.

The heat capacity is defined as

\[C_V := \left(\frac{\partial e}{\partial T}\right)_V = T \left(\frac{\partial S}{\partial T}\right)_V\]

which provides a natural means by which to relate the entropy/energy to the temperature and form a complete equation of state. However, there are specific requirements placed on the functional form of the heat capacity described below.

The Maxwell relations are the consequence of the requirement that mixed partial second derivatives of the thermodynamic potentials should be independent of the order of differentiation. This concept can be extended to third derivatives (as Davis does in his Appendix B here) to produce the relationship,

\[\frac{V}{C_V^2}\left(\frac{\partial C_V}{\partial V}\right)_S = \left(\frac{\partial \Gamma}{\partial S}\right)_V.\]

This is often referred to as a “compatibility condition” (see also Menikoff) and provides an important connection between the Gruneisen parameter, \(\Gamma\), and the constant volume heat capacity, \(C_V\), in the context of a complete equation of state developed from the internal energy. Importantly, if the Gruneisen form forces us to assume that the Gruneisen parameter is a sole function of density/volume, then the implication is that the heat capacity *must* then be a sole function of entropy.

Again from Menikoff, it can be shown that even in the simplest possible case of a constant heat capacity, the temperature is related to the energy through

\[T(\rho, e) = T_0 \phi(\rho) + \frac{e - e_{S,0}(\rho)}{C_V}\]

where \(T_0\) represents a reference temperature, \(e_{S,0}\) is the energy along the isentrope that passes through the reference temperature, and \(\phi(\rho)\) is an integrating factor given by

\[\phi(\rho) = \exp\left(\int\limits_{\rho_0}^{\rho} \rho \Gamma (\rho) \mathrm{d}\rho \right).\]

As an EOS of a Mie-Gruneisen form becomes more complicated with more complex functional forms for \(\Gamma\) and the reference curves, the task of calculating a thermodynamically consistent temperature becomes more complicated.

Available EOS Information and 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, Specific Internal Energy
\(P(\rho, e)\)	Pressure
\(S(\rho, e)\)	Specific Entropy
\(e(\rho, T)\)	Specific Internal Energy	Density, Temperature
\(P(\rho, T)\)	Pressure
\(S(\rho, T)\)	Specific Entropy
\(\rho(P, T)\) [1]	Density	Pressure, [1] Temperature [1]
\(e(P, T)\) [1]	Specific Internal Energy
\(C_V(\rho, T)\)	Constant Volume Specific Heat Capacity	Density, Temperature
\(C_V(\rho, e)\)		Density, Specific Internal Energy
\(B_S(\rho, T)\)	Isentropic Bulk Modulus	Density, Temperature
\(B_S(\rho, e)\)		Density, Specific Internal Energy
\(\Gamma(\rho, T)\)	Gruneisen Parameter	Density, Temperature
\(\Gamma(\rho, e)\)		Density, Specific Internal Energy

[1] (1,2,3,4)

Note: Using pressure and temperature as independent variables is fraught since both pressure and energy are often multi-valued in density for many EOS due to the presence of phase changes (especially tabular EOS). For analytic EOS in singularity-eos where there is not an analytic inversion to pressure-temperature space, a root-find is typically used that uses the density at standard temperature and pressure (STP) as an initial guess. The notable exceptions to this are the spiner EOS that allow the introduction of an initial guess via the lambda function argument.

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 and entopry are always specific since we are working in terms of density (the inverse of specific volume).

Entropy availability

For an arbitrary equation of state, a change in entropy in terms of temperature and volume is given by

\[\Delta S = \int\limits_{S_0}^{S} \mathrm{d}S = \int\limits_{T_0}^{T} \left(\frac{\partial S}{\partial T}\right)_V \mathrm{d}T + \int\limits_{V_0}^{V} \left(\frac{\partial S}{\partial V}\right)_T \mathrm{d}V,\]

which can be simplified using a definition of the heat capacity and a Maxwell relation to become

\[\Delta S = \int\limits_{T_0}^{T} \frac{C_V}{T} \mathrm{d}T + \int\limits_{V_0}^{V} \left(\frac{\partial P}{\partial T}\right)_V \mathrm{d}V.\]

Similarly, expressing the entropy in terms of energy and volume yields

\[\Delta S = \int\limits_{e_0}^{e} \frac{1}{T(\rho, e)} \mathrm{d}e + \int\limits_{V_0}^{V} \frac{P(\rho, e)}{T(\rho, e)} \mathrm{d}V\]

after substituting for the appropriate derivatives using the first law of thermodynamics.

Importantly for the analytic EOS, these integrals will tend to diverge as the temperature and volume approach zero if the heat capacity does not also approach zero. This necessitates appropriate choices for the reference states \(T_0\) and \(V_0\).

Note

All EOS objects will expose the functions EntropyFromDensityInternalEnergy () and EntropyFromDensityTemperature() but many EOS cannot currently calculate entropy, either because the EOS is not thermodynamically consistent or because the feature has not yet been implemented. In these cases, the use of these functions will abort the calculation or raise an exception depending on the execution space (host or device). The cases where this occurs are noted below.

Nomenclature Disambiguation

The 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/cm3	1	1 mg/mm3
\(e\)	erg/g	10-12 Terg/g	10-10 J/mg
\(S\)	erg/g-K	10-12 Terg/g-K	10-10 J/mg-K
\(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: sometimes 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. 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\]

Although this formulation 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.

The entropy for an ideal gas is given by

\[S = C_V \ln\left(\frac{T}{T_0}\right) + \Gamma C_V \ln\left(\frac{\rho_0} {\rho}\right),\]

Note

The entropy diverges to negative infinity at absolute zero due to the constant heat capacity assumption. Care should be taken when using temperatures significantly below that of the reference state.

we have assumed that the entropy is zero at the reference state given by \(T_0\) and \(\rho_0\). By default, \(T_0 = 298\) K and the reference density is given by

\[\rho_0 = \frac{P_0}{\Gamma C_V T_0},\]

where \(P_0 = 1\) bar.

It should be noted that this equation diverges as the temperature approaches zero or the density approaches infinity. From a physical perspective, this is a limitation of the constant heat capacity assumption and would be remedied if the heat capacity went to zero at absolute zero. However, this represents a significant deviation from true ideal gas physics so we do not include it here.

The settable parameters for this EOS are the Gruneisen parameter and specific heat capacity. Optionally, the reference state for the entropy calculation can be provided by setting both the reference density and temperature.

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)

Optionally, the reference temperature and density for the entropy calculation, can be provided in the constructor via EntropyT0 and EntropyRho0:

IdealGas(Real gm1, Real Cv, Real EntropyT0, Real EntropyRho0)

Note that these parameters are provided solely for the entropy calculation. When these values are not set, they will be the same as those returned by the ValuesAtReferenceState() function. However, if the entropy reference conditions are given, the return values of the ValuesAtReferenceState() function will not be the same.

Gruneisen EOS

Warning

Entropy is not available for this 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 the inconsisetency in the temperature, we have made the choice not to expose the entropy for this EOS. Requesting an entropy value will result in an error.

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. Note that the parameter \(s_1\) is related to the fundamental derivative of shock physics as shown by Wills.

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)

Extendended Vinet EOS

The extended Vinet EOS is a full EOS, extended in both temperature and density from the Vinet universal EOS for solids (also called Rose cold curve). It is expected to work well in compression but is untested in expansion. It is published in Appendix 2 in J. Appl. Phys. 119, 015904 (2016).

While the Mie-Gruneisen EOS is based on a Hugoniot as reference curve, the Vinet is based on an isotherm:

\[P(\rho,T) = P_{ref}(\rho) + \alpha_0 B_0 (T - T_{ref})\]

where the reference isotherm is

\[\begin{split}P_{ref}(X)&=\frac{3 B_0}{X^2} Z \exp[\eta_0 Z] \left( 1 + \sum_{n=2}^N d_n Z^n \right) \, , \\ X &= \left( \frac{\rho_0}{\rho} \right)^{1/3} \\ Z &= 1-X\end{split}\]

Note that \(P_{ref}=0\) when \(\rho = \rho_0\), the reference state on the reference isotherm is always at ambient pressure. However, the reference isotherm is not necessarily at room temperature.

It can be shown that \(B_0\) is the isothermal bulk modulus, and \(\alpha_0\) the thermal expansion coefficient, at the reference state, and that

\[\eta_0 = \frac{3}{2}\left[ \left[ \frac{\partial B}{\partial P}\right]_0 -1\right] \, .\]

By assuming that also the constant volume heat capacity is a constant, \({C_V}_0\), an entropy can be derived

\[S(V,T) = S_0 + \alpha_0 B_0 (V - V_0) + {C_V}_0 \ln \frac{T}{T_{ref}}\]

and from that a thermodynamic consistent energy

\[\begin{split}E(X,T) =& 9 \frac{B_0 V_0}{{\eta_0}^2}\left(f_0 - \exp[\eta_0 Z] \left (f_0 - \eta_0 Z \left(f_0 + \sum_{n=1}^N f_n Z^n \right)\right)\right) \\ & - \alpha_0 B_0 V_0 (1-X^3) T_{ref} + {C_V}_0 (T - T_{ref}) + E_0\end{split}\]

where the energy coefficients \(f_n\) are determined from the pressure coefficients \(d_n\), \(n\geq 2\), by

\[\begin{split}f_N &= d_N \\ f_n &= d_n - \frac{n+2}{\eta_0} f_{n+1} \\ d_0 &= 1.0 \\ d_1 &= 0.0\end{split}\]

Note

The entropy diverges to negative infinity at absolute zero due to the constant heat capacity assumption. Care should be taken when using temperatures significantly below that of the reference state.

The constructor for the Vinet EOS has the signature

Vinet(const Real rho0, const Real T0, const Real B0, const Real BP0, const Real A0,
           const Real Cv0, const Real E0, const Real S0, const Real *expconsts)

where rho0 is \(\rho_0\), T0 is \(T_{ref}\), B0 is \(B_0\), BP0 is \((\partial B/\partial P)_0\), A0 is \(\alpha_0\), Cv0 is \({C_V}_0\), E0 is \(E_0\), S0 is \(S_0\), and expconsts is a pointer to the constant array of length 39 containing the expansion coefficients \(d_2\) to \(d_{40}\). Expansion coefficients not used should be set to 0.0.

JWL EOS

Warning

Entropy is not available for this 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

Warning

Entropy is not yet available for this EOS

The Davis reactants EOS uses an isentrope passing through a reference state and as the reference curve and then 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 Gruneisen parameter is given 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\) is a dimensionless parameter and \(y = 1 - \rho0/\rho\). Along an isentrope, the Gruneisen parameter can be expressed as

\[\Gamma_S(\rho) = \frac{\rho}{T} \left(\frac{\partial T}{\partial \rho}\right)_S,\]

which, upon integration can produce the temperature along the reference isentrope:

\[\begin{split}T_{S,0}(\rho) = \begin{cases} T_0\left(\frac{\rho}{\rho_0}\right)^{\Gamma_0} & \rho < \rho_0 \\ T_0\exp\left(-Zy\right)\left(\frac{\rho}{\rho_0}\right)^{\Gamma_0 + Z} & \rho >= \rho_0 \end{cases}\end{split}\]

where \(T_{S,0}\) is the temperature along the reference isentrope, \(S = S_0\).

Using the fact that the heat capacity can be expressed as

\[C_V = T\left( \frac{\partial S}{\partial T} \right)_V,\]

the temperature off of the reference isoentrope can be integrated from this identity to yield

\[T(\rho, S) = T_{S,0}(\rho) \left( \frac{C_V(S)}{C_{V,0}} \right)^{\frac{1}{\alpha}},\]

Now requiring that the entropy go to zero at absolute zero in accordance with the Nernst postulate and the third law of thermodynamics, the entropy can be expressed as a function of temperature and density such that

\[S(\rho, T) = \frac{C_{V,0}}{\alpha} \left( \frac{T}{T_{S,0}(\rho)} \right)^\alpha.\]

The \(e(\rho, P)\) formulation can now be 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)\).

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

\[\begin{split}P_S(\rho) = P_0 + \frac{\rho_0 A^2}{4B} \begin{cases} \exp \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},\]

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

Warning

Entropy is not yet available for this 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

Warning

Entropy is not yet available for this 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.

The StellarCollapse model, if used alone, also provides several additional functions of interest for those running, e.g., supernova simulations:

void MassFractionsFromDensityTemperature(const Real rho, const Real temperature, Real &Xa, Real &Xn, Real &Xp, Real &Abar, Real &Zbar, Real *lambda = nullptr) const

which returns the mass fractions for alpha particles, Xa, heavy ions Xh, neutrons Xn, and protons Xp, as well as the average atomic mass Abar and atomic number Zbar for heavy ions, assuming nuclear statistical equilibrium.

In addition, the user may query the bounds of the table via the functions lRhoMin(), lRhoMax(), lTMin(), lTMax(), TMin(), TMax(), YeMin(), YeMax(), sieMin(), and sieMax(), which all return a Real number. The l prefix indicates log base 10.

EOSPAC EOS

Warning

Entropy is not yet available for this 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.

Note for performance reasons this EOS uses a slightly different vector API. See EOSPAC Vector Functions for more details.