Improved Enthalpy-based Equation of State and application to Solid and Porous Lithium Deuteride

1. Introduction

Equation of state (EOS) of materials is essential for the description of shock wave propagation [1]. The EOS defines a relation between pressure P, specific volume V and temperature T in a material in thermodynamic equilibrium [2], and it is generally formulated using V and T as independent variables [3]. The choice of independent variables is purely guided by convenience, in fact, the variables P and T are found more appropriate for describing hydrodynamic phenomena in porous and granular materials. This is mainly due to the ease in treatment of the shock induced compaction phase in these materials; it is particularly useful in dealing regions with multi-valued pressure versus volume variations. There have been significant progress in developing such EOS formulations for porous materials [4,5], which also have firm theoretical basis [6]. The choice of P and T variables is, perhaps, most appropriate for describing shock wave propagation in heterogeneous mixtures, as pressure P and material velocity $\overrightarrow{U}$ relax quite fast to steady values after shock traversal, in comparison to thermal energy and hence temperature T[7]. Procedures to extract parameters of the EOS model (with $P-T$ variables) from Hugoniot data of solid as well as porous metals have also been attempted [8,9].

Authors of this paper pursued an approach for using P and T as independent variables via the enthalpy-based equation of state (EEOS) [10], starting from the thermodynamic definition of the enthalpy-parameter $\chi$. This approach yielded the correct formulation for explicit accounting of contributions of electrons in EEOS [10]. It also yielded simple modeling of shock wave properties in solid and porous binary mixtures [11] and multi-component mixtures [12]. In addition, together with a suitable model for the fluid phase of the material, it was possible to evaluate the Hugoniot of highly porous Cu (porosity up to 10), whose final states are in the expanded volume region [13]. Furthermore, EEOS could be implemented within a solver for Euler equations of fluid flow [14]. This solver, which employs a second order accurate flux-corrected-transport algorithm (FCT) [15], provided pressure versus particle-velocity $(P-U)$ Hugoniot curves in excellent agreement with experimental data for different solid and porous materials. The present paper is mainly concerned with (i) introducing explicit temperature dependence in the enthalpy-parameter, (ii) correct implementation of compaction modeling in the $P-\alpha$ model, (iii) a new set of formulas for thermal pressure and energy of electrons (based on Thomas-Fermi model), and (iv) applications to solid as well as porous lithium-6 deuteride (Li₆D) and natural lithium deuteride (LiD-N). New schemes in the hydrodynamic simulations using the EEOS, to generate $P-U$ curves for solid and porous Li₆D, are also discussed. The last item offers the possibility of prediction of $P-U$ data prior to experiments, as main details of the experimental setup could be implemented in numerical simulations.

The paper is organized as follows. In Section II, the general enthalpy-based EOS, together with temperature dependent $\overline{\chi }$ is briefly discussed. The last feature follows from using $\chi$ versus P variations of different $P-V$ curves and then using a non-linear interpolation formula. This development also shows that there is no need to explicitly introduce electron contributions in the EEOS, as done earlier [10,13]; electronic effects need only to be introduced in the computation of $\chi -P$ curves. Modeling of the compaction phase and the modified $P-\alpha$ model are covered in Section III. Applications and comparisons of the results of the model to LiD-N and Li₆D are given in Section IV. Issues arising in the implementation of EEOS in the FCT solver for obtaining $P-U$ shock wave data (of Li₆D), are covered in section V. Finally, Section VI summaries the work. The basic formulations of zero-temperature isotherm (Vinet EOS supplemented with quantum statistical model (QSM)), ion thermodynamic functions including Debye-Eintein model, melting, rotation and vibration contributions of LiD molecule, and its dissociation are given in the Appendix. Accurate new fits to Thomas-Fermi model functions for electron thermal pressure and energy, are also given there.

2. Enthalpy-Based EOS

Development of EEOS starts with the thermodynamic definitions of enthalpy-parameter $\chi$ and constant-pressure specific heat $C_P$. These are given by ${(\partial V/\partial T)}_P=(\chi /P)C_P$ and ${(\partial H/\partial T)}_P=C_P$, respectively, where $H=E+PV$ and E denote the specific enthalpy and internal energy, respectively. These parameters contain thermal contributions of ions and electrons, and so can be separated as $\chi C_P={\chi }_i{C}_{Pi}+{\chi }_e{C}_{Pe}$, where ${\chi }_i$ and ${\chi }_e$ denote the ionic and electronic enthalpy parameters. Similarly, $C_P=C_{Pi}+C_{Pe}$, where $C_{Pi}$ and $C_{Pe}$ represent the ionic and electronic specific heats, respectively. These separations are not absolutely essential in the formulation as shown below. Now, integration from 0 to T of the basic definitions yield

$$\begin{array}{ccc}\hfill V& =& V_c\left(P\right)+{\displaystyle \frac{1}{P}}{\displaystyle \underset{0}{\overset{T}{\int }}}\chi (P,\tau )C_P(P,\tau )d\tau \hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill H& =& H_c\left(V\right)+{\displaystyle \underset{0}{\overset{T}{\int }}}{C}_P(P,\tau )d\tau .\hfill \end{array}$$

(2)

Here, $H=H_c+H_{th}$ is total enthalpy; the two terms in the sum represent zero-temperature and thermal contributions, respectively. The reference state parameters $V_c\left(P\right)$ and $H_c\left(P\right)$, which refer to $T=0$ and given pressure P, arise as integration constants. Now, Equation (1) is rewritten as

$$\begin{array}{ccc}\hfill V& =& V_c\left(P\right)+{\displaystyle \frac{1}{P}}\overline{\chi }(P,T)(H-H_c)\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill \overline{\chi }(P,T)& =& \left[{\displaystyle \underset{0}{\overset{T}{\int }}}\chi (P,\tau )C_P(P,\tau )d\tau \right][{\displaystyle \underset{0}{\overset{T}{\int }}}{C}_P(P,\tau )d\tau {]}^{-1}.\hfill \end{array}$$

(4)

where $\overline{\chi }(P,T)$ is introduced as an average enthalpy parameter. It should be noted that Equation (3) is exact and, importantly, it shows that $\overline{\chi }(P,T)$ should be temperature dependent in general. However, along a specific $P-V$ curve, T is expressible as a function of P by inverting the EOS: $P=P(V,T)$. Then $\overline{\chi }(P,T)$, along that curve, reduces to a function of P only. It is denoted as ${\chi }_h\left(P\right)$ with the subscript h indicating the specific curve. It should be emphasized that this argument holds only for ${\chi }_h\left(P\right)$ along a specific $P-V$ curve like the isotherm, isentrope, principal Hugoniot, etc. By using the same arguments for ${\chi }_i$, ${\chi }_e$, $C_{Pi}$ and $C_{Pe}$, where the subscripts $iande$ denote ion and electron components, ${\chi }_h\left(P\right)$ is expressible as a weighted average ${\chi }_h\left(P\right)=({\chi }_i{C}_{Pi}+{\chi }_e{C}_{Pe})/C_P$. But it is clear that ${\chi }_h\left(P\right)$ is directly computed if the electron components are added in the functions $P(V,T)$ and $E(V,T)$.

2.1. Specifications for EOS

A three component description, using the zero-temperature isotherm and energy $(P_c,E_c)$, the ionic components of thermal energy and pressure $(P_{ti},E_{ti})$ and the electron components $(P_{te},E_{te})$, is needed in any EOS specification [3,16]. These aspects are briefly discussed below, and in more detail in the Appendix.

The zero-temperature isotherm of Vinet et al [17,18] is known to be quite accurate up to about 1 TPa pressure [19], and is given by $P_v=3B_0{(V/V_{c0})}^{-2}(1-(V/Vc0))exp[(1-(V/V_{c0}))(B_0^{\prime }-1)3/2]$ (see Appendix A.1). The volume $V_{c0}$ is to be determined such that total pressure is ${10}^{-4}$ GPa at ambient temperature $T_0$ (300 K) and volume $V_0$. For very high pressures, energy and pressure must approach those of an electron gas around the nucleus, as given by the quantum statistical model (QSM) [20]. This model accounts for exchange and correlation effects in addition to electron density gradients [21]. Electron pressure in QSM is analytically expressed as $P_q\left(V\right)=({\hslash }^2/2m_e)[(2/5){\left(3{\pi }^2\right)}^{2/3}{N}_s^{5/3}-(2/a_B)(13/36){(3/\pi )}^{1/3}{N}_s^{4/3}]$ (see Appendix A.1). An interpolation procedure [22] to smoothly go over from Vinet’s model to QSM uses a adjustable volume parameter $V_m$ such that the corrected pressure $P_c=E_v$ and $E_c=E_v$ for $V\ge V_m$. These for smaller volumes are obtained by imposing continuity of $P_c\left(V\right)$ and its first two derivatives at $V_m$ (see Appendix A.1). This procedure gives a smooth transition from Vinet’s model to the QSM [23], and is expected to be better than arbitrary modification of $P_v\left(V\right)$ and $E_v\left(V\right)$ at low volumes.

The thermal components of ionic energy and pressure (per gram) are expressed as $E_{ti}(V,T)=E_S+N{k}_{B}T(E_0+E_{\psi }+E_1)$ and $P_{ti}(V,T)=P_S+(1/V)(2{\Gamma }_{th}-2/3)N{k}_{B}T(E_0+E_{\psi }+E_1)$ (see Appendix A.2). Here the original Jonson’s model for metals [24] is generalized so that the solid-energy $E_S=E_d+E_e$ consists of both the Debye ($E_d$ acoustic modes) and Einstein ($E_e$ optical modes) contributions [25]. Similar is the case for pressure $P_S=P_d+P_e$. The components $E_0$ and $E_{\psi }$ (non-zero for $T\ge T_M$, melting temperature) accounts for changes in specific heat upon melting and heat of fusion [16,24]. An additional term $E_1$ ($\ne 0$ for $T\ge T_M$) is added to account for rotation, vibration and dissociation of LiD molecule. Three-term fits [26] to the Grüneisen parameters ${\Gamma }_d$ and ${\Gamma }_e$ are employed in the present model, and specific heat weighting provides the average ${\Gamma }_{th}$ (see Appendix A.3). The adequacy of these models are checked against experimental data on LID-N (natural) in the section dealing with results (see Section 4.1).

It is well known that electron thermal properties are provided by the finite-temperature TF model [27]. A new numerical method to solve the TF equation [28], and thereafter Brachman’s differential equation for electron energy [29,30], has been developed to generate extensive tables for $P_{te}$ and energy $E_{te}$ as a part of the present work (see Appendix A.4). The scaled variables $t=T{Z}^{-4/3}$ and $R{Z}^{1/3}$, where Z and R denote atomic number and cell radius, respectively, varied over ${10}^{-4}\le T{Z}^{-4/3}\le {10}^3$ and $0.9\le R{Z}^{1/3}\le 20$, in these tables. Thereafter, rational function approximations in the form $E_{te}=Z^{7/3}[e_1{t}^2+e_2{t}^3]{[1+e_3{t}^{e_4}+t^2]}^{-1}$ and $P_{te}=Z^{10/3}[p_1{t}^2+p_2{t}^3]{[1+p_3{t}^{p_4}+t^2]}^{-1}$ are derived, which have the correct low and high temperature limiting forms [28]. Dependence of $E_{te}$ and $P_{te}$ on cell radius is incorporated in the fitting parameters $e_n$ and $p_n$, for $1\le n\le 4$ (see Appendix A.4). Unlike the analytical fits available in the literature [31,32], these new expressions are more accurate and also applicable in expanded volume regions.

2.2. Temperature Dependent $\overline{\chi }$

The formulation of thermodynamic functions, outlined above, is appropriate for computing the parameters ${\chi }_h$ and $C_P$ along any arbitrary $P-V$ curve. As the expressions quoted are derivable from respective free energies, they provide consistent thermodynamic results. The constant-volume specific heat $C_V={(\partial E/\partial T)}_V$, Grüneisen parameter $\Gamma =V{(\partial P/\partial T)}_V/C_V$, and isothermal bulk modulus $K_T=-V{(\partial P/\partial V)}_T$ are readily determined on the specified $P-V$ curve. Temperature T on the curve is obtained by inverting the equation $P=P(V,T)$ for given values of V and P. This inversion procedure for T is avoidable for the case of Hugoniot or isentrope, as T is determined by solving the temperature-equation [13] in an iterative manner (see Section 4 and Equation (10)). Then, the thermodynamic relations $C_P/C_V$ = $K_S/K_T$ = $1+(T/V){\Gamma }^2{C}_V/K_T$ provide $C_P$ and isentropic bulk modulus $K_S$. It is important to emphasis that all these quantities are obtained explicitly in terms of, say, V and T on the given $P-V$ curve, and the approximation that $\Gamma$ depends only on V is unnecessary. Finally, the defining relation, ${\chi }_h=P\Gamma /K_S$, provides the enthalpy parameter. Note that it also yields the equation $K_S/K_T$ = ${[1-\left(C_{V}T\right)(\Gamma /V)({\chi }_h/P)]}^{-1}$, which is useful whenever ${\chi }_h$ is available.

A simple scheme to generate solid Hugoniot is to employ the parameters C and S in the linear relation $U_s=C+S{U}_p$ between shock speed $U_s$ and particle speed $U_p$. Then the Rankine-Hugoniot jump conditions [1] provide volume $V\left(P\right)$ on the Hugoniot explicitly as [14]

$$V_h=V_0-V_0{\displaystyle \frac{\beta \left(P\right)}{S}}[1-{(1-{\beta }^{-2})}^{1/2}];\beta \left(P\right)=1+{\displaystyle \frac{C^2}{2S{V}_0(P-P_0)}}$$

(5)

Here $P_0$ and $V_0={\rho }_0^{-1}$ are the initial pressure, volume and density. The parameter C is same as adiabatic sound speed ( bulk modulus $B_0={\rho }_0{C}^2$) while S is related to the pressure derivative of bulk modulus as $B_0^{\prime }=(4S-1)$. This equation is used to generate three $P-V$ and $T-P$ curves with the fixed parameters ${\rho }_0=0.8$ g/cm³, $P_0={10}^{-4}$ GPa, $C=0.66$ cm/s, which are appropriate to Li₆D. The three curves in Figure 1 are obtained with the values $S=1.03$ (curve-1), $S=1.31$ (curve-2) and $S=1.51$ (curve-3). The other two (curve-0 and curve-H) correspond to zero-temperature isotherm and principal Hugoniot (S=1.11). The $V-P$ and $T-P$ curves corresponding the four cases (curves 1, 2, 3 and H) are also shown in the insets in this figure. The values of ${\chi }_h$ versus P obtained using the EOS model outlined above (see Appendix also) are shown in the main figure. Numerical fits of the ${\chi }_h$ data to the expression ${\xi }_h=(w_{1}P+w_2{P}^{3/2}+{\chi }_{\infty }{P}^2){(1+w_{3}P+w_4{P}^{3/2}+P^2)}^{-1}$ are shown as continuous lines in the figure. The fit constants $w_1,w_2,\dots w_4$ are provided in Table 1 for all the five curves. The term ${\chi }_{\infty }=2/5$ in this expression is the limiting value of ${\chi }_h$ corresponding to mono-atomic ideal gas $[\Gamma /(1+\Gamma \left)\right]$ with $\Gamma =2/3$. The fitting function is chosen so as to yield this limit for high pressures. The temperature data is fitted to the form ${\tau }_h=(T_0+t_{1}P+t_2{P}^2+t_3{P}^3){(1+t_4{P}^2)}^{-1}$ ($T_0=300$ K is ambient temperature). Table 1 also contains values of these fitting constants.

The formulas for ${\xi }_h\left(P\right)$ and ${\tau }_h\left(P\right)$ are used to develop an interpolation scheme for computing $\overline{\chi }(P,T)$. First of all, note that in the vicinity of any of the five curves, the expression ${\chi }_h(P,T)=\{{\xi }_0+{\chi }_{\infty }{\beta }_h\left(P\right)(T/{\tau }_h)\}/\{1+{\beta }_h\left(P\right)(T/{\tau }_h)\}$, where ${\xi }_0$ is its value at $T=0$, provides a temperature variation. This function tends to the asymptotic limit ${\chi }_{\infty }$ as $T\to \infty$. The choice of the function ${\beta }_h\left(P\right)=({\chi }_{\infty }-{\xi }_h)/({\xi }_h-{\xi }_0)$ ensures that ${\chi }_h(P,T)$ reduces to ${\xi }_h\left(P\right)$ on the curve-h. Thus the interpolation tends to ${\chi }_{\infty }$ as $P\to \infty$ as well. Now the function $T/{\tau }_j$ is generalized to the Lagrange polynomial, viz. ${\theta }_j\left(T\right)=(T/{\tau }_j){\Pi }_{k\ne j}(T-{\tau }_k)/({\tau }_j-{\tau }_k)$ and $\overline{\chi }(P,T)$ is expressed as

$$\overline{\chi }(P,T)={\displaystyle \frac{{\xi }_0+{\chi }_{\infty }{\sum }_{1\le j\le H}{\beta }_j\left(P\right){\theta }_j\left(T\right)}{1+{\sum }_{1\le j\le H}{\beta }_j\left(P\right){\theta }_j\left(T\right)}}.$$

(6)

On any of the five curves, say curve-h, the parameter ${\theta }_j$ reduces to ${\delta }_{hj}$ and hence $\overline{\chi }(P,T)$ would yield the fitted function ${\xi }_h$.

3. Modeling Compaction Phase

A parameter that is used to characterize porous materials is its initial porosity, defined as ${\alpha }_0=V_{00}/V_0$, where $V_{00}$ and $V_0$ are, respectively, the initial specific volumes of the porous and the solid materials. Carroll and Holt [33] considered the compaction of a spherical pore, surrounded by an in-compressible matrix material, due to an applied time-dependent pressure P on the outer surface of the matrix. Time dependence of P is characterized by the linear ramp rate $\dot{P}$. These authors derived a formula for the variation of porosity function with pressure, and it is expressed as $\alpha \left(P\right)={\alpha }_0$ for $0\le P\le P_{cri}$ and $\alpha \left(P\right)=(1-exp[-3P/2Y]{)}^{-1}$ for $P_{cri}<P<\infty$. Here, Y is the strength of the matrix material and $P_{cri}=(2/3)Yln({\alpha }_0/({\alpha }_0-1))$ is the elastic critical pressure of the porous material. Compaction process is initiated when P exceeds $P_{cri}$. This formula agrees well with numerical results of the collapse model. Their analysis also explained the validity of the relation $V=\alpha V_m$, where V and $V_m$ represents the specific volumes of the porous and matrix materials. However, it is found that there are some differences between numerical results for pressure P in the porous material and those provided by the $P-\alpha$ model formula, viz. $P=P_m/\alpha$, for $\dot{P}>0.01GPa/ns$. This relation is also expressed as $PV=P_m{V}_m$, where $PV$ and $P_m{V}_m$ represent thermal energy in the porous and matrix materials, respectively. Thus the range of validity of the relation $V=\alpha V_m$, $\left(withrespecttothevaluesof\dot{P}\right)$, is found to be much more than that of $P=P_m/\alpha$. Another assumption in the $P-\alpha$ model is that the zero-temperature energy $E_c\left(V_0\right)$ is unaltered due to the presence of pores, but pressure increases during the compaction process. As the solid matrix is assumed to be in-compressible during compaction, $V_m$ is the same as $V_0$. So, a simple way to extend Equation (1) to porous materials is to compute $\alpha$ from Carol-Holt relation for a specified value of P and obtain the zero-temperature volume for porous material as $V_c^{\prime }=\alpha V_c$. A more accurate approach is discussed next.

3.1. Low Pressure Compaction

First step in this detailed approach to model compaction process is to rewrite Equation (1) for the matrix material as

$$V_m=V_c\left(P_m\right)+{\displaystyle \frac{\overline{\chi }(P_m,T)}{1-\overline{\chi }(P_m,T)}}{\displaystyle \frac{1}{P_m}}(E_m-E_c\left(P_m\right)).$$

(7)

Here the definitions $H=E+PV$ and $H_c=E_c+P{V}_c$ are used. The subscript m in $V_m$, $E_m$ and $P_m$ denotes that these variables pertain to the solid matrix. Note that Equation (7) is the $P-V-T$ isotherm of the solid matrix. Now, using the relations $P=P_m/\alpha$, $E=E_m$ and $V=\alpha V_m$ yields the equation

$$V=\alpha V_c\left(\alpha P\right)+{\displaystyle \frac{\overline{\chi }(\alpha P,T)}{1-\overline{\chi }(\alpha P,T)}}{\displaystyle \frac{1}{P}}(E-E_c\left(\alpha P\right))$$

(8)

Starting with a specified value of P, and Carol-Holt formula for $\alpha \left(P\right)$, Equation (8) provides the $P-V-T$ isotherm of the porous material. All the aspects of the $P-\alpha$ model are present in this equation. A similar equation, however, using an approximate Mie-Grüneisen relation (which assumes $\Gamma V$ is a constant), has been derived elsewhere [34]. But Equation (8) is more general because accurate expressions for $\overline{\chi }(P_m,T)$ and $E_c\left(V_m\right)$ of the solid matrix can be employed. The advantage of using P and T as independent variables is most evident in obtaining Equation (8) for porous materials.

The compaction curve along the Hugoniot follows by substituting the Hugoniot equation for energy, $E-E_0=(1/2)(P+P_0)(V_{00}-V)$ into Equation (8) (see also Equation (9) in Section 4). Results obtained in this manner at 300 K and experimental data [35] on porous Cu are compared in Figure 2. Excellent agreement obtained establishes the accuracy of Equation (8). Since temperature involved is very low, ${\chi }_0$(P) along the zero-temperature isotherm of Vinet et al is used for the computations. It is verified that the isotherm of Li et al [36] also produces indistinguishable results. All the data needed for Cu are tabulated in the Appendix (see Table A1 and Table A2) .

3.2. Menikoff’s Modified $P-\alpha$ Model

The assumption of $P-\alpha$ model that internal energy does not change during compaction is thermodynamically inconsistent. The modification proposed by Menikoff [37] starts with the Helmholtz free energy model of the porous material $F(V,T,\varphi )=F_m(V_m,T)+B\left(\varphi \right)$. Here $F_m(V_m,T)=F_m(\varphi V,T)$ is the free energy of the solid matrix and $\varphi =V_m/V={\alpha }^{-1}$ is its volume fraction at an intermediate compaction; ${\varphi }_0=V_m/V_{00}={\alpha }_0^{-1}$ being the initial volume fraction. The term $B\left(\varphi \right)$, arising purely due to the pores, in the expression for free energy may be related to the surface energy $E_{sur}\left(\varphi \right)$ of pores as $B\left(\varphi \right)=E_{sur}\left({\varphi }_0\right)-E_{sur}\left(\varphi \right)$, so that $B\left(\varphi \right)$ is an increasing function of $\varphi$ and it vanishes at ${\varphi }_0$. In general, $\varphi$ is an independent variable characterizing the porous state, and hence pressure in the porous material is $P=-{[\partial F(V,T,\varphi )/\partial V]}_{T,\varphi }=\varphi P_m$, which is the basic equation of $P-\alpha$ model. Internal energy of the porous material is $E=E_m(V_m,T)+B\left(\varphi \right)$, as entropy ${[-\partial F(V,T,\varphi )/\partial T]}_{V,\varphi }$ is the same as that of the solid matrix.

Now, if the material relaxes very fast to equilibrium conditions after attaining pressure P and temperature T, then the volume fraction also corresponds to equilibrium conditions, and can be denoted as ${\varphi }_e$. It should be determined from the condition ${[\partial F/\partial \varphi ]}_e=0$ (i.e., ${\varphi }_e$ should correspond to minimum of F). This yields the relation $P_m{V}_m={[\varphi \partial B\left(\varphi \right)/\partial \varphi ]}_e$, which implies a relationship of the form ${\varphi }_e=\psi \left(P_m{V}_m\right)$ at equilibrium conditions, where $\psi$ is a suitable, usually empirical, function. Note that in Carroll-Holt model, $\varphi$ (or $\alpha$) corresponds to a time-dependent pressure P, and is dependent only on $P_m$. However, volume at equilibrium, $V_m$ (in $\psi \left(P_m{V}_m\right))$ is determined in terms of $P_m$ (and T) via the EOS of the solid matrix. The differential equation $P_m{V}_m={\psi }^{-1}\left({\varphi }_e\right)={\varphi }_e\partial B\left({\varphi }_e\right)/\partial {\varphi }_e$ is integrated starting form ${\varphi }_0$ to obtain $B\left({\varphi }_e\right)$. The Carroll-Holt compaction function is rewritten as $P=P_{cri}+(2Y/3)ln[(1-{\varphi }_0)/(1-\varphi )]$ for $P\ge P_{cri}$. The assumption that this is also valid in equilibrium condition, yields the differential equation $\partial B\left({\varphi }_e\right)/\partial {\varphi }_e=P_{cri}{V}_m/{\varphi }_e^2+(2V_{m}Y/3)(1/{\varphi }_e^2)ln[(1-{\varphi }_0)/(1-{\varphi }_e)]$. As the matrix material is in-compressible during compaction [33], $V_m=V_0$. This reduces the equation for $B\left({\varphi }_e\right)$ purely in terms of matrix material constants. Thus the improved $P-\alpha$ model also provides changes in internal energy during compaction, thereby removing the inconsistency. The differential equation is integrated along the low pressure compaction curve of Cu considered earlier (see Section 3.1). The variation of $B\left({\varphi }_e\right)$ versus P (as ${\varphi }_e$ is related to $P_m$ and hence P) is shown in in the inset figure of Figure 2. Energy on the Hugoniot at 2.5 GPa is found to be $E_h=-5170(J/g)$, so the value of $B=25(J/g)$, seen from the inset figure at the same pressure, is negligible. Thus, the original $P-\alpha$ model (within its range of validity) is adequate for all practical applications.

4. Application to Lithium Deuteride

Substitution of the expression for energy along the Hugoniot, $E=E_0+(1/2)(P+P_0)(V_{00}-V)$, into Equation (8) yields the volume-pressure curve (on the Hugoniot)

$$V={\displaystyle \frac{1-{\overline{\chi }}_m}{1-{\overline{\chi }}_m+{\overline{\chi }}_m^{*}}}\alpha \left(P\right)V_{cm}+{\displaystyle \frac{{\overline{\chi }}_m}{1-{\overline{\chi }}_m+{\overline{\chi }}_m^{*}}}\left(\right(V_{00}^{*}+{\displaystyle \frac{1}{P\alpha \left(P\right)}}(E_{00}-E_c\left(V_{cm}\right)))$$

(9)

where ${\overline{\chi }}_m^{*}=({\overline{\chi }}_m/2\alpha )(1+P_0/P)$. Similarly, $V_{00}^{*}=(V_{00}/2\alpha )(1+P_0/P)$. The subscript m on ${\overline{\chi }}_m$ and $V_{cm}$ indicates that these quantities should be computed at pressure $P_m=P\alpha \left(P\right)$ in the matrix material. Thus, note that ${\overline{\chi }}_m^{*}$ and $V_{00}^{*}$ contain the parameter $\alpha \left(P\right)$ in their present definitions. In the equations used by earlier authors [4,5,13], the parameter $\alpha \left(P\right)$ appears only as $\left[\alpha \left(P\right)V_c\left(P\right)\right]$. Thus, quantities like $V_c$, $\overline{\chi }$, $E_c$ etc. occurring in Equation (9) were computed at P in those papers, and not at $P_m$. Of course, the differences in above equation are significant only for higher porosity, like ${\alpha }_0\ge 2$, and in the compaction phase. Nevertheless, it should be noted that Equations (8) and (9) follow from correct application of the $P-\alpha$ theory. It is not straight forward to implement this model when V and T are treated as independent variables. For high values of P, Equation (9) shows that $V/V_{00}\to 0.5{\overline{\chi }}_m/(1-0.5{\overline{\chi }}_m)$ as the parameters $\alpha \to 1$ and $V_{cm}\to 0$. Note that this is just the mono-atomic ideal gas limit as ${\overline{\chi }}_m\to 2/5$.

Volume in Equation (9) is determined together with an equation for temperature T on the Hugoniot, as ${\overline{\chi }}_m$ depends explicitly on T. The differential equation for T is [14]

$${\displaystyle \frac{d}{dP}}T-{\displaystyle \frac{{\overline{\chi }}_m}{P}}T\left(P\right)={\displaystyle \frac{0.5}{C_{Pm}}}((V_{00}-V)+(P-P_{00}){\displaystyle \frac{dV}{dP}})$$

(10)

Here $C_{Pm}$ is the constant-pressure specific heat, at pressure $P_m$, of the matrix material. Using a forward finite difference scheme, a discrete form of Equation (10) is obtained as

$$(1-\Delta P{\displaystyle \frac{{\overline{\chi }}_{mk}}{P_k}})T_k=T_{k-1}+{\displaystyle \frac{0.5}{C_{Pmk}}}((V_{00}-V_k)\Delta P+(P_k-P_{00})(V_k-V_{k-1}))$$

(11)

Here the subscript k on different variables denotes their respective values at $P_k=P_{00}+k\Delta P$ where $\Delta P$ is an increment in pressure, and $k=1,2,3,\dots$. Finally, note that Equations (9) and (11) are evaluated iteratively at each $P_k$ to obtain $V_k$. Starting with a guess value of $T_k$, say $T_{k-1}$, one iteration is complete on obtaining $V_k$, using Equation (9), and then updating $T_k$, via Equation (11). These steps are repeated until values of $T_k$ converge within a prescribed error (say, ${10}^{-3}$). Centered values of parameters in Equation (11) are introduced after first iteration, which converge in few (3-4) steps. Completion of the whole process for $k=1,2,3,\dots$, provides volume and temperature along the Hugoniot in a consistent manner. The method is applied to (natural) LiD and ${Li}_{6}D$ and compared with experimental data below.

4.1. Isotherms of Natural LiD

Before considering shock wave data on (natural) LiD-N, the ion model and zero-temperature isotherms (discussed in Appendix) are evaluated using experimental data. Employing in situ high-pressure and high-temperature neutron diffraction techniques (for pressures up to 4.1 GPa and temperatures 1280 K), thermodynamic properties of LiD-N (${\rho }_0=0.8837{\mathrm{gm}/\mathrm{cm}}^3$) are available in the literature [38]. Analysis of the pressure-volume data at 300 K, using the finite-strain EOS $P=3f{(1+2f)}^{5/2}{B}_0(1-2\zeta f)$, where $\zeta =(3/4)(4-B_0^{\prime })$ and $f=(1/2)[{(V_0/V)}^{2/3}-1]$, authors of this work deduced $B_0=31.5$ GPa assuming that with $B_0^{\prime }=4$. Similar analysis of the data at 500 K shows that $B_0=26.1$ GPa [38]. These data are re-interpreted in the present paper using the zero-temperature isotherm of Vinet et al [17,18] (see Appendix A.1). It is known [39] that this is an excellent model for LiD-N (see the inset in Figure 3). Einstein-Debye model, which takes into account the optical and acoustic branches of lattice vibrations (see Appendix A.2), is used as the thermal model. The results so obtained for the two isotherms of LiD-N are compared with experimental data [38] in Figure 3. Solid lines are from present computations while symbols denote experimental data. Agreement to the data at 500 K needs explicit accounting of optical branch while the acoustic branch alone is found adequate for 300 K data. The parameters used for all calculations are listed in the Appendix (see Table A1 and Table A2).

4.2. Principal Hugoniot of natural LiD

A series of shock wave experiments driven by impact with velocities in the range $17-32$ km/s are reported on single crystals of LiD-N with initial density around ${\rho }_0=0.886g/c{m}^3$ [40]. Pressure versus volume on the principal Hugoniot (up to 570 GPa) and behind re-shock states (up to 920 GPa) are available for comparison with results of EEOS formalism. These results significantly extend those of earlier LASL shock wave experimental data (up to 45 GPa) (see Section 4.3) available in the literature [41]. Results of computations are shown in Figure 4, and compared with both sets of data (filled circles and filled stars for high and low pressures, respectively). Excellent comparison shows the accuracy of the EEOS model up to pressures of about 600 GPa. The agreement between computed results and data is similar to that obtained via ab-initio molecular dynamic (AIMD) calculations [40]. The inset graph in Figure 4 shows computed temperatures on the principal Hugoniot versus pressure and experimental data (filled circles). To account for the contribution of dissociation of LiD molecule to ionic energy, a term of the form $E_v=(3/2)(T_v/T){(Exp[T_v/T]-1)}^{-1}$, with $T_v=2500°K$, is added so that the total ionic energy is $3K_{B}T$ (per LiD molecule) in the high temperature limit (see Appendix A.2). Although there are only four data points, the agreement with the results of EEOS is satisfactory.

Next, the results of EEOS model are compared with experimental data [40] on re-shock Hugoniot of LiD-N. The same formula in Equation (9), with initial conditions on the principal Hugoniot, generates the re-shock states. The first-shock and re-shock states (filled circles), shown in Figure 5, agree well with computed results (continuous lines). This comparison further affirms the accuracy of the model up to about 900 GPa. Note that higher densities (and lower temperatures), in comparison to that on the principal Hugoniot, are obtained on the re-shock states.

4.3. Solid and Porous ${Li}_{6}D$ Hugoniot

Extensive compilation of results of shock wave experiments from LASL, for a variety of solid and porous materials, are available [41]. These data are primarily based on measurements of fluid velocity $U_p$, shock velocity $U_s$ and free surface velocity $U_{fs}$. Many materials follow a linear relation $U_s=C+S{U}_p$ beyond the elastic pressure regime (see Table A1). Tables of P versus V, deduced from measurements, in the data files are directly comparable with computed results shown in Figure 6. Data considered for comparison correspond to solid LiD-N , solid ${Li}_{6}D$ and two porous cases of ${Li}_{6}D$ (with initial porosity ${\alpha }_0=1.212$ and $1.781$, respectively) (see the legends in the figure). Successive graphs in this figure are shifted upwards by 10 units for clarity. There is good agreement between computed results and the data. Similar comparison is also found for other porosity values ${\alpha }_0=$ 1.052, 1.081, 1.379 and 1.568.

The increase in pressure P due to increase in the initial volume $V_0$ is clearly evident from these results. The Hugoniot relations [1], which imply $P=(U_s^2/V_0)(1-V/V_0)$, show that P increases with $V_0$ for the same shock wave energy density $U_s^2/V_0$ and final volume V. This feature is also clear from the results of solid LiD-N and ${Li}_{6}D$. But, shock compression of porous materials involves initial compaction, leading to pore closure, followed by compression of the matrix material [34]. The microscopic mechanisms involved during the process of pore collapse are: (i) initial shock loading of solid matrix, (ii) shock break out at the pore interfaces, (iii) collision energy dissipation of the ejecta, and (iv) consequent local heating. Thus, higher temperatures and pressures are produced on shock loading of porous materials [10]. Furthermore, there is anomalous expansion in the shocked state, for higher initial porosity. For instance, when initial porosity ${\alpha }_0$ of Cu is more than 2, the final volume increases with pressure, instead of decreasing. Then P taken as a function of V becomes double valued. In this situation, the choice of P and T as independent variables in thermodynamic modeling, and hence EEOS, becomes preferable.

Next, Hugoniot of ${Li}_{6}D$ at extreme pressures (up to $\sim {10}^6$ GPa) is considered, where multi-valued pressure versus volume arises. Results of extensive computations of the principal Hugoniot of ${Li}_{6}D$ up to pressures about $2.5\times {10}^5$ GPa is reported [42]. These authors used the Khon-Sham density-functional molecular dynamics (KSMD) for temperatures in the range 0.5 to 25 eV, and extended it with orbital-free density-functional molecular dynamics (OFMD) up to about ${10}^3$ eV. The OFMD scheme is based on Thomas-Fermi-Dirac model. Both methods are matched in the intermediate range of 10 to 20 eV. Detailed results for pressure P and temperature T versus density $\rho$ on the principal Hugoniot of the solid (${\rho }_0=0.8{g/cm}^3$), and P versus $\rho$ data for 4 porous-cases (${\rho }_0=0.71,0.67,0.58and0.45{g/cm}^3$) are provided.

The results computed using EEOS together with data obtained via KSMD (filled stars for P > 100 GPa) and OFMD (filled circles) for solid ${Li}_{6}D$ are shown in Figure 7. LASL experimental data [41] up to 45 GPa (filled stars) are also shown. Initial build up followed by exponential increase (note the Log-scale) and thereafter bending to lower densities are the characteristic features of the $P-\rho$ curve. The bending, and consequent double-valued pressures, is due to ions and electrons in ${Li}_{6}D$ approaching the behavior of ideal fluids. Maximum compression $(\rho /{\rho }_0)$ achieved is just $4.271$ which may be compared with 4 for mono-atomic ideal gas. This increase arises from absorption of energy from the shock for dissociation of ${Li}_{6}D$ molecule and thermal ionization. At the extreme limit of pressure, $P\sim 5\times {10}^5$ GPa, the shocked state would contain ions of Li₆, D and electrons. Almost all features of the complete Hugoniot are brought out by the EEOS model. Temperature T along the Hugoniot is compared with KSMD (filled stars) and OFMD (filled circles) data in the inset figure. Although there are some differences, main features of $T-\rho$ curve, particularly temperature $5\times {10}^5$K reached at maximum compression $\rho /{\rho }_0=4.271$, and bending of the curve are reproduced by EEOS model.

Finally, the results based on OFMD scheme for solid and four porous cases are compared with those of EEOS model in Figure 8. The comparison is similar to that of pure solid case. As intial density decreases from ${\rho }_0=0.8{g/cm}^3$, maximum density reached decreases, although compression factor $(\rho /{\rho }_0)$ remain the same. For the case of ${\rho }_0=0.45{g/cm}^3$ maximum density reached is ${\rho }_0=1.879{g/cm}^3$.

5. Implementation of EEOS in Euler Solver

The Euler equations describing the conservation laws of fluid flow are non-linear partial differential equations, in $\overrightarrow{r}$ and time t, for density $\rho$, fluid velocity $\overrightarrow{u}$ and total energy $E_t=E+0.5\times \overrightarrow{u}·\overrightarrow{u}$. These are expressible in a fixed (Eulerian) reference frame or a moving (Lagrangian) frame [1]. For the cases of all one-dimensional geometries (planar, cylindrical and spherical) the integral form of general conservation law is written as:

$${\displaystyle \frac{\partial }{\partial t}}{\displaystyle \underset{V\left(t\right)}{\int }}GdV+{\displaystyle \underset{\partial V\left(t\right)}{\int }}G(u-u^g)·d\overrightarrow{A}+{\displaystyle \underset{V\left(t\right)}{\int }}[{\displaystyle \frac{1}{r^{n-1}}}{\displaystyle \frac{\partial }{\partial r}}\left(r^{n-1}B\right)-C{\displaystyle \frac{\partial D}{\partial r}}]dV-{\displaystyle \underset{V\left(t\right)}{\int }}QdV=0.$$

(12)

Here, $V\left(t\right)$ is the time-dependent volume element (arising due to moving grid in Lagrangian frame) and $\partial V\left(t\right)$ its surface element. The spatial coordinate is r, and u and $u^g$ represent the fluid velocity and grid velocity along the r-coordinate. The grid velocity $u^g$ assumes values 0 and u in Eulerian and Lagrangian schemes, respectively. Volume and surface elements in one-dimension are: $dV=2^{n-1}\pi r^{n-1}dr$ and $d|\overrightarrow{A}|=2^{n-1}\pi r^{n-1}$, respectively. The index n takes values 1, 2 and 3 in planar, cylindrical, and spherical geometry, respectively. Here, G stands for $\rho$, $\rho u$ or $\rho E_t$ corresponding to the three Euler equations for mass, momentum and total energy [1]. The internal source function (third term in Equation (12)) is explicitly written in terms of coefficients B, C and D. Note that $B=Pu$ for the energy equation, but zero for mass and momentum equations. Similarly, $C=-1$ and $D=P$ for the momentum equation, but these are zero for the other two cases. These and source rate Q are functions of r and t. An EOS which provides pressure in either form $P=P(\rho ,E)$ or $V=V(P,E)$ is required to close the system of equations. It is sufficient to have an algorithm to treat Equation (12) for solving Euler equations in the order $\rho \to \rho u\to \rho E_t$. In fact, this form is much more general and a large class of problems is mapped into it by defining the coefficients appropriately.

The one dimensional FCT algorithm, developed by Boris, Book and co-workers [15], solves a system of conservation equations. This algorithm consists of a convective diffusive stage which is followed by an anti-diffusive stage. The fluid variables are updated using a three point spatial difference scheme during the former stage in every time step. Some amount of diffusive transport is allowed in this stage so that numerical stability and positivity are preserved. However, the latter stage of the algorithm involves addition of anti-diffusion fluxes for reducing the smearing effects of diffusion. The important aspect of the method is that the anti-diffusion fluxes are corrected such that no new overshoots in fluid variables are generated. Integration in time is improved using a two-step scheme wherein time-centered source terms are generated before the full time-step integration. The FCT algorithm is coded as FORTRAN subroutines [15] which are easily used with a driver program and an EOS package. The algorithm has the flexibility that either the Eulerian scheme or moving Lagrangian scheme is used in the simulation. The method is also discussed in detail, and applied to simulate several strong-shock problem [43], and the results so obtained are found to be quite satisfactory. Next, the updating equations for P and T, while using EEOS, are discussed.

5.1. Updating P and T

After every hydrodynamic sweep (in time), spatial profiles of P and T in the meshes are to be updated. The information available from the FCT algorithm for this purpose are the increments $\Delta \rho$ or $\Delta V$ and $\Delta E$ at all meshes. Assumption of local thermodynamic equilibrium yields the relations

$$\begin{array}{ccc}\hfill \Delta V& =& -{\displaystyle \frac{V}{K_T}}\Delta P+C_P{\displaystyle \frac{\overline{\chi }}{P}}\Delta T,\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill T\Delta S& =& -T{C}_P{\displaystyle \frac{\overline{\chi }}{P}}\Delta P+C_P\Delta T.\hfill \end{array}$$

(14)

Here, Equation (13) follows by considering $V=V(P,T)$ (as in EEOS) and the definition $K_T=-V{(\partial P/\partial V)}_T$ of iso-thermal bulk modulus. A basic relation involving Grüneisen parameter $\Gamma$, $C_V$ and $K_T$, viz. $\Gamma C_V/K_T=\overline{\chi }{C}_P/P={(\partial V/\partial T)}_P$, is also employed. Similarly, Equation () arises by taking $S=S(P,T)$ and using the definition $C_P=T{(\partial S/\partial T)}_P={(\partial H/\partial T)}_P$. This equation can be rewritten using the `TdS’ relation $T\Delta S=dE+P\Delta V$ (which assumes energy conservation and reversible changes), and also Equation (13). The resulting $(2\times 2)$ matrix equation, which provides $\Delta P$ and $\Delta T$ is given by

$$\begin{array}{c}\hfill \left(\begin{array}{cc}-V/K_T\hfill & C_P\overline{\chi }/P\hfill \\ PV/K_T-T{C}_P\overline{\chi }/P\hfill & C_P(1-\overline{\chi })\hfill \end{array}\right)\left(\begin{array}{c}\Delta P\\ \Delta T\end{array}\right)=\left(\begin{array}{c}\Delta V\\ \Delta E\end{array}\right)\end{array}$$

(15)

The vector on the RHS comes from the hydrodynamic cycle. All the terms (like $K_T$, $\overline{\chi }$, $C_P$, etc.) occurring in the matrix elements are readily computed using the variables $V_1$, $T_1$, $P_1$ available at the beginning of the cycle. The EEOS model provides $V_1$ in term of the known values of $P_1$ and $T_1$. After solving the $(2\times 2)$ system, $P_2=P_1+\Delta P$ and $T_2=T_1+\Delta T$ provide estimates of P and T at the end of the cycle. Another iteration, after updating the matrix elements as averages of beginning and end values, is done to obtain a time-centered algorithm.

A root-finding method [14] to determine $P_2$ and $T_2$ is feasible by rewriting Equations (8) and (13) as

$$\begin{array}{ccc}\hfill V_2& =& \alpha V_c\left(\alpha P\right)+{\displaystyle \frac{{\overline{\chi }}^{*}}{P}}(E_2-E_c\left(\alpha P\right)),\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill \Delta V& =& C_P{\displaystyle \frac{\overline{\chi }}{\overline{P}}}(T_2-T_1)-{\displaystyle \frac{\overline{V}}{K_T}}(P_2-P_1).\hfill \end{array}$$

(17)

Here $V_2=V_1+\Delta V$ and $E_2=E_1+\Delta E$ are provided from hydrodynamic cycle. Further, ${\overline{\chi }}^{*}=\overline{\chi }(\alpha P,T)/(1-\overline{\chi }(\alpha P,T))$, and $\overline{P}$ and $\overline{V}$ denote averages of beginning and end values. Starting with $T_1$, first Equation (16) is solved via root-finding method to get $P_2$. Thereafter Equation () directly provides $T_2$. A second iteration is done to redefine other parameters $(\overline{\chi },C_P,K_T,etc.)$ as centered values, as explained earlier.

5.2. Shock Pressure Versus Fluid Velocity in ${Li}_{6}D$

As an application of EEOS package in FCT solver, a configuration of ${Li}_{6}D$ slabs, similar to that used in impact experiments, is considered next. In impact experiments, a projectile (also called impactor) impacts on a target plate and launches two shock waves, one into the plate and the other back to the projectile. After an initial acceleration, fluid particles in the target attain a constant velocity, called fluid (or particle) velocity $U_p$. The shock wave moves in the target with its own velocity denoted as shock velocity $U_s$. Measurement of $U_p$ and $U_s$ in a series of impact experiments provide the principal Hugoniot of the target material.

Simulations in ${Li}_{6}D$ used a configuration consisting of 5 mm thick target, impacted with 3 mm thick impactor. Impact velocities used are up to about 12 km/s. In the first set of simulations, both the impactor and target are chosen as solid ${Li}_{6}D$ (density ${\rho }_0=0.8g/cc$), because inelastic collision between them (same material) yields a fluid velocity which is half of the impact velocity. At time $t=0$, all the mesh-interfaces in the impactor are assigned the impact velocity. Meshes of width $0.01mm$ and time steps of ${10}^{-5}\mu s$ are used in each calculation. The resulting hydrodynamic motion in the target is followed up to about 100 $ns$, thereby yielding steady values of pressure and fluid velocity in a single simulation. A series of such simulations provide the P versus $U_p$ Hugoniot shown in Figure 9 (curve-1), which is compared with the LASL experimental data (symbols) [41]. Excellent agreement obtained with experimental data demonstrates the adequacy of EEOS and its integration with the FCT code. The same procedure is then repeated with targets of two initial porosity specified via ${\alpha }_0$. The results obtained (curve-2) and (curve-3), respectively, for ${\alpha }_0=1.211$ and ${\alpha }_0=1.781$, are also compared with experimental data in Figure 9. Good agreement is found in these cases as well. It should be noted that there is significant attenuation of the shock wave in porous targets, and so P and $U_p$ values have to be picked up, from the simulation data, from meshes close to the impactor-target interface.

6. Summary

The main aim in this paper is to provide a consolidated account of the different aspects of the generalized enthalpy-based EOS. Utility of this scheme (which uses $PandT$ as independent variables) over the conventional approach (with $VandT$ as independent variables), in dealing with multi-valued P versus V variations is brought out. Several new ideas such as: (i) an average temperature-dependent enthalpy parameter $\overline{\chi }$ (see Equation ()), (ii) explicit interpolation formula for $\overline{\chi }$ (see Equation (6)), (iii) correct use of the $P-\alpha$ model for treating compaction (see Equation (8)), and (iv) updating P and T in the EEOS implementation in the FCT solver for Euler equations (see Section 5.1), are discussed. Results of applications to solid (natural) LiD, and ${Li}_{6}D$ with porosity up to ${\alpha }_0=1.781$, are provided (see Section 4.1, Section 4.2, Section 4.3). The results on the principal Hugoniot of ${Li}_{6}D$ extend up to extreme pressures $\sim 2.5\times {10}^5$ GPa. These results of EEOS show good agreement with experimental data (up to ∼ 900 GPa), and also computed data using ab initio molecular dynamics (up to ∼$2.5\times {10}^5$ GPa ). New fits for computing electron-properties, together with brief outlines of zero temperature and ionic EOS, are presented in the Appendix.

Appendix A.1. Zero Temperature Isotherm

The zero temperature isotherm (for pressure, energy and bulk modulus) of Vinet et al [17,18] is

$$\begin{array}{ccc}\hfill P_v& =& 3B_0{X}^{-2}(1-X)exp\left(\eta (1-X)\right)\hfill \\ \hfill E_v& =& -E_b+(9/{\eta }^2)B_0{V}_0[1-Yexp\left(\eta (1-X)\right)]\hfill \\ \hfill B_v& =& B_0{X}^{-2}(1+(\eta X+1)(1-X))exp\left(\eta (1-X)\right)\hfill \end{array}$$

(A1)

where $X={(V/V_{c0})}^{1/3}$. Further, $Y=1-\eta (1-X)$ and $\eta =(3/2)(B_0^{\prime }-1)$, and $B_0$, $B_0^{\prime }$ and $E_b$ are, respectively, the bulk modulus, its pressure derivative and cohesive energy at $V_{c0}$. Numerical inversion of Equation (A1) to get $V_c\left(P\right)$, for a given pressure P, is effected using an iteration method (starting with $X=1$) after rewriting it as $X=(X^{-1/3}+[P/\left(3B_0\right)]exp[-\eta (1-X^{1/3})]{)}^{-3/2}$. Note the finite limit of $E_v$ as $V\to 0$. But energy and pressure must approach those of QSM [20] in this limit. Pressure in QSM is

$$P_q\left(V\right)={\displaystyle \frac{{\hslash }^2}{2m_e}}[{\displaystyle \frac{2}{5}}{\left(3{\pi }^2\right)}^{2/3}{N}_s^{5/3}-{\displaystyle \frac{2}{a_B}}{\displaystyle \frac{13}{36}}{(3/\pi )}^{1/3}{N}_s^{4/3}],N_s\left(V\right)={\displaystyle \frac{Z}{V}}Exp[-\alpha ({\displaystyle \frac{V}{V_{c0}}}{)}^{1/3}-\beta \left({\displaystyle \frac{V}{V_{c0}}}{)}^{2/3}\right]$$

(A2)

Here, $\alpha =0.1935Z^{[0.495-0.039\left({log}_{10}Z\right)]}(R/a_B)$ while $\beta =[0.068+0.078\left({log}_{10}Z\right)-0.086{\left({log}_{10}Z\right)}^2]{(R/a_B)}^2$. Further, ℏ is reduced Planck’s constant, $m_e$ electron mass, $a_B$ Bohr radius, Z atomic number and R the Wigner-Seitz cell radius. Energy per atom, $E_q\left(V\right)$, is obtained by integrating the equation $P=-dE/dV$ from a suitable initial volume, say $V_0$. In an interpolation procedure [22] employed to smoothly go over from Vinet model to QSM, a volume $V_q$ is chosen such that the corrected pressure and energy are, respectively, $P_c\left(V\right)=P_v\left(V\right)$ and $E_0\left(V\right)=E_v\left(V\right)$ for $V\ge V_q$. These are expressed as $E_c\left(V\right)=[E_q\left(V\right)-E_q\left(V_q\right)]B_i\left(V\right)+E_v\left(V_q\right)$, and similarly $P_c\left(V\right)=P_q\left(V\right)B_i\left(V\right)+[E_q\left(V\right)-E_q\left(V_q\right)]B_i^{\prime }\left(V\right)$ for $V\le V_q$. The parameters in interpolation function $B_i\left(V\right)=[1+b_{1}V+b_2{V}^{4/3}+b_3{V}^{5/3}]$ are determined such that $P_c\left(V\right)$ and its first two derivatives are continuous at $V_q$. The constants for Li₆D are $b_1=-200.81$, $b_2=889.27$, $b_3=-1038.42$ with $V_q=0.03132$, which is (approximately) the zero-temperature volume at the turning point in the Hugoniot (at $\sim 13.5$ TPa pressure, see Figure 7.

The zero temperature enthalpy parameter (see Section 2.2) is given as ${\chi }_c=P{\Gamma }_{th}\left(V_c\right)/B_0\left(V_c\right)$, where ${\Gamma }_{th}$ is the average Grüneisen parameter (see Appendix A.3).

Appendix A.2. Ionic EOS

Johnson’s model [24] for ionic thermal energy is based on the temperature-variation of constant-volume specific heat $C_V$ from the $T^3$ law to $3N{k}_B$ above Debye’s temperature and, finally, to $(3/2)N{k}_B$, where N is the number of particles per gram and $k_B$ Boltzmann’s constant. Johnson adds an extra contribution $3N{k}_B(T/T_M)$ in the interval $T_M$ to $1.2T_M$ (where $T_M$ denotes the melting temperature) to account for the heat of fusion. The decrease of $C_V$ from its value $3N{k}_B$ at $T_M$ to $9N{k}_B/4$ at $5T_M$ and thereafter a linear variation in $ln\left(T\right)$ to the ideal gas value $(3/2)N{k}_B$ are built in to the model. The original model is modified so as to include Debye model for acoustic modes and Einstein model for optical modes of lattice vibrations. Then, $C_V$ reduces from $6N{k}_B$ to ideal gas value of $3N{k}_B$ (for Li and D atoms) after the melting transition. Further, it is also necessary to add rotation and vibration contributions of LiD molecule after melting, and its dissociation into atoms. Thus, specific thermal internal energy and thermal pressure are given by

$$\begin{array}{ccc}\hfill E_{it}(V,T)& =& [E_S(V,T)-E_S(V,0)]+N{k}_{B}T(E_0+E_{\psi }+E_1)\hfill \end{array}$$

(A3)

$$\begin{array}{ccc}\hfill P_{it}(V,T)& =& [P_S(V,T)-P_S(V,0)]+{\displaystyle \frac{1}{V}}(2{\Gamma }_{th}-2/3)N{k}_{B}T(E_0+E_{\psi }+E_1)\hfill \\ \hfill E_d(V,T)& =& N{k}_{B}T[{\displaystyle \frac{9}{8}}{\displaystyle \frac{{\theta }_d}{T}}+3I_d];I_d=3({\displaystyle \frac{T}{{\theta }_d}}{)}^3{\displaystyle \underset{0}{\overset{{\theta }_d/T}{\int }}}{\displaystyle \frac{z^3}{e^z-1}}dz;P_d={\displaystyle \frac{{\Gamma }_{ac}}{V}}{E}_d(V,T)\hfill \\ \hfill E_e(V,T)& =& 3({\mathcal{N}}_m-1)\left(N{k}_{B}T\right)[{\displaystyle \frac{{\theta }_e}{2T}}+{\displaystyle \frac{{\theta }_e}{T}}{\{exp\left({\displaystyle \frac{{\theta }_e}{T}}\right)-1\}}^{-1}];P_e={\displaystyle \frac{{\Gamma }_{op}}{V}}{E}_e(V,T)\hfill \end{array}$$

(A4)

Here $E_S=E_d+E_e$ and $P_S=P_d+P_e$ are the sums of Debye’s model and Einstein’s model specific internal energy and pressure, respectively. Further ${\mathcal{N}}_m$ is the number of atoms (two) in the LiD molecule so that $3({\mathcal{N}}_m-1)$ is the number of vibration modes, $\psi =T/T_M$ is scaled temperature, and ${\theta }_d$ and ${\theta }_e$ are the characteristic temperatures. The energy parameters $E_0$, $E_{\psi }$, and $E_1$, which are fitted functions of $\psi$ satisfying the constraints on $C_V$, are given by

$$\begin{array}{ccc}\hfill E_0& =& [3/2+3({\mathcal{N}}_m-1)](f_1+f_2);for1\le \psi <\infty ;\hfill \\ \hfill E_{\psi }& =& {\displaystyle \frac{3}{2}}(\psi -1/\psi ),for1\le \psi \le 1.2,E_{\psi }=0.66/\psi ,for1.2\le \psi <\infty \hfill \\ \hfill f_1& =& -1+a_4({\psi }^{-3/2}-{\psi }^{-2}/2);f_2={\displaystyle \frac{2}{3}}{a}_2[a_{3}y+{\psi }^{1-y}){\left({\psi }^y{(a_3+{\psi }^{1-y})}^2\right]}^{-1},\hfill \\ \hfill y& =& [201{\{1600a_1^2+2398(4a_1+5)\}}^{1/2}-40(5-197a_1)]{\left[3980(4a_1+5)\right]}^{-1}\hfill \\ \hfill a_1& =& -5.7,a_2={\displaystyle \frac{3}{2}}{(1+a_3)}^3{\left[a_3(1-y)(a_{3}y+2-y)\right]}^{-1},a_3=200,a_4=-{\displaystyle \frac{8}{5}}(a_1+a_2/(1+a_3))\hfill \end{array}$$

(A5)

The contribution to specific internal energy due to heat of fusion is ${ϵ}_{\psi }$, and $E_0$ describes the variation of $C_V$ after melting. The factor $f_1+f_2$, which varies in (0, -1) for $1\le \psi <\infty$, provides the interpolation from $6N{k}_B$ (as the acoustic and optical modes contribute $3N{k}_B$ each) to $3N{k}_B/2$, mentioned above. Thus, these terms treat LiD as a single unit through the melting transition. The additional term $\left(N{k}_{B}T\right)E_1$, which accounts for rotation, vibration and dissociation of free LiD molecules, is modeled as

$$E_1=[1-{\alpha }_d\left(T\right)]tanh\left[{\displaystyle \frac{\psi -1}{\omega }}\right][1+{\displaystyle \frac{{\theta }_{vib}}{T}}{\{exp\left({\displaystyle \frac{{\theta }_{vib}}{T}}\right)-1\}}^{-1}]+{\alpha }_d\left(T\right)[{\displaystyle \frac{3}{2}}+{\displaystyle \frac{{\theta }_{dis}}{T}}],1\le \psi <\infty$$

(A6)

where ${\alpha }_d\left(T\right)$ is the dissociation fraction. Here it is assumed that $N{k}_{B}T[3/2+{\theta }_{dis}/T]$ is the extra energy acquired by the atoms after dissociation, over and above that of the kinetic energy $(3/2)N{k}_{B}T$ present in the molecular term $E_0$. Before dissociation, the molecules have energy $N{k}_{B}T[1+{\displaystyle \frac{{\theta }_{vib}}{T}}{\{exp\left({\displaystyle \frac{{\theta }_{vib}}{T}}\right)-1\}}^{-1}]$ due to rotation and vibration. Classical energy of rotation $N{k}_{B}T$ is assumed; further ${\theta }_{dis}$ and ${\theta }_{vib}$ denote the dissociation and vibration temperature, respectively. The factor $tanh\left[{\displaystyle \frac{\psi -1}{\omega }}\right]$, with $\omega \sim 2$, provides smooth addition of these contributions. It is noted that the factor $(2{\Gamma }_{th}-2/3)$ in $P_{it}(V,T)$ facilitates smooth approach of ${\Gamma }_{th}$ to its ideal gas limit (2/3).

The constant-volume specific heat of the ionic model (omitting the contribution from $E_{\psi }$) is given by the sum $C_{Vi}=C_d+C_e+C_0+C_1$, where the terms are

$$\begin{array}{ccc}\hfill C_d& =& 3N{k}_B[4I_d-3{\displaystyle \frac{{\theta }_d}{T}}{\{exp\left({\displaystyle \frac{{\theta }_d}{T}}\right)-1\}}^{-1}];\hfill \\ \hfill C_e& =& 3({\mathcal{N}}_m-1)\left(N{k}_B\right)[{\left({\displaystyle \frac{{\theta }_e}{T}}\right)}^{2}exp\left({\displaystyle \frac{{\theta }_e}{T}}\right){\{exp\left({\displaystyle \frac{{\theta }_e}{T}}\right)-1\}}^{-1}];\hfill \\ \hfill C_0& =& N{k}_B[3/2+3({\mathcal{N}}_m-1)](f_3+f_4);f_3=-1-a_4({\psi }^{-3/2}-{\psi }^{-2})/2;\hfill \\ \hfill f_4& =& {(1+a_3)}^3[a_{3}y+(2-y){\psi }^{1-y}]{\left[(a_{3}y+2-y){\psi }^y{(a_3+{\psi }^{1-y})}^3\right]}^{-1},\hfill \\ \hfill C_1& =& N{k}_B[1-\alpha \left(T\right)]tanh\left[{\displaystyle \frac{\psi -1}{\omega }}\right][1+{\left({\displaystyle \frac{{\theta }_{vib}}{T}}\right)}^{2}exp\left({\displaystyle \frac{{\theta }_{vib}}{T}}\right){\{exp\left({\displaystyle \frac{{\theta }_{vib}}{T}}\right)-1\}}^{-2}]+N{k}_B\alpha \left(T\right)\left({\displaystyle \frac{3}{2}}\right).\hfill \end{array}$$

(A7)

The basic approach [44] to obtain ${\alpha }_d$ is to minimize the total free energy of the mixture (LiD, Li and D), including the entropy of mixing. This yields ${\alpha }_d=\sqrt{\mathcal{K}/(1+\mathcal{K})}$ where $\mathcal{K}=exp[(F_{LiD}-F_{Li}-F_D)/k_{B}T]$ with $F_{LiD}$, etc. denoting the free energy of the components of the mixture. Substitution of the expressions of these free energies [1] readily yields

$$\mathcal{K}={\displaystyle \frac{V}{N_{LiD}}}[{\displaystyle \frac{2\pi m_{Li}{m}_D{k}_{B}T}{m_{LiD}{h}^2}}{]}^{3/2}{\displaystyle \frac{1}{Z_r{Z}_v}}{\displaystyle \frac{g_{Li}{g}_D}{g_{LiD}}}exp[-U/k_{B}T]$$

(A8)

Here, U is the dissociation energy, $N_{LiD}$ the number of un-dissociated molecules per gram, $g_D=g_{Li}=2$ and $g_{LiD}=1$ the multiplicities of electronic states, $m_{LiD}$, etc. the particle masses, and h Planck’s constant. It is adequate to use the classical rotational partition function $Z_r=\left(k_{B}T\right)(8{\pi }^{2}I/h^2)$, where I is the moment of inertia of the molecule [1]. However, for the vibration partition function $Z_v$, a finite number of vibration states, $N_d\sim U/h\nu$ where $\nu$ denotes the vibration frequency, are considered so that $Z_r={\sum }_{k=1}^{N_d}exp[-k(h\nu /k_{B}T)]$. For the numerical work, $h\nu /k_B=2500$ K and $U/k_B=29000$ K are employed. Interaction between the particles is neglected as the dissociation temperature is high.

Appendix A.3. Ionic Grüneisen Parameter

Ionic Grüneisen parameters $({\Gamma }_{th},{\Gamma }_{ac},{\Gamma }_{op})$, Debye / Einstein temperatures $({\theta }_d,{\theta }_e)$ and melting temperature $\left(T_M\right)$ are also needed in the models. It is assumed that ${\Gamma }_{th}={\Gamma }_{ac}={\Gamma }_{op}$. The formula due to Preston $\mathit{et}\mathit{al}$[26] is ${\Gamma }_{th}\left(V\right)=(1/2)+{\gamma }_1{(V/V_0)}^{1/3}+{\gamma }_2{(V/V_0)}^q$, which has the correct asymptotic behavior for strong compression. The parameters ${\gamma }_1,{\gamma }_2$, and q are determined using experimental data on ${\Gamma }_{th}$ at $T=300K$, at zero-pressure melting point, and asymptotic ($V\to 0$) approximation for free electron states [26]. Expression for ${\theta }_d\left(V\right)$, which follow from Debye-Grüneisen law, is ${\theta }_d\left(V\right)={\theta }_{d0}({\displaystyle \frac{V_0}{V}}{)}^{1/2}exp[3{\gamma }_1(V_0^{1/3}-V^{1/3})+{\displaystyle \frac{{\gamma }_2}{q}}(V_0^q-V^q)]$. Similarly, Lindemann’s law gives $T_M\left(V\right)=T_{M0}({\displaystyle \frac{V_r}{V}}{)}^{1/3}exp[6{\gamma }_1(V_r^{1/3}-V^{1/3})+{\displaystyle \frac{2{\gamma }_2}{q}}(V_r^q-V^q)]$. Values of reference parameters ${\theta }_{d0}$, $T_{M0}$ and $V_r$, which correspond to $P=0$, are given in Table A1 and Table A2 . The expression for ${\theta }_e\left(V\right)$ is the same as that of ${\theta }_d\left(V\right)$, but with ${\theta }_{e0}$ as the reference parameter.

Table A1. Parameters for EOS.

Material	${\mathit{\rho }}_0$	${\mathit{B}}_0$	${\mathit{B}}_0^{\prime }$	${\mathit{E}}_{\mathit{c}\mathit{o}\mathit{h}}$	${\mathit{C}}_{\mathit{V}0}$	${\mathit{\Gamma }}_0$	C	S
	g cm−3	GPa		k J g−1	J g −1 K−1		k m s−1
Cu	8.93	134.8	5.19	5.302	0.3851	2.14	4.14	1.408
LiD-N	0.886	32.1†	3.53†	10.14	3.6312	1.04	0.610	1.101
Li6D	0.8	33.1	3.73	10.14	3.0932	1.22	0.659	1.095

† 300 K values [38].

Table A1. Parameters for EOS.

Material	${\mathit{\rho }}_0$	${\mathit{B}}_0$	${\mathit{B}}_0^{\prime }$	${\mathit{E}}_{\mathit{c}\mathit{o}\mathit{h}}$	${\mathit{C}}_{\mathit{V}0}$	${\mathit{\Gamma }}_0$	C	S
	g cm−3	GPa		k J g−1	J g −1 K−1		k m s−1
Cu	8.93	134.8	5.19	5.302	0.3851	2.14	4.14	1.408
LiD-N	0.886	32.1†	3.53†	10.14	3.6312	1.04	0.610	1.101
Li6D	0.8	33.1	3.73	10.14	3.0932	1.22	0.659	1.095

† 300 K values [38].

Table A2. Parameters for EOS (continued).

Material	Y	${\mathit{V}}_{\mathit{r}}$	${\mathit{T}}_{\mathit{M}0}$	${\mathit{\theta }}_{\mathit{d}0}$	${\mathit{\theta }}_{\mathit{e}0}$	${\mathit{\gamma }}_1$	${\mathit{\gamma }}_2$	q
	GPa	cm3 g−1	K	K	K
Cu	0.448	0.1194	1358	343		0.901	0.7818	4.7
LiD-N	0.001	1.185	588	1030	812	0.1	0.5740	5.1
Li6D	0.001	1.312	529	995	812	0.1	0.4476	5.1

Table A2. Parameters for EOS (continued).

Material	Y	${\mathit{V}}_{\mathit{r}}$	${\mathit{T}}_{\mathit{M}0}$	${\mathit{\theta }}_{\mathit{d}0}$	${\mathit{\theta }}_{\mathit{e}0}$	${\mathit{\gamma }}_1$	${\mathit{\gamma }}_2$	q
	GPa	cm3 g−1	K	K	K
Cu	0.448	0.1194	1358	343		0.901	0.7818	4.7
LiD-N	0.001	1.185	588	1030	812	0.1	0.5740	5.1
Li6D	0.001	1.312	529	995	812	0.1	0.4476	5.1

Appendix A.4. Electron-EOS-Fits

The thermal electron properties are accurately provided by the finite-temperature TF model. Numerical solution of the nonlinear TF equation [28] are used for this purpose. Using the results of extensive calculations, covering the range ${10}^{-4}\le T{Z}^{-4/3}\le {10}^3$ and $0.9\le R{Z}^{1/3}\le 20$, rational function approximations for $P_e(R,T)$ and $E_{et}(R,T)$ are obtained as

$$E_{et}(R,T)=Z^{7/3}{\displaystyle \frac{e_1{t}^2+e_2{t}^3}{1+e_3{t}^{e_4}+t^2}};P_{et}(R,T)=Z^{10/3}{\displaystyle \frac{p_1{t}^2+p_2{t}^3}{1+p_3{t}^{p_4}+t^2}};$$

(A9)

Here, $e_n$ and $p_n$, for $1\le n\le 4$, are functions of the scaled cell radius $R{Z}^{1/3}$ and $t=T/Z^{4/3}$ is the scaled temperature. The rational functions reduce to asymptotic form $\sim t^2$ and $\sim t$ in the limits $t\to 0$ and $t\to \infty$, respectively. Therefore, it is evident that $e_2=3/2$ and $p_2={\left(ZV\right)}^{-1}$. The parameters $e_1$ and $p_1$ are already determined by Gilvarry using solutions of temperature-perturbed TF equation [47,48]. These are given by

$$\begin{array}{ccc}\hfill e_1& =& P_0(15/2)\left(ZV\right)(\sigma +2\tau +3\omega )\xi {\left[Z^{4/3}\right]}^2;p_1=P_0(5/2)(\sigma +2\tau )\xi {\left[Z^{4/3}\right]}^2\hfill \\ \hfill \sigma & =& {\chi }_b/{\varphi }_b;\omega ={\chi }_i^{\prime }/\left(x_b^{1/2}{\varphi }_b^{5/2}\right);\tau ={(x_b/{\varphi }_b)}^2\hfill \\ \hfill {\varphi }_b& =& [0.48075x_b^{2/2}+0.06934x_b^{4/2}+0.0097x_b^{5/2}+0.0033704x_b^{6/2}{]}^{-1}\hfill \\ \hfill {\chi }_b& =& -0.3205x_b^3-0.04021x_b^4-0.002519x_b^5,{\chi }_i^{\prime }=[-0.005805x_b^{-0.228}-0.3755x_b^{-1}-3.12x_b^{-2}{]}^{-1}\hfill \end{array}$$

(A10)

Here $P_0$ is pressure at zero-temperature, $x_b=(R{Z}^{1/3}/{\mu }^{*})$ and $\xi =({\pi }^2/8){({\mu }^{*}{Z}^{-4/3}/e^2)}^2$. Further, ${\mu }^{*}=Z^{1/3}\mu =a_0{(9{\pi }^2/128)}^{1/3}$ is the TF screening length, $a_0$ is the Bohr radius. Remaining parameters are

$$\begin{array}{ccc}\hfill e_3& =& {\displaystyle \frac{-9.46531+2.1505x_b-0.33598x_b^2+0.00357102x_b^3-0.000045534x_b^4}{1-0.829696x_b+0.00159347x_b^2}}\hfill \\ \hfill e_4& =& {\displaystyle \frac{0.952738+0.01242x_b+0.00373893x_b^2}{1+0.0171359x_b+0.00734102x_b^2}}\hfill \\ \hfill p_3& =& {\displaystyle \frac{-12.795+3.98404x_b-0.659487x_b^2+0.000101473x_b^3-0.0000148037x_b^4}{1-0.828017x_b-0.00295839x_b^2}}\hfill \\ \hfill p_4& =& {\displaystyle \frac{0.735021-0.0368263x_b+0.0144916x_b^2-0.0000886233x_b^3}{1-0.172122x_b+0.0251956x_b^2}}\hfill \end{array}$$

(A11)

Gilvarry’s fits to zero-temperature pressure is $P_0=[Z^{10/3}{e}^2/10\pi {\mu }^{*4}]{\left\{{\sum }_{j=2}^6{c}_j{x}_b^{j/2+1}\right\}}^{-5/2}$ where $c_2=0.48075,c_3=0,c_4=0.06934,c_5=0.0097,$ and $c_6=0.0033704$ are constants.