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

TL;DR

This paper evaluates various multicomponent flow models and interface capturing schemes for accurately simulating spherical bubble collapse, highlighting the advantages of pressure-disequilibrium models over traditional approaches.

Contribution

It compares five- and six-equation models and interface-capturing methods, demonstrating the superiority of pressure-disequilibrium models in spherical bubble dynamics simulations.

Findings

Traditional 5-equation model is inadequate for bubble collapse.

Augmented 5-equation model better captures bubble dynamics.

Pressure-disequilibrium models show versatility in 3D bubble collapse simulations.

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 general, 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], the 5-equation model of Kapila et al. [2], and the 6-equation model of Saurel et al. [3] 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 of Allaire et al. [1] for spherical bubble collapse problems and explain the corresponding advantages of the augmented model of Kapila et al. [2] 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 pressure-disequilibrium models. Lastly, the performance of pressure disequilibrium model for representing a three-dimensional spherical bubble collapse for different bubble interior/exterior pressure ratios is evaluated for different numerical methods. Pathologies associated with each factor and their origins are identified and discussed.

Tables

Table 1. Table 1: Nominal initial conditions for the cases simulated.

Case	$p_{\infty }$ [Pa]	$p_b$ [Pa]	$p_{\infty }/p_b$
1: Low-pressure-ratio	${10}^5$	${10}^4$	10
2: High-pressure-ratio	$5\times {10}^6$	$3550$	1427

Table 2. Table 2: Interface thickness T / R 0 𝑇 subscript 𝑅 0 T/R_{0} at specific times and for the cases as labeled. Case 1 : N R 0 = 25 subscript 𝑁 subscript 𝑅 0 25 N_{R_{0}}=25 and t = 2 ​ t c 𝑡 2 subscript 𝑡 𝑐 t=2t_{c} ; case 2 : N R 0 = 50 subscript 𝑁 subscript 𝑅 0 50 N_{R_{0}}=50 and t = 1.35 ​ t c 𝑡 1.35 subscript 𝑡 𝑐 t=1.35t_{c} Here, T 𝑇 T is computed via the number of cells satisfying 0.01 ⩽ α ⩽ 0.99 0.01 𝛼 0.99 0.01\leqslant\alpha\leqslant 0.99 .

Case	MUSCL2	MUSCL2 + THINC	WENO3	WENO3 + THINC
1: $p_{\infty }/p_b=10$	0.32	0.16	0.44	0.16
2: $p_{\infty }/p_b=1427$	0.36	0.1	0.5	0.08

Equations

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

$$\displaystyle\mathbf{q}=\left[\begin{array}[]{c}\alpha_{1}\\ \alpha_{1}\rho_{1}\\ \alpha_{2}\rho_{2}\\ \rho\mathbf{u}\\ \rho E\end{array}\right],\quad\mathbf{F}=\left[\begin{array}[]{c}\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}\\ \left(\rho E+p\right)\mathbf{u}\end{array}\right],\quad\mathbf{h}=\left[\begin{array}[]{c}-\alpha_{1}\\ 0\\ 0\\ \mathbf{0}\\ 0\end{array}\right],\quad\mathbf{r}=\left[\begin{array}[]{c}0\\ 0\\ 0\\ \mathbf{0}\\ 0\end{array}\right],$$

$$E=e+\frac{1}{2}\Vert \mathbf{u}{\Vert }^2,$$

$$e=\sum _{k=1}^2{Y}_k{e}_k\left({\rho }_k,p\right).$$

$$Y_k=\frac{{\alpha }_k{\rho }_k}{\rho }.$$

$$p_g=({\gamma }_g-1){\rho }_g{e}_g,$$

$$p_l=({\gamma }_l-1){\rho }_l{e}_l-{\gamma }_l{\pi }_{\infty },$$

$$\rho =\sum _{k=1}^2{\alpha }_k{\rho }_k\text{and}p=\sum _{k=1}^2{\alpha }_k{p}_k,$$

$$\rho c^2=\sum _{k=1}^2\frac{{\alpha }_k{\rho }_k{c}_k^2}{\beta \left({\gamma }_k-1\right)},\beta =\sum _{k=1}^2\frac{{\alpha }_k}{{\gamma }_k-1},$$

$$\displaystyle\mathbf{q}=\left[\begin{array}[]{c}\alpha_{1}\\ \alpha_{1}\rho_{1}\\ \alpha_{2}\rho_{2}\\ \rho\mathbf{u}\\ \rho E\end{array}\right],\quad\mathbf{F}=\left[\begin{array}[]{c}\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}\\ \left(\rho E+p\right)\mathbf{u}\end{array}\right],\quad\mathbf{h}=\left[\begin{array}[]{c}-\alpha_{1}\\ 0\\ 0\\ \mathbf{0}\\ 0\end{array}\right],\quad\mathbf{r}=\left[\begin{array}[]{c}K\nabla\cdot\mathbf{u}\\ 0\\ 0\\ \mathbf{0}\\ 0\end{array}\right],$$

$$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}},$$

$$\frac{1}{\rho c^2}=\sum _{k=1}^2\frac{{\alpha }_k}{{\rho }_k{c}_k^2},$$

$$\displaystyle\mathbf{q}=\left[\begin{array}[]{c}\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{array}\right],\quad\mathbf{F}=\left[\begin{array}[]{c}\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{array}\right],\quad\mathbf{h}=\left[\begin{array}[]{c}-\alpha_{1}\\ 0\\ 0\\ \mathbf{0}\\ \alpha_{1}p_{1}\\ \alpha_{2}p_{2}\end{array}\right],\quad\mathbf{r}=\left[\begin{array}[]{c}\mu\delta p\\ 0\\ 0\\ \mathbf{0}\\ -\mu p_{I}\delta p\\ \mu p_{I}\delta p\end{array}\right],$$

$$p_I=\frac{z_2{p}_1+z_1{p}_2}{z_1+z_2},$$

$$\delta p=p_1-p_2,$$

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

$$c^2=\sum _{k=1}^2{Y}_k{c}_k^2.$$

$${\mathbf{q}}_i^{n+1}={\mathbf{q}}_i^n-\frac{\Delta t}{V_i}\left(\sum _{s=1}^N{A}_s{\mathbf{F}}_s^{\star }\cdot {\mathbf{n}}_s+\mathbf{h}\left({\mathbf{q}}_i^n\right)\sum _{s=1}^N{A}_s{\mathbf{u}}_s^{\star }\cdot {\mathbf{n}}_s\right),$$

$${\mathbf{q}}_i^{n+\frac{1}{2}}$$

$${\mathbf{q}}_i^{n+1}$$

$${\mathbf{q}}_i^{(1)}$$

$${\mathbf{q}}_i^{(2)}$$

$${\mathbf{q}}_i^{n+1}$$

$$p(R)=p_{\infty }+\frac{R_0}{R}\left(p_b-p_{\infty }\right).$$

$${\alpha }_g=\frac{1}{2}\left[1-\text{tanh}\left(\frac{R-R_0}{2D}\right)\right],$$

$$R={\left(\frac{3V_b}{4\pi }\right)}^{\frac{1}{3}},\text{where}{V}_b=\sum _{i=1}^N{\alpha }_{g,i}{V}_{c,i}$$

$$t_c=0.915R_0\sqrt{\frac{{\rho }_l}{p_{\infty }}}$$

$$\epsilon =\frac{1}{N_t}\sum _{i=0}^{N_t}\frac{\Vert R(t_i)-R_{KM}(t_i)\Vert }{R_{KM}(t_i)},$$

$$\Psi =\frac{{\pi }^{\frac{1}{3}}{\left(6V_b^{\prime }\right)}^{\frac{2}{3}}}{A_b},$$

$$p=\{\begin{array}{ll}{p}_b& \text{mbox}for0\le R\le R_0,\\ p_{\infty }& \text{mbox}otherwise.\end{array}$$

Full text

††journal: J. Comp. Phys.

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 [4, 5, 6, 7, 8] and polymeric coatings and biofouling in industry [9]. Preliminary studies have shown that bubble dynamics are sensitive to the properties of these materials [10, 11], 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.

Diffuse interface (interface-capturing) methods appear to be well-suited for this problem when compared to other interface tracking and capturing schemes [12, 13]. Such methods combine a multicomponent flow model with shock-capturing finite volume methods. 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. Herein, 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] has been widely used and can faithfully represent shock-induced collapses [14, 15, 16, 17] and droplet atomization [18, 19]. Unfortunately, this model cannot predict the collapse time and minimum radius of the Rayleigh collapse problem [20, 21]. This problem can be averted via the thermodynamically consistent model of Kapila et al. [2] [20, 21], which includes a term ($K\nabla\cdot\mathbf{u}$) in the volume-fraction evolution equation to represent compressibility in mixture regions. Unfortunately, this additional term leads to numerical instabilities during strong compression and expansion near the interface [3, 22]. Instead, we propose using a pressure-disequilibrium model [3], 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. [2] under mesh refinement, and while it has been utilized for cavitating flows [3, 23], detonating flows [24], surface-tension driven flows [25], droplet atomization [26, 27], and fracture and fragmentation in ductile materials [28, 29], it has not been applied to bubble dynamics or particularly to collapsing bubbles.

The multicomponent flow models are solved using shock-capturing finite-volume schemes and Riemann solvers for fluxes [30, 31]. High-order spatial reconstructions, such as MUSCL [32, 31, 33] and WENO [34, 35, 14, 20], are often used, along with their variants WENO-Z [36], WENO-CU6 [37, 38], and TENO [39]. Herein, we will consider MUSCL and the WENO of Jiang and Shu [34], coupled with the HLLC approximate Riemann solver [31, 14, 3] as standard approaches for solving the multicomponent flow equations. Following usual procedure, these are coupled to total-variation-diminishing time integrators as an attempt to suppress spurious oscillations at material interfaces under refinement [30, 40, 31, 41].

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

[TABLE]

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]

We first consider the mechanical-equilibrium model of Allaire et al. [1], which we call the 5-equation model. For a two-phase flow, we have

[TABLE]

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

[TABLE]

where $e$ is the mixture specific internal energy

[TABLE]

In (24), $e_{k}$ is defined via an equation of state and $Y_{k}$ are the mass fractions

[TABLE]

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

[TABLE]

where $\gamma_{g}=1.4$, and the liquid is modeled by the stiffened-gas equation of state

[TABLE]

where $\gamma_{l}$ and $\pi_{\infty}$ are case-specific model parameters [42]. The mixture quantities are

[TABLE]

and

[TABLE]

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, 25].

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

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

[TABLE]

where only $\mathbf{r}$ is different from (22). Here, $K$ is

[TABLE]

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

[TABLE]

which is also the Wood speed of sound [43, 44].

2.3 Pressure-disequilibrium model of Saurel et al. [3]

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

[TABLE]

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

[TABLE]

where ${z_{k}}={\rho_{k}}{c_{k}}$ is the acoustic impedance of the phase $k$, and

[TABLE]

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

[TABLE]

We note that (80) 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. [3]).

The mixture speed of sound is defined according to

[TABLE]

After applying the infinite pressure-relaxation procedure detailed in section 3.2, the effective mixture speed of sound matches (52). We will discuss the influence of sound speed for interface problems in section 5.3.2.

3 Numerical methods

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

[TABLE]

where $n$ is the time-step index. Here, we label this basic first-order-accurate finite-volume scheme as FV1. At the volume–volume interfaces, the associated Riemann problem is computed using the HLLC approximate solver [14, 3, 31], giving the flux tensor and flow-velocity vector $\mathbf{F}^{\star}_{s}$ and $\mathbf{u}^{\star}_{s}$, respectively. The solution of (82) is restricted by the usual CFL criterion.

3.1 Spatial and time reconstruction

Herein, we utilize both MUSCL and WENO spatial reconstructions. We use the second-order-accurate MUSCL scheme of Schmidmayer et al. [33] (labeled here as MUSCL2) with two-step time integration

[TABLE]

where $\mathbf{L}$ is the numerically approximated fluxes and non-conservative terms. The first step is a prediction for the second step and the usual piece-wise linear MUSCL reconstruction [31] is used on the primitive variables. The monotonized central (MC) [40] 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 [3, 23, 24, 25, 33, 27, 45, 29].

The WENO scheme is either third- or fifth-order accurate (labeled here as WENO3 and WENO5, respectively) and reconstructs the primitive variables [14]. In this case, the time derivative is computed via the third-order TVD Runge–Kutta algorithm [41]

[TABLE]

This method has previously been implemented for the pressure-equilibrium models of Allaire et al. [1] [14, 15, 19, 18, 16] and of Kapila et al. [2] [20, 21, 46]; here, we also utilize it for the pressure-disequilibrium model of Saurel et al. [3].

3.2 Pressure-relaxation procedure

The pressure-disequilibrium model (77) requires pressure-relaxation to converge to a single, equilibrium pressure. We use the infinite-relaxation procedure of Saurel et al. [3]. At each time step it solves the non-relaxed, hyperbolic equations ($\mu\to 0$) using (82), then relaxes the disequilibrium pressures for $\mu\to+\infty$. The latter 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. [2] (50). When multi-stage time integration is used, the relaxation procedure is 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] and Kapila et al. [2] 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 [47, 20]:

[TABLE]

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}\text{,}\mathrm{Pa}$$ [42, 48, 49, 23, 50, 15].

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\leqslant 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 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 few grid cells. The smearing procedure is commonly employed in multicomponent models when fifth-order WENO reconstruction is used [51, 52, 53, 35, 20, 14, 54, 15, 21, 55, 16, 22], 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 [20]

[TABLE]

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 multicomponent models. We define an effective bubble radius, $R$, as

[TABLE]

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

[TABLE]

is the nominal total collapse time from its initial (maximum) radius $R_{0}$ [47]. In our implementation, we compute about $69\times 10^{3}$ and $18\times 10^{3}$ 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. [20].

Figure 2 shows that in both pressure-ratio cases, only the model with $K\nabla\cdot\mathbf{u}$ agrees with a nominal exact solution following the Keller–Miksis equation [56]; a compressible form of the Rayleigh–Plesset equation. In this case, the initial interface smearing does not affect this agreement.

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. [20], who attributed the better results to enforcement of the second law of thermodynamics [2]. We seek here an alternative explanation in terms of the dynamics.

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 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. Here, we consider the 6-equation model as a potential 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 analytic solution of the Keller–Miksis equations [56], 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.

5.3 Numerical schemes for the 6-equation model

The 6-equation model can be solved via a number of 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 [57] slope limiter, the Minmod [58, 31] 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 WENO5 scheme.

To confirm that numerical dissipation is the cause of the discrepancy between the results of the MUSCL2, WENO3 and WENO5 schemes, we consider the spatial convergence of the numerical methods. In Figure 5 (b), this convergence is presented in terms of the discrete $L_{2}$ error $\varepsilon$ as

[TABLE]

where $N_{t}$ is the number of time steps in the temporal window $t\in[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 [30, 59, 31]. The WENO5 method has the smallest $\varepsilon$, 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 WENO5 scheme has no apparent consequence on simulation accuracy.

For the flow configurations we consider, the spherical bubble interface is known to be physically stable [60, 47], 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 [61]

[TABLE]

which is the ratio of the surface area of a sphere with the same volume as the bubble $V_{b}^{\prime}$, to the surface area of the bubble $A_{b}$. By this definition, a spherical shape has $\Psi=1$ and distorted shapes have $\Psi<1$. We define the bubble as the region with $\alpha_{g}\geq 0.5$ and its surface is the isosurface of $\alpha_{g}=0.5$. We compute $V_{b}^{\prime}$ 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 WENO5 scheme maintains sphericity during the entire collapse–rebound process. The MUSCL2 and 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_{\text{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 [62]. 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 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 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.

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

[TABLE]

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 (52), 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 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 [63, 64], would benefit from keeping interfaces as sharp as possible.

We buttress the 6-equation model and the MUSCL2 and WENO3 schemes with the THINC interface-sharpening method [65] as an illustration of the behavior of interface-sharpening methods for spherical bubble dynamics. We use THINC, rather than anti-diffusion [66] or regularization methods [20], 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 note that our implementation of the THINC method faithfully represented the solutions to the 1D shock tube problem of Appendix A and 2D shock–bubble gas–gas interaction problems (such as that of Shyue and Xiao [65] and Deng et al. [67]), and so we can proceed with a faithful comparison for collapsing spherical bubbles. Further, we noticed that the THINC method has numerical instabilities when coupled to the WENO5 scheme; this seems to be a result of significant interface sharpening, which were previously shown to trigger instabilities for this method, and so we do not consider WENO5 herein.

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 $\varepsilon$ for the low-pressure-ratio case and $59\%$ smaller $\varepsilon$ 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 shapes with and without THINC. We see that the THINC 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. In 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 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 methods 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] is unable to accurately represent a spherical bubble collapse and demonstrated how the additional $K\nabla\cdot\mathbf{u}$ term introduced by Kapila et al. [2] is required to ensure good agreement with the Keller–Miksis solution [54]. Since the 5-equation model with $K\nabla\cdot\mathbf{u}$ is known to produce instabilities in some numerical experiments [3, 22], 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\mathbf{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, and despite the fact that 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 be helpful for maintaining 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 [65] 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 are 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.

Appendix A Water-air shock tube

We consider a water-air shock tube, whose initial configuration is shown in Figure 15 [3, 27, 15]. 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}\text{,}\mathrm{Pa}$$ [68, 3, 25].

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\text{,}\mathrm{\SIUnitSymbolMicro 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 [3, 50]. The problem setup is shown in Figure 17; there is a uniform initial pressure $p=$10^{5}\text{,}\mathrm{Pa} and densities $\rho_{l}=$10^{3}\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3} and $\rho_{g}=$1\text{,}\mathrm{kg}\text{,}{\mathrm{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}\text{\,}\mathrm{Pa}.

Figure 18 shows the results of the primitive variables at $t=$1.85\text{,}\mathrm{ms}$$ for the vacuum 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.