Kevin Schmidmayer
email: kevinsch@caltech.edu

Spencer H Bryngelson
email: spencer@caltech.edu

Tim Colonius
email: colonius@caltech.edu

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

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

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 5equation 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 numerical methods. Pathologies associated with each factor and their origins are identified and discussed.

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 [START_REF] Brennen | Cavitation in medicine[END_REF][START_REF] Dollet | Bubble dynamics in soft and biological matter[END_REF][START_REF] Estrada | High strain-rate soft material characterization via inertial cavitation[END_REF][START_REF] Pan | Phenomenology of bubble-collapse-driven penetration of biomaterial-surrogate liquid-liquid interfaces[END_REF][START_REF] Oguri | Cavitation bubble nucleation induced by shock-bubble interaction in a gelatin gel[END_REF] and polymeric coatings and biofouling in industry [START_REF] Turangan | Numerical studies of cavitation erosion on an elastic-plastic material caused by shock-induced bubble collapse[END_REF]. Preliminary studies have shown that bubble dynamics are sensitive to the properties of these materials [START_REF] Barajas | The effects of heat and mass diffusion on freely oscillating bubbles in a viscoelastic, tissue-like medium[END_REF][START_REF] Gaudron | Bubble dynamics in a viscoelastic medium with nonlinear elasticity[END_REF], 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, threedimensional (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 [START_REF] Turangan | Numerical studies of cavitation erosion on an elastic-plastic material caused by shock-induced bubble collapse[END_REF][START_REF] Ball | Shock-induced collapse of a cylindrical air cavity in water: a free-lagrange simulation[END_REF][START_REF] Turangan | Free-lagrange simulations of the expansion and jetting collapse of air bubbles in water[END_REF], front-tracking methods [START_REF] Glimm | Conservative front tracking with improved accuracy[END_REF][START_REF] Terashima | A front-tracking/ghost-fluid method for fluid interfaces in compressible flows[END_REF], and level-set/ghost fluid methods [START_REF] Pan | Phenomenology of bubble-collapse-driven penetration of biomaterial-surrogate liquid-liquid interfaces[END_REF][START_REF] Pan | A conservative interface-interaction method for compressible multi-material flows[END_REF][START_REF] Hu | A conservative interface method for compressible flows[END_REF][START_REF] Han | Adaptive multi-resolution method for compressible multiphase flows with sharp interface model and pyramid data structure[END_REF][START_REF] Chang | Direct numerical simulation of interfacial instabilities: a consistent, conservative, all-speed, sharp-interface method[END_REF]. 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 [START_REF] Shahbazi | High-order finite difference scheme for compressible multi-component flow computations[END_REF], 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 [START_REF] Mirjalili | Comparison between the diffuse interface and volume of fluid methods for simulating two-phase flows[END_REF] 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. [START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF] and Massoni et al. [START_REF] Massoni | 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[END_REF] has been widely used and can faithfully represent shock-induced collapses [START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flows[END_REF][START_REF] Beig | Maintaining interface equilibrium conditions in compressible multiphase flows using interface capturing[END_REF][START_REF] Veilleux | Transient cavitation in pre-filled syringes during autoinjector actuation[END_REF][START_REF] Xiang | Numerical investigation on the interaction of planar shock wave with an initial ellipsoidal bubble in liquid medium[END_REF] and droplet atomization [START_REF] Meng | Numerical simulations of the early stages of high-speed droplet breakup[END_REF][START_REF] Meng | Numerical simulation of the aerobreakup of a water droplet[END_REF]. Unfortunately, this model cannot predict the collapse time and minimum radius of the Rayleigh collapse problem [START_REF] Tiwari | A diffuse interface model with immiscibility preservation[END_REF][START_REF] Rasthofer | Large scale simulation of cloud cavitation collapse[END_REF]. This problem can be averted via the thermodynamically consistent model of Kapila et al. [START_REF] Kapila | Two-phase modeling of DDT in granular materials: Reduced equations[END_REF] [START_REF] Tiwari | A diffuse interface model with immiscibility preservation[END_REF][START_REF] Rasthofer | Large scale simulation of cloud cavitation collapse[END_REF], which includes a term (K∇ • u) in the volume-fraction evolution equation that follows from an asymptotic expansion of the total-disequilibrium model of Baer and Nunziato [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF]. This additional term represents compressibility in mixture regions, though unfortunately it leads to numerical instabilities during strong compression and expansion near the interface [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF][START_REF] Beig | Temperatures produced by inertially collapsing bubbles near rigid surfaces[END_REF]. 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 [START_REF] Maso | Definition and weak stability of nonconservative products[END_REF]. Instead, we propose using a pressure-disequilibrium model [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF], which relaxes the phase-specific pressures algorithmically at each time step and averts the stability issues of the K∇ • u term. This model theoretically converges to the mechanical-equilibrium model of Kapila et al. [START_REF] Kapila | Two-phase modeling of DDT in granular materials: Reduced equations[END_REF] under mesh refinement, and while it has been utilized for cavitating flows [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF][START_REF] Petitpas | Diffuse interface models for high speed cavitating underwater systems[END_REF], detonating flows [START_REF] Petitpas | Modelling detonation waves in condensed energetic materials: Multiphase CJ conditions and multidimensional computations[END_REF], surfacetension driven flows [START_REF] Schmidmayer | A model and numerical method for compressible flows with capillary effects[END_REF], droplet atomization [START_REF] Schmidmayer | Simulation de l'atomisation d'une goutte par un écoulement à grande vitesse[END_REF][START_REF] Schmidmayer | Adaptive Mesh Refinement algorithm based on dual trees for cells and faces for multiphase compressible flows[END_REF], and fracture and fragmentation in ductile materials [START_REF] Favrie | Solid-fluid diffuse interface model in cases of extreme deformations[END_REF][START_REF] Ndanou | Multi-solid and multi-fluid diffuse interface model: Applications to dynamic fracture and fragmentation[END_REF], 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 [START_REF] Pelanti | A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves[END_REF], 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 [START_REF] Godunov | A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics[END_REF][START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF]. High-order spatial reconstructions, such as MUSCL [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF][START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method[END_REF][START_REF] Schmidmayer | ECOGEN: An open-source tool for multiphase, compressible, multiphysics flows[END_REF] and WENO [START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flows[END_REF][START_REF] Tiwari | A diffuse interface model with immiscibility preservation[END_REF][START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF][START_REF] Shukla | An interface capturing method for the simulation of multi-phase compressible flows[END_REF], are often used, along with their variants WENO-Z [START_REF] Zhang | A weighted essentially non-oscillatory numerical scheme for a multi-class traffic flow model on an inhomogeneous highway[END_REF], WENO-CU6 [START_REF] Hu | An adaptive central-upwind weighted essentially nonoscillatory scheme[END_REF][START_REF] Hu | Scale separation for implicit large eddy simulation[END_REF], and TENO [START_REF] Fu | A family of high-order targeted ENO schemes for compressible-fluid simulations[END_REF]. Herein, we will consider MUSCL and the WENO of Jiang and Shu [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF], coupled with the HLLC approximate Riemann solver [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF][START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flows[END_REF][START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF] 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 [START_REF] Godunov | A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics[END_REF][START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF][START_REF] Van Leer | Towards the ultimate conservative difference scheme III. Upstreamcentered finite-difference schemes for ideal compressible flow[END_REF][START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF].

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∇ • 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.

Multicomponent flow models

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

∂q ∂t + ∇ • F (q) + h (q) ∇ • u = r (q) , (1)
where q is the state vector, F is the flux tensor, u is the velocity field, and h and r are nonconservative 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.

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. [START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF] and Massoni et al.

[2], which we call the 5-equation model. For a two-component flow, it is

q =            α 1 α 1 ρ 1 α 2 ρ 2 ρu ρE            , F =            α 1 u α 1 ρ 1 u α 2 ρ 2 u ρu ⊗ u + pI (ρE + p) u            , h =            -α 1 0 0 0 0            , r =            0 0 0 0 0            , (2)
where ρ, u, and p are the mixture density, velocity, and pressure, respectively, and α k is the volume fraction, for which k indicates the phase index. The mixture total energy is

E = e + 1 2 u 2 , ( 3
)
where e is the mixture specific internal energy

e = 2 k=1 Y k e k (ρ k , p) . (4)
In ( 4), e k is defined via an equation of state (EOS) and Y k are the mass fractions

Y k = α k ρ k ρ . (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 = (γ g -1)ρ g e g , (6)
where γ g = 1.4, and the liquid is modeled by the stiffened-gas EOS

p l = (γ l -1)ρ l e l -γ l π ∞ , (7)
where γ l and π ∞ are model parameters [START_REF] Métayer | Elaborating equations of state of a liquid and its vapor for two-phase flow models[END_REF]. We note that the stiffened-gas EOS is not a complete EOS but is reasonably accurate at the conditions considered here [START_REF] Métayer | Elaborating equations of state of a liquid and its vapor for two-phase flow models[END_REF]. More general EOS could easily be substituted in our framework. The mixture quantities are

ρ = 2 k=1 α k ρ k and p = 2 k=1 α k p k , (8)
and

ρc 2 = 2 k=1 α k ρ k c 2 k β (γ k -1) , β = 2 k=1 α k γ k -1 , ( 9
)
where c is the mixture speed of sound, and c k and γ 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 [START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF][START_REF] Massoni | 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[END_REF][START_REF] Schmidmayer | A model and numerical method for compressible flows with capillary effects[END_REF].

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

The thermodynamically consistent mechanical-equilibrium model of Kapila et al. [START_REF] Kapila | Two-phase modeling of DDT in granular materials: Reduced equations[END_REF],

which we call the 5-equation model with K∇ • u, has

q =            α 1 α 1 ρ 1 α 2 ρ 2 ρu ρE            , F =            α 1 u α 1 ρ 1 u α 2 ρ 2 u ρu ⊗ u + pI (ρE + p) u            , h =            -α 1 -K 0 0 0 0            , r =            0 0 0 0 0            , (10)
where only h is different from [START_REF] Massoni | 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[END_REF]. Here, K is

K = ρ 2 c 2 2 -ρ 1 c 2 1 ρ 2 c 2 2 α 2 + ρ 1 c 2 1 α 1 , (11)
and K∇ • u represents expansion and compression of each phase in mixture regions. In this case, the mixture speed of sound follows from

1 ρc 2 = 2 k=1 α k ρ k c 2 k , (12)
which is also the Wood speed of sound [START_REF] Wood | A textbook of sound[END_REF][START_REF] Wallis | One-dimensional two-phase flow[END_REF].

Pressure-disequilibrium model of Saurel et al. [4]

The pressure-disequilibrium model of Saurel et al. [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF], which we call the 6-equation model, is expressed as

q =               α 1 α 1 ρ 1 α 2 ρ 2 ρu α 1 ρ 1 e 1 α 2 ρ 2 e 2               , F =               α 1 u α 1 ρ 1 u α 2 ρ 2 u ρu ⊗ u + pI α 1 ρ 1 e 1 u α 2 ρ 2 e 2 u               , h =               -α 1 0 0 0 α 1 p 1 α 2 p 2               , r =               µδp 0 0 0 -µp I δp µp I δp               , (13)
where r represents the relaxation of pressures between the phases with coefficient µ. The interfacial pressure is

p I = z 2 p 1 + z 1 p 2 z 1 + z 2 , (14)
where z k = ρ k c k is the acoustic impedance of the phase k, and

δp = p 1 -p 2 , (15)
is the pressure difference between the two phases. Since p 1 = 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

∂ρE ∂t + ∇ • [(ρE + p) u] = 0. ( 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. [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF]).

The mixture speed of sound is defined according to

c 2 = 2 k=1 Y k c 2 k . (17)
After applying the infinite pressure-relaxation procedure detailed in section 3.3, the effective mixture speed of sound matches [START_REF] Gaudron | Bubble dynamics in a viscoelastic medium with nonlinear elasticity[END_REF]. We will discuss the influence of sound speed for interface problems in section 5.3.2.

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 q on a computational cell i with volume V i and surface

A with normal unit vector n is given by the explicit finite-volume Godunov [START_REF] Godunov | A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics[END_REF] scheme

q n+1 i = q n i - ∆t V i N s=1 A s F s • n s + h (q n i ) N s=1 A s u s • n s , ( 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 [START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flows[END_REF], 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 [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF][START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flows[END_REF][START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF], giving the flux tensor and flow-velocity vector F s and u s , respectively. The solution of ( 18) is restricted by the usual CFL criterion.

MUSCL scheme

We use the second-order-accurate MUSCL scheme of Schmidmayer et al. [START_REF] Schmidmayer | ECOGEN: An open-source tool for multiphase, compressible, multiphysics flows[END_REF] (labeled here as MUSCL2) with two-step time integration

q n+ 1 2 i = q n i + 1 2 ∆tL (q n i ) , (19)
q n+1 i = q n i + ∆tL q n+ 1 2 i , (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 [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF] is used on the primitive variables. The monotonized central (MC) [START_REF] Van Leer | Towards the ultimate conservative difference scheme III. Upstreamcentered finite-difference schemes for ideal compressible flow[END_REF] 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].

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 [START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF] q

(1)

i = q n i + ∆tL (q n i ) , (21) q
(2)

i = 3 4 q n i + 1 4 q (1) i + 1 4 ∆tL q (1) i , (22)
q n+1 i = 1 3 q n i + 2 3 q (2) i + 2 3 ∆tL q (2) i . (23)
This method has previously been implemented for the pressure-equilibrium models of Allaire et al. [START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF] 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. [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF].

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. [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF]. At each time step it solves the non-relaxed, hyperbolic equations (µ → 0) using ( 18), then relaxes the disequilibrium pressures for µ → +∞. Specifically, after manipulation of the equations realized in agreement with thermodynamic considerations and saturation constraint ( k α k = 1) , ( 13) is solved with respect to

k (αρ) k v k (p) = 1, (24)
where (αρ) 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) = p 0 k + γ k π ∞k + p(γ k -1) γ k (p + π ∞k ) v 0 k , (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 = ρe - k α k γ k π ∞k γ k -1 k α k γ k -1 , ( 26
)
where ρe is determined using [START_REF] Terashima | A front-tracking/ghost-fluid method for fluid interfaces in compressible flows[END_REF]. 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. [START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF], Massoni et al. [START_REF] Massoni | 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[END_REF] and Kapila et al. [START_REF] Kapila | Two-phase modeling of DDT in granular materials: Reduced equations[END_REF] for the spherical-bubble-collapse cases we consider subsequently.

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 ∞ according to the Rayleigh-Plesset equation [START_REF] Tiwari | A diffuse interface model with immiscibility preservation[END_REF][START_REF] Brennen | Cavitation and bubble dynamics[END_REF]: In the following, this pressure initialization is labeled as initial interface equilibrium with Ṙ0 = 0. We consider cases with both modest and high initial pressure ratios, as shown in Table 1. The water is parameterized by γ l = 2.35 and π ∞ = 10 9 Pa [START_REF] Beig | Maintaining interface equilibrium conditions in compressible multiphase flows using interface capturing[END_REF][START_REF] Petitpas | Diffuse interface models for high speed cavitating underwater systems[END_REF][START_REF] Pelanti | A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves[END_REF][START_REF] Métayer | Elaborating equations of state of a liquid and its vapor for two-phase flow models[END_REF][START_REF] Métayer | Modeling evaporation fronts with reactive Riemann solvers[END_REF][START_REF] Saurel | Modelling phase transition in metastable liquids: application to cavitating and flashing flows[END_REF].

p(R) = p ∞ + R 0 R (p b -p ∞ ) . (27)
Case p ∞ [Pa] p b [Pa] p ∞ /p b 1: Low-pressure-ratio 10 5 10 4 10
2: High-pressure-ratio 5 × 10 6 3550 1427

Table 1: Nominal initial conditions for the cases simulated.

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 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 [START_REF] Tiwari | A diffuse interface model with immiscibility preservation[END_REF]

α g = 1 2 1 -tanh R -R 0 2D , ( 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 = 3V b 4π 1 3
, where

V b = N i=1 α g,i V c,i (29)
is the total volume of the gas phase, N is the total number of grid cells, and α 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.915R 0 ρ l p ∞ ( 30
)
is the nominal total collapse time from its initial (maximum) radius R 0 [START_REF] Brennen | Cavitation and bubble dynamics[END_REF]. In our implementation, we compute about 69 × 10 3 and 18 × 10 3 time steps per t c for cases 1 and 2, respectively.

Results

Effect of K∇ • u on the 5-equation model

We first reconsider the behavior and influence of the K∇ • u term from the 5-equation models on the spherical-bubble-collapse problem using the WENO5 scheme as previously presented by Tiwari et al. [START_REF] Tiwari | A diffuse interface model with immiscibility preservation[END_REF].

Figure 2 shows that in both pressure-ratio cases, only the model with K∇ • u agrees with a semi-analytical solution following the Keller-Miksis equation [START_REF] Keller | Bubble oscillations of large amplitude[END_REF]; a compressible form of the Rayleigh-Plesset equation. The Keller-Miksis equation is based on an asymptotic (a)

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t >

(b)

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > 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∇ • u to represent spherical bubble collapse was previously observed by Tiwari et al. [START_REF] Tiwari | A diffuse interface model with immiscibility preservation[END_REF], who attributed the better results of the model of Kapila et al. [START_REF] Kapila | Two-phase modeling of DDT in granular materials: Reduced equations[END_REF] 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∇ • u along a radial coordinate at two instances in time. For t ≈ 0, the initial interface smearing results in a mixture region at the interface,

for which K = 0, but K∇ • u ≈ 0 because ∇ • u ≈ 0.
However, during the collapse, ∇ • u = 0 and the fluid volume fractions are modified. We see that K∇ • 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∇ • u term. This keeps the interface relatively sharp and results in the larger interface velocity observed in t = 0.85t c < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > t ⇡ 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > (a)

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t >

(b)

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > 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∇ • u, since this term accounts for expansion and compression in mixture regions.

Comparison of the 5-and 6-equation models

While the 5-equation model with K∇ • 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∇ • u term [START_REF] Maso | Definition and weak stability of nonconservative products[END_REF] as discussed in section 1. Here, we consider the 6-equation model as a potential (a)

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t >

(b)

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > (a)

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t >

(b)

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t >

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.

Spherical bubble collapse with initial interface equilibrium

We first consider the case of initial interface equilibrium, Ṙ0 = 0. Figure 5 (a) shows the interface evolution for p ∞ /p b = 10, spatial resolution N R 0 = 25 and for the FV1, MUSCL2, WENO3 and WENO5 schemes. In addition to the MC [START_REF] Gammie | HARM: A numerical scheme for general relativistic magnetohydrodynamics[END_REF] slope limiter, the Minmod [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF][START_REF] Sweby | High resolution schemes using flux limiters for hyperbolic conservation laws[END_REF] 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 ε as

ε = 1 N t Nt i=0 R(t i ) -R KM (t i ) R KM (t i ) , ( 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 [START_REF] Godunov | A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics[END_REF][START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF][START_REF] Van Leer | Towards the ultimate conservative difference scheme. II. Monotonicity and conservation combined in a second-order scheme[END_REF]. The WENO5 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 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 [START_REF] Brennen | Cavitation and bubble dynamics[END_REF][START_REF] Frost | Effects of ambient pressure on the instability of a liquid boiling explosively at the superheat limit[END_REF], 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 [START_REF] Wadell | Volume, shape, and roundness of quartz particles[END_REF] Ψ = π

1 3 (6V b ) 2 3 A b , (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 WENO5 scheme maintains sphericity during the entire collapserebound process, similar to that observed by Tiwari et al. [START_REF] Tiwari | A diffuse interface model with immiscibility preservation[END_REF], Whereas the MUSCL2 and

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t >

(b)

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > 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 ∞ /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 [START_REF] Fuster | Liquid compressibility effects during the collapse of a single cavitating bubble[END_REF]. 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 ε over the temporal window t ∈ [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 ε 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 of the bubble surface are shown Figure 8.

Almost no grid-based artifacts on the bubble surface are visible until t ≈ t c for all numerical methods, at which point Ψ decreases significantly. Compared to the low-pressure-ratio case, the interface evolves more quickly and all methods conserve sphericity for t t c . However, after the collapse, significant distortions are visible and Ψ 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 Ψ 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.

c water < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > c air < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t >

Spherical bubble collapse with initial interface disequilibrium

Lastly, we consider the case of initial interface disequilibrium, and thus Ṙ0 (t = 0) = 0.

We enforce this by setting the internal and external interface pressures to different values as

p =      p b for 0 ≤ R ≤ R 0 , p ∞ otherwise. ( 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 (α = 0 or 1). As shown in Figure 10, the Wood speed of sound [START_REF] Gaudron | Bubble dynamics in a viscoelastic medium with nonlinear elasticity[END_REF], which is also the speed of sound of the 5-equation model with K∇ • u, varies in this region and is much less than that of either of the pure phases (α = 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 α and pressure p both modify the pressure-wave behaviors uniquely, though both ultimately pollute the bubble dynamics.

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 pressuredisequilibrium events, such as a collapsing ellipsoidal bubble near a wall [START_REF] Pishchalnikov | Experimental observations and numerical modeling of lipid-shell microbubbles with calcium-adhering moieties for minimally-invasive treatment of urinary stones[END_REF][START_REF] Pishchalnikov | High-speed video microscopy and numerical modeling of bubble dynamics near a surface of urinary stone[END_REF], would

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t >

(b)

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > 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 ε for the lowpressure-ratio case and 59% smaller ε 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 lowpressure-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 [START_REF] Ii | An interface capturing method with a continuous function: The THINC method with multi-dimensional reconstruction[END_REF]. 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 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 pressuredisequilibrium events.

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. [START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF] and Massoni et al. [START_REF] Massoni | 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[END_REF] is unable to accurately represent a spherical bubble collapse and demonstrated how the additional K∇ • u term naturally present in the model of Kapila et al. [START_REF] Kapila | Two-phase modeling of DDT in granular materials: Reduced equations[END_REF] is required to ensure good agreement with the Keller-Miksis solution [START_REF] Tiwari | Growth-and-collapse dynamics of small bubble clusters near a wall[END_REF]. Since the 5-equation model with K∇ • u is known to produce instabilities in some numerical experiments [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF][START_REF] Beig | Temperatures produced by inertially collapsing bubbles near rigid surfaces[END_REF], we investigated the 6-equation pressuredisequilibrium 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∇ • 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 surrogatetruth 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 [START_REF] Shyue | An Eulerian interface sharpening algorithm for compressible twophase flow: The algebraic THINC approach[END_REF] 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.

computations utilized the Extreme Science and Engineering Discovery Environment, which is supported by the National Science Foundation grant number CTS120005.

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.

Air, p b Water, p ∞ R 0 Figure 1 :

 01 Figure 1: Problem configuration for a collapsing spherical bubble.

Figure 2 :

 2 Figure 2: Radial bubble-wall evolution for (a) p ∞ /p b = 10 with N R0 = 25 and (b) p ∞ /p b = 1427 with N R0 = 50. Solutions are computed using the 5-equation models with WENO5 as well as the Keller-Miksis equation.

Figure 3 :

 3 Figure 3: Liquid volume-fraction, pressure, and K∇ • u for varying radial position R and the case p ∞ /p b = 10 and N R0 = 25 using the 5-equation model with K∇ • u. Times (a) t ≈ 0 and (b) t = 0.85t c are shown.

Figure 2 .

 2 Figure 2.

Figure 4 :

 4 Figure 4: Radial bubble-wall evolution for (a) p ∞ /p b = 10 and (b) p ∞ /p b = 1427. Solutions are computed using WENO5.

Figure 5 :

 5 Figure 5: Radial bubble-wall evolution for initial interface equilibrium and p ∞ /p b = 10. (a) Results for all schemes and flux limiters as labeled (fixed resolution N R0 = 25) and (b) spatial convergence of the numerical methods.

Figure 6 :

 6 Figure 6: Evolution of the bubble sphericity for p ∞ /p b = 10 and N R0 = 25. Nominal bubble shapes as represented by α = 0.5 isosurfaces are also shown for times t = 0.7t c , t (R = R min ), and 2t c .

Figure 7 :

 7 Figure 7: Radial bubble-wall evolution for initial interface equilibrium and p ∞ /p b = 1427. (a) Results (fixed resolution N R0 = 50) and (b) spatial convergence of the MUSCL2, WENO3 and WENO5 schemes.

Figure 8 :

 8 Figure 8: Bubble sphericity evolution for p ∞ /p b = 1427 and N R0 = 50. Nominal bubble shapes (α = 0.5) are also shown for times t = 0.5t c , 1.02t c , and 1.35t c .

Figure 9 :

 9 Figure9: Bubble sphericity and associated bubble interface for varying mesh resolution N R0 at time t = 1.35t c for the p ∞ /p b = 1427 case and MUSCL2 method. For N R0 = 100 and 150 the adaptive-mesh-refinement technique of Schmidmayer et al.[START_REF] Schmidmayer | Adaptive Mesh Refinement algorithm based on dual trees for cells and faces for multiphase compressible flows[END_REF] is used to minimize computational expense.

Figure 10 :

 10 Figure 10: Wood speed of sound for water-air mixture. Here, c water = 1625 m/s and c air = 350 m/s.

Figure 11 :

 11 Figure 11: Radial bubble-wall evolution for initial interface disequilibrium, p ∞ /p b = 1427 and N R0 = 50. Solutions are computed using the 6-equation model, and the Keller-Miksis result is shown as surrogate truth.

Figure 13 :Table 2 :

 132 Figure 13: Radial bubble-wall evolution for initial interface equilibrium and numerical methods as labeled; (a) p ∞ /p b = 10 and N R0 = 25, (b) p ∞ /p b = 1427 and N R0 = 50.

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

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > (a)

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > (b) < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > (c) < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > 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 interfacesharpening method [START_REF] Shyue | An Eulerian interface sharpening algorithm for compressible twophase flow: The algebraic THINC approach[END_REF] 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 [START_REF] So | Anti-diffusion interface sharpening technique for two-phase compressible flow simulations[END_REF] or regularization methods [START_REF] Tiwari | A diffuse interface model with immiscibility preservation[END_REF], 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 [START_REF] Shyue | An Eulerian interface sharpening algorithm for compressible twophase flow: The algebraic THINC approach[END_REF] and Deng et al. [START_REF] Deng | High fidelity discontinuity-resolving reconstruction for compressible multiphase flows with moving interfaces[END_REF]), and so we can proceed with a faithful comparison for collapsing spherical bubbles.

MUSCL2

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > MUSCL2 + THINC

< l a t e x i t s h a 1 _ b a s e 6 4 = " B q h u K Z 6 y A z 5 j Z S p x T 9 q O 3 8 d t

r B f r 3 f q Y t 5 a s Y q Y K f 2 B 9 / g B T n J M + < / l a t e x i t > t = 0.7t c < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > t(R = R min ) < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > t = 2t c < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > WENO3 < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > WENO3 + THINC

z s e 8 t e D k M 4 f w B 8 7 n D y f R j R s = < / l a t e x i t > t = 0.5t c < l a t e x i t s h a 1 _ b a s e 6 4 = " n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > t = 1.35t c < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > WENO3 < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > t = 1.02t c < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > WENO3 + THINC

S P m t p 6 h O P K S e e P j L G x 1 r p 4 l 4 g d f m

MUSCL2

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > b)

< l a t e x i t s h a 1 _ b a s e 6 4 = " K C H Z q B A + 5 d k P n r 4 F + 0 f J e j a 0 Z n We simulate the flow in the shock tube using both the 5-equation with K∇ • u and 6equation 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 [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF][START_REF] Pelanti | A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves[END_REF]. The problem setup is shown in Figure 17; there is a uniform initial pressure p = 10 5 Pa and densities ρ l = 10 3 kg/m 3 and ρ 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 γ l = 4.4 and π ∞ = 6 × 10 8 Pa.

Figure 18 shows the results of the primitive variables at t = 1.85 ms for the vacuum (a)

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t >

(b)

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > (c)

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t >

(d)

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > (b)

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > (c)

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t >

(d)

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t >