# An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics

Kevin Schmidmayer, Spencer Bryngelson, Tim Colonius

► To cite this version:

Kevin Schmidmayer, Spencer Bryngelson, Tim Colonius. An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics. Journal of Computational Physics, 2020, 402, pp.109080. 10.1016/j.jcp.2019.109080. hal-03649385

HAL Id: hal-03649385

https://hal.science/hal-03649385

Submitted on 22 Apr 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
An assessment of multicomponent flow models and interface
capturing schemes for spherical bubble dynamics

Kevin Schmidmayer*, Spencer H. Bryngelson, Tim Colonius

Division of Engineering and Applied Science, California Institute of Technology
1200 E. California Blvd., Pasadena, CA 91125, USA

Abstract

Numerical simulation of bubble dynamics and cavitation is challenging; even the seemingly
simple problem of a collapsing spherical bubble is difficult to compute accurately with a gen-
eral, three-dimensional, compressible, multicomponent flow solver. Difficulties arise due to
both the physical model and the numerical method chosen for its solution. We consider the 5-
equation model of Allaire et al. [1] and Massoni et al. [2], the 5-equation model of Kapila et al.
[3], and the 6-equation model of Saurel et al. [4] as candidate approaches for spherical bubble
dynamics, and both MUSCL and WENO interface-capturing methods are implemented and
compared. We demonstrate the inadequacy of the traditional 5-equation model for spherical
bubble collapse problems and explain the corresponding advantages of the augmented model
of Kapila et al. [3] for representing this phenomenon. Quantitative comparisons between the
augmented 5-equation and 6-equation models for three-dimensional bubble collapse problems
demonstrate the versatility of the pressure-disequilibrium model. Lastly, the performance
of the pressure-disequilibrium model for representing a three-dimensional spherical bubble
collapse for different bubble interior/exterior pressure ratios is evaluated for different numer-
ical methods. Pathologies associated with each factor and their origins are identified and
discussed.

Keywords: Bubble dynamics, interface-capturing schemes, diffuse-interface method,
multiphase flow, compressible flow

<footnote>*Corresponding author</footnote>

<footnote>Email addresses: kevinsch@caltech.edu (Kevin Schmidmayer), spencer@caltech.edu (Spencer H. Bryngelson), colonius@caltech.edu (Tim Colonius)</footnote>
# 1. Introduction

Amongst other features, cavitation involves the growth and collapse of a gas bubble in a liquid. Many applications require a detailed understanding of this process in or near soft materials, including biological tissues for medical purposes [5–9] and polymeric coatings and biofouling in industry [10]. Preliminary studies have shown that bubble dynamics are sensitive to the properties of these materials [11, 12], motivating a comprehensive multi-scale theory capable of predicting complex bubble cavitation.

Before considering the viscoelasticity of soft materials, accurate algorithms for bubble dynamics in Newtonian liquids must be developed. Indeed, even the seemingly simple problem of a collapsing spherical bubble is challenging to compute accurately with general, three-dimensional (3D), fully-compressible computational methods for a significant range of bubble/ambient pressure ratios (and thus interface Mach numbers). Here, we use this problem as a case study for the ability of a physical model, and its coupled numerical method, to predict bubble dynamics generally.

Computational models for multi-component flows can generally be categorized as either interface tracking and interface capturing. Interface-tracking methods treat material interfaces as sharp features in the flow, and thus fluids and interfacial physics are modeled separately. Such methods are generally subdivided into free-Lagrange methods [10, 13, 14], front-tracking methods [15, 16], and level-set/ghost fluid methods [8, 17–20]. Unfortunately, interface-tracking schemes do not naturally enforce discrete or total conservation properties. Herein, we instead use interface-capturing methods because the discrete conservation properties they naturally achieve are of principal interest to many of the key phenomenologies outlined above. Such methods combine a multicomponent flow model with shock-capturing finite volume method. Their discrete-level conservation allows the compressibility of all phases and mixtures to be represented on the computational grid and interfaces appear and vanish naturally, irrespective of their corresponding density ratio. We note that finite difference shock-capturing methods can also be used and fulfill these requirements [21], though
they do not generally improve performance when compared to finite volume implementations. We also prefer interface capturing here as they are generally more efficient than interface-tracking methods [22] and their complexity does not increase with simulation dimensionality or the number of components. Here, we assess the difficulties that arise during a spherical bubble collapse from the physical multicomponent flow model and its coupling to the numerical method.

The mechanical-equilibrium multicomponent model of Allaire et al. [1] and Massoni et al. [2] has been widely used and can faithfully represent shock-induced collapses [23–26] and droplet atomization [27, 28]. Unfortunately, this model cannot predict the collapse time and minimum radius of the Rayleigh collapse problem [29, 30]. This problem can be averted via the thermodynamically consistent model of Kapila et al. [3] [29, 30], which includes a term ($K\nabla \cdot \mathbf{u}$) in the volume-fraction evolution equation that follows from an asymptotic expansion of the total-disequilibrium model of Baer and Nunziato [31]. This additional term represents compressibility in mixture regions, though unfortunately it leads to numerical instabilities during strong compression and expansion near the interface [4, 32]. This is because the equations cannot be cast in a purely hyperbolic form, and thus a unique solution to the associated Riemann problem does not exist [33]. Instead, we propose using a pressure-disequilibrium model [4], which relaxes the phase-specific pressures algorithmically at each time step and averts the stability issues of the $K\nabla \cdot \mathbf{u}$ term. This model theoretically converges to the mechanical-equilibrium model of Kapila et al. [3] under mesh refinement, and while it has been utilized for cavitating flows [4, 34], detonating flows [35], surface-tension driven flows [36], droplet atomization [37, 38], and fracture and fragmentation in ductile materials [39, 40], it has not been applied to bubble dynamics or particularly to collapsing bubbles. We note that other, more modern pressure-disequilibrium models are also viable [41], though they do not directly address challenges associated with simulating bubble cavitation.

The multicomponent flow models are solved using shock-capturing finite-volume schemes and Riemann solvers for the fluxes [42, 43]. High-order spatial reconstructions, such as
MUSCL [43–45] and WENO [23, 29, 46, 47], are often used, along with their variants WENO-Z [48], WENO-CU6 [49, 50], and TENO [51]. Herein, we will consider MUSCL and the WENO of Jiang and Shu [46], coupled with the HLLC approximate Riemann solver [4, 23, 43] as standard approaches for solving the multicomponent flow equations. Following the usual procedure, these are coupled to total-variation-diminishing time integrators as an attempt to suppress spurious oscillations at material interfaces under refinement [42, 43, 52, 53].

We first present the diffuse-interface multicomponent models in section 2. The numerical methods we employ to solve the resulting equations are outlined in section 3. The setup of the spherical-bubble-collapse problems we consider are presented in section 4. In section 5.1 we demonstrate and explain the utility of the $K\nabla \cdot \mathbf{u}$ term in the mechanical-equilibrium models. The convergence and behavior of this improved equilibrium model and the usual pressure-disequilibrium model are studied in section 5.2 for the collapse and rebound of spherical bubbles. Artifacts of the numerical methods we consider are examined in section 5.3, including an investigation of interface sharpening techniques in section 5.4. Finally, the pathologies identified are discussed in section 6.

## 2. Multicomponent flow models

The compressible multicomponent flow models we present can all be written as

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

where $\mathbf{q}$ is the state vector, $\mathbf{F}$ is the flux tensor, $\mathbf{u}$ is the velocity field, and $\mathbf{h}$ and $\mathbf{r}$ are non-conservative quantities we describe subsequently. We only consider mechanical-equilibrium models that formally conserve mass, momentum, and total energy, and neglect the effects of viscosity, phase change and surface tension.
## 2.1. Mechanical-equilibrium model of Allaire et al. [1] and of Massoni et al. [2]

We first consider the mechanical-equilibrium model of Allaire et al. [1] and Massoni et al. [2], which we call the 5-equation model. For a two-component flow, it is

$$ \mathbf{q} = \begin{bmatrix} \alpha_1 \\ \alpha_1 \rho_1 \\ \alpha_2 \rho_2 \\ \rho \mathbf{u} \\ \rho E \end{bmatrix}, \quad \mathbf{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 \mathbf{h} = \begin{bmatrix} -\alpha_1 \\ 0 \\ 0 \\ \mathbf{0} \\ 0 \end{bmatrix}, \quad \mathbf{r} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ \mathbf{0} \\ 0 \end{bmatrix}, \quad (2) $$

where $\rho$, $\mathbf{u}$, and $p$ are the mixture density, velocity, and pressure, respectively, and $\alpha_k$ is the volume fraction, for which $k$ indicates the phase index. The mixture total energy is

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

where $e$ is the mixture specific internal energy

$$ e = \sum_{k=1}^{2} Y_k e_k (\rho_k, p). \quad (4) $$

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

$$ Y_k = \frac{\alpha_k \rho_k}{\rho}. \quad (5) $$

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

$$ p_g = (\gamma_g - 1) \rho_g e_g, \quad (6) $$

where $\gamma_g = 1.4$, and the liquid is modeled by the stiffened-gas EOS

$$ p_l = (\gamma_l - 1) \rho_l e_l - \gamma_l \pi_{\infty}, \quad (7) $$

where $\gamma_l$ and $\pi_{\infty}$ are model parameters [54]. We note that the stiffened-gas EOS is not a complete EOS but is reasonably accurate at the conditions considered here [54]. More general EOS could easily be substituted in our framework. The mixture quantities are

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

$$ \rho c^2 = \sum_{k=1}^{2} \frac{\alpha_k \rho_k c_k^2}{\beta (\gamma_k - 1)}, \quad \beta = \sum_{k=1}^{2} \frac{\alpha_k}{\gamma_k - 1}, \qquad (9) $$

where c is the mixture speed of sound, and $c_k$ and $\gamma_k$ are the speed of sound and polytropic coefficient of phase k. We note that while this model conserves mass, momentum, and total energy, it does not strictly obey the second law of thermodynamics [1, 2, 36].

## 2.2. Mechanical-equilibrium model of Kapila et al. [3]

The thermodynamically consistent mechanical-equilibrium model of Kapila et al. [3], which we call the 5-equation model with $K\nabla \cdot \mathbf{u}$, has

$$ \mathbf{q} = \begin{bmatrix} \alpha_1 \\ \alpha_1 \rho_1 \\ \alpha_2 \rho_2 \\ \rho \mathbf{u} \\ \rho E \end{bmatrix}, \quad \mathbf{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 \mathbf{h} = \begin{bmatrix} -\alpha_1 - K \\ 0 \\ 0 \\ \mathbf{0} \\ 0 \end{bmatrix}, \quad \mathbf{r} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ \mathbf{0} \\ 0 \end{bmatrix}, \qquad (10) $$

where only $\mathbf{h}$ is different from (2). Here, 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}}, \qquad (11) $$

and $K\nabla \cdot \mathbf{u}$ represents expansion and compression of each phase in mixture regions. In this case, the mixture speed of sound follows from

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

which is also the Wood speed of sound [55, 56].
## 2.3. Pressure-disequilibrium model of Saurel et al. [4]

The pressure-disequilibrium model of Saurel et al. [4], which we call the 6-equation model, is expressed as

$$ \mathbf{q} = \begin{bmatrix} \alpha_1 \\ \alpha_1\rho_1 \\ \alpha_2\rho_2 \\ \rho\mathbf{u} \\ \alpha_1\rho_1 e_1 \\ \alpha_2\rho_2 e_2 \end{bmatrix}, \quad \mathbf{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} \\ \alpha_1\rho_1 e_1\mathbf{u} \\ \alpha_2\rho_2 e_2\mathbf{u} \end{bmatrix}, \quad \mathbf{h} = \begin{bmatrix} -\alpha_1 \\ 0 \\ 0 \\ \mathbf{0} \\ \alpha_1 p_1 \\ \alpha_2 p_2 \end{bmatrix}, \quad \mathbf{r} = \begin{bmatrix} \mu\delta p \\ 0 \\ 0 \\ \mathbf{0} \\ -\mu p_I \delta p \\ \mu p_I \delta p \end{bmatrix}, \quad (13) $$

where $\mathbf{r}$ represents the relaxation of pressures between the phases with coefficient $\mu$. The interfacial pressure is

$$ p_I = \frac{z_2 p_1 + z_1 p_2}{z_1 + z_2}, \quad (14) $$

where $z_k = \rho_k c_k$ is the acoustic impedance of the phase $k$, and

$$ \delta p = p_1 - p_2, \quad (15) $$

is the pressure difference between the two phases. Since $p_1 \neq p_2$ 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. \quad (16) $$

We note that (16) is redundant when the internal energy equations are also computed. However, in practice we include it in our computations to ensure that the total energy is numerically conserved, and thus preserve a correct treatment of shock waves (more details can be found in Saurel et al. [4]).

The mixture speed of sound is defined according to

$$ c^2 = \sum_{k=1}^{2} Y_k c_k^2. \quad (17) $$
After applying the infinite pressure-relaxation procedure detailed in section 3.3, the effective mixture speed of sound matches (12). We will discuss the influence of sound speed for interface problems in section 5.3.2.

### 3. Numerical methods

We solve (1) 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 [42] scheme

$$ \mathbf{q}_i^{n+1} = \mathbf{q}_i^n - \frac{\Delta t}{V_i} \left( \sum_{s=1}^{N} A_s \mathbf{F}_s^* \cdot \mathbf{n}_s + \mathbf{h}(\mathbf{q}_i^n) \sum_{s=1}^{N} A_s \mathbf{u}_s^* \cdot \mathbf{n}_s \right), \quad (18) $$

where $n$ is the time-step index. The relaxation terms, if any, are then solved using the procedure detailed in 3.3 to complete the time-step integration. Here, we label this basic first-order-accurate finite-volume scheme as FV1. We also utilize both MUSCL and WENO spatial reconstructions of the primitive state variables; these are presented in the following subsections. We note that reconstructing the conservative variables instead leads to spurious oscillations near material interfaces [23], 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 [4, 23, 43], giving the flux tensor and flow-velocity vector $\mathbf{F}_s^*$ and $\mathbf{u}_s^*$, respectively. The solution of (18) is restricted by the usual CFL criterion.

#### 3.1. MUSCL scheme

We use the second-order-accurate MUSCL scheme of Schmidmayer et al. [45] (labeled here as MUSCL2) with two-step time integration

$$ \mathbf{q}_i^{n+\frac{1}{2}} = \mathbf{q}_i^n + \frac{1}{2} \Delta t \mathbf{L}(\mathbf{q}_i^n), \quad (19) $$

$$ \mathbf{q}_i^{n+1} = \mathbf{q}_i^n + \Delta t \mathbf{L} \left( \mathbf{q}_i^{n+\frac{1}{2}} \right), \quad (20) $$
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 for the second step and the usual piece-wise linear MUSCL reconstruction [43] is used on the primitive variables. The monotonized central (MC) [52] slope limiter is employed as an attempt to minimize interface diffusion and its behavior is investigated in section 5.3. This method has been previously implemented for the pressure-disequilibrium model [4, 34–36, 38, 40, 45, 57].

## 3.2. WENO scheme

We also implement third- and fifth-order accurate WENO schemes for comparison purposes (labeled here as WENO3 and WENO5, respectively). In this case, the time derivative is computed via the third-order TVD Runge–Kutta algorithm [53]

$$ \mathbf{q}_i^{(1)} = \mathbf{q}_i^n + \Delta t \mathbf{L} (\mathbf{q}_i^n), \qquad (21) $$

$$ \mathbf{q}_i^{(2)} = \frac{3}{4}\mathbf{q}_i^n + \frac{1}{4}\mathbf{q}_i^{(1)} + \frac{1}{4}\Delta t \mathbf{L} (\mathbf{q}_i^{(1)}), \qquad (22) $$

$$ \mathbf{q}_i^{n+1} = \frac{1}{3}\mathbf{q}_i^n + \frac{2}{3}\mathbf{q}_i^{(2)} + \frac{2}{3}\Delta t \mathbf{L} (\mathbf{q}_i^{(2)}). \qquad (23) $$

This method has previously been implemented for the pressure-equilibrium models of Allaire et al. [1] and of Massoni et al. [2] [23–25, 27, 28], and of Kapila et al. [3] [29, 30, 58]; here, we also utilize it for the pressure-disequilibrium model of Saurel et al. [4].

## 3.3. Pressure-relaxation procedure

The pressure-disequilibrium model (13) requires stiff pressure relaxation to converge to a single, equilibrium pressure. We use the infinite-relaxation procedure of Saurel et al. [4]. At each time step it solves the non-relaxed, hyperbolic equations ($\mu \to 0$) using (18), then relaxes the disequilibrium pressures for $\mu \to +\infty$. Specifically, after manipulation of the equations realized in agreement with thermodynamic considerations and saturation constraint ($\sum_k \alpha_k = 1$), (13) is solved with respect to

$$ \sum_k (\alpha \rho)_k v_k (p) = 1, \qquad (24) $$
where $(\alpha\rho)_k$ are constant during the relaxation process and $v_k(p)$ are the specific volumes of each state determined with the help of the EOS as

$$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, \quad (25)$$

where superscript 0 indicates the hyperbolic step index. We ultimately solve (24) using the Newton-Raphson method to find the relaxed pressure and the phase densities and volume fractions are determined for the next step.

This relaxation procedure is combined with a re-initialization procedure to ensure the conservation of total energy, and thus converges to the mechanical-equilibrium model of Kapila et al. [3] (10). This follows from the EOS and the mixture total-energy conservation law as

$$p = \frac{\rho e - \sum_k \frac{\alpha_k \gamma_k \pi_{\infty k}}{\gamma_k - 1}}{\sum_k \frac{\alpha_k}{\gamma_k - 1}}, \quad (26)$$

where $\rho e$ is determined using (16). Once $p$ is determined, the internal energies of the phases are reinitialized using their respective EOS. When multi-stage time integration is used, these procedures are performed at each stage. Thus, there is only one pressure at the end of each stage and the reconstructed variables are the same for all models. As a result, simulations of the pressure-disequilibrium model are only about 5% more expensive than the models of Allaire et al. [1], Massoni et al. [2] and Kapila et al. [3] for the spherical-bubble-collapse cases we consider subsequently.

## 4. Setup of the spherical-bubble-collapse problem

As a step towards understanding the practical differences between the presented models and methods, we consider the behavior of a collapsing spherical bubble. The problem setup is shown in Figure 1. We initialize the bubble with radius $R_0$ and the computational domain has size $L = 320R_0$, which is sufficiently large to avoid boundary effects. Initially, the bubble has a uniform internal pressure $p_b$, and the exterior pressure increases gradually up to the far-field pressure $p_\infty$ according to the Rayleigh-Plesset equation [29, 59]:

$$p(R) = p_{\infty} + \frac{R_0}{R} (p_b - p_{\infty}). \quad (27)$$
<figure><img src="image_2.png" /><figcaption>Figure 1: Problem configuration for a collapsing spherical bubble.</figcaption></figure>

In the following, this pressure initialization is labeled as initial interface equilibrium with $\dot{R}_0 = 0$. We consider cases with both modest and high initial pressure ratios, as shown in Table 1. The water is parameterized by $\gamma_l = 2.35$ and $\pi_\infty = 10^9$ Pa [24, 34, 41, 54, 60, 61].

<table><caption>Table 1: Nominal initial conditions for the cases simulated.</caption><thead><tr><th>Case</th><th>$p_\infty$ [Pa]</th><th>$p_b$ [Pa]</th><th>$p_\infty/p_b$</th></tr></thead><tbody><tr><td>1: Low-pressure-ratio</td><td>$10^5$</td><td>$10^4$</td><td>10</td></tr><tr><td>2: High-pressure-ratio</td><td>$5 \times 10^6$</td><td>3550</td><td>1427</td></tr></tbody></table>

We simulate the flow on a cubical, rectilinear grid with $N_{R_0}$ nodes in each coordinate direction per initial bubble radius near the bubble ($R \le 1.5R_0$); far from the bubble ($R > 1.5R_0$), the grid is stretched nonuniformly to accommodate the large computational domain $L$. To reduce the computational cost, one octant of the domain is computed, with symmetry boundary conditions mimicking the bubble dynamics in neighboring regions. We performed two simulations for each pressure ratio, one without mesh stretching and another with the complete physical domain (no symmetry boundary conditions), and compared them against the simulations presented hereafter to confirm that our results are insensitive to both of these procedures.

When using the WENO5 method, the bubble interface is smeared in the radial direction over a few grid cells. The smearing procedure is commonly employed in models when fifth-order WENO reconstruction is used [23–25, 29, 30, 32, 47, 62–66], as it appears that unphysical oscillations or numerical instabilities can occur without it. The initial interface smearing
procedure we employ involves smearing the volume fraction across the interface using an hyperbolic tangent function [29]

$$ \alpha_g = \frac{1}{2} \left[ 1 - \tanh \left( \frac{R - R_0}{2D} \right) \right], \qquad (28) $$

where D is the characteristic length of the corresponding computational cell; the conservative variables then follow from simple mixture relations, allowing thermodynamic consistency. The physical artifacts associated with this procedure are discussed in section 5.3.2.

In the following, we use the radial bubble-wall evolution to compare the performance of the three different models. We define an effective bubble radius, R, as

$$ R = \left( \frac{3V_b}{4\pi} \right)^{\frac{1}{3}}, \quad \text{where} \quad V_b = \sum_{i=1}^{N} \alpha_{g,i} V_{c,i} \qquad (29) $$

is the total volume of the gas phase, 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 (30) $$

is the nominal total collapse time from its initial (maximum) radius $R_0$ [59]. In our implementation, we compute about 69 × 10³ and 18 × 10³ time steps per $t_c$ for cases 1 and 2, respectively.

# 5. Results

## 5.1. Effect of $K\nabla \cdot \mathbf{u}$ on the 5-equation model

We first reconsider the behavior and influence of the $K\nabla \cdot \mathbf{u}$ term from the 5-equation models on the spherical-bubble-collapse problem using the WENO5 scheme as previously presented by Tiwari et al. [29].

Figure 2 shows that in both pressure-ratio cases, only the model with $K\nabla \cdot \mathbf{u}$ agrees with a semi-analytical solution following the Keller–Miksis equation [67]; a compressible form of the Rayleigh–Plesset equation. The Keller–Miksis equation is based on an asymptotic
<figure><img src="image_3.png" /><figcaption>Figure 2: Radial bubble-wall evolution for (a) $p_\infty/p_b = 10$ with $N_{R_0} = 25$ and (b) $p_\infty/p_b = 1427$ with $N_{R_0} = 50$. Solutions are computed using the 5-equation models with WENO5 as well as the Keller–Miksis equation.</figcaption></figure>

expansion in Mach number which also assumes that the bubble remains spherical. Its use here is predicated on the idea that errors measured relative to it are 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. Furthermore, in this case, the initial interface smearing does not affect the agreement with the Keller–Miksis solution.

The inability of the 5-equation model without $K\nabla \cdot \mathbf{u}$ to represent spherical bubble collapse was previously observed by Tiwari et al. [29], who attributed the better results of the model of Kapila et al. [3] to the enforcement of the second law of thermodynamics. Herein, we seek an alternative explanation in terms of the bubble dynamics themselves.

Figure 3 shows key quantities and $K\nabla \cdot \mathbf{u}$ along a radial coordinate at two instances in time. For $t \approx 0$, the initial interface smearing results in a mixture region at the interface, for which $K \neq 0$, but $K\nabla \cdot \mathbf{u} \approx 0$ because $\nabla \cdot \mathbf{u} \approx 0$. However, during the collapse, $\nabla \cdot \mathbf{u} \neq 0$ and the fluid volume fractions are modified. We see that $K\nabla \cdot \mathbf{u}$ is positive and is larger on the liquid side of the interface. As a result, the liquid volume fraction increases faster, particularly on the liquid side of the interface, than it would without the $K\nabla \cdot \mathbf{u}$ term. This keeps the interface relatively sharp and results in the larger interface velocity observed in
<figure><img src="image_4.png" /><figcaption>Figure 3: Liquid volume-fraction, pressure, and $K\nabla \cdot \mathbf{u}$ for varying radial position $R$ and the case $p_\infty/p_b = 10$ and $N_{R_0} = 25$ using the 5-equation model with $K\nabla \cdot \mathbf{u}$. Times (a) $t \approx 0$ and (b) $t = 0.85t_c$ are shown.</figcaption></figure>

Figure 2.

This behavior can be explained via the bubble pressure evolution. Initially, the pressure is small inside the bubble and increases gradually outside of it. During the collapse, the bubble pressure increases, which reduces the bubble volume due to compression of the gas. This is dynamically coupled to the interface and intensifies the collapse. Further, the pressure always increases in the radial direction. Thus, the gas in the mixture region is more highly compressed on the water side of the bubble interface, and so its volume fraction decreases more rapidly. This process also intensifies the bubble collapse. Note that this second effect is not present into the 5-equation model without $K\nabla \cdot \mathbf{u}$, since this term accounts for expansion and compression in mixture regions.

## 5.2. Comparison of the 5- and 6-equation models

While the 5-equation model with $K\nabla \cdot \mathbf{u}$ can accurately represent spherical bubble dynamics in some cases, it is also often numerically unstable. This is a result of significant compression and expansion near the interface, which can occur during strong shock or expansion waves and cannot be easily treated due to the non-conservative nature of the $K\nabla \cdot \mathbf{u}$ term [33] as discussed in section 1. Here, we consider the 6-equation model as a potential
<figure><img src="image_5.png" /><figcaption>Figure 4: Radial bubble-wall evolution for (a) $p_\infty/p_b = 10$ and (b) $p_\infty/p_b = 1427$. Solutions are computed using WENO5.</figcaption></figure>

solution to this issue; under infinite pressure-relaxation, it theoretically converges to the 5-equation model with $K\nabla \cdot \mathbf{u}$. However, when discretized, the equation sets are different and equivalence has neither been demonstrated for high-order schemes, such as the WENO5 method we consider, nor for the challenging spherical-bubble-collapse test problems. To test our implementation and confirm their convergence to one another, comparison between these methods are presented for shock tube and vacuum problems in Appendix A and B, and consider the collapse of a spherical air bubble in water next.

Simulation results for nominal low- and high-pressure ratios are shown in Figure 4 for both the 5-equation with $K\nabla \cdot \mathbf{u}$ and 6-equation models, following the previous subsections. Both methods agree closely for both cases with the semi-analytic solution of the Keller-Miksis equations [67], which are initialized at equilibrium with $\dot{R}_0 = 0$. The high-pressure-ratio case of Figure 4 (b) also shows the spatial convergence of the models.

Since the 6-equation model has closely matched the 5-equation model with $K\nabla \cdot \mathbf{u}$ for all three of our challenging test cases, we consider it a potential surrogate to the 5-equation model that does not inherit its stability issue. Next, we investigate the behavior of the 6-equation model when solved by numerical schemes of different character and accuracy.
<figure><img src="image_6.png" /><figcaption>Figure 5: Radial bubble-wall evolution for initial interface equilibrium and $p_\infty/p_b = 10$. (a) Results for all schemes and flux limiters as labeled (fixed resolution $N_{R_0} = 25$) and (b) spatial convergence of the numerical methods.</figcaption></figure>

### 5.3. Numerical schemes for the 6-equation model

The 6-equation model can be solved via many different interface-capturing numerical methods. We compute its solution using the methods described in section 3 for a collapsing spherical bubble of varying initial pressure ratio and interface states as a critical assessment of the viability of the numerical schemes for cavitating flows.

#### 5.3.1. Spherical bubble collapse with initial interface equilibrium

We first consider the case of initial interface equilibrium, $\dot{R}_0 = 0$. Figure 5 (a) shows the interface evolution for $p_\infty/p_b = 10$, spatial resolution $N_{R_0} = 25$ and for the FV1, MUSCL2, WENO3 and WENO5 schemes. In addition to the MC [68] slope limiter, the Minmod [43, 69] limiter is implemented for the MUSCL2 scheme (note that the MC limiter is used when not specified for the MUSCL2 scheme). Here, the slope limiters attempt to reduce numerical dissipation of the scheme. As mentioned in section 4, the interface is initially smeared for the WENO5 cases to guarantee numerical stability. However, the MUSCL2 and WENO3 schemes do not require this procedure to remain stable, and thus all interfaces are kept sharp at the grid level in these cases. In Figure 5 (a) we see that the MC slope limiter performs significantly better than the Minmod limiter and similarly for the WENO3 scheme, although
the corresponding results are still less accurate than those of the WEN05 scheme.

To confirm that numerical dissipation is the cause of the discrepancy between the results of the MUSCL2, WEN03 and WEN05 schemes, we consider the spatial convergence of the numerical methods. In Figure 5 (b), this convergence is presented in terms of the discrete L₂ error ε as

$$ \epsilon = \frac{1}{N_t} \sum_{i=0}^{N_t} \frac{\|R(t_i) - R_{KM}(t_i)\|}{R_{KM}(t_i)}, \qquad (31) $$

where N_t is the number of time steps in the temporal window t ∈ [0, 2t_c], and R(t_i) and R_{KM}(t_i) are the bubble radius at time t_i of our simulations and the Keller–Miksis solution, respectively. We see that all methods converge at first order, matching the expected rate for the numerical solution of flows with discontinuities [42, 43, 70]. The WEN05 method has the smallest ε, and so we conclude that for small initial pressure ratios higher-order reconstructions have smaller errors as they suppress numerical diffusion. In this case, the interface smearing procedure we employ for the WEN05 scheme has no apparent consequence on simulation accuracy.

For the flow configurations we consider, the spherical bubble interface is known to be physically stable [59, 71], and so non-spherical interfaces are an artifact of the numerical method; we use this property to assess the performance of the numerical methods. The bubble sphericity is computed as [72]

$$ \Psi = \frac{\pi^{\frac{1}{3}} (6V'_{b})^{\frac{2}{3}}}{A_{b}}, \qquad (32) $$

which is the ratio of the surface area of a sphere with the same volume as the bubble V'_b, to the surface area of the bubble A_b. By this definition, a spherical shape has Ψ = 1 and distorted shapes have Ψ < 1. We define the bubble as the region with α_g ≥ 0.5 and its surface is the isosurface of α_g = 0.5. We compute V'_b and A_b using high-order interpolation of the data.

Sphericity and bubble shape evolution for the small pressure ratio case are shown in Figure 6. We see that the WEN05 scheme maintains sphericity during the entire collapse–rebound process, similar to that observed by Tiwari et al. [29], Whereas the MUSCL2 and
<figure><img src="image_7.png" /><figcaption>Figure 6: Evolution of the bubble sphericity for $p_\infty/p_b = 10$ and $N_{R_0} = 25$. Nominal bubble shapes as represented by $\alpha = 0.5$ isosurfaces are also shown for times $t = 0.7t_c$, $t(R = R_{\min})$, and $2t_c$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Radial bubble-wall evolution for initial interface equilibrium and $p_\infty/p_b = 1427$. (a) Results (fixed resolution $N_{R_0} = 50$) and (b) spatial convergence of the MUSCL2, WENO3 and WENO5 schemes.</figcaption></figure>

WENO3 schemes develop grid-specific artifacts, which are visible beginning at $t = 0.7t_c$; these are presumably due to anisotropic dispersion on the grid with faster propagation of the interface along the Cartesian coordinate directions. By the time of minimum radius $t(R = R_{min})$, the bubble shape is significantly distorted, and at $t = 2t_c$ distortions are still visible.

The radial bubble-wall evolution and convergence results for the larger pressure ratio $p_\infty/p_b = 1427$ are shown in Figure 7. In Figure 7 (a), we only show the Keller–Miksis solution until $t = 1.05t_c$, just after the minimum bubble radius is achieved, since the subsequent rebounds for large pressure ratios are well-known to be physically inaccurate [73]. We see that MUSCL2 is marginally more accurate at predicting the minimum bubble radius and collapse time than the WENO5 method. This seems to be a result of two factors; first, the interface moves more quickly for larger pressure ratios and thus, the MUSCL2 results are less polluted by numerical diffusion over the significantly fewer time steps to reach collapse than were required for the low-pressure-ratio case; second, the initial smearing introduced for the WENO5 method results in an initial diffusion greater than that what ultimately develops during MUSCL2 and WENO3 simulations.

In Figure 7 (b) we plot the observed spatial convergence of the numerical schemes. Here, we only compute $\varepsilon$ over the temporal window $t \in [0, 1.05t_c]$, commensurate with the physical
<figure><img src="image_9.png" /><figcaption>Figure 8: Bubble sphericity evolution for $p_\infty/p_b = 1427$ and $N_{R_0} = 50$. Nominal bubble shapes ($\alpha = 0.5$) are also shown for times $t = 0.5t_c$, $1.02t_c$, and $1.35t_c$.</figcaption></figure>

accuracy of the Keller–Miksis solution over this interval. We again observe approximately first-order convergence for all numerical methods we consider. However, in this case, MUSCL2 has the smallest error $\varepsilon$ and WENO5 the largest. Again, this appears to be a result of the dissipation introduced by the initial smearing procedure used for the WENO5 simulations.

The bubble sphericity and illustrations of the bubble surface are shown in Figure 8. Almost no grid-based artifacts on the bubble surface are visible until $t \approx t_c$ for all numerical methods, at which point $\Psi$ decreases significantly. Compared to the low-pressure-ratio case, the interface evolves more quickly and all methods conserve sphericity for $t \lesssim t_c$. However, after the collapse, significant distortions are visible and $\Psi$ does not reach unity for any of the methods. Furthermore, we see that the WENO5 method results in stronger distortions than the MUSCL2 or WENO3 schemes immediately after the collapse. For larger $t$, the WENO3 scheme
<figure><img src="image_10.png" /><figcaption>Figure 9: Bubble sphericity and associated bubble interface for varying mesh resolution $N_{R_0}$ at time $t = 1.35t_c$ for the $p_\infty/p_b = 1427$ case and MUSCL2 method. For $N_{R_0} = 100$ and 150 the adaptive-mesh-refinement technique of Schmidmayer et al. [38] is used to minimize computational expense.</figcaption></figure>

develops further distortions, eventually reaching similar $\Psi$ values as WENO5 result, whereas the MUSCL2 scheme maintains sphericity after the initial collapse.

For $N_{R_0} = 50$, the minimum radius is about $0.09R_0$, which corresponds to about 4.5 cells per bubble radius in each direction and seemingly leads to a significant amount of anisotropy. Thus, we investigate the effect of mesh resolution on the bubble shape in Figure 9. We see that sphericity indeed improves with increasing the mesh resolution, Note that we only shows results for MUSCL2, but similar behavior is expected for the WENO schemes.

For initial interface equilibrium, we conclude that the WENO5 scheme converges more quickly and can better maintain sphericity when the pressure ratio is relatively small, and thus the maximum interface velocity is much smaller than the Mach number. However, when the pressure ratio is much larger, and so the interface velocity exceeds the Mach number, all the schemes show similar performance, with MUSCL2 only modestly outperforming the others.
<figure><img src="image_11.png" /><figcaption>Figure 10: Wood speed of sound for water-air mixture. Here, $c_{water} = 1625$ m/s and $c_{air} = 350$ m/s.</figcaption></figure>

### 5.3.2. Spherical bubble collapse with initial interface disequilibrium

Lastly, we consider the case of initial interface disequilibrium, and thus $\dot{R}_0(t = 0) \neq 0$. We enforce this by setting the internal and external interface pressures to different values as

$$p = \begin{cases} p_b & \text{for } 0 \le R \le R_0, \\ p_\infty & \text{otherwise.} \end{cases} \qquad (33)$$

This condition represents the discontinuities present, for example, during bubble wall impact. The other initial conditions are identical to previous test cases and thus of section 4.

With the initial interface smearing employed for the WENO5 scheme (or indeed after a sufficient number of time steps for any scheme due to numerical diffusion) the interface has non-negligible thickness, giving rise to a mixture region ($\alpha \neq 0$ or 1). As shown in Figure 10, the Wood speed of sound (12), which is also the speed of sound of the 5-equation model with $K\nabla \cdot \mathbf{u}$, varies in this region and is much less than that of either of the pure phases ($\alpha = 0$ and 1). After the pressure relaxation procedure, the effective speed of sound of the 6-equation model also converges to the Wood speed of sound (see Appendix B).

Figure 11 shows that the smearing procedure employed to keep the WENO5 scheme stable results in an inaccurate solution for the collapse of a bubble in initial pressure disequilibrium. We also see that the MUSCL2 and WENO3 schemes behave similarly when the interface is initially
<figure><img src="image_12.png" /><figcaption>Figure 11: Radial bubble-wall evolution for initial interface disequilibrium, $p_\infty/p_b = 1427$ and $N_{R_0} = 50$. Solutions are computed using the 6-equation model, and the Keller–Miksis result is shown as surrogate truth.</figcaption></figure>

smeared, though this procedure is not required for numerical stability in these cases; for both schemes, the non-smeared cases agree closely with the Keller–Miksis dynamics.

The poor performance associated with the initial interface smearing procedure appears to be due to a wave-trapping phenomenon that results from a lower mixture sound speed, reducing the initial interface velocity. This is illustrated in Figure 12 for the MUSCL2 method. Pressure contours are shown in the $t-R$ space for three degrees of initial smearing (a)–(c). When the interface is not smeared, the pressure waves travel at the pure-phase speed of sound. However, when either the volume fraction or both the volume fraction and mixture pressure are spatially smeared, these waves evolve in a more complex manner due to the reduced sound speeds within the interface mixture region. Pressure waves that escape the mixture region again travel at the liquid speed of sound. The difference between Figure 12 (b) and (c) shows that smearing of the volume fraction $\alpha$ and pressure $p$ both modify the pressure-wave behaviors uniquely, though both ultimately pollute the bubble dynamics.

### 5.4. Interface-sharpening techniques for collapsing spherical bubbles

The numerical dissipation inherent in any interface capturing scheme will eventually smear even initially sharp interfaces. Thus, problems involving multiple interface pressure-disequilibrium events, such as a collapsing ellipsoidal bubble near a wall [74, 75], would
<figure><img src="image_13.png" /><figcaption>Figure 12: Pressure in the time and radial-direction space for initial interface disequilibrium, $p_\infty/p_b = 1427$ and $N_{R_0} = 50$. Solutions are computed using the 6-equation model with the MUSCL2 scheme for three initial configurations: (a) no initial smearing, (b) smearing on only the volume fraction $\alpha$, and (c) smearing on both $\alpha$ and pressure $p$.</figcaption></figure>

benefit from keeping interfaces as sharp as possible. As an attempt to accomplish this, we buttress the 6-equation model and the MUSCL2 and WENO3 schemes with the THINC interface-sharpening method [76] implemented via directional splitting. We note that we observed numerical instabilities when coupling this THINC method to the WENO5 scheme; this seems to be a result of significant interface sharpening when compared to the numerical diffusion of the scheme, which was previously shown to trigger instabilities for this method, and so we do not consider WENO5 herein. Further, we use THINC, rather than anti-diffusion [77] or regularization methods [29], as its conservative property matches the conservative methods we already employ and the others displayed significant numerical instabilities for our methods (particularly for high-pressure-ratio cases). We confirm that our implementation of this method faithfully represents the solutions of the 1D shock tube problem of appendix A and 2D shock-bubble gas-gas interaction problems (such as that of Shyue and Xiao [76] and Deng et al. [78]), and so we can proceed with a faithful comparison for collapsing spherical bubbles.
<figure><img src="image_14.png" /><figcaption>Figure 13: Radial bubble-wall evolution for initial interface equilibrium and numerical methods as labeled; (a) $p_\infty/p_b = 10$ and $N_{R_0} = 25$, (b) $p_\infty/p_b = 1427$ and $N_{R_0} = 50$.</figcaption></figure>

<table><caption>Table 2: Interface thickness $T/R_0$ at specific times and for the cases as labeled. Case 1: $N_{R_0} = 25$ and $t = 2t_c$; case 2: $N_{R_0} = 50$ and $t = 1.35t_c$. Here, $T$ is computed via the number of cells satisfying $0.01 \le \alpha \le 0.99$.</caption><thead><tr><th>Case</th><th>MUSCL2</th><th>MUSCL2 + THINC</th><th>WENO3</th><th>WENO3 + THINC</th></tr></thead><tbody><tr><td>1: $p_\infty/p_b = 10$</td><td>0.32</td><td>0.16</td><td>0.44</td><td>0.16</td></tr><tr><td>2: $p_\infty/p_b = 1427$</td><td>0.36</td><td>0.1</td><td>0.5</td><td>0.08</td></tr></tbody></table>

Figure 13 shows the radial bubble-wall evolution for the low- and high-pressure-ratio interface-pressure equilibrium cases we considered in section 5.3.1. When compared to non-THINC-equipped methods, the THINC results have about 44% larger error $\epsilon$ for the low-pressure-ratio case and 59% smaller $\epsilon$ for the high-pressure-ratio case; however, in all cases the error is already relatively small. Despite having an inconsistent effect on the error, the THINC scheme does keep the bubble interface sharper, as shown in Table 2.

Figure 14 compares bubble shape results for methods with and without THINC. We see that the THINC-coupled method results in significantly less spherical shapes for the low-pressure-ratio cases, though for the high-pressure-ratio cases the shapes are nearly the same. We note that the directional splitting we use could be the source of this discrepancy, though it is unclear if intrinsically multi-dimensional THINC methods will improve this [79]. In
<figure><img src="image_15.png" /><figcaption>Figure 14: Nominal bubble shapes ($\alpha = 0.5$) for times and methods as labeled: (a) $p_\infty/p_b = 10$ and $N_{R_0} = 25$, (b) $p_\infty/p_b = 1427$ and $N_{R_0} = 50$.</figcaption></figure>
general, we see that the THINC method is better behaved when coupled to the MUSCL2,
rather than the WENO3, scheme.

Thus, we conclude that the directionally-split THINC method did not reliably improve
the accuracy of our results, and in most cases disturbed the interface sphericity. As such,
it only offers a partial solution when considering collapsing bubbles with multiple pressure-
disequilibrium events.

**6. Discussion and conclusion**

We analyzed the ability of diffuse-interface models and their associated numerical meth-
ods to represent the collapse and rebound of spherical gas bubbles in a liquid. We confirmed
that the 5-equation model of Allaire et al. [1] and Massoni et al. [2] is unable to accurately
represent a spherical bubble collapse and demonstrated how the additional $K\nabla \cdot u$ term nat-
urally present in the model of Kapila et al. [3] is required to ensure good agreement with the
Keller-Miksis solution [65]. Since the 5-equation model with $K\nabla \cdot u$ is known to produce
instabilities in some numerical experiments [4, 32], we investigated the 6-equation pressure-
disequilibrium model as a potential surrogate. We observed good agreement between these
models for challenging test problems, including a 1D water-air shock tube, a 1D vacuum
developing in a water-air mixture, and the collapse of a 3D spherical bubble. Thus, the
6-equation model is a good candidate to remedy the stability issues of the 5-equation model
with the $K\nabla \cdot u$ source term.

We also considered the behavior and pathologies of the 6-equation model when coupled to
MUSCL and WENO numerical methods for a collapsing spherical bubble. We first analyzed
bubbles at initial interface pressure equilibrium. For this, the bubble interface evolution
of the WENO5-based solution more closely matched the associated Keller-Miksis surrogate-
truth solution than did the MUSCL2 and WENO3 schemes for relatively small pressure ratios.
This was due to the more substantial numerical diffusion intrinsic to the lower-order schemes,
even though the WENO5 scheme required an initially smeared interface to maintain simulation
stability. When the initial pressure ratio was larger, all three methods showed similar results,
quickly converging to the Keller-Miksis solution. Further, we noticed that the relatively
small bubble size at the collapse time resulted in significantly distorted interface shapes. However, these shapes were shown to be more spherical for finer spatial meshes. Thus, an adaptive-mesh-refinement technique would help maintain bubble interface sphericity at the same computational cost as a uniform mesh near the bubble.

When the bubble interface was in initial disequilibrium, we saw that the smearing procedure implemented for the WENO5 method precluded an accurate solution for large pressure ratios. This was a result of the relatively large degree of initial diffusion, which produced a mixture region with a much smaller speed of sound that polluted the dynamics. We also noted that the numerical dissipation inherent in any interface capturing scheme will eventually smear even initially sharp interfaces and, therefore, these schemes would benefit from keeping interfaces as sharp as possible. Interface-sharpening techniques are one way to minimize this dissipation, and we surveyed the THINC method [76] for the same spherical bubble collapse problems. While the THINC method did keep the interfaces sharper, in most cases it further disturbed the interface sphericity. Additionally, we did not observe a consistent increase in simulation accuracy. Thus, further investigation and possibly method improvement is required to maintain surface sharpness while guaranteeing a conservative behavior and numerical stability.

Ultimately, we saw that WENO-based schemes were preferable for bubble dynamics that involve small pressure ratios, and thus slower interface dynamics, and the MUSCL and WENO-based schemes performed similarly for large pressure ratios and thus fast interface speeds. Thus, the WENO5 scheme is generally preferred, except in cases involving interface pressure discontinuities, for which the interface smearing required to keep the scheme stable pollutes the dynamics. As such, the instability of high-order WENO schemes for interface problems warrants future attention.

## Acknowledgments

The authors would like to thank Dr. Mauro Rodriguez, Prof. Eric Johnsen, and Dr. Shahaboddin Alahyari Beig for fruitful discussions. This work was supported by the Office of Naval Research under grant numbers N0014-18-1-2625 and N0014-17-1-2676, and associated
computations utilized the Extreme Science and Engineering Discovery Environment, which
is supported by the National Science Foundation grant number CTS120005.

References

[1] G. Allaire, S. Clerc, S. Kokh, A five-equation model for the simulation of interfaces between compressible fluids, J. Comp. Phys. 181 (2002) 577-616.

[2] J. Massoni, R. Saurel, B. Nkonga, R. Abgrall, 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 45 (2002) 1287–1307.

[3] A. Kapila, R. Menikoff, J. Bdzil, S. Son, D. Stewart, Two-phase modeling of DDT in granular materials: Reduced equations, Phys. Fluids 13 (2001) 3002–3024.

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

[5] C. E. Brennen, Cavitation in medicine, Interface Focus 5 (2015) 20150022.

[6] B. Dollet, P. Marmottant, V. Garbin, Bubble dynamics in soft and biological matter, Annu. Rev. Fluid Mech. 51 (2019) 331–355.

[7] J. B. Estrada, C. Barajas, D. L. Henann, E. Johnsen, C. Franck, High strain-rate soft material characterization via inertial cavitation, Journal of the Mechanics and Physics of Solids 112 (2018) 291–317.

[8] S. Pan, S. Adami, X. Hu, N. A. Adams, Phenomenology of bubble-collapse-driven penetration of biomaterial-surrogate liquid-liquid interfaces, Phys. Rev. Fluids 3 (2018) 114005.
[9] R. Oguri, K. Ando, Cavitation bubble nucleation induced by shock-bubble interaction in a gelatin gel, Phys. Fluids 30 (2018) 051904.

[10] C. K. Turangan, G. J. Ball, A. R. Jamaluddin, T. G. Leighton, Numerical studies of cavitation erosion on an elastic-plastic material caused by shock-induced bubble collapse, Proc. Royal Soc. A 473 (2017) 20170315.

[11] C. Barajas, E. Johnsen, The effects of heat and mass diffusion on freely oscillating bubbles in a viscoelastic, tissue-like medium, J. Acoust. Soc. Am. 141 (2017) 908-918.

[12] R. Gaudron, M. Warnez, E. Johnsen, Bubble dynamics in a viscoelastic medium with nonlinear elasticity, J. Fluid Mech. 766 (2015) 54-75.

[13] G. Ball, B. Howell, T. Leighton, M. Schofield, Shock-induced collapse of a cylindrical air cavity in water: a free-lagrange simulation, Shock Waves 10 (2000) 265-276.

[14] C. Turangan, A. Jamaluddin, G. Ball, T. Leighton, Free-lagrange simulations of the expansion and jetting collapse of air bubbles in water, J. Fluid Mech. 598 (2008) 1-25.

[15] J. Glimm, X. Li, Y. Liu, Z. Xu, N. Zhao, Conservative front tracking with improved accuracy, SIAM J. Num. Anal. 41 (2003) 1926-1947.

[16] H. Terashima, G. Tryggvason, A front-tracking/ghost-fluid method for fluid interfaces in compressible flows, J. Comp. Phys. 228 (2009) 4012-4037.

[17] S. Pan, L. Han, X. Hu, N. A. Adams, A conservative interface-interaction method for compressible multi-material flows, J. Comp. Phys. 371 (2018) 870-895.

[18] X. Hu, B. Khoo, N. Adams, F. Huang, A conservative interface method for compressible flows, J. Comp. Phys. 219 (2006) 553-578.

[19] L. Han, X. Hu, N. Adams, Adaptive multi-resolution method for compressible multi-phase flows with sharp interface model and pyramid data structure, J. Comp. Phys. 262 (2014) 131-152.
[20] C.-H. Chang, X. Deng, T. Theofanous, Direct numerical simulation of interfacial instabilities: a consistent, conservative, all-speed, sharp-interface method, J. Comp. Phys. 242 (2013) 946–990.

[21] K. Shahbazi, High-order finite difference scheme for compressible multi-component flow computations, Comp. Fluids 190 (2019) 425–439.

[22] S. Mirjalili, C. B. Ivey, A. Mani, Comparison between the diffuse interface and volume of fluid methods for simulating two-phase flows, Int. J. Multiphase Flow 116 (2019) 221–238.

[23] V. Coralic, T. Colonius, Finite-volume WENO scheme for viscous compressible multi-component flows, J. Comp. Phys. 274 (2014) 95–121.

[24] S. A. Beig, E. Johnsen, Maintaining interface equilibrium conditions in compressible multiphase flows using interface capturing, J. Comp. Phys. 302 (2015) 548–566.

[25] J.-C. Veilleux, K. Maeda, T. Colonius, J. E. Shepherd, Transient cavitation in pre-filled syringes during autoinjector actuation, Int. Symp. Cav. (2018).

[26] G. Xiang, B. Wang, Numerical investigation on the interaction of planar shock wave with an initial ellipsoidal bubble in liquid medium, AIP Adv. 8 (2018) 075128.

[27] J. Meng, T. Colonius, Numerical simulations of the early stages of high-speed droplet breakup, Shock Waves (2014) 1–16.

[28] J. Meng, T. Colonius, Numerical simulation of the aerobreakup of a water droplet, J. Fluid Mech. 835 (2018) 1108–1135.

[29] A. Tiwari, J. Freund, C. Pantano, A diffuse interface model with immiscibility preservation, J. Comp. Phys. 252 (2013) 290–309.

[30] U. Rasthofer, F. Wermelinger, P. Hadiijdoukas, P. Koumoutsakos, Large scale simulation of cloud cavitation collapse, Procedia Comput. Sci. 108 (2017) 1763–1772.
[31] M. Baer, J. Nunziato, A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials, Int. J. of Multiphase Flows 12 (1986) 861–889.

[32] S. A. Beig, B. Aboulhasanzadeh, E. Johnsen, Temperatures produced by inertially collapsing bubbles near rigid surfaces, J. Fluid Mech. 852 (2018) 105–125.

[33] G. Dal Maso, P. Le Floch, F. Murat, Definition and weak stability of nonconservative products, J. Math. Pures App. 74 (1995) 483–548.

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

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

[36] K. Schmidmayer, F. Petitpas, E. Daniel, N. Favrie, S. Gavrilyuk, A model and numerical method for compressible flows with capillary effects, J. Comp. Phys. 334 (2017) 468–496.

[37] K. Schmidmayer, Simulation de l'atomisation d'une goutte par un écoulement à grande vitesse, Ph.D. thesis, Aix-Marseille, 2017.

[38] K. Schmidmayer, F. Petitpas, E. Daniel, Adaptive Mesh Refinement algorithm based on dual trees for cells and faces for multiphase compressible flows, J. Comp. Phys. 388 (2019) 252–278.

[39] N. Favrie, S. L. Gavrilyuk, R. Saurel, Solid-fluid diffuse interface model in cases of extreme deformations, J. Comp. Phys. 228 (2009) 6037–6077.

[40] S. Ndanou, N. Favrie, S. Gavrilyuk, Multi-solid and multi-fluid diffuse interface model: Applications to dynamic fracture and fragmentation, J. Comp. Phys. 295 (2015) 523–555.
[41] M. Pelanti, K.-M. Shyue, A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves, J. Comp. Phys. 259 (2014) 331–357.

[42] S. K. Godunov, A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics, Mat. Sb. 89 (1959) 271–306.

[43] E. Toro, Riemann solvers and numerical methods for fluid dynamics, Springer Verlag, Berlin, 1997.

[44] B. Van Leer, Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method, J. Comp. Phys. 32 (1979) 101–136.

[45] K. Schmidmayer, F. Petitpas, S. Le Martelot, E. Daniel, ECOGEN: An open-source tool for multiphase, compressible, multiphysics flows, Submitted (2019).

[46] G.-S. Jiang, C.-W. Shu, Efficient implementation of weighted ENO schemes, J. Comp. Phys. 126 (1996) 202–228.

[47] R. K. Shukla, C. Pantano, J. B. Freund, An interface capturing method for the simulation of multi-phase compressible flows, J. Comp. Phys. 229 (2010) 7411–7439.

[48] P. Zhang, S. C. Wong, C.-W. Shu, A weighted essentially non-oscillatory numerical scheme for a multi-class traffic flow model on an inhomogeneous highway, J. Comp. Phys. 212 (2006) 739–756.

[49] X. Hu, Q. Wang, N. A. Adams, An adaptive central-upwind weighted essentially non-oscillatory scheme, J. Comp. Phys. 229 (2010) 8952–8965.

[50] X. Hu, N. A. Adams, Scale separation for implicit large eddy simulation, J. Comp. Phys. 230 (2011) 7240–7249.

[51] L. Fu, X. Y. Hu, N. A. Adams, A family of high-order targeted ENO schemes for compressible-fluid simulations, J. Comp. Phys. 305 (2016) 333–359.
[52] B. Van Leer, Towards the ultimate conservative difference scheme III. Upstream-centered finite-difference schemes for ideal compressible flow, J. Comp. Phys. 23 (1977) 263-275.

[53] S. Gottlieb, C.-W. Shu, Total variation diminishing Runge-Kutta schemes, Math. Comput. 67 (1998) 73-85.

[54] O. Le Métayer, J. Massoni, R. Saurel, Elaborating equations of state of a liquid and its vapor for two-phase flow models, Int. J. Therm. Sci. 43 (2004) 265-276.

[55] A. Wood, A textbook of sound, G. Bell and Sons LTD, London (1930).

[56] G. B. Wallis, One-dimensional two-phase flow, McGraw-Hill, 1969.

[57] K. Schmidmayer, A. Marty, F. Petitpas, E. Daniel, ECOGEN, an open source tool dedicated to multiphase compressible multiphysics flows, Int. Conf. Appl. Aero. (2018).

[58] M. Rodriguez, E. Johnsen, A high-order accurate, five-equations compressible multi-phase approach for viscoelastic fluids and solids with relaxation and elasticity, J. Comp. Phys. (2018).

[59] C. E. Brennen, Cavitation and bubble dynamics, Oxford University Press, 1995.

[60] O. Le Métayer, J. Massoni, R. Saurel, Modeling evaporation fronts with reactive Riemann solvers, J. Comp. Phys. 205 (2005) 567-610.

[61] R. Saurel, F. Petitpas, R. Abgrall, Modelling phase transition in metastable liquids: application to cavitating and flashing flows, J. Fluid Mech. 607 (2008) 313-350.

[62] E. Johnsen, T. Colonius, Implementation of WENO schemes in compressible multicomponent flow problems, J. Comp. Phys. 219 (2006) 715-732.

[63] E. Johnsen, T. Colonius, Numerical simulations of non-spherical bubble collapse, J. Fluid Mech. 629 (2009) 231-262.
[64] J. H. J. Niederhaus, J. A. Greenough, J. G. Oakley, D. Ranjan, M. Anderson, R. Bonazza, A computational parameter study for the three-dimensional shock-bubble interaction, J. Fluid Mech. 594 (2008) 85-124.

[65] A. Tiwari, C. Pantano, J. Freund, Growth-and-collapse dynamics of small bubble clusters near a wall, J. Fluid Mech. 775 (2015) 1-23.

[66] G. Xiang, B. Wang, Numerical study of a planar shock interacting with a cylindrical water column embedded with an air cavity, J. Fluid Mech. 825 (2017) 825-852.

[67] J. B. Keller, M. Miksis, Bubble oscillations of large amplitude, J. Acoust. Soc. Am. 68 (1980) 628-633.

[68] C. F. Gammie, J. C. McKinney, G. Tóth, HARM: A numerical scheme for general relativistic magnetohydrodynamics, Astrophys. J. 589 (2003) 444.

[69] P. K. Sweby, High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM J. Num. Anal. 21 (1984) 995-1011.

[70] B. Van Leer, Towards the ultimate conservative difference scheme. II. Monotonicity and conservation combined in a second-order scheme, J. Comp. Phys. 14 (1974) 361-370.

[71] D. Frost, B. Sturtevant, Effects of ambient pressure on the instability of a liquid boiling explosively at the superheat limit, J. Heat Trans. 108 (1986) 418-424.

[72] H. Wadell, Volume, shape, and roundness of quartz particles, J. Geol. 43 (1935) 250-280.

[73] D. Fuster, C. Dopazo, G. Hauke, Liquid compressibility effects during the collapse of a single cavitating bubble, J. Acoust. Soc. Am. 129 (2011) 122-131.

[74] Y. A. Pishchalnikov, W. Behnke-Parks, K. Maeda, T. Colonius, M. Mellema, M. Hopcroft, A. Luong, S. Wiener, M. L. Stoller, T. Kenny, D. J. Laser, Experimental observations and numerical modeling of lipid-shell microbubbles with calcium-adhering
moieties for minimally-invasive treatment of urinary stones, in: Proc. Meet. Acoust.,
volume 35, 2018, p. 020008.

[75] Y. A. Pishchalnikov, W. M. Behnke-Parks, K. Schmidmayer, K. Maeda, T. Colonius, T. W. Kenny, D. J. Laser, High-speed video microscopy and numerical modeling of bubble dynamics near a surface of urinary stone, J. Acoust. Soc. Am. 146 (2019) 516–531.

[76] K. Shyue, F. Xiao, An Eulerian interface sharpening algorithm for compressible two-phase flow: The algebraic THINC approach, J. Comp. Phys. 268 (2014) 326–354.

[77] K. K. So, X. Y. Hu, N. A. Adams, Anti-diffusion interface sharpening technique for two-phase compressible flow simulations, J. Comp. Phys. 231 (2012) 4304–4323.

[78] X. Deng, S. Inaba, B. Xie, K.-M. Shyue, F. Xiao, High fidelity discontinuity-resolving reconstruction for compressible multiphase flows with moving interfaces, J. Comp. Phys. 371 (2018) 945–966.

[79] S. Ii, K. Sugiyama, S. Takeuchi, S. Takagi, Y. Matsumoto, F. Xiao, An interface capturing method with a continuous function: The THINC method with multi-dimensional reconstruction, J. Comp. Phys. 231 (2012) 2328–2358.

[80] A. Murrone, H. Guillard, A five equation reduced model for compressible two phase flow problems, J. Comp. Phys. 202 (2005) 664–698.
<figure><img src="image_16.png" /><figcaption>Figure 15: Shock tube problem setup and initial and boundary conditions.</figcaption></figure>

## Appendix A Water-air shock tube

We consider a water-air shock tube, whose initial configuration is shown in Figure 15 [4, 24, 38]. The domain has length $L$, the initial discontinuity is located at $L^* = L/7$ and $10^3$ nodes are used. Here, the water has stiffened-gas parameters $\gamma_l = 4.4$ and $\pi_\infty = 6 \times 10^8$ Pa [4, 36, 80].

We simulate the flow in the shock tube using both the 5-equation with $K\nabla \cdot \mathbf{u}$ and 6-equation models and the WENO5 numerical scheme. A uniform and one-dimensional mesh of $10^3$ nodes is used. Results for the primitive variables at $t = 241$ µs are shown in Figure 16. A rightward shock wave propagates into the air, followed by a contact discontinuity, observable in (a) and (b), that delimits the interface between the two phases; left-going expansion waves propagate into the water. We observe good agreement between the numerical implementations of both models and the exact solution. Indeed, differences can only be seen at the tail of the expansion waves and near the contact discontinuity. Note that these differences diminish with increasing resolution.

## Appendix B Vacuum generation into a water-air mixture

We consider a vacuum generation into a water-air mixture [4, 41]. The problem setup is shown in Figure 17; there is a uniform initial pressure $p = 10^5$ Pa and densities $\rho_l = 10^3$ kg/m$^3$ and $\rho_g = 1$ kg/m$^3$, and a flow is generated by the initial discontinuity in velocity. Again, $10^3$ nodes are used and the water has stiffened-gas parameters $\gamma_l = 4.4$ and $\pi_\infty = 6 \times 10^8$ Pa.

Figure 18 shows the results of the primitive variables at $t = 1.85$ ms for the vacuum
<figure><img src="image_17.png" /><figcaption>Figure 16: Water-air shock-tube interaction problem at $t = 241$ µs. Numerical and exact solutions are as labeled above.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 17: Problem setup for a vacuum generation into a water-air mixture.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 18: Cavitating water-air mixture problem at $t = 1.85$ ms. Numerical and exact solutions are as labeled above.</figcaption></figure>
problem using the same methods and computational parameterization as Appendix A. The discontinuity in velocity generates left- and right-going expansion waves, and thus generates a $p = 0$ vacuum in the center of the domain. Mixture compressibility ensures that the water volume fraction, and thus the mixture density, decreases in the vacuum region. We observe good agreement between the numerical simulations and exact solution. However, the 6-equation model generally performs better, with no pressure oscillations at the head of the expansion waves.