# Modelling interactions between waves and diffused interfaces

Kevin Schmidmayer, Joris Cazé, Fabien Petitpas, Eric Daniel, Nicolas Favrie

► To cite this version:

Kevin Schmidmayer, Joris Cazé, Fabien Petitpas, Eric Daniel, Nicolas Favrie. Modelling interactions between waves and diffused interfaces. International Journal for Numerical Methods in Fluids, 2023, 95 (2), pp.215-241. 10.1002/fld.5142. hal-03387818v2

HAL Id: hal-03387818

https://hal.science/hal-03387818v2

Submitted on 24 Aug 2022

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
ARTICLE TYPE

# Modelling interactions between waves and diffused interfaces

Kevin Schmidmayer*<fnref n="1" /> | Joris Cazé<fnref n="1" />,<fnref n="2" /> | Fabien Petitpas<fnref n="1" /> | Éric Daniel<fnref n="1" /> | Nicolas Favrie<fnref n="1" />

Correspondence

*Kevin Schmidmayer, Email:
kevin.schmidmayer@univ-amu.fr /
kevin.schmidmayer@gmail.com

Abstract

When simulating multiphase compressible flows using the diffuse-interface methods, the test cases presented in the literature to validate the modellings with regard to interface problems are always textbook cases: interfaces are sharp and the simulations therefore easily converge to the exact solutions. In real problems, it is rather different because the waves encounter moving interfaces which consequently have already undergone the effects of numerical diffusion. Numerical solutions resulting from the interactions of waves with diffused interfaces have never been precisely investigated and for good reasons, the results obtained are extremely dependent on the model used. Precisely, well-posed models present similar and important issues when such an interaction occurs, coming from the appearance of a wave-trapping phenomenon. To circumvent those issues, we propose to use a thermodynamically-consistent pressure-disequilibrium model with finite, instead of infinite, pressure-relaxation rate to overcome the difficulties inherent in the computation of these interactions. Because the original method to solve this model only enables infinite relaxation, we propose a new numerical method allowing infinite as well as finite relaxation rates. Solutions of the new modelling are examined and compared to literature, in particular we propose the study of a shock on a water-air interface, but also for problems of helium-air and water-air shock tubes, spherical and non-spherical bubble collapses.

KEYWORDS:

Diffuse-interface method, multiphase, compressible, interface, interaction, relaxation

1 | INTRODUCTION

Understanding and controlling complex and physically rich flows, such as multiphase compressible flows, which are often unsteady, are of great importance in various fields such as aeronautics, automotive, aerospace, nuclear energy, naval and also medicine. Good examples are cavitating flows such as the ones encountered around naval propellers where cavitation pockets form at the vicinity of the blades and lead to a turbulent bubbly flow in the wake<fnref n="1" />, or in biomedical applications such as in lithotripsy (treatment for kidney stones)<fnref n="2" /> or, recently, histotripsy (non-invasive treatment for cancers)³ where cavitation bubbles violently collapse and interact with biomaterials.

A possible way of investigation is to use benefits of numerical simulations to analyse in details these flows where experiments often lack of insight and of spatial and temporal resolution. There are numerous techniques to model multiphase compressible flows, all having their pros and cons⁴. Herein we are focusing on diffuse-interface methods, where interfaces are not explicitly tracked but allowed to diffuse numerically⁵,⁶. These methods are particularly interesting because they are able to deal with

<footnote n="1">¹Aix Marseille Univ, CNRS, IUSTI,
Marseille, France</footnote>

<footnote n="2">²Centre National d'Etudes Spatiales, Paris,
France</footnote>
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).

In this context, the test cases presented in the literature to validate the modellings with regard to interface problems are always textbook cases: interfaces are sharp and the simulations therefore easily converge to the exact solutions. In real problems (multi-dimensions, multiple interactions between waves), it is rather different because the waves encounter moving interfaces which consequently have already undergone the effects of numerical diffusion. To our knowledge, numerical solutions resulting from the interactions of waves with diffused interfaces have never been precisely investigated and for good reasons, the results obtained are extremely dependent on the model used. Thus herein, a model is selected and a new numerical method is proposed to overcome the difficulties inherent in the computation of these interactions.

We focus here on interface problem modelling such that it is not necessary to consider velocity disequilibrium between phases. Among single-velocity models, the model of Allaire et al.⁷ and Massoni et al.⁸ has been widely used, such as for shock-induced collapses⁹,¹⁰ and droplet atomization¹¹,¹². However, this model, as mentioned by Allaire et al.⁷ themselves, cannot be mathematical proven with an increasing entropy in mixture regions (i.e. where the interface is diffused) under isobaric closure with a general equation of state. An isothermal closure solves this issue but it is not physically acceptable under pressure disequilibrium. Further, this model cannot predict the collapse time and minimum radius of the Rayleigh collapse of a bubble¹³,¹⁴,¹⁵. Another candidate is the thermodynamically-consistent mechanical-equilibrium model of Kapila et al.¹⁶,¹³,¹⁴,¹⁵. This model includes additional non-conservative terms in the volume-fraction equations that reflect the differences in phase compressibilities in mixture regions. These non-conservative terms unfortunately lead to numerical instabilities during strong compression and expansion near the interface¹⁷. Pressure-disequilibrium models¹⁷,¹⁸ can also be used together with an algorithmic relaxation at infinite rate of the phase pressures at each time step to avert the stability issues of preceding single-pressure model. These modellings theoretically and numerically converge to the mechanical-equilibrium model of Kapila et al.¹⁶.

Although these thermodynamically-consistent, mechanical-equilibrium modellings¹⁶,¹⁷,¹⁸ have been utilized for cavitating flows¹⁷,¹⁹, detonating flows²⁰, surface-tension driven flows²¹,²², droplet atomization²³,²⁴, fracture and fragmentation in ductile materials²⁵,²⁶, bubble dynamics or particularly to collapsing bubbles¹⁵,²⁷,²⁸, they all present similar and important issues when a wave interacts with a diffused interface. Indeed, a wave-trapping phenomenon appears that results from a lower mixture sound speed, Wood's²⁹, than the pure phases. One can note that such issues also appear for models with additional thermal equilibrium (not considered herein) since the behaviour of the mixture sound speed is similar to Wood's.

To circumvent those issues, we propose to use the pressure-disequilibrium model using internal energies¹⁷ with finite, instead of infinite, pressure-relaxation rate, which therefore does not inherit the lower mixture sound speed. However, the original method¹⁷ to solve this model works with an energy-correction procedure that only enables infinite relaxation. Hence, we propose a new numerical method acting directly on the energy fluxes and allowing infinite as well as finite relaxation rates. In addition to its use to solve wave-trapping phenomena, such modelling could also be useful in compaction context for porous materials, where the solid and fluid phases can behave in pressure disequilibrium³⁰. This modelling is implemented in ECOGEN³¹, an open-source tool for multiphase, compressible, multiphysics flows, and it can directly work or can be extended to numerics, such as adaptive mesh refinement²³, or to physics, such as viscosity and surface tension²¹,²², respectively. We also note that finite relaxation rates have already been considered in previous work. In particular, Re and Abgrall³² proposed a pressure-based method for weakly compressible two-phase flows under a Baer–Nunziato-type model³³, which also considers velocity disequilibrium. And Chiocchetti and Müller³⁴ proposed an integration technique for the finite relaxation rates for the Baer–Nunziato model and the pressure-disequilibrium model based on total energies¹⁸, but they only considered mixture test cases (no pure phases). Hence, to our knowledge, previous work did not considered the problem of waves travelling through an interface.

The problem of wave interactions in multiphase flow models is first presented in Section 2. We highlight the differences between interactions with sharp and diffused interfaces. We select the pressure-disequilibrium model using internal energies and detail it in Section 3. The numerical method we employ to solve the resulting equations, including the new energy-correction and pressure-relaxation procedures, are outlined in Section 4. Solutions of the new modelling are examined and compared to literature in Section 5 for problems of helium–air and water–air shock tubes, shock on water–air interface, spherical and nonspherical bubble collapses. Note that the expression modelling indicates herein a combination of the mathematical model and the numerical scheme used to solve it. Finally, a discussion on the finite pressure relaxation rate and a detailed conclusion are given in Sections 6 and 7, respectively.
<figure><img src="image_2.png" /><figcaption>FIGURE 1 Problem configuration for the problem of a shock on a water-air interface.</figcaption></figure>

## 2 | PROBLEM: INTERACTIONS BETWEEN WAVES AND DIFFUSED INTERFACES

In order to assess the interactions between waves and diffused interfaces from well-posed models, we present a simplified test where a wave, here a weak shock propagating in water, interacts with a water-air interface (high density jumps). We compare the results when the water-air interface is initially sharp with those obtained when it is artificially diffused. This last situation is representative of a wave interacting with an interface in motion and occurs many times in multi-dimensional simulations. To do so, a shock tube of 1 m long is presented with initially three states (Figure 1). A shock at $x_s = 0.3$ m is separating water in a shocked and non-shocked states. The water-air interface is positioned at $x_i = 0.7$ m. Precise information on this test case can be found in the result section 5.3.

In both cases (initially sharp or diffused interface), simulations are performed with the modellings of Saurel et al.¹⁷ and of Pelanti & Shyue¹⁸ corresponding to two different formulations of a well-posed pressure-disequilibrium model. These two formulations are solved in the context of infinite pressure-relaxation rate. Results are compared with the exact solution in Figure 2. We observe that both modellings present excellent results when the shock interacts with the sharp interface. However, when the interface is diffused, the solutions obtained are seriously polluted by the interaction. Indeed, for these modellings, the speed of sound in the artificial mixture (corresponding to the diffused interface location) follows the Wood speed of sound²⁹,³⁵

$$ \frac{1}{\rho c^2} = \sum_{k=1}^{N} \frac{\alpha_k}{\rho_k c_k^2}, \qquad (1) $$

which presents significantly lower speeds than the ones of the pure phases. Figure 3 shows its behaviour in function of the volume fraction for a standard water-air mixture. These lower speeds affect significantly the wave propagation through the interface and therefore induces a complete reshaping of the waves, here the transmitted and reflected shocks, and with particular transitions between the different states.

One should note that such behaviour is not so pronounced if there is not a clear discontinuity of bulk modulus between the phases, e.g. helium and air. As a matter of fact, the Wood sound speed does not manifest lower speeds than the ones of the pure phases in this framework.

While this wave-interface interaction demonstrates the need of improvement regarding the modelling in order to circumvent problems shown in the mixture region, the literature flourish of interface sharpening techniques, such as³⁶,¹³,³⁷ to name a few. Having an interface as sharp as possible only partially alleviates the problem because these methods only show more disturbance of the flow variables near the interface for steeper interfaces and inevitably diffusion still occurs. It is also important to note that these methods are often non-conservatives which poses other issues. In addition, high-order methods, such as WENO³⁸,³⁹,⁴⁰,⁹,¹³,¹⁵,⁴¹ or ADER⁴²,²², require initial smearing of the interfaces to avoid stability issues. In conclusion, act on the sharpness of the interface is not the solution. And because all the well-posed models show similar behaviour, we instead propose to act on the pressure relaxation of the pressure-disequilibrium model. This leads us to propose a new method allowing finite relaxation rates while guaranteeing conservation.
<figure><img src="image_3.png" /><figcaption>FIGURE 2 Solution of the shock on an initially sharp (left) and diffused (right) water-air interface for the Saurel et al.¹⁷ and Pelanti & Shyue¹⁸ modellings (infinite relaxations). The exact solution is also plotted for comparison. (b) and (d) are magnified views of (a) and (c), respectively.</figcaption></figure>

## 3 | MULTIPHASE FLOW MODEL

Compressible single-velocity multiphase flow models can be written in a general form as

$$ \frac{\partial \mathbf{q}}{\partial t} + \nabla \cdot \mathbf{F}(\mathbf{q}) + \mathbf{h}(\mathbf{q}) \nabla \cdot \mathbf{u} + \mathbf{g}(\mathbf{q}) \nabla \cdot \mathbf{m}(\mathbf{q}) = \mathbf{r}(\mathbf{q}), \qquad (2) $$

where $\mathbf{q}$ is the state vector, $\mathbf{F}$ is the flux tensor, $\mathbf{u}$ is the mixture velocity field, and $\mathbf{h}, \mathbf{g}, \mathbf{m}$ and $\mathbf{r}$ are non-conservative quantities we describe subsequently. We recall that the models mentioned herein are in velocity equilibrium and formally ensured the principles of conservation of mass, momentum, and total energy. They also respect the second law of thermodynamics and are hyperbolic with eigenvalues either equal to $u$ or $u \pm c$, where $u$ is the velocity in the $x$-direction and $c$ is the mixture speed of sound. We select the pressure-disequilibrium model based on internal energies and we present it below. Whereas, for completeness, the mechanical-equilibrium model of Kapila et al.¹⁶ and the pressure-disequilibrium model based on total energies¹⁸,⁴³, which we call UEqTotE, are recalled in A and B, respectively.
<figure><img src="image_4.png" /><figcaption>FIGURE 3 Wood and frozen speeds of sound for standard water–air mixture.</figcaption></figure>

The pressure-disequilibrium model based on internal energies is expressed for N phases as

$$ \mathbf{q} = \begin{bmatrix} \alpha_k \\ \alpha_k \rho_k \\ \rho \mathbf{u} \\ \alpha_k \rho_k e_k \end{bmatrix}, \quad \mathbf{F} = \begin{bmatrix} \alpha_k \mathbf{u} \\ \alpha_k \rho_k \mathbf{u} \\ \rho \mathbf{u} \otimes \mathbf{u} + p \mathbf{I} \\ \alpha_k \rho_k e_k \mathbf{u} \end{bmatrix}, \quad \mathbf{h} = \begin{bmatrix} -\alpha_k \\ 0 \\ \mathbf{0} \\ \alpha_k p_k \end{bmatrix}, \quad \mathbf{r} = \begin{bmatrix} \delta p_k \\ 0 \\ \mathbf{0} \\ -p_I \delta p_k \end{bmatrix}, \qquad (3) $$

$\mathbf{g} = \mathbf{m} = [0, 0, 0, 0]^T$, where $\alpha_k$, $\rho_k$, $p_k$ and $e_k$ are the volume fraction, density, pressure and internal energy of each phase, respectively, and for which k indicates the phase index. The mixture density and pressure are

$$ \rho = \sum_{k=1}^{N} \alpha_k \rho_k \quad \text{and} \quad p = \sum_{k=1}^{N} \alpha_k p_k, \qquad (4) $$

while the mixture total energy is

$$ E = e + \frac{1}{2} \|\mathbf{u}\|^2, \qquad (5) $$

where e is the mixture specific internal energy

$$ e = \sum_{k=1}^{N} Y_k e_k (\rho_k, p_k). \qquad (6) $$

In (6), $e_k(\rho_k, p_k)$ is defined via an equation of state (EOS) and $Y_k$ are the mass fractions

$$ Y_k = \frac{\alpha_k \rho_k}{\rho}. \qquad (7) $$

Herein, we will consider two-phase mixtures of gas (g) and liquid (l), for which the gas is modeled by the ideal-gas EOS

$$ p_g = \rho_g(\gamma_g - 1)(e_g - e_{g,ref}), \qquad (8) $$

and the liquid is modeled by the stiffened-gas (SG) EOS

$$ p_l = \rho_l(\gamma_l - 1)(e_l - e_{l,ref}) - \gamma_l \pi_{\infty}, \qquad (9) $$

where $\gamma$, $e_{ref}$ and $\pi_{\infty}$ are model parameters$^{44}$. Other EOS, e.g. NASG$^{45}$ or Mie–Grüneisen, could easily be substituted in our framework. r represents the relaxation of pressures between the phases with

$$ \delta p_k = \sum_{j \neq k}^{N} \mu_{k,j} (p_k - p_j), \qquad (10) $$
where j are phases different from k and $\mu_{k,j}$ are the pressure-relaxation coefficients related to the k-j interactions. Herein, the pressure-relaxation coefficient $\mu$ is considered the same for each phase combination. The interfacial pressure is defined as

$$p_I = \frac{\sum_k^N (p_k \sum_{j \neq k}^N z_j)}{(N-1) \sum_k^N z_k}, \qquad (11)$$

where $z_k = \rho_k c_k$ and $c_k$ are the acoustic impedance and speed of sound of the phase k, respectively.

Since pressures are in disequilibrium here, the total energy equation of the mixture is replaced by the internal-energy equation for each phase. Nevertheless, conservation of the mixture total energy can be written in its usual form

$$\frac{\partial \rho E}{\partial t} + \nabla \cdot [(\rho E + p) \mathbf{u}] = 0. \qquad (12)$$

We note that (12) is redundant when the internal energy equations are also computed. However, in practice, we include it in our computations (detailed in Section 4.3) to ensure that the total energy is numerically conserved, and thus preserve a correct treatment of shock waves.

The mixture speed of sound, also called frozen speed of sound, is defined according to

$$c^2 = \sum_{k=1}^{N} Y_k c_k^2. \qquad (13)$$

Its behaviour is also shown in Figure 3 in function of the volume fraction for a standard water-air mixture. Note that after applying an infinite pressure relaxation, the model converge to the mechanical-equilibrium model of Kapila et al.¹⁶ (A1) and the effective mixture speed of sound matches Wood's (1).

# 4 | NUMERICAL METHODS

In this section is presented the numerical method including the general scheme, the new total-energy correction scheme and the pressure-relaxation procedures for infinite and finite rates. Note that the method, except the new correction scheme, can also be applied to solve the UEqTotE model (B3) of Pelanti & Shyue¹⁸.

## 4.1 | General scheme

We solve (2) numerically using a splitting procedure between the left-hand-side terms associated with the flow and the right-hand-side terms associated with our relaxation procedure. First, the time evolution of $\mathbf{q}$ on a computational cell $i$ with volume $V_i$ and surface $A$ with normal unit vector $\mathbf{n}$ is given by the explicit finite-volume Godunov⁴⁶ scheme

$$\mathbf{q}_i^{n+1} = \mathbf{q}_i^n - \frac{\Delta t}{V_i} \left( \sum_{s=1}^{N_s} A_s \mathbf{F}_s^* \cdot \mathbf{n}_s + \mathbf{h}(\mathbf{q}_i^n) \sum_{s=1}^{N_s} A_s \mathbf{u}_s^* \cdot \mathbf{n}_s + \mathbf{g}(\mathbf{q}_i^n) \sum_{s=1}^{N_s} A_s \mathbf{m}_s^* (\mathbf{u}_s^*) \cdot \mathbf{n}_s \right), \qquad (14)$$

where $n$ is the time-step index. The relaxation terms are then solved using one of the procedures detailed in 4.4 to complete the time-step integration. We also utilize MUSCL spatial reconstructions of the primitive state variables presented in 4.2. We note that reconstructing the conservative variables instead leads to spurious oscillations near material interfaces⁹, and using a characteristic-based reconstruction in our implementation significantly increases computational costs but does not improve results. At the volume-volume interfaces, the associated Riemann problem is computed using the HLLC approximate solver¹⁷,⁴⁷, giving the flux tensor $\mathbf{F}_s^*$, the flow-velocity vector $\mathbf{u}_s^*$ and non-conservative vector $\mathbf{m}_s^*$. Note that in the case of the volume-fraction equation, $\mathbf{u}_s^*$ is depreciated to the speed of the contact discontinuity $\mathbf{s}_s^*$ to ensure a correct treatment of this transport equation (here without the relaxation terms). The solution of (14) is restricted by the usual CFL criterion.
## 4.2 | MUSCL scheme

We use the second-order-accurate MUSCL scheme of Schmidmayer et al. 31 with two-step time integration

$$ q_i^{n+\frac{1}{2}} = q_i^n + \frac{1}{2} \Delta t \mathbf{L}(q_i^n), \qquad (15) $$

$$ q_i^{n+1} = q_i^n + \Delta t \mathbf{L} \left( q_i^{n+\frac{1}{2}} \right), \qquad (16) $$

where the operator L is the numerically approximated fluxes and non-conservative terms, function of the state vector q at different time stages. The first step is a prediction and the usual piece-wise linear MUSCL reconstruction 47 is used on the primitive variables. The Minmod 47,48, van Leer 49, monotonized central (MC) 50 and THINC (only for interface) 37 slope limiters are employed to minimize diffusion and are specified for each test case. The corrections based on total energy and relaxations, detailed below, are applied at each stage.

## 4.3 | Total-energy correction schemes

To ensure the conservation of total energy, a procedure correcting the non-conservative terms of the pressure-disequilibrium model based on internal energies (3) is required and it uses the mixture total-energy relation (12).

The method proposed in Saurel et al. 17 corrects the total energy after the relaxation procedure and imposes an infinite relaxation. In the following, we call this method PUEq.

Herein, we propose a new method correcting the total energy before the relaxation procedure and therefore allowing finite or infinite relaxations. This procedure takes place during the flux computation of the hyperbolic step. In the following, we call this new method UEq.

Let us first express the discretized equation on internal energy of each phase

$$ (\alpha_k \rho_k e_k)_i^{n+1} = (\alpha_k \rho_k e_k)_i^n + F_{\alpha_k \rho_k e_k}^*, \qquad (17) $$

where

$$ F_{\alpha_k \rho_k e_k}^* = -\frac{\Delta t}{V_i} \left( \sum_{s=1}^{N_s} A_s (\alpha_k \rho_k e_k \mathbf{u})_s^* \cdot \mathbf{n}_s + (\alpha_k \rho_k)^n \sum_{s=1}^{N_s} A_s \mathbf{u}_s^* \cdot \mathbf{n}_s \right). \qquad (18) $$

Hence, the sum of internal energies can be expressed as

$$ \sum_{k=1}^{N} (\alpha_k \rho_k e_k)_i^{n+1} = \sum_{k=1}^{N} (\alpha_k \rho_k e_k)_i^n + F_{\rho e}^*, \qquad (19) $$

where

$$ F_{\rho e}^* = \sum_{k=1}^{N} F_{\alpha_k \rho_k e_k}^*. \qquad (20) $$

To ensure conservation of total energy, we must satisfy that the sum of internal energies equals the mixture internal energy ρe obtained from the mixture total-energy equation (12) and momentum (3)

$$ \sum_{k=1}^{N} (\alpha_k \rho_k e_k)_i^n = (\rho e)_i^n, \quad \sum_{k=1}^{N} (\alpha_k \rho_k e_k)_i^{n+1} = (\rho e)_i^{n+1} \qquad (21) $$

with

$$ \rho e = \rho \left( E - \frac{1}{2} \| \mathbf{u} \|^2 \right). \qquad (22) $$

Defining the variation (flux) of mixture internal energy as

$$ \Delta \rho e = (\rho e)_i^{n+1} - (\rho e)_i^n, \qquad (23) $$

one can write

$$ \sum_{k=1}^{N} (\alpha_k \rho_k e_k)_i^{n+1} = \sum_{k=1}^{N} (\alpha_k \rho_k e_k)_i^n + \Delta \rho e, \qquad (24) $$
and therefore one must satisfy

$$F_{pe}^{*} = \Delta \rho e. \qquad (25)$$

To do so, we propose the following correction scheme

$$(\alpha_k \rho_k e_k)_i^{n+1} = (\alpha_k \rho_k e_k)_i^n + F_{\alpha_k \rho_k e_k}^* + \alpha_k^n \epsilon, \qquad (26)$$

where

$$\epsilon = \Delta \rho e - \sum_{k=1}^{N} F_{\alpha_k \rho_k e_k}^{*}, \qquad (27)$$

and where the mixture internal-energy variation can be computed by

$$\Delta \rho e = F_{\rho E}^{*} - \frac{1}{2} \left( \frac{\| (\rho \mathbf{u})_i^{n+1} \|^2}{\rho_i^{n+1}} - \frac{\| (\rho \mathbf{u})_i^n \|^2}{\rho_i^n} \right), \qquad (28)$$

with

$$F_{\rho E}^{*} = -\frac{\Delta t}{V_i} \sum_{s=1}^{N_s} A_s [(\rho E + p) \mathbf{u}_s^{*} \cdot \mathbf{n}_s]. \qquad (29)$$

This assures the conservation of total energy

$$\sum_{k=1}^{N} (\alpha_k \rho_k e_k)_i^{n+1} = \sum_{k=1}^{N} (\alpha_k \rho_k e_k)_i^n + \sum_{k=1}^{N} (F_{\alpha_k \rho_k e_k}^* + \alpha_k^n \epsilon), \qquad (30)$$

$$\sum_{k=1}^{N} (\alpha_k \rho_k e_k)_i^{n+1} = \sum_{k=1}^{N} (\alpha_k \rho_k e_k)_i^n + \Delta \rho e. \qquad (31)$$

Note that one can use the mass fractions $Y_k^n$ instead of the volume fractions $\alpha_k^n$ in (26) and that we experienced very similar results.

## 4.4 | Pressure-relaxation procedures

The pressure-disequilibrium model (3) can be solved with either a finite or an infinite relaxation procedure for $\mathbf{r}$. At each time step we solve the non-relaxed, hyperbolic equations ($\mu \to 0$) using (14), then we solve the system of ordinary differential equations (ODE)

$$\frac{\partial \mathbf{q}}{\partial t} = \mathbf{r}(\mathbf{q}), \qquad (32)$$

that relaxes the disequilibrium pressures for a given $\mu$ or $\mu \to +\infty$. When multi-stage time integration is used, these procedures are performed at each stage.

### 4.4.1 | Infinite relaxation

In the case of the infinite relaxation (stiff pressure relaxation), $\mu \to +\infty$, the model converges to a single, equilibrium pressure. We use the infinite-relaxation procedure of Saurel et al.¹⁷.

Specifically, after manipulations of the internal-energy equations and the mixture mass equation, the mixture energy equation can be written for two phases as

$$e - e^0 + (\hat{p}_{I,1} - \hat{p}_{I,2}) (Y_1 v_1 - Y_1 v_1^0) = 0, \qquad (33)$$

where superscript 0 indicates the hyperbolic step index, $v_k$ are the specific volumes of each state and $\hat{p}_{I,k} = \frac{1}{v_k - v_k^0} \int_0^{\Delta t} p_I \frac{\partial v_k}{\partial t} dt$. In order that the mixture energy conservation be fulfilled it is necessary that $\hat{p}_{I,1} = \hat{p}_{I,2} = \hat{p}_I$. A possible estimate compatible with the entropy inequality is $\hat{p}_I = p$, the relaxed pressure. This leads to solve the following system of equations for N phases

$$e_k (p, v_k) - e_k^0 (p_k^0, v_k^0) + p (v_k - v_k^0) = 0, \qquad (34)$$
which involves $N + 1$ unknowns, $v_k$ and $p$. Its closure is achieved using the saturation constraint ($\sum_k \alpha_k = 1$) and this system can be replaced by a single equation to solve with a single unknown ($p$)

$$ \sum_k (\alpha \rho)_k v_k (p) = 1, \qquad (35) $$

where $(\alpha \rho)_k$ are constant during the relaxation process and $v_k(p)$ are determined with the help of the EOS. For the general case of SG EOS, the energy equations become

$$ v_k(p) = \frac{p_k^0 + \gamma_k \pi_{\infty,k} + p(\gamma_k - 1)}{\gamma_k(p + \pi_{\infty,k})} v_k^0. \qquad (36) $$

We ultimately solve (35) using the Newton-Raphson method to find the relaxed pressure. Then the phase densities and volume fractions are determined.

We recall that after applying the infinite pressure-relaxation procedure, the model converges to the mechanical-equilibrium model of Kapila et al.¹⁶ (A1) and the effective mixture speed of sound matches (1).

### 4.4.2 | Finite relaxation

For a given and finite $\mu$, system (32) of ODE can be generalized and simplified by using the pressure equations instead of the energy equations

$$ \frac{\partial p_k}{\partial t} = - \frac{\rho_k c_{I,k}^2}{\alpha_k} \mu \Delta p_k, \qquad (37) $$

where

$$ c_{I,k}^2 = \frac{\left( \frac{p_I}{\rho_k^2} - \frac{\partial e_k}{\partial \rho_k} \right)_{p_k}}{\left( \frac{\partial e_k}{\partial p_k} \right)_{\rho_k}} \quad \text{and} \quad \Delta p_k = \sum_{j \neq k}^{N} p_k - p_j. \qquad (38) $$

$c_{I,k}^2$ is expressed for ideal and SG EOS in C. In addition, the following thermodynamical constraints have to be verified at each time step:

* Volume fractions conserved between [0, 1].

* Pressures can only converge to the same value and cannot cross each other, i.e. the relaxed pressure of the initially (before relaxation) higher-pressure fluid cannot be inferior to the relaxed pressure of the lower-pressure fluid, and inversely.

This system of ODE can be solved using your favourite ODE solver, such as LSODA⁵¹,⁵². However, even though these types of solvers are very good from a mathematical point of view, they do not take into account thermodynamical constraints. Therefore, for a more controlled and faster procedure, we propose a classic, first-order, explicit, Euler scheme with time-step subdivisions. The number of subdivisions is adapted at each time step to verify the constraints.

#### Scheme

The scheme to integrate the relaxation terms is as follow

$$ \begin{align} \alpha_k^{n+1} &= \alpha_k^n + \Delta t_l \mu \Delta p_k^n, \notag \\ \rho_k^{n+1} &= \frac{\alpha_k^n \rho_k^n}{\alpha_k^{n+1}}, \tag{39} \\ p_k^{n+1} &= p_k^n - \Delta t_l \frac{\rho_k^n (c_{I,k}^n)^2}{\alpha_k^n} \mu \Delta p_k^n, \notag \end{align} $$

where $\Delta t_l$ is the local, Euler, subdivided time step respecting the constraints detailed hereafter and initialized by the remaining time necessary to complete the global time step $\Delta t$.
## Constraint on volume fractions

The volume fractions have to be conserved between [0, 1]. Therefore, the following algorithm can be applied for each phase to determined $\Delta t_l$

$$ 
\begin{array}{l}
\text{if } (\mu \Delta p_k^n > \varepsilon) \quad \Delta t_{\max} = \min \left( \Delta t_{\max}, \frac{1 - \alpha_k^n}{\mu \Delta p_k^n} \right), \\
\text{else if } (\mu \Delta p_k^n < -\varepsilon) \quad \Delta t_{\max} = \min \left( \Delta t_{\max}, \frac{-\alpha_k^n}{\mu \Delta p_k^n} \right), \\
\Delta t_{\max} = 0.5 \Delta t_{\max}, \\
\text{if } (\Delta t_{\max} < \Delta t_l) \quad \Delta t_l = \Delta t_{\max},
\end{array}
\qquad (40) $$

where $\varepsilon$ is a small value (typically $10^{-10}$) and $\Delta t_{\max}$ is initialized with a significant value.

## Constraint on pressures

The condition is such that the pressures can only relax at a maximum where they are all equals $p_k^{n+1} = p^{n+1}$. For two phases in interaction (subscripted as 0 and 1), one can therefore determine the maximal time step:

$$ \frac{\partial p_1}{\partial t} = - \frac{\rho_1 c_{I,1}^2}{\alpha_1} \mu (p_1 - p_2) \quad \text{and} \quad \frac{\partial p_2}{\partial t} = \frac{\rho_2 c_{I,2}^2}{\alpha_2} \mu (p_1 - p_2), \qquad (41) $$

hence

$$ \frac{\partial p_1 - p_2}{\partial t} = -\mu (p_1 - p_2) \rho \bar{c}^2, \qquad (42) $$

where

$$ \rho \bar{c}^2 = \frac{\rho_1 c_{I,1}^2}{\alpha_1} + \frac{\rho_2 c_{I,2}^2}{\alpha_2} \qquad (43) $$

is the mixture bulk modulus. Once discretized, we obtain

$$ (p_1 - p_2)^{n+1} - (p_1 - p_2)^n = -\Delta t_{\max} \mu (p_1 - p_2)^n (\rho \bar{c}^2)^n. \qquad (44) $$

Finally, with $(p_1 - p_2)^{n+1} = 0$, the following algorithm can be applied for each combination of two phases to determined $\Delta t_l$

$$ 
\begin{array}{l}
\Delta t_{\max} = \frac{1}{\mu (\rho \bar{c}^2)^n}, \\
\Delta t_{\max} = 0.5 \Delta t_{\max}, \\
\text{if } (\Delta t_{\max} < \Delta t_l) \quad \Delta t_l = \Delta t_{\max}.
\end{array}
\qquad (45) $$

## Case of a unique pressure

During this procedure, if the pressures are completely relaxed, i.e. a unique pressure for all phases, we terminate the Euler scheme and we perform from the initial state an infinite relaxation to guarantee a unique pressure and better estimate the solution. This also assures a faster computation.

## 4.5 | Modelling summary

Finally, the steps involved in the first-order method can be summarized as

1. Hyperbolic step using (14) combined with the total-energy correction scheme (26) for the internal-energy equations.

2. Infinite- or finite-relaxation procedure (detailed in 4.4).

We recall that the total-energy correction scheme and relaxation procedures are applied at each stage for a higher-order method.

# 5 | RESULTS

In the following, the new numerical method, UEq, is highlighted in bold to distinguish it from existing methods we compare against: Saurel et al.¹⁷ and Pelanti & Shyue¹⁸,⁴³. We remind that we call the latter PUEq and UEqTotE, respectively. The
<figure><img src="image_5.png" /><figcaption>**FIGURE 4** Problem configuration for the helium–air and water–air shock-tube problem.</figcaption></figure>

modellings are tested in the limit of infinite and null pressure relaxations, even if the latter has no physical significance. Indeed, an infinite relaxation rate is physically accurate when the relaxation time to obtain equilibrium is smaller than the time-step. While a null relaxation rate means there is not any equilibrium reached after an infinite time.

Before getting into the detailed test cases, we note that in the context of finite relaxation, we experienced through our tests different results when computing the **UEq** model (3) with or without the non-conservative terms $\alpha_k p_k \nabla \cdot \mathbf{u}$ of the internal-energy equations. Details and comparisons for interface problems are given in the following and we chose the notation where “**UEq**” stands for the **UEq** model without computing the non-conservative terms and “**UEq NC**” stands for the **UEq** model computing those terms. Note that we observed extremely minor differences when infinite relaxation was computed and both configurations are therefore not presented. Moreover, we recall that the scheme is still fully conservative with regard to the energy, whether we compute or not those non-conservative terms. In addition, note for the UEqTotE modelling that we experienced valid results only if all the non-conservative terms are computed, whatever the relaxation process.

The EOS used in the following test cases are summarized in Table 1. Note that two sets of parameters are given for water as they are respectively more accurate, with regard to real fluid behaviour, in the flow range of the corresponding test cases.

<table><caption>**TABLE 1** EOS of the different fluids used in the different test cases.</caption><thead><tr><th>Fluid</th><th>EOS type</th><th>γ [-]</th><th>π<sub>∞</sub> [Pa]</th><th>e<sub>ref</sub> [J/kg]</th></tr></thead><tbody><tr><td>Air</td><td>IG</td><td>1.4</td><td>—</td><td>0</td></tr><tr><td>Helium</td><td>IG</td><td>1.666667</td><td>—</td><td>0</td></tr><tr><td>Water1 (shock-tube tests)</td><td>SG</td><td>4.4</td><td>6x10<sup>8</sup></td><td>0</td></tr><tr><td>Water2 (bubble tests)</td><td>SG</td><td>2.35</td><td>1x10<sup>9</sup></td><td>0</td></tr></tbody></table>

Tests cases on transport are not presented since we did not observe any fluctuations (excellent results) from the different modellings.

The first test case presented is a helium–air shock tube and is typical for testing numerical methods. We are going to show that no matter what numerics we use this test case always goes well. On the other hand, if a significantly different phase bulk modulus is encountered (true two-phase interface), such as for water–air interfaces presented hereafter, numerics matter.

## 5.1 | Helium–Air shock tube

This test case aims to present the behaviour of the modellings in shock-tube configurations with low density jumps. To do so, a helium–air shock tube of 1 m long is presented with the initial discontinuity (diaphragm) between the two fluids positioned at $x = 0.7$ m (Figure 4). Helium at $p_{He} = 10^6$ Pa and $\rho_{He} = 0.5$ kg/m$^3$ is initially present in the left region, while air at $p_{air} = 10^5$ Pa and $\rho_{air} = 1$ kg/m$^3$ is present in the right region. Initial velocities are nulls. Non-reflecting boundary conditions, the van Leer and THINC limiters and $10^3$ cells are used.

Figure 5 shows the ability of the different modellings to produce very reliable solutions for low density jumps when infinite relaxation is used. Without any relaxation, as shown in Figure 6, the modellings also produce reliable solutions though with a
<figure><img src="image_6.png" /><figcaption>**FIGURE 5** Helium–Air shock-tube solution at *t* = 0.1 ms for the proposed **UEq** modelling as well as the PUEq and UEqTotE modellings with infinite relaxation. The exact solution is also plotted as reference.</figcaption></figure>

relatively small pressure peak at the interface location. We also note that **UEq** and UEqTotE give slightly better results than the **UEq NC** modelling.

In addition, we recall that a null pressure relaxation rate, even if it has no physical significance, is a limit test case and since the shock-tube problem is a problem between pure phases, relaxation between the phases should theoretically not even be considered.

## 5.2 | Water–Air shock tube

This test case aims to present the behaviour of the modellings in shock-tube configurations with high density jumps. To do so, a water–air shock tube of 1 m long is presented with the initial discontinuity (diaphragm) between the two fluids positioned at *x* = 0.7 m (Figure 4). Water at *p*<sub>water</sub> = 10<sup>9</sup> Pa and *ρ*<sub>water</sub> = 10<sup>3</sup> kg/m<sup>3</sup> is initially present in the left region, while air at *p*<sub>air</sub> = 10<sup>5</sup> Pa and *ρ*<sub>air</sub> = 50 kg/m<sup>3</sup> is present in the right region. Initial velocities are nulls. Non-reflecting boundary conditions, the van Leer and THINC limiters and 10<sup>3</sup> cells are used.

Figure 7 shows the ability of the different modellings to produce very reliable solutions for high density jumps when infinite relaxation is used. Though one should note a relatively small fluctuation on pressure and velocity near the tail of the expansion waves (typically observed in the literature).
<figure><img src="image_7.png" /><figcaption>FIGURE 6 Helium–Air shock-tube solution at t = 0.1 ms for the proposed UEq and UEq NC modellings as well as the UEqTotE modelling without relaxation (μ = 0). The exact solution is also plotted as reference. Since pressures are in disequilibrium when the relaxation rate is not infinite, the mixture pressure p = ΣNk=1 αkpk is presented.</figcaption></figure>

Figure 8 shows the solutions given by the different modellings without relaxation. We observe that the **UEq** modelling presents good results with a small density undershot near the interface location. While the **UEq NC** and UEqTotE modellings, which are mathematically identical (though not the same computed models and numerical methods), present very close but wrong results. Indeed, only the expansion waves and interface location seem respected while the shocked states and shock position are incorrect and present non-expected fluctuations. This behaviour is explained by the values taken from the non-conservative terms α_k p_k ∇·**u** of the internal-energy equations. Unlike the previous test case with helium and air, those terms can be about three orders of magnitude higher in the current configuration (higher pressure). Even though the scheme ensures mixture energy conservation, those values produce an inadequate distribution of the energy between the phases. One should note that the pressure relaxation acts as a recalibration, and therefore as a correction, of the energy distribution. This explains why **UEq** and **UEq NC** give almost identical results when infinite relaxation is performed.

We conclude from this type of problem that there is no need to relax with our modelling to obtain good results, whereas it is necessary for the other modellings to correct the energy distribution produced by the schemes. Based on this conclusion, the results presented in Figure 2 and the fact that we want to do finite relaxations, the **UEq NC** and UEqTotE modellings are therefore no longer presented in what follows.

## 5.3 | Shock on water-air interface

This test case is the one presented in Section 2 and it aims to present the behaviour of the modellings when a wave, here a weak shock, interacts with a sharp or diffuse interface (high density jump). We recall that to do so, a shock tube of 1 m long is presented with initially three states (Figure 1). A shock at $x_s = 0.3$ m is separating water in a shocked and non-shocked states. The water-air interface is positioned at $x_i = 0.7$ m. The initial conditions are given in Table 2. Note that this shocked state can be obtained from an initial shock tube with a high-pressure state of $p = 10^7$ Pa, $\rho = 10^3$ kg/m$^3$ and $u = 0$ m/s. Non-reflecting boundary conditions, the van Leer and THINC limiters and $10^3$ cells are used.

Three relaxation rates are investigated. The first one is infinite, to match with the PUEq modelling. The second is null ($\mu = 0$), that is to say no pressure relaxation. The third is finite with a value to be defined using the following approximations.

Similarly to what is done by Petitpas et al.²⁰ in detonation context, we propose an approximation of Eq. (42) as

$$ \frac{1}{\tau} = \mu \rho \bar{c}^2, \qquad (46) $$
<figure><img src="image_8.png" /><figcaption>FIGURE 7 Water–Air shock-tube solution at $t = 241$ µs for the proposed UEq modelling as well as the PUEq and UEqTotE modellings with infinite relaxation. The exact solution is also plotted as reference. The Water1 SG EOS is used for water.</figcaption></figure>

<table><caption>TABLE 2 Initial conditions of the three states of the problem of shock on water–air interface.</caption><thead><tr><th>State position</th><th>Fluid</th><th>p [Pa]</th><th>ρ [kg/m³]</th><th>u [m/s]</th></tr></thead><tbody><tr><td>Left</td><td>Shocked water</td><td>5.0421806x10<sup>6</sup></td><td>1001.8658</td><td>3.0337923</td></tr><tr><td>Middle</td><td>Water</td><td>10<sup>5</sup></td><td>1000</td><td>0</td></tr><tr><td>Right</td><td>Air</td><td>10<sup>5</sup></td><td>1</td><td>0</td></tr></tbody></table>

allowing to link the pressure relaxation coefficient $\mu$ with the mixture bulk modulus (43) and a characteristic time $\tau$. Since we are considering the case of a wave travelling through an interface, we propose to define the characteristic time as the time for a wave to cross entirely the interface, leading to

$$ \tau = \frac{D_I}{\tilde{c}}, \qquad (47) $$

where $D_I$ is the interface thickness and $\tilde{c}$ is the average speed of sound we are expecting.
<figure><img src="image_9.png" /><figcaption>FIGURE 8 Water–Air shock-tube solution at t = 241 µs for the proposed UEq and UEq NC modellings as well as the UEqTotE modelling without relaxation. The exact solution is also plotted as reference. The Water1 SG EOS is used for water.</figcaption></figure>

In the context of this test case with the diffused-interface configuration (detailed below), we estimate $D_I \approx 0.02$ m and $\tilde{c} \approx 10^3$ m/s, leading to $\tau \approx 2 \times 10^{-5}$ s. Considering $\alpha = 0.5$, the mixture bulk modulus can be estimated to $\rho\bar{c}^2 \approx 5 \times 10^9$ Pa, finally leading to $\mu \approx 10^{-5}$ Pa$^{-1}$s$^{-1}$, which is taken as finite value for this test case.

### 5.3.1 | Sharp-interface configuration

Figure 9 shows the behaviour of the interaction of the shock wave with an initially sharp interface ($\alpha = 0$ or 1). The UEq modelling is presented with infinite and finite relaxations. We observe that all the solutions are in very good agreement with the exact solution. Even though there is a relatively small pressure peak at the interface location when a finite relaxation is used.

### 5.3.2 | Diffused-interface configuration

To initially have a diffused interface, we employ an initial interface smearing procedure involving smearing the volume fraction across the interface using an hyperbolic tangent function

$$ \alpha_{\text{water}} = \frac{1}{2} \left[ 1 - \tanh \left( \frac{x - x_i}{h} \right) \right], \qquad (48) $$

where $h$ is the characteristic length of the smearing, here $h = 3$ mm. The conservative variables then follow from simple mixture relations, allowing thermodynamic consistency. The smeared interface is observable on Figure 10 (a).

Figure 10 shows the behaviour of the interaction of the shock wave with this initially diffused interface and this for different relaxation rates.

First, as observed in Section 2, the solutions obtained with an infinite relaxation (PUEq and UEq) are seriously polluted by the interaction of the shock with the diffused interface. Indeed, we recall that for these modellings the speed of sound in a mixture (here interface) follows the Wood speed of sound which presents significantly lower speeds than the ones of the pure phases (see Figure 3). This affects the wave propagation through the interface and therefore induces a complete reshaping of the waves, here the transmitted and reflected shocks, and with particular transitions between the different states.

Second, we observe a relatively good agreement of the UEq modelling without any relaxation. Even though relatively small perturbations are noticeable. Furthermore, one can note that the transmitted shock wave is in advance compared to its theoretical location. Indeed, when $\mu = 0$, the speed of sound in a mixture (interface) now follows the frozen speed of sound which is almost equal to the speed of sound of the higher-density fluid in almost all the range of the volume fraction (see Figure 3). Therefore, the shock travels faster when crossing the interface: it travels almost like it is travelling only in water and not in air for most of the interface thickness.
<figure><img src="image_10.png" /><figcaption>**FIGURE 9** Solution of the shock on an initially sharp water-air interface at $t = 600$ µs for the proposed UEq modelling with infinite or finite relaxations. The exact and PUEq solutions are also plotted for comparison. The Water1 SG EOS is used for water. (d) is a magnified view of (c).</figcaption></figure>

Third, when picking a finite relaxation rate of $\mu = 10^{-5}$ s⁻¹, we remark fewer perturbations than the $\mu = 0$ case and the shock location is also improved. Indeed, we tend to have an interface speed of sound varying between the Wood and frozen speeds of sound (though closer to the frozen one), which slows down the shock.

## 5.4 | Spherical bubble collapse

This test case aims to present the behaviour of different relaxation rates for UEq during the spherical collapse of a bubble in a free field. To reduce the computational cost, a one-dimensional domain of 3 mm long is used with spherical axi-symmetry to mimic a three-dimensional bubble (Figure 11). The two-phase, axi-symmetric system of equations solved is presented in Appendix D. The domain contains high-pressure water at $p_\infty = 50.6625 \times 10^5$ Pa (50 atm) and $\rho_l = 10^3$ kg/m³. At its origin is positioned an air bubble of initial radius $R_0 = 0.1$ mm. Its state is equivalent to saturated vapor with $p_b = 3550$ Pa ($p_\infty/p_b = 1427$) and $\rho_b = 0.027$ kg/m³. Initial velocities are nulls. One could note that this configuration enforces an initial interface disequilibrium. A non-reflecting boundary condition is used at the far field limit while a symmetry boundary condition is used at the origin of the domain. The mesh contains 150 cells from 0 to 0.3 mm, which corresponds to 100 cells per bubble diameter, and then the grid is stretched non-uniformly to accommodate the large computational domain. The MC and THINC limiters are used.
<figure><img src="image_11.png" /><figcaption>FIGURE 10 Solution of the shock on an initially diffused water-air interface at $t = 600 \ \mu s$ for the proposed UEq modelling with infinite or finite relaxations. The exact and PUEq solutions are also plotted for comparison. The Water1 SG EOS is used for water. (d) is a magnified view of (c).</figcaption></figure>

In the following, we use the radial bubble-wall evolution to compare the performance of the different modellings. Because we are in one dimension of space, we define an effective bubble radius, $R$, as

$$R = \sum_{i=1}^{N} \alpha_{g,i} V_{c,i}, \qquad (49)$$

where $N$ is the total number of grid cells, and $\alpha_{g,i}$ and $V_{c,i}$ are the gas volume fraction and the volume of cell $i$, respectively. The radial bubble-wall evolution is presented in a non-dimensionalized form where

$$t_c = 0.915 R_0 \sqrt{\frac{\rho_l}{p_\infty}} \qquad (50)$$

is the nominal total collapse time, also called Rayleigh collapse time, from its initial (maximum) radius $R_0$⁵³.

We also compute a semi-analytical solution following the Keller–Miksis equation⁵⁴; a compressible form of the Rayleigh–Plesset equation. The Keller–Miksis equation is based on an asymptotic expansion in Mach number of the Bernoulli equation which also assumes that the bubble remains spherical. Its use here is predicated on the idea that errors measured relative to it are
<figure><img src="image_12.png" /><figcaption>FIGURE 11 Problem configuration for a collapsing spherical bubble.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 12 Radial bubble-wall evolution for spherical bubble collapse with $p_\infty/p_b = 1427$ and $N_{R_0} = 50$ cells. Solutions are computed using the PUEq and proposed **UEq** modellings as well as the Keller–Miksis equation. The Water2 SG EOS is used for water. (b) is a magnified view of (a).</figcaption></figure>

larger than any errors associated with the asymptotic expansion and presumption of sphericity inherent to it. This assumption is borne out in the results presented below.

In Figure 12 is shown the solutions obtained from the PUEq and **UEq** modellings as well as from the Keller–Miksis equation. Note that the solutions are only shown until $t = 1.05t_c$, just after the minimum bubble radius is reached, since the subsequent rebounds for large pressure ratios for the Keller–Miksis equation are well-known to be physically inaccurate⁵⁵. We remark that **UEq** with infinite relaxation performs as well as PUEq. Whereas the smaller the finite relaxation rate, the higher the minimum bubble radius and therefore the lower the accuracy. Furthermore, the maximum pressure being directly related to the minimum radius, the maximum pressure for smaller finite relaxation rate is reduced.

We conclude that in the scenario, where we are not looking at waves crossing the interface, it is highly preferable to compute an infinite relaxation rate.

## 5.5 | Non-spherical bubble collapse

This test case aims to present again the behaviour of different relaxation rates for **UEq** but here for a more complex flow: non-spherical bubble collapse. To reduce the computational cost, a two-dimensional domain (12 mm x 12 mm) is used with cylindrical axi-symmetry to mimic a three-dimensional bubble. The two-phase, axi-symmetric system of equations solved is presented in Appendix D. Figure 13 shows the flow configuration considered.

The initial bubble is spherical with radius $R_0 = 0.5$ mm and stand-off distance, distance from the wall to bubble centre, $S = 3R_0/2$. The bubble is filled with non-condensable gas of initial pressure $p_b = 3550$ Pa and density $\rho_b = 0.027$ kg/m³.
<figure><img src="image_14.png" /><figcaption>**FIGURE 13** Problem configuration for a non-spherical bubble collapse near a wall. Two dimensions are used with cylindrical axi-symmetry.</figcaption></figure>

Furthermore, the collapse dynamics are also only weakly sensitive to the internal bubble pressure when the driving pressure differences are large$^{27}$.

The bubble is surrounded by water with a density of $\rho_l = 10^3$ kg/m$^3$ and varying pressure

$$p(r, t = 0) = p_{\infty} + \frac{R_0}{r} (p_b - p_{\infty}) \quad \text{for } r > R_0, \qquad (51)$$

where r is the radial coordinate with origin at the bubble centre. Initial velocities are null. This initialization matches the pressure distribution predicted by the Rayleigh equation for the Besant problem$^{53,56}$. For the configuration considered, it provides a valid approximation of the realistically evolving pressure field and suppresses the formation of spurious pressure waves due to pressure jumps$^{13,57}$. Further, it has been established that this approximation evolves towards an exact solution of the Besant problem within a very short time$^{58}$. We use $p_{\infty} = 50.6625 \times 10^5$ Pa (50 atm), which matches the previous test case and serves as a representation of actual applications involving liquid cavitation.

Unlike the spherical-bubble-collapse test case presented above, one could note that this configuration involves an initial interface equilibrium. Non-reflecting boundary conditions are used at the far field limits while symmetry and wall boundary conditions are used on the bottom and left boundaries, respectively. The mesh contains 400 cells from 0 to 2 mm in the x-direction and 250 cells from 0 to 1.25 mm in the y-direction, which corresponds to 200 cells per bubble diameter, and then the grid is stretched non-uniformly to accommodate the large computational domain. The MC and THINC limiters are used.

The total simulation time is 1.15t* where

$$t^* = R_0 \sqrt{\frac{\rho_l}{\Delta p}} \qquad (52)$$

is an estimate of the collapse time of a bubble collapse near a solid wall$^{59}$, where $\Delta p \equiv p_{\infty} - p_b$ is the driving pressure difference. The wall has a retarding effect on the collapse and thus $t^*$ is longer than the Rayleigh collapse time for spherical collapses ($t_c = 0.915t^*$).

Figure 14 visualizes the flow of the collapsing bubble using a similar setup than what is done by Trummler et al.$^{28}$, meaning using the pressure field $p$ and numerical schlieren $\Phi^{60}$ as

$$\Phi = \exp \left( - \frac{k |\nabla \rho|}{\max |\nabla \rho|} \right), \qquad (53)$$

where $k = 400$ is used to ensure waves in the liquid are visible$^{61,12}$. However herein, because we present results with finite relaxation rates, the flow encounters negative mixture pressures (coming from the liquid). Hence, to avoid any visualization failure when rendering the log-scale pressure field, we present the following non-dimensional pressure field

$$\frac{p + \pi_{\infty, \text{water}}}{p_{\infty}}. \qquad (54)$$

Figure 14 shows results for six selected times and for infinite and finite relaxations. Note that the $\mu = 0$ simulation, which rather represents a limit test of the modelling (no communication between the phases with regard to energy) than a physical
behaviour, was unstable for this test case and is therefore not presented. In the three presented configurations, we observe the non-spherical collapse of the bubble with several characteristic stages:

1. (i), the bubble slowly collapses with the part of the interface further from the wall moving faster and in direction to the wall.

2. (ii), a jet is forming until the two ends of the interface interact.

3. (iii), the jet propagates in direction to the wall and a violent shock wave is emitted.

4. (iv) and (v), the waves emitted in direction to the wall slowly dissipate due to the three-dimensional propagation and then they are reflected by the wall.

5. (vi), those reflected waves travel back to the bubble and interact with the now diffused interface.

As expected from the results of the spherical collapse, we remark that the higher the relaxation rate, the smaller the minimum bubble volume. However, it appears that the released energy between the three cases is similar since they all produce an equivalent pressure peak, about $2 \times 10^8$ Pa, on the wall at the location of the axis of symmetry (Figure 15).

Properly evaluate the results is a difficult task considering the complicated wave patterns in such an advanced test case. However, considering the results obtained on the spherical collapse and the validation of the modelling with infinite relaxation rate made in the literature$^{28,27,62,15}$, we can assume the results with $\mu = \infty$ are the closest to a physical behaviour until the interaction of the wall-reflected waves with the bubble. Indeed, based on the results shown in Section 5.3, this interaction is certainly not properly handled because the interface is now diffused. It is important to note that this situation could also happen in numerous situations, such as when there is a cloud of bubbles where numerous bubbles collapse and emit waves which are therefore impinging on diffused interfaces$^{58,57}$ or when two bubbles interact and collapse near a wall$^{63}$.

# 6 | DISCUSSION ON FINITE PRESSURE RELAXATION RATE

Compared to other models and as stated in the introduction, the Kapila et al.$^{16}$ model and its corresponding numerical solutions$^{17,18}$ have many mathematical and thermodynamical advantages leading to very good results for numerous problems. However, as shown herein, they lack reliability for the problem of the interaction of a wave with a diffused interface. Make use of a finite relaxation rate instead of an infinite one is one solution. In this context, the question of the choice of the rate value arises and the flow topology, i.e. wave travelling through an interface, bubble collapse, cavitation or any kind of mixture phenomenon, governs what should be taken.

For the peculiar case of a wave interacting with a diffused interface, the target is an effective speed of sound within the mixture region (interface) evolving linearly with the volume fraction between the two pure-phase speeds of sound. To obtain such results, we proposed to use Eq. (43), (46) and (47) for which parameters are the thickness of the interface and the speed of sound of each pure phase. One should note that the thickness of the interface may evolve with the flow.

For a resolved bubble collapse, we showed that an infinite relaxation rate produces the best results, at least until the interaction of the wall-reflected waves with the bubble where inaccuracies are expected. We note that during the collapse, no waves travel through the interface but rather a compression is formed at the interface vicinity, leading to the pressure peak, shock wave and jet. These are phenomenologically different, which explains the difference in relaxation rate.

While in the context of mixtures, the choice of the rate value could be evaluated by using for example the collapse time of gas pores within porous media subjected to sudden increase in pressure$^{64}$, or by using experimental information on the negative pressure signal from expansion waves travelling through a liquid before cavitation (appearance of a gas phase and therefore with positive pressure) is reached$^{65}$.

Finally, we propose herein a numerical context in which one can choose a pressure relaxation rate, even if this choice is in itself a challenging task. In particular because the value of the relaxation rate must be adaptive to the local flow topology.

# 7 | CONCLUSION

We propose herein a new method to compute the pressure-disequilibrium model (3) where the correction of the internal energies of the phases is done, unlike the previous method of Saurel et al.$^{17}$, directly on the fluxes (before the relaxations). This allows
<figure><img src="image_15.png" /><figcaption>**FIGURE 14** Numerical schlieren (bottom) and log-scale pressure fields (top) of an air bubble collapsing onto a wall at selected times (i)–(vi) for three relaxation coefficients. Gas volume fraction $\alpha_g$ is shown as a shaded area of decreasing opacity with decreasing $\alpha_g$ (bottom), while the $\alpha_g = 0.5$ bubble interface is shown as a solid curve (top). The initial bubble interface ($t = 0$) is also shown with a red dash line at time (i). The Water2 SG EOS is used for water.</figcaption></figure>

to compute infinite and finite pressure relaxations. This new modelling, called **UEq** herein, has been validated on numerous test cases: transport, classical shock tubes, shock on interface, spherical and non-spherical bubble collapses. When infinite relaxation is used, the proposed modelling (**UEq**) presents almost exactly the same solutions than the PUEq modelling¹⁷ while being more generalized (works for finite relaxations). We also propose a methodology to integrate the finite relaxation terms which optimizes the sub-time steps. The latter being determined under specific conditions on volume fractions and pressures.
<figure><img src="image_16.png" /><figcaption>**FIGURE 15** Pressure on the wall at the axis of symmetry for the non-spherical collapse test case and for three relaxation coefficients. The Water2 SG EOS is used for water.</figcaption></figure>

The influence of computing or not the non-conservative terms $\alpha_k p_k \nabla \cdot \mathbf{u}$ of the internal-energy equations of (3) has been investigated. We remind that whether we compute or not the terms, the scheme is conservative. If an infinite relaxation is done, results are almost identical. Whereas if a finite relaxation is done, another situation arises. When interface problems are simulated, the solutions are definitely better when the non-conservative terms are not computed, especially when the terms have significant values. Indeed, the latter produce an inadequate distribution of the energy between the phases. Results were presented to exhibit this behaviour and we conclude it is better not to compute those terms.

Comparisons against the model of Pelanti & Shyue<sup>18,43</sup> (B3), using total energies instead of internal energies and herein called UEqTotE, were also made. In all the investigated cases, UEqTotE presented similar results than the **UEq NC** modelling (taking into account the non-conservative terms). Meaning it worked well but for the water-air shock-tube case where $\mu = 0$ was considered. Indeed, like for **UEq NC**, the non-conservative terms produced an inadequate distribution of the energy between the phases. One may note that the original wave-propagation scheme proposed in<sup>18,43</sup> has not been tested and this could maybe improve the results when $\mu = 0$.

Through our different tests, we have shown that waves behave differently when they form at the interface (mixture) compared to when they pass through it. Indeed, waves interacting with a diffused interface have shown issues when infinite relaxation is used. Solutions are seriously polluted due to the mixture speed of sound following Wood's relation (1) for which significantly lower speeds than the ones of the pure phases are present. In this context, it has been shown that using a finite relaxation rate is a good candidate to avoid such problems and to allow good wave transmission. For the particular case of waves crossing an interface, we also proposed a methodology to determine the value of the finite rate $\mu$.

When considering bubble collapse, spherical or not, it appears that an infinite rate is preferred to match physical phenomena during the collapse phase. However, interactions of the waves, potentially produced by the bubbles, with diffused interfaces (here bubble wall) are to avoid if infinite relaxation rate is kept. Note that diffused interfaces obviously appear during the course of a simulation but can also simply be present because an initial smearing is used to stabilize the high-order method<sup>39,13,9,15</sup>.

Finally, in this paper, we first raised the problem of interaction of waves with diffused interfaces and its consequences on solutions. We then proposed a solution as an initial step toward better wave propagation in the context of diffused interfaces.

## ACKNOWLEDGMENTS

K.S. and N.F. acknowledge support from the A*MIDEX and the ANR under grants ANR-11-LABX-0092, ANR-11-IDEX-0001-02 and ANR-ASTRID project SNIP ANR-19-ASTR-0016-01. J.C., F.P. and E.D. acknowledge support from the CNES under grant RT-CT-2310000-2001-CNES.
References

1. Russell PS, Giosio DR, Venning JA, Pearce BW, Brandner PA, Ceccio S. Microbubble generation from condensation and turbulent breakup of sheet cavitation. In: 31st Symposium on Naval Hydrodynamics. ; 2016: 11–16.

2. Brennen CE. Cavitation in medicine. *Interface Focus* 2015; 5(5): 20150022.

3. Khokhlova VA, Fowlkes JB, Roberts WW, et al. Histotripsy methods in mechanical disintegration of tissue: Towards clinical applications. *Int. J. hyperthermia* 2015; 31(2): 145–162.

4. Drui F. Modélisation et simulation Eulériennes des écoulements diphasiques à phases séparées et dispersées: développement d'une modélisation unifiée et de méthodes numériques adaptées au calcul massivement parallèle. PhD thesis. Université Paris-Saclay (ComUE), Gif-sur-Yvette; 2017.

5. Abgrall R, Karni S. Computations of compressible multifluids. *J. Comp. Phys.* 2001; 169(2): 594–623.

6. Saurel R, Abgrall R. A simple method for compressible multifluid flows. *SIAM J. Sci. Comp.* 1999; 21(3): 1115–1145.

7. Allaire G, Clerc S, Kokh S. A five-equation model for the simulation of interfaces between compressible fluids. *J. Comp. Phys.* 2002; 181: 577–616.

8. Massoni J, Saurel R, Nkonga B, Abgrall R. Proposition de méthodes et modèles Eulériens pour les problèmes à interfaces entre fluides compressibles en présence de transfert de chaleur: Some models and Eulerian methods for interface problems between compressible fluids with heat transfer. *Int. J. Heat and Mass Transfer* 2002; 45(6): 1287–1307.

9. Coralic V, Colonius T. Finite-volume WENO scheme for viscous compressible multicomponent flows. *J. Comp. Phys.* 2014; 274: 95–121.

10. Beig SA, Johnsen E. Maintaining interface equilibrium conditions in compressible multiphase flows using interface capturing. *J. Comp. Phys.* 2015; 302: 548–566.

11. Meng J, Colonius T. Numerical simulations of the early stages of high-speed droplet breakup. *Shock Waves* 2014; 1–16.

12. Meng J, Colonius T. Numerical simulation of the aerobreakup of a water droplet. *J. Fluid Mech.* 2018; 835: 1108–1135.

13. Tiwari A, Freund J, Pantano C. A diffuse interface model with immiscibility preservation. *J. Comp. Phys.* 2013; 252: 290–309.

14. Rasthofer U, Wermelinger F, Hadiijdoukas P, Koumoutsakos P. Large scale simulation of cloud cavitation collapse. *Procedia Comput. Sci.* 2017; 108: 1763–1772.

15. Schmidmayer K, Bryngelson SH, Colonius T. An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics. *J. Comp. Phys.* 2020; 402: 109080.

16. Kapila A, Menikoff R, Bdzil J, Son S, Stewart D. Two-phase modeling of DDT in granular materials: Reduced equations. *Phys. Fluids* 2001; 13: 3002–3024.

17. Saurel R, Petitpas F, Berry R. Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures. *J. Comp. Phys.* 2009; 228(5): 1678–1712.

18. Pelanti M, Shyue KM. A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves. *J. Comp. Phys.* 2014; 259: 331–357.

19. Petitpas F, Massoni J, Saurel R, Lapebie E, Munier L. Diffuse interface models for high speed cavitating underwater systems. *Int. J. Multiphase Flow* 2009; 35(8): 747–759.

20. Petitpas F, Saurel R, Franquet E, Chinnayya A. Modelling detonation waves in condensed energetic materials: Multiphase CJ conditions and multidimensional computations. *Shock Waves* 2009; 19(5): 377–401.
21. Schmidmayer K, Petitpas F, Daniel E, Favrie N, Gavrilyuk S. A model and numerical method for compressible flows with capillary effects. *J. Comp. Phys.* 2017; 334: 468–496.

22. Chiocchetti S, Peshkov I, Gavrilyuk S, Dumbser M. High order ADER schemes and GLM curl cleaning for a first order hyperbolic formulation of compressible flow with surface tension. *J. Comp. Phys.* 2021; 426: 109898.

23. Schmidmayer K, Petitpas F, Daniel E. Adaptive Mesh Refinement algorithm based on dual trees for cells and faces for multiphase compressible flows. *J. Comp. Phys.* 2019; 388: 252–278.

24. Dorschner B, Biasiori-Poulanges L, Schmidmayer K, El-Rabii H, Colonius T. On the formation and recurrent shedding of ligaments in droplet aerobreakup. *J. Fluid Mech.* 2020; 904: A20.

25. Favrie N, Gavrilyuk SL, Saurel R. Solid-fluid diffuse interface model in cases of extreme deformations. *J. Comp. Phys.* 2009; 228(16): 6037–6077.

26. Ndanou S, Favrie N, Gavrilyuk S. Multi-solid and multi-fluid diffuse interface model: Applications to dynamic fracture and fragmentation. *J. Comp. Phys.* 2015; 295: 523–555.

27. Pishchalnikov YA, Behnke-Parks WM, Schmidmayer K, et al. High-speed video microscopy and numerical modeling of bubble dynamics near a surface of urinary stone. *J. Acoust. Soc. Am.* 2019; 146: 516–531.

28. Trummler T, Bryngelson SH, Schmidmayer K, Schmidt SJ, Colonius T, Adams NA. Near-surface dynamics of a gas bubble collapsing above a crevice. *J. Fluid Mech.* 2020; 899: A16.

29. Wood AB. A textbook of sound. *G. Bell and Sons LTD, London* 1930.

30. Favrie N, Schmidmayer K, Massoni J. A multiphase irreversible-compaction model for granular-porous materials. *Continuum Mech. Thermodyn.* 2021.

31. Schmidmayer K, Petitpas F, Le Martelot S, Daniel E. ECOGEN: An open-source tool for multiphase, compressible, multiphysics flows. *Comp. Phys. Com.* 2020; 251: 107093.

32. Re B, Abgrall R. A pressure-based method for weakly compressible two-phase flows under a Baer–Nunziato type model with generic equations of state and pressure and velocity disequilibrium. *Int. J. Numer. Meth. Fl.* 2022.

33. Baer MR, Nunziato JW. A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials. *Int. J. of Multiphase Flow* 1986; 12: 861–889.

34. Chiocchetti S, Müller C. A Solver for Stiff Finite-Rate Relaxation in Baer–Nunziato Two-Phase Flow Models. In: Springer. 2020 (pp. 31–44).

35. Wallis GB. *One-dimensional two-phase flow*. McGraw-Hill, 1969.

36. So KK, Hu XY, Adams NA. Anti-diffusion interface sharpening technique for two-phase compressible flow simulations. *J. Comp. Phys.* 2012; 231(11): 4304–4323.

37. Shyue K, Xiao F. An Eulerian interface sharpening algorithm for compressible two-phase flow: The algebraic THINC approach. *J. Comp. Phys.* 2014; 268: 326–354.

38. Jiang GS, Shu CW. Efficient implementation of weighted ENO schemes. *J. Comp. Phys.* 1996; 126(1): 202–228.

39. Johnsen E, Colonius T. Implementation of WENO schemes in compressible multicomponent flow problems. *J. Comp. Phys.* 2006; 219(2): 715–732.

40. Shukla RK, Pantano C, Freund JB. An interface capturing method for the simulation of multi-phase compressible flows. *J. Comp. Phys.* 2010; 229(19): 7411–7439.

41. Bryngelson SH, Schmidmayer K, Coralic V, Meng JC, Maeda K, Colonius T. MFC: An open-source high-order multi-component, multi-phase, and multi-scale compressible flow solver. *Comp. Phys. Com.* 2020; 266: 107396.
42. Schwartzkopff T, Dumbser M, Munz CD. Fast high order ADER schemes for linear hyperbolic equations. *J. Comp. Phys.* 2004; 197(2): 532–539.

43. Pelanti M, Shyue KM. A numerical model for multiphase liquid-vapor-gas flows with interfaces and cavitation. *Int. J. Multiphase Flow* 2019; 113: 208–230.

44. Le Métayer O, Massoni J, Saurel R. Elaborating equations of state of a liquid and its vapor for two-phase flow models. *Int. J. Therm. Sci.* 2004; 43: 265–276.

45. Le Métayer O, Saurel R. The Noble-Abel Stiffened-Gas equation of state. *Phys. Fluids* 2016; 28(4): 046102.

46. Godunov SK. A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics. *Mat. Sb.* 1959; 89(3): 271–306.

47. Toro E. *Riemann solvers and numerical methods for fluid dynamics*. Berlin: Springer Verlag. 1997.

48. Sweby PK. High resolution schemes using flux limiters for hyperbolic conservation laws. *SIAM J. Num. Anal.* 1984; 21(5): 995–1011.

49. Van Leer B. Towards the ultimate conservative difference scheme. II. Monotonicity and conservation combined in a second-order scheme. *J. Comp. Phys.* 1974; 14(4): 361–370.

50. Van Leer B. Towards the ultimate conservative difference scheme III. Upstream-centered finite-difference schemes for ideal compressible flow. *J. Comp. Phys.* 1977; 23(3): 263–275.

51. Hindmarsh AC. ODEPACK, a systematized collection of ODE solvers. *Scientific computing* 1983: 55–64.

52. Petzold L. Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations. *SIAM J. Sci. Stat. Comput.* 1983; 4(1): 136–148.

53. Brennen CE. *Cavitation and bubble dynamics*. Oxford University Press. 1995.

54. Keller JB, Miksis M. Bubble oscillations of large amplitude. *J. Acoust. Soc. Am.* 1980; 68(2): 628–633.

55. Fuster D, Dopazo C, Hauke G. Liquid compressibility effects during the collapse of a single cavitating bubble. *J. Acoust. Soc. Am.* 2011; 129(1): 122–131.

56. Besant W. *A Treatise on Hydrostatics and Hydrodynamics*. Bell, London: Deighton. 1859.

57. Tiwari A, Pantano C, Freund J. Growth-and-collapse dynamics of small bubble clusters near a wall. *J. Fluid Mech.* 2015; 775: 1–23.

58. Rasthofer U, Wermelinger F, Karnakov P, Šukys J, Koumoutsakos P. Computational study of the collapse of a cloud with 12 500 gas bubbles in a liquid. *Phys. Rev. Fluids* 2019; 4(6): 063602.

59. Plesset MS, Chapman RB. Collapse of an initially spherical vapour cavity in the neighbourhood of a solid boundary. *J. Fluid Mech.* 1971; 47(2): 283–290.

60. Quirk J, Karni S. On the dynamics of a shock-bubble interaction. *J. Fluid Mech.* 1996; 318: 129–163.

61. Johnsen E. *Numerical simulations of non-spherical bubble collapse*. PhD thesis. California Institute of Technology, Pasadena; 2007.

62. Beig SA, Aboulhasanzadeh B, Johnsen E. Temperatures produced by inertially collapsing bubbles near rigid surfaces. *J. Fluid Mech.* 2018; 852: 105–125.

63. Alahyari Beig S. *A Computational Study of the Inertial Collapse of Gas Bubbles Near a Rigid Surface*. PhD thesis. University of Michigan, Ann Arbor; 2018.

64. Favrie N, Schmidmayer K, Massoni J. A multiphase irreversible-compaction model for granular-porous materials. *Continuum Mech. Thermodyn.* 2022; 34.1: 217–241.
65. Veilleux JC. *Pressure and stress transients in autoinjector devices*. PhD thesis. California Institute of Technology, Pasadena; 2019.

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

# APPENDIX

## A MECHANICAL-EQUILIBRIUM MODEL OF KAPILA ET AL.

The mechanical-equilibrium model of Kapila et al.¹⁶ is expressed for a two-phase flow (N = 2) as

$$ q = \begin{bmatrix} \alpha_1 \\ \alpha_1 \rho_1 \\ \alpha_2 \rho_2 \\ \rho \mathbf{u} \\ \rho E \end{bmatrix}, \quad F = \begin{bmatrix} \alpha_1 \mathbf{u} \\ \alpha_1 \rho_1 \mathbf{u} \\ \alpha_2 \rho_2 \mathbf{u} \\ \rho \mathbf{u} \otimes \mathbf{u} + p \mathbf{I} \\ (\rho E + p) \mathbf{u} \end{bmatrix}, \quad h = \begin{bmatrix} -\alpha_1 - K \\ 0 \\ 0 \\ \mathbf{0} \\ 0 \end{bmatrix}, \quad (A1) $$

$r = g = m = [0, 0, 0, 0, 0]^T$, K is

$$ K = \frac{\rho_2 c_2^2 - \rho_1 c_1^2}{\frac{\rho_2 c_2^2}{\alpha_2} + \frac{\rho_1 c_1^2}{\alpha_1}}, \quad (A2) $$

and $K \nabla \cdot \mathbf{u}$ represents expansion and compression of each phase in mixture regions. For this model, the mixture speed of sound follows Wood's²⁹.

## B PRESSURE-DISEQUILIBRIUM MODEL BASED ON TOTAL ENERGIES

The pressure-disequilibrium model based on total energies¹⁸,⁴³, which we call UEqTotE, is expressed for N phases as

$$ q = \begin{bmatrix} \alpha_k \\ \alpha_k \rho_k \\ \rho \mathbf{u} \\ \alpha_k \rho_k E_k \end{bmatrix}, \quad F = \begin{bmatrix} \alpha_k \mathbf{u} \\ \alpha_k \rho_k \mathbf{u} \\ \rho \mathbf{u} \otimes \mathbf{u} + p \mathbf{I} \\ \alpha_k \rho_k E_k \mathbf{u} \end{bmatrix}, \quad h = \begin{bmatrix} -\alpha_k \\ 0 \\ \mathbf{0} \\ -Y_k p + \alpha_k p_k \end{bmatrix}, \quad (B3) $$

$$ g = \begin{bmatrix} 0 \\ 0 \\ \mathbf{0} \\ Y_k \end{bmatrix}, \quad m = \begin{bmatrix} 0 \\ 0 \\ \mathbf{0} \\ p \mathbf{u} \end{bmatrix}, \quad r = \begin{bmatrix} \delta p_k \\ 0 \\ \mathbf{0} \\ -p_I \delta p_k \end{bmatrix}, $$

where

$$ E_k = e_k + \frac{1}{2} \| \mathbf{u} \|^2. \quad (B4) $$

For this model, pressures are also in disequilibrium and the total energy equation of the mixture is replaced by the total-energy equation for each phase. Nevertheless, conservation of the mixture total energy (12) can be recovered. However, unlike model (3), there is no need to include (12) in our computations to numerically ensure conservation of the total energy. Indeed, the sum of the non-conservative terms of each phase cancels. Further, the mixture speed of sound is the same than for model (3).
C PHASE INTERFACIAL SOUND SPEED

We recall that for a finite pressure relaxation, the pressure equations are solved

$$ \frac{\partial p_k}{\partial t} = - \frac{\rho_k c_{I,k}^2}{\alpha_k} \mu \Delta p_k, \qquad (C5) $$

where

$$ c_{I,k}^2 = \frac{\left( \frac{p_I}{\rho_k^2} - \frac{\partial e_k}{\partial \rho_k} \right)_{p_k}}{\left( \frac{\partial e_k}{\partial p_k} \right)_{\rho_k}} \quad \text{and} \quad \Delta p_k = \sum_{j \neq k}^{N} (p_k - p_j), \qquad (C6) $$

and where $c_{I,k}^2$ needs to be expressed for specific EOS.

In the context of ideal and SG⁴⁴ EOS, one can write the internal energies as

$$ e_k = \frac{p_k + \gamma_k \pi_{\infty,k}}{(\gamma_k - 1) \rho_k}. \qquad (C7) $$

Hence

$$ \left. \frac{\partial e_k}{\partial \rho_k} \right|_{p_k} = - \frac{p_k + \gamma_k \pi_{\infty,k}}{(\gamma_k - 1) \rho_k^2} \quad \text{and} \quad \left. \frac{\partial e_k}{\partial p_k} \right|_{\rho_k} = \frac{1}{(\gamma_k - 1) \rho_k}. \qquad (C8) $$

This leads to

$$ c_{I,k}^2 = \frac{(\gamma_k - 1) p_I + p_k + \gamma_k \pi_{\infty,k}}{\rho_k}. \qquad (C9) $$

D TWO-PHASE, PRESSURE-DISEQUILIBRIUM MODEL WITH AXI-SYMMETRY

In the context of cylindrical or spherical axi-symmetry, the pressure-disequilibrium model based on internal energies¹⁷ for two phases is written as follows. Note that the axi-symmetric terms herein are solved as source terms.

Cylindrical coordinates

$$ \left\{ \begin{aligned} & \frac{\partial \alpha_1}{\partial t} + u \frac{\partial \alpha_1}{\partial r} + v \frac{\partial \alpha_1}{\partial z} &&= \mu (p_1 - p_2), \\ & \frac{\partial \alpha_1 \rho_1}{\partial t} + \frac{\partial \alpha_1 \rho_1 u}{\partial r} + \frac{\partial \alpha_1 \rho_1 v}{\partial z} &&= \frac{\alpha_1 \rho_1 u}{r}, \\ & \frac{\partial \alpha_2 \rho_2}{\partial t} + \frac{\partial \alpha_2 \rho_2 u}{\partial r} + \frac{\partial \alpha_2 \rho_2 v}{\partial z} &&= \frac{\alpha_2 \rho_2 u}{r}, \\ & \frac{\partial \rho u}{\partial t} + \frac{\partial \rho u^2 + p}{\partial r} + \frac{\partial \rho u v}{\partial z} &&= \frac{\rho u^2}{r}, \\ & \frac{\partial \rho v}{\partial t} + \frac{\partial \rho u v}{\partial r} + \frac{\partial \rho v^2 + p}{\partial z} &&= \frac{\rho u v}{r}, \\ & \frac{\partial \alpha_1 \rho_1 e_1}{\partial t} + \frac{\partial \alpha_1 \rho_1 e_1 u}{\partial r} + \frac{\partial \alpha_1 \rho_1 e_1 v}{\partial z} + \alpha_1 p_1 \left( \frac{\partial u}{\partial r} + \frac{\partial v}{\partial z} \right) &&= - \frac{\alpha_1 (\rho_1 e_1 + p_1) u}{r} - \mu p_I (p_1 - p_2), \\ & \frac{\partial \alpha_2 \rho_2 e_2}{\partial t} + \frac{\partial \alpha_2 \rho_2 e_2 u}{\partial r} + \frac{\partial \alpha_2 \rho_2 e_2 v}{\partial z} + \alpha_2 p_2 \left( \frac{\partial u}{\partial r} + \frac{\partial v}{\partial z} \right) &&= - \frac{\alpha_2 (\rho_2 e_2 + p_2) u}{r} - \mu p_I (p_1 - p_2), \end{aligned} \right. \qquad (D10) $$

with $\mathbf{u} = (u, v)^T$ in two dimensions and with the total-energy equation

$$ \frac{\partial \rho E}{\partial t} + \frac{\partial (\rho E + p) u}{\partial r} + \frac{\partial (\rho E + p) v}{\partial z} = \frac{(\rho E + p) u}{r}. \qquad (D11) $$
Spherical coordinates

$$ \left\{ \begin{aligned} \frac{\partial \alpha_1}{\partial t} + u \frac{\partial \alpha_1}{\partial r} &= \mu (p_1 - p_2), \\ \frac{\partial \alpha_1 \rho_1}{\partial t} + \frac{\partial \alpha_1 \rho_1 u}{\partial r} &= \frac{2 \alpha_1 \rho_1 u}{r}, \\ \frac{\partial \alpha_2 \rho_2}{\partial t} + \frac{\partial \alpha_2 \rho_2 u}{\partial r} &= \frac{2 \alpha_2 \rho_2 u}{r}, \\ \frac{\partial \rho u}{\partial t} + \frac{\partial \rho u^2 + p}{\partial r} &= \frac{2 \rho u^2}{r}, \\ \frac{\partial \alpha_1 \rho_1 e_1}{\partial t} + \frac{\partial \alpha_1 \rho_1 e_1 u}{\partial r} + \alpha_1 p_1 \frac{\partial u}{\partial r} &= - \frac{2 \alpha_1 (\rho_1 e_1 + p_1) u}{r} - \mu p_I (p_1 - p_2), \\ \frac{\partial \alpha_2 \rho_2 e_2}{\partial t} + \frac{\partial \alpha_2 \rho_2 e_2 u}{\partial r} + \alpha_2 p_2 \frac{\partial u}{\partial r} &= - \frac{2 \alpha_2 (\rho_2 e_2 + p_2) u}{r} - \mu p_I (p_1 - p_2), \end{aligned} \right. \quad (D12) $$

with $u = u$ in one dimension and with the total-energy equation

$$ \frac{\partial \rho E}{\partial t} + \frac{\partial (\rho E + p) u}{\partial r} = \frac{2 (\rho E + p) u}{r}. \quad (D13) $$