# A model and numerical method for compressible flows with capillary effects

Kevin Schmidmayer, Fabien Petitpas, Eric Daniel, Nicolas Favrie, Sergey

Gavrilyuk

► To cite this version:

Kevin Schmidmayer, Fabien Petitpas, Eric Daniel, Nicolas Favrie, Sergey Gavrilyuk. A model and numerical method for compressible flows with capillary effects. Journal of Computational Physics, 2017, 334, pp.468-496. 10.1016/j.jcp.2017.01.001. hal-01297432v3

HAL Id: hal-01297432

https://hal.science/hal-01297432v3

Submitted on 26 Jan 2017

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.

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

Distributed under a Creative Commons Attribution 4.0 International License
A model and numerical method for compressible flows with
capillary effects

Kevin Schmidmayer$^{a,*}$, Fabien Petitpas$^a$, Eric Daniel$^a$, Nicolas Favrie$^a$, Sergey Gavrilyuk$^a$

$^a$Aix Marseille Univ, CNRS, IUSTI, Marseille, France

**Abstract**

A new model for interface problems with capillary effects in compressible fluids is presented
together with a specific numerical method to treat capillary flows and pressure waves propa-
gation. This new multiphase model is in agreement with physical principles of conservation
and respects the second law of thermodynamics. A new numerical method is also proposed
where the global system of equations is split into several submodels. Each submodel is hy-
perbolic or weakly hyperbolic and can be solved with an adequate numerical method. This
method is tested and validated thanks to comparisons with analytical solutions (Laplace
law) and with experimental results on droplet breakup induced by a shock wave.

*Keywords:* diffuse interface, Godunov type methods, hyperbolic systems, multiphase flows, shock waves, surface tension

# 1. Introduction

The breakup of liquid droplets induced by high speed flows has a wide range of engi-
neering and scientific applications and has given rise to a large number of publications. In
some cases, this phenomenon causes damages as for example when droplets are impacting
aircrafts in supersonic flight causing erosion of its surface (Engel [8], Joseph et al. [23], Igra
and Takayama [20, 21]). Studying of droplets behavior in a high speed flow may also be
encountered when security issues are considered as, for example, for shock wave attenuation
(Chauvin et al. [4, 5]). Other applications can be found in explosive science or in combus-
tion systems where a liquid jet is atomized (Welch and Boyle [52], Meng and Colonius [33],
Devassy et al. [7]). Detailed reviews on droplet breakup can be found in Pilch and Erdman
[43], Wierzba and Takayama [53], Hsiang and Faeth [18].

Concerning numerical simulations, the breakup study is usually focused on the first stages
of droplet deformation when Richtmyer-Meshkov and/or Rayleigh-Taylor instabilities appear
(Yang *et al.* [55], Quirk and Karni [44], Layes and Le Metayer [28], Meng and Colonius [33]),
but not on the further stages when capillary and/or viscous effects become significant.

<footnote>*Corresponding author</footnote>

<footnote>Email addresses: kevin.schmidmayer@univ-amu.fr (Kevin Schmidmayer),
fabien.petitpas@univ-amu.fr (Fabien Petitpas), eric.daniel@univ-amu.fr (Eric Daniel),
nicolas.favrie@univ-amu.fr (Nicolas Favrie), sergey.gavrilyuk@univ-amu.fr (Sergey Gavrilyuk)</footnote>
In the last decades, several theoretical studies have been performed to treat capillary effects in multiphase flows. The seminal work of Brackbill et al. [3] succeeded in transforming a surface force into a volume force, quite easy to treat as a source term in a multiphase flow model. The surface tension volume force is expressed thanks to a color function $\tilde{c}(x)$. This approach has been used in Chen and Doolen [6], Sussman et al. [50], Gueyffier et al. [16], Osher and Fedkiw [36, 37], Tryggvason et al. [51], Périgaud and Saurel [39], Le Martelot et al. [29] where capillary effects are added into the momentum and the energy equations.

The aim of this work is to develop a mathematical model for fluid flows with capillary effects that is hyperbolic, verifies conservation principles and entropy inequality together with a suitable numerical method capable to treat the effect of the flow on the droplet from the short time scale when the shock wave interacts with the droplet to the long time scale when capillary effects become significant. We focus in this study on multiphase compressible fluid flows only. Viscous and heat conduction are not taken into account and will be a part of future works. Some ideas on the treatment of heat conduction in multiphase compressible flows can already be found in [40].

Section 2 presents the Brackbill et al. [3] method to treat the surface tension and a review of existing models with a conservative form of the capillary terms. In Section 3, the new model with capillary effects is presented. The model is in agreement with the conservation principles and with the second law of thermodynamics. It is shown that the model is weakly hyperbolic. It has two sound characteristics associated with the classical compression waves and two new sound characteristics associated with the capillary effects. However, for multiple contact characteristics one eigenvector is always missing. Section 4 is devoted to the building of a numerical method able to solve capillary terms in a conservative manner. The method is based on split models that are separately hyperbolic or weakly hyperbolic. These submodels are solved thanks to adequate numerical schemes. Section 5 presents the validation of the method on 2D test cases. It shows that the model and the numerical method are able to treat accurately both capillary effects and shock wave propagation. Quantitative comparisons are done with other methods based on source terms integration to show the importance of the conservative formulation. To illustrate the capabilities of the model, the aerodynamic breakup of a water column induced by a shock wave is numerically solved and is compared with experiments. In Appendix, the model derivation is given.

## 2. Compressible two-phase capillary flows: state of the art

### 2.1. Surface tension force and color function

The main difficulty in modeling the capillary effects is about considering a surface force in numerical models that solve volume average quantities. The seminal work of Brackbill et al. [3], called CSF (Continuum Surface Force) method, succeeded to do it by using a color function, $\tilde{c}(x)$. Thanks to this function, the surface tension volume force is then expressed:

$$F_v(x) = \sigma\kappa(x) \frac{\nabla\tilde{c}(x)}{[\tilde{c}]},$$
where $\sigma$ is the surface tension coefficient and $\kappa(\mathbf{x})$ the local curvature of the interface defined by:

$$ \kappa(\mathbf{x}) = -\nabla \cdot \mathbf{n}(\mathbf{x}), $$

where $\mathbf{n}(\mathbf{x})$ is the normal vector to the interface between the both phases:

$$ \mathbf{n}(\mathbf{x}) = \frac{\nabla \tilde{c}(\mathbf{x})}{\|\nabla \tilde{c}(\mathbf{x})\|} $$

The color function $\tilde{c}(\mathbf{x})$ allows locations of the different fluids and the interface. $\tilde{c}(\mathbf{x})$ is defined as:

$$ \tilde{c}(\mathbf{x}) = \begin{cases} c_1 & \text{in fluid 1,} \\ c_2 & \text{in fluid 2,} \\ c_1 \le \tilde{c}(\mathbf{x}) \le c_2 & \text{in the transition region.} \end{cases} \qquad (1) $$

In the transition region $\tilde{c}(\mathbf{x})$ is given by interpolation, meaning that the interface has a non zero thickness. $[\tilde{c}] = c_2 - c_1$ is the jump of the color function.

It is assumed that the color function obeys a transport equation [3]:

$$ \frac{\partial \tilde{c}(\mathbf{x})}{\partial t} + \mathbf{u}_I \cdot \nabla \tilde{c}(\mathbf{x}) = 0, $$

where $\mathbf{u}_I$ is the interface velocity.

Numerical results using this force can be found in [6, 29, 36, 37, 50, 51]. In these references, the surface tension force is treated as source terms in the momentum and the energy equations. Nevertheless, this treatment of capillary effects violates conservation principles.

## 2.2. Review of existing compressible models with capillary effects

Two family of methods are available to treat interface problems.

* The first family of methods considers interfaces as sharp. Sharp interfaces can be obtained using interface-tracking methods, where usually a level set function tracks the interface (Osher and Sethian [38]). However, such formulations often involve slight modifications of the governing equations. For example a pressure evolution equation can replace the energy equation around interfaces (Karni [25, 26]). In the Ghost Fluid Method (Fedkiw et al. [12]) and its simplified version (Koren et al. [27]), thermodynamically similar variables are added across interfaces to complete stencils. These methods do not generate spurious oscillations at interfaces, however, they are basically not conservative (Liu et al. [30]). Even if progress has been recently done concerning this aspect (Hu et al. [19], Luo et al. [31], Han et al. [17], Schranner et al. [49]), that makes them less desirable for problems where shock waves are involved.

* The second family of methods is described below and is called diffuse interface methods. In this class of methods, interfaces are not explicitly tracked but allowed to diffuse numerically (Abgrall and Karni [1], Saurel and Abgrall [45]). These methods are particularly interesting because they are able to deal with dynamic appearance and
disappearance of interfaces. Moreover, this is also the only class of models where the thermodynamics of mixture cells is well defined, thanks to a specific equation of state for each phase (liquid or gas).

The study of capillary effects within the framework of the diffuse interface methods is based on the generalization of the Allaire et al. [2] model. This model originally does not include capillary effects. The second law of thermodynamics is verified only if the condition of thermal equilibrium $T = T_k$, with $k = \{1, 2\}$, is retained. Périgaud and Saurel [39] extended this model by including the capillary effects. As a result, the surface tension volume force appears as a flux term in the momentum equation as well as the work of this force $F_v(x) \cdot u$ in the total energy equation. In this reference, as Gueyffier et al. [16] did within the incompressible flows framework, a conservative formulation was obtained:

$$ 
\left\{ 
\begin{aligned}
\frac{\partial \alpha_1}{\partial t} + \mathbf{u} \cdot \nabla \alpha_1 &= 0, \\
\frac{\partial \alpha_k \rho_k}{\partial t} + \nabla \cdot (\alpha_k \rho_k \mathbf{u}) &= 0, \\
\frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot \left( \rho \mathbf{u} \otimes \mathbf{u} + P \bar{I} - \sigma \left( \| \nabla \alpha_1 \| \bar{I} - \frac{\nabla \alpha_1 \otimes \nabla \alpha_1}{\| \nabla \alpha_1 \|} \right) \right) &= 0, \\
\frac{\partial \rho E + \varepsilon_\sigma}{\partial t} + \nabla \cdot \left( \mathbf{u} (\rho E + \varepsilon_\sigma + P) - \sigma \left( \| \nabla \alpha_1 \| \bar{I} - \frac{\nabla \alpha_1 \otimes \nabla \alpha_1}{\| \nabla \alpha_1 \|} \right) \right) \cdot \mathbf{u} &= 0,
\end{aligned}
\right. 
$$

where $\alpha_k$ and $\rho_k$ are the volume fraction and the density of phase k. $\rho$, $\mathbf{u}$, $P$, $E = e + \frac{1}{2} \|\mathbf{u}\|^2$ and $e$ are respectively the mixture variables for density, velocity, pressure, total energy and internal energy. The specific mixture internal energy is defined as $e = \sum_k Y_k e_k (\rho_k, P)$ and each fluid is governed by its own equation of state (EOS) $e_k = e_k (\rho_k, P)$. Introducing the capillary effects in fluxes (conservative formulation) leads to a new term in the mixture total energy equation, $\varepsilon_\sigma = \sigma \|\nabla \alpha_1\|$ which is defined as a capillary potential energy. Note that $\alpha_1$ obeys a transport equation and plays the role of the color function (1) because $\mathbf{u}_I = \mathbf{u}$ when the velocity equilibrium is considered. In [39], the authors used a pressure equilibrium closure, allowing the resolution of interface problems dedicated to high speed flows. Nevertheless, this closure does not allow to respect the second law of thermodynamics for the mixture.

Le Martelot et al. [29] corrected this drawback by including the capillary effects in a velocity, pressure and temperature equilibrium system of equations by considering the mass fraction as the color function: the second law of thermodynamics is then respected. However, the characteristic velocities were estimated under assumption that the local curvature is fixed. Even if such an approach is physically reasonable, it does not allow us to rigorously estimate the propagation speeds. Also, the thermal equilibrium closure is reasonable for the description of the boiling phenomenon, but it is too drastic for droplet breakup study under high speed flows.

It becomes then obvious that the introduction of a new model having clearly defined sound speeds with the capillary effects in conservative form that can still be compatible with the mixture entropy increase without the temperature equilibrium assumption is necessary.
### 3. Mechanical equilibrium model with capillary effects

#### 3.1. The model

The new model we propose is reminiscent that of the Kapila et al. model [24]. This last model has been shown to be suitable to treat interface problems between compressible fluids. Multiple extensions of this model have been developed to solve several concrete problems dealing with phase transition (Massoni et al. [32], Saurel et al. [47], cavitation (Petitpas et al. [41]), detonation in high energetic materials (Petitpas et al. [42]), solid-fluid interaction and compaction of granular media (Favrie and Gavrilyuk [10, 9]) and low Mach number flows (Murrone and Guillard [34]). It considers compressible two-phase flows in mechanical equilibrium (pressure and velocity equilibrium). In this model, the volume fraction obeys the following equation coming from the pressure equilibrium condition:

$$ \frac{\partial \alpha_1}{\partial t} + \mathbf{u} \cdot \nabla \alpha_1 = K \nabla \cdot \mathbf{u}, \qquad (2) $$

where the term $K \nabla \cdot \mathbf{u}$ accounts for the differences in the acoustic behaviour of both phases. K is given by:

$$ K = \frac{\rho_2 a_2^2 - \rho_1 a_1^2}{\frac{\rho_2 a_2^2}{\alpha_2} + \frac{\rho_1 a_1^2}{\alpha_1}}, $$

$a_k$ being the speed of sound of phase k.

The capillary effects are then added in conservative form. Because the color function is a purely geometric variable, a supplementary equation for the color function is added. It will be shown in the following that the introduction of this supplementary equation is an important ingredient to perform the hyperbolicity study. Applying the Hamilton principle for the model derivation (see Appendix for details), the model becomes:

$$ \left\{ \begin{aligned} \frac{\partial \alpha_1}{\partial t} + \mathbf{u} \cdot \nabla \alpha_1 - K \nabla \cdot \mathbf{u} &= 0, \\ \frac{\partial \alpha_k \rho_k}{\partial t} + \nabla \cdot (\alpha_k \rho_k \mathbf{u}) &= 0, \\ \frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u} + P \bar{I} + \bar{\Omega}) &= 0, \\ \frac{\partial \rho E + \varepsilon_\sigma}{\partial t} + \nabla \cdot (\rho E + \varepsilon_\sigma + P) \mathbf{u} + \bar{\Omega} \cdot \mathbf{u} &= 0, \\ \frac{\partial c}{\partial t} + \mathbf{u} \cdot \nabla c &= 0, \end{aligned} \right. \qquad (3) $$

with $\bar{\Omega}$ being the capillary tensor given by:

$$ \bar{\Omega} = -\sigma \left( \| \nabla c \| \bar{I} - \frac{\nabla c \otimes \nabla c}{\| \nabla c \|} \right). \qquad (4) $$

The capillary energy is equal to $\varepsilon_\sigma = \sigma \| \nabla c \|$ and the color function is normalized by its jump: $c = \tilde{c}/[\tilde{c}]$.

The surface tension terms do not affect the pressure and the entropy equations. The pressure equation evolution then reads:

$$ \frac{dP}{dt} + \rho a^2 \nabla \cdot (\mathbf{u}) = 0. $$
where a is the Wood mixture speed of sound [54]:

$$a^2 = \left( \rho \sum_k \frac{\alpha_k}{\rho_k a_k^2} \right)^{-1} . \qquad (5)$$

As in Kapila et al. [24], the entropy equations remain unchanged in continuous motion:

$$\frac{ds_k}{dt} = 0,$$

with the material derivative operator $d(\cdot)/dt = \partial(\cdot)/\partial t + \mathbf{u} \cdot \nabla(\cdot)$. The mixture entropy is equation obviously assured:

$$\frac{ds}{dt} = \sum_k \frac{dY_k s_k}{dt} = 0.$$

## 3.2. Hyperbolicity

In this section, the study of the hyperbolicity of system (3) is done. Thanks to the rotational invariance of the equations, the study can be reduced to the study of only 1D equation.

### 3.2.1. Primitive form

The model needs to be transformed into a vector form:

$$\frac{\partial \mathbf{W}}{\partial t} + \bar{A}(\mathbf{W}) \frac{\partial \mathbf{W}}{\partial x} = 0. \qquad (6)$$

The color function is first rewritten by taking its gradient to obtain a conservative equation for $\mathbf{w} = \nabla c$:

$$\frac{\partial \mathbf{w}}{\partial t} + \nabla (\mathbf{u} \cdot \mathbf{w}) = 0.$$

Or, in an equivalent form:

$$\frac{\partial \mathbf{w}}{\partial t} + \left( \frac{\partial \mathbf{w}}{\partial x} \right)^T \cdot \mathbf{u} + \left( \frac{\partial \mathbf{u}}{\partial x} \right)^T \cdot \mathbf{w} = 0. \qquad (7)$$

Since $\mathbf{w}$ is a gradient, we have:

$$\mathrm{curl}(\mathbf{w}) = \mathbf{0}, \qquad (8)$$

i.e.

$$\left( \frac{\partial \mathbf{w}}{\partial x} \right)^T = \left( \frac{\partial \mathbf{w}}{\partial x} \right).$$

Thus, equation (7) with the constraint (8) becomes:

$$\frac{\partial \mathbf{w}}{\partial t} + \left( \frac{\partial \mathbf{w}}{\partial x} \right) \cdot \mathbf{u} + \left( \frac{\partial \mathbf{u}}{\partial x} \right)^T \cdot \mathbf{w} = \mathbf{0}. \qquad (9)$$
Constraint (8) also appears in solid mechanics where the formulation of the extended system under a “curl” constraint is important in the hyperbolicity study (see Ndanou et al. [35]). Equation (9) may be developed in three dimensions, where we denote $\mathbf{u} = (u, v, w)^T$, $\mathbf{w} = (w_1, w_2, w_3)^T$ and $\mathbf{x} = (x, y, z)^T$:

$$ \left\{ \begin{aligned} \frac{\partial w_1}{\partial t} + \frac{\partial w_1}{\partial x} u + \frac{\partial w_1}{\partial y} v + \frac{\partial w_1}{\partial z} w + \frac{\partial u}{\partial x} w_1 + \frac{\partial v}{\partial x} w_2 + \frac{\partial w}{\partial x} w_3 &= 0, \\ \frac{\partial w_2}{\partial t} + \frac{\partial w_2}{\partial x} u + \frac{\partial w_2}{\partial y} v + \frac{\partial w_2}{\partial z} w + \frac{\partial u}{\partial y} w_1 + \frac{\partial v}{\partial y} w_2 + \frac{\partial w}{\partial y} w_3 &= 0, \\ \frac{\partial w_3}{\partial t} + \frac{\partial w_3}{\partial x} u + \frac{\partial w_3}{\partial y} v + \frac{\partial w_3}{\partial z} w + \frac{\partial u}{\partial z} w_1 + \frac{\partial v}{\partial z} w_2 + \frac{\partial w}{\partial z} w_3 &= 0. \end{aligned} \right. \quad (10) $$

In the case where all the variables depend only on $(t, x)$, (10) reduces to:

$$ \left\{ \begin{aligned} \frac{\partial w_1}{\partial t} + \frac{\partial w_1}{\partial x} u + \frac{\partial u}{\partial x} w_1 + \frac{\partial v}{\partial x} w_2 + \frac{\partial w}{\partial x} w_3 &= 0, \\ \frac{\partial w_2}{\partial t} + \frac{\partial w_2}{\partial x} u &= 0, \\ \frac{\partial w_3}{\partial t} + \frac{\partial w_3}{\partial x} u &= 0. \end{aligned} \right. $$

The vector $\mathbf{w}$ is introduced to rewrite components for the capillary tensor $\overline{\Omega}$ (4) in the $x$-direction:

$$ \Omega_{11} = \sigma \left( \frac{-w_2^2 - w_3^2}{\sqrt{w_1^2 + w_2^2 + w_3^2}} \right), $$

$$ \Omega_{12} = \sigma \left( \frac{w_1 w_2}{\sqrt{w_1^2 + w_2^2 + w_3^2}} \right), $$

$$ \Omega_{13} = \sigma \left( \frac{w_1 w_3}{\sqrt{w_1^2 + w_2^2 + w_3^2}} \right). $$

Finally, the system rewritten in vector form reads for two phases:

$$ \left\{ \begin{aligned} \frac{\partial \alpha_1}{\partial t} + u \frac{\partial \alpha_1}{\partial x} - K \frac{\partial u}{\partial x} &= 0, \\ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + \frac{1}{\rho} \left( \frac{\partial P}{\partial x} + \frac{\partial \Omega_{11}}{\partial w_1} \frac{\partial w_1}{\partial x} + \frac{\partial \Omega_{11}}{\partial w_2} \frac{\partial w_2}{\partial x} + \frac{\partial \Omega_{11}}{\partial w_3} \frac{\partial w_3}{\partial x} \right) &= 0, \\ \frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + \frac{1}{\rho} \left( \frac{\partial \Omega_{12}}{\partial w_1} \frac{\partial w_1}{\partial x} + \frac{\partial \Omega_{12}}{\partial w_2} \frac{\partial w_2}{\partial x} + \frac{\partial \Omega_{12}}{\partial w_3} \frac{\partial w_3}{\partial x} \right) &= 0, \\ \frac{\partial w}{\partial t} + u \frac{\partial w}{\partial x} + \frac{1}{\rho} \left( \frac{\partial \Omega_{13}}{\partial w_1} \frac{\partial w_1}{\partial x} + \frac{\partial \Omega_{13}}{\partial w_2} \frac{\partial w_2}{\partial x} + \frac{\partial \Omega_{13}}{\partial w_3} \frac{\partial w_3}{\partial x} \right) &= 0, \\ \frac{\partial P}{\partial t} + u \frac{\partial P}{\partial x} + \rho a^2 \frac{\partial u}{\partial x} &= 0, \\ \frac{\partial w_1}{\partial t} + u \frac{\partial w_1}{\partial x} + w_1 \frac{\partial u}{\partial x} + w_2 \frac{\partial v}{\partial x} + w_3 \frac{\partial w}{\partial x} &= 0, \\ \frac{\partial w_2}{\partial t} + u \frac{\partial w_2}{\partial x} &= 0, \\ \frac{\partial w_3}{\partial t} + u \frac{\partial w_3}{\partial x} &= 0, \\ \frac{\partial s_1}{\partial t} + u \frac{\partial s_1}{\partial x} &= 0, \\ \frac{\partial s_2}{\partial t} + u \frac{\partial s_2}{\partial x} &= 0, \\ \frac{\partial Y_1}{\partial t} + u \frac{\partial Y_1}{\partial x} &= 0, \\ \frac{\partial c}{\partial t} + u \frac{\partial c}{\partial x} &= 0. \end{aligned} \right. $$
### 3.2.2. Eigenvalues

The vector $\mathbf{W}$ and the matrix $\overline{\mathbf{A}}$ in (6) are then defined by:

$$ \mathbf{W} = [\alpha_1, u, v, w, P, w_1, w_2, w_3, s_1, s_2, Y_1, c]^T, $$

$$ \overline{\mathbf{A}}(\mathbf{W}) = \begin{bmatrix} u & -K & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & u & 0 & 0 & \frac{1}{\rho} & \frac{1}{\rho} \frac{\partial \Omega_{11}}{\partial w_1} & \frac{1}{\rho} \frac{\partial \Omega_{11}}{\partial w_2} & \frac{1}{\rho} \frac{\partial \Omega_{11}}{\partial w_3} & 0 & 0 & 0 \\ 0 & 0 & u & 0 & 0 & \frac{1}{\rho} \frac{\partial \Omega_{12}}{\partial w_1} & \frac{1}{\rho} \frac{\partial \Omega_{12}}{\partial w_2} & \frac{1}{\rho} \frac{\partial \Omega_{12}}{\partial w_3} & 0 & 0 & 0 \\ 0 & 0 & 0 & u & 0 & \frac{1}{\rho} \frac{\partial \Omega_{13}}{\partial w_1} & \frac{1}{\rho} \frac{\partial \Omega_{13}}{\partial w_2} & \frac{1}{\rho} \frac{\partial \Omega_{13}}{\partial w_3} & 0 & 0 & 0 \\ 0 & \rho a^2 & 0 & 0 & u & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & w_1 & w_2 & w_3 & 0 & u & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & u & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & u & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & u & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & u & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & u \end{bmatrix} $$

The determinant of matrix $\overline{\mathbf{A}}$ can be calculated:

$$ \det(\overline{\mathbf{A}}(\mathbf{W}) - \lambda \overline{\mathbf{I}}) = (u - \lambda)^8 [(u - \lambda)^4 + L(u - \lambda)^2 + M], $$

where:

$$ L = -a^2 - \frac{1}{\rho} \left( w_1 \frac{\partial \Omega_{11}}{\partial w_1} + w_2 \frac{\partial \Omega_{12}}{\partial w_1} + w_3 \frac{\partial \Omega_{13}}{\partial w_1} \right), $$

$$ M = \frac{a^2}{\rho} \left( w_2 \frac{\partial \Omega_{12}}{\partial w_1} + w_3 \frac{\partial \Omega_{13}}{\partial w_1} \right). $$

8 real eigenvalues are straightforward found:

$$ \lambda_{1,2,3,4,5,6,7,8} = u. $$

The 4 other eigenvalues are the roots of the quadratic equation:

$$ X^2 + LX + M = 0, \tag{11} $$

with $X = (u - \lambda)^2$. If the discriminant $\Delta$ is positive, the roots of equation (11) are real and complex numbers otherwise. Defining the components of the normal vector $n_k = w_k / \|\mathbf{w}\|$, after some calculations the discriminant is:

$$ \Delta = \left( a^2 + \frac{\sigma \|\mathbf{w}\|}{\rho} (n_2^2 + n_3^2) \right)^2 - 4a^2 \frac{\sigma \|\mathbf{w}\|}{\rho} (n_2^2 + n_3^2)^2. $$

But:

$$ (n_2^2 + n_3^2) \le 1. $$
<figure><img src="image_2.png" /><figcaption>Figure 1: Model (3) admits two types of acoustic waves.</figcaption></figure>

Hence:

$$ \begin{aligned} \Delta &\ge \left( a^2 + \frac{\sigma \|\mathbf{w}\|}{\rho} (n_2^2 + n_3^2) \right)^2 - 4a^2 \frac{\sigma \|\mathbf{w}\|}{\rho} (n_2^2 + n_3^2) \\ &= \left( a^2 - \frac{\sigma \|\mathbf{w}\|}{\rho} (n_2^2 + n_3^2) \right)^2 \ge 0. \end{aligned} $$

The roots of the quadratic equation (11) are then all real. Thus, the corresponding eigenvalues of model (3) are real, explicit and given by:

$$ \lambda_{9,10} = u \pm a_s, \quad (12) $$

$$ \lambda_{11,12} = u \pm a_c, \quad (13) $$

where:

$$ a_s^2 = \frac{a^2 + b + \sqrt{(a^2 + b)^2 - 4a^2b(n_2^2 + n_3^2)}}{2}, $$

$$ a_c^2 = \frac{a^2 + b - \sqrt{(a^2 + b)^2 - 4a^2b(n_2^2 + n_3^2)}}{2}, $$

a is the Wood mixture speed of sound previously defined (5) and:

$$ b = \frac{\sigma \|\mathbf{w}\|}{\rho} (n_2^2 + n_3^2). $$

Finally, there are 8 multiple eigenvalues $\lambda = u$ and 4 eigenvalues corresponding to sound waves (12) and capillary waves (13) (see Figure 1).

The system is hyperbolic if the multiple eigenvalues $\lambda = u$ have exactly 8 linearly independent eigenvectors. One can prove that it is not the case, one eigenvector is missing (see Appendix), so the system is only weakly hyperbolic. The whole system of the eigenvectors will not be given because this equilibrium system will not be numerically solved in this form (see Section 4 for discussion of the numerical method).
# 4. Numerical resolution of model (3)

## 4.1. Basic ideas

The numerical resolution of model (3) represents a challenge regarding the two following points:

* The first difficulty is due to the $K \nabla \cdot \mathbf{u}$ term in the volume fraction evolution equation. This is the precious ingredient leading to the respect of the mixture entropy equation when a pressure equilibrium assumption is retained. Nevertheless, the presence of this non-conservative term considerably complicates the numerical method which crucially depends on the choice of appropriate Rankine-Hugoniot relations (see [14] and [48] for details). This is a reason why a non-equilibrium pressure model is preferred with a pressure relaxation term instead of the non-conservative term $K \nabla \cdot \mathbf{u}$ in the volume fraction equation. This model is presented below.

* The second difficulty is in the simultaneous treatment of 5 waves which are present in the model. This difficulty is circumvented by the use of split models that are proven to be hyperbolic and weakly hyperbolic respectively.

### 4.1.1. Pressure relaxation model with capillary effects

Because the non-conservative volume fraction equation presents a major problem regarding numerical resolution, the following weakly hyperbolic non-equilibrium pressure model is proposed:

$$ 
\left\{ 
\begin{aligned}
\frac{\partial \alpha_1}{\partial t} + \mathbf{u} \cdot \nabla \alpha_1 &= \mu (P_1 - P_2), \\
\frac{\partial \alpha_1 \rho_1}{\partial t} + \nabla \cdot (\alpha_1 \rho_1 \mathbf{u}) &= 0, \\
\frac{\partial \alpha_2 \rho_2}{\partial t} + \nabla \cdot (\alpha_2 \rho_2 \mathbf{u}) &= 0, \\
\frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot \left( \rho \mathbf{u} \otimes \mathbf{u} + P \bar{I} + \bar{\Omega} \right) &= 0, \\
\frac{\partial \alpha_1 \rho_1 e_1}{\partial t} + \nabla \cdot (\alpha_1 \rho_1 e_1 \mathbf{u}) + \alpha_1 P_1 \nabla \cdot \mathbf{u} &= -\mu P_I (P_1 - P_2), \\
\frac{\partial \alpha_2 \rho_2 e_2}{\partial t} + \nabla \cdot (\alpha_2 \rho_2 e_2 \mathbf{u}) + \alpha_2 P_2 \nabla \cdot \mathbf{u} &= \mu P_I (P_1 - P_2), \\
\frac{\partial c}{\partial t} + \mathbf{u} \cdot \nabla c &= 0,
\end{aligned}
\right. 
\quad (14) 
$$

where $\mu$ is the pressure relaxation coefficient, $P_I = \frac{Z_2 P_1 + Z_1 P_2}{Z_1 + Z_2}$ (see [46] for details) and $Z_k = \rho_k a_k$ is the acoustic impedance of the phase $k$. The mixture pressure is given by:

$$ P = \alpha_1 P_1 + \alpha_2 P_2. $$

Due to the condition $P_1 \neq P_2$ in this model, the total energy equation of the mixture is replaced by the internal energy equation for each phase. Nevertheless, the mixture total energy equation of the system can be written in usual form:

$$ \frac{\partial \rho E + \varepsilon_{\sigma}}{\partial t} + \nabla \cdot \left( (\rho E + \varepsilon_{\sigma} + P) \mathbf{u} + \bar{\Omega} \cdot \mathbf{u} \right) = 0. \quad (15) $$

The equation (15) is redundant when both phasic internal energy equations are solved, but it will appear to be an important ingredient for numerical method to ensure the energy conservation and to preserve a correct treatment of shock waves.
Following Section 3.2, one can prove that model (14) is also weakly hyperbolic.

One can note that the surface tension effects are missing in the phasic energy equations since it is only a mixture characteristic.

The entropy equations for system (14) can also be written:

$$ \alpha_1 \rho_1 T_1 \frac{ds_1}{dt} = \mu(P_1 - P_2)^2 \frac{Z_1}{Z_1 + Z_2}, $$

$$ \alpha_2 \rho_2 T_2 \frac{ds_2}{dt} = \mu(P_1 - P_2)^2 \frac{Z_2}{Z_1 + Z_2}, $$

that insures the mixture entropy ($s = Y_1s_1 + Y_2s_2$) increases.

Model (3) can be recovered as the asymptotic limit of Model (14) when the pressure of each phase tends to be equal. A special splitting procedure will be done for the numerical resolution of model (14).

### 4.1.2. Splitting procedure

Model (14) without the relaxation terms is split in two submodels. The first submodel does not take into account the surface tension terms. The hyperbolicity is then easily verified. The second submodel contains the only capillary terms and is proven weakly hyperbolic. Such an approach was first proposed in Favrie et al. [11] where the splitting procedure was used to separate the treatment of longitudinal and shear waves in hyperelasticity. The submodels are presented below only in the x-direction.
## Hyperbolic submodel 1

The first submodel is similar to that presented in [48] with additional decoupled equations for the gradient of the color function:

$$ 
\left\{
\begin{aligned}
& \frac{\partial \alpha_1}{\partial t} + u \frac{\partial \alpha_1}{\partial x} &&= 0, \\
& \frac{\partial \alpha_1 \rho_1}{\partial t} + \frac{\partial \alpha_1 \rho_1 u}{\partial x} &&= 0, \\
& \frac{\partial \alpha_2 \rho_2}{\partial t} + \frac{\partial \alpha_2 \rho_2 u}{\partial x} &&= 0, \\
& \frac{\partial \rho u}{\partial t} + \frac{\partial \rho u^2 + \alpha_1 P_1 + \alpha_2 P_2}{\partial x} &&= 0, \\
& \frac{\partial \rho v}{\partial t} + \frac{\partial \rho u v}{\partial x} &&= 0, \\
& \frac{\partial \rho w}{\partial t} + \frac{\partial \rho u w}{\partial x} &&= 0, \\
& \frac{\partial \alpha_1 \rho_1 e_1}{\partial t} + \frac{\partial \alpha_1 \rho_1 e_1 u}{\partial x} + \alpha_1 P_1 \frac{\partial u}{\partial x} &&= 0, \\
& \frac{\partial \alpha_2 \rho_2 e_2}{\partial t} + \frac{\partial \alpha_2 \rho_2 e_2 u}{\partial x} + \alpha_2 P_2 \frac{\partial u}{\partial x} &&= 0, \\
& \frac{\partial w_1}{\partial t} + \frac{\partial w_1 u}{\partial x} &&= 0, \\
& \frac{\partial w_2}{\partial t} + u \frac{\partial w_2}{\partial x} &&= 0, \\
& \frac{\partial w_3}{\partial t} + u \frac{\partial w_3}{\partial x} &&= 0.
\end{aligned}
\right.
\quad (16) $$

This system describes only the transport and the compression waves. The equation for $w_1$ is taken in conservative form to let the possibility to consider weak solutions. The other terms in this equation will be treated in the second submodel.

Model (16) may be rewritten in a vector form (6) with:

$$ \mathbf{W} = [\alpha_1, s_1, s_2, u, v, w, P_1, P_2, w_1, w_2, w_3]^T $$

$$ \bar{\bar{A}}(\mathbf{W}) = \begin{bmatrix}
u & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & u & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & u & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
\frac{P_1-P_2}{\rho} & 0 & 0 & u & 0 & 0 & \frac{\alpha_1}{\rho} & \frac{\alpha_2}{\rho} & 0 & 0 \\
0 & 0 & 0 & 0 & u & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & u & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & \rho_1 a_1^2 & 0 & 0 & u & 0 & 0 & 0 \\
0 & 0 & 0 & \rho_2 a_2^2 & 0 & 0 & 0 & u & 0 & 0 \\
0 & 0 & 0 & w_1 & 0 & 0 & 0 & 0 & u & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & u \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & u
\end{bmatrix} $$

The eigenvalues of the system are:

$$ \lambda_{1,2,3,4,5,6,7,8,9} = u, $$
$$ \lambda_{10} = u - a_f, $$

$$ \lambda_{11} = u + a_f, $$

where $a_f$ is the frozen mixture sound speed:

$$ a_f^2 = Y_1 a_1^2 + Y_2 a_2^2. $$

The hyperbolicity of this first submodel is proven in [48].

### Weakly hyperbolic submodel 2

The second submodel is:

$$ \left\{ \begin{aligned} \frac{\partial \alpha_1}{\partial t} &= 0, \\ \frac{\partial \alpha_1 \rho_1}{\partial t} &= 0, \\ \frac{\partial \alpha_2 \rho_2}{\partial t} &= 0, \\ \frac{\partial \rho u}{\partial t} + \left( \frac{\partial \Omega_{11}}{\partial w_1} \frac{\partial w_1}{\partial x} + \frac{\partial \Omega_{11}}{\partial w_2} \frac{\partial w_2}{\partial x} + \frac{\partial \Omega_{11}}{\partial w_3} \frac{\partial w_3}{\partial x} \right) &= 0, \\ \frac{\partial \rho v}{\partial t} + \left( \frac{\partial \Omega_{12}}{\partial w_1} \frac{\partial w_1}{\partial x} + \frac{\partial \Omega_{12}}{\partial w_2} \frac{\partial w_2}{\partial x} + \frac{\partial \Omega_{12}}{\partial w_3} \frac{\partial w_3}{\partial x} \right) &= 0, \\ \frac{\partial \rho w}{\partial t} + \left( \frac{\partial \Omega_{13}}{\partial w_1} \frac{\partial w_1}{\partial x} + \frac{\partial \Omega_{13}}{\partial w_2} \frac{\partial w_2}{\partial x} + \frac{\partial \Omega_{13}}{\partial w_3} \frac{\partial w_3}{\partial x} \right) &= 0, \\ \frac{\partial \alpha_1 \rho_1 e_1}{\partial t} &= 0, \\ \frac{\partial \alpha_2 \rho_2 e_2}{\partial t} &= 0, \\ \frac{\partial w_1}{\partial t} + w_2 \frac{\partial v}{\partial x} + w_3 \frac{\partial w}{\partial x} &= 0, \\ \frac{\partial w_2}{\partial t} &= 0, \\ \frac{\partial w_3}{\partial t} &= 0. \end{aligned} \right. \quad (17) $$

This second system describes the capillary effects. Also, the non-conservative product in the equation for $w_1$ is well defined because $w_2$ and $w_3$ are continuous through the shock.

Model (17) may be rewritten in a vector form (6) with:

$$ \mathbf{W} = [\alpha_1, s_1, s_2, u, v, w, P_1, P_2, w_1, w_2, w_3]^T $$

$$ \bar{\bar{A}}(\mathbf{W}) = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{\rho} \frac{\partial \Omega_{11}}{\partial w_1} & \frac{1}{\rho} \frac{\partial \Omega_{11}}{\partial w_2} & \frac{1}{\rho} \frac{\partial \Omega_{11}}{\partial w_3} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{\rho} \frac{\partial \Omega_{12}}{\partial w_1} & \frac{1}{\rho} \frac{\partial \Omega_{12}}{\partial w_2} & \frac{1}{\rho} \frac{\partial \Omega_{12}}{\partial w_3} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{\rho} \frac{\partial \Omega_{13}}{\partial w_1} & \frac{1}{\rho} \frac{\partial \Omega_{13}}{\partial w_2} & \frac{1}{\rho} \frac{\partial \Omega_{13}}{\partial w_3} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & w_2 & w_3 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix} $$

The eigenvalues of the system are:

$$ \lambda_{1,2,3,4,5,6,7,8,9} = 0, $$
$$ \lambda_{10} = -(n_2^2 + n_3^2) \sqrt{\frac{\sigma \|\mathbf{w}\|}{\rho}}, $$

$$ \lambda_{11} = (n_2^2 + n_3^2) \sqrt{\frac{\sigma \|\mathbf{w}\|}{\rho}}. $$

The eigenvalues of the system are all real but, as in the case of model (3), there are not as much linearly independent eigenvectors as eigenvalues: one eigenvector is always missing. So, the second submodel is weakly hyperbolic.

## 4.2. Numerical Method

Finally, the numerical method is presented as a 3-step method. Each step is successively performed in order to circumvent specific numerical problems:

* First, the hyperbolic non-equilibrium pressure model (16) is solved using a Godunov-type method.

* Second, model (17) is solved. A specific attention is paid to the choice for the flux terms in order to ensure the momentum and energy conservation.

* Third, a relaxation procedure leads to the pressure equilibrium.

The chain of these three steps is equivalent to solve model (3). Each step of the method is presented in details hereafter.

The full system of equations is first rewritten in the following vector form:

$$ \frac{\partial \mathbf{U}}{\partial t} + \frac{\partial (\mathbf{F}_h^x(\mathbf{U}) + \mathbf{F}_c^x(\mathbf{U}))}{\partial x} + \frac{\partial (\mathbf{F}_h^y(\mathbf{U}) + \mathbf{F}_c^y(\mathbf{U}))}{\partial y} + \frac{\partial (\mathbf{F}_h^z(\mathbf{U}) + \mathbf{F}_c^z(\mathbf{U}))}{\partial z} + \mathbf{H}_{nc} \nabla \cdot \mathbf{u} = \mathbf{H}_{relax}. $$

The vector $\mathbf{U}$ contains the unknown quantities defined in the system:

$$ \mathbf{U} = [\alpha_1, \alpha_1 \rho_1, \alpha_2 \rho_2, \rho u, \rho v, \rho w, \alpha_1 \rho_1 e_1, \alpha_2 \rho_2 e_2, c, \rho E + \varepsilon_\sigma]^T $$

The vectors $\mathbf{F}_h^\beta(\mathbf{U})$, $\mathbf{F}_c^\beta(\mathbf{U})$, $\mathbf{H}_{nc}$ and $\mathbf{H}_{relax}$, with $\beta = \{x, y, z\}$, contain respectively the hydrodynamic fluxes, the capillary fluxes, the non-conservative terms and the relaxation terms:
$$ \mathbf{F}_h^x (\mathbf{U}) = \begin{bmatrix} \alpha_1 u \\ \alpha_1 \rho_1 u \\ \alpha_2 \rho_2 u \\ \rho u^2 + P \\ \rho uv \\ \rho uw \\ \alpha_1 \rho_1 e_1 u \\ \alpha_2 \rho_2 e_2 u \\ cu \\ (\rho E + P) u \end{bmatrix} $$

$$ \mathbf{F}_h^y (\mathbf{U}) = \begin{bmatrix} \alpha_1 v \\ \alpha_1 \rho_1 v \\ \alpha_2 \rho_2 v \\ \rho uv \\ \rho v^2 + P \\ \rho vw \\ \alpha_1 \rho_1 e_1 v \\ \alpha_2 \rho_2 e_2 v \\ cv \\ (\rho E + P) v \end{bmatrix} $$

$$ \mathbf{F}_h^z (\mathbf{U}) = \begin{bmatrix} \alpha_1 w \\ \alpha_1 \rho_1 w \\ \alpha_2 \rho_2 w \\ \rho uw \\ \rho vw \\ \rho w^2 + P \\ \alpha_1 \rho_1 e_1 w \\ \alpha_2 \rho_2 e_2 w \\ cw \\ (\rho E + P) w \end{bmatrix} $$

$$ \mathbf{F}_c^x (\mathbf{U}) = \begin{bmatrix} 0 \\ 0 \\ 0 \\ \Omega_{11} \\ \Omega_{12} \\ \Omega_{13} \\ 0 \\ 0 \\ 0 \\ \varepsilon_\sigma u + \Omega_{11}u + \Omega_{12}v + \Omega_{13}w \end{bmatrix} $$

$$ \mathbf{F}_c^y (\mathbf{U}) = \begin{bmatrix} 0 \\ 0 \\ 0 \\ \Omega_{21} \\ \Omega_{22} \\ \Omega_{23} \\ 0 \\ 0 \\ 0 \\ \varepsilon_\sigma u + \Omega_{21}u + \Omega_{22}v + \Omega_{23}w \end{bmatrix} $$

$$ \mathbf{F}_c^z (\mathbf{U}) = \begin{bmatrix} 0 \\ 0 \\ 0 \\ \Omega_{31} \\ \Omega_{32} \\ \Omega_{33} \\ 0 \\ 0 \\ 0 \\ \varepsilon_\sigma u + \Omega_{31}u + \Omega_{32}v + \Omega_{33}w \end{bmatrix} $$
$$ \mathbf{H}_{nc} = \begin{bmatrix} -\alpha_1 \\ 0 \\ 0 \\ 0 \\ 0 \\ \alpha_1 P_1 \\ \alpha_2 P_2 \\ -c \\ 0 \end{bmatrix} $$

$$ \mathbf{H}_{relax} = \begin{bmatrix} \mu(P_1 - P_2) \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ -\mu P_I (P_1 - P_2) \\ \mu P_I (P_1 - P_2) \\ 0 \\ 0 \end{bmatrix} $$

It is important to note that the additional equation for the mixture total energy has been added to the system for numerical solution purposes. This equation is obviously in agreement with the complete system (14) and will be necessary to correct the energy conservation equation during the final relaxation step.

The unknown vector $\mathbf{U}^{n+1}$ is obtained from the initial condition $\mathbf{U}^n$ by application of the three successive operators according to the sequence:

$$ \mathbf{U}^{n+1} = L_{relax} L_{cap} L_{hyper} (\mathbf{U}^n). $$

Each step of the numerical method corresponds to the application of one of the three operators detailed below in a cartesian 2D framework.

### 4.2.1. Hyperbolic operator

The application of the first operator $L_{hyper}$ corresponds to the resolution of the hyperbolic submodel (16) using a Godunov-type method [15] extended to an high-order scheme with a MUSCL-Hancock procedure. The solution for this step is given for the cell $(i, j)$ by:

$$ \mathbf{U}_{i,j}^{\text{hyper}} = \mathbf{U}_{i,j}^{n} - \mathbf{G}\Delta t, $$

$$ \mathbf{G} = \left\{ \begin{array}{l} \frac{1}{\Delta x} \left( \mathbf{F}_{h; i+\frac{1}{2}, j}^* \left( \mathbf{U}_{i,j,R}^{n+\frac{1}{2}}, \mathbf{U}_{i+1,j,L}^{n+\frac{1}{2}} \right) - \mathbf{F}_{h; i-\frac{1}{2}, j}^* \left( \mathbf{U}_{i-1,j,R}^{n+\frac{1}{2}}, \mathbf{U}_{i,j,L}^{n+\frac{1}{2}} \right) \right) \\ + \frac{1}{\Delta y} \left( \mathbf{F}_{h; i,j+\frac{1}{2}}^* \left( \mathbf{U}_{i,j,T}^{n+\frac{1}{2}}, \mathbf{U}_{i,j+1,B}^{n+\frac{1}{2}} \right) - \mathbf{F}_{h; i,j-\frac{1}{2}}^* \left( \mathbf{U}_{i,j-1,T}^{n+\frac{1}{2}}, \mathbf{U}_{i,j,B}^{n+\frac{1}{2}} \right) \right) \\ + \frac{\mathbf{H}_{nc; i,j}^n}{\Delta x} \left( u_{i+\frac{1}{2}, j}^* \left( \mathbf{U}_{i,j,R}^{n+\frac{1}{2}}, \mathbf{U}_{i+1,j,L}^{n+\frac{1}{2}} \right) - u_{i-\frac{1}{2}, j}^* \left( \mathbf{U}_{i-1,j,R}^{n+\frac{1}{2}}, \mathbf{U}_{i,j,L}^{n+\frac{1}{2}} \right) \right) \\ + \frac{\mathbf{H}_{nc; i,j}^n}{\Delta y} \left( v_{i,j+\frac{1}{2}}^* \left( \mathbf{U}_{i,j,T}^{n+\frac{1}{2}}, \mathbf{U}_{i,j+1,B}^{n+\frac{1}{2}} \right) - v_{i,j-\frac{1}{2}}^* \left( \mathbf{U}_{i,j-1,T}^{n+\frac{1}{2}}, \mathbf{U}_{i,j,B}^{n+\frac{1}{2}} \right) \right) \end{array} \right\} . $$

The superscript * represents the solution of the Riemann problem on the corresponding cell boundary using the extrapolated values to the cell boundary $\mathbf{U}^{n+\frac{1}{2}}$ in the case of high order method (the subscripts R, L, T, B represent respectively the right, left, top and bottom neighbouring cells). These fluxes can be computed by any Riemann solver. Here a Harten-Lax-van Leer Contact (HLLC) approximate Riemann solver is used. Details of this method can be found in Saurel et al. [48]. Since the time marching scheme is explicit, the time step obeys a classical Courant-Friedrichs-Lewy (CFL) criterion.
### 4.2.2. Capillary operator

The second operator $L_{cap}$ introduces capillary terms and corresponds to the resolution of submodel (17):

$$ \mathbf{U}^{cap} = L_{cap}(\mathbf{U}^{hyper}). $$

Application of the operator $L_{cap}$ to vector $\mathbf{U}^{hyper}$ is done in the finite volume framework:

$$ \mathbf{U}_{i,j}^{cap} = \mathbf{U}_{i,j}^{hyper} - \Delta t \left\{ \frac{1}{\Delta x} \left( \mathbf{F}_{c; i+\frac{1}{2}, j}^{x} \left( \mathbf{U}_{i+\frac{1}{2}, j}^{hyper} \right) - \mathbf{F}_{c; i-\frac{1}{2}, j}^{x} \left( \mathbf{U}_{i-\frac{1}{2}, j}^{hyper} \right) \right) + \frac{1}{\Delta y} \left( \mathbf{F}_{c; i, j+\frac{1}{2}}^{y} \left( \mathbf{U}_{i, j+\frac{1}{2}}^{hyper} \right) - \mathbf{F}_{c; i, j-\frac{1}{2}}^{y} \left( \mathbf{U}_{i, j-\frac{1}{2}}^{hyper} \right) \right) \right\}. $$

The only equations affected by the capillary effects are the momentum and the total energy equations. These equations are developed in two dimensions:

$$ (\rho u)_{i,j}^{cap} = (\rho u)_{i,j}^{hyper} - \sigma \Delta t \left\{ -\frac{\|\mathbf{w}\|_{i+\frac{1}{2},j}^{hyper} - \|\mathbf{w}\|_{i-\frac{1}{2},j}^{hyper}}{\Delta x} + \frac{\left(\frac{w_1^2}{\|\mathbf{w}\|}\right)_{i+\frac{1}{2},j}^{hyper} - \left(\frac{w_1^2}{\|\mathbf{w}\|}\right)_{i-\frac{1}{2},j}^{hyper}}{\Delta y} \right\}, $$

$$ (\rho v)_{i,j}^{cap} = (\rho v)_{i,j}^{hyper} - \sigma \Delta t \left\{ -\frac{\|\mathbf{w}\|_{i,j+\frac{1}{2}}^{hyper} - \|\mathbf{w}\|_{i,j-\frac{1}{2}}^{hyper}}{\Delta y} + \frac{\left(\frac{w_1 w_2}{\|\mathbf{w}\|}\right)_{i+\frac{1}{2},j}^{hyper} - \left(\frac{w_1 w_2}{\|\mathbf{w}\|}\right)_{i-\frac{1}{2},j}^{hyper}}{\Delta x} \right\}, $$

$$ (\rho E + \varepsilon_{\sigma})_{i,j}^{cap} = (\rho E + \varepsilon_{\sigma})_{i,j}^{hyper} $$

$$ -\sigma \Delta t \left\{ \frac{\left(\frac{w_1^2 u}{\|\mathbf{w}\|}\right)_{i+\frac{1}{2},j}^{hyper} - \left(\frac{w_1^2 u}{\|\mathbf{w}\|}\right)_{i-\frac{1}{2},j}^{hyper}}{\Delta x} + \frac{\left(\frac{w_1 w_2 v}{\|\mathbf{w}\|}\right)_{i+\frac{1}{2},j}^{hyper} - \left(\frac{w_1 w_2 v}{\|\mathbf{w}\|}\right)_{i-\frac{1}{2},j}^{hyper}}{\Delta y} + \frac{\left(\frac{w_2^2 v}{\|\mathbf{w}\|}\right)_{i,j+\frac{1}{2}}^{hyper} - \left(\frac{w_2^2 v}{\|\mathbf{w}\|}\right)_{i,j-\frac{1}{2}}^{hyper}}{\Delta y} \right\}. $$

The different capillary terms at the cell boundary are obtained by the means of an arithmetic average of quantities of neighbouring cells. The vector $\mathbf{w}$ requires derivatives of the color function which are computed by using second-order finite difference approximations.

### 4.2.3. Pressure relaxation operator

To go through the operators chain, the solution at time $n + 1$ is obtained by a pressure relaxation algorithm and corrects the components of $\mathbf{U}^{cap}$:

$$ \mathbf{U}^{n+1} = L_{relax}(\mathbf{U}^{cap}). $$

The details about the pressure relaxation algorithm as well as the correction procedure used to guarantee total energy conservation can be found in Saurel et al. [48].
## 5. Numerical results and validations

In this section, 2D test cases are proposed to prove that the model and the numerical method are able to treat accurately both capillary effects and shock wave propagation.

In each presented cases, the equation of state (EOS) for the air obeys to the ideal gas law:

$$P_{air} = (\gamma_{air} - 1) \rho_{air} e_{air},$$

with $\gamma_{air} = 1.4$.

The liquid obeys the stiffened gas EOS:

$$P_{liquid} = (\gamma_{liquid} - 1) \rho_{liquid} e_{liquid} - \gamma_{liquid} P_{\infty,liquid},$$

where the stiffened gas EOS parameters are:

* $\gamma_{liquid} = 2.1$ and $P_{\infty,liquid} = 1.10^6 Pa$ for the capillary verification test cases of Section 5.1 and 5.2.

* $\gamma_{liquid} = 4.4$ and $P_{\infty,liquid} = 6.10^8 Pa$ for water in the compressible validation test case of Section 5.3.

### 5.1. Static capillary effects validation

2D tests are proposed to study the ability of the method to treat accurately capillary effects. The proposed method described in the present paper is compared with previous existing methods that use source terms integration. A particular attention is paid to the convergence, the numerical stability, the conservative property, the accuracy with respect to the Laplace pressure law and the dynamics of capillary flows. It is shown that the present method is able to simulate capillary effects in presence of pressure waves. In the following, “conservative method” will always mean the numerical method of Section 4 to distinguish it from the source terms integration method that do not conserve neither the momentum nor the energy.

The first test cases are static tests and they consist in finding the equilibrium state of a cylindrical column of liquid (2D liquid droplet) placed in air. The initial density in air is equal to $1 \text{ kg.m}^{-3}$ and $1000 \text{ kg.m}^{-3}$ in the liquid. The radius of the liquid column $R$ is set to $0.15 \text{ m}$ and the surface tension coefficient is equal to $800 \text{ N.m}^{-1}$. These unrealistic values are chosen to magnify the model properties. The pressure is initially uniform in the whole domain ($75 \text{ cm} \times 75 \text{ cm}$) and is set to $1 \text{ bar}$ (the initial conditions are presented in Figure 2).

#### 5.1.1. Accuracy regarding Laplace law and mesh convergence

The following test cases treating mesh convergence and accuracy regarding the Laplace law are performed using outgoing pressure wave boundary conditions. It consists in imposing the Neumann boundary conditions for the pressure and imposing in-outgoing conditions expressed in terms of the Riemann invariants.
<figure><img src="image_3.png" /><figcaption>Figure 2: Sketch of the initial conditions for the simulations of a cylindrical liquid column (2D liquid droplet) placed in air.</figcaption></figure>

Because of the capillary effects, the pressure of the liquid column converges to satisfy the Laplace law. The expression of the theoretical pressure jump of Laplace in 2 dimensions is:

$$[P] = \frac{\sigma}{R},$$

where $[P]$ expresses the pressure jump between inside and outside the droplet, here $P_{liquid} - P_{air}$.

The pressure convergence is analyzed using the relative residual:

$$\epsilon = Max \left( \frac{|P_{i,j}^n - P_{i,j}^{n-1}|}{P_{i,j}^n} \right).$$

The convergence is considered to be reached when the criterion $\epsilon \le 1.10^{-4}$ is verified. An example of evolution of this relative residual with a mesh of 120x120 cells is presented in Figure 3 for both conservative method and source terms integration method. It is clear that the pressure convergence is obtained faster using the conservative method than using the source terms integration one.

The pressure profiles through the liquid column is presented in Figure 4 for the mesh with 120x120 computational cells. Pressure profiles are presented for both methods after the pressure convergence time mentioned previously. The new conservative method has a better accuracy than the source terms integration one.

The mesh convergence study is performed on 4 different meshes containing 60x60 up to 120x120 computational cells. For each chosen mesh, the relative residual $\epsilon$ is used to match
<figure><img src="image_4.png" /><figcaption>Figure 3: Relative residual $\epsilon$ function of the physical time for the 2D liquid column test case with a mesh of 120x120 cells. Source terms integration (red) and conservative (blue) methods. Convergence is reached after 2 seconds with the conservative method compared to more than 6 seconds with the source terms integration method.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Pressure jump in a centered section of a steady 2D liquid column in air. $\sigma = 800\,\text{N.m}^{-1}$, $R = 0.1496\,\text{m}$ and so [$P$] = 5347 Pa. Theoretical (black lines) and converged simulations results with source terms integration (red crosses) and conservative (blue dots) methods are shown at times 6.68s and 2.07s, respectively.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Relative pressure error $\xi$ averaged in time when the convergence is reached versus the number of cells in one direction $N$ for the 2D liquid column test case. Grid convergence toward the theoretical pressure jump of Laplace is shown for the conservative method proposed in the paper. The bars show the domain of variation of $\xi$ in time after the convergence criteria is reached, associated to each mesh. The interval of variation always decreases with the cells refinement.</figcaption></figure>

pressure convergence. Then the relative pressure error $\xi$ is computed by:

$$ \xi = \frac{|[P_{num}] - [P_{th}]|}{[P_{th}]}, $$

where $[P_{num}]$ is the pressure jump between the averaged pressure inside the column ($P_{num} = \sum \alpha_{liquid}P_{liquid} / \sum \alpha_{liquid}$ when $\alpha_{liquid} \ge 0.99$) and the pressure of the air $P_{air}$. The theoretical pressure jump $[P_{th}]$ is calculated using the following relation:

$$ [P_{th}] = \sigma \sqrt{\frac{\pi}{S}}, $$

where $S$ represents the 2D volume of the droplet. It is determined doing the sum of the volume fraction of liquid over the whole domain multiplied by the cell volume: $S = \sum \alpha_{liquid}S_{cell}$. The relative pressure error $\xi$ is averaged in time and plotted versus grid resolution in Figure 5 for the conservative method proposed in the paper. Grid convergence is validated and the variation of $\xi$ over time, when convergence is reached, decreases with the cells refinement.

### 5.1.2. Stability

Regarding the stability of the method, contours of the liquid column are presented in Figure 6. Wall boundary conditions are used in this example on a mesh containing 120x120
<figure><img src="image_7.png" /><figcaption>Figure 6: Magnified view of column contours through color function ($c = 0.5$). Initial contour (black contour) together with cartesian mesh is on the left hand side; Contours after convergence is reached: source terms integration method is in the middle (red contour) and conservative method is on the right hand side (blue contour) for physical times about 3.68s and 2.32s, respectively. With source terms integration method, the position of the column as already changed when convergence is reached.</figcaption></figure>

cells. The left hand side of the figure shows the initial contour together with the mesh. Effect of the mesh on the initial circle is noticeable. The two other pictures show contour positions immediately after the convergence is reached, corresponding to a physical time of about 3.68s for the source term integration method and 2.32s for the conservative method. In this test case, if the source term integration method is used, one can observe a slight move of the gravity center of the column before a converged solution is reached. In fact the gravity center starts to move at a physical time about 2.57s. At the converged time (3.68s), the position of the gravity center is still incorrect (picture in the middle of Figure 6). This drift may be considered as a visible result of the non-conservation of the momentum. This is crippling, for example, when studying interaction between a droplet and a shock wave (see the last subsection): the droplet must be at rest before the interaction with the shock wave. With the conservative method of the present paper, it becomes possible to consider the treatment of this kind of problem.

Nevertheless, the gravity center starts to move after a long time in the two methods. In this example, with the conservative method, it is about 4.52s so well after the convergence is reached (2.32s). But these displacements cannot be compared with the drift observed for the source terms integration method. In other words, one can consider that the new method maintains the position of the gravity center at the converged solution. For longer time, with the accumulation of rounding errors, some spurious motions can be seen as traces of numerical instabilities. They do not have the same origin as the drift observed with the source terms integration method as soon as the calculation begins.

### 5.1.3. *Energy conservation*

The conservative property of the new model is verified by plotting the evolution of the relative error for the total mixture energy over the whole domain for the liquid column using wall boundary conditions with a mesh containing 120x120 cells (Figure 7). The energy is obviously conserved with the conservative method but not with the source terms integration
<figure><img src="image_8.png" /><figcaption>Figure 7: Evolution of the relative error for the total mixture energy over the whole domain for the liquid column with wall boundary conditions. Expression for the evolution of the relative error is: $100(E(t) - E(t_0))/E(t_0)$. Source terms integration (red) and conservative (blue) methods.</figcaption></figure>

one. Moreover, for the source terms integration method, the error does not saturate and then still increases when the convergence is reached. It is also noticeable that the difference in energies appears even initially and comes from the capillary energy $\epsilon_σ$ which does not appear in the source term form expression of the total mixture energy equation. Note that its value is here: $\epsilon_σ = 2.10^7 J.m^{-3}$.

## 5.2. Dynamic capillary effects validation

The next two test cases for the dynamic verification of surface tension are slightly different from the previous ones.

### 5.2.1. Recovering a circle shape

In this test case we are looking for the shape recovering and the volume evolution of an initially square liquid column, with $L_{square} = 0.2m$, placed in an atmospheric air with the same pressure $P_{air}$. Because of the capillary terms the solution should evolve to a circle shape of the liquid interface and the pressure should verify the Laplace jump: $P_{final} = P_{air} + \sigma/R$. The successive dynamic stages leading to the converged circle shape are presented in Figure 8.

Moreover, the mass conservation implies a variation of the liquid column volume. This variation is following an isentropic transformation leading to the final volume:

$$V_{final} = V_{initial} \left( \frac{P_{initial} + P_{\infty}}{P_{final} + P_{\infty}} \right)^{\frac{1}{\gamma}}$$
<figure><img src="image_9.png" /><figcaption>Figure 8: Magnified schlieren images of the mixture density of the 2D liquid column under a square shape in an air environment converging to a circle shape. Initially $\sigma = 800N.m^{-1}$ and $L_{square} = 0.2m$. Simulation at different times with a mesh of 150x150 cells: (a) initialization, (b) 5ms, (c) 15ms, (d) 20ms, (e) 30ms, (f) 35ms, (g) 45ms, (h) 50ms, (i) 105ms and (j) 3000ms.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Evolution of the relative variation of the 2D volume of the liquid column ((V - Vfinal)/Vfinal) for 3 different meshes. The final isentropic 2D volume (S = 0.039878m²) is numerically recovered with an error of 0.011% for the mesh with 200x200 cells.</figcaption></figure>

Initially the area of the square column is $V_{initial} = 0.04m^2$ and it converges to the theoretical value of $V_{final} = 0.039878m^2$ with a precision equal to 0.011% at time $t = 3s$ with a mesh made of 200x200 cells. In Figure 9, the relative variation of the liquid 2D volume ((V - Vfinal)/Vfinal) is plotted versus time for 3 different meshes. Grid convergence toward the theoretical value is observed.

### 5.2.2. Uniform velocity flow

In this section, we focus on the ability of the method to maintain a cylindrical column with a correct pressure jump in a uniform velocity flow. A sketch of the initial condition is presented in Figure 10. It consists in a liquid column with an initial radius of $R = 0.15m$ moving in a uniform velocity field ($u = 50m.s^{-1}$). At the initialization, the liquid column is in equilibrium with the ambient air. Simulations are performed on 4 different meshes containing 180x60 up to 360x120 computational cells for a physical domain of $2.25m \times 0.75m$.

In Figure 11, the relative pressure error $\xi$ defined previously is averaged in time and plotted versus grid resolution. Mesh convergence is observed and results show that the Laplace pressure jump is well verified.

During the simulations, the position of the liquid column is in perfect agreement with the theoretical one (Figure 12). It is then important to note that surface tension effects have influence only on the shape of the column, since initially there is a grid impact (see Figure 6 for an example), but not on the position. This validates the method for this simple advection test case.
<figure><img src="image_11.png" /><figcaption>Figure 10: Sketch of the initial conditions for the simulation of a cylindrical liquid column (2D liquid droplet) in air, coupled with advection of the whole domain.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11: Relative pressure error $\xi$ averaged in time versus the number of cells in the y-direction $N_y$ for the 2D liquid column test case with surface tension and advection of the whole domain. Mesh convergence is observed. The results show that the Laplace pressure jump is well verified using the new method.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 12: Evolution of the position of the liquid column along the *x* axis minus its initial position for the 2D liquid column test case with surface tension and advection of the whole domain. Numerical results using the new method are in perfect agreements with the theoretical one.</figcaption></figure>

### 5.3. Shock wave interaction with a water column

We focus here on the early stages of droplet breakup in the high speed flow behind a shock wave.

The numerical simulations are performed in two dimensions to simulate the interaction of an air flow and a water column. Comparisons with Igra and Takayma experiments [21, 22] are carried out.

A cylindrical water column with initial diameter of $D = 6.4\,\text{mm}$ is exposed to a shock wave at Mach number 1.3 in atmospheric air (see Figure 13 for initialization sketch). The initial densities are $\rho_{\text{air}} = 1.2\,\text{kgm}^{-3}$ and $\rho_{\text{water}} = 1000\,\text{kg.m}^{-3}$. The corresponding initial Weber number ($We = \rho_{\text{air}}u^2D/\sigma$) in these conditions is 3690. 2D computations are performed on a cartesian mesh containing 3200x1200 cells representing a physical domain of 220mm x 82.5mm. Shocked air is entering at the left. Outgoing pressure wave boundary condition is used at the right and the Newmann boundary conditions is used at the top and the bottom.

The results are presented at three different stages. Each stage is representative of a typical physical effect:

* The first stage is presented in Figure 14 where the flow dynamics is governed by the wave pattern. This is a totally compressible stage corresponding to the shock-column interaction for a large Weber number. Comparisons with experiments are presented in Figure 15.

* After the shock wave propagation, the column shape is changing due to inertial effects.
<figure><img src="image_14.png" /><figcaption>Figure 13: Sketch of the initialization for interaction of a shock wave with a water column.</figcaption></figure>

Liquid-gas interface is stretched and the first liquid filaments are appearing. The interface location is compared with experiments in Figure 16.

* At the third stage, the breakup begins because of high inertia of the water column. This effect has already been observed in bubble breakup (Yang et al. [55], Quirk and Karni [44], Layes and Le Metayer [28]). The Weber number is lower (due to diameter change and also because the velocity gap is reducing) and the capillary effects tend to tear the filaments out the main column (last pictures of Figure 17).

In Figure 18, results obtained for a flow without capillary effects (left picture) and a flow with capillary effects (right picture) are compared at time $t = 1.2ms$. One can see that the capillary effects begin to have influence on the flow topology, particularly on the filaments shape where the surface tension effects tend to break the filaments to form little droplets.

## Conclusion

A new multiphase model treating interface problems and capillary effects has been derived. Hyperbolicity study has been completely done and a splitting numerical method that guarantees conservation of the mass, the momentum and the energy has been built. Comparison with “source terms integration method” and with experiments have shown the advantages of the new model and numerical method. The simulation of an aerodynamic breakup of a water column induced by a shock wave showed that the method is able to treat accurately both pressure waves interaction and capillary effects in the same formulation. Future works will include viscous effects, heat conduction effects and massive simulations including a large number of droplets in 3D.
<figure><img src="image_15.png" /><figcaption>Figure 14: Magnified view of the early stages of the aerodynamic breakup of the water column in the high-speed airstream behind the shock wave of Mach number 1.3 in atmospheric air propagating from the left to the right. 2D simulation with a mesh of 3200x1200 cells and schlieren images of the mixture density are given at the following times: (a) 0µs, (b) 10µs, (c) 20µs, (d) 24µs, (e) 37µs and (f) 55µs.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 15: Magnified view of the early stages of the aerodynamic breakup of the water column in the high-speed airstream behind the shock wave of Mach number 1.3 in atmospheric air propagating from the left to the right. 2D simulation on the left with a mesh of 3200x1200 cells and a schlieren image of the mixture density. Experiment on the right from Igra and Takayama [22] with an interferogram on the upper half and an unreconstructed hologram on the lower half. Results are shown at times 15µs (top) and 55µs (bottom).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 16: Magnified view of the aerodynamic breakup of the water column in the high-speed airstream behind the shock wave of Mach number 1.3 in atmospheric air propagating from the left to the right. 2D simulation on the left with a mesh of 3200x1200 cells and a schlieren image of the mixture density. Experiment on the right from Igra and Takayama [22]. Results are shown at times 200µs (top) and 300µs (bottom).</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 17: Magnified view of the aerodynamic breakup of the water column in the high-speed airstream behind the shock wave of Mach number 1.3 in atmospheric air propagating from the left to the right. 2D simulation with a mesh of 3200x1200 cells and schlieren images of the mixture density (upper half) coupling with volume fraction of water visualization (lower half, in white the air and in blue the water) are given at the following times: (a) 160µs, (b) 320µs, (c) 420µs, (d) 500µs, (e) 660µs, (f) 820µs, (g) 1000µs and (h) 1200µs.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 18: Magnified view of the aerodynamic breakup of the water column in the high-speed airstream behind the shock wave of Mach number 1.3 in atmospheric air propagating from the left to the right. 2D simulations respectively without (left) and with (right) capillary effects taken into account, with a mesh of 3200x1200 cells and a schlieren image of the mixture density. Results are shown at time 1.2ms.</figcaption></figure>

# Appendix I: Variational principle for compressible mixtures with capillary effects

Consider a continuum characterized by the Lagrangian:

$$L = \rho \left( \frac{\|\mathbf{u}\|^2}{2} - e \right) - \sigma \|\nabla c\|.$$

Here $\rho = \alpha_1\rho_1 + \alpha_2\rho_2$ is the mixture density, $e = Y_1e_1 + Y_2e_2$ is the mixture specific energy, $Y_i = \alpha_i\rho_i/\rho$ are the mass fractions, $\sigma$ is the surface tension coefficient, $c$ is the color function.

Consider the Hamilton action:

$$a = \int_{t_1}^{t_2} \int_{D} L dtdD, \qquad (18)$$

where $t_1$ and $t_2$ are fixed time instants, and $D$ is a material fluid volume.

The governing equations are extremal curves of the Hamilton action under the following constraints:

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0,$$

$$\frac{\partial Y_1}{\partial t} + \mathbf{u} \cdot \nabla (Y_1) = 0,$$

$$\frac{\partial s_1}{\partial t} + \mathbf{u} \cdot \nabla (s_1) = 0,$$
$$ \begin{aligned} & \frac{\partial s_2}{\partial t} + \mathbf{u} \cdot \nabla (s_2) = 0, \\ & \frac{\partial c}{\partial t} + \mathbf{u} \cdot \nabla (c) = 0. \end{aligned} \qquad (19) $$

The Eulerian variations of the unknown variables in terms of virtual displacement $\delta \mathbf{x}$ are given by (see Gavrilyuk [13] for details):

$$ \begin{aligned} & \delta \rho = -\nabla \cdot (\rho \delta \mathbf{x}), \\ & \delta Y_1 = -\nabla Y_1 \cdot \delta \mathbf{x}, \\ & \delta s_1 = -\nabla s_1 \cdot \delta \mathbf{x}, \\ & \delta s_2 = -\nabla s_2 \cdot \delta \mathbf{x}, \\ & \delta c = -\nabla c \cdot \delta \mathbf{x}, \\ & \delta \mathbf{u} = \frac{\partial \delta \mathbf{x}}{\partial t} + \frac{\partial \delta \mathbf{x}}{\partial \mathbf{x}} \cdot \mathbf{u} - \frac{\partial \mathbf{u}}{\partial \mathbf{x}} \cdot \delta \mathbf{x}. \end{aligned} $$

Taking the variation of the Hamilton action (18) under the boundary conditions:

$$ \begin{aligned} & \delta \mathbf{x}|_{\partial D} = 0, \\ & \delta \mathbf{x}|_{t=t_1} = 0, \\ & \delta \mathbf{x}|_{t=t_2} = 0, \end{aligned} $$

one can obtain the momentum equation:

$$ \rho \left( \frac{\partial \mathbf{u}}{\partial t} + \frac{\partial \mathbf{u}}{\partial \mathbf{x}} \cdot \mathbf{u} \right) + \nabla P = -\sigma (\nabla \cdot \mathbf{n}) \nabla c, $$

with $\mathbf{n} = \nabla c / \|\nabla c\|$.

The equation is exactly the same as in Brackbill et al. [3] for the one component case. It can be written in conservative form:

$$ \frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot \left( \rho \mathbf{u} \otimes \mathbf{u} + P \bar{\bar{I}} + \bar{\bar{\Omega}} \right) = 0, $$

where the capillary tensor $\bar{\bar{\Omega}}$ is given by:

$$ \bar{\bar{\Omega}} = -\sigma \| \nabla c \| (\bar{\bar{I}} - \mathbf{n} \otimes \mathbf{n}). $$

The variation of the Hamilton action with respect to the volume fraction gives as the pressure equilibrium condition:

$$ P_1 = P_2 = P = \alpha_1 P_1 + \alpha_2 P_2. $$

resulting to the non-conservative equation (2) for $\alpha_1$ (see Gavrilyuk [13] for details).
Using the constraints on the entropies and mass fractions of each phase, the mixture entropy equation is given by:

$$ \frac{\partial s}{\partial t} + \mathbf{u} \cdot \nabla (s) = 0. $$

The mixture entropy, mass and momentum equations admit the following mixture total energy equation:

$$ \frac{\partial \rho E}{\partial t} + \nabla \cdot ((\rho E + P) \mathbf{u}) = - (\nabla \cdot \overline{\Omega}) \cdot \mathbf{u}. $$

To transform the previous equation to a fully conservative form, the following relation is used:

$$ (\nabla \cdot \overline{\Omega}) \cdot \mathbf{u} = \nabla \cdot (\overline{\Omega} \cdot \mathbf{u}) - \overline{\Omega} : \frac{\partial \mathbf{u}}{\partial \mathbf{x}}, $$

where the right term still needs to be transformed. It can be developed as:

$$ \sigma \left( \| \nabla c \| \overline{I} - \frac{\nabla c \otimes \nabla c}{\| \nabla c \|} \right) : \frac{\partial \mathbf{u}}{\partial \mathbf{x}} = \sigma \left( \| \nabla c \| \nabla \cdot \mathbf{u} - \left( \nabla c \cdot \frac{\partial \mathbf{u}}{\partial \mathbf{x}} \right) \cdot \frac{\nabla c}{\| \nabla c \|} \right). \quad (20) $$

Using the material derivative equation of $\| \nabla c \|$:

$$ \frac{d \| \nabla c \|}{dt} = \frac{d \nabla c}{dt} \cdot \frac{\nabla c}{\| \nabla c \|}, \quad (21) $$

with the material derivative operator $d(\cdot)/dt = \partial(\cdot)/\partial t + \mathbf{u} \cdot \nabla(\cdot)$ and the Schwarz theorem, we have:

$$ \frac{d\nabla c}{dt} = \nabla \left(\frac{dc}{dt}\right) - \nabla c \cdot \frac{\partial \mathbf{u}}{\partial \mathbf{x}}. \quad (22) $$

Hence with equations (19) and (22), equation (21) is rewritten as:

$$ \frac{d \| \nabla c \|}{dt} = - \left( \nabla c \cdot \frac{\partial \mathbf{u}}{\partial \mathbf{x}} \right) \cdot \frac{\nabla c}{\| \nabla c \|}. $$

Thus relation (20) is now developed as:

$$ \sigma \left( \| \nabla c \| \overline{I} - \frac{\nabla c \otimes \nabla c}{\| \nabla c \|} \right) : \frac{\partial \mathbf{u}}{\partial \mathbf{x}} = \sigma \left( \| \nabla c \| \nabla \cdot \mathbf{u} + \frac{d \| \nabla c \|}{dt} \right) \\ = \frac{\partial \sigma \| \nabla c \|}{\partial t} + \nabla \cdot (\sigma \| \nabla c \| \mathbf{u}). $$

Finally, the energy conservation law is:

$$ \frac{\partial \rho E + \varepsilon_{\sigma}}{\partial t} + \nabla \cdot \left( (\rho E + \varepsilon_{\sigma} + P) \mathbf{u} + \overline{\Omega} \cdot \mathbf{u} \right) = 0, $$

with the capillary energy term $\varepsilon_{\sigma} = \sigma \| \nabla c \|$.
## Appendix II: Eigenvectors of model (3)

Remind that model (3) can be written under a vector form (6) with the vector $\mathbf{W}$ and the matrix $\bar{\mathbf{A}}$ defined by:

$$ \mathbf{W} = [\alpha_1, u, v, w, P, w_1, w_2, w_3, s_1, s_2, Y_1, c]^T, $$

$$ \bar{\mathbf{A}}(\mathbf{W}) = \begin{bmatrix} u & -K & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & u & 0 & 0 & \frac{1}{\rho} & \frac{1}{\rho} \frac{\partial \Omega_{11}}{\partial w_1} & \frac{1}{\rho} \frac{\partial \Omega_{11}}{\partial w_2} & \frac{1}{\rho} \frac{\partial \Omega_{11}}{\partial w_3} & 0 & 0 & 0 \\ 0 & 0 & u & 0 & 0 & \frac{1}{\rho} \frac{\partial \Omega_{12}}{\partial w_1} & \frac{1}{\rho} \frac{\partial \Omega_{12}}{\partial w_2} & \frac{1}{\rho} \frac{\partial \Omega_{12}}{\partial w_3} & 0 & 0 & 0 \\ 0 & 0 & 0 & u & 0 & \frac{1}{\rho} \frac{\partial \Omega_{13}}{\partial w_1} & \frac{1}{\rho} \frac{\partial \Omega_{13}}{\partial w_2} & \frac{1}{\rho} \frac{\partial \Omega_{13}}{\partial w_3} & 0 & 0 & 0 \\ 0 & \rho a^2 & 0 & 0 & u & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & w_1 & w_2 & w_3 & 0 & u & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & u & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & u & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & u & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & u & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & u \end{bmatrix} $$

The corresponding eigenvalues are real, explicit and given by:

$$ \lambda_{1,2,3,4,5,6,7,8} = u, $$

$$ \lambda_{9,10} = u \pm a_s, $$

$$ \lambda_{11,12} = u \pm a_c, $$

where:

$$ a_s^2 = \frac{a^2 + b + \sqrt{(a^2 + b)^2 - 4a^2b(n_2^2 + n_3^2)}}{2}, $$

$$ a_c^2 = \frac{a^2 + b - \sqrt{(a^2 + b)^2 - 4a^2b(n_2^2 + n_3^2)}}{2}, $$

$a$ is the Wood mixture speed of sound previously defined (5) and:

$$ b = \frac{\sigma \| \mathbf{w} \|}{\rho} (n_2^2 + n_3^2). $$

The system is hyperbolic if the multiple eigenvalues $\lambda = u$ have exactly 8 linearly independent eigenvectors $\mathbf{R}$. 5 of them are straightforward:
$$ \mathbf{R}_1 = \begin{bmatrix} 1 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} $$

$$ \mathbf{R}_2 = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 1 \\ 0 \\ 0 \\ 0 \end{bmatrix} $$

$$ \mathbf{R}_3 = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 1 \\ 0 \\ 0 \end{bmatrix} $$

$$ \mathbf{R}_4 = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 1 \\ 0 \end{bmatrix} $$

$$ \mathbf{R}_5 = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 1 \end{bmatrix} $$

And the 3 last eigenvectors for $\lambda = u$ have to be found using the following reduced matrix:

$$ \bar{\bar{a}} = \begin{bmatrix} u & 0 & 0 & \frac{1}{\rho} & \frac{\partial \Omega_{11}}{\rho \partial w_1} & \frac{\partial \Omega_{11}}{\rho \partial w_2} & \frac{\partial \Omega_{11}}{\rho \partial w_3} \\ 0 & u & 0 & 0 & \frac{\partial \Omega_{12}}{\rho \partial w_1} & \frac{\partial \Omega_{12}}{\rho \partial w_2} & \frac{\partial \Omega_{12}}{\rho \partial w_3} \\ 0 & 0 & u & 0 & \frac{\partial \Omega_{13}}{\rho \partial w_1} & \frac{\partial \Omega_{13}}{\rho \partial w_2} & \frac{\partial \Omega_{13}}{\rho \partial w_3} \\ \rho a^2 & 0 & 0 & u & 0 & 0 & 0 \\ w_1 & w_2 & w_3 & 0 & u & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & u & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & u \end{bmatrix} $$

The corresponding reduced eigenvectors $\mathbf{r}$ are the solutions of $(\bar{\bar{a}} - u\bar{\bar{I}}) \mathbf{r} = \mathbf{0}$, where $\mathbf{r}$ is written as:

$$ \mathbf{r} = [r_1, r_2, r_3, r_4, r_5, r_6, r_7]^T $$

3 distinct systems of equations results from this equality:

$$ \rho a^2 r_1 = 0, $$
$$w_1 r_1 + w_2 r_2 + w_3 r_3 = 0,$$

$$\begin{bmatrix} \frac{1}{\rho} r_4 \\ 0 \\ 0 \end{bmatrix} + \bar{C} \begin{bmatrix} r_5 \\ r_6 \\ r_7 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix}$$

where $\bar{C}$ is the matrix with the only capillary terms:

$$\bar{C} = \begin{bmatrix} \frac{1}{\rho} \frac{\partial \Omega_{11}}{\partial w_1} & \frac{1}{\rho} \frac{\partial \Omega_{11}}{\partial w_2} & \frac{1}{\rho} \frac{\partial \Omega_{11}}{\partial w_3} \\ \frac{1}{\rho} \frac{\partial \Omega_{12}}{\partial w_1} & \frac{1}{\rho} \frac{\partial \Omega_{12}}{\partial w_2} & \frac{1}{\rho} \frac{\partial \Omega_{12}}{\partial w_3} \\ \frac{1}{\rho} \frac{\partial \Omega_{13}}{\partial w_1} & \frac{1}{\rho} \frac{\partial \Omega_{13}}{\partial w_2} & \frac{1}{\rho} \frac{\partial \Omega_{13}}{\partial w_3} \end{bmatrix}$$

A solution of the first two equations is: $r_1 = 0$, $r_2 = -\eta w_3$ and $r_3 = \eta w_2$, where $\eta$ is any real constant. If the matrix $\bar{C}$ is invertible, the components $r_5$, $r_6$ and $r_7$ are uniquely defined for a given $r_4$. First, the derivatives of the capillary tensor terms are:

$$\frac{\partial \Omega_{11}}{\partial w_1} = \frac{\sigma}{\rho} n_1 (n_2^2 + n_3^2),$$

$$\frac{\partial \Omega_{11}}{\partial w_2} = \frac{\sigma}{\rho} n_2 (n_2^2 + n_3^2 - 2),$$

$$\frac{\partial \Omega_{11}}{\partial w_3} = \frac{\sigma}{\rho} n_3 (n_2^2 + n_3^2 - 2),$$

$$\frac{\partial \Omega_{12}}{\partial w_1} = \frac{\sigma}{\rho} n_2 (n_2^2 + n_3^2),$$

$$\frac{\partial \Omega_{12}}{\partial w_2} = \frac{\sigma}{\rho} n_1 (n_1^2 + n_3^2),$$

$$\frac{\partial \Omega_{12}}{\partial w_3} = \frac{\sigma}{\rho} n_1 n_2 n_3,$$

$$\frac{\partial \Omega_{13}}{\partial w_1} = \frac{\sigma}{\rho} n_3 (n_2^2 + n_3^2),$$

$$\frac{\partial \Omega_{13}}{\partial w_2} = \frac{\sigma}{\rho} n_1 n_2 n_3,$$

$$\frac{\partial \Omega_{13}}{\partial w_3} = \frac{\sigma}{\rho} n_1 (n_1^2 + n_2^2).$$

And then the determinant of the matrix $\bar{C}$ is:

$$\det(\bar{C}) = \left(\frac{\sigma}{\rho}\right)^3 n_1(n_2^2(n_1^2 + n_2^2)^3 + n_3^2(n_1^2 + n_3^2)^3 \\ + 6n_1^4 n_2^2 n_3^2 + 9n_1^2 n_2^4 n_3^2 + 4n_2^6 n_3^2 + 9n_1^2 n_2^2 n_3^4 + 6n_2^4 n_3^4 + 4n_2^2 n_3^6).$$
The determinant is non-zero as long as $n_1$ is non-vanishing, which is numerically always the case. Thus:

$$ \begin{bmatrix} r_5 \\ r_6 \\ r_7 \end{bmatrix} = -r_4 \bar{C}^{-1} \begin{bmatrix} \frac{1}{\rho} \\ 0 \\ 0 \end{bmatrix} $$

Finally the resulting reduced eigenvectors are:

$$ \mathbf{r} = \begin{bmatrix} 0 \\ -\eta w_3 \\ \eta w_2 \\ r_4 \\ -r_4 \bar{C}^{-1} \begin{bmatrix} \frac{1}{\rho} \\ 0 \\ 0 \end{bmatrix} \end{bmatrix} $$

Thus, only a two-parameter family of eigenvectors exists. Linearly independent eigenvectors can be given as:

$$ \mathbf{R}_6 = \begin{bmatrix} 0 \\ 0 \\ -w_3 \\ w_2 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} $$

$$ \mathbf{R}_7 = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 1 \\ -\bar{C}^{-1} \begin{bmatrix} \frac{1}{\rho} \\ 0 \\ 0 \end{bmatrix} \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} $$

and thus 1 eigenvector is missing to have a hyperbolic model. Model (3) is then weakly hyperbolic.

## Acknowledgments

Authors are particularly grateful to Dan Igra for providing experimental data and to Sarah Hank for fruitful discussions.

## References

[1] R. Abgrall and S. Karni. Computations of compressible mutlifluids. *Journal of Computational Physics*, 169(2):594–623, 2001.
[2] G. Allaire, S. Clerc, and S. Kokh. A five-equation model for the simulation of interfaces between compressible fluids. *J. Comp. Phys.*, 181:577–616, 2002.

[3] J. Brackbill, D. Kothe, and C. Zemach. A continuum method for modelling surface tension. *Journal of Computational Physics*, 100:335–354, 1992.

[4] A. Chauvin, E. Daniel, A. Chinnayya, J. Massoni, and G. Jourdan. Shock waves in sprays: numerical study of secondary atomization and experimental comparison. *Shock Waves*, pages 1–13, 2015.

[5] A. Chauvin, G. Jourdan, E. Daniel, L. Houas, and R. Tosello. Experimental investigation of the propagation of a planar shock wave through a two-phase gas-liquid medium. *Physics of Fluids (1994-present)*, 23(11):113301, 2011.

[6] S. Chen and G.D. Doolen. Lattice boltzmann method for fluid flows. *Annual review of fluid mechanics*, 30(1):329–364, 1998.

[7] B.M. Devassy, C. Habchi, and E. Daniel. Atomization modelling of liquid jets using a two-surface-density approach. *Atomization and Sprays*, 25(1):47–80, 2015.

[8] O.G. Engel. Fragmentation of waterdrops in the zone behind an air shock. *J. Res. Natl. Bur. Stand*, 60(3):245–280, 1958.

[9] N. Favrie and S.L. Gavrilyuk. Diffuse interface model for compressible fluid–compressible elastic–plastic solid interaction. *Journal of Computational Physics*, 231(7):2695–2723, 2012.

[10] N. Favrie and S.L. Gavrilyuk. Dynamic compaction of granular materials. *Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences*, 469(2160):20130214, 2013.

[11] N. Favrie, S.L. Gavrilyuk, and S. Ndanou. A thermodynamically compatible splitting procedure in hyperelasticity. *Journal of Computational Physics*, 270:300–324, 2014.

[12] R.P. Fedkiw, T. Aslam, B. Merriman, and S. Osher. A non oscillatory Eulerian approach to interfaces in multimaterial flows (The Ghost Fluid Method). *Journal of Computational Physics*, 152:457–492, 1999.

[13] S.L. Gavrilyuk. Multiphase flow modeling via Hamilton's principle. In the book: F. Dell'Isola and S.L. Gavrilyuk, Eds. *Variational models and methods in solid and fluid mechanics, CISM Courses and Lectures*, volume 535. Springer, 2012.

[14] S.L. Gavrilyuk and R. Saurel. Rankine-Hugoniot relations for shocks in heterogeneous mixtures. *Journal of Fluid Mechanics*, 575:495–507, 2007.

[15] S.K. Godunov. A finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics. *Math Sb.*, 47:357–393, 1959.
[16] D. Gueyffier, L. Li, A. Nadim, R. Scardovelli, and S. Zaleski. Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows. *Journal of Computational Physics*, 152:423-456, 1999.

[17] L.H. Han, X.Y. Hu, and N.A. Adams. Scale separation for multi-scale modeling of free-surface and two-phase flows with the conservative sharp interface method. *Journal of Computational Physics*, 280:387-403, 2015.

[18] L-P. Hsiang and G.M. Faeth. Near-limit drop deformation and secondary breakup. *International Journal of Multiphase Flow*, 18(5):635-652, 1992.

[19] X.Y. Hu, B.C. Khoo, N.A. Adams, and F.L. Huang. A conservative interface method for compressible flows. *Journal of Computational Physics*, 219(2):553-578, 2006.

[20] D. Igra and K. Takayama. Investigation of aerodynamic breakup of a cylindrical water droplet. *Atomization and Sprays*, 11(2):167-185, 2001.

[21] D. Igra and K. Takayama. Numerical simulation of shock wave interaction with a water column. *Shock Waves*, 11(3):219-228, 2001.

[22] D. Igra and K. Takayama. A study of shock wave loading on a cylindrical water column. *Report of the Institute of Fluid Science*, Tohoku University, 13:19-36, 2001.

[23] D.D. Joseph, J. Belanger, and G.S. Beavers. Breakup of a liquid drop suddenly exposed to a high-speed airstream. *International Journal of Multiphase Flow*, 25(6):1263-1303, 1999.

[24] A. Kapila, R. Menikoff, J. Bdzil, S. Son, and D. Stewart. Two-phase modeling of DDT in granular materials: Reduced equations. *Physics of Fluids*, 13:3002-3024, 2001.

[25] S. Karni. Multicomponent flow calculations by a consistent primitive algorithm. *Journal of Computational Physics*, 112(1):31-43, 1994.

[26] S. Karni. Hybrid multifluid algorithms. *SIAM Journal on Scientific Computing*, 17(5):1019-1039, 1996.

[27] B. Koren, M.R. Lewis, E.H. van Brummelen, and B. van Leer. Riemann-problem and level-set approaches for homotropic two-fluid computations. *Journal of Computational Physics*, 181:654-674, 2002.

[28] G. Layes and O. Le Metayer. Quantitative numerical and experimental studies of the shock accelerated heterogeneous bubbles motion. *Physics of Fluids*, 19(042105), 2007.

[29] S. Le Martelot, R. Saurel, and B. Nkonga. Towards the direct numerical simulation of nucleate boiling flows. *International Journal of Multiphase Flow*, 66:62-78, 2014.

[30] T.G. Liu, B.C. Khoo, and K.S. Yeo. Ghost fluid method for strong shock impacting on material interface. *Journal of Computational Physics*, 190(2):651-681, 2003.
[31] J. Luo, X.Y. Hu, and N.A. Adams. A conservative sharp interface method for incompressible multiphase flows. *Journal of Computational Physics*, 284:547-565, 2015.

[32] J. Massoni, R. Saurel, B. Nkonga, and R. Abgrall. Proposition de methodes et modeles Euleriens pour les problemes a interfaces entre fluides compressibles en presence de transfert de chaleur. *Int. J. Heat and Mass Transfer*, 45:1287-1307, 2002.

[33] J.C. Meng and T. Colonius. Numerical simulations of the early stages of high-speed droplet breakup. *Shock Waves*, pages 1-16, 2014.

[34] A. Murrone and H. Guillard. Behavior of upwind scheme in the low mach number limit: Iii. preconditioned dissipation for a five equation two phase model. *Computers & Fluids*, 37(10):1209-1224, 2008.

[35] S. Ndanou, N. Favrie, and S.L. Gavrilyuk. Criterion of hyperbolicity in hyperelasticity in the case of the stored energy in separable form. *Journal of Elasticity*, 115(1):1-25, 2014.

[36] S. Osher and R. Fedkiw. Level set methods: An overview and some recent results. *Journal of Computational Physics*, 169:463-502, 2001.

[37] S. Osher and R. Fedkiw. *Level set methods and dynamic implicit surfaces*, volume 153. Springer Science & Business Media, 2006.

[38] S. Osher and J.A. Sethian. Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations. *Journal of computational physics*, 79(1):12-49, 1988.

[39] G. Périgaud and R. Saurel. A compressible flow model with capillary effects. *J. Comp. Phys.*, 209:139-178, 2005.

[40] F. Petitpas and S. Le Martelot. A discrete method to treat heat conduction in compressible two-phase flows. *Computational Thermal Sciences: An International Journal*, 6(3), 2014.

[41] F. Petitpas, J. Massoni, R. Saurel, E. Lapebie, and L. Munier. Diffuse interface models for high speed cavitating underwater systems. *Int. J. of Multiphase Flows*, 35(8):747-759, 2009.

[42] F. Petitpas, R. Saurel, E. Franquet, and A. Chinnayya. Modelling detonation waves in condensed energetic materials: Multiphase CJ conditions and multidimensional computations. *Shock waves*, 19(5):377-401, 2009.

[43] M. Pilch and C.A. Erdman. Use of breakup time data and velocity history data to predict the maximum size of stable fragments for acceleration-induced breakup of a liquid drop. *International Journal of Multiphase Flow*, 13(6):741-757, 1987.
[44] J.J. Quirk and S. Karni. On the dynamics of a shock-bubble interaction. *Journal of Fluid Mechanics*, 318:129-163, 1996.

[45] R. Saurel and R. Abgrall. A simple method for compressible multifluid flows. *SIAM J. Sci. Comp.*, 21(3):1115-1145, 1999.

[46] R. Saurel, S.L. Gavrilyuk, and F. Renaud. A multiphase model with internal degrees of freedom: Application to shock-bubble interaction. *Journal of Fluid Mechanics*, 495:283-321, 2003.

[47] R. Saurel, F. Petitpas, and R. Abgrall. Modelling phase transition in metastable liquids: application to cavitating and flashing flows. *Journal of Fluid Mechanics*, 607:313-350, 2008.

[48] R. Saurel, F. Petitpas, and R.A. Berry. Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures. *Journal of Computational Physics*, 228(5):1678-1712, 2009.

[49] F.S. Schranner, X.Y. Hu, and N.A. Adams. On the convergence of the weakly compressible sharp-interface method for two-phase flows. *Journal of Computational Physics*, 324:94-114, 2016.

[50] M. Sussman, P. Smereka, and S. Osher. A level set approach for computing solutions to incompressible two-phase flow. *Journal of Computational Physics*, 114(1):146-159, 1994.

[51] G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber, J. Han, S. Nas, and Y-J. Jan. A front-tracking method for the computations of multiphase flow. *Journal of Computational Physics*, 169(2):708-759, 2001.

[52] P. Welch and P. Boyle. New turbines to enable efficient geothermal power plants. *Geothermal Resources Council Transactions*, 33:765-772, 2009.

[53] A. Wierzba and K. Takayama. Experimental investigation on liquid droplet breakup in a gas stream. *Rept. Inst. High Speed Mech. (Tohoku Univ.)*, 53(382):1-99, 1987.

[54] A.B. Wood. A textbook of sound. *G. Bell and Sons LTD, London*, 1930.

[55] J. Yang, T. Kubota, and E.E. Zukoski. A model for characterization of a vortex pair formed by shock passage over a light-gas inhomogeneity. *Journal of Fluid Mechanics*, 258:217-244, 1994.