Spencer H Bryngelson

Kevin Schmidmayer

Vedran Coralic

Jomela C Meng

Kazuki Maeda

Tim Colonius
email: colonius@caltech.edu

MFC: An open-source high-order multi-component, multi-phase, and multi-scale compressible flow solver

Keywords: computational fluid dynamics, multi-phase flow, diffuse-interface method, compressible flow, ensemble averaging, bubble dynamics Program

published or not. The documents may come

Introduction

The multi-component flow code (MFC) is an open source high-order solver for multi-phase and multi-component flows. Such flows are central to a wide range of engineering problems. For example, cavitating flow phenomena are of critical importance to the development of artificial heart valves and pumps [START_REF] Brennen | Cavitation in medicine[END_REF], minimizing injury due to blast trauma [START_REF] Laksari | Computational simulation of the mechanical response of brain tissue under blast loading[END_REF][START_REF] Proud | The high-strain rate loading of structural biological materials[END_REF], and improving shockand burst-wave lithotripsy treatments [START_REF] Coleman | Acoustic cavitation generated by an extracorporeal shockwave lithotripter[END_REF][START_REF] Pishchalnikov | Cavitation bubble cluster activity in the breakage of kidney stones by lithotripter shockwaves[END_REF][START_REF] Ikeda | Cloud cavitation control for lithotripsy using high intensity focused ultrasound[END_REF]. Bubble cavitation is also pervasive in flows around hydrofoils, submarines, and high-velocity projectiles [START_REF] Saurel | Modelling phase transition in metastable liquids: Application to cavitating and flashing flows[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], during underwater explosions [START_REF] Etter | Underwater acoustic modeling and simulation[END_REF][START_REF] Kedrinskii | Negative pressure profile in cavitation zone at underwater explosion near free surface[END_REF], and within pipe systems and hydraulic machinery [START_REF] Streeter | Transient cavitating pipe flow[END_REF][START_REF] Weyler | An investigation of the effect of cavitation bubbles on the momentum loss in transient pipe flow[END_REF]. Unfortunately, cavitation in these settings is usually detrimental, causing noise and material deterioration [START_REF] Streeter | Transient cavitating pipe flow[END_REF][START_REF] Weyler | An investigation of the effect of cavitation bubbles on the momentum loss in transient pipe flow[END_REF]. Other cases of interest include the breakup of liquid droplets and jets [START_REF] Meng | Numerical simulation of the aerobreakup of a water droplet[END_REF][START_REF] Meng | Numerical simulations of the early stages of high-speed droplet breakup[END_REF][START_REF] Schmidmayer | A model and numerical method for compressible flows with capillary effects[END_REF], erosion of aircraft surfaces during supersonic flight [START_REF] Engel | Fragmentation of waterdrops in the zone behind an air shock[END_REF][START_REF] Joseph | Breakup of a liquid drop suddenly exposed to a high-speed airstream[END_REF], shock-wave attenuation of nuclear blasts [START_REF] Chauvin | Shock waves in sprays: numerical study of secondary atomization and experimental comparison[END_REF], and needle-free injection for drug delivery [START_REF] Tagawa | Needle-free injection into skin and soft matter with highly focused microjets[END_REF][START_REF] Veilleux | Pressure and stress transients in autoinjector devices[END_REF].

Robust simulation of multi-component compressible flow phenomena requires a numerical method that maintains discrete conservation, suppresses oscillations near discontinuities, and preserves numerical stability. For such simulations to be computationally efficient, the method used should also be high-order accurate away from discontinuities. Schemes that can potentially achieve these requirements can be classified as either interface-tracking or interface-capturing [START_REF] Fuster | A review of models for bubble clusters in cavitating flows[END_REF]. Examples of interface-tracking methods include free-Lagrange [START_REF] Ball | Shock-induced collapse of a cylindrical air cavity in water: a free-Lagrange simulation[END_REF][START_REF] Turangan | Numerical studies of cavitation erosion on an elastic-plastic material caused by shock-induced bubble collapse[END_REF], front-tracking [START_REF] Glimm | Conservative front tracking with improved accuracy[END_REF][START_REF] Cocchi | A Riemann problem based method for the resolution of compressible multimaterial flows[END_REF], and level-set/ghost fluid schemes [START_REF] Abgrall | Computations of compressible multifluids[END_REF][START_REF] Liu | Ghost fluid method for strong shock impacting on material interface[END_REF][START_REF] Liu | A conservative modification to the ghost fluid method for compressible multiphase flows[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 multi-phase 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]; these methods differ from interface-capturing as they treat material interfaces as sharp features in the flow. This allows interfacial fluids to have differing equations of state and ensures that interfacial physics are straightforward to implement. Unfortunately, they do not naturally enforce discrete or total conservation, though there have been some recent attempts to partially include these properties [START_REF] Denner | Pressure-based algorithm for compressible interfacial flows with acoustically-conservative interface discretisation[END_REF][START_REF] Fuster | An all-Mach method for the simulation of bubble dynamics problems in the presence of surface tension[END_REF]. Interface-capturing methods instead treat interfaces as discontinuities in material properties via advected volume fractions. Such methods are generally more efficient than interface-tracking schemes [START_REF] Mirjalili | Comparison between the diffuse interface and volume of fluid methods for simulating two-phase flows[END_REF] and can achieve discrete conservation by simply solving the governing equations in conservative form. However, they also smear material interfaces via numerical diffusion [START_REF] Johnsen | Preventing numerical errors generated by interface-capturing schemes in compressible multi-material flows[END_REF][START_REF] Perigaud | A compressible flow model with capillary effects[END_REF][START_REF] Shyue | A fluid-mixture type algorithm for compressible multicomponent flow with van der waals equation of state[END_REF]. In these smeared regions, we must ensure that the mixture properties are treated in a thermodynamically-and numerically-consistent way, such that spurious oscillations (and other physical inconsistencies) are avoided in the presence of high density contrasts between materials. Fortunately, such methods have been developed and prove to be a robust treatment of interfacial dynamics [START_REF] Johnsen | Preventing numerical errors generated by interface-capturing schemes in compressible multi-material flows[END_REF][START_REF] Johnsen | Implementation of weno schemes in compressible multicomponent flow problems[END_REF][START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flow problems[END_REF]. Still, the lack of a coherent material interface means that interfacial physics are more challenging to implement than interface-tracking methods, though conservative treatments do exist for both interfacial heat transfer [START_REF] Massoni | Some models and Eulerian methods for interface problems between compressible fluids with heat transfer[END_REF] and capillary effects [START_REF] Schmidmayer | A model and numerical method for compressible flows with capillary effects[END_REF][START_REF] Perigaud | A compressible flow model with capillary effects[END_REF][START_REF] Meng | Numerical simulations of droplet aerobreakup[END_REF].

Here, we choose to use an interface-capturing scheme because computational efficiency and discrete conservation are of principal importance to many problems of interest. The interfacecapturing schemes we use follow from the so-called 5-and 6-equation models, which are known to be sufficient for representing a wide range of flow phenomenologies [START_REF] Massoni | Some models and Eulerian methods for interface problems between compressible fluids with heat transfer[END_REF][START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF][START_REF] Kapila | Two-phase modeling of DDT in granular materials: Reduced equations[END_REF][START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF]. They are complemented by a discretely-conservative numerical method that solves the conservative form of the compressible flow equations. To maintain a non-oscillatory behavior near material interfaces, the material advection equations are formulated in a quasi-conservative form [START_REF] Abgrall | How to prevent pressure oscillations in multicomponent flow calculations: a quasi conservative approach[END_REF]. The equations of motion are then closed by a thermodynamically consistent set of mixture rules [START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF]. The governing equations are solved with a shock-capturing finite-volume method. Specifically, we adopt a WENO spatial reconstruction that can achieve high-order accuracy while maintaining non-oscillatory behavior near material interfaces [START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flow problems[END_REF][START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF][START_REF] Balsara | Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy[END_REF][START_REF] Henrick | Mapped weighted essentially non-oscillatory schemes: achieving optimal order near critical points[END_REF]. This scheme is then coupled with an HLL-type approximate Riemann solver [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF][START_REF] Toro | Restoration of the contact surface in the HLL-Riemann solver[END_REF] and a total-variation-diminishing (TVD) time stepper [START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF].

MFC is, of course, not the only viable option for simulation of compressible multi-phase flows. For example, ECOGEN [START_REF] Schmidmayer | an open-source tool dedicated to multiphase compressible multiphysics flows[END_REF] offers a pressure-disequilibrium-based interface-capturing scheme that is well suited for the same problems reached by MFC. However, ECOGEN is built upon an intrinsically low-order MUSCL scheme that can inhibit both efficient simulation and physically fidelity when compared to the WENO schemes we use here, especially when augmented with our high-order cell-average approximations. This is demonstrated in section 5.4 for an isentropic vortex problem and in section 5.3 and Schmidmayer et al. [START_REF] Schmidmayer | An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics[END_REF] for cavitating gas bubbles. In pursuit of this, MFC was also constructed with a phase-averaged flow model that represent unresolved multi-phase dynamics at the sub-grid level. More mature CFD solvers such as OpenFOAM are also available. OpenFOAM natively supports finite-volume methods for multi-component flows [START_REF] Weller | A tensorial approach to computational continuum mechanics using object-oriented techniques[END_REF], though higherorder methods and interface-capturing models are only available via links to external forked projects. MFC instead offers an integrated approach that avoids any conflicts from such libraries. Finally, the parallel I/O file systems we employ ensures that the MFC architecture can scale up to the largest modern HPC systems.

Herein, we describe version 1.0 of MFC. In section 2 we present an overview of what is included in the MFC package, including its organization, features, and logistics. The physical models we use are presented in section 3, including the associated mixture rules, governing equations, equations of state, and our implementation of an ensemble-averaged bubbly flow model. The numerical methods used to solve the associated equations are presented in section 4. A series of test cases simulated using MFC are discussed in section 5; these verify and validate our method. These are complemented by a set of illustrative example cases that further demonstrate MFC's capabilities in section 6 and parallel benchmarking in section 7, which analyzes the performance of MFC on large scale computing clusters. Section 8 concludes our presentation of MFC.

Overview and features

Package, installation, and testing

MFC is available at https://mfc-caltech.github.io. The source code is written using Fortran 90 with MPI bindings used for parallel communication and Python scripts are used to generate input files. Installation of MFC requires the FFTW package for cylindrical coordinate treatment [START_REF] Frigo | FFTW: An adaptive software architecture for the FFT[END_REF], and, optionally, Silo and its dependencies for post-treatment of data files [START_REF] Miller | Silo-a mesh and field I/O library and scientific database[END_REF]. We only provide the FFTW package with MFC to keep the package relatively small.

The MFC package includes several directories; their organization and descriptions are shown in table 1. New users should consult the CONFIGURE and INSTALL files for instructions on how to compile MFC. In brief, the user must ensure that Python and an MPI Fortran compiler are loaded and that the FFTW package can be located by Makefile; this can be done by pointing Makefile.user to the correct location or by installing the included FFTW distribution in the installers directory. Once the software has been built, the test target of Makefile should be called; it runs multiple tests (which are located in the tests directory) to ensure that MFC is operating as intended.

Features

MFC uses a fully parallel environment via message passing interfaces (MPI), the performance of which is the subject of section 7. Computationally, it includes structured Cartesian and cylindrical grids with non-uniform mesh stretching available; characteristic-based Thompson, periodic, and free-slip boundary conditions have also been implemented. be used with a flexible number of components, as discussed in section 3. Ensemble-averaged dilute bubbly flow modeling is also available, including options for Gilmore and Keller-Miksis single-bubble models (see section 3.3). The numerical methods are discussed in section 4; they include 1-, 3-, and 5-th order accurate WENO reconstructions on optionally the primitive, conservative, or characteristic variables. Within each finite volume, high-order evaluation is available for the cell-averaged variables via Gaussian quadrature for multi-dimensional problems. The shock-capturing schemes are paired with either HLL, HLLC, or exact Riemann solvers. For time-stepping, 1-5-th order accurate Runge-Kutta methods are available. These features will, of course, evolve and expand with time.

Of great practical importance are the user interfaces we utilize. MFC features Python input scripts, which operate via dictionaries to automatically write input files that are read via Fortran namelists. Additionally, the file system and data formats were selected to enable large-scale parallel simulations. Specifically, we use the Lustre file system to generate and read restart files; it can support gigabyte-per-second-scale IO operations and petabyte-scale storage requirements [START_REF] Halbwachs | The synchronous data flow programming language lustre[END_REF], which ensures that the MFC can utilize the full capabilities of modern HPC systems. We also utilize HDF5 Silo databases, which keeps the file structure compact and enables parallel visualization. We have confirmed that MFC works as expected on various high-performance computing platforms, including modern SGI-and Dell-based supercomputers.

Software structure 2.3.1. Pre-processing

The pre-processor generates initial conditions and spatial grids from the physical patches specified in the Python input file and exports them as binary files to be read by the simulator. Specifically, this involves allocating and writing either a Cartesian or cylindrical mesh, with the option of mesh stretching, according to the input parameters. The specified physical variables for each patch are transformed into their conservative form and written in a manner consistent with the mesh. The pre-processor is comprised of individual Fortran modules that read input values and export mesh and initial condition files, assign then distribute global variables via MPI, perform variable transformations, generate grids, parse and assign patch types, and check that specified input variables are physically consistent and that specified options do not contradict each other.

Simulation

The simulator, given the initial-condition files generated by the pre-processor, solves the corresponding governing flow equations with the specified boundary conditions using our interfacecapturing numerical method. Simulations are conducted for the number of time steps indicated. The simulator exports run-time information, restart files that can be used to either restart the simulation or post-process the associated data, and, optionally, human-readable output data. The structure of the simulator follows that of pre-processor, with individual Fortran modules conducting each software component; this includes reading and exporting data and grid files, performing Fourier transforms, assigning and distributing global variables via MPI, performing variable transformations, computing time and spatial derivatives using WENO and the Riemann solver specified, computing boundary values, including ensemble-averaged bubbly flow physics, and checking that the input variables are valid.

Post-processing

The post-processor reads simulation data and exports HDF5/Silo databases that include variables and derived variables, as specified in the input file. Since the simulator can export human-readable data, post-processing is not essential for the usage of MFC, but is a useful tool, especially for large or parallel data structures. Specifically, the post-process component of the MFC reads the restart files exported by the simulator at distinct time intervals and computes the necessary derived quantities. The HDF5 database is then generated and exported, and can be readily viewed using, for example, VisIt [START_REF] Childs | VisIt: An end-user tool for visualizing and analyzing very large data[END_REF] or Paraview [START_REF] Ahrens | Paraview: An end-user tool for large data visualization[END_REF]. Again, individual Fortran modules perform the associated tasks, including reading data, parameter conversion, assigning and distributing global MPI variables, computing Fourier transforms, exporting HDF5 Silo databases, and checking that the input parameters are consistent.

Description of input/output files

We next describe the contents of a case-specific directory and its logistics. The specific file structure is shown in table 2. The Python script input.py is used to generate the input files ( * .inp) for the source codes and execute an MFC component (one of pre process, simulation, or post process) either in the active window or as a submitted batch script. This Python file contains the input parameters available for the MFC. MFC, depending upon the component used and options selected, will generate several files and directories. If enabled, the run time.inf file will be generated by the simulator and includes details about the current time step, simulation time, and stability criterion. Directories that contain binary restart data and output files for visualization or further post-treatment are generated, again depending upon the specified options, by the simulator and post-processor.

Physical model and governing equations

The mechanical-equilibrium compressible multi-component flow models we use can be written as

∂q ∂t + ∇ • F (q) + h (q) ∇ • u = s (q) , (1)
where q is the state vector, F is the flux tensor, u is the velocity field, and h and s are nonconservative quantities we describe subsequently.

5-equation model

We first introduce our implementation of 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]. Our multi-component implementation can be used for N k components, though we present a two-component (N k = 2) configuration here for demonstration purposes. It consists of five partial differential equations as

q =       α 1 α 1 ρ 1 α 2 ρ 2 ρu ρE       , F =       α 1 u α 1 ρ 1 u α 2 ρ 2 u ρu ⊗ u + pI -T (ρE + p) u -T • u       , h =       -α 1 -K 0 0 0 0       , s = 0, (2)
where ρ, u, and p are the mixture density, velocity, and pressure, respectively, α k is the volume fraction of component k, and T is the viscous stress tensor

T = 2η D - 1 3 (∇ • u)I , (3)
where η is the mixture shear viscosity and

D = 1 2 ∇u + (∇u) (4)
is the strain rate tensor. The mixture total and internal energies are E = e + u 2 /2 and

e = N k k=1 Y k e k (ρ k , p) , (5)
respectively, where Y k = α k ρ k /ρ are the mass fractions of each component. We close (5) using the stiffened-gas equation of state, which is chosen for its ability to faithfully model both liquids and gases [START_REF] Menikoff | The Riemann problem for fluid-flow of real materials[END_REF]; for component k it is

p k = (γ k -1)ρ k e k -γ k π ∞,k , (6)
where γ is the specific heat ratio and π ∞ is the liquid stiffness (gases have π ∞ = 0) [START_REF] Métayer | Elaborating equations of state of a liquid and its vapor for two-phase flow models[END_REF]. For liquids, these are usually interpreted as fitted parameters from shockwave-Hugoniot data [START_REF] Marsh | LASL Shock Hugoniot Data, Los Alamos Series on Dynamic Material Properties[END_REF][START_REF] Gojani | Shock Hugoniot and equations of states of water, castor oil, and aqueous solutions of sodium chloride, sucrose and gelatin[END_REF]. The speed of sound of each component is then

c k = γ k (p k + π ∞,k ) ρ k . (7)
The K term in h of (2) represents expansion and compression in mixture regions. For an

N k = 2 configuration it is K = ρ 2 c 2 2 -ρ 1 c 2 1 ρ 2 c 2 2 α 2 + ρ 1 c 2 1 α 1 (8)
and the mixture speed of sound c follows as the so-called Wood speed of sound [START_REF] Wood | A textbook of sound[END_REF][START_REF] Wallis | One-dimensional two-phase flow[END_REF]

1 ρc 2 = N k k=1 α k ρ k c 2 k . (9)
Ultimately, the equations are closed by the usual set of mixture rules

1 = N k k=1 α k , ρ = N k k=1 α k ρ k , ρe = N k k=1 α k ρ k e k , and η = N k k=1 α k η k . ( 10
)
We note that the models 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 | Some models and Eulerian methods for interface problems between compressible fluids with heat transfer[END_REF] do not include the K term in (2) and thus do not strictly obey the second-law of thermodynamics, nor reproduce the correct mixture speed of sound [START_REF] Pelanti | A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves[END_REF]. While MFC also supports these models, accurately representing the sound speed is known to be important for some problems, such as the cavitation of gas bubbles [START_REF] Schmidmayer | An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics[END_REF]. However, it is also known that the K term can result in numerical instabilities for problems with strong compression or expansion in mixture regions due to its non-conservative nature [START_REF] Schmidmayer | An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics[END_REF]. Thus, the decision of what 5-equation model to use (if any) is problem dependent and left to the user.

6-equation model

While the 5-equation model described in section 3.1 is efficient and represents the correct physics, the K∇ • u term that makes the model thermodynamically consistent can sometimes introduce numerical instabilities [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF][START_REF] Schmidmayer | An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics[END_REF]. In such cases, a pressure-disequilibrium is preferable [START_REF] Schmidmayer | An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics[END_REF]. We also support the 6-equation 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 for a two-component configuration 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 -T α 1 ρ 1 e 1 u α 2 ρ 2 e 2 u         , h =         -α 1 0 0 0 α 1 p 1 α 2 p 2         , s =         µδp 0 0 0 -µp I δp -α 1 T 1 : ∇u µp I δp -α 2 T 2 : ∇u         , (11)
where T k are the component-specific viscous stress tensors and the other terms of s represent the relaxation of pressures between components with coefficient µ. The interfacial pressure is

p I = z 2 p 1 + z 1 p 2 z 1 + z 2 , ( 12
)
where

z k = ρ k c k is the acoustic impedance of component k and δp = p 1 -p 2 , (13)
is the pressure difference. Since p 1 = p 2 , the total energy equation of the mixture is replaced by the internal-energy equation for each component. The mixture speed of sound is defined according to

c 2 = 2 k=1 Y k c 2 k , (14)
though after applying the numerical infinite pressure-relaxation procedure detailed in section 4.6 the effective mixture speed of sound matches (9).

Bubbly flow model

Implementation

MFC includes support for the ensemble-phase-averaged bubbly flow model of Zhang and Prosperetti [START_REF] Zhang | Ensemble phase-averaged equations for bubbly flows[END_REF], and our implementation of it matches that of Bryngelson et al. [START_REF] Bryngelson | A quantitative comparison of phase-averaged models for bubbly, cavitating flows[END_REF]. The bubble population has void fraction α b , which is assumed to be small, and the carrier components have mixture pressure p l . The equilibrium radii of the bubble population are represented discretely as R o , which are N bin bins of an assumed log-normal PDF with standard deviation σ p [START_REF] Colonius | Statistical equilibrium of bubble oscillations in dilute bubbly flows[END_REF]. The instantaneous bubble radii are a function of these equilibrium states as

R(R o ) = {R 1 , R 2 , . . . , R N bin }.
The total mixture pressure is modified as

p = (1 -α b )p l + α b R 3 p bw R 3 -ρ R 3 Ṙ2 R 3 , ( 15
)
where Ṙ are the bubble radial velocities and p bw are the bubble wall pressures. Overbars • denote the usual moments with respect to the log-normal PDF. The bubble void fraction is advected as

∂α b ∂t + u • ∇α b = 3α b R 2 Ṙ R 3 , (16)
and the bubble dynamic variables are evolved as

∂nφ ∂t + ∇ • (nφu) = n φ, (17)
where φ ≡ R, Ṙ, p b , m v (see section 3.3.2) and n is the conserved bubble number density per unit volume

n = 3 4π α b R 3 . ( 18
)

Single-bubble dynamics

A partial differential equation following ( 17) is evolved for each bin representing equilibrium radius R o . These equations assume that each bubble evolves, without interaction with its neighbors, in an otherwise uniform flow whose properties are dictated by the local mixture-averaged flow quantities [START_REF] Ando | Effects of polydispersity in bubbly flows[END_REF]. We also assume that the bubbles remain spherical, maintain a uniform internal pressure, and do not break-up, or coalesce. Our model includes the thermal effects, viscous and acoustic damping, and phase change. The bubble radial accelerations R are computed by the Keller-Miksis equation [START_REF] Keller | Bubble oscillations of large amplitude[END_REF]:

R R 1 - Ṙ c b + 3 2 Ṙ2 1 - Ṙ 3c b = p bw -p l ρ 1 + Ṙ c b + R ṗbw ρc b , (19)
where c b is the usual speed of sound associated with the bubble and

p bw = p b - 4µ Ṙ R - 2σ R (20)
is the bubble wall pressure, for which p b is the internal bubble pressure, σ is the surface tension coefficient, and µ is the liquid viscosity. The evolution of p b is evaluated using the model of Ando [START_REF] Ando | Effects of polydispersity in bubbly flows[END_REF]:

ṗb = 3γ b R R v T bw ṁv -Ṙp b + γ b -1 γ b λ bw ∂T ∂r r=w , ( 21
)
where T is the temperature, λ is the thermal conductivity, R v is the gas constant and γ b is the specific heat ratio of the gas. Mass transfer of the bubble contents follows the reduced model of Preston et al. [START_REF] Preston | A reduced-order model of diffusion effects on the dynamics of bubbles[END_REF] as

ṁv = Dρ bw 1 -χ vw ∂χ v ∂r w . ( 22
)

Solution method

Our numerical scheme generally follows that of Coralic and Colonius [START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flow problems[END_REF]. The spatial discretization of (1) in three-dimensional Cartesian coordinates is

∂q ∂t + F x (q) ∂x + F y (q) ∂y + F z (q) ∂z = s(q) -h(q)∇ • u, (23)
where F x i are the i ∈ (x, y, z)-direction flux vectors and the treatment of ∇ • u is discussed later.

Treatment of spatial derivatives

We use a finite volume method to treat the spatial derivatives of [START_REF] Ball | Shock-induced collapse of a cylindrical air cavity in water: a free-Lagrange simulation[END_REF]. The finite volumes are

I i,j,k = [x i-1/2 , x i+1/2 ] × [y j-1/2 , y j+1/2 ] × [z k-1/2 , z k+1/2 ]. (24)
We spatially integrate [START_REF] Ball | Shock-induced collapse of a cylindrical air cavity in water: a free-Lagrange simulation[END_REF] within each cell-centered finite volume as

dq i,j,k dt = 1 ∆x i [F x i-1/2,j,k -F x i+1/2,j,k ]+ 1 ∆y j [F y i,j-1/2,k -F y i,j+1/2,k ]+ 1 ∆z k [F z i,j,k-1/2 -F z i,j,k+1/2 ] + s(q i,j,k ) -h(q i,j,k )(∇ • u) i,j,k , (25)
where

q i,j,k = 1 V i,j,k I i,j,k q(x, y, z, t) dxdydz, ( 26
)
s i,j,k = 1 V i,j,k I i,j,k s(x, y, z, t) dxdydz, ( 27
)
F i+1/2,j,k = 1 ∆y j ∆z k A i+1/2,j,k F (x, y, z, t) dydz, (28)
are cell-volume and face averages, for which ∆x i = x i+1/2 -x i-1/2 are the mesh spacings (∆y and ∆z have the same form), and V i,j,k and A i+1/2,j,k are the cell volumes and face areas. MFC can approximate this equation using high-order quadratures as opposed to simple cellcentered averages. This approach computes the flux surface integrals and source terms using a two-point, fourth-order, Gaussian quadrature rule, e.g.

F i+1/2,j,k = 1 4 2 m=1 2 l=1 F (q(x i+1/2 , y j l , z km )), (29)
where l and m are the Gaussian quadrature point indices and

y j l = y j + (2l + 1) ∆y j 2 √ 3 and z km = z k + (2m + 1) ∆z k 2 √ 3 . ( 30
)
The divergence terms are treated using a midpoint rule

(∇ • u) i,j,k = 1 ∆x i (u i+1/2,j,k -u i-1/2,j,k ) + 1 ∆y j (v i,j+1/2,k -v i,j-1/2,k ) + 1 ∆z k (w i,j,k+1/2 -w i,j,k-1/2 ), (31)
where u = {u, v, w} are the cell-averaged velocity components computed analogous to [START_REF] Liu | A conservative modification to the ghost fluid method for compressible multiphase flows[END_REF].

To avoid spurious oscillations at material interfaces, we ultimately evaluate the fluxes by reconstructing the primitive or characteristic variables at the cell faces via a 5th-order-accurate WENO scheme [START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flow problems[END_REF] (though reconstruction of the conservative variables and lower-order WENO schemes are also supported). This allows us to apply a Riemann solver, so

F i+1/2,j,k = 1 4 2 m=1 2 l=1 F (q L i+1/2,j,k , q R i+1/2,j,k ), ( 32
)
where F is the numerical flux function of the Riemann solver. We use the HLLC approximate Riemann solver to compute the fluxes [START_REF] Toro | Restoration of the contact surface in the HLL-Riemann solver[END_REF], though other Riemann solvers are also supported.

Limiters for improved numerical stability 4.2.1. Volume fraction limiting

Following our mixture rules [START_REF] Etter | Underwater acoustic modeling and simulation[END_REF], the volume fractions are physically required to sum to unity. However, the accumulation of numerical errors can preclude this if the first N k -1 volume fractions exceed unity, and thus one of the volume fractions must be negative [START_REF] Meng | Numerical simulations of droplet aerobreakup[END_REF]. This is of course unphysical and leads to other numerical issues, such as complex speeds of sound. To treat this, we impose the volume fraction mixture rule by limiting each volume fraction as 0 ≤ α k ≤ 1 for all components k, then rescaling them as

α k = α k N k k=1 α k for k = 1, . . . , N k . ( 33
)
We note that we only use this limiting when computing the mixture properties, and do not alter them otherwise as to avoid polluting the mass conservation properties of the method.

Flux limiting

The WENO schemes we utilize are, in general, not TVD. This can be problematic when WENO cannot form a smooth stencil to reconstruct, and can lead to numerical instabilities even when the usual CFL criteria are met. MFC supports advective flux limiting to treat this issue, which improves stability, though it also increases numerical dissipation and thus smearing of material interfaces. However, we use the gradient of the local volume fraction χ to minimize this effect, which localizes the limiter to non-smooth regions of the flow. In one dimension of our dimensional-splitting procedure this yields

χ i =        α i -α i-1 α i+1 -α i if u * ≥ 0, α i+2 -α i+1 α i+1 -α i if u * < 0, ( 34
)
where i is the spatial index and u * is the local velocity computed by the Riemann solver. Ultimately, the modified flux is a combination of a low-and high-order accurate flux approximation. The low-order flux is chosen to be equivalent to a first-order WENO reconstruction and the high-order flux is the Riemann flux from the WENO reconstruction. Specifically, MFC supports the minmod, MC, ospre, superbee, Sweby, van Albada, or van Leer flux limiters, each of which is a function of the volume fraction gradient χ. Further details of our numerical implementation are located in Meng [START_REF] Meng | Numerical simulations of droplet aerobreakup[END_REF].

Cylindrical coordinate considerations

MFC also supports the use of cylindrical coordinates. They are formulated in similar fashion to [START_REF] Ball | Shock-induced collapse of a cylindrical air cavity in water: a free-Lagrange simulation[END_REF], though with the addition of an additional set of source terms on the right-hand-side associated with the 1/r cylindrical terms of the divergence operator. Our implementation was detailed by Meng [START_REF] Meng | Numerical simulations of droplet aerobreakup[END_REF] and does not use an alternative treatment of the WENO weights in the azimuthal direction, such as would be required to ensure high-order accuracy away from discontinuities. While such methods for this exist, they are generally complex and suffer from numerical stability issues [START_REF] Li | WENO schemes for cylindrical and spherical geometry[END_REF][START_REF] Mignone | High-order conservative reconstruction schemes for finite volume methods in cylindrical and spherical coordinates[END_REF][START_REF] Wang | High-order schemes for cylindrical/spherical geometries with cylindrical/spherical symmetry[END_REF]. The cylindrical coordinate treatment also requires velocity gradients, rather than just velocities, at cell boundaries. For this, we use a second-order-accurate finite-difference and averaging procedure to obtain the necessary velocity gradients at these points. Finally, the coordinate singularity at r → 0 is treated using the method of Mohseni and Colonius [START_REF] Mohseni | Numerical treatment of polar coordinate singularities[END_REF], which performs differentiation in the radial direction via a redefinition of the radial coordinate; in our implementation, we place the singularity at the finite-volume cell boundary (rather than the center). We note that this implementation requires an even number of grid cells in the azimuthal coordinate direction. A consequence of the cylindrical coordinate system is that the grid cells near the radial center are much smaller than those far away, which restricts the global CFL criterion. Following Mohseni and Colonius [START_REF] Mohseni | Numerical treatment of polar coordinate singularities[END_REF], we address this issue by using a spectral filter, which filters the high-frequency components of the solution near the centerline and relaxes the CFL criterion. Our implementation was verified by Meng [START_REF] Meng | Numerical simulations of droplet aerobreakup[END_REF] to be second-order accurate away from discontinuities via a propagating spherical pressure pulse, while the construction of the viscous stress tensor was verified using the method of manufactured solutions.

One-way acoustic wave generation

The source term s of ( 1) and ( 23) can be augmented with additional terms s s (t) for the generation of one-way acoustic waves, for examples to model an ultrasound transducer immersed in the flow. Following Maeda and Colonius [START_REF] Maeda | A source term approach for generation of one-way acoustic waves in the Euler and Navier-Stokes equations[END_REF], these take the form

s s (t) = Γ Ω Γ (ξ, t)δ(X(ξ, t) -x)dξ (35)
where Ω Γ is (possibly time dependent) forcing, Γ is the surface the forcing acts upon, δ is the Dirac delta function, and X maps the ξ coordinate to physical space x. In our numerical framework this is represented at cell I i,j,k as

s s i,j,k (t) = M m=1 Ω Γ (ξ m , t)δ h (|X(ξ m , t) -x i,j,k |)∆ξ m , (36)
where δ h is the discrete delta function operator, ∆ξ k are the sizes of the discrete patch or line the forcing is applied to, and M are the number of such patches. For example, in two dimensions the one-dimensional line operator follows as

δ h (h) = 1 2πσ 2 e -1 2 h 2 σ 2 , ( 37
)
where σ = 3∆ and ∆ is the largest mesh spacing. For a single-component, two-dimensional problem (note that there is no volume-fraction advection equation in this case) the forcing is expressed as

Ω Γ (t) = f (t)[1/c o , cos(θ), sin(θ), c 2 o /(γ -1)], (38)
where f (t) is the time-dependent pulse amplitude, c o is the speed of sound, and θ is the forcing direction as measured from the first-coordinate axis.

Time integration

Once the spatial derivatives have been approximated, ( 23) becomes a semi-discrete system of ordinary differential equations in time. We treat the temporal derivative using a Runge-Kutta time-marching scheme for the state variables. To achieve high-order accuracy and avoid spurious oscillations, we use the third-order-accurate total variation diminishing scheme of Gottlieb and Shu [START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF]:

q (1) i,j,k = q n i,j,k + ∆tL(q n i,j,k ), q (2) i,j,k = 3 4 q n i,j,k + 1 4 q (1) i,j,k + 1 4 ∆tL(q (1) i,j,k ), ( 39
)
q n+1 i,j,k = 1 3 q n i,j,k + 2 3 q (2) i,j,k + 2 3 ∆tL(q (2) i,j,k ),
where ( 1) and ( 2) are intermediate time-step stages, L represents the right hand side of ( 25), and n is the time-step index. We note that MFC also supports Runge-Kutta schemes of orders 1-5.

Pressure-relaxation procedure

The pressure-disequilibrium model ( 11) requires a pressure-relaxation procedure to converge to an 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 first-order-accurate explicit time step integration and a re-initialization procedure that ensures total energy conservation at the discrete level. After this, the disequilibrium pressures are relaxed as µ → +∞. This procedure is performed at each Runge-Kutta stage, and so there is a unique pressure at the end of each stage and the 5-and 6-equation models reconstruct the same variables. As a result, simulations of the pressure-disequilibrium model are only modestly more expensive than the 5-equation models (about 5% for spherical bubble collapses [START_REF] Schmidmayer | An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics[END_REF]).

Simulation verification and validation

We next present several test cases that validate and verify MFC's capabilities. These include one-, two-, and three-dimensional cases that span a wide variety of flow problems.

Shock-bubble interaction

We first consider a Mach 1.22 shock wave impinging on a 5 cm diameter spherical helium bubble in air. This problem was investigated via experimental methods by Haas and Sturtevant [START_REF] Haas | Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities[END_REF] and has been previously used as a validation case for multi-component flow simulations [START_REF] Terashima | A front-tracking/ghost-fluid method for fluid interfaces in compressible flows[END_REF][START_REF] Fedkiw | A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF][START_REF] Hu | An interface interaction method for compressible multifluids[END_REF]. Our simulations are performed using an axisymmetric configuration; further simulation specifications can be found in Coralic and Colonius [START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flow problems[END_REF].

Visualizations of the shock impinging the bubble and subsequent breakup and vortex ring production are shown in figure 1. We see that the simulation results qualitatively match those of the experiment. Importantly, no spurious oscillations can be seen in the numerical schlieren images, despite their sensitivity to small density differences. We quantitatively compare our results to the experiment by considering the velocity of key flow features: Haas and Sturtevant [START_REF] Haas | Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities[END_REF] measured the velocity of the incident, reflected, and transmitted shocks, as well as the up and down-stream interfaces and jets. Our simulation results are within 10% of the experiments for all cases and are generally within about 5% of the experimental means. Our results are also consistent with those computed independently via the level set [START_REF] Hejazialhosseini | High order finite volume methods on wavelet-adapted grids with local time-stepping on multicore architectures for the simulation of shock-bubble interactions[END_REF] and diffuse-interface methods [START_REF] So | Anti-diffusion interface sharpening technique for two-phase compressible flow simulations[END_REF], including the Kelvin-Helmholtz instability that develops along the bubble interface (see figure 1 (b)-(e)).

Shock-droplet/cylinder interaction

We next consider air shocks interacting with liquid media in two and three dimensions. These problems are more computationally challenging, primarily due to the larger density ratio. The first case we analyze consists of a Mach 1.47 shock impinging a 4.8 mm diameter liquid water cylinder. The simulation parameterization can be found in Meng and Colonius [START_REF] Meng | Numerical simulations of the early stages of high-speed droplet breakup[END_REF]. Figure 2 shows a visualization of experimental and our numerical results as the shock passes over the liquid cylinder. At early times it is difficult to assess the cylinder's deformation, so we instead compare the primary and secondary waves that are generated. In figure 2 (a) we see that the primary wave system, including the incident and reflected shock, have the same locations for both experimental and numerical results. The secondary wave system is generated when the Mach stems on both sides of the cylinder converge to the rear stagnation point; in figure 2 (b) we see that these also match closely.

We also consider the breakup of a spherical water droplet due to a helium shock (Mach 0.59 observed in the post-shock flow), following the experimental conditions of Theofanous et al. [START_REF] Theofanous | The physics of aerobreakup. part II. Viscous liquids[END_REF]. A full exposition of the simulation conditions can be found in Meng and Colonius [START_REF] Meng | Numerical simulation of the aerobreakup of a water droplet[END_REF]. Figure 3 shows their experimental image and volume fraction isosurfaces and sliced isocontours from our simulations. We use a small α l = 0.01 value for the isosurface of figure 3 (b) for comparison purposes since images from experiments are often obscured by the fine mist generated. While it is challenging to obtain accurate timing data from the experiments, a qualitative agreement between experiment and simulation are still observed for the shear-induced entrainment of the droplet.

Spherical bubble dynamics

Numerical simulation of cavitating spherical gas bubbles is challenging because mixture-region compressibility must be properly treated, discrete conservation must be enforced, and sphericity should be maintained in the presence of large density and pressure ratios. We consider a collapsing and rebounding air bubble in water at 10 times higher pressure as a test of the capabilities of MFC. Specific simulation specifications were presented in Schmidmayer et al. [START_REF] Schmidmayer | An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics[END_REF]; further, we include a guided description of the simulation setup for this problem in the example cases/3D sphbubcollapse directory of the MFC package (see table 1). [START_REF] Brennen | Cavitation and bubble dynamics[END_REF], then rebounds, as expected. The radius R of our simulations is computed from the gas volume by assuming the shape is nearly spherical; this closely matches the solution expected following the Keller-Miksis equation [START_REF] Keller | Bubble oscillations of large amplitude[END_REF]. We assess and compare the quality of our simulation with ECOGEN [START_REF] Schmidmayer | an open-source tool dedicated to multiphase compressible multiphysics flows[END_REF] via computation of the bubble sphericity during the collapse-rebound process. Figure 4 (b) shows this sphericity ψ of the bubble, which is defined following Wadell [START_REF] Wadell | Volume, shape, and roundness of quartz particles[END_REF] and a value of 1 indicates a spherical bubble. We see that the WENO scheme used by MFC can better maintain sphericity than a MUSCL scheme of ECOGEN formulated for the same diffuse-interface model [START_REF] Schmidmayer | An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics[END_REF], and is thus preferable for this problem.

Isentropic and Taylor-Green vortices

We next consider two-and three-dimension vortex problems as a means to verify our solutions of the flow equations. The two-dimension problem we consider is the evolution of a steady, inviscid, isentropic, ideal-gas vortex. This problem has been used previously to assess the convergence properties of high-order WENO schemes for smooth solutions to the Euler equations [START_REF] Balsara | Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy[END_REF][START_REF] Titarev | Finite-volume WENO schemes for three-dimensional conservation laws[END_REF]. We use it here to verify that we obtain high-order accuracy away from shocks and material interfaces; details regarding the numerical parameters and exact problem formulation can be found in Coralic and Colonius [START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flow problems[END_REF]. Figure 5 (a) shows the density error for both low-and high-order finite volume cell-averaging (following section 4). Since the solution should be steady, the error is computed as the deviation from the initial condition after 1 dimensionless time unit as a function of the grid size, for which N is the number of finite volumes in one spatial direction. The convergence is 2nd-and 5th-order accurate for 2nd-and 4th-order-accurate cell-averaging schemes, respectively. Thus, we conclude that for multi-dimensional problems the 4th-order-accurate cell averaging we employ is required to achieve 5th-order accuracy associated with the WENO reconstructions.

We use the three-dimensional Taylor-Green vortex problem of Brachet et al. [START_REF] Brachet | Small-scale structure of the Taylor-Green vortex[END_REF] to study the production of small length scales, including vortex stretching and dissipation. The simulation parameters again follow from Coralic and Colonius [START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flow problems[END_REF]. Figure 5 (b) shows the dimensionless dissipation rate of the kinetic energy ε in dimensionless time t, as computed over the entire computational domain. We see that the vortex stretching grows until t ≈ 5, after which the effects of viscous dissipation begin to dominate. MFC results closely match those of the direct solution computed by Brachet et al. [START_REF] Brachet | Small-scale structure of the Taylor-Green vortex[END_REF] using spectral methods.

Further verification

MFC has been verified using several other problems, including several one-dimensional test cases. For example, of great importance are the development of spurious oscillations at material interfaces, which can significantly pollute simulation quality. To determine if such oscillations appear when using our method, Coralic and Colonius [START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flow problems[END_REF] considered the advection of an isolated air-water interface at constant velocity in a periodic domain. It was shown that when conservative variables are reconstructed, the interface is corrupted by spurious oscillations and the pressures can even become negative. Whereas when primitive variables are reconstructed, their character is oscillation free for both velocity and pressure down to round off error.

It is also challenging to predict the correct position and speed of waves that emit from shockinterface interactions. While the shock-bubble interaction problem considered in section 5.1 showed that our method can approximate these quantities when comparing to experiments, it is helpful to consider a similar problem in one-dimension where exact solutions are available. Following Liu et al. [START_REF] Liu | Ghost fluid method for strong shock impacting on material interface[END_REF], Coralic and Colonius [START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flow problems[END_REF] used the methods employed by MFC to analyze a Mach 8.96 helium shock wave impinging an air interface. It was shown that the numerical results quantitatively matches the associated exact solution, correctly identifying the position and speed of all waves in the problem while avoiding any spurious oscillations. Of similar character is the gas-liquid shock tube problem of Cocchi et al. [START_REF] Cocchi | Treatment of interface problems with Godunov-type schemes[END_REF], which has been used as a model for underwater explosions. For this, Coralic and Colonius [START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flow problems[END_REF] also showed that the numerical solution matches the exact one and correctly identifies the position and speed of all waves.

Finally, the ensemble-averaged bubbly flow model introduced in section 3.3 was verified by simulating a weak acoustic pulse impinging a dilute bubble screen and comparing to the linearized bubble dynamic results of Commander and Prosperetti [START_REF] Commander | Linear pressure waves in bubbly liquids: Comparison between theory and experiments[END_REF]. We saw that the measured phase speed and acoustic attenuation, computed via the method of Ando [START_REF] Ando | Effects of polydispersity in bubbly flows[END_REF] and Bryngelson et al. [START_REF] Bryngelson | A quantitative comparison of phase-averaged models for bubbly, cavitating flows[END_REF], match the expected results. Further, Bryngelson et al. [START_REF] Bryngelson | A quantitative comparison of phase-averaged models for bubbly, cavitating flows[END_REF] showed that this method quantitatively matches the volume-averaged formulation of the same problem [START_REF] Fuster | Modelling bubble clusters in compressible liquids[END_REF][START_REF] Maeda | Eulerian-Lagrangian method for simulation of cloud cavitation[END_REF].

Illustrative examples

Shock-bubble dynamics in a vessel phantom

We demonstrate the capabilities of MFC by first considering the shock-induced collapse of a gas bubble inside a deformable vessel. This is closely related to the vascular injury that can occur during shock-wave lithotripsy treatments [START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flow problems[END_REF][START_REF] Coralic | Shock-induced collapse of a bubble inside a deformable vessel[END_REF]. This problem is particularly challenging because it involves lareg density-, pressure-, and viscosity-ratios as well as components with significantly different equation of state parameters. Specifically, we consider a 20 µm diameter air bubble centered in a 26 µm diameter cylindrical vessel filled with water and surrounded by 10% gelatin (see figure 6 (a)). The problem is initialized via a 40 MPa shock wave impinging the side of the vessel. Details on the specific numerical setup can be found in Coralic [START_REF] Coralic | Simulation of shock-induced bubble collapse with application to vascular injury in shockwave lithotripsy[END_REF]. We visualize bubble dynamics and subsequent impingement and deformation of the vessel in figure 6. We see that by t = 39 ns the bubble shape is asymmetric and the vessel contracts due to the incoming shock, after t = 46 ns the bubble surface has gained an inflection point and by t = 66 ns it impinges the vessel surface and becomes mushroom-shaped. Importantly, all surfaces remain smooth and free of spurious oscillations, despite the large density-and viscosity-ratios.

Vocalizing humpback whales

We next demonstrate the utility, flexibility, and robustness of the ensemble-averaged bubbly flow model described in section 3.3. For this, we model the humpback whale bubble-net feeding process [START_REF] Hain | Feeding behavior of the humpback whale, megaptera novaeangliae, in the western north Atlantic[END_REF]. Specifically, multiple humpback whales vocalize towards an annular bubbly region called a bubble net, which is modeled using the acoustic source terms of section 4.4 and bubbly regions described using the phase-averaged model. We show the acoustics associated with the periodic excitation of the bubble net in figure 7 for both two-and three-dimensional configurations. In both cases, we see that the impedance associated with the relatively dilute bubble net (void fraction 10 -4 ) effectively shields the core region from the vocalizations. Indeed, it is anticipated that the whales use their nets as a tool for corralling their prey into this relatively quiet, compact region. We also see that the curved material interfaces remain smooth and free of spurious oscillations.

Parallel performance benchmarks

It is important to ensure that our parallel implementation can utilize modern, large computer resources. To do this, we benchmark MFC's parallel architecture via the usual scaling and speedup tests. These tests were carried out using 288 compute nodes, each containing two 2.6 GHz six-core AMD Opteron processors. ) measures how well the computational load is balanced across the available cores by measuring solution time while fixing the grid size distributed to each core and varying the number of cores (thus changing the overall problem size), and the speedup test (figure 8 (c)) measures how the solution time increases with respect to serial computation as the number of cores varies. Thus, speedup is defined by the ratio of time cost of a parallel simulation with a certain of cores to a serial computation. The strong scaling and speedup tests are carried out on a 500 3 grid, while a constant load of 50 3 cells per core is maintained during the weak scaling test. For all tests, MFC performs very near the ideal threshold until the number of cores is 4096, at which point the results deviate slightly from ideal.

Conclusions

We presented MFC, an open-source tool capable of simulating multi-component, multi-phase, and multi-scale flows. It uses state-of-the-art diffuse-interface models, coupled with high-order interfacecapturing and Riemann solvers, to represent multi-material dynamics. MFC includes a variety of flow models and numerical methods, including spatial and temporal orders of accuracy, that are useful when considering the computational requirements of challenging open problems. It also includes options for additional physics and modeling techniques, including a sub-grid ensemble-averaged bubbly flow model.

We also described the requirements to build MFC and its design. This included external open-source software libraries that are readily available online. MFC was divided into three main components that initialize and simulate the flow, then process the exported simulation data. Each of these components is modular, and thus can be readily modified by new developers. They are coupled together via an intuitive input Python script that automatically generates the required Fortran input files and executes the software component. The exported simulation files can be readily analyzed or treated via parallel post-processing.

Finally, we presented a comprehensive set of validations, verifications, and illustrative examples. Validation was performed via comparisons to expected bubble dynamics and shock-bubble, shockdroplet, shock-water-cylinder experiments, while verification was obtained via numerical experiments involving isentropic and Taylor-Green vortices, as well as advected-and interface-interaction problems. The capabilities and fidelity of MFC were also illustrated using challenging studies of shock-bubble-viscous-vessel-wall interaction in application to shock-wave lithotripsy and acousticbubble-net dynamics in application to feeding humpback whales. A set of performance benchmarks also showed that MFC was able to perform near the ideal threshold of parallel computational efficiency.

Figure 1 :

 1 Figure 1: Comparison between (i) experimental shadowgraphs of Haas and Sturtevant [80] and (ii) numerical schlieren visualizations [84] using MFC at select times (a)-(e) as labeled. Experimental images are ©Cambridge University Press 1987.

  (a) t = 16 µs (b) t = 32 µs(i) Experiment (ii) MFC (i) Experiment (ii) MFC

Figure 2 :

 2 Figure 2: Comparison between (i) holographic interferograms [87] and (ii) numerical schlieren visualizations using MFC at select times (a) and (b) as labeled. The experimental images are reprinted from Igra and Takayama [87].

Figure 3 :

 3 Figure 3: Comparison of (a) experimental [88] and (b) numerical water droplets isosurface α l = 0.01. (c) shows isocontours of α l ranging from 0.01 to 0.99.

Figure 4 :

 4 Figure 4: Evolution of (a) the dimensionless bubble radius R and (b) its sphericity ψ. In (b) the nominal bubble shapes, represented by α l = 0.5 isosurfaces, are shown at select times for MFC (WENO) and ECOGEN (MUSCL).

Figure 4 (

 4 Figure4(a) shows the evolution of the bubble radius; it reaches a minimum near the nominal Rayleigh collapse time t c[START_REF] Brennen | Cavitation and bubble dynamics[END_REF], then rebounds, as expected. The radius R of our simulations is computed from the gas volume by assuming the shape is nearly spherical; this closely matches the solution expected following the Keller-Miksis equation[START_REF] Keller | Bubble oscillations of large amplitude[END_REF]. We assess and compare the quality of our simulation with ECOGEN[START_REF] Schmidmayer | an open-source tool dedicated to multiphase compressible multiphysics flows[END_REF] via computation of the bubble sphericity during the collapse-rebound process. Figure4(b) shows this sphericity ψ of the bubble, which is defined following Wadell[START_REF] Wadell | Volume, shape, and roundness of quartz particles[END_REF] and a value of 1 indicates a spherical bubble. We see that the WENO scheme used by MFC can better maintain sphericity than a MUSCL scheme of ECOGEN formulated for the same diffuse-interface model[START_REF] Schmidmayer | An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics[END_REF], and is thus preferable for this problem.

Figure 5 :

 5 Figure 5: (a) L∞ density error associated with an isentropic steady vortex for the cell-averaged scheme accuracy labeled. (b) Non-dimensional kinetic energy dissipation rate ε associated with the three-dimensional Taylor-Green vortex problem; the MFC solution is compared to the direct spectral solution of Brachet et al. [92].

  (a) t = 0 ns (b) t = 39 ns (c) t = 46 ns (d) t = 58 ns

Figure 6 :

 6 Figure 6: Temporal snapshots (a)-(d) that show the bubble collapse and vessel wall deformation. The shock impinges the vessel from right to left. The bubble and vessel walls are illustrated via their 0.5 volume-fraction isosurfaces.

Figure 7 :

 7 Figure 7: Visualizations of four model humpback whales vocalizing towards an annular bubble net in (a) three-and (b) two-dimensions. The acoustics are shown via isocontours of pressure and the annular bubble net is shaded.

Figure 8 :

 8 Figure 8: MFC performance benchmarks: (a) strong scaling, (b) weak scaling, and (c) speedup tests.

Figure 8

 8 Figure 8 shows parallel performance benchmarks of MFC. The strong scaling test of figure 8 (a) measures how the solution time varies for a fixed problem size as the number of computing cores

Table 1 :

 1 The 5-and 6-equation flow models can Descriptions of the files and directories included in MFC.[component] is one of pre process, simulation, or post process.

	Name	Description
	example cases/	Example and demonstration cases
	doc/	Documentation
	installers/	Package installers: Includes FFTW
	lib/	Libraries
	src/	Source code
	master scripts/ Python modules and dictionaries, optional source files
	pre process code/ Generates initial conditions and grids
	simulation code/ Flow solver
	post process code/ Processes simulation data
	tests/	Test cases to ensure software is operating as intended
	AUTHORS	List of contributors and their contact information
	CONFIGURE	Package configuration guide
	COPYRIGHT	Copyright notice
	INSTALL	Installation guide
	LICENSE	The GNU public license file
	Makefile.in	Makefile input; generally does not need to be modified
	Makefile.user	User inputs for compilation; requires attention from user
	Makefile	Targets: all (default), [component], test, clean
	RELEASE	Release notes

Table 2 :

 2 Input/output files in the case-specific directory.

Acknowledgments

The authors are grateful for the suggestions of Dr. Benedikt Dorschner when making MFC open source. This work was supported in part by multiple past grants from the US National Institute of Health (NIH), the US Office of Naval Research (ONR), and the US National Science Foundation (NSF), as well as current NIH grant number 2P01-DK043881 and ONR grant numbers N0014-17-1-2676 and N0014-18-1-2625. The computations presented here utilized the Extreme Science and Engineering Discovery Environment, which is supported under NSF grant number CTS120005. K.M. acknowledges support from the Funai Foundation for Information Technology via the Overseas Scholarship.

link to code/repository: https://mfc-caltech.github.io Operating systems under which the program has been tested: UNIX,