<figure><img src="image_1.png" /></figure>

HAL
open science

# A three-phase equation of state for shock-induced phase transitions in tin

Corentin Stéphan, Hervé Jourdren, Nicolas Seguin

► To cite this version:

Corentin Stéphan, Hervé Jourdren, Nicolas Seguin. A three-phase equation of state for shock-induced phase transitions in tin. 2025. hal-05254026v1

HAL Id: hal-05254026

https://hal.science/hal-05254026v1

Preprint submitted on 18 Sep 2025 (v1), last revised 8 Jan 2026 (v2)

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
A THREE-PHASE EQUATION OF STATE FOR SHOCK-INDUCED
PHASE TRANSITIONS IN TIN

CORENTIN STÉPHAN<fnref n="1" />,<fnref n="2" />, HERVÉ JOURDREN<fnref n="2" />,<fnref n="3" /> AND NICOLAS SEGUIN<fnref n="1" />,<fnref n="4" />

**Abstract.** This paper addresses the construction of an equation of state considering three phases of tin for shock-induced phase transitions. First we present general results on such an equation of state which ensure the existence and uniqueness of the problem of maximization on the mixture entropy. Next, a multiphase equation of state under *strict* thermodynamic equilibrium is obtained using a combination of tabulation and Newton-Raphson iterations. Using this equation of state calibrated with standard static and dynamic experimental data, *negative* values of the fundamental derivative appear in the $(\beta/\gamma)$ mixture zone. A major consequence is the occurrence of composite waves for shock-induced $(\beta \to \gamma)$ phase transition. We construct the five self-similar reference solutions, using Hugoniot relations and Riemann invariants. Finally, numerical simulations are provided, illustrating these self-similar regimes of the $(\beta \to \gamma)$ phase transition at thermodynamic equilibrium for increasing piston velocities.

**Résumé.** Cet article aborde la construction d'une équation d'état prenant en compte trois phases de l'étain pour les transitions de phase induites par choc. Tout d'abord, nous présentons des résultats généraux sur une telle équation d'état, qui garantissent l'existence et l'unicité du problème de maximisation de l'entropie de mélange. Ensuite, une équation d'état multiphase à l'équilibre thermodynamique *strict* est obtenue en utilisant une méthode combinant tabulation et itérations de Newton-Raphson. En utilisant cette équation d'état calibrée sur des données expérimentales statiques et dynamiques de la littérature, des valeurs *négatives* de la dérivée fondamentale apparaissent dans la zone de mélange $(\beta/\gamma)$. Une conséquence importante est l'apparition d'ondes composites pour des transitions de phase $(\beta \to \gamma)$ sous choc. Nous construisons les cinq solutions auto-similaires de référence, en utilisant les relations de Rankine-Hugoniot et les invariants de Riemann. Enfin, des simulations numériques sont présentées, illustrant ces régimes auto-similaires de la transition de phase $(\beta \to \gamma)$ à l'équilibre thermodynamique pour des vitesses de piston croissantes.

<footnote n="1">¹ IMAG, Université de Montpellier, CNRS, Montpellier, France; e-mail: coretin.stephan@umontpellier.fr</footnote>

<footnote n="2">² CEA, DAM, DIF, F-91297 Arpajon, France</footnote>

<footnote n="3">³ Université Paris-Saclay, CEA DAM DIF, Laboratoire en Informatique Haute Performance pour le Calcul et la simulation, F-91297 Arpajon, France</footnote>

<footnote n="4">⁴ IMAG, Inria, Univ Montpellier, CNRS, Montpellier, France; e-mail: nicolas.seguin@inria.fr</footnote>
INTRODUCTION

Phase transitions appear in many applications and models can vary considerably depending on the material and experimental conditions. In the case of high dynamic loads, solids have to be considered as compressible materials, so that different solid phases of the material may appear. This question has been investigated for a long time, for instance in [8, 13, 36]. In this paper, we concentrate on solid $\beta$, solid $\gamma$, and liquid phases of tin. In particular, the phase diagram includes one triple point, which will be our reference to construct the equation of state describing the thermodynamic behavior of the material, including all the phase transitions.

Here, the dynamics of the material is governed by the Euler equations, written in time-mass Lagrangian coordinates $(t, m)$:

$$ \left\{ \begin{array}{l} \partial_t \tau - \partial_m u = 0, \\ \partial_t u + \partial_m P = 0, \\ \partial_t e + \partial_m (uP) = 0, \end{array} \right. \qquad (1) $$

where $\tau$, $u$ and $e$ respectively are the specific volume, the velocity and the specific total energy, while the pressure $P$ is linked to the other variables by

$$ \left\{ \begin{array}{l} \varepsilon = e - \frac{1}{2}|u|^2, \\ P = P(\tau, \varepsilon), \end{array} \right. \qquad (2) $$

where $\varepsilon$ is the specific internal energy. Such a model, based on mass, momentum and energy conservation is able to reproduce shock waves, but also dynamic phase transition fronts, when the pressure law is sufficiently realistic.

Equations of state which include phase transition do not fulfill classical convexity assumptions, leading to composite waves, so that the analysis and the approximation go beyond the standard framework, see for instance [18, 26].

When using classical Euler equations (1), all the thermodynamic information is encoded by function $P(\tau, \varepsilon)$ which has to provide appropriate values for separated phases but also good approximations of the transitions between phases. The construction of such pressure law is far from been obvious, direct interpolation of experimental measurements leads in general to numerical instabilities and possible crash of the simulation code. A common way of deriving a multiphase equation of state is to construct first a mixture entropy which is the weighted sum of the entropies of each phase, the weights being the mass fractions. To deduce the multiphase equation of state, it suffices to maximize this mixture entropy over all possible mixtures, which is equivalent to impose the thermodynamic equilibrium, see for instance [9, 14, 18] and reference therein. In this paper, we first sum up this process of construction in the case of an arbitrary number of phases, and propose a simple assumption which ensures the compatibility of the single phase entropies so that the thermodynamic equilibrium can be explicitly described. This assumption also enables a full description of triple points. Here, each phase is described by a complete Mie–Grüneisen equation of state [15, 22]. In practice, all this derivation is very sensitive with respect to the coefficients of the single phase equations of state, so that the whole process of construction of the three-phase equation of state has to be carefully done with respect to experimental data, see [31] for more details on the calibration and the data sets we used.

During numerical simulations with the Euler equations (1), the multiphase equation of state may be called several times by cell and time step, so that the associated procedure has to be very efficient. On the other hand, in order to obtain results of arbitrary precision, it seems unavoidable to solve
the maximization process when phase transitions are present. In order to balance the accuracy and the computational cost of the method, we propose a hybrid technique, using a tabulated partition of the phase space which ensures the direct use of single phase Mie-Grüneisen equations of state when only one phase is present, while the maximization problem is solved only in the vicinity of phase transition regions by a Newton-Raphson method initialized with a coarse grid tabulated data. (Note that this problem can be explicitly solved in the zone which corresponds to the triple point.) As a result, the overall cost of calculation is more than acceptable, since only cells near phase transition fronts require the use of the Newton-Raphson method.

The outline of this paper is the following. In the first section, we first recall the theory of general mixture of immiscible phases, showing in particular how the phase space is partitioned according to the maximum reached by the mixture entropy. We provide detailed description of the equilibrium mixture zones, including the triple point, by the use of hypothesis 1.2 between the single phase equations of state. In a second stage, we describe the hybrid method for the implementation of the equation of state: how the lines which separate the different zones are constructed and discretized, and how tabulation is used to call the maximization problem only on the mixture zones. In Section 3, we provide the three Mie-Grüneisen equations of state, with explicit values of all the coefficients for the $\beta$, $\gamma$ and liquid phases of tin. This allows the construction of exact self-similar solutions composed by constant states separated by shock waves, phase transition fronts, and isentropic compression waves. In the last section, some numerical results are given, using first- and second-order accurate Godunov-type schemes, showing expected profiles, with respect to exact solutions.

# 1. CONSTRUCTION OF A MULTIPHASE EQUATION OF STATE AT THERMODYNAMIC EQUILIBRIUM

We consider a continuous material and want to describe the possible phase transitions according to its specific volume $\tau$ and specific internal energy $\varepsilon$. Each possible phase is described by an associated equation of state, and we invoke the thermodynamic equilibrium to construct the complete equation of state, which includes phase transitions. We provide first a general framework for an arbitrary number of phases, under a compatibility assumption on the single phase equations of state. In the second part, we study the particular case of tin, with just three possible phases: solid $\beta$ and $\gamma$, and liquid.

## 1.1. The continuous material as a mixture of phases

The continuous compressible material is considered as a mixture of $n$ immiscible phases, $n \ge 2$, see for instance [6]. Each phase is represented by a number $\alpha \in [1, n]$. We define $\tau_\alpha$ and $\varepsilon_\alpha$ the specific volume and the specific internal energy of $\alpha$ phase, and the space of admissible states as

$$ \mathcal{F} = \{(\tau, \varepsilon) | \tau > 0, T > 0\}, \qquad (3) $$

where $T$ is the temperature, given later thanks to the equation of state (5).

For the moment, we assume that all the phases are defined for all values of $(\tau, \varepsilon)$ in $\mathcal{F}$, which is of course not the case from the physical point of view. The question is to identify the phases which are present from the knowledge of the specific volume of the mixture $\tau$, and its specific internal energy $\varepsilon$, leading to the so-called phase diagram. This is done by the computation of the thermodynamical equilibrium.
We assume that the thermodynamic behavior of each $\alpha$ phase of the material is described by a complete equation of state so that we can define its specific entropy

$$ (\tau, \varepsilon) \mapsto s_{\alpha}(\tau, \varepsilon) \qquad (4) $$

satisfying the fundamental relation of thermodynamics

$$ T_{\alpha}ds_{\alpha} = d\varepsilon_{\alpha} + P_{\alpha}d\tau_{\alpha}, \qquad (5) $$

meaning that pressure $P_\alpha$ and temperature $T_\alpha$ write

$$ P_{\alpha} = \frac{\partial_{\tau} s_{\alpha}}{\partial_{\varepsilon} s_{\alpha}}, \quad T_{\alpha} = \frac{1}{\partial_{\varepsilon} s_{\alpha}}, \qquad (6) $$

while the Gibbs free energy $G_\alpha$ is defined by $G_\alpha = \varepsilon_\alpha + P_\alpha\tau_\alpha - T_\alpha s_\alpha$. We also assume that the specific entropy $s_\alpha \in C^2(\mathcal{F})$ is strictly concave for each $\alpha \in [1, n]$, in the sense that the Hessian matrix $D^2s_\alpha$ is negative definite on $\mathcal{F}$.

We now consider a mixture of different phases of the material. We note $\tau$ the global specific volume, and $\varepsilon$ the global specific internal energy of the material, assuming $(\tau, \varepsilon) \in \mathcal{F}$. Since we consider immiscible phases, we have

$$ \begin{pmatrix} \tau \\ \varepsilon \end{pmatrix} = \sum_{\alpha=1}^{n} \lambda_{\alpha} \begin{pmatrix} \tau_{\alpha} \\ \varepsilon_{\alpha} \end{pmatrix} \quad \text{with} \quad \sum_{\alpha=1}^{n} \lambda_{\alpha} = 1, \qquad (7) $$

where $\lambda_\alpha \in [0, 1]$ is the mass fraction of $\alpha$ phase. In the case when all phase states $(\tau_\alpha, \varepsilon_\alpha)_{1\le\alpha\le n}$ are different, the mass fractions can be seen as barycentric coordinates:

$$ \lambda_{\alpha} = \frac{\tau - \sum_{\theta \neq \alpha} \lambda_{\theta} \tau_{\theta}}{\tau_{\alpha}} = \frac{\varepsilon - \sum_{\theta \neq \alpha} \lambda_{\theta} \varepsilon_{\theta}}{\varepsilon_{\alpha}}. \qquad (8) $$

One can also define two other fractions of a $\alpha$ phase, namely the volume fraction $z_\alpha$ and the energy fraction $\psi_\alpha$, by

$$ \begin{pmatrix} z_{\alpha} \\ \psi_{\alpha} \end{pmatrix} = \lambda_{\alpha} \begin{pmatrix} \tau_{\alpha} / \tau \\ \varepsilon_{\alpha} / \varepsilon \end{pmatrix}. \qquad (9) $$

Here again, these fractions belong to $[0, 1]$ and their sum is equal to 1. It is worth to note that the following equivalences hold: for any $(\tau_\alpha, \varepsilon_\alpha)$ and $(\tau, \varepsilon)$ in $\mathcal{F}$,

$$ \lambda_{\alpha} = 0 \iff z_{\alpha} = 0 \iff \psi_{\alpha} = 0. \qquad (10) $$

We now consider a given global state $(\tau, \varepsilon) \in \mathcal{F}$ of the mixture, and define the set of admissible mixtures by

$$ \mathcal{D}(\tau, \varepsilon) = \{ (\lambda_{\alpha}, \tau_{\alpha}, \varepsilon_{\alpha})_{1\le\alpha\le n} \in ([0, 1] \times \mathcal{F})^n \mid \text{equalities (7) hold} \}. \qquad (11) $$

From now on, we assume that $(\tau, \varepsilon) \in \mathcal{F}$ and $(\lambda_\alpha, \tau_\alpha, \varepsilon_\alpha)_{1\le\alpha\le n} \in \mathcal{D}(\tau, \varepsilon)$, and introduce the notation of the triplet $W_\alpha = (\lambda_\alpha, \tau_\alpha, \varepsilon_\alpha) \in [0, 1] \times \mathcal{F}$. By additivity, the specific entropy of mixture is defined as

$$ s : (\tau, \varepsilon, (W_{\alpha})_{1\le\alpha\le n}) \mapsto \sum_{\alpha=1}^{n} \lambda_{\alpha} s_{\alpha}(\tau_{\alpha}, \varepsilon_{\alpha}). \qquad (12) $$
Following for instance [9], one can show that the specific mixture entropy at a fixed value of $(\tau, \varepsilon)$: $(W_\alpha)_{1\le\alpha\le n} \mapsto s(\tau, \varepsilon, (W_\alpha)_{1\le\alpha\le n})$ is concave but not strictly concave in $D(\tau, \varepsilon)$.

## 1.2. The mixture at thermodynamic equilibrium

We can now study the thermodynamic equilibrium in order to deduce the phases which remain present and their repartition in terms of mass fraction. To do so, it is sufficient to maximize the mixture entropy $s$ on the set of all physical mixtures described by equations (7) and the set $D(\tau, \varepsilon)$.

At this stage, it is important to note that the thermodynamic equilibrium and the associated phase diagram will be described in terms of variables $\tau$ and $\varepsilon$, while it may be simpler to use the variables $P$ and $T$, see [10]. The reason is that we intend to include the resulting equilibrium equation of state in Euler equations, where $\tau$ and $\varepsilon$ are primitive variables, see (1)-(2).

**Definition 1.1** (Thermodynamic equilibrium). For all $(\tau, \varepsilon) \in \mathcal{F}$, the equilibrium entropy $s^{eq}$ is defined by

$$s^{eq} : (\tau, \varepsilon) \mapsto \max_{(W_\alpha)_{1\le\alpha\le n} \in D(\tau, \varepsilon)} s(\tau, \varepsilon, (W_\alpha)_{1\le\alpha\le n}). \quad (13)$$

The equilibrium entropy $s^{eq}$ describes the thermodynamic behavior of the continuous material. As a consequence, the thermodynamic equilibrium associated with a global state $(\tau, \varepsilon)$ is defined by any admissible mixture $(W_\alpha^*)_\alpha$ that maximizes the mixture entropy $s$.

Here are some properties of the equilibrium entropy $s^{eq}$, see for instance the Faccanioni's thesis [9]:

* $s^{eq} \in C^1(\mathcal{F})$,

* $s^{eq}$ is concave but not strictly,

* the graph of $s^{eq}$ is the concave hull of $(\tau, \varepsilon) \mapsto \max_{\alpha \in [1, n]} s_\alpha(\tau, \varepsilon)$.

Let us now describe the mixture at equilibrium, and in particular identify the phases which are present. To do so, additional assumptions are needed since many degenerate cases can occur. In [9], the case of two phases is studied and compatibility conditions are introduced between the phase entropies. In the case of a general number of phases $n \ge 3$, it seems difficult to extend these compatibility conditions, and we prefer to introduce a condition using the tangent plane to the graph of $s^{eq}$, which implicitly constrains phase entropies $s_\alpha$ between each other.

For any given state $w = (\tau, \varepsilon) \in \mathcal{F}$, we consider the tangent plane to the graph of $s^{eq}$ at point $w$, defined by the affine function

$$\forall \tilde{w} \in \mathcal{F}, \quad \phi_w^{eq}(\tilde{w}) = s^{eq}(w) + \nabla s^{eq}(w)(\tilde{w} - w).$$

Let $A(w) \subset [1, n]$ be the set of phases whose graph of entropy $s_\alpha$ intersects this tangent plane, that is to say

$$A(w) := \{\alpha \in [1, n] \text{ s.t. } \exists w_\alpha^* \in \mathcal{F}, s_\alpha(w_\alpha^*) = \phi_w^{eq}(w_\alpha^*)\}. \quad (14)$$

For each $\alpha \in A(w)$ phase, the (strict) concavity of $s_\alpha$ implies the uniqueness of the intersection point $w_\alpha^*$. For the rest, we explicit the dependence w.r.t. $w$, i.e. we use the notation $w_\alpha^*(w)$ for these intersection points. In order to eliminate the pathological cases and focus on realistic situations, we make the following assumption:

**Hypothesis 1.2.** For all $w \in \mathcal{F}$, $\#A(w) = \dim C(w) + 1$, where $C(w)$ the set of states $\tilde{w} \in \mathcal{F}$ such that $s^{eq}(\tilde{w}) = \phi_w^{eq}(\tilde{w})$.
Since $\mathcal{C}(w) \subset \mathcal{F}$, $\dim \mathcal{C}(w) \le 2$. Moreover, $\mathcal{C}(w)$ is a closed convex set, which includes the states $(w_{\alpha}^{*}(w))_{\alpha \in A(w)}$. Thanks to the assumption, $\mathcal{C}(w)$ is uniquely defined by the states $(w_{\alpha}^{*}(w))_{\alpha \in A(w)}$ which are its vertices, and is a simplex: $\mathcal{C}(w)$ is a point if $\#A(w) = 1$, a segment if $\#A(w) = 2$ and a triangle if $\#A(w) = 3$.

We are now in position to describe the mixture of phases at thermodynamic equilibrium.

**Theorem 1.3.** Consider a state $w = (\tau, \varepsilon) \in \mathcal{F}$, the associated set $A(w)$ and states $(w_{\alpha}^{*}(w))_{\alpha \in A(w)}$ defined by (14).

(1) Let $(\lambda_{\alpha}^{*}(w))_{\alpha \in A}$ be the (unique) barycenter coordinates of $w$ in $\mathcal{C}(w)$. The partial mixture described by $(\lambda_{\alpha}^{*}(w), w_{\alpha}^{*}(w))_{\alpha \in A}$ provides a maximizer of problem (13), i.e.

$$s^{eq}(w) = \sum_{\alpha \in A} \lambda_{\alpha}^{*}(w) s_{\alpha}(w_{\alpha}^{*}(w)), \qquad (15)$$

so that any maximizer of (13) can be written as

$$\begin{cases} (\lambda_{\alpha}^{*}(w), w_{\alpha}^{*}(w)) & \text{if } \alpha \in A(w), \\ (0, \tilde{w}_{\alpha}) & \text{else,} \end{cases}$$

where $\tilde{w}_{\alpha}$ can be any state of $\mathcal{F}$.

(2) The phases in $A(w)$ are in thermodynamic equilibrium, i.e. if we define

$$(P^{eq}, T^{eq}, G^{eq})(w) = \left( \frac{\partial_{\tau} s^{eq}(w)}{\partial_{\varepsilon} s^{eq}(w)}, \frac{1}{\partial_{\varepsilon} s^{eq}(w)}, \varepsilon + P^{eq}_{\tau} - T^{eq}(w) s^{eq}(w) \right)$$

then

$$(P, T, G)_{\alpha}(w_{\alpha}^{*}(w)) = (P^{eq}, T^{eq}, G^{eq})(w) \quad \forall \alpha \in A(w). \qquad (16)$$

*Proof.* (1) By definition of the set $\mathcal{C}(w)$, $w \in \mathcal{C}(w)$. As a result, since $\mathcal{C}(w)$ is a simplex, the barycenter coordinates of $w$ with respect to the vertices of $\mathcal{C}(w)$ exist and are unique, i.e. the decomposition

$$w = \sum_{\alpha \in A} \lambda_{\alpha}^{*}(w) w_{\alpha}^{*}(w),$$

holds and is unique. Then, the properties of $\phi_{w}^{eq}$ and the definition of the $w_{\alpha}^{*}$ directly lead to (15). The general form of the maximizers is a straightforward consequence of (15), and the fact that there is no maximizer with another form is deduced from the uniqueness of the $(w_{\alpha}^{*})_{\alpha \in A(w)}$ and $(\lambda_{\alpha}^{*})_{\alpha \in A(w)}$.

(2) The equality of pressures and temperatures are deduced from the equalities

$$\nabla s_{\alpha}(w_{\alpha}^{*}(w)) = \nabla s^{eq}(w) \quad \forall \alpha \in A(w).$$

On the other hand, since $s^{eq}$ is affine in $\mathcal{C}(w)$,

$$s^{eq}(w) - s^{eq}(w_{\alpha}^{*}(w)) = \nabla s^{eq}(w) (w - w_{\alpha}^{*}(w)).$$

Rearranging the terms and using the previous equalities lead to

$$G_{\alpha}(\tau_{\alpha}^{*}, \varepsilon_{\alpha}^{*}) = \varepsilon_{\alpha}^{*} + P_{\alpha}(\tau_{\alpha}^{*}, \varepsilon_{\alpha}^{*})\tau_{\alpha}^{*} - T_{\alpha}(\tau_{\alpha}^{*}, \varepsilon_{\alpha}^{*})s_{\alpha}(\tau_{\alpha}^{*}, \varepsilon_{\alpha}^{*}) = G^{eq}(w) \quad \forall \alpha \in A(w),$$
with the definition of $G^{eq}(w)$ provided in the theorem. □

**Remark 1.4.** In this presentation, $w \in \mathcal{F}$ is fixed, and then we deduce the phases which are present. One can also reverse the point of view. Let us consider $a \subset [1, n]$ and define $\mathcal{F}_a = \{w \in \mathcal{F}, A(w) = a\}$. This set can be described by the use of its bounds: they can be bounds of $\mathcal{F}$, but also common bounds with sets $\mathcal{F}_{a'}$, with $a \cap a' \neq \emptyset$. Doing so with all possible $a$, one can retrieve the phase diagram in the $(\tau, \varepsilon)$ plane. In the $(P, T)$ plane, the thermodynamical equilibrium and assumption 1.2 lead to two-dimensional sets for single phases, separated by curves. By assumption 1.2, curves correspond to the coexistence of two phases and intersection points to local triple points. Note that, still using assumption 1.2, a triple point in the $(P, T)$ plane corresponds to a triangle in the $(\tau, \varepsilon)$ plane.

**Remark 1.5.** In our description of the thermodynamical equilibrium, some of barycenter coordinates in (15) can vanish: $\lambda_{\alpha}^*(w)$ belongs to $[0, 1]$ and the bounds of this interval can be reached. As a result, any phase may belong to $A(w)$ but with a mass fraction equal to zero, so that $A(w)$ is not exactly the set of phases which are really present, i.e. with a positive mass fraction. This is due to the geometric description we choose and the fact that $C(w)$ is a closed set.

### 1.3. Construction of a three-phase diagram

We apply the previous results to the case of three possible phases ($n = 3$). In reference to the application investigated later, the set of the three phases is $\mathcal{P} = \{\beta, \gamma, \ell\}$.

**Hypothesis 1.6.** There exists a unique triple point, which is described in the $(\tau, \varepsilon)$ plane by a triplet of states $(w_{\alpha}^t)_{\alpha \in \mathcal{P}} \in \mathcal{F}^3$, which are not aligned in $\mathcal{F}$, such that

$$A(w) = \{\beta, \gamma, \ell\} \qquad (17)$$

for all $w$ in the closed triangle $\mathcal{T}$ whose vertices are $w_{\beta}^t$, $w_{\gamma}^t$ and $w_{\ell}^t$.

The phase diagram includes three lines of phase-change (noted here $\beta/\gamma$, $\beta/\ell$ and $\gamma/\ell$) with shapes as in figure 1.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. Illustration of a three-phase diagram in the $(P, T)$ plane and $(\tau, \varepsilon)$ plane.</figcaption></figure>
### 1.3.1. The triple point

The triple point is a point in the $(P, T)$ diagram such that the three phases coexist in thermodynamical equilibrium. Consequently, hypothesis 1.6 implies the existence of the triplet $(T^t, P^t, G^t)$ such as

$$ egin{cases} P_{\beta}((\tau, \varepsilon)_{\beta}^{t}) = P_{\gamma}((\tau, \varepsilon)_{\gamma}^{t}) = P_{\ell}((\tau, \varepsilon)_{\ell}^{t}) = P^{t} \\ T_{\beta}((\tau, \varepsilon)_{\beta}^{t}) = T_{\gamma}((\tau, \varepsilon)_{\gamma}^{t}) = T_{\ell}((\tau, \varepsilon)_{\ell}^{t}) = T^{t} \\ G_{\beta}((\tau, \varepsilon)_{\beta}^{t}) = G_{\gamma}((\tau, \varepsilon)_{\gamma}^{t}) = G_{\ell}((\tau, \varepsilon)_{\ell}^{t}) = G^{t}. \end{cases} \quad (18) $$

Indeed all the points $(\tau, \varepsilon)$ living in $\mathcal{T}$ have for unique maximizer the triplet $(\tau_{\alpha}^{t}, \varepsilon_{\alpha}^{t})_{\alpha \in \mathcal{P}}$ and their barycentric coordinates in this system.

### 1.3.2. Definition of the phase-change lines

The phase-change lines separate two regions of pure phases $a, b \in \mathcal{P}$, $a \neq b$, in the $(P, T)$ plane. They are bounded on one side by the triple point $(P^t, T^t)$. In order to properly define these lines, we make the following assumption (which is satisfied in the case of tin, see figure 1):

**Hypothesis 1.7.** The phase-change lines are monotone functions w.r.t. $P$ in the $(P, T)$ plane.

We are then able to consider the function

$$ F_{a/b} : P^* \mapsto \left( \begin{array}{c} T \\ (\tau_a, \varepsilon_a) \\ (\tau_b, \varepsilon_b) \end{array} \right) \quad (19) $$

defined by the thermodynamic equilibrium between phases $a$ and $b$, i.e.

$$ egin{cases} P_a(\tau_a, \varepsilon_a) = P_b(\tau_b, \varepsilon_b) = P^* \\ T_a(\tau_a, \varepsilon_a) = T_b(\tau_b, \varepsilon_b) = T \\ G_a(\tau_a, \varepsilon_a) = G_b(\tau_b, \varepsilon_b) = G \end{cases} \quad (20) $$

where $(\tau_a, \varepsilon_a)(P^*), (\tau_b, \varepsilon_b)(P^*) \in \mathcal{F}$.

We are now in position to define the different phase-change sets according to the phase plane we consider:

* In the $(P, T)$ plane, the phase-change line is the graph of the first component of $F_{a/b}: \mathcal{L}_{a/b}$ is the set of points $(P, (F_{a/b})^{(1)}(P))$ where $P < P^t$ if $a = \ell$ and $b = \beta$ (or the contrary), and $P > P^t$ otherwise, see figure 1-left.

* In the $(\tau, \varepsilon)$ plane, the zone of coexistence of the two phases $a$ and $b$ is bounded by the edge $[(\tau_a, \varepsilon_a), (\tau_b, \varepsilon_b)]$ of $\mathcal{T}$ and two curves, defined by

$$ \mathcal{L}_{a/b}^{a} = \{(\tau_a, \varepsilon_a) = (F_{a/b})^{(2)}(P), P\} \quad (21) $$

$$ \mathcal{L}_{a/b}^{b} = \{(\tau_b, \varepsilon_b) = (F_{a/b})^{(3)}(P), P\} \quad (22) $$

where $P < P^t$ if $a = \ell$ and $b = \beta$ (or the contrary), and $P > P^t$ otherwise, see figure 1-right.
## 2. NUMERICAL IMPLEMENTATION OF THE EQUATION OF STATE

Using the theoretical framework of the previous section and some additional monotonicity assumptions, we construct and describe the phase diagrams in the $(P, T)$ and $(\tau, \varepsilon)$ planes. Considering a tabulation of this latter plane, a hybrid approach is proposed, enabling the computation of the thermodynamic equilibrium only in the multiphase regions.

### 2.1. Computation of the phase diagram

#### 2.1.1. The triple point and associated triangle

Let us introduce the function

$$F^t: ((\tau_\beta, \varepsilon_\beta), (\tau_\gamma, \varepsilon_\gamma), (\tau_\ell, \varepsilon_\ell)) \mapsto \begin{pmatrix} P_\beta - P_\gamma \\ P_\beta - P_\ell \\ T_\beta - T_\gamma \\ T_\beta - T_\ell \\ G_\beta - G_\gamma \\ G_\beta - G_\ell \end{pmatrix}. \qquad (23)$$

Then the unique solution $((\tau, \varepsilon)_\beta^t, (\tau, \varepsilon)_\gamma^t, (\tau, \varepsilon)_\ell^t)$ of equalities (18) verifies

$$F^t((\tau, \varepsilon)_\beta^t, (\tau, \varepsilon)_\gamma^t, (\tau, \varepsilon)_\ell^t) = 0. \qquad (24)$$

The solution provides the triangle $\mathcal{T}$ and the triple point $(P^t, T^t)$. We can solve this system with a classical Newton-Raphson algorithm.

#### 2.1.2. Discretization of the phase-change lines in the $(P, T)$ plane

Once the triple point is obtained, we discretize the phase-change lines given by the equations (19) and (20). The two phases in presence are noted $a, b \in \mathcal{P}$. As we can see on the diagram in the $(P, T)$ plane (figure 1), the phase-change lines are monotone functions of $P$ then the curves can be parametrized by the pressure. The pressure varies in an interval whose triple point must be one of the end points, the other endpoint has to be chosen according to the minimal/maximal possible values in the applications.

In practice, we consider a monotone discretization of this pressure segment by $N$ points: $(P^k)_{1 \le k \le N}$, with $P_1 = P^t$. This sequence is decreasing if $a = \ell$ and $b = \beta$ (or the contrary), increasing otherwise. For each of these pressures, we introduce a function $F_{P^k}$ such that

$$F_{P^k}: ((\tau_a, \varepsilon_a), (\tau_b, \varepsilon_b)) \mapsto \begin{pmatrix} P_a(\tau_a, \varepsilon_a) - P^k \\ P_b(\tau_b, \varepsilon_b) - P^k \\ T_a(\tau_a, \varepsilon_a) - T_b(\tau_b, \varepsilon_b) \\ G_a(\tau_a, \varepsilon_a) - G_b(\tau_b, \varepsilon_b) \end{pmatrix}. \qquad (25)$$

Therefore, for each $P^k$, one obtains two points, $(\tau, \varepsilon)_a^k$ and $(\tau, \varepsilon)_b^k$, which belong to the phase-change lines and verify

$$F_{P^k}((\tau, \varepsilon)_a^k, (\tau, \varepsilon)_b^k) = 0 \quad 1 \le k \le N. \qquad (26)$$
System (26) can be solved with a Newton-Raphson algorithm. For $k=1$, the solution is known and is the associated vertex of the triangle $\mathcal{T}$, i.e.

$$((\tau, \varepsilon)_a^1, (\tau, \varepsilon)_b^1) = ((\tau_a, \varepsilon_a)^t, (\tau_b, \varepsilon_b)^t).$$

For $k \ge 2$, the initial data of the Newton-Raphson algorithm is set to $((\tau, \varepsilon)_a^{k-1}, (\tau, \varepsilon)_b^{k-1})$, since we know that these curves are smooth. At the end, we obtain sets of points belonging to the phase-change lines in the $(\tau, \varepsilon)$ plane

$$\{(\tau, \varepsilon)_a^k, 1 \le k \le N\} \subset \mathcal{L}_{a/b}^a, \qquad (27)$$

$$\{(\tau, \varepsilon)_b^k, 1 \le k \le N\} \subset \mathcal{L}_{a/b}^b, \qquad (28)$$

and sets of points belonging to phase-change lines in $(P, T)$ plane

$$\{(P^k, T_a((\tau, \varepsilon)_a^k)), 1 \le k \le N\} \subset \mathcal{L}_{a/b}. \qquad (29)$$

From the numerical point of view, these sets of points will be used to construct piecewise linear approximations of the different phase-change lines.

**Remark 2.1.** We describe the phase diagram and saturation curves in the $(\tau, \varepsilon)$ plane, since these variables are those which are provided when the equation of state is used in compressible hydrodynamics codes. It is then very different from classical literature on equations of state for liquid-vapor phase transition where the $(P, T)$ plane is preferred, see for instance [28] and [10]. In the latter reference, only the saturation curves in the $(P, T)$ plane are approximated and phase-change lines in the $(\tau, \varepsilon)$ are deduced directly using pure phase equations of state. As it will be presented in the next section, our hybrid method is constructed such that the dependence on the discretization of the pressure segment is minimized, if not independent. In particular, thorough consistency tests have been performed and are described in [31].

### 2.1.3. Related phase diagram in the $(\tau, \varepsilon)$ plane

In the $(P, T)$ plane, multiphase regions are only curves and points. As a consequence, if $\mathcal{D}_{P,T}^\alpha$ denotes the domain of $\alpha$ phase, $\alpha \in \mathcal{P}$, the whole domain writes

$$\mathcal{D}_{P,T} = \bigcup_{\alpha \in \mathcal{P}} \overline{\mathcal{D}_{P,T}^\alpha}. \qquad (30)$$

As mentioned above, the phase-change lines are approximated by piecewise linear curves so that the domain actually are polygonal sets. For simplicity, we will use the same notations for the exact domains and their polygonal approximations.

Let us now describe the phase diagram in the $(\tau, \varepsilon)$ plane, according to the number of coexisting phases:

* Three-phase zone: this is by definition the set which corresponds to the triple point, and is the triangle $\mathcal{T}$ whose vertices are $(\tau, \varepsilon)_\beta^t$, $(\tau, \varepsilon)_\gamma^t$, and $(\tau, \varepsilon)_\ell^t$.

* Two-phase zones: let $a$ and $b$ the two different phases, then the associated set of coexistence is fully defined by its bounds, which are the segment between $(\tau, \varepsilon)_a^t$ and $(\tau, \varepsilon)_b^t$, the curve $\mathcal{L}_{a/b}^a$ and the curve $\mathcal{L}_{a/b}^b$. It is noted $\mathcal{D}_{\tau, \varepsilon}^{a,b}$ and there are three two-phase zones.
* Single-phase zones: let $\alpha \in \mathcal{P}$. The associated set, denoted $\mathcal{D}_{\tau, \varepsilon}^{\alpha}$, is fully defined by its bounds which are the curves $\mathcal{L}_{a/b}^{a}$ and $\mathcal{L}_{a/b}^{b}$, with $a, b \in \mathcal{P}$ such that $\alpha, a$ and $b$ are distinct from each other.

To sum up, the $(\tau, \varepsilon)$ plane is partitioned into seven zones, whose indices are in the set

$$ \mathcal{B} = \{\beta, \gamma, \ell, \beta/\gamma, \beta/\ell, \gamma/\ell, t\}. $$

In practice we do not consider all possible positive values of specific volume and temperature, since pressures and temperatures have to stay in the domain of validity of the Mie-Grüneisen EOS to be considered later. In order to identify in a unique manner the index in $\mathcal{B}$ for each $(\tau, \varepsilon)$, we assume that $\mathcal{T}$ is closed and that the single-phase zones are open sets. We represent in figure 2 the phase diagram in the $(P, T)$ plane and in the $(\tau, \varepsilon)$ plane.

<figure><img src="image_3.png" /><figcaption>FIGURE 2. The seven zones of phase diagrams in planes $(P, T)$ and $(\tau, \varepsilon)$.</figcaption></figure>

## 2.2. A hybrid implementation of the equation of state

Let us present now the way the equation of state is implemented, that is to say how to compute the equilibrium variables for any state $(\tau, \varepsilon) \in \mathcal{F}$.

### 2.2.1. Partition of the phase diagram

A naïve method would be to solve the full thermodynamical equilibrium system. However, it would be very expensive and the fact that all the phases are not always present may lead to numerical difficulties. Therefore, in order to avoid unnecessary calculations, we do not use the same calculation method for all zones, but rather perform off-line computations to partition the phase diagram into the seven regions described above. First, we precompute the vertices of the triangle $\mathcal{T}$ by solving (24) and the phase-change lines by the procedure described in section 2.1.2. Let us recall that the phase-change lines are discretized by piecewise linear curves, the seven zones to investigate actually are polygons, so that it is easy to determine the zone to which a point $(\tau, \varepsilon)$ belongs. Moreover, in the triangle $\mathcal{T}$, the pressure is constant and known, while in the single-phase zones, we can directly use the associated single-phase equation of state.

It remains to detail the computation of the pressure for $(\tau, \varepsilon)$ in some two-phase region, say $\mathcal{D}_{\tau, \varepsilon}^{a,b}$. We search $(\lambda_a, \tau_a, \varepsilon_a)$ and $(\lambda_b, \tau_b, \varepsilon_b)$ in $\mathcal{D}(\tau, \varepsilon)$ which satisfy the equalities $(P, T, G)(\tau_a, \varepsilon_a) =$
$(P, T, G)(\tau_b, \varepsilon_b)$. Since $\lambda = \lambda_b = 1 - \lambda_a$, we introduce the function $H_{a,b}^{\tau, \varepsilon}$ defined by

$$H_{a,b}^{\tau, \varepsilon} : (\lambda, (\tau_a, \varepsilon_a), (\tau_b, \varepsilon_b)) \mapsto \begin{pmatrix} P_a(\tau_a, \varepsilon_a) - P_b(\tau_b, \varepsilon_b) \\ T_a(\tau_a, \varepsilon_a) - T_b(\tau_b, \varepsilon_b) \\ G_a(\tau_a, \varepsilon_a) - G_b(\tau_b, \varepsilon_b) \\ (1-\lambda)\tau_a + \lambda\tau_b - \tau \\ (1-\lambda)\varepsilon_a + \lambda\varepsilon_b - \varepsilon \end{pmatrix}. \quad (31)$$

As a consequence, the thermodynamic equilibrium satisfies the system

$$H_{a,b}^{\tau, \varepsilon}(\lambda, (\tau, \varepsilon)_a, (\tau, \varepsilon)_b) = 0 \quad (32)$$

which is solved in practice with a classical Newton-Raphson method.

### 2.2.2. Optimization of the equation of state computation

In order to go further in the optimization of the code, we use a tabulation of the phase diagram in variable $(\tau, \varepsilon)$. The goal is to accelerate the method which enables to identify the zone to which a state $(\tau, \varepsilon)$ belongs. To do so, the full domain $\mathcal{F}$ is discretized by a uniform Cartesian grid, with indices $i$ for the $\tau$ direction, and $j$ for the $\varepsilon$ direction, and denote $R_{ij}$ the cells of this grid. For each cell $R_{ij}$, we save the zones that intersect it, and their number, i.e.

$$Q_{ij} := \{\alpha \in \mathcal{B} \mid \mathcal{D}_{\tau, \varepsilon}^{\alpha} \cap R_{ij} \neq \emptyset\}, \quad n_{ij} := \text{cardinal } Q_{ij}. \quad (33)$$

This way, we only have to find the cell that contains the state $(\tau, \varepsilon)$, which can be done by simple integer divisions (since the grid is a uniform Cartesian grid), and after the thermodynamic variables can be obtained using the calculations described previously.

### 2.2.3. The full hybrid algorithm

Let us now describe the full algorithm to compute the equilibrium pressure and temperature from a given state $(\tau, \varepsilon)$, separating the preprocessing and the remaining part of the method.

#### PREPROCESSING

Compute the vertices $(\tau, \varepsilon)_\beta^t$, $(\tau, \varepsilon)_\gamma^t$, and $(\tau, \varepsilon)_\ell^t$ of the triangle $\mathcal{T}$, using (23)-(24)

Compute the phase-change lines $\mathcal{L}_{a/b}^a$ and $\mathcal{L}_{a/b}^b$, $a, b \in \mathcal{P}$, $a \neq b$, following section 2.1.2

Construct the grid $(R_{ij})_{i,j}$ and compute $Q_{ij}$ and $n_{ij}$ for all $(i, j)$, using (33)

The purpose of the following algorithm is to compute the equilibrium pressure and temperature from any state $(\tau, \varepsilon) \in \mathcal{F}$.

#### ALGORITHM FOR THE EQUATION OF STATE

Let $(\tau, \varepsilon) \in \mathcal{F}$ be given

Find $(i, j)$ such that $(\tau, \varepsilon) \in R_{ij}$

if $n_{ij} = 1$ then

▷ The cell contains one zone

Let $\alpha \in \mathcal{B}$ such that $Q_{ij} = \{\alpha\}$

if $\alpha = \{\beta, \gamma, \ell\}$ then

▷ $(\tau, \varepsilon)$ is in a single-phase zone

$$(P, T)(\tau, \varepsilon) = (P, T)_\alpha(\tau, \varepsilon)$$

else

if $\alpha \in \{\beta/\gamma, \beta/\ell, \gamma/\ell\}$ then

▷ $(\tau, \varepsilon)$ is in a two-phase zone
Solve the two-phase equilibrium (31)-(32)

else

$$ (P, T) = (P^t, T^t) $$

end if

end if

else

if $n_{ij} = 2$ then

if $\{\beta/\gamma, \beta/\ell, \gamma/\ell\} \cap Q_{ij} \neq \emptyset$ then

Solve the two-phase equilibrium (31)-(32)

else

Determine the zone which contains $(\tau, \varepsilon)$

Compute $(P, T)$ as in the case $n_{ij} = 1$

end if

else

Determine the zone which contains $(\tau, \varepsilon)$

Compute $(P, T)$ as in the case $n_{ij} = 1$

end if

end if

▷ $(\tau, \varepsilon)$ is in the three-phase zone $\mathcal{T}$

▷ The cell contains two zones

▷ (⋆)

▷ (‡)

▷ The cell contains more than two zones

▷ (‡)

In the case $(\star)$, the resolution of the two-phase equilibrium (31)-(32) may lead to wrong mass fractions if the state $(\tau, \varepsilon)$ belongs either to a single-phase case or to the triangle $\mathcal{T}$. In both cases, the solution of two-phase equilibrium (31)-(32) is replaced by the appropriate computation of $(P, T)$. One may also remark that the triple point is computed only once, during the preprocessing step, since the resolution of system (23)-(24) can be expensive.

In the cases $(\ddagger)$, the zone is obtained by the use of a classical detection method, since the zones are polygonal sets.

Most of the cells of the grid contain only one zone so that, on the majority of the domain, the step of detection of the zones is avoided. Moreover a fine tabulation (for instance a 900 × 900 grid) will increase the proportion of tabulation elements with only one zone.

**Remark 2.2.** Using 64-bit IEEE, the adimensional convergence criteria is set to $10^{-11}$ on the variation of the mass fractions in the Newton–Raphson algorithm considering *exact* analytical expression of the Jacobian matrix.

**Remark 2.3.** As mentioned above, the tabulation is used to detect the number of phases present in each cell $R_{ij}$, but is also used to store local estimations of the equilibrium values of the variables, which serve as initializations of the Newton–Raphson method.

**Remark 2.4.** We assume here that the tabulation of the phase diagram is fine enough, so that the computation of the cardinals $n_{ij}$ is accurate, despite the fact that the phase-change lines are approximated by piecewise linear curves. As already mentioned in remark 2.1, since the thermodynamic equilibrium is computed as soon as the point $(\tau, \varepsilon)$ is close to a two-phase zone, the approximation of the phase-change lines does not lead to numerical difficulties. More details are given in [31].

### 2.2.4. Thermodynamic consistency

The following formula resulting from the Maxwell's rule (Schwarz theorem)

$$ T^{eq} \frac{\partial P^{eq}}{\partial T^{eq}}|_{\tau} = P^{eq} + \frac{\partial \varepsilon^{eq}}{\partial \tau}|_{T} \qquad (34) $$
is verified up to machine precision (see [31] for more details).

### 2.2.5. Computation of the sound speed and fundamental derivative

In addition to pressure and temperature, other quantities are needed when implementing numerical schemes.

The speed of sound depends on the derivative of the pressure with respect to $\tau$ at constant entropy $s$.

$$c^2 = -\tau^2 \left( \frac{\partial P}{\partial \tau} \right)_s . \qquad (35)$$

In contrast to a pure phase zone, the speed of sound cannot be calculated with an explicit formula in two-phase mixtures. We use finite differences to approximate the derivative of pressure. In the triangle $\mathcal{T}$ corresponding to three-phase mixtures, the pressure is constant. Therefore, the derivative of pressure, as well as the speed of sound, is zero.

**Pure phase zone.** The sound speed and the fundamental derivative can be expressed thanks to partial derivatives with respect to $\tau$ and $\varepsilon$. We thus obtain analytical expressions for the pure phases:

$$c^2(\tau, \varepsilon) = -\tau^2 (\partial_{\tau} P - P \partial_{\varepsilon} P), \qquad (36)$$

$$\mathcal{G}(\tau, \varepsilon) = -\frac{\tau}{2} \left( \frac{\partial^2 P}{\partial \tau^2} \right)_s = -\frac{\tau}{2} \frac{\partial^2 P}{\partial \tau} - \partial_{\tau} P \partial_{\varepsilon} P + P \left[ (\partial_{\varepsilon} P)^2 - 2 \partial_{\tau \varepsilon}^2 P \right] + P^2 \partial_{\varepsilon}^2 P. \qquad (37)$$

When we consider a single phase with the invertible Heuzé potential [15], we have $\partial_{\varepsilon\varepsilon}^2 P = \partial_{\tau\varepsilon}^2 P = 0$, and $\mathcal{G}$ simplifies:

$$\mathcal{G}(\tau, \varepsilon) = -\frac{\tau}{2} \frac{\partial_{\tau}^2 P - \partial_{\tau} P \partial_{\varepsilon} P + P(\partial_{\varepsilon} P)^2}{\partial_{\tau} P - P \partial_{\varepsilon} P}. \qquad (38)$$

**Two-phase zone.** The derivative in (35) is approximated by a second order approximation as

$$c^2(\tau, \varepsilon) = -\tau^2 \left( \frac{P(\tau + \delta\tau, \varepsilon + \delta\varepsilon) - P(\tau - \delta\tau, \varepsilon - \delta\varepsilon)}{2\delta\tau} \right) + \mathcal{O}((\delta\tau)^2)$$

where $\delta\varepsilon = -P\delta\tau$, and $\delta\tau/\tau = 10^{-3}$.

In the same way, the calculation of the fundamental derivative $\mathcal{G}$ uses the following formula for the second-order derivative:

$$\left(\frac{\partial^2 P}{\partial \tau^2}\right)_s = \frac{P(\tau + \delta\tau, \varepsilon + \delta\varepsilon) + P(\tau - \delta\tau, \varepsilon - \delta\varepsilon) - 2P(\tau, \varepsilon)}{\delta\tau^2} + \mathcal{O}((\delta\tau)^2).$$

## 3. DEFINITION AND ASSEMBLY OF THE THREE SINGLE-PHASE EOS

### 3.1. Choice of a Mie-Grüneisen formulation for each phase

We have chosen to model each phase of our three-phase EOS with a same Mie-Grüneisen formulation, which writes in incomplete EOS form as

$$P(\tau, \varepsilon) = P_{\text{ref}}(\tau) + \frac{\Gamma(\tau)}{\tau}(\varepsilon - \varepsilon_{\text{ref}}(\tau)), \qquad (39)$$

where
* $\varepsilon_{\text{ref}}(\tau)$ and $P_{\text{ref}}(\tau)$ are reference curves,

* $\Gamma(\tau) = -\tau(\partial_{\varepsilon}P)_{\tau}$ is the Grüneisen coefficient which is a function only of $\tau$.

We consider the case where the two following assumptions are satisfied:

* the fraction $\Gamma(\tau)/\tau = \frac{\Gamma_0}{\tau_0}$ is constant,

* the reference curve is an isentrope i.e. $\varepsilon'_{\text{ref}} = -P_{\text{ref}}$.

The processes of expansion of the incomplete Mie-Grüneisen EOS allow to access to temperature and entropy [15, 25]. Using Heuzé's formalism, let us introduce the Debye temperature

$$ \Theta(\tau) = \Theta_0 e^{\frac{\Gamma_0}{\tau_0}(\tau_0 - \tau)} $$

where $\Theta_0$ is a constant of the material. The complete EOS then writes

$$ s(\tau, \varepsilon) = s_r + C_{v0} \ln \left( u_r + \frac{\varepsilon - \varepsilon_{\text{ref}}(\tau)}{C_{v0}\Theta(\tau)} \right) \quad (40) $$

$$ T(\tau, \varepsilon) = u_r \Theta(\tau) + \frac{1}{C_{v0}} (\varepsilon - \varepsilon_{\text{ref}}(\tau)) \quad (41) $$

where

* $(\tau_0, \varepsilon_0)$ is the reference state, and $P_0 = P(\tau_0, \varepsilon_0)$, $T_0 = T(\tau_0, \varepsilon_0)$, $s_0 = s(\tau_0, \varepsilon_0)$. The reference curve passes through the reference state $\varepsilon_{\text{ref}}(\tau_0) = \varepsilon_0$ and $P_{\text{ref}}(\tau_0) = P_0$.

* The constants of integration are defined as $u_r = T_0/\Theta_0$ and $s_r = s_0 - C_{v0} \ln(u_r)$.

### 3.1.1. Use of Heuzé's invertible potential

We choose to use the Heuzé's invertible potential [17], defined by

$$ P_{\text{ref}}(\tau) = \frac{K_0}{N_0 + 1} \left( e^{\frac{N_0+1}{\tau_0}(\tau_0 - \tau)} - 1 \right) \quad (42) $$

$$ \varepsilon_{\text{ref}}(\tau) = \varepsilon_0 - \int_{\tau_0}^{\tau} P_{\text{ref}}(\tilde{\tau}) d\tilde{\tau} = \frac{K_0 \tau_0}{(N_0 + 1)} \left( \frac{e^{\frac{N_0+1}{\tau_0}(\tau_0 - \tau)} - 1}{N_0 + 1} - \frac{\tau_0 - \tau}{\tau_0} \right) \quad (43) $$

where

* $K_0$ is the isothermal bulk modulus i.e. $K_0 = -\tau(\partial_{\tau}P)_T$ at $(\tau_0, \varepsilon_0)$.

* $N_0$ is the derivative of the bulk modulus $K_0$ versus pressure.

**Remark 3.1.** Let us list some properties of the Heuzé's invertible potential:

* The pressure of the reference state is always zero i.e. $P_0 = P_{\text{ref}}(\tau_0) = 0$.

* The strict decrease of $P_{\text{ref}}$ implies the bijectivity of the function $(\tau, \varepsilon) \mapsto (P(\tau, \varepsilon), T(\tau, \varepsilon))$. This property is necessary for the construction of the phase diagram.

* The positivity of $\frac{\Gamma_0}{\tau_0}$ and the convexity of the potential $P_{\text{ref}}$ implies the convexity of the EOS i.e. for all isentrope curve $\tau \mapsto \varepsilon_s$, the pressure along the isentrope curve $P(\tau, \varepsilon_s)$ is convex.

The definition of the complete EOS of each phase needs a reference state $(\tau_0, \varepsilon_0, s_r, u_r)$ and five material constants $(\Gamma_0, C_{v0}, K_0, N_0, \Theta_0)$.
## 3.2. A set of EOS coefficients for three-phase tin

The five material constants by phase have been calibrated against static and dynamic experimental data [1–4, 23, 24, 30, 33–35] also considered in [7, 29].

The EOS of each phase depends also on a reference state defined by four parameters ($\tau_0, \varepsilon_0, s_r, u_r$).

The phase $\beta$ is our reference phase which means that the reference states of the other phases depend on the $\beta$ phase. The choice of the reference state of the $\gamma$ and $\ell$ phase will determine the value of the triple point as well as the phase-change lines, so they are of paramount importance in creating the equation of state and are based on experimental values.

For $\gamma$ and $\ell$ phases, we make the assumption that the coefficient $u_r$ is equal to that of $\beta$ phase i.e.

$$u_{r,\beta} = u_{r,\gamma} = u_{r,\ell}. \qquad (44)$$

### 3.2.1. Reference state of the reference $\beta$ phase

The reference state is the standard state of the material at $P_{0,\beta} = P_{\beta}(\tau_{0,\beta}, \varepsilon_{0,\beta}) = 0$ GPa and $T_{0,\beta} = T_{\beta}(\tau_{0,\beta}, \varepsilon_{0,\beta}) = 300$ K the ambient temperature. We can deduce $s_r$ and $u_r$ from the others parameters. Now we have completely characterized the $\beta$ phase, let us determine the reference states of the $\gamma$ and $\ell$ phases.

### 3.2.2. Assembly of phases and reference states of the $\gamma$ and $\ell$ phases

In order to ensure a restitution of the phase-change lines between the three phases, the reference state of the $\alpha$ phase, $\alpha \in \{\gamma, \ell\}$, is deduced from a point of the phase-change line ($P_m, T_m$) on which the slope ($\frac{dP}{dT}$)$_{P_m, T_m}$ and the jump of specific volume ($\Delta\tau$)$_{P_m, T_m}$ are known.

We show here how to deduce the values of ($\tau_0, \varepsilon_0, s_r$)$_\alpha$, $\alpha \in \{\gamma, \ell\}$, thanks to the equation of state of $\beta$ phase and the values of parameters

$$(P_m, T_m, (dP/dT)_{P_m, T_m}, (\Delta\tau)_{P_m, T_m}), \quad (K_0, N_0, \Gamma_0, C_{v0}, \Theta_0)_\alpha, \quad \alpha \in \{\gamma, \ell\}.$$

Note that $u_{r,\alpha} = u_{r,\beta}$ is known.

Let us introduce the predecessors of the reference point ($P_m, T_m$) for the $\beta$ and $\alpha \in \{\gamma, \ell\}$ phases, i.e.

$$(P_m, T_m) = (P_{\beta}, T_{\beta})(\tilde{\tau}_{\beta}, \tilde{\varepsilon}_{\beta}) = (P_{\alpha}, T_{\alpha})(\tilde{\tau}_{\alpha}, \tilde{\varepsilon}_{\alpha}).$$

At this point, the equality of the Gibbs free energies and the Clausius-Clapeyron relation write

$$G_{\beta}(\tilde{\tau}_{\beta}, \tilde{\varepsilon}_{\beta}) = G_{\alpha}(\tilde{\tau}_{\alpha}, \tilde{\varepsilon}_{\alpha}), \qquad \left(\frac{dP}{dT}\right)_{P_m, T_m} = \frac{s_{\alpha}(\tilde{\tau}_{\alpha}, \tilde{\varepsilon}_{\alpha}) - s_{\beta}(\tilde{\tau}_{\beta}, \tilde{\varepsilon}_{\beta})}{(\Delta\tau)_{P_m, T_m}}. \qquad (45)$$

The state $(\tilde{\tau}_{\beta}, \tilde{\varepsilon}_{\beta})$ is determined since we have access to the equation of state of the $\beta$ phase. Therefore the values of $(\tilde{\tau}_{\alpha}, s_{\alpha}(\tilde{\tau}_{\alpha}, \tilde{\varepsilon}_{\alpha}), \tilde{\varepsilon}_{\alpha})$ are determined thanks to equalities (45) as follows

$$\tilde{\tau}_{\alpha} = \tilde{\tau}_{\beta} + (\Delta\tau)_{P_m, T_m}, \\ s_{m,\alpha} = s_{\alpha}(\tilde{\tau}_{\alpha}, \tilde{\varepsilon}_{\alpha}) = s_{\beta}(\tilde{\tau}_{\beta}, \tilde{\varepsilon}_{\beta}) + (\Delta\tau)_{P_m, T_m} \left(\frac{dP}{dT}\right)_{P_m, T_m}, \\ \tilde{\varepsilon}_{\alpha} = \tilde{\varepsilon}_{\beta} - P_m (\Delta\tau)_{P_m, T_m} + T_m (\Delta\tau)_{P_m, T_m} \left(\frac{dP}{dT}\right)_{P_m, T_m}.$$
Finally $(\tau_0, \varepsilon_0, s_r)_\alpha, \alpha \in \{\gamma, \ell\}$, are determined by the resolution of the following system with a Newton-Raphson method:

$$ (P, T, s)_\alpha(\tilde{\tau}_\alpha, \tilde{\varepsilon}_\alpha) = (P, T, s)(\tilde{\tau}_\alpha, \tilde{\varepsilon}_\alpha; \tau_{0,\alpha}, \varepsilon_{0,\alpha}, s_{r,\alpha}) = (P_m, T_m, s_{m,\alpha}). \quad (46) $$

### 3.2.3. Resulting values of the parameters

We now provide a full set of parameters for the three-phase equation of state of tin, see table 1. The calibration procedure, which relies on nested optimization algorithms, is detailed in [31].

The parameters of the nine first lines of table 1 are independent of the four last lines, with values extracted from [22] for phase $\beta$. The value of the triple point

$$ (P, T)^{\text{t}} = (3 \text{ GPa}, 580 \text{ K}) $$

has been imposed at machine precision. Note that the values of phase-change points $(P_m, T_m)_\gamma = (9.4 \text{ GPa}, 300 \text{ K})$ and $(P_m, T_m)_\ell = (0 \text{ GPa}, 505 \text{ K})$ directly appear in the set of coefficients.

<table><caption>TABLE 1. Parameters of the three-phase EOS using Heuzé's invertible potential for each phase.</caption><thead><tr><th colspan="2">Phase</th><th>β</th><th>γ</th><th>liquid</th></tr></thead><tbody><tr><td colspan="2">Phase ref</td><td>-</td><td>β</td><td>β</td></tr><tr><td>K<sub>0</sub></td><td>GPa</td><td>53.5971</td><td>34.2332</td><td>45.5888</td></tr><tr><td>N<sub>0</sub></td><td></td><td>5.53799</td><td>6.26166</td><td>5.22616</td></tr><tr><td>Γ<sub>0</sub></td><td></td><td>2.27</td><td>1.9182</td><td>2.32028</td></tr><tr><td>C<sub>v0</sub></td><td>J/kg/K</td><td>225</td><td>228.166</td><td>224.912</td></tr><tr><td>Θ<sub>0</sub></td><td>K</td><td>250.</td><td>300.</td><td>350.</td></tr><tr><td>T&#773;<sub>m</sub></td><td>K</td><td>-</td><td>300.</td><td>505.</td></tr><tr><td>P<sub>m</sub></td><td>GPa</td><td>-</td><td>9.4</td><td>0.</td></tr><tr><td>Δτ</td><td>m<sup>3</sup>/kg</td><td>-</td><td>-2.23 × 10<sup>-6</sup></td><td>3.86988 × 10<sup>-6</sup></td></tr><tr><td>dP/dT</td><td>Pa/K</td><td>-</td><td>-1.46437 × 10<sup>7</sup></td><td>3.39364 × 10<sup>7</sup></td></tr><tr><td>ρ<sub>0</sub> = 1/τ<sub>0</sub></td><td>kg/m<sup>3</sup></td><td>7287.</td><td>7114.45</td><td>7033.36</td></tr><tr><td>ε<sub>0</sub></td><td>J/kg</td><td>0.</td><td>32491.3</td><td>94742.7</td></tr><tr><td>s<sub>r</sub></td><td>J/kg/K</td><td>-38.2875</td><td>42.6973</td><td>172.852</td></tr><tr><td>u<sub>r</sub></td><td></td><td>1.2</td><td>1.2</td><td>1.2</td></tr></tbody><tfoot><tr><td colspan="5" style="text-align: right;">ref states</td></tr></tfoot></table>

In complement we provide in figure 3 some illustrations of the calibration with respect to experimental data. Comparisons can be done with recent publications [29] and [7]. We also present in figure 4 the diagrams with isovalues for speed of sound and fundamental derivative. One can check in particular the very good behavior near phase-change lines of both quantities.

More numerical results are provided in [31], such as thermodynamic consistency tests and numbers of iterations of the Newton-Raphson algorithm with respect to the size of the tabulation.

## 4. ANALYTICAL AND NUMERICAL RESOLUTION OF RIEMANN PROBLEMS

As usual with phase transition [9, 18, 26], this equation of state has special features so that the characteristic fields are neither genuinely nonlinear nor linearly degenerate. This leads to composite waves, that we describe here in some particular, but representative, cases.
<figure><img src="image_4.png" /><figcaption>FIGURE 3. EOS calibration to experimental data: (A) Phase diagram with experimental triple point and $\beta/\gamma$ and $\beta/\ell$ transitions, and two other calibration points used on the Hugoniot, (B) Isobaric expansion, (C), (D) Isotherm 300 K, (E), (F), (G), (H) Principal Hugoniot.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Isovalues in the plane $(\tau, \varepsilon)$ obtained by the hybrid method.</figcaption></figure>
## 4.1. Reference solutions of Riemann problems

### 4.1.1. Euler equations

We are interested in solutions of the one-dimensional Euler equations, written in time-mass Lagrangian coordinates $(t, m)$:

$$ \left\{ \begin{array}{l} \partial_t \tau - \partial_m u = 0, \\ \partial_t u + \partial_m P = 0, \\ \partial_t e + \partial_m (uP) = 0, \end{array} \right. \qquad (47) $$

where $\tau$, $u$ and $e$ respectively are the specific volume, the velocity and the specific total energy, while the pressure $P$ is linked to the other variables by

$$ \left\{ \begin{array}{l} e = \varepsilon + \frac{1}{2}|u|^2, \\ P = P(\tau, \varepsilon), \end{array} \right. \qquad (48) $$

where $\varepsilon$ is the specific internal energy. Let introduce the notation $U = (\tau, u, e)^t$ and $f(U) = (-u, P, Pu)^t$ such as (47) rewrites

$$ \partial_t U + \partial_m f(U) = 0. $$

We note the domain of states $\Omega = \{(\tau, u, e) \in \mathbb{R}^3 | \tau > 0, u \in \mathbb{R}, T > 0\}$. We recall the notation $w = (\tau, \varepsilon)$ defined in section 1.2. We allow ourselves to use $U$ and $(w, u)$ for the vector of unknowns.

### 4.1.2. Wave curve

The purpose of this paragraph is to determine the set of states $U$ that can be connected to $U_0$ by discontinuities and isentropic compression waves, allowing phase transitions.

**Definition 4.1** (Admissible discontinuities). A single discontinuity that belongs to the same characteristic family between two different states $U_a, U_b \in \Omega$ is admissible if the two states verify

* The Rankine-Hugoniot jump relation, i.e.

$$ \exists \sigma(U_a, U_b) \in \mathbb{R} \text{ such as } \sigma(U_a, U_b)(U_a - U_b) = f(U_a) - f(U_b), \qquad (49) $$

where $\sigma(U_a, U_b)$ is the Lagrangian speed of the discontinuity.

* The Oleinik-Liu entropy condition [21], i.e. for all state $\tilde{U} \in \Omega$ verifying the jump relation (49) with $U_b$ such as $\tilde{\tau} \in (\min(\tau_a, \tau_b), \max(\tau_a, \tau_b))$, inequality

$$ \sigma(\tilde{U}, U_b) \leq \sigma(U_a, U_b) \qquad (50) $$

holds.

The rest of the discussion is independent of the velocity $u$. Indeed for all states $U(w, u)$ and $U_0(w_0, u_0)$ that are connected through a discontinuity, the speed $u$ is uniquely determined by $w, w_0$ and $u_0$ thanks to the jump relation (49) (see [11]). Moreover, we will focus on symmetric initial data of Riemann problems, in compressive situations, with fixed initial thermodynamic states $w_0 = (\tau_0, \varepsilon_0)$ such that the tin is in the $\beta$ phase. Using the Galilean invariance of the model, it is equivalent to describe the wave curve associated with wave $-c/\tau$ starting from the left state and moving rightwards or the wave curve associated with wave $+c/\tau$ starting from the right state and moving leftwards. Due to the compressive configuration we investigate, we end up in both cases
with decreasing specific volumes $\tau$ through discontinuities, and the Hugoniot curve $H(w_0)$ crosses two phase-change zones $(\beta/\gamma)$ and $(\gamma/\ell)$.

From the Rankine-Hugoniot jump relation (49) we define

* the Rayleigh line between $U_a$ and $U_b$ as

$$R(w_a, w_b) = \left\{ \tilde{w} \in \mathcal{F} \mid r_{w_b, w_a}(\tilde{w}) := \frac{P(\tilde{w}) - P(w_b)}{\tilde{\tau} - \tau_b} - \frac{P(w_a) - P(w_b)}{\tau_a - \tau_b} = 0 \right\}, \quad (51)$$

* and the Hugoniot curve centered at $U_b$ as

$$H(w_b) = \left\{ \tilde{w} \in \mathcal{F} \mid \tilde{\tau} < \tau_b \text{ and } h_{w_b}(\tilde{w}) := \tilde{\varepsilon} - \varepsilon_b + \frac{1}{2}(P(\tilde{\tau}, \tilde{\varepsilon}) + P(w_b))(\tilde{\tau} - \tau_b) = 0 \right\}. \quad (52)$$

**Remark 4.2.** The Oleinik-Liu entropy condition can be interpreted graphically in the $(\tau, P)$ plane thanks to the Rayleigh line: restricted to $[\tau_a, \tau_b]$, the Rayleigh line $R(\tilde{w}, w_b)$ is below the Rayleigh line $R(w_a, w_b)$.

When the fundamental derivative is positive, classical compressive shock waves are admissible. In the case of negative values, rarefaction shocks occur and it is impossible to construct compressive shock waves which satisfy the Oleinik-Liu criteria. The description of this specific behavior goes back to the note [19] published by Jouguet in 1904, see also [20, 32]. As a result, since negative values of fundamental derivative appear, isentropic compression waves have to be considered. For a precise definition of smooth isentropic waves, we refer to [11] and use the following properties:

**Proposition 4.3** (Isentropic compression wave). Let $U_a, U_b \in \Omega$ two states. If there exists a k-isentropic compression ($k=1$ or $3$) wave between the states $U_a$ and $U_b$ then the Riemann invariants at $U_a$ and $U_b$ are equal i.e.

$$u_b = \begin{cases} u_a - \psi_a(\tau_b) & \text{1-isentropic compression} \\ u_a + \psi_a(\tau_b) & \text{3-isentropic compression} \end{cases} \quad \text{with } \psi_a(\tau) = - \int_{\tau_a}^{\tau} \frac{c}{\tilde{\tau}} d\tilde{\tau} |_{s(w_a)} \\ s(w_b) = s(w_a).$$

Then in the $(\tau, \varepsilon)$-plane, the isentropic compression curve corresponds to the isentrope passing through $w_a$ and will be noted $I(w_a)$.

We defined the wave curve $S(w_0)$ such as the set of thermodynamic states which can be connected to $w_0$ through admissible discontinuities and isentropic compression waves with $\tau < \tau_0$.

We now describe a whole continuum of solutions, parameterized by the difference between the specific volume of the initial state and that of the end state. We do not provide the construction details of these solutions here. Indeed, it would be necessary to develop a really tedious argumentation, far too long for this publication, see for instance [26, 27].

**Single shock in $\beta$ phase:** we note $w_1(w_0)$ the point at the boundary between the $\beta$ phase and the two-phase zone $(\beta/\gamma)$, defined as $w_1 = (\tau_1, \varepsilon_1) := H(w_0) \cap \mathcal{L}_{\beta, \gamma}^{\beta}$. The interval on $\tau$ is $[\tau_1, \tau_0]$. On this interval for all $w \in H(w_0)$, $(w, w_0)$ verifies the discontinuity admissibility condition. Then, the curve $S(w_0)$ is exactly $H(w_0)$ on the interval $[\tau_0, \tau_1]$.

**Composite wave (shock + isentropic compression) in the mixture zone $(\beta/\gamma)$:** we note $w_2(w_0)$ the point at the boundary between the $\gamma$ phase and the two-phase zone $(\beta/\gamma)$, defined as $w_2 = (\tau_2, \varepsilon_2) := I(w_1) \cap \mathcal{L}_{\beta, \gamma}^{\gamma}$. The interval on $\tau$ is $[\tau_2, \tau_1]$. On this interval, one can
<figure><img src="image_6.png" /><figcaption>FIGURE 5. Representation of the wave curve $S(w_0)$ in the $(\tau, P)$ plane. It is composed of 5 parts: $H(w_0)$, $I(w_1)$, $J(w_1)$, $H(w_1)$, and again $H(w_0)$.</figcaption></figure>

check that admissible discontinuities from $w_0$ and any state of $H(w_0)$ cannot appear. It is the same if two shock waves are considered, a first one between $w_0$ and $w_1$ and a second one from $w_1$ and any state of $H(w_1)$. As a result, the admissible solution is composed by a shock wave between $w_0$ and $w_1$ and an isentropic compression wave between $w_1$ and any state of $I(w_1)$. Note that these two waves are separated by a constant state, namely $w_1$.

**Composite wave (shock + isentropic compression + shock) in the $\gamma$ phase:** We note $w_3(w_0)$ the point of intersection between the tangent of the isentrope $I(w_1)$ at point $w_1$ and $H(w_1)$. As a consequence, the derivative of $I(w_1)$ at $w_1$ in the $(\tau, p)$-plane $(\partial_\tau p|_s = -(\frac{c}{\tau})^2)$ is equal to the shock velocity between states $w_1$ and $w_3$.

The solution for a given $\tau \in [\tau_3, \tau_2]$ is composed by a shock wave between $w_0$ and $w_1$, followed by an isentropic compression wave from $w_1$ to some point $w_* \in I(w_1)$, and a shock wave from $w_*$ and $w$, so that $w \in H(w_*)$. The state $w_*$ is the unique point such that $\partial_\tau p|_s(w_*) = \sigma(w_*, w)$, which means that the end of the isentropic compression wave between $w_1$ and $w_*$ is sticked to the shock wave between $w_*$ and $w$.

**Composite wave (shock+shock) in $\gamma$ phase:** We note $w_4(w_0)$ the unique point of intersection between the Rayleigh curve $R(w_1, w_0)$ and the Hugoniot curve $H(w_0)$ such as $\tau_4 < \tau_3$. The interval on $\tau$ is $[\tau_4, \tau_3]$. On this interval the Hugoniot curve $H(w_0)$ is below the Rayleigh line $R(w_1, w_0)$ so for all $w \in H(w_0)$, the pair $(w, w_0)$ does not satisfy the Oleinik–Liu criteria. The admissible solution between $w$ and $w_0$ is composed of the two admissible shocks $(w, w_1)$ and $(w_1, w_0)$. The curves $S(w_0)$ and $H(w_1)$ coincide.

**Single shock in $\gamma$ phase, the mixture zone ($\gamma/\ell$) and $\ell$ phase:** The interval on $\tau$ is $[0, \tau_4]$. As in the first interval, for all $w \in H(w_0)$, $(w, w_0)$ verifies the admissible condition. Then the curves $S(w_0)$ and $H(w_0)$ coincide.

Let us define $w_5(w_0) \in S(w_0)$ the point of entry into the two-phase zone ($\gamma/\ell$) i.e. $w_5 = (\tau_5, \varepsilon_5) := H(w_0) \cap \mathcal{L}_{\gamma, \ell}^{\gamma}$, and $w_6(w_0) \in S(w_0)$ the point of exit from the biphasic zone ($\gamma/\ell$) i.e. $w_6 = (\tau_6, \varepsilon_6) := H(w_0) \cap \mathcal{L}_{\gamma, \ell}^{\ell}$. Then on the interval $[\tau_5, \tau_4]$, $w \in S(w_0)$ is in $\gamma$ phase, on the interval $[\tau_6, \tau_5]$, $w \in S(w_0)$ is in mixture ($\gamma/\ell$) and on the interval $[0, \tau_6]$, $w \in S(w_0)$ is in $\ell$ phase.
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Graph of the wave curve $S(w_0)$ ($w_0$ such as $P(w_0) = 0$ GPa and $T(w_0) = 300$ K) and the 6 characteristic points ($w_k$)$_{1 \le k \le 6}$ in the phase diagram $(P,T)$.</figcaption></figure>

## 4.2. Symmetric compressive Riemann problems

We focus on initial configurations which write as

$$ U(m, 0) = \begin{cases} \tilde{U}_0 = (\tau_0, u, \varepsilon_0 + \frac{1}{2}u^2) & \text{if } m < 0, \\ U_0 = (\tau_0, -u, \varepsilon_0 + \frac{1}{2}u^2) & \text{if } m > 0. \end{cases} \qquad (53) $$

where $u \in \mathbb{R}_+$. We will keep $w_0$ constant and investigate the possible solutions letting $u$ vary.

Let us defined the state $U_i = (\tau_i, u_i, \varepsilon_i)$ which is located at the interface $m = 0$. It is connected to $U_0$ by waves associated with the eigenvalue $c/\tau$ and to $\tilde{U}_0$ by waves associated with the eigenvalue $-c/\tau$. The solution is symmetric, so that $u_i = 0$ and $w_i \in S(w_0)$.

We parameterize the solutions by the initial speed $u$. For the states which are connected by a single shock wave, i.e. states $w_k$ with $k \in \{1, 4, 5, 6\}$, we use the jump relation (49), so that

$$ u_k^2 = 2(\varepsilon_k - \varepsilon_0) - 2P_0(\tau_k - \tau_0). \qquad (54) $$

On the other hand, the initial speed $u_2$ needed to reach the thermodynamic state $w_2$,

$$ u_2 = u_1 + \int_{\tau_1}^{\tau_2} \frac{c}{\tau} d\tau \qquad (55) $$

Using the previous description of $S(w_0)$, one can check that four types of solutions can appear. Let us describe them with respect to the initial velocity $u$.

$0 < u < u_1$: There is no phase-change and $w \in H(w_0)$. The solution consists of three constant states separated by two single shocks.

The solution writes

$$ U(m,t) = \begin{cases} (\tau_0, -\text{sign}(m)u, \varepsilon_0 + \frac{1}{2}u^2)^t & \text{if } \sigma < |m|/t, \\ (\tau_i, 0, \varepsilon_i)^t & \text{if } 0 \le |m|/t < \sigma. \end{cases} \qquad (56) $$
<figure><img src="image_8.png" /><figcaption>FIGURE 7. Reference solution in the (m, t) diagram of a simple shock configuration.</figcaption></figure>

where $\sigma$ is the speed shock between $U_i$ and $U_0$. See Figure 7.

$u_1 < u < u_2$: The symmetric solution reveals two composite waves (shock+isentropic compression) for the 1-wave and the 3-wave as shown in Figure 8.

<figure><img src="image_9.png" /><figcaption>FIGURE 8. Reference solution in the (m, t) diagram of a composite wave (shock+isentropic compression) configuration.</figcaption></figure>

The solution writes

$$ U(m,t) = \begin{cases} (\tau_0, -\operatorname{sign}(m)u, \varepsilon_0 + \frac{1}{2}u^2)^t & \text{if } \sigma(U_0, U_i) < |m|/t, \\ (\tau_1, -\operatorname{sign}(m)(u-u_1), \varepsilon_1 + \frac{1}{2}(u-u_1)^2)^t & \text{if } \frac{c}{\tau}(w_1) < |m|/t < \sigma(U_0, U_i), \\ (\tau_R, -\operatorname{sign}(m)(u-u_1-u_R), \varepsilon_R + \frac{1}{2}(u-u_1-u_R)^2)^t(\frac{m}{t}) & \text{if } \frac{c}{\tau}(w_i) < |m|/t < \frac{c}{\tau}(w_1), \\ (\tau_i, 0, \varepsilon_i)^t & \text{if } 0 < |m|/t < \frac{c}{\tau}(w_i). \end{cases} \tag{57} $$

where $(\tau_R, u_R, \varepsilon_R + \frac{1}{2}u_R^2)^t(\frac{m}{t})$ is the isentropic compression wave between $U_1$ and $U_i$.

$u_2 < u < u_3$: The symmetric solution reveals two composite waves (shock+isentropic compression+shock) as shown in Figure 9.
<figure><img src="image_10.png" /><figcaption>FIGURE 9. Reference solution in the (m, t) diagram of a composite wave (shock+isentropic compression+shock) configuration.</figcaption></figure>

The solution writes

$$ U(m,t) = \begin{cases} (\tau_0, -\operatorname{sign}(m)u, \varepsilon_0 + \frac{1}{2}u^2)^t & \text{if } \sigma_1 < |m|/t, \\ (\tau_1, -\operatorname{sign}(m)(u-u_1), \varepsilon_1 + \frac{1}{2}(u-u_1)^2)^t & \text{if } \frac{c}{\tau}(w_1) < |m|/t < \sigma_1, \\ (\tau_R, -\operatorname{sign}(m)(u-u_1-u_R), \varepsilon_R + \frac{1}{2}(u-u_1-u_R)^2)^t(\frac{m}{t}) & \text{if } \frac{c}{\tau}(w_*) < |m|/t < \frac{c}{\tau}(w_1), \\ (\tau_i, 0, \varepsilon_i)^t & \text{if } 0 \le |m|/t < \frac{c}{\tau}(w_*) \end{cases} \tag{58} $$

where $\sigma_1 = \sigma(U_0, U_1)$, and $(\tau_R, u_R, \varepsilon_R + \frac{1}{2}u_R^2)(\frac{m}{t})^t$ is the isentropic compression wave between $U_1$ and $U_*$, such as $w_* \in I(w_1)$ as described in the description above.

### $u_3 < u < u_4$: The thermodynamic state $w$ belongs to $H(w_1)$ and does not lie in the $\beta$ phase. The solution is a double shock for the left-facing wave and the right-facing wave as shown in Figure 10.

<figure><img src="image_11.png" /><figcaption>FIGURE 10. Reference solution in the (m, t) diagram of a double shock configuration.</figcaption></figure>

The solution writes

$$ U(m,t) = \begin{cases} (\tau_0, -\operatorname{sign}(m)u, \varepsilon_0 + \frac{1}{2}u^2)^t & \text{if } \sigma_1 < |m|/t, \\ (\tau_1, -\operatorname{sign}(m)(u-u_1), \varepsilon_1 + \frac{1}{2}(u-u_1)^2)^t & \text{if } \sigma < |m|/t < \sigma_1, \\ (\tau_i, 0, \varepsilon_i)^t & \text{if } 0 \le |m|/t < \sigma. \end{cases} \tag{59} $$
u_4 < u: The thermodynamic state w belongs to H(w_0). The solution is similar to the case 0 < u < u_1 however the (β/γ) phase transition occur in the front wave. If u > u_6, the (β/γ) and (γ/ℓ) phase transitions occur in the front wave.

4.3. Reference values for the compressive wave curve

The values of the seven notable states of the wave curve that appear in the solution of the specific Riemann problems (53), see Figures 5 and 6, are given in Table 2.

<table><caption>TABLE 2. Reference values of the compressive wave curve.</caption><thead><tr><th></th><th>U<sub>0</sub></th><th>U<sub>1</sub></th><th>U<sub>2</sub></th><th>U<sub>3</sub></th><th>U<sub>4</sub></th><th>U<sub>5</sub></th><th>U<sub>6</sub></th></tr><tr><th>Phase</th><th>β</th><th>β</th><th>γ</th><th>γ</th><th>γ</th><th>γ</th><th>ℓ</th></tr></thead><tbody><tr><th>τ (m<sup>3</sup>/kg)(&times;10<sup>-5</sup>)</th><td>13.723</td><td>12.341</td><td>12.005</td><td>11.991</td><td>11.546</td><td>9.2667</td><td>9.0397</td></tr><tr><th>ε (J/kg)(&times;10<sup>4</sup>)</th><td>7.1673&times;10<sup>-14</sup></td><td>5.3566</td><td>8.0904</td><td>8.2219</td><td>13.29</td><td>129.2</td><td>168.29</td></tr><tr><th>P (GPa)</th><td>0</td><td>7.7492</td><td>8.5197</td><td>8.6216</td><td>12.206</td><td>57.986</td><td>71.866</td></tr><tr><th>T (K)</th><td>300</td><td>405.97</td><td>358.28</td><td>359.34</td><td>402.79</td><td>2853.6</td><td>3289.4</td></tr><tr><th>u (m/s)</th><td>-</td><td>327.31</td><td>378.15</td><td>382.57</td><td>515.55</td><td>1607.5</td><td>1834.6</td></tr><tr><th>s (J/kg/K)</th><td>2.7349</td><td>19.345</td><td>19.345</td><td>19.57</td><td>31.753</td><td>407.53</td><td>486.71</td></tr></tbody></table>

The values of speed of sound and fundamental derivative are provided in table 3. Since these two quantities are discontinuous through phase-change lines, two values must be given: one for the pure phase and one for the mixture. Moreover, the fundamental derivative G is negative in the mixture zone (β/γ), see Table 4, but remains positive in the mixture zone (γ/ℓ), see Table 5.

<table><caption>TABLE 3. Speed of sound and fundamental derivative.</caption><thead><tr><th></th><th></th><th>U<sub>0</sub></th><th>U<sub>1</sub></th><th>U<sub>2</sub></th><th>U<sub>3</sub></th><th>U<sub>4</sub></th><th>U<sub>5</sub></th><th>U<sub>6</sub></th></tr><tr><th>Phase</th><th></th><th>β</th><th>β</th><th>γ</th><th>γ</th><th>γ</th><th>γ</th><th>ℓ</th></tr></thead><tbody><tr><th>c (m/s) (pure phase)</th><td></td><td>2712</td><td>3394</td><td>3170</td><td>3179</td><td>3437.6</td><td>5061.5</td><td>5141.2</td></tr><tr><th>c (m/s) (mixture of phases)</th><td></td><td>-</td><td>1948</td><td>1747</td><td>-</td><td>-</td><td>4797</td><td>5059</td></tr><tr><th>G (pure phase)</th><td></td><td>2.134</td><td>1.914</td><td>3.117</td><td>3.114</td><td>2.993</td><td>2.339</td><td>1.926</td></tr><tr><th>G (mixture of phases)</th><td></td><td>-</td><td>-1.019</td><td>-0.417</td><td>-</td><td>-</td><td>1.125</td><td>1.041</td></tr></tbody></table>

For the set of numerical values in Tables 2, 3, 4 and 5, the points linked by one or several discontinuities to w₀, i.e. w₁, w₃, w₄, w₅ and w₆ are computed up to machine precision by the Newton-Raphson method thanks to the Rankine-Hugoniot relations. The isentrope which starts from the point I(w₁) necessary to compute the point w₂ and all points in the mixture zone (β/γ) is discretized with a step of specific volume given by Δτ/τ = 10⁻⁶.

4.4. Basic 1D Lagrangian Godunov-type scheme

We consider a space domain discretized initially by the interfaces $(x_{i+1/2}^0)_{i \in \mathbb{Z}}$ defined as $x_{i+1/2}^0 = i\Delta x$. The interfaces of the mass domain $(m_{i+1/2})_{i \in \mathbb{Z}}$ are defined by $m_{i+1/2} = i\Delta x/\tau_0(m_{i+1/2})$ which are constant over time. We note the cell mass $\Delta m_i = m_{i+1/2} - m_{i-1/2}$. (In the case we are interested in, $\tau_0$ is constant over the mass domain.)
<table><caption>TABLE 4. Mixture zone (β/γ) on the isentropic compression wave.</caption><thead><tr><th>λ<sub>γ</sub></th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th><th>0.6</th><th>0.7</th><th>0.8</th><th>0.9</th></tr></thead><tbody><tr><th>τ (m<sup>3</sup>/kg) (×10<sup>-5</sup>)</th><td>12.31</td><td>12.27</td><td>12.24</td><td>12.21</td><td>12.17</td><td>12.14</td><td>12.11</td><td>12.07</td><td>12.04</td></tr><tr><th>ε (J/kg) (×10<sup>4</sup>)</th><td>5.618</td><td>5.883</td><td>6.150</td><td>6.420</td><td>6.693</td><td>6.967</td><td>7.245</td><td>7.525</td><td>7.806</td></tr><tr><th>P (GPa)</th><td>7.832</td><td>7.914</td><td>7.994</td><td>8.073</td><td>8.150</td><td>8.227</td><td>8.302</td><td>8.376</td><td>8.448</td></tr><tr><th>T (K)</th><td>401.0</td><td>396.1</td><td>391.2</td><td>386.4</td><td>381.6</td><td>376.8</td><td>372.1</td><td>367.5</td><td>362.8</td></tr><tr><th>u (m/s)</th><td>332.6</td><td>337.8</td><td>343.0</td><td>348.2</td><td>353.3</td><td>358.3</td><td>363.4</td><td>368.3</td><td>373.3</td></tr><tr><th>s (J/kg/K)</th><td>19.34</td><td>19.34</td><td>19.34</td><td>19.34</td><td>19.34</td><td>19.34</td><td>19.34</td><td>19.34</td><td>19.34</td></tr><tr><th>c (m/s)</th><td>1926</td><td>1904</td><td>1883</td><td>1862</td><td>1842</td><td>1822</td><td>1803</td><td>1784</td><td>1766</td></tr><tr><th>G</th><td>-0.93</td><td>-0.85</td><td>-0.77</td><td>-0.71</td><td>-0.65</td><td>-0.59</td><td>-0.54</td><td>-0.49</td><td>-0.45</td></tr></tbody></table>

<table><caption>TABLE 5. Mixture zone (γ/ℓ) on the principle Hugoniot.</caption><thead><tr><th>λ<sub>ℓ</sub></th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th><th>0.6</th><th>0.7</th><th>0.8</th><th>0.9</th></tr></thead><tbody><tr><th>τ (m<sup>3</sup>/kg) (×10<sup>-5</sup>)</th><td>9.244</td><td>9.221</td><td>9.198</td><td>9.175</td><td>9.153</td><td>9.13</td><td>9.108</td><td>9.085</td><td>9.062</td></tr><tr><th>ε (J/kg) (×10<sup>4</sup>)</th><td>133.1</td><td>137</td><td>140.9</td><td>144.8</td><td>148.7</td><td>152.6</td><td>156.5</td><td>160.4</td><td>164.3</td></tr><tr><th>P (GPa)</th><td>59.42</td><td>60.85</td><td>62.26</td><td>63.66</td><td>65.05</td><td>66.43</td><td>67.80</td><td>69.16</td><td>70.52</td></tr><tr><th>T (K)</th><td>2903</td><td>2950</td><td>2997</td><td>3042</td><td>3086</td><td>3129</td><td>3171</td><td>3211</td><td>3251</td></tr><tr><th>u (m/s)</th><td>1592</td><td>1614</td><td>1636</td><td>1658</td><td>1679</td><td>1700</td><td>1721</td><td>1741</td><td>1762</td></tr><tr><th>s (J/kg/K)</th><td>416.3</td><td>424.9</td><td>433.3</td><td>441.5</td><td>449.5</td><td>457.2</td><td>464.9</td><td>472.3</td><td>479.6</td></tr><tr><th>c (m/s)</th><td>4834</td><td>4868</td><td>4899</td><td>4928</td><td>4956</td><td>4980</td><td>5003</td><td>5024</td><td>5043</td></tr><tr><th>G</th><td>1.125</td><td>1.125</td><td>1.121</td><td>1.116</td><td>1.108</td><td>1.098</td><td>1.087</td><td>1.073</td><td>1.058</td></tr></tbody></table>

The initial solution is $U_i^0 = \frac{1}{m_i} \int_{m_{i-1/2}^0}^{m_{i+1/2}^0} U_0(\tilde{m}) d\tilde{m}$ with $U_0$ the initial condition (53). Then at time $t^{n+1} = t^n + \Delta t^n$, the approximation $\Delta m_i U_i^n$ of the integral of $U(m, t^n)$ verifies

$$U_i^{n+1} = U_i^n - \frac{\Delta t^n}{\Delta m_i} (F_{i+1/2}^n - F_{i-1/2}^n) \quad (60)$$

where $F_{i+1/2}^n = (u^*, P^*, u^*P^*)_{i+1/2}^t$ is the numeral flux.

The scheme consists of the five steps.

1. **EOS call:** using the multiphase EOS, the pressure and Lagrangian speed of sound are computed according to

$$P_i^n = P(\tau_i^n, \varepsilon_i^n), \\ (\rho c)_i^n = (\rho c)(\tau_i^n, \varepsilon_i^n).$$

2. **Computation of the numerical flux at interfaces:** the approximations of the pressure $P_{i+1/2}^*$ and material velocity $u_{i+1/2}^*$ at interface $m_{i+1/2}$ and time $t^n$ are computed by linearizing the Riemann invariants. Along the characteristic curves $C_\pm$ defined by $dm = \pm(\rho c)dt$, relations $dp \pm (\rho c)du = 0$ lead to

$$(P_{i+1/2}^* - P_i^n) + (\rho c)_i^n (u_{i+1/2}^* - u_i^n) = 0, \\ (P_{i+1/2}^* - P_{i+1}^n) + (\rho c)_{i+1}^n (u_{i+1/2}^* - u_{i+1}^n) = 0.$$
This linear system gives the flux formula, the so-called Godunov's acoustic Lagrangian Riemann solver [12]

$$u_{i+1/2}^* = \frac{(\rho c)_i^n u_i^n + (\rho c)_{i+1}^n u_{i+1}^n}{(\rho c)_i^n + (\rho c)_{i+1}^n} + \frac{1}{(\rho c)_i^n + (\rho c)_{i+1}^n} (P_i^n - P_{i+1}^n),$$

$$p_{i+1/2}^* = \frac{(\rho c)_{i+1}^n P_i^n + (\rho c)_i^n P_{i+1}^n}{(\rho c)_i^n + (\rho c)_{i+1}^n} + \frac{(\rho c)_i^n (\rho c)_{i+1}^n}{(\rho c)_i^n + (\rho c)_{i+1}^n} (u_i^n - u_{i+1}^n).$$

3. **Computation of the time step:** the time step is computed according to the CFL stability condition

$$\Delta t^n = CFL \min_i \frac{\Delta m_i}{((\rho c)_i^n)}, \quad \text{with } 0 < CFL \le 1/2. \qquad (61)$$

4. **Computation of the solution at time $t^{n+1}$:**

$$\tau_i^{n+1} = \tau_i^n + \frac{\Delta t^n}{\Delta m_i}(u_{i+1/2}^* - u_{i-1/2}^*),$$

$$u_i^{n+1} = u_i^n - \frac{\Delta t^n}{\Delta m_i}(P_{i+1/2}^* - P_{i-1/2}^*),$$

$$e_i^{n+1} = e_i^n - \frac{\Delta t^n}{\Delta m_i} (P_{i+1/2}^* u_{i+1/2}^* - P_{i-1/2}^* u_{i-1/2}^*).$$

5. **Actualisation of the grid:**

$$x_{i+1/2}^{n+1} = x_{i+1/2}^n + \Delta t^n u_{i+1/2}^*.$$

The CFL condition (61) in step 2 is rather rough. It could be much refined to guarantee the positivity of the density and the specific energy, and discrete entropy inequalities is sufficient to recover solutions selected with the Oleinik-Liu criteria, but in all the numerical tests we performed, it turns out that the CFL condition (61) has been sufficient to recover the expected solutions, even in presence of composite waves.

### 4.5. Examples of numerical simulations

These test cases are obtained using initial data (53) for specific values of the initial material velocity $u > 0$. The domain is centered in zero and large of 10 mm. Initially, the material is in the $\beta$ phase, with $P_0 = P(w_0) = 0$ GPa and $T_0(w_0) = 300$ K and the initial configuration writes as

$$U(m, 0) = \begin{cases} (\tau_0, u, \varepsilon_0 + \frac{1}{2}u^2) & \text{if } -5.10^{-3} \le x \le 0, \\ (\tau_0, -u, \varepsilon_0 + \frac{1}{2}u^2) & \text{if } 0 < x \le 5.10^{-3}. \end{cases} \qquad (62)$$

The time step is computed with $CFL = 0.45$ in (61).
4.5.1. *Shock+shock configuration*

The initial material velocity $\mathbf{u} = 400 \ \mathrm{m/s}$ is such that $u_3 < \mathbf{u} < u_4$. The two waves of the symmetric solution are both composed of two shocks. We display the reference solution obtained above along with the numerical results using meshes of 250, 1000 and 4000 cells at the final time $1.2 \ \mu\mathrm{s}$, see Figures 11 and 12. The dissipation of this first-order Godunov-type scheme is significant

<figure><img src="image_12.png" /><figcaption>FIGURE 11. Pressure profiles for a $\gamma$ phase shock in a double shock configuration using the Godunov-type scheme. Top: full domain. Bottom: zoom of the left-facing composite wave.</figcaption></figure>

so 250 cells is not enough to capture accurately the intermediate $w_1$ level. With the increase in the number of cells, the simulation approaches the reference solution. The phase transition between the $\beta$ phase and the $\gamma$ phase is clearly visible through the shock front in Figure 12.

4.5.2. *Shock+isentropic compression configuration*

With the initial speed $\mathbf{u} = u_2 = 378.15 \ \mathrm{m/s}$, the two waves of the symmetric solution are composed of a shock followed by an isentropic compression up to full $\gamma$ phase transition.
<figure><img src="image_13.png" /><figcaption>FIGURE 12. Mass fractions profiles of the three phases for a $\gamma$ phase shock in a double shock configuration using the Godunov-type scheme.</figcaption></figure>

We display in Figures 13 and 14 the reference solution at the final time 1.2 µs and the numerical results using meshes of 250, 1000 and 4000 cells. With 250 cells, the mesh is not fine enough to detect if the profile of the second wave is an isentropic compression or a shock. The profile of the isentropic compression wave clearly appears only for 4000 cells. The (smooth) phase transition between $\beta$ and $\gamma$ phases occurs only through the isentropic compression waves, see Figure 14.

### 4.5.3. Self-similar regimes of the $\beta/\gamma$ phase transition

Using previous solutions and Galilean invariance we plot six different solutions corresponding to increasing velocities $u$ in relation to the notable speeds listed in Table 2. To display the different types of solutions on a single figure, we shift the profiles such that the fluid is initially at rest.

Due to the non-convexity of the EOS and negative values of the fundamental derivative, the various self-similar regimes of the ($\beta \to \gamma$) transition appear for increasing velocity $u$:

(a) $u = u_1 = 327.3 \text{ m/s}$: single shock wave in $\beta$ phase.

(b) $u_1 = 327.3 \text{ m/s} < u = 368.3 \text{ m/s} < u_2 = 378.1 \text{ m/s}$: shock wave in $\beta$ phase + $\beta/\gamma$ isentropic compression.

(c) $u_2 = 378.3 \text{ m/s} < u = (u_2 + u_3)/2 < u_3 = 382.5 \text{ m/s}$: shock wave in $\beta$ phase + $\beta/\gamma$ isentropic compression + shock wave from $\beta/\gamma$ region to $\gamma$ phase.

(d) $u = u_3 = 515.5 \text{ m/s}$: shock wave in $\beta$ phase + shock wave from $\beta$ phase to $\gamma$ phase.
<figure><img src="image_14.png" /><figcaption>FIGURE 13. Pressure profiles using the Godunov-type scheme in a shock+isentropic compression case. Top: full domain. Bottom: zoom on the left-facing composite wave.</figcaption></figure>

(e) $u_3 = 382.5 \text{ m/s} < u = (u_3 + u_4)/2 < u_4 = 515.5 \text{ m/s}$:
shock wave in $\beta$ phase + shock wave from $\beta$ phase to $\gamma$ phase.

(f) $u_4 = 515.5 \text{ m/s} < u = 1.1 \times u_4$:
single shock wave from $\beta$ phase to $\gamma$ phase.

Even with a fine 2 $\mu$m space resolution, the first-order Godunov-type scheme does not provide sufficiently accurate profiles and the successive configurations are difficult to distinguish. However, using a second-order accurate GAD Godunov-type scheme with a minmod limiter, see [16] for more details, the different wave configurations are well reproduced, see figure 15.

Due to the negative values of fundamental derivative in the mixture $\beta/\gamma$ at equilibrium, the solution (shock wave+isentropic compression) appears on a large range of velocities of 50.8 m/s. On the other hand, the solutions of type (b) (shock+isentropic compression+shock) appears on a very small range of velocities of 4.4 m/s.
<figure><img src="image_15.png" /><figcaption>FIGURE 14. Mass fractions profiles of the 3 phases of the numerical result using a mesh of 4000 cells.</figcaption></figure>

Note that quite significant numerical dissipation can be observed at the foot of the isentropic compression. This is in part due to the large difference of wave speed between the first shock wave and the isentropic compression wave (recall that sound speeds in mixtures are very low compared to their pure phase values). This kind of solutions with composite waves represents a real numerical challenge and constitute valuable benchmarks for hydrodynamic schemes well-adapted to multiphase EOS for metals with phase transition.

## 5. CONCLUSION

In this article, we have discussed the modeling of shock-induced phase transitions focusing on the solid phase $\beta$, solid phase $\gamma$ and liquid phase of tin: a three-phase equation of state in *strict* thermodynamic equilibrium has been constructed, self-similar reference solutions established, illustrated by numerical simulations of compressive configurations.

We have introduced hypothesis 1.2, which guarantees the existence of a triple point in the $(P, T)$ phase diagram. Using a complete description of the phase diagram in the $(\tau, \epsilon)$ plane, a hybrid numerical method (tabulation+Newton-Raphson iterations) has been proposed, enabling the computation of the thermodynamic equilibrium at *machine precision*.

Using a Mie-Grüneisen formalism with Heuzé's invertible potential for each phase, the parameters of each equation of state have been calibrated with static and dynamic experimental data of the literature. Associated diagrams for the speed of sound and the fundamental derivative have been provided, showing the very good behavior of these variables through the phase-change lines. One can also check that the fundamental derivative becomes negative in the $(\beta/\gamma)$ mixture zone.

By a careful study of the wave curves, a set of solutions has been described in symmetric compressive cases, involving shock waves and isentropic compression waves, with phase transitions. The negative values of the fundamental derivative lead to complex successions of waves, providing challenging numerical benchmarks for hydrodynamic codes. First- and second-order accurate Godunov-type Lagrangian schemes have finally been used, showing their ability to capture the physically admissible solutions.
<figure><img src="image_16.png" /><figcaption>FIGURE 15. The five self-similar regimes at thermodynamic equilibrium of shock-induced ($\beta \to \gamma$) phase transition. Top: pressure profiles at $t = 1.2 \ \mu s$ with space resolution 2 $\mu m$ (GAD scheme, CFL= 0.3). Bottom: zoom of the above figure.</figcaption></figure>

We should also mention that this work has been the subject of more in-depth studies, see [31]. In particular, several experimental measurements of high-speed impacts point to the presence of phases that are not at thermodynamic equilibrium over non-negligible durations. It is therefore necessary to introduce kinetic phase transitions to model this behavior. Again, promising results in this direction are described in [31], using relaxation models on mass fractions equations.

**Acknowledgments.** The second author (H.J.) thanks Olivier Heuzé for the numerous fruitful discussions over the years.

## REFERENCES

[1] Boris Batokovich Alchagirov and Asiyat Maskhutovna Chochaeva. Temperature dependence of the density of liquid tin. *High Temperature*, 38(1):44–48, 2000.
[2] L. V. Al'Tshuler, A. A. Bakanova, I. P. Dudoladov, E. A. Dynin, R. F. Trunin, and B. S. Chekin. Shock adiabatic curves of metals: New data, statistical analysis, and general laws. *Journal of Applied Mechanics and Technical Physics*, 22(2):145-169, 1981.

[3] L. V. Al'Tshuler, A. A. Bakanova, and R. F. Trunin. Shock adiabats and zero isotherms of seven metals at high pressures. *Sov. Phys. JETP*, 15(1):65-75, 1962.

[4] L. V. Al'Tshuler, K. Krupnikov, B. N. Ledenev, V. I. Zhuchikhin, and M. I. Brazhnik. Dynamic compressibility and equation of state of iron under high pressure. *Soviet Physics JETP*, 34(7), 1958.

[5] François Bouchut. *Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources*. Frontiers in Mathematics. Birkhäuser Verlag, Basel, 2004.

[6] H. B. Callen. *Thermodynamics and an introduction to thermostatistics*. second ed., John Wiley and Sons, 1985.

[7] G. A. Cox and M. A. Christie. Fitting of a multiphase equation of state with swarm intelligence. *Journal of Physics: Condensed Matter*, 27(40):405201, sep 2015.

[8] G. E. Duvall and R. A. Graham. Phase transitions under shock-wave loading. *Rev. Mod. Phys.*, 49:523-579, 1977.

[9] Gloria Faccanoni. *Étude d'un modèle fin de changement de phase liquide-vapeur. Contribution à l'étude de la crise d'ébullition*. Phd thesis, Ecole Polytechnique X, November 2008. Thèse en cotutelle entre l'École Polytechnique et l'Université di Trento (Italie).

[10] Gloria Faccanoni, Samuel Kokh, and Grégoire Allaire. Modelling and simulation of liquid-vapor phase transition in compressible flows based on thermodynamical equilibrium. *ESAIM: M2AN*, 46(5):1029-1054, 2012.

[11] Edwige Godlewski and Pierre-Arnaud Raviart. *Numerical approximation of hyperbolic systems of conservation laws*, volume 118 of *Applied Mathematical Sciences*. Springer-Verlag, New York, 1996.

[12] S. K. Godunov. A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics. *Mat. Sb. (N.S.)*, 47 (89):271-306, 1959.

[13] D. B. Hayes. Wave propagation in a condensed medium with N transforming phases: Application to solid-I-solid-II-liquid bismuth. *J. Appl. Phys.*, 46(8):3438-3443, 1975.

[14] Philippe Helluy and Nicolas Seguin. Relaxation models of phase transition flows. *ESAIM: M2AN*, 40(2):331-352, 2006.

[15] Olivier Heuzé. General form of the Mie-Grüneisen equation of state. *Comptes Rendus. Mécanique*, 340(10):679-687, 2012.

[16] O. Heuzé, S. Jaouen, and H. Jourdren. Dissipative issue of high-order shock capturing schemes with non-convex equations of state. *Journal of Computational Physics*, 228:833-860, 2009.

[17] Olivier Heuzé and Damian Swift. Analysis and modeling of laser ramps and shocks in titanium and zirconium with phase transitions. *AIP Conference Proceedings*, 1426(1):1541-1544, 03 2012.

[18] Stéphane Jaouen. *Étude mathématique et numérique de stabilité pour des modèles hydrodynamiques avec transition de phase*. Phd thesis, Université Pierre et Marie Curie, 2001.

[19] Émile Jouguet. Remarques sur la loi adiabatique d'Hugoniot. *Comptes rendus de l'Académie des Sciences*, 139:786-789, 1904.

[20] L.D. Landau. On shock waves. *J. Phys. USSR*, 6, 1942.

[21] Tai-Ping Liu. The entropy condition and the admissibility of shocks. *Journal of Mathematical Analysis and Applications*, 53(1):78-88, 1976.

[22] Catherine Mabire. *Transformation polymorphique et fusion de l'étain sous choc dans la gamme 0-100 GPa. Etude expérimentale et modélisation*. PhD thesis, Université de Poitiers, 1999.

[23] Stanley P. Marsh. *LASL shock Hugoniot data*, volume 5. Univ of California Press, 1980.

[24] R. G. McQueen and S. P. Marsh. Equation of State for Nineteen Metallic Elements from Shock-Wave Measurements to Two Megabars. *Journal of Applied Physics*, 31(7):1253-1269, 07 1960.

[25] Ralph Menikoff. Complete Mie-Grüneisen equation of state. Technical report, Los Alamos National Laboratory (LANL), Los Alamos, NM (United States), 2016.

[26] Ralph Menikoff and Bradley J. Plohr. The Riemann problem for fluid flow of real materials. *Rev. Mod. Phys.*, 61:75-130, 1989.

[27] Christian Merkle and Christian Rohde. The sharp-interface approach for fluids with phase change: Riemann problems and ghost fluid techniques. *M2AN Math. Model. Numer. Anal.*, 41(6):1089-1123, 2007.

[28] Andrés Piña-Martinez, Romain Privat, and Jean-Noël Jaubert. Use of 300,000 pseudo-experimental data over 1800 pure fluids to assess the performance of four cubic equations of state: SRK, PR, tc-RK, and tc-PR. *AIChE Journal*, 68(2):e17518, 2022.
[29] Daniel A. Rehn, Carl W. Greeff, Leonid Burakovsky, Daniel G. Sheppard, and Scott D. Crockett. Multiphase tin equation of state using density functional theory. Phys. Rev. B, 103:184102, May 2021.

[30] Ashkan Salamat, Richard Briggs, Pierre Bouvier, Sylvain Petitgirard, Agnès Dewaele, Melissa E. Cutler, Furio Corà, Dominik Daisenberger, Gaston Garbarino, and Paul F. McMillan. High-pressure structural transformations of sn up to 138 gpa: Angle-dispersive synchrotron x-ray diffraction study. Phys. Rev. B, 88:104104, Sep 2013.

[31] Corentin Stéphan. *Schémas hydrodynamiques d'ordre élevé et cinétique de changement de phase*. PhD thesis, Université de Montpellier, 2025.

[32] Philip A. Thompson. A fundamental derivative in gasdynamics. *Physics of Fluids*, 14(9):1843–1849, September 1971.

[33] Y. S. Touloukian, R. K. Kirby, R. E. Taylor, and P. D. Desai. *Thermal expansion: metallic elements and alloys*. Springer, 1975.

[34] John M. Walsh, Melvin H. Rice, Robert G. McQueen, and Frederick L. Yarger. Shock-wave compressions of twenty-seven metals. equations of state of metals. Phys. Rev., 108:196–216, Oct 1957.

[35] Lianwen Wang, Qiang Wang, Aiping Xian, and Kunquan Lu. Precise measurement of the densities of liquid bi, sn, pb and sb. *Journal of Physics: Condensed Matter*, 15(6):777, feb 2003.

[36] Ya.B. Zel'dovitch and Yu.P. Raizer. *Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena*. Academic Press, 1967.