# A seven-equation diffused interface method for resolved multiphase flows

Achyut Panchal  Spencer H. Bryngelson  Suresh Menon  School of Aerospace
Engineering,
Georgia Institute of Technology,
Atlanta, GA, 30332.  School of Computational Science and Engineering,
Georgia Institute of Technology,
Atlanta, GA, 30332.

###### Abstract

The seven-equation model is a compressible multiphase formulation that allows
for phasic velocity and pressure disequilibrium. These equations are solved
using a diffused interface method that models resolved multiphase flows. Novel
extensions are proposed for including the effects of surface tension,
viscosity, multi-species, and reactions. The allowed non-equilibrium of
pressure in the seven-equation model provides numerical stability in strong
shocks and allows for arbitrary and independent equations of states. A
discrete equations method (DEM) models the fluxes. We show that even though
stiff pressure- and velocity-relaxation solvers have been used, they are not
needed for the DEM because the non-conservative fluxes are accurately modeled.
An interface compression scheme controls the numerical diffusion of the
interface, and its effects on the solution are discussed. Test cases are used
to validate the computational method and demonstrate its applicability. They
include multiphase shock tubes, shock propagation through a material
interface, a surface-tension-driven oscillating droplet, an accelerating
droplet in a viscous medium, and shock–detonation interacting with a deforming
droplet. Simulation results are compared against exact solutions and
experiments when possible.

##  1 Introduction

Multiphase flows are characterized by more than one phase, each with its own
material properties. Such flows are encountered in many situations, including
liquid fuel-powered devices [1, 2], solid explosives [3], pharmaceutical
sprays [4], and more. The design and engineering of these applications require
computational modeling. From this perspective, there are two sets of methods:
resolved approaches where the multiphase entities (MPE) are fully-resolved on
the computational grid [5, 6], and dispersed approaches where the MPEs are
unresolved on the grid and instead treated as point particles [7, 8, 9, 10].

This work focuses on resolved multiphase flow modeling, which can be further
segmented into sharp- or diffused-interface methods (SIM or DIM). The SIM
naturally represents materials as infinitely sharp [11, 12]. Specialized
techniques are used to advect the interface, including the level-set method
(LS) [11], the volume of fluid approach (VOF) [13], and the coupled LS–VOF
method (CLSVOF) [14]. Example methods for handling jumps at the interface are
continuum surface forces (CSF) [15] and ghost fluid methods (GFM) [16]. SIM
are particularly attractive for their ease in specifying jump conditions and
efficient interface tracking. However, they have historically been limited by
low-speed flows, low density–viscosity ratios [17], and mass-conservation
errors [11]. We note that some of these drawbacks are changing via bespoke
methods [18, 19, 20].

Diffused interface methods (DIM) are hallmarked by the Baer–Nunziato (BN)
“seven-equation” formulation of a compressible and mixture conservative flow
[21, 6]. The equations follow from an averaging procedure to the phasic
governing flow equations. The derived partial differential equation (PDE)
solves for phasic volume fractions and phase-averaged quantities. This
formulation evolves the conserved quantities on the grid, usually via
numerical methods that allow for sharp gradients [22, 23]. As a result, the
interface naturally diffuses across multiple grid cells due to numerical
diffusion [6]. Besides their conservation properties, DIMs also allow for the
creation or destruction of interfaces, common during supercavitation and
flashing [24, 25, 26, 27].

The first numerical method for a single-component inviscid seven-equation
model was developed by Saurel and Abgrall [28]. They used a two-step approach,
where the two-phase equations were first solved using a hyperbolic solver,
followed by a relaxation step to handle phase disequilibrium. After this
initial effort, reduced forms of the Baer–Nunziato equations were developed to
relieve the otherwise stiff pressure- and velocity-relaxation procedures.
These methods are broadly used today. For example, Kapila et al. [29]
developed a reduced five-equation model with single velocity and pressure, and
its use has been widespread [30, 31, 32, 33, 34]. This is partly due to the
simplicity of including additional physical phenomena like surface tension
[35], mass-transfer [27], viscous effects [36], flashing [37], multiple
components [38], and cavitation [39].

The drawback to the five-equation approach is that a single-pressure requires
establishing a mixture EOS with a non-monotonic sound speed [40]. This results
in numerical instability near strong shocks, difficulty maintaining volume
fraction-positivity, and inaccurate wave transmission through the interface
[41]. The five-equation model by Allaire et al. [30] is relatively robust, but
it cannot model the creation and destruction of the interface due to the
absence of a volume fraction source term. A six-equation model was developed
to resolve this, which used a single velocity but phase-independent pressures
[41]. Extensions of this work include phase-change effects [42, 43, 44, 45],
viscous effects [46], and arbitrary numbers of components [38, 47].

Herein, we use the more extensive seven-equation model. The seven-equation
model maintains phase-independent pressures so that it can represent with
strong shocks and large density ratios:
$\rho^{(2)}/\rho^{(1)},p^{(2)}/p^{(1)}\gg 1$. There is no need to formulate a
mixture EOS, so one can use phase-independent, arbitrary EOSs. It also retains
phase-independent velocities, which allows representation of dispersed
unresolved MPEs [8]. Considering this, the various features of SIM and DIM
variants are summarized in Table 1.

After the initial development of the seven-equation model by Saurel and
Abgrall [28], there have been relatively limited extensions and application
use-cases [48, 49, 50, 51, 42]. This appears to be due to the complexity of
formulating an efficient numerical method. In this work, we take a step
towards addressing this problem. We start with the baseline numerical approach
of Abgrall and Saurel [48] and provide novel extensions to develop a robust
computational method for the seven-equation model that can handle strong
shocks, arbitrary EOS, viscous effects, surface tension, multiple species, and
reactions [52].

Novel contributions of this work:

  * 1.

Various numerical methods have been developed for solving the BN PDE [28, 48,
49, 50, 51]. Still, the discrete equations method (DEM) [48] accurately models
non-conservative terms at the interface. For 1D Riemann problems, DEM relieves
the need to use a stiff pressure–velocity relaxation solver for pure fluid
interfaces [48]. This effect is further studied in this work for more complex
problems.

  * 2.

Surface tension is required for modeling droplet breakup and deformation. Past
works included it in the five-equation model [35] and a path-conservative
framework for the seven-equation model [51] via source terms. We include it
within the DEM approach via appropriate modifications to the Riemann solver.

  * 3.

Viscous effects are necessary for applications such as drag computation of a
deforming droplet or droplets interacting with turbulence. These have been
previously considered in the DEM framework [36] for the five-equation model.
This work extends it to the seven-equation model.

  * 4.

The developed computational framework can deal with arbitrary and phase-
independent EOS due to the pressure non-equilibrium, and this is demonstrated
by using calorically perfect gas (CPG), thermally perfect gas (TPG), stiffened
gas (SG), and Mie–Grüneisen (MG) EOS.

  * 5.

Being able to handle multi-species mixtures for reactive flows, e.g., droplets
in a detonation wave [53]. Each phase in this model includes an arbitrary
number of species and reactions.

  * 6.

The DIM interface can numerically diffuse across multiple cells. Considering
this, an interface compression scheme, initially developed for the five-
equation model [54] is modified to be used with the seven-equation model.

Table 1: Various numerical methods for modeling resolved MPEs and their features. P: Possible; NP: Not possible; PNA: Possible, not attempted. Bold text refers to cases considered in this work. |  SIM |  5-eq DIM |  6-eq DIM |  7-eq DIM
---|---|---|---|---
$\rho^{(2)}/\rho^{(1)}$ $\gg 1,\quad$ $p^{(2)}/p^{(1)}$ $\gg 1$ |  Limited [18, 19, 20, 55] |  Limited robustness [29] |  P [41] |  P [28, 48]
Arbitrary EOS |  P [16] |  NP |  PNA |  PNA
Viscous |  P [55] |  P [36] |  P [38] |  PNA
Surface tension |  P [56] |  P [35] |  P [57] |  P [51]
Flashing / trans-critical |  NP |  P [37, 39] |  P [42, 58, 43] |  P [42]
Evaporation |  P [59, 18] |  P [27] |  P [42, 58, 43] |  P [42]
Dispersed |  NP |  NP |  NP |  P [28]
Multi-species / reaction |  P [18] |  PNA |  PNA |  PNA
Interface Compression |  Not needed |  Used [54, 60] |  Used [61] |  PNA

Herein, Section 2 describes the seven-equation mathematical model, Section 3
formulates the numerical method, and Section 4 shows verification and
validation cases.

##  2 Mathematical Formulation

###  2.1 Fine-grained description

We start from a fine-grained description of a multiphase flow before BN
averaging. An indicator function for phase $\kappa=1,2$ is defined as
$I^{\kappa}(x_{\textrm{i}},t)$, which is $1$ if the location $x_{\textrm{i}}$
is in phase $\kappa$ at time $t$, and $0$ otherwise. An evolution equation for
$I^{\kappa}$ is

$\frac{\partial I^{\kappa}}{\partial t}+u^{I}_{\textrm{i}}\frac{\partial
I^{\kappa}}{\partial x_{\textrm{i}}}=0\quad(\textrm{i}=1,2,3),$ (2.1)

where $u^{I}_{\textrm{i}}$ is the interface velocity. The single-phase
Navier–Stokes equations accompany these for each phase as

$\frac{\partial U}{\partial t}+\frac{\partial F_{\textrm{i},inv}}{\partial
x_{\textrm{i}}}+\frac{\partial F_{\textrm{i},visc}}{\partial
x_{\textrm{i}}}=\dot{S}_{\kappa}\quad({\textrm{i}}=1,2,3),$ (2.2)

where,

$U=\begin{pmatrix}1,\\\ \rho,\\\ \rho u_{\textrm{j}},\\\ \rho E,\\\ \rho
Y_{\textrm{m}}\end{pmatrix},\,F_{{i},inv}=\begin{pmatrix}0,\\\ \rho
u_{\textrm{i}},\\\ \rho u_{\textrm{i}}u_{\textrm{j}},\\\ \rho
Eu_{\textrm{i}},\\\ \rho
Y_{\textrm{m}}u_{\textrm{i}}\end{pmatrix}+\begin{pmatrix}0,\\\ 0,\\\
p\delta_{\textrm{ij}},\\\ pu_{\textrm{i}},\\\
0\end{pmatrix},\,F_{{i},visc}=\begin{pmatrix}0,\\\ 0,\\\
\tau_{\textrm{i}\textrm{j}},\\\
\tau_{\textrm{i}\textrm{j}}u_{\textrm{j}}+q_{\textrm{i}},\\\ \rho
Y_{\textrm{m}}v_{\textrm{im}}\end{pmatrix}.$ (2.3)

where $U$, $F_{\textrm{i},inv}$, and $F_{\textrm{i},visc}$ are solution
variables, inviscid fluxes and viscous fluxes, respectively. The total flux is
$F_{\textrm{i}}=F_{\textrm{i},inv}+F_{\textrm{i},visc}$, and $\rho$,
$u_{\textrm{j}},\ (\textrm{j}=1,2,3)$, $p$, and $Y_{\textrm{m}},\
(\textrm{m}=1,\dots,N_{s,\kappa})$ represent density, velocity, pressure, and
mass-fractions for $\textrm{m}^{\textrm{th}}$ species, respectively. Note that
$U(0)=1$ is redundant but corresponds to a volume fraction evolution equation
after the averaging operation. We will show this later. The number of species
that are tracked, $N^{\kappa}_{s}$, can be different for different phases.

The total energy $E$ is related to the internal energy
$e=h_{0,mix}+C_{v,mix}T$ as $E=e+1/2\ u_{\textrm{i}}u_{\textrm{i}}$, where
$C_{v,mix}$ is the mixture constant volume heat capacity and $h_{0,mix}$ is
the mixture enthalpy of formation per mass. Pressure $p$ is computed as
$P^{\kappa}(\rho,e,Y_{\textrm{m}})$ via an EOS. (EOS) The mixture formation
enthalpy per mass is computed as
$h_{0,mix}=\sum_{m=1}^{N_{s}^{\kappa}}h_{0,\textrm{m}}^{\kappa}Y_{\textrm{m}}$
with $h_{0,\textrm{m}}^{\kappa}$ being the formation enthalpy per mass of
${\textrm{m}}^{\textrm{th}}$ species. The mixture heat capacities are computed
as

$\displaystyle
C_{v,mix}=\sum_{\textrm{m}=1}^{N_{s}^{\kappa}}C_{v,\textrm{m}}^{\kappa}Y_{\textrm{m}}\quad\text{and}\quad
C_{p,mix}=\sum_{\textrm{m}=1}^{N_{s}^{\kappa}}C_{p,\textrm{m}}^{\kappa}Y_{\textrm{m}}.$
(2.4)

The mixture specific enthalpy is $h_{mix}=h_{0,mix}+C_{p,mix}T$. The species-
specific heat capacities are $C_{p,\textrm{m}}^{\kappa}$ and
$C_{v,\textrm{m}}^{\kappa}$, and they can be curve-fits of the local
temperature if needed. Both phases have independent EOS that determine
$P^{\kappa}$ and thermodynamic properties $h_{0,\textrm{m}}^{\kappa}$,
$C_{p,\textrm{m}}^{\kappa}$, $C_{v,\textrm{m}}^{\kappa}$.

The viscous stress term $\tau_{\textrm{i}\textrm{j}}$ is computed from the
shear stress tensor

$\displaystyle\mathcal{T}^{\kappa}(u_{\textrm{i}})=\mu^{\kappa}\ {\partial
u_{\textrm{i}}}/{\partial x_{j}}-2/3\mu^{\kappa}\
\delta_{\textrm{i}\textrm{j}}{\partial u_{\textrm{k}}}/{\partial
x_{\textrm{k}}}.$ (2.5)

The heat flux is modeled as

$\displaystyle q_{\textrm{i}}=\lambda^{\kappa}\ {\partial T}/{\partial
x_{\textrm{i}}}+\sum_{\textrm{m}=1}^{N_{s}^{\kappa}}h_{\textrm{m}}^{\kappa}Y_{\textrm{m}}v_{\textrm{im}}$
(2.6)

where $v_{\textrm{im}}$ are species diffusive velocities, and $\mu^{\kappa}$
and $\lambda^{\kappa}$ are viscosity and thermal conductivity. Similar to the
thermodynamic properties, the transport properties $\mathcal{T}^{\kappa}$,
$\mu^{\kappa}$, $\lambda^{\kappa}$, and species diffusivities are independent
of the phases. The same form of the shear stress tensor $\mathcal{T}^{\kappa}$
is used for both the phases (corresponding to a Newtonian fluid). However, a
different form can also be used, for example, to model a solid material.

The source terms $\dot{S}$ are
$(0,0,g_{\textrm{i}},g_{\textrm{i}}u_{\textrm{i}},\dot{\omega}_{\textrm{m}})^{\text{T}}$,
where $g_{\textrm{i}}$ represents gravity and
$\dot{\omega}_{\textrm{m}}^{\kappa}$ are the chemical reaction rates. We do
not consider gravity hereon but include it for completeness. In the presence
of chemical reactions, we also have the relations

$\displaystyle\sum_{\textrm{{m}}}\dot{\omega}_{\textrm{m}}^{\kappa}=0\quad\text{and}\quad\sum_{\textrm{{m}}}Y_{\textrm{m}}v_{\textrm{mi}}=0.$
(2.7)

The kinetics that models the reaction rates
$\dot{\omega}_{\textrm{m}}^{\kappa}$ are also different and independent
between the phases.

These equations are complemented by the jump conditions across the phases as

$\left[(F_{\textrm{i}}-u^{I}_{\textrm{i}}U)\right]=\mathcal{J}^{U}_{\textrm{i}}\quad(\textrm{i}=1,2,3),$
(2.8)

where

$\mathcal{J}^{U}_{\textrm{i}}=\begin{pmatrix}0,\\\ 0,\\\ \sigma\beta
n_{\textrm{i}}n_{\textrm{j}},\\\ \sigma\beta u^{I}_{\textrm{i}},\\\
0\end{pmatrix},$ (2.9)

and $[f]$ is the jump across the interface for any arbitrary field variable
$f$, e.g., $f^{(2)}-f^{(1)}$. Here, since surface tension is considered,
$n_{i}$ is the interface surface normal pointing from the liquid ($\kappa=2$)
to the gas ($\kappa=1$) phase, and $\sigma$ and $\beta=\partial
n_{\textrm{i}}/\partial x_{\textrm{i}}$ are the surface tension and the
interface curvature, respectively. These jump conditions provide boundary
conditions at the phase interfaces.

In the absence of any mass transfer across the phase interface and the
presence of a thermal equilibrium due to heat transfer, these jump conditions
transition into

$[u_{\textrm{i}}]=0,\quad[p]=\sigma\beta-\left[\tau_{\textrm{ij}}n_{\textrm{i}}n_{\textrm{j}}\right],\quad[T]=0.$
(2.10)

One would also maintain a Gibbs free energy equilibrium across the interface
if mass transfer across the phase interface had been considered [42, 58, 44].
In the absence of mass-transfer the interface velocity $u^{I}_{\textrm{i}}$ is
$u_{\textrm{i}}$ because $[u_{\textrm{i}}]=0$.

###  2.2 Averaged equations

To derive the coarse-grained Baer–Nunziato formulation from the fine-grained
description, an averaging operator $\langle\cdot\rangle$ is applied. The
coarse-grained quantities are defined as $\mkern
1.5mu\overline{\mkern-1.5mu\alpha\vphantom{b}\mkern-1.5mu}\mkern
1.5mu^{\kappa}=\langle I^{\kappa}\rangle$ (volume fraction) and $\mkern
1.5mu\overline{\mkern-1.5muf\vphantom{b}\mkern-1.5mu}\mkern
1.5mu^{\kappa}\mkern
1.5mu\overline{\mkern-1.5mu\alpha\vphantom{b}\mkern-1.5mu}\mkern
1.5mu^{\kappa}=\langle I^{\kappa}f\rangle$ (phase-averaged quantities). Here
the operator $\overline{\cdot}$ is spatial filtering as a result of the phase-
averaging. However, since both the flow features and the interface are
considered resolved in this work, $\mkern
1.5mu\overline{\mkern-1.5muf\vphantom{b}\mkern-1.5mu}\mkern 1.5mu^{\kappa}$ is
represented as $f^{\kappa}$ from here onward. The generalized seven-equation
formulation is derived from Eq. (2.1) and Eq. (2.2), by applying
$\langle\cdot\rangle$ as

$\begin{split}\underbrace{\frac{\partial\alpha^{\kappa}{U}^{\kappa}}{\partial
t}}_{\textrm{I}}+\underbrace{\frac{\partial{\alpha}^{\kappa}{F}^{\kappa}_{\textrm{i},inv}}{\partial
x_{\textrm{i}}}}_{\textrm{II}}+\underbrace{\frac{\partial{\alpha}^{\kappa}{F}^{\kappa}_{\textrm{i},visc}}{\partial
x_{\textrm{i}}}}_{\textrm{III}}=&\underbrace{\left\langle\left(F^{\kappa}_{\textrm{i},inv}-u_{\textrm{i}}^{\kappa}U^{\kappa}\right)\frac{\partial
I^{\kappa}}{\partial x_{\textrm{i}}}\right\rangle}_{\textrm{IV}}\\\
&+\underbrace{\left\langle F^{\kappa}_{\textrm{i},visc}\frac{\partial
I^{\kappa}}{\partial
x_{\textrm{i}}}\right\rangle}_{\textrm{V}}+\underbrace{\left\langle
U^{\kappa}(u^{\kappa}_{\textrm{i}}-u^{I}_{\textrm{i}})\frac{\partial
I^{\kappa}}{\partial
x_{\textrm{i}}}\right\rangle}_{\textrm{VI}}+\underbrace{\alpha^{\kappa}\dot{S}^{\kappa}}_{\textrm{VII}}.\end{split}$
(2.11)

Here, term I is the time-derivative of flow-field variables, terms II and III
are conservative flues, terms IV and V are non-conservative terms that are
active only at the multiphase interfaces, and they are responsible for
interphase exchanges. Term VI is mass transfer across the phase boundary, and
they are neglected here but will be a part of future work. Term VII are source
terms.

Regardless of whether the size of the MPE is larger or smaller than the
filter-size of the averaging operator $\langle\cdot\rangle$ (implicitly same
as the finite-volume grid), the same Eq. (2.11) can be used for modeling
either dispersed or resolved multiphase flows. Numerical modeling of terms I,
II, III, and VII is straightforward and they remain the same for either
resolved or dispersed flows. IV and V are unclosed, and their modeling
requires a specific focus. In the absence of mass transfer, these are
typically modeled for the seven-equation model as [28, 51, 36, 42]

$\begin{gathered}\left(IV+V\right)^{(1)}=I^{(1)}+R^{(1)}_{p}+R^{(1)}_{u}+R^{(1)}_{T},\\\
I^{(1)}=\begin{pmatrix}-u_{\textrm{i}}^{I}\frac{\partial\alpha^{(1)}}{\partial
x_{\textrm{j}}},\\\ 0,\\\
((p^{I}-\sigma\beta\delta^{(1)})\delta_{\textrm{ij}}+\tau_{\textrm{ij}}^{I})\frac{\partial\alpha^{(1)}}{\partial
x_{\textrm{j}}},\\\
\left((p^{I}-\sigma\beta\delta^{(1)})u_{\textrm{i}}^{I}+\tau_{\textrm{ij}}^{I}u^{I}_{\textrm{i}}+q_{\textrm{i}}^{I}\right)\frac{\partial\alpha^{(1)}}{\partial
x_{\textrm{i}}},\\\ 0\end{pmatrix},\\\
R^{(1)}_{p}=\begin{pmatrix}\theta^{p}\Delta p,\\\ 0\\\ 0,\\\
(p^{I}-\sigma\beta\delta^{(1)})\theta^{p}\Delta p,\\\
0\end{pmatrix},R^{(1)}_{u}=\begin{pmatrix}0,\\\ 0,\\\ \theta^{u}\Delta
u_{i},\\\ u_{\textrm{i}}^{I}\theta^{u}\Delta u_{\textrm{i}}\\\
0\end{pmatrix},R^{(1)}_{T}=\begin{pmatrix}\frac{1}{\kappa^{I}}\theta^{T}\Delta
T,\\\ 0,\\\ 0,\\\ \theta^{T}\Delta T,\\\ 0\end{pmatrix}\end{gathered}$ (2.12)

where $\Delta p=p^{(1)}-p^{(2)}+\sigma\beta$, $\Delta
u_{i}=u^{(1)}_{i}-u^{(2)}_{i}$, and $\Delta T=T^{(1)}-T^{(2)}$ are differences
of pressure, velocity, and temperature between the phases. Furthermore,
$p^{I}$, $u_{{i}}^{I}$ (modeling discussed in [28]), $\tau_{{ij}}^{I}$,
$q_{{i}}^{I}$ (modeling discussed in [36]), and $\kappa^{I}$ (modeling
discussed in [42]) are interface-averaged terms within the filter
$\langle\cdot\rangle$, and since they are unknown, they typically have to be
modeled as listed in the references listed. For instance, $p^{I}$,$u^{I}$ have
been modeled in the past as $p^{(1)}$ and $u^{(2)}$ (with $\kappa=1$ being the
gas phase and $\kappa=2$ as the liquid phase) for gas-liquid flows. In this
work, however, the discrete equations method (DEM) is used [48], and its use
circumvents the need for this modeling by solving Riemann problems at each
discrete interface instead. A similar form as in Eq. (2.12) is also used for
the phase $(2)$, and it is not repeated here for brevity. As noted before,
assuming that the phase $(1)$ is a gas-phase the phase $(2)$ is a liquid-
phase, the surface tension force $\sigma\beta$ acts from the gas-phase towards
the liquid-phase, and as a result, in Eq. (2.12), $\delta^{(1)}=1$ for the
gas-phase, whereas $\delta^{(2)}=0$ for the liquid-phase [51].

The terms containing $\theta^{p}$, $\theta^{u}$, and $\theta^{T}$ are referred
to as the relaxation terms for pressure, velocity, and temperature,
respectively. The pressure relaxation results in volume exchange between the
phases, the velocity relaxation results in interphase momentum exchanges,
e.g., drag, and the temperature relaxation results in heat exchange between
the phases. For resolved interfaces, these are typically modeled as
$\theta^{p},\theta^{u},\theta^{T},\theta^{f}\rightarrow\infty$, whereas, for
dispersed phases, they are modeled using finite-time processes such as drag
and heat transfer [28]. Enforcing this condition analytically would have
resulted in a reduced five-/six-equation model [29, 41]. Past works [28, 42]
had to use a stiff relaxation solver with the seven-equation model to enforce
$\theta^{p},\theta^{u}\rightarrow\infty$. However, the use of DEM circumvented
this as it was shown for 1D canonical cases in [48]. This is further confirmed
for more complex test cases in this work (Section 4 for more details).

###  2.3 Material properties

The governing equations are accompanied by an equation of state ($P^{\kappa}$)
and material properties for each phase. For demonstration, stiffened gas (SG,
includes parameters $p_{0}^{\kappa}$ and $\gamma^{\kappa}$) [28], calorically
perfect gas (CPG, has parameter $\gamma^{\kappa}$), thermally perfect gas
(TPG, requires property curve-fits for individual species), and Mie–Grüneisen
(MG, requires properties for individual species) [62] EOS are used in this
work. These and the corresponding parameters are well-described in the
available literature, and their forms are not repeated here for brevity. The
corresponding parameters for various materials are described separately for
each test as required in Section 4.

A Sutherland transport model is used for the gas-phase viscosity
($\mu^{(1)}$). The liquid- viscosity ($\mu^{(2)}$) and surface tension
$\sigma$ are assumed to be constants, although it is possible to use curve-
fits for those as well. The thermal conductivity $\lambda^{\kappa}$ is
computed from the viscosity $\mu^{\kappa}$ using a constant Prandtl number
($Pr=0.739$) as $\lambda^{\kappa}=C_{p,mix}^{\kappa}\mu^{\kappa}/Pr$. The
species diffusion is modeled using a unity Lewis number ($Le$) assumption,
with
$Le^{\kappa}_{\textrm{m}}=\lambda^{\kappa}/\left(\rho^{\kappa}D^{\kappa}_{\textrm{m}}C_{p,\textrm{m}}^{\kappa}\right)=1$.
For the reactive simulations shown in this work, following previous work on
hydrogen-air detonations [63], a 7-step, 6-species chemical mechanism [64] is
used to model finite-rate kinetics.

##  3 Numerical Method

A three-dimensional finite-volume computational framework is used to solve Eq.
(2.11). The equations are implemented in LESLIE, a parallel multi-physics
solver, which has been used extensively in the past for multiphase reacting
flow simulations [8, 65, 66]. The numerical approach here primarily follows
published literature [48, 28, 51, 54, 36], and every detail is not repeated
for brevity. However, the following novel contributions are made, and details
about them are provided in the following sections.

  * 1.

Implementation of the DEM [48] for multi-block structured curvilinear grids.

  * 2.

Modification in HLLC [67] Riemann solver within the DEM to include surface
tension.

  * 3.

Extension of viscous effects within the DEM [36] for multiple dimensions.

  * 4.

Improved robustness of the stiff relaxation solver [28] for dealing with
arbitrary EOS.

  * 5.

Application of an interface compression scheme [54] for the seven-equation
model.

###  3.1 Baseline approach

A multi-block structured grid is used, wherein $\mathcal{C}_{ijk}$ is a
computational cell with $i$, $j$, and $k$ as the indices in the three
curvilinear directions. The corresponding cell volume is denoted as $V_{ijk}$,
and $S_{l,ijk}$ is used to indicate the cell faces. Note that these $i,\ j,\
k$ are different from the i, j, k that we used in Section 2 for the Einstein
notation. Let $(\xi,\eta,\zeta)$ are the computational directions of
increasing $i$, $j$ and $k$, and $\partial\xi/\partial x_{i}$,
$\partial\eta/\partial x_{i}$, $\partial\zeta/\partial x_{i}$ are the
corresponding grid metrics. Quantities computed at the face-centers are
denoted with sub-scripts $i+1/2$, $j+1/2$ and $k+1/2$. For example,
$S_{i+\frac{1}{2}jk}$ and ${\mathcal{N}}_{\textrm{i},i+\frac{1}{2}jk}$ are the
face-area and the face-normal, respectively, between the cells
$\mathcal{C}_{i+1jk}$ and $\mathcal{C}_{ijk}$.

A volume-integrated form of Eq. (2.11) is numerically solved on this grid. The
cell-centered quantity to be updated at each time-step is
$\left(\alpha^{\kappa}U^{\kappa}\right)_{ijk}$, and it is also referred to as
$Q_{ijk}$. The corresponding time-update is denoted as
$d\left(\alpha^{\kappa}U^{\kappa}\right)_{ijk}$ or $dQ_{ijk}$. The increment
$dQ$ (or $d\left(\alpha_{\kappa}U_{\kappa}\right)$) has contributions due to
various terms in Eq. (2.11) and they are denoted using subscripts as
$dQ_{\textrm{II}},dQ_{\textrm{III}}$, etc.

The computation is advanced in time using an explicit two-stage MacCormack
approach. The time-step $\Delta t$ is dynamically computed considering
convective, acoustic, and diffusive propagation speeds of both the phases as

$\Delta
t^{\kappa}_{ijk}=CFL\cdot{V}_{ijk}\bigg{/}\left(\frac{1}{6}\sum_{l=1}^{6}{\mid
v^{\kappa}_{l,ijk}\mid
S_{{l,ijk}}+c^{\kappa}\overline{S}_{ijk}+\frac{2\gamma^{\kappa}\nu^{\kappa}}{Pr}^{\kappa}\frac{{\overline{S}_{ijk}}^{2}}{V_{ijk}}}\right),\quad\Delta
t=\underset{ijk}{\textrm{min}}\ \textrm{min}\left(\Delta t^{(1)}_{ijk},\Delta
t^{(2)}_{ijk}\right)$ (3.1)

where, $S_{l,ijk}$ are face area and
$v^{\kappa}_{l,ijk}=\sum_{\textrm{i}=1}^{3}u^{\kappa}_{\textrm{i},ijk}n_{\textrm{i},ijk}$
is the normal velocity on the $l^{\textrm{th}}$ face of the cell
$\mathcal{C}_{ijk}$. The averaged face area for cell $\mathcal{C}_{ijk}$ is
given by $\overline{S}_{ijk}=\frac{1}{6}\sum_{l=1}^{6}S_{l,ijk}$, and
$\gamma^{\kappa}$, $\nu^{\kappa}=\mu^{\kappa}/\rho^{\kappa}$, and
$Pr^{\kappa}$ are the heat capacity, the kinematic viscosity and the Prandtl
number, respectively. The $CFL$ is chosen as $0.5$ for stability.

###  3.2 Inviscid fluxes

In the finite-volume scheme, the inviscid increments due to the terms $II$ and
$IV$ in Eq. (2.11) are provided as

$d\left(\alpha^{\kappa}U^{\kappa}\right)_{ijk,II/IV}=d\left(\alpha^{\kappa}U^{\kappa}\right)_{ijk,II/IV}^{i}+d\left(\alpha^{\kappa}U^{\kappa}\right)_{ijk,II/IV}^{j}+d\left(\alpha^{\kappa}U^{\kappa}\right)_{ijk,II/IV}^{k},$
(3.2)

where the superscripts $i,j,k$ denote contributions from each direction. The
fluxes in each direction are then computed using the DEM. Their computation is
illustrated here only for the $i$-direction, but the same forms are used in
the other directions as well.

####  3.2.1 Discrete equations method (DEM)

Instead of directly solving for the conservative and the non-conservative
terms, DEM divides the problem into three Riemann problems at each cell-face.
A schematic for this is shown in Fig. 1(a), where the three Riemann problems
at each face are shown corresponding to three possibilities of the two-phase
boundaries, i.e., (1)–(1) (R1), (1)–(2) (R2) and (2)–(2) (R3). Note that
piece-wise constant fields are shown in each cell here for illustration.
However, we go up to second-order accuracy with piece-wise linear fields for
the numerical method implemented. On each face, left- and right-interpolations
are handled using a first/second-order MUSCL approach [23], and the
interpolated values are denoted with subscripts $i+1/2,l$, $i+1/2,r$, etc.

(a) Riemann sub-problems

(b) Fluxes via DEM

Figure 1: Schematic showing the three Riemann problems (a) and the flux
computation (b) with the DEM. Piecewise constant values are shown within the
cell for illustration. However, the numerical scheme allows for second-order
accuracy via piecewise linear interpolations. The Gray/white-colored region
shows phases (2)/(1).

####  3.2.2 Riemann solution

Solutions to the Riemann sub-problems (1)–(1), (1)–(2), and (2)–(2) are now
needed. An HLLC approximate Riemann solver [67] has been used traditionally
[28, 48]. However, modifications to it are required, including the surface
tension. Note that for a (1)–(2) Riemann problem, a discontinuity in pressure
would be maintained across the contact surface such as
$p^{(2)}-p^{(1)}=p^{\sigma}$, whereas, the solution should stay unaffected for
(1)–(1) or (2)–(2) problems. A typical 1D Riemann problem is depicted in Fig.
2, which could be either of the three sub-problems R1, R2, or R3. The state
vector $\mathcal{U}=(1,\rho,\rho u_{\textrm{i}},\allowbreak\rho E,\rho
Y_{\textrm{m}})^{T}$ is set equal to either $U^{(1)}$ or $U^{(2)}$ (on both
the left and the right side of the face) depending on whether it belong to the
phase (1) or (2). In the case where the left and the right state vectors
$\mathcal{U}^{L}$ or $\mathcal{U}^{R}$ belong to different phases, i.e.,
(1)–(2) Riemann problem, the contact wave separates the phases as the problem
evolves.

$\mathcal{U}^{R}$$\mathcal{U}^{L}$$c^{*}(\mathcal{U}^{L},\mathcal{U}^{R})$$\mathcal{U}^{R*}$$\mathcal{U}^{L*}$ContactShock/expansionShock/expansion
Figure 2: A sample Riemann problem is shown here, where $\mathcal{U}^{L}$ and
$\mathcal{U}^{R}$ are the initial left- and right-states, and
$\mathcal{U}^{L*}$, $\mathcal{U}^{R*}$ are the generated intermediate states.
The waves on the extreme right/left could be a shock or an expansion wave
(solid line), and the wave in the middle is a contact wave (dashed line),
whose speed is denoted using $c^{*}$.

A sharp jump is assumed across both the left and the right shock/expansion
waves in the HLLC approximate Riemann solver [67], and these waves travel with
the speeds $c^{L}$ and $c^{R}$, respectively. These are the acoustic
propagation speeds corresponding to the state vectors $\mathcal{U}^{L}$ and
$\mathcal{U}^{R}$, respectively. Considering this, we have the following
relations across these jumps.

$\mathcal{F}^{K*}-\mathcal{F}^{K}=c^{K}(\mathcal{U}^{K*}-\mathcal{U}^{K}),\
K\in\\{L,R\\},\quad\mathcal{F}^{R*}-\mathcal{F}^{L*}+\mathcal{F}^{\sigma}=c^{*}(\mathcal{U}^{R*}-\mathcal{U}^{L*}).$
(3.3)

where $\mathcal{F}=\left(0,\rho v_{n},\rho
u_{\textrm{i}}v_{n}+pn_{\textrm{i}},\rho Ev_{n}+pv_{n},\rho
Y_{\textrm{m}}v_{n}\right)$ is the flux vector, and that is defined for each
state in the Riemann problem. As noted before,
$v_{n}=\sum_{\textrm{i}=1}^{3}u_{\textrm{i}}n_{\textrm{i}}$ is the velocity
normal to the corresponding cell-face. The term
$\mathcal{F}^{\sigma}=\left(0,0,p^{\sigma}n_{\textrm{i}},p^{\sigma}c^{*},0\right)$
is the flux due to the surface tension at the interface boundary (at the
contact wave in this case).

The pressures and the normal velocities across the interface (contact wave)
are related as $p^{L*}=p^{R*}+p^{\sigma},\quad v^{L*}_{n}=v^{R*}_{n}$, where,
$p^{\sigma}=n^{LR}\sigma\beta$ is the pressure jump across the interface due
to the surface tension. The normal vector $n^{LR}$ is $1$ if the left phase is
the liquid phase and the right phase is the gas phase, i.e., (2)–(1), $-1$ if
the left phase is the gas phase and the right phase is the liquid-phase, i.e.,
(1)–(2), and $0$ if the phases on both the side of the interface are the same,
i.e., (1)–(1) or (2)–(2). Combining these equations and following the
derivation in Toro [67] with these additional modifications due to the surface
tension force, the speed of the contact wave and the intermediate states are
given as

$c^{*}\left(\mathcal{U}^{L},\mathcal{U}^{R}\right)=\frac{p^{R}-p^{L}+p^{\sigma}+\rho^{L}v^{L}_{n}(c^{L}-v^{L}_{n})-\rho^{R}v^{R}_{n}(c^{R}-v^{R}_{n})}{\rho^{L}(c^{L}-v^{L}_{n})-\rho^{R}(c^{R}-v^{R}_{n})},\quad\mathcal{F}^{*}=\rho^{K}\left(\frac{c^{K}-v^{K}_{n}}{c^{K}-c^{*}}\right)$
(3.4)
$\mathcal{U}^{K*}\left(\mathcal{U}^{L},\mathcal{U}^{R}\right)=\begin{pmatrix}1,\\\
\mathcal{F}^{*},\\\
\mathcal{F}^{*}\left(u_{\textrm{i}}^{K}+(c^{*}-v^{K}_{n})n_{\textrm{i}}\right),\\\
\mathcal{F}^{*}E^{K}+\mathcal{F}^{*}\left(c^{*}-v^{K}_{n}\right)\left(c^{*}+\frac{p^{K}}{\rho^{K}(c^{K}-v^{K}_{n})}\right),\\\
\mathcal{F}^{*}Y_{\textrm{m}}^{K}\end{pmatrix}.$ (3.5)

Once the intermediate states are known, the conservative Riemann flux at the
cell face is computed as

$\begin{split}F^{cons}_{inv}\left(\mathcal{U}^{L},\mathcal{U}^{R}\right)=\begin{cases}\mathcal{F}(\mathcal{U}^{L})\quad\textrm{if
}0\leq c^{L}\\\ \mathcal{F}(\mathcal{U}^{L*})\quad\textrm{if }c^{L}\leq 0\leq
c^{*}\\\ \mathcal{F}(\mathcal{U}^{R*})\quad\textrm{if }c^{*}\leq 0\leq
c^{R}\\\ \mathcal{F}(\mathcal{U}^{R})\quad\textrm{if }c^{R}\leq 0\\\
\end{cases}\end{split}$ (3.6)

Note that the conventional single-phase Riemann solution is recovered by
letting $p^{\sigma}=0$.

####  3.2.3 Inviscid conservative flux (Terms $II$)

Based on the Riemann solutions that are constructed for the sub-problems
(1)–(1), (1)–(2), and (2)–(2), we can now compute the inviscid conservative
fluxes. To illustrate this process with DEM, a schematic is shown in Fig. 1
(b). For this example, during the Riemann problem evolution, the interface
that was initially at $i+1/2$ has now moved within the cell on the right based
on the speed of the contact wave $c^{*}$. $F^{cons}_{inv}$ are the
conservative fluxes for the Riemann problems (1)–(1), (1)–(2), and (2)–(2)
that act at the $i+1/2$ cell-face.

Considering the solutions from the three Riemann problems, the averaged
i-directional conservative fluxes (term II in Eq. (2.11)) at face $i+1/2$ for
phase (1) are given using the DEM as

$\begin{split}\left\langle
I^{(1)}\left(F^{(1)}_{inv}\right)\right\rangle_{i+\frac{1}{2}}=&\quad\textrm{min}\left(\alpha^{(1)}_{i+\frac{1}{2},l},\alpha^{(1)}_{i+\frac{1}{2},r}\right)F^{cons}_{inv}\left(U_{i+\frac{1}{2},l}^{(1)},U_{i+\frac{1}{2},r}^{(1)}\right)\\\
&+\textrm{max}\left(\alpha^{(1)}_{i+\frac{1}{2},l}-\alpha^{(1)}_{i+\frac{1}{2},r},0\right)\left(\beta_{i+\frac{1}{2}}^{(1,2)}\right)^{+}F^{cons}_{inv}\left(U_{i+\frac{1}{2},l}^{(1)},U_{i+\frac{1}{2},r}^{(2)}\right)\\\
&+\textrm{max}\left(\alpha^{(2)}_{i+\frac{1}{2},l}-\alpha^{(2)}_{i+\frac{1}{2},r},0\right)\left(-\beta_{i+\frac{1}{2}}^{(2,1)}\right)^{+}F^{cons}_{inv}\left(U_{i+\frac{1}{2},l}^{(2)},U_{i+\frac{1}{2},r}^{(1)}\right),\end{split}$
(3.7)

where $F^{cons}_{inv}$ are the fluxes based on the Riemann solutions as
discussed earlier, and their coefficients are the probabilities of
encountering a gas-gas, a liquid-liquid, a gas-liquid, or a liquid-gas
interface on a fine-grained discrete level. Further details and justifications
are provided in the original DEM work [48], and they are not repeated here for
brevity. Coefficients $\beta_{i+{1}/{2}}^{(p,q)}$ are computed as
$\textrm{sign}\left(c^{*}\left(U_{i+1/2,l}^{p},U_{i+1/2,r}^{q}\right)\right)$,
and $\left(\beta\right)^{+}=\mathrm{max}\left(\beta,0\right)$.

After computing the fluxes, the increment
$d\left(\alpha^{\kappa}U^{\kappa}\right)$ due to the $i-$directional
conservative inviscid flux (term II in Eq. (2.11)) is computed as

$\begin{split}Vd\left(\alpha^{\kappa}U^{\kappa}\right)_{ijk,II}^{i}=&S_{i+\frac{1}{2}jk}\left\langle
I^{\kappa}\left(F^{\kappa}_{inv}\right)\right\rangle_{i+\frac{1}{2}jk}-S_{i+\frac{1}{2}jk}\left\langle
I^{\kappa}\left(F^{\kappa}_{inv}\right)\right\rangle_{i-\frac{1}{2}jk}.\end{split}$
(3.8)

####  3.2.4 Inviscid non-conservative flux (Terms $IV$)

The non-conservative fluxes are referred to as the ones that act at a (1)-(2)
interface (contact) but within a finite-volume cell once the multiphase
interface has moved into it, e.g., the flux $F^{lag}_{L/R}$ acting on the
interface that has moved into the cell $i+1$ in Fig. 1 (b). The subscripts
$L/R$ correspond to whether the flux is acting on the material on the left or
right side of the interface. Following [48], second-order $i-$directional
update for phase $(1)$ due to the inviscid non-conservative flux (term IV in
Eq. (2.11)) is given as

$\begin{split}Vd\left(\alpha^{(1)}U^{(1)}\right)_{ijk,IV}^{i}&=\textrm{max}\left(\alpha^{(1)}_{i+\frac{1}{2},l}-\alpha^{(1)}_{i+\frac{1}{2},r},0\right)\left(-\beta_{i+\frac{1}{2}}^{(2,1)}\right)^{+}F^{lag}_{L,inv}\left(U_{i+\frac{1}{2},l}^{(1)},U_{i+\frac{1}{2},r}^{(2)}\right)S_{i+\frac{1}{2}}\\\
&-\textrm{max}\left(\alpha^{(2)}_{i+\frac{1}{2},l}-\alpha^{(2)}_{i+\frac{1}{2},r},0\right)\left(-\beta_{i+\frac{1}{2}}^{(2,1)}\right)^{+}F^{lag}_{R,inv}\left(U_{i+\frac{1}{2},l}^{(2)},U_{i+\frac{1}{2},r}^{(1)}\right)S_{i+\frac{1}{2}}\\\
&+\textrm{max}\left(\alpha^{(1)}_{i-\frac{1}{2},l}-\alpha^{(1)}_{i-\frac{1}{2},r},0\right)\left(\beta_{i+\frac{1}{2}}^{(1,2)}\right)^{+}F^{lag}_{L,,inv}\left(U_{i+\frac{1}{2},l}^{(1)},U_{i+\frac{1}{2},r}^{(2)}\right)S_{i-\frac{1}{2}}\\\
&-\textrm{max}\left(\alpha^{(2)}_{i-\frac{1}{2},l}-\alpha^{(2)}_{i-\frac{1}{2},r},0\right)\left(\beta_{i-\frac{1}{2}}^{(2,1)}\right)^{+}F^{lag}_{R,inv}\left(U_{i-\frac{1}{2},l}^{(2)},U_{i-\frac{1}{2},r}^{(1)}\right)S_{i-\frac{1}{2}}\\\
&+\textrm{max}\left(\alpha^{(1)}_{i-\frac{1}{2},r}-\alpha^{(1)}_{i+\frac{1}{2},l},0\right)F^{lag}_{L}\left(U_{i}^{(1)},U_{i}^{(2)}\right)\frac{1}{2}\left(S_{i+\frac{1}{2}}+S_{i-\frac{1}{2}}\right)\\\
&-\textrm{max}\left(\alpha^{(2)}_{i-\frac{1}{2},r}-\alpha^{(2)}_{i+\frac{1}{2},l},0\right)F^{lag}_{R}\left(U_{i}^{(2)},U_{i}^{(1)}\right)\frac{1}{2}\left(S_{i+\frac{1}{2}}+S_{i-\frac{1}{2}}\right)\\\
&+\langle
N_{int}\rangle\left(F_{R}^{lag}\left(U_{i}^{(2)},U_{i}^{(1)}\right)-F_{L}^{lag}\left(U_{i}^{(1)},U_{i}^{(2)}\right)\right)\frac{1}{2}\left(S_{i+\frac{1}{2}}+S_{i-\frac{1}{2}}\right).\end{split}$
(3.9)

where $F^{lag}_{L/R,inv}(U^{L},U^{R})$ are fluxes between the phases at the
interface, and they are computed as part of the previously defined (1)-(2)
Riemann solution as

$\begin{split}F^{lag}_{L,inv}(U^{L},U^{R})=\mathcal{F}(\mathcal{U}^{L*})-c^{*}\left(\mathcal{U}^{L},\mathcal{U}^{L}\right)\mathcal{U}^{L*},\\\
F^{lag}_{R,inv}(U^{L},U^{R})=\mathcal{F}(\mathcal{U}^{R*})-c^{*}\left(\mathcal{U}^{L},\mathcal{U}^{R}\right)\mathcal{U}^{R*}\end{split}$
(3.10)

Note that based on the approximate Riemann solution derived earlier, we have
$F^{lag}_{L}-F^{lag}_{R}=F^{\sigma}$ for flows with surface tension.

The first two terms in Eq. (3.9) correspond to the non-conservative fluxes
that act within the cell $i$ as a result of a multiphase interface that would
have entered through the cell-face $i+1/2$. The following two terms correspond
to an interface that would have entered through the cell-face $i-1/2$. The
subsequent two terms are for second-order accuracy when piecewise linear
interpolations are used within a cell. In the last term, $\langle
N_{int}\rangle$ is the average number of internal interfaces within the cell,
and this term was related to the relaxation processes within a cell, which
would lead to local equilibrium. For resolved interfaces, this could be
modeled via the stiff relaxation solvers as described later in Section 3.5.
The reader can find more details about the DEM elsewhere [48], and they are
not repeated here for brevity.

###  3.3 Viscous fluxes

Terms III and V in Eq. (2.11) are the viscous conservative and non-
conservative updates. Abgrall and Rodio [36] extended the original DEM
framework [48] for viscous fluxes. However, they only did this for 1D flows,
and here it is also used for 2D/3D. The viscous updates
$d\left(\alpha^{\kappa}U^{\kappa}\right)_{ijk,III}$ and
$d\left(\alpha_{\kappa}U_{\kappa}\right)_{ijk,V}$ follow the same formula as
Eq. (3.2) for their directional contributions from $i/j/k$. The flux
computations are described here only for the $i$-direction, but the other
directions use a similar form.

####  3.3.1 Single-phase viscous fluxes

A Riemann problem is inherently defined for 1D inviscid flows, and therefore,
it is not possible to completely replicate the inviscid computations shown
earlier for the viscous terms. The computation of single-phase viscous fluxes
$F^{\kappa}_{i,visc}$ is discussed first in this section, and their
compilation into the two-phase conservative and non-conservative terms is
described in the following sections.

Centered differencing and centered averaging are used for computing the terms
$\partial_{x_{\textrm{i}}}f\,\vline_{\ i+1/2jk}$, $f_{i+1/2jk}$, etc. The
physical derivatives $\partial_{x_{\textrm{i}}}f$, e.g.,
$\partial_{x_{\textrm{j}}}u_{\textrm{i}}$, $\partial_{x_{\textrm{i}}}T$, etc.,
are obtained from the computational derivatives as

$\frac{\partial f}{\partial x_{\textrm{i}}}=\frac{\partial
f}{\partial\zeta}\frac{\partial\zeta}{\partial x_{\textrm{i}}}+\frac{\partial
f}{\partial\eta}\frac{\partial\eta}{\partial x_{\textrm{i}}}+\frac{\partial
f}{\partial\xi}\frac{\partial\xi}{\partial x_{\textrm{i}}},$ (3.11)

where $\partial_{\zeta}f$, $\partial_{\eta}f$, and $\partial_{\xi}f$ are
computed from the flow-variables as

$\begin{split}&\quad\quad\frac{\partial f}{\partial\zeta}\vline_{\
i+\frac{1}{2}jk}=f_{i+1jk}-f_{ijk},\\\ &\frac{\partial
f}{\partial\eta}\vline_{\
i+\frac{1}{2}jk}=\frac{1}{2}\left(f_{i+\frac{1}{2}j+1k}-f_{i+\frac{1}{2}j-1k}\right),\\\
&\frac{\partial f}{\partial\xi}\vline_{\
i+\frac{1}{2}jk}=\frac{1}{2}\left(f_{i+\frac{1}{2}jk+1}-f_{i+\frac{1}{2}jk-1}\right).\end{split}$
(3.12)

In addition, any face-averaged quantities, e.g., $\alpha^{\kappa}$,
$\mu^{\kappa}$, etc., are computed as

$f_{i+\frac{1}{2}jk}=\frac{1}{2}(f_{i+1jk}+f_{ijk}).$ (3.13)

With these, the $i$-directional single-phase viscous fluxes are computed as

$\begin{split}F^{\kappa}_{visc,i+\frac{1}{2}jk}&=F^{\kappa}_{\textrm{i},visc}\left(\mu^{\kappa}_{i+\frac{1}{2}jk},\dots,\frac{\partial
u^{\kappa}_{\textrm{i}}}{\partial x_{i}}\vline\
_{i+\frac{1}{2}jk},\dots\right)\cdot
n_{\textrm{i},i+\frac{1}{2}jk}.\end{split}$ (3.14)

####  3.3.2 Viscous conservative flux (Terms $III$)

For the DEM computation of the viscous conservative fluxes, similar forms as
Eq. (3.7) and Eq. (3.8) are used. However, instead of using the
$F^{cons}_{inv}$ from the Riemann solutions as before, the terms
$F^{\kappa}_{visc}$ at the cell-faces (e.g., $i+1/2,j,k$) are now computed
using Eq. (3.14). The DEM viscous fluxes are then given as

$\begin{split}\left\langle
I^{(1)}\left(F^{(1)}_{visc}\right)\right\rangle_{i+\frac{1}{2}}=&\textrm{min}\left(\alpha^{(1)}_{i+\frac{1}{2},l},\alpha^{(1)}_{i+\frac{1}{2},r}\right)F_{visc}^{(1)}\\\
&+\textrm{max}\left(\alpha^{(1)}_{i+\frac{1}{2},l}-\alpha^{(1)}_{i+\frac{1}{2},r},0\right)\left(\beta_{i+\frac{1}{2}}^{(1,2)}\right)^{+}F_{visc}^{(1)}\\\
&+\textrm{max}\left(\alpha^{(2)}_{i+\frac{1}{2},l}-\alpha^{(2)}_{i+\frac{1}{2},r},0\right)\left(-\beta_{i+\frac{1}{2}}^{(2,1)}\right)^{+}F_{visc}^{(1)},\end{split}$
(3.15)

and a similar form is used for $\kappa=2$ as well.

####  3.3.3 Viscous non-conservative (Terms $V$)

The non-conservative viscous terms follow the same form as Eq. (3.9), however,
instead of $F^{lag}_{L/R,inv}$ coming from the 1D inviscid Riemann problems,
$F^{lag}_{L/R,visc}$ are modeled as a function of the single-phase viscous
terms $F^{\kappa}_{visc}$ computed at cell-faces (e.g., $i+1/2,j,k$) as

$\begin{gathered}F^{lag}_{L,visc}(U^{(1)}_{i+\frac{1}{2},l},U^{(2)}_{i+\frac{1}{2},l})=F^{lag}_{R,visc}(U^{(2)}_{i+\frac{1}{2},l},U^{(1)}_{i+\frac{1}{2},l})=F^{(1)}_{visc,i+\frac{1}{2}},\\\
F^{lag}_{L,visc}(U^{(1)}_{i},U^{(2)}_{i})=F^{lag}_{R,visc}(U^{(2)}_{i},U^{(1)}_{i})=\frac{S_{i+\frac{1}{2}}F^{(1)}_{visc,i+\frac{1}{2}}+S_{i-\frac{1}{2}}F^{(1)}_{visc,i-\frac{1}{2}}}{S_{i+\frac{1}{2}}+S_{i-\frac{1}{2}}}.\\\
\end{gathered}$ (3.16)

where $\kappa=1$ is assumed as the gas-phase, and therefore, it is assumed
that $F^{I}_{visc}=F^{(1)}_{visc}$. This is similar to the previous viscous
DEM approach of Abgrall and Rodio [36].

###  3.4 Curvature computation

The interface curvature $\beta$ is required to compute
$p^{\sigma}=\sigma\beta$ at the cell-centers $(i,j,k)$. Once the curvature is
computed at the cell centers, it is then interpolated at the cell faces using
the MUSCL approach in the same way as the other field quantities for flux
computation. The curvature for a field $f$ is defined as

$\beta=\frac{\partial}{\partial x_{i}}\left(\frac{\partial f}{\partial
x_{i}}\bigg{/}\bigg{|}\frac{\partial f}{\partial x_{i}}\bigg{|}\right).$
(3.17)

A previous approach for curvature computation is followed that used a vertex-
based staggering [51, 35], which is illustrated in 2D. The gradient of $f$ at
a grid vertex is first computed as

$\begin{split}\frac{\partial f}{\partial
x_{1}}_{i+\frac{1}{2},j+\frac{1}{2}}&=\frac{1}{2}\left(\frac{f_{i+1,j+1}-f_{i,j+1}}{\Delta
x}+\frac{f_{i+1,j}-f_{i,j}}{\Delta x}\right),\\\ \frac{\partial f}{\partial
x_{2}}_{i+\frac{1}{2},j+\frac{1}{2}}&=\frac{1}{2}\left(\frac{f_{i+1,j+1}-f_{i+1,j}}{\Delta
x}+\frac{f_{i,j+1}-f_{i,j}}{\Delta x}\right).\end{split}$ (3.18)

where $\Delta x$ and $\Delta y$ are the grid sizes in the $i$ and the $j$
directions. Based on this gradient, a unit vector normal to the interface is
computed as

$n_{\textrm{i},i+\frac{1}{2},j+\frac{1}{2}}=\begin{cases}\frac{\frac{\partial
f}{\partial
x_{\textrm{i}}}_{i+\frac{1}{2},j+\frac{1}{2}}}{\bigg{|}\bigg{|}\frac{\partial
f}{\partial x_{\textrm{i}}}_{i+\frac{1}{2},j+\frac{1}{2}}\bigg{|}\bigg{|}},\ \
\ \ \ \text{if}\ \ \ \ \ \bigg{|}\bigg{|}\frac{\partial f}{\partial
x_{\textrm{i}}}_{i+\frac{1}{2},j+\frac{1}{2}}\bigg{|}\bigg{|}>0,\\\ 0,\ \ \ \
\ \ \text{otherwise},\end{cases}$ (3.19)

where, $||\cdot||$ is the magnitude of the vector. The normal vector is then
interpolated at the cell-faces as

$n_{\textrm{i},i+\frac{1}{2},j}=\frac{1}{2}\left(n_{\textrm{i},i+\frac{1}{2},j+\frac{1}{2},}+n_{\textrm{i},i+\frac{1}{2},j-\frac{1}{2}}\right),$
(3.20)

and finally, the curvature is computed at the cell-centers as

$\kappa_{i,j}=\frac{n_{1,i+\frac{1}{2},j}-n_{1,i-\frac{1}{2},j,}}{\Delta
x}+\frac{n_{2,i,j+\frac{1}{2},}-n_{i,j-\frac{1}{2}}}{\Delta y}.$ (3.21)

The previous approach [51] used $f=\alpha^{(1)}$ in Eq. (3.17) as that would
result in the interface curvature at the droplet surface. However,
$\alpha^{(1)}$ varies sharply from 0 to 1 across an interface, and this can
result in numerical errors. To address this, in this work, $f=\phi$ is used,
which is defined as

$\phi=\frac{\left(\alpha^{(2)}\right)^{m}}{\left(\alpha^{(2)}\right)^{m}+\left(1-\alpha^{(2)}\right)^{m}}\quad\textrm{for}\
m<1.$ (3.22)

Shukla et al. [54] used this form of $\phi$ for their interface compression
scheme, and it is used for the curvature computation in this work. The field
$\phi$ is smoother as compared to $\alpha^{(2)}$, and $m$ is a parameter that
decides its smoothness. Note that $m=1$ would result in $\phi=\alpha^{(2)}$.
Finally, one would compute the curvature everywhere in the flow-field using
Eq. (3.21). However, to reduce any numerical noise, it is only used in the
regions where $0.99>\alpha^{(1)}>0.01$, and it is set to zero otherwise.

To verify the curvature computation, a circular drop with diameter $d_{0}=9.6$
mm is placed at the center ($x_{0}=0,y_{0}=0$) of a square domain of size
$[6d_{0},6d_{0}]$. The computational domain is uniformly divided into
$240\times 240$ cells. The volume fraction for the droplet is initialized via

$f=\sqrt{\left(\frac{x-x_{0}}{d_{0}/2}\right)^{2}+\left(\frac{y-y_{0}}{d_{0}/2}\right)^{2}},\quad\text{and}\quad\alpha^{(2)}=\frac{1}{2}\tanh\left(\frac{(f-1)d_{0}}{\delta
2\Delta x}\right)+\frac{1}{2}.$ (3.23)

where $\delta$ is the thickness over which the numerical interface is diffused
initially. The corresponding volume fraction fields for $\delta=0.5$ and
$\delta=2.5$, as well as the numerically computed curvature are shown in Fig.
3 along the center-line. The exact curvature for a droplet with diameter
$d_{0}$ would be $2/d_{0}$ (shown by the solid horizontal line in Fig. 3).
However, since the interface is diffused, if we were to compute the curvature
using Eq. (3.17) and Eq. (3.23), then it would be $2/(d_{0}f)$ along the
radius of the droplet (shown via dashed line in Fig. 3).

A comparison of the curvature computed using Eq. (3.21) and Eq. (3.22) at $m$
in Fig. 3 shows that for $\delta=2.5$, when the interface is diffused across
8-10 cells, $m=0.1,0.5,1$, provide a good match against $2/d_{0}f$. The
surface tension only acts at the sharp interface where the curvature is
$2/d_{0}$ (horizontal line). Indeed, a numerically diffused volume fraction
could result in further errors. Therefore, the curvature computation must work
for sharp interfaces and should compute the curvature as $2/d_{0}$. For
$\delta=0.5$, the interface is diffused over 2-3 cells. For this case,
choosing $m=1$ (i.e., $f=\alpha^{(2)}$, which was used in earlier works [51])
results in about $50$% errors, but with $m=0.5$, the errors are $<5$%.

Choosing $m=0.1$ leads to slightly lower errors as compared to $m=0.5$.
However, the test considered here is an ideal scenario where the volume
fraction is defined using Eq. (3.23), and as a result, $f$ is monotonic. In
reality, there could be numerical oscillations in the $\alpha^{(2)}$ field,
and those, when coupled with the high gradients of $\partial f/\partial\alpha$
near $\alpha=0,1$, could result in a non-monotonic behavior of $f$ for a lower
$m$ such as 0.1. This could lead to further errors in the computation, and
considering this, $m=0.5$ is chosen as a balance between the two effects in
this work.

Figure 3:  The curvature and the volume fraction are shown for a circular drop
along the centerline initialized as Eq. (3.23). The computed curvature using
Eq. (3.21) and Eq. (3.22) is shown for $m=0.1$, $m=0.5$ and $m=1.0$. The exact
curvature along the droplet radius is shown as $2/d_{0}f$, and the exact
droplet curvature at the interface is shown as $2/d_{0}$.

###  3.5 Relaxation process

The quantities $R_{u}^{\kappa}$, $R_{p}^{\kappa}$ and $R_{T}^{\kappa}$ in Eq.
(2.12) are relaxation terms. Physically, they represent the interface jump
conditions of Eq. (2.9) and Eq. (2.10) for a resolved interface. Since
$\theta^{p},\theta^{u},\theta^{T}\rightarrow\infty$ for the seven-equation
model, the relaxation processes are stiff, and they are applied separately
from the inviscid and the viscous updates through a Strang splitting [28].
Correspondingly, the flow solution is updated as
$Q^{n+1}=\mathcal{R}\mathcal{L}^{c}\mathcal{R}\mathcal{L}^{p}(Q^{n})$, where
$\mathcal{R}$ represents the stiff relaxation operator, $\mathcal{L}^{p}$ and
$\mathcal{L}^{c}$ correspond to the predictor and the corrector steps of the
time-update, and the superscript $n$ is used to represent the temporal
iteration number. The pressure, velocity and the temperature relaxation
operators are represented as $\mathcal{R}_{p}$, $\mathcal{R}_{u}$,
$\mathcal{R}_{T}$, respectively, and
$\mathcal{R}=\mathcal{R}_{T}\mathcal{R}_{u}\mathcal{R}_{p}$.

In terms of the numerical modeling of the interphase exchanges through the
DEM, the terms with $\langle N_{int}\rangle$ in Eq. (3.9) correspond to the
relaxation terms as shown by Abgrall and Saurel [48]. Here, $\langle
N_{int}\rangle$ represents the number of unresolved interfaces within a
finite-volume cell, and this would typically be finite for dispersed flows
where the droplets are unresolved. Still, it is supposed to be negligible for
the resolved interfaces such as those considered here. As a result of this,
while using the DEM, there should be no need to use a separate stiff
relaxation solver for resolved interfaces, and this is further confirmed in
this work (see Section 4).

Interestingly, a recent work [68] showed that it is not necessary to model the
non-conservative terms $I^{\kappa}$ in Eq. (2.12) if the relaxation terms
$R^{\kappa}$ are modeled consistently. This suggests that both the non-
conservative terms $I^{\kappa}$ and the relaxation terms $R^{\kappa}$ model
the same interphase exchange effects but differently. While Schmidmayer et al.
[68] were able to model the interphase exchanges only via $R^{\kappa}$, our
results here show that it is possible to model them only via the non-
conservative terms $I^{\kappa}$ with the use of the DEM.

However, to provide a consistent comparison against the previous works that
did use a relaxation solver [28], and because the seven-equation model would
asymptote to the widely used five-equation model in the stiff pressure and
velocity relaxation limit, the operators $\mathcal{R}_{p}$ and
$\mathcal{R}_{u}$ are implemented in the current computational framework. The
stiff temperature relaxation operator $\mathcal{R}_{T}$ is not included in
this work, so the thermal exchanges through the interface are only via the
resolved interphase exchange terms $I^{\kappa}$ in Eq. (3.9) via the DEM.

The stiff velocity relaxation operator $\mathcal{R}_{u}$ modifies the phase-
averaged momentum and energy ($\rho^{\kappa}\alpha^{\kappa}u_{i}^{\kappa}$,
$\rho^{\kappa}\alpha^{\kappa}E^{\kappa}$) while keeping the mixture momentum
and energy conserved. The algorithm for this is described elsewhere [28] and
it is not repeated here. The stiff pressure relaxation $\mathcal{R}_{p}$
results in local expansion/compression of the phases via a change in the
volume fraction $\alpha^{\kappa}$, and corresponding pressure work terms are
included in the energy $\rho^{\kappa}\alpha^{\kappa}E^{\kappa}$ (see Eq.
(3.9)). The algorithm matches that of [28], but with the following
improvements:

  * 1.

The viscous and the surface tension forces are included in the pressure jump
condition as shown in Eq. (2.10). As identified earlier [36],
$\tau_{\textrm{i}\textrm{j}}^{(1)}n_{\textrm{i}}n_{\textrm{j}}=\tau_{\textrm{i}\textrm{j}}^{(2)}n_{\textrm{i}}n_{\textrm{j}}$
is maintained by the viscous DEM, and therefore, the pressure jump that would
have to be applied at the interface via the relaxation process is
$p^{(2)}-p^{(1)}=\sigma\beta$. We can combine the volume fraction and the
energy equations, resulting in

$\frac{\partial\alpha^{(1)}\rho^{(1)}E^{(1)}}{\partial
t}=-\left(p^{(1)}-\sigma\beta\right)\frac{\partial\alpha^{(1)}}{\partial
t}\quad\text{and}\quad\frac{\partial\alpha^{(2)}\rho^{(2)}E^{(2)}}{\partial
t}=p^{(1)}\frac{\partial\alpha^{(1)}}{\partial t}.$ (3.24)

An iterative method similar to [28] is used to solve these equations and
compute an updated volume fraction $\alpha^{(1)}$ and energies $E^{\kappa}$ to
obtain $p^{(2)}-p^{(1)}=\sigma\beta$. Unlike previous works without the
surface tension, the energy exchange terms do not sum to zero, and the
difference between the two, $\sigma\beta\theta^{p}\Delta p$, is a source/sink
due to the interfacial surface tension energy.

  * 2.

Gradient descent is used as an iterative scheme. However, an adaptive stepping
makes the algorithm robust for the various EOSs considered in this work and
for a wide range of pressures. Specifically, it is used for the gradient
descent that would reduce the step size if either the pressures computed via
the EOS (either $\kappa=1,\,2$) become negative during the iterations or the
updated volume fraction goes beyond $[\epsilon,1-\epsilon]$ with
$\epsilon=10^{-10}$.

###  3.6 Interface compression

The interface numerically diffuses across multiple cells as the simulation
evolves. Several techniques have been proposed to alleviate this.
Modifications to the MUSCL face-reconstruction [69] or the limiter [61] and
applications of higher-order schemes [70, 47] have been useful for DIM.
However, the DEM used in this work is a specialized technique, and we cannot
directly apply it. Considering this, an interface compression scheme is
borrowed here [54] which works independently of the other parts of the solver.
Jain et al. [71] recently conducted a one-to-one comparison of various
interface compression techniques and concluded that the choice would often
depend on the modeling problem of interest. Their coupling with the seven-
equation model can be considered in the future, but it is not addressed here.
The employed iterative technique for interface compression requires solving
the following PDE in a pseudo-time $\tau$:

$\frac{\partial\alpha^{(1)}}{\partial\tau}=n_{\textrm{i}}\frac{\partial}{\partial
x_{\textrm{i}}}\left(\epsilon_{h}\ \vline\frac{\partial\alpha^{(1)}}{\partial
x_{\textrm{i}}}\vline-\alpha^{(1)}(1-\alpha^{(1)})\right)=n_{\textrm{i}}\frac{\partial
f_{comp}}{\partial x_{\textrm{i}}}.$ (3.25)

The right-hand side of this equation contains compression and diffusion terms
and the interface normal, which balance each other when a steady-state is
reached in the pseudo-time. This would result in a consistent thickness across
the diffused interface, which would be a function of $\epsilon_{h}$. Shukla et
al. [54] used this approach for a five-equation DIM, where they had to use an
additional correction for the mixture density. This approach is not required
for the seven-equation model as the densities for the phases are computed
separately.

In curvilinear coordinates Eq. (3.25) is written as

$\frac{\partial\alpha^{(1)}}{\partial\tau}=\left(n_{\textrm{i}}\frac{\partial\xi}{\partial
x_{\textrm{i}}}\frac{\partial
f_{comp}}{\partial\xi}+n_{\textrm{i}}\frac{\partial\eta}{\partial
x_{\textrm{i}}}\frac{\partial
f_{comp}}{\partial\eta}+n_{\textrm{i}}\frac{\partial\zeta}{\partial
x_{\textrm{i}}}\frac{\partial f_{comp}}{\partial\zeta}\right),$ (3.26)

where, $\alpha^{(1)}$ is defined at the cell-centers and updates to it in the
pseudo-time are denoted as $d\alpha_{ijk}^{\tau}$. First, $f_{comp}$ in Eq.
(3.25) is computed at the cell-faces. To compute $f_{comp,i+\frac{1}{2}jk}$ at
the cell-faces, a centered averaging is used for $\alpha^{(1)}$ following Eq.
(3.13), and the derivatives $\partial\alpha^{(1)}/\partial x_{\textrm{i}}$ are
computed using Eq. (3.11) and Eq. (3.12). Instead of computing the derivatives
$\partial\alpha^{(1)}/\partial x_{\textrm{i}}$ directly with $\alpha^{(1)}$,
for numerical stability, they are computed as $\partial\phi/\partial
x_{\textrm{i}}$ with the smooth field $\phi$ defined as previously in Eq.
(3.22) [54].

As the next step, $\partial{f_{comp}}/\partial\xi$, etc., are obtained at the
cell-centers via a centered differencing such as

$\displaystyle\frac{\partial f_{comp}}{\partial\xi}\ \vline_{\
ijk}=f_{comp,i+\frac{1}{2}jk}-f_{comp,i-\frac{1}{2}jk}.$ (3.27)

Similar forms are used for the $j$ and the $k$-directions. Last, for the
updates $d\alpha_{ijk}^{\tau}$ in Eq. (3.26), the interface normal

$\displaystyle n_{\textrm{i}}=\frac{\partial\phi/\partial
x_{\textrm{i}}}{\lvert\partial\phi/\partial x_{\textrm{i}}\rvert}$ (3.28)

is computed at the cell-centers via Eq. (3.11), and the corresponding cell-
centered derivatives are computed as ${\partial\phi}/{\partial\zeta}\
\vline_{\ i+\frac{1}{2}jk}={1}/{2}\left(\phi_{i+1jk}-\phi_{i-1jk}\right)$
(shown here representatively in the $i$-direction).

While applying the interface compression, the solution is evolved in the
pseudo-time using explicit Euler until a steady-state solution is reached. The
pseudo-time-step $d\tau$ is taken to be 0.2 for numerical stability [54]. The
compression parameter $\epsilon_{h}$ is taken as $V_{ijk}^{1/3}/2$, which
would result in the interface diffusing across 2-3 cells.

Eq. (3.25) is not in a strictly conservative form. As a result, it was
observed that when these pseudo-iterations are conducted for a very long time,
the solution initially converges to a steady-state, but later it diverges.
This is shown in Fig. 4 for example. For this test, a circular drop with
diameter $d_{0}=9.6$ mm is placed at the center ($x_{0}=0,y_{0}=0$) of a
square domain of size $[2d_{0},2d_{0}]$. The computational domain is uniformly
divided into $30\times 30$ cells. The volume fraction for the droplet is
initialized using Eq. (3.23) with $\delta=0.2$ and $\delta=2.5$. The
convergence is measured using an $L_{2}$ norm of $\alpha^{(1)}$ between the
current and the previous iteration. For both the cases, initially, the
iterations converge (in $\approx$10-30 steps), and a steady-state solution is
reached. However, after running this for a longer time, i.e. $>10^{3}$
iterations, the solution diverges. The intermediate steady-state solution
shows the interface diffused across 2–3 cells.

The solution to this problem in a standalone computation of the interface
compression technique is straightforward. We could use an absolute or a
relative tolerance limit to stop the iterations. However, while the
compression scheme is being used as part of the DIM simulations, this is not
possible. In addition to the interface compression, the flow fields are also
being updated every time step, and even if the interface compression update
was to be applied only once after every certain number (called $n^{comp}$) of
simulation time-steps, the multiphase simulations can go on for an arbitrarily
long number of iterations that could eventually cause this divergence. As a
solution to this, after every $n^{comp}$ simulation time-steps, the interface
compression scheme is iterated for at least 3 pseudo-time-steps, and the
volume fraction field is updated only if the $L_{2}$ norms during this
monotonically decrease. If not, the solution is reverted to its original state
without any interface compression updates. The numerical results show that the
interface compression did not cause any divergence by using this approach even
when we ran the simulations for an arbitrarily long time (up to $10^{7}$ flow-
field iterations). The interface compression was still active but acted only
when necessary, and it maintained the interface thickness diffused only across
2-3 cells. $n^{comp}=100$ is used in this work. The results for only a
circular drop are shown here, but the interface compression scheme has also
been validated for more complex shapes, e.g., ellipse, two lobes, etc. Those
results are not shown here for brevity.

(a) Initial $\delta=0.2$

(b) Initial $\delta=2.5$

Figure 4: L2 norms and volume fraction fields are shown with the interface
compression algorithm iterations.

###  3.7 Boundary Conditions

Seven different treatments are used in this work at the computational domain
boundaries: slip walls, no-slip walls, subsonic inflow, subsonic outflow,
supersonic inflow, supersonic outflow, periodic. For the no-slip walls, the
wall-normal and the tangential velocities of both the phases are set to zero
at the boundary face. Only the wall-normal component is set to zero for the
slip walls, whereas a zero gradient Neumann condition is applied for the
tangential component. All the other quantities, i.e., the temperature, the
species mass-fractions, and the volume fraction, are set using a zero gradient
Neumann condition at the walls.

The eigenvalues of the seven-equation system being solved here are
$\\{u^{(1)}-c^{(1)},u^{(1)},u^{(1)}+c^{(1)},u^{(2)}-c^{(2)},u^{(2)},u^{(2)}+c^{(2)},u^{I}\\}$
[28]. Correspondingly, seven waves are propagating in this hyperbolic system
of equations. The inflow/outflow treatment is dependent on the wave-
propagation directions at the boundaries. For supersonic flows, i.e.,
$u^{(1)},u^{(2)}>c^{(1)}$ and $u^{(1)},u^{(2)}>c^{(2)}$, all the waves
propagate in the same direction. Therefore, a supersonic inflow, where all the
waves enter the domain, is set using a Dirichlet condition. A supersonic
outflow, where all the waves leave the domain, is set using a zero-gradient
Neumann condition.

For subsonic inflow/outflows, where certain waves are entering the domain and
certain that are leaving the domain, characteristic-based boundary conditions
[72] are used. These were initially developed for multi-species single-phase
flows, but the same form is used here for the two-phase flows. Since the wave-
propagation speeds considered within this single-phase framework only
correspond to the gas-phase wave speeds
$\\{u^{(1)}-c^{(1)},u^{(1)},u^{(1)}+c^{(1)}\\}$, it is implicitly assumed that
the other phase properties, as well as the volume fraction, are propagating at
speed $u^{(1)}$ as if they are passive scalars. This is not a valid approach
in general. Still, it is retained here, considering that only a single phase
(typically the gas) is present at the inflow/outflow boundaries for the
configurations studied. The development of a characteristic-based boundary
condition for fully compressible two-phase flows is the subject of future
work.

##  4 Results

Various features of the seven-equation DIM framework are evaluated in this
section. Table 2 lists the tests that are considered in this work and the
features that they validate. For the seven-equation model, the phase-averaged
flow fields for both the phases (i.e., $\rho^{\kappa}$, $p^{\kappa}$, etc.)
are defined everywhere in the computational domain. For instance, even in a
region where only the phase (1) is present, $\rho^{(2)}$, $p^{(2)}$ etc. are
still mathematically defined. The corresponding volume fraction of the phase
(2) would be very small, i.e., $\alpha^{(2)}=\epsilon$, and this would result
$f^{mix}=\sum_{\kappa=1}^{2}\alpha^{\kappa}f^{\kappa}\approx f^{(1)}$, so the
phase (2) values, i.e., $\rho^{(2)}$, $p^{(2)}$, etc., aren’t expected to
affect the mixture flow-fields with $\epsilon\rightarrow 0$, but they still
have to be defined and initialized. To avoid a division by zero and to make
sure that the numerical solution does not get corrupted, the corresponding
$\epsilon$ has to be larger than the machine precision (double precision used
here, $10^{-15}$), and in this work it is chosen as $10^{-6}$ unless
explicitly mentioned otherwise.

Table 2: List of verification and validation test cases for the seven-equation DIM.  Feature |  Test-case and approach
---|---
Baseline solver |  Periodic pulse convection
- Order of scheme
- Non-disturbing condition [28, 48]
$p^{(2)}/p^{(1)},\rho^{(2)}/\rho^{(1)}>>1$ |  Multiphase shock tube
- Comparison against exact solution [48]
Arbitrary EOS & shocks |  Shock transmission through
- Air-Aluminum (CPG/SG) [73]
- Air-HMX (CPG/MIEG) interface
- Comparison against an exact solution
Surface tension |  $\Delta p$ in a circular drop [51, 35]
- Comparison against exact jump
- Analysis of parasitic currents
Oscillating elliptical droplet
- Compared against exact frequency [51, 35]
Viscous effects |  Two-phase lid-driven cavity
- Comparison against exact solution [74]
Cylindrical and deforming droplet drag
- Comparison against $C_{d}$ correlations [75]
Droplet deformation |  Shock-droplet interaction
- Comparison against experiments [76, 77, 31]
Reactive flows |  Detonation-droplet interaction demonstration

###  4.1 Convection

A 1D periodic computational domain of length $L=1$ m is used for this test.
The volume fraction $\alpha^{(1)}$ is initialized (at $t=t_{0}$) as a sine
function varying between 0 to 1 as $\alpha^{(1)}(x)=(1/4)\sin\left(x/(2\pi
L)\right)+1/2$. It is capped between $\epsilon$ and $1-\epsilon$ to prevent
any division by zero in the code. Volume fraction for the other phase is set
as $\alpha^{(2)}=1-\alpha^{(1)}$. Uniform pressure and velocity are specified
as $p^{(1)}=p^{(2)}=10^{5}$ Pa and $u^{(1)}=u^{(2)}=100$ m/s. The densities
are set as $\rho^{(1)}=1.0$ kg/m3 and $\rho^{(2)}=10^{3}$ kg/m3. The EOS for
phase (1) is specified as CPG, with $\gamma^{(1)}=1.4$, and for phase (2) it
is specified as SG with $\gamma^{(2)}=4.4$ and $p_{0}^{(2)}=6\times 10^{8}$.
The surface tension, the viscous terms, and the interface compression are
absent for this test, and a single specie is used for both the phases. Four
different grids are tested with the number of cells varying as
$N_{x}=20,\,50,\,100,\,200$ as well as the first- and the second-order
accurate schemes. Simulations are conducted both with and without stiff
relaxation.

The volume fraction pulse is allowed to convect in the periodic domain for one
flow-through time ($t=t_{f}$). The solution at the end is compared against the
initial solution ($t=t_{0}$), which serves as a ‘truth’ solution for this
case. The ‘truth’ solutions are denoted with a superscript ‘$truth$’. The
$L_{2}$ error norms for any field-variable $f$ are computed as

$||f||_{2}=\sqrt{\frac{1}{N_{x}}\sum_{i_{x}=1}^{N_{x}}\left({f(i_{x},t=t_{f})-f^{truth}(i_{x})}\right)}.$
(4.1)

As noted in earlier studies [48], a “non-disturbing” condition for this case
means that regardless of the volume fraction in the domain, since the
simulation starts with a uniform pressure/velocity, one should maintain a
uniform pressure/velocity at all times with machine accuracy. Both the
pressure and the velocity maintain uniform values with
$||p^{\kappa}||_{2}/p^{(1)}<10^{-15}$ and
$||u^{\kappa}||_{2}/u^{(1)}<10^{-15}$. The volume fraction $L_{2}$ error norms
are plotted in Fig. 5. The volume fraction suffers numerical diffusion.
However, the scheme order of accuracy is maintained as expected. The
differences in the reference and the numerical order of accuracy slopes are
less than 5% suggesting a reasonable match. The results shown here are without
using the stiff relaxation solver. Since the uniform pressure/velocity
condition is maintained between both the phases anyway, using the stiff
relaxation solver did not make any difference. Those results are not shown
here for brevity.

Figure 5: Normalized $L_{2}$ error norms for the 1D pulse convection test are
shown with different grid sizes and numerical scheme orders. The solid lines
have $1/1$ and $1/2$ slopes for reference.

###  4.2 Multiphase shock tube

The objective of this test is to evaluate the method’s ability to deal with
sharp jumps across a material interface [48]. A 1D computational domain of
length $L=1$ m is used. A phase interface is initialized at $x/L=0.8$. The
flow variables at the initial time ($t_{0}$) are set as
$\alpha^{(1)}=\epsilon$, $p^{(2)}=p^{(1)}=2\times 10^{8}$ for $x/L<0.8$, and
$\alpha^{(2)}=\epsilon$, $p^{(2)}=p^{(1)}=10^{5}$ for $x/L\geq 0.8$. The
initial velocities are $u^{(1)}=u^{(2)}=0$ everywhere. The initial densities
are set as $\rho^{(1)}=50$ and $\rho^{(2)}=1000$ everywhere. Based on this,
the pressure and the velocity equilibrium are already maintained at $t_{0}$.
The phase (1) is modeled as CPG with $\gamma^{(1)}=1.4$ to mimic compressed
air, whereas the phase (2) is modeled as SG with $\gamma^{(2)}=4.4$,
$p_{0}^{(2)}=6\times 10^{8}$ to mimic liquid water. Supersonic outflows are
used at both ends of the tube to let the waves pass without reflection.
Surface tension, viscous terms, and interface compression are absent for this
test. Four different grids with the number of cells
$N_{x}=200,\,500,\,1000,\,2000$ are tested with both the first- and the
second-order-accurate schemes.

The simulation is conducted for 0.2 ms, and the final solution is compared
against an analytical ‘truth’ solution computed via an exact Riemann solver as
shown in Fig. 6. A shock propagates on the right side of the interface (in
air), and a rarefaction travels towards the left (in water). The contact wave,
and therefore, the interface, also moves slightly towards the right. Note that
the flow variables for both phases ($\kappa=1,\,2$) are evolved for the seven-
equation model. However, for the resolved interface simulations that are
considered here, only the mixture properties computed as $f^{mix}$ are of
relevance, and only those are compared. The simulation results match the exact
solution with increasing accuracy with an increase in the grid resolution and
the order of the scheme (see Fig. 6). For instance, the normalized $L_{2}$
error norm of the pressure ($||p^{mix}||_{2}/p^{(2)}_{x/L<0.8}$) reduces from
$0.07$ to $0.03$ going from $N_{x}=200$ to $N_{x}=1000$ for the first order
scheme, and almost a similar improvement in the accuracy is obtained for
$N_{x}=200$ by going from the first to the second order scheme.

Results with and without the stiff relaxation (SR) solver match each other
closely at all grid resolutions and scheme orders. For example, the normalized
$L_{2}$ norm of the difference in $p^{mix}$ with and without the SR for
$N_{x}=200$ and the second-order-accurate scheme is 0.002. This shows that the
use of the DEM relieves the need to use a separate stiff relaxation solver.
The DEM accurately models the interphase exchanges. The pressure and the
velocity equality between the phases are obtained within the numerically
diffused interface ([48] for further details) even without enforcing it. When
a stiff pressure relaxation is applied the pressure and the velocity equality
are enforced everywhere, and therefore, $p^{mix}=p^{1}=p^{2}$,
$u^{mix}=u^{1}=u^{2}$. However, without it, $p^{1}$, $p^{2}$ and $u^{1}$,
$u^{2}$ can be different away from the numerically diffused interface, but
$f^{mix}\approx f^{\kappa}$ is still be maintained in the pure fluid regions
as $\alpha^{\kappa}\approx 1$.

Figure 6: Results for a multiphase shock tube using the seven-equation DIM are
compared against an exact solution.

###  4.3 Shock interaction with material interface

The passage of a shock through a material interface is simulated in this test.
A 1D computational domain of length $L=1$ m is used. The region $x/L<0.5$ is
filled with air ($\alpha^{(2)}=\epsilon$, CPG, $\gamma^{(1)}=1.4$), whereas,
$x/L>0.5$ is Aluminum ($\alpha^{(1)}=\epsilon$, SG, $p_{0}^{(2)}=21.13\times
10^{9}$, $\gamma^{(2)}=3.8$). A shock of strength Mach 2 is initialized in the
air at $x/L=0.3$ to propagate towards the material interface. The pre-shock
conditions are set to $p^{(1)}=p^{(2)}=10^{5}$ Pa, $u^{(1)}=u^{(2)}=0$ m/s,
and $\rho^{(1)}=1.2$ kg/m3, $\rho^{(2)}=2784$ kg/m3. The post-shock conditions
for $x/L<0.3$ are computed using the Rankine-Hugoniot relations for air as
$p^{(1)}=4.56\times 10^{5}$ Pa, $u^{(1)}=429.0$ m/s, and $\rho^{(1)}=3.211$
kg/m3.

Only pure air is present for $x/L<0.3$ in reality. However, we also need to
specify the flow fields for phase (2) for the seven-equation numerical
modeling. The different initial conditions (IC) attempt to understand their
effects on the final solution.

  * IC-1

: $\epsilon=10^{-3}$, $p^{(2)}=p^{(1)}$ and $u^{(2)}=u^{(1)}$ for $x/L<0.3$

  * IC-2

: $\epsilon=10^{-6}$, $p^{(2)}=p^{(1)}$ and $u^{(2)}=u^{(1)}$ for $x/L<0.3$

  * IC-3

: $\epsilon=10^{-9}$, $p^{(2)}=p^{(1)}$ and $u^{(2)}=u^{(1)}$ for $x/L<0.3$

  * IC-4

: $\epsilon=10^{-3}$, $p^{(2)}=10^{5}$Pa and $u^{(2)}=0$ for $x/L<0.3$.

For the first three IC, the pressure and the velocity equilibria are already
attained between the phases at $t_{0}$. However, this can be a potential
problem since the pre/post-shock conditions that are specified correspond to a
Mach 2 shock in the air, and they are not for the liquid, which still has a
volume fraction of $\epsilon$ in this region. This would create a Riemann
problem for the liquid resulting in artificial waves even before the shock
hits the material interface. IC-4 aims at solving this problem by initializing
phase (2) with uniform pressure/velocity even though phase (1) is initialized
as a shock. This means that the SR cannot be enforced for this, but there
would not be any artificial waves at the start.

The results at $t=0.3$ ms, after the shock has imparted the material
interface, are analyzed and compared against an analytical solution [73] in
Fig. 7. Depending on the impendence of the materials and the shock strength,
the imparted shock would reflect/transmit as either a shock or a rarefaction.
In this case, both the reflected and the transmitted waves are shocks. Due to
the shock impact, the material interface also moves, although for this case,
its velocity is only 0.095 m/s, and it is not noticeable from these plots.

In terms of the effect of the IC, due to the inconsistent pre/post-shock
conditions as discussed before and a higher $\epsilon$, as expected, IC-1
shows significant errors resulting from the artificial waves created in the
phase (2) near the initial shock location ($x/L=0.3$). The errors are already
present at t=0.04 ms, well before the shock hits the interface, and therefore,
the results at later times are also corrupted (not shown here). Both IC-2 and
IC-3 show a significant improvement over IC-1. Since $\epsilon$ is smaller for
them, even though the artificial waves are created in phase (2) initially,
their effects on the mixture quantities, i.e., $p^{mix}$, $u^{mix}$ seem to be
minimal. IC-4 also seems to be a valid initialization approach for the seven-
equation model and captures the shock transmission/reflection through the
material interface, even with $\epsilon=10^{-3}$ as there are no artificial
waves due to the initialization anymore. However, as noted before, this IC
cannot be used with an SR solver since it is not enforced at the $t_{0}$.

The results shown here are only for the second-order scheme. They improve
accuracy with grid refinement while comparing against the reference exact
solution (except for IC-1). Like the multiphase shock tube test earlier, the
results with and without a stiff relaxation (SR) show a similar trend and
reasonably match. Another test, where a shock initiated in the Aluminum enters
the air, was also simulated. For this, the transmitted wave is a shock,
whereas the reflected wave is a rarefaction due to the lower impendence of the
air. These results are not shown here for brevity.

Next, to demonstrate the ability of this method to handle arbitrary EOS, the
Al that was previously modeled using SG is now replaced with HMX, which is
modeled using an MG EOS with the same parameters as used in an earlier work
[62]. Unfortunately, an analytical solution to compare is not available here,
but the obtained results are shown in Fig. 8. The results look qualitatively
similar to the Al/air test, where the imparted shock resulted in a transmitted
and a reflected shock at the interface. These results are obtained with the
IC-2 approach, and they prove the robustness of this method to handle
arbitrary EOS and large pressure and density ratios.

Figure 7: The $p^{mix}$ is shown for the 1D shock interacting with Air/Al
material interface. The vertical dashed line shows the location of the
material interface. The results shown here are with the second-order scheme.
Figure 8: The $p^{mix}$ is shown for the 1D shock interacting with Air/HMX
material interface. The results are at $t=0.38$ ms. The vertical dashed line
shows the location of the material interface. The results shown here are with
the second-order scheme.

###  4.4 Surface tension effects

Several tests are conducted for validating the surface tension modeling
following Nguyen and Dumbser [51]. A 2D square computational domain of width
$L=1$ m is filled with air for the first test. A circular liquid droplet of
radius $r_{0}=0.1549$ m is placed at its center. Slip walls are used as
boundary conditions. As a result of the droplet curvature, higher pressure
should be maintained within the droplet. The objective of the first test is to
evaluate the ability of the DEM surface tension modifications to handle and
maintain this pressure jump. Instead of computing the curvature from the
volume fraction field, it is set to $\beta=1/r_{0}$ everywhere. The volume
fractions inside and outside the droplet are initialized as
$\alpha^{(1)}=\epsilon$ and $\alpha^{(2)}=\epsilon$. The velocities are
initialized as zero everywhere. The pressures are set as $p^{(1)}=10^{3}$ Pa
and $p^{(2)}=p^{(1)}+\sigma\beta$ everywhere, and this would result in a jump
in $p^{mix}$ at the droplet interface. In terms of the DEM flux computation,
this is appropriate since this means that the pressure jump is only present at
a phase (1)–(2) interface, where we included modifications for the surface
tension. Whereas the pressure is uniform across the (1)–(1) and the (2)–(2)
interfaces. The gas and the liquid phases are modeled using CPG and SG,
respectively, with $\gamma^{(1)}=1.4$, $\gamma^{(2)}=4.4$,
$p_{0}^{(2)}=6\times 10^{8}$. The surface tension is set to $\sigma=342$ N/m.
The gas-phase and the liquid-phase densities are $\rho^{(1)}=1.0$ kg/s and
$\rho^{(2)}=100$ kg/s, respectively. This test is conducted both with and
without the stiff relaxation solver. However, the pressure difference
$p^{(2)}-p^{(1)}$ is already maintained to $p_{\sigma}=\sigma\beta$ from the
start, which would mean that the relaxation solver does not have to act.

Three different grids, $200\times 200$, $400\times 400$, and $800\times 800$
and both first and the second-order schemes are tested. The simulations evolve
until a physical time of 0.085 s, and at the end of it, the $L_{2}$ error
norms for $p^{mix}$, $u^{mix}$ and $\alpha^{(2)}$ are computed considering the
initial conditions as the truth. For all the grids, for all the orders of the
scheme, with or without the stiff relaxation solver, and with or without the
interface compression, the normalized errors are maintained to machine
accuracy ($<10^{-15}$), suggesting that the DEM modifications for the surface
tension can handle the surface jump condition exactly.

For the second test, the same setup is used. However, the curvature is
computed from the volume fraction field as described in Section 3.4 instead of
setting it as a constant. The pressure within the domain is initially set as
$p^{(1)}=p^{(2)}=10^{3}$ Pa, and as a result of the surface tension force, it
is expected to rise within the droplet (ideally to $\sigma\beta=2207$ Pa). The
simulations are run for 0.034 s ($\sim$40,000 steps for the 800$\times$800
grid), and the pressure jump within the droplet at the end of it is measured
by integrating over the entire computational domain $\Omega$ as

$\Delta
p=\int_{\Omega}p^{(2)}\alpha^{(2)}\bigg{/}\int_{\Omega}\alpha^{(2)}-\int_{\Omega}p^{(1)}\alpha^{(1)}\bigg{/}\int_{\Omega}\alpha^{(1)}.$
(4.2)

This is compared against the exact pressure difference based on Young-Laplace
relation as $\Delta p^{err}=|\Delta p-\sigma\beta|$. Even though the
background flow should have been at rest, parasitic currents are observed due
to the numerical errors in the curvature computation. They are consistent with
the previous works [51]. These are quantified using the kinetic energy as

$KE^{mix}=\int_{\Omega}\left(\frac{1}{2}\alpha^{(1)}\rho^{(1)}{u^{(1)}}^{2}+\frac{1}{2}\alpha^{(2)}\rho^{(2)}{u^{(2)}}^{2}\right)\bigg{/}\int_{\Omega}\left(\alpha^{(1)}+\alpha^{(2)}\right).$
(4.3)

Values of both the $\Delta p^{err}$ and the $KE^{mix}$ are reported in Table 3
for various grids and with/without the stiff relaxation (SR) and with/without
the interface compression (comp.). Only the second-order scheme results are
shown here. However, the results with the first-order scheme are also
consistent with the conclusions here. The lowest errors are observed for the
case without the SR but with the interface compression. The errors in the
computed pressure jump are less than 2% for all the grids for this case, and
the parasitic velocities are $<0.2$ m/s. The parasitic currents and the
$p^{mix}$ within the droplet for this case are shown in Fig. 9 for reference.

Regardless of the interface compression, the errors with the SR are at least
an order of magnitude higher than those without the SR. Further insights into
this are provided below. The errors do not necessarily reduce with the grid,
and this could be due to the contrasting effects, as we would compute the
curvature more accurately for the finer grids. However, the coarser grids
would help the parasitic currents to diffuse numerically. The errors with the
interface compression are also significantly lower than those without it since
the droplet interface becomes irregular due to the parasitic currents. Still,
the interface compression helps keep it regular, resulting in a more accurate
curvature computation.

Table 3: Errors $\Delta p^{err}$ [Pa] and $KE^{mix}$ [kg m2/s2] for the circular droplet surface tension test are shown for various grids and various simulation options. Case | 200$\times$200 | 400$\times$400 | 800$\times$800
---|---|---|---
| $\Delta p^{err}$ | $KE^{mix}$ | $\Delta p^{err}$ | $KE^{mix}$ | $\Delta p^{err}$ | $KE^{mix}$
w/o SR, w/o comp. | $1.06\times 10^{2}$ | $7.99\times 10^{-6}$ | $2.64\times 10^{2}$ | $5.92\times 10^{-5}$ | $2.09\times 10^{2}$ | $4.35\times 10^{-5}$
w SR, w/o comp. | $8.84\times 10^{2}$ | $1.58\times 10^{-4}$ | $1.18\times 10^{3}$ | $1.56\times 10^{-3}$ | $1.25\times 10^{3}$ | $2.18\times 10^{-3}$
w/o SR, w comp. | $3.52\times 10^{1}$ | $2.26\times 10^{-6}$ | $1.49\times 10^{1}$ | $2.06\times 10^{-6}$ | $5.86\times 10^{1}$ | $2.10\times 10^{-6}$
w SR, w comp. | $1.44\times 10^{2}$ | $1.68\times 10^{-4}$ | $1.01\times 10^{3}$ | $4.59\times 10^{-4}$ | $4.66\times 10^{2}$ | $1.49\times 10^{-3}$

To understand the reasons for the large errors with the SR, $p^{\kappa}$,
$\rho^{\kappa}$, and $\alpha^{(2)}$ are plotted along the centerline for the
800$\times$800 test with and without the SR in Fig. 10. As noted before, the
DEM modifications in the surface tension are designed so that a pressure jump
is handled across the (1)-(2) interface, whereas the (1)–(1) and the (2)–(2)
interfaces should behave as before. As expected, the case with the SR
maintains pressure equality away from the interface and enforces a pressure
jump between $p^{(1)}$ and $p^{(2)}$ within the numerically diffused
interface. However, because of this, pressure jumps are present in both
$p^{(1)}$ and $p^{(2)}$ near the boundaries of the diffused interface. They
result in pressure jumps at the (1)–(1) and the (2)–(2) interfaces within the
DEM, possibly causing these larger errors.

On the other hand, for the case without the SR, $p^{(1)}$ stays constant, and
only the $p^{(2)}$ changes within the numerically diffused interface. The
liquid pressure $p^{(2)}$ reduces almost to the same value as $p^{(1)}$ right
before the interface ($x/d_{0}>0.5$) and then increases back due to the
surface tension terms at the interface. The initial reduction in it is still
in the gas phase (i.e., $\alpha^{(2)}$ is small), and therefore, it does not
adversely affect the overall mixture ($mix$) flow fields. The surface tension
is accurately captured without the SR within the numerically diffused
interface due to the DEM modifications in the (1)–(2) Riemann problem. The
errors mentioned above due to the SR are also apparent in the density fields,
as $\rho^{(1)}$ in the presence of the SR increases by almost $100$% within
the droplet. In contrast, the oscillations in both $\rho^{(1)}$ and
$\rho^{(2)}$ are $<0.5$% for the case without the SR.

Figure 9: Velocity vectors and $p^{mix}$ contours are shown for the
800$\times$800 cylindrical droplet case without the stiff relaxation and
interface compression at t=0.085 s. The droplet interface $\alpha^{(2)}=0.5$
is shown via a dashed line. Figure 10: $p^{(1)}$, $p^{(2)}$, $\rho^{(1)}$,
$\rho^{(1)}$, and $\alpha^{(2)}$ are plotted along the centerline for the
800$\times$800 test with and without the SR.

For the third test, instead of starting from a circular droplet, an elliptical
drop is used with the following shape [35, 51]

$\frac{(x-x_{0})^{2}}{0.2^{2}}+\frac{(y-y_{0})^{2}}{0.12^{2}}=1.$ (4.4)

The ellipsoidal droplet starts converting into a circular shape because of the
initial higher/lower curvature on the horizontal/vertical ends of the droplet.
The potential energy stored in the droplet’s surface tension converts into
kinetic energy. The droplet oscillates from an initial ellipse into an
intermediate circular shape and back to an ellipse. The kinetic energy
oscillates during this process, and we can compute its frequency against an
analytical value based on the Rayleigh formula [78]

$\omega=(o^{3}-o)\frac{\sigma}{(\rho^{(1)}+\rho^{(2)})r_{0}^{3}}.$ (4.5)

where $o$ is the oscillation mode. Corresponding to $o=2$ and $r_{0}=0.1549$
m, an oscillation frequency is determined as $t_{osc}=2\pi/\omega=0.085$ s.

Based on the conclusions from the previous test, this test is simulated
without the SR. The evolution of the droplet shape at various times is shown
Fig. 11, and corresponding $KE^{mix}$ is shown in Fig. 12. Similar to the
previous test, interface compression is necessary for accurate curvature and
surface tension computation. As expected, the initial elliptical droplet
converts to a circular droplet by $t/t_{osc}=0.3$. This is accompanied by an
increase in the $KE^{mix}$. The droplet regains an elliptical shape by
$t/t_{osc}=0.5$ and $KE^{mix}$ shows reduces to zero. Ideally, this process
should continue in the absence of any viscous effects. However, because of the
numerical diffusion here, the peaks of $KE^{mix}$ reduce with time. The
400$\times$400 grid only captures two peaks of $KE^{mix}$, whereas the finer
grids can capture more. The computed oscillation frequency matches the exact
value with $<$15% accuracy.

Figure 11: Oscillating ellipsoidal droplet at various times for the
800$\times$800 grid without the interface compression. Figure 12: $KE^{mix}$
[kg m2/s2] for an oscillating ellipsoidal droplet. The dashed vertical lines
correspond to the exact oscillation period.

###  4.5 Viscous effects

To verify the implementation of the viscous terms, a lid-driven cavity with a
Reynolds number (Re) of 100 is simulated first. Initially, the fluid is at
rest in a square domain of size $L$. The top wall is accelerated to a velocity
$u^{\infty}$ at $t=0$, and all the other walls stay at rest. The flow within
the domain starts moving and rotating due to the no-slip walls and the viscous
effects. The simulations are run until a steady-state solution is achieved.
Single-phase simulations [79] for this case have been established before, and
the results match with an exact solution [74]. Two-phase results are obtained
here with the DIM, where both the phases are set as CPG and with the same
properties. The volume fraction $\alpha^{(2)}$ is initialized corresponding to
a circular droplet within the domain using Eq. (3.23). Both the single-phase
and the two-phase simulation results with and without the SR match the exact
solution as shown in Fig. 13, verifying the implementation of the viscous
terms. These results are with the second-order scheme.

Figure 13: Steady-state velocities are compared against an exact solution for
the lid-driven cavity test.

Next, to validate the viscous effects in two-phase flows, the drag of a
droplet in a viscous medium is computed. A 2D domain of length $L_{x}=15r_{0}$
and width $L_{y}=10r_{0}$ is used, where $r_{0}=0.5$ mm is the radius of an
initially circular droplet initially placed at $(x,y)=(4r_{0},5r_{0})$. The
gas and the liquid phases are modeled using CPG and SG, respectively, with
$\gamma^{(1)}=1.4$, $\gamma^{(2)}=4.4$, $p_{0}^{(2)}=6\times 10^{8}$. The
velocity field surrounding the droplet is initiated using a potential flow
solution, using a stream function

$\psi=U_{0}r\rm{sin}(\theta)[1-R_{0}^{2}/r^{2}],\quad
u^{(1)}_{1}=u^{(2)}_{1}=-\partial\psi/\partial y,\quad
u^{(1)}_{2}=u^{(2)}_{2}=\partial\psi/\partial x.$ (4.6)

where $U_{0}$ is the free-stream velocity, and $\theta=1^{0}$ is used to
trigger vortex shedding. The velocities inside the droplet are initialized as
zero. The dynamic viscosity for the gas phase is set as $\mu^{(1)}=1.8\times
10^{-5}$ kg/(m-s), and for the liquid it is $\mu^{(1)}=1.14\times 10^{-3}$
kg/(m-s). The surface tension is specified as $\sigma=0.0728$ N/m. The
Reynolds number ($Re$) is defined as $Re=2\rho^{(1)}r_{0}U_{0}/\mu_{g}$, and
the Weber number ($We$) is defined as $We=2\rho^{(1)}r_{0}U_{0}^{2}/\sigma$.

Two conditions are simulated following a previous work using the five-equation
model [80]. One with $U_{0}=5$ m/s, that corresponds to $Re=81.3$ and
$We=0.1$, and another with $U_{0}=50$ m/s, that corresponds to $Re=813$ and
$We=10$. A time-scale for normalization is defined as $t_{0}=r_{0}/U_{0}$.
Three grids are used with the number of cells across the droplet diameter
$N_{x}=20,50,100$. The simulations are conducted with and without the SR and
with and without the interface compression. The interface diffuses irregularly
without interface compression, similar to the previous surface tension test.
The boundary and shear layers that are supposed to grow surrounding the sharp
gas-liquid interface are wrongly predicted. Therefore, only the results with
the interface compression are discussed further in this section.

The flow fields for the two simulated conditions are shown in Fig. 14. The
droplet at $We=0.1$ maintains a circular shape, but it deforms for $We=10$.
The low $Re$ case develops a separation bubble behind the droplet, and the
flow field remains relatively steady. On the other hand, the high $Re$ case
demonstrates a vortex shedding process. The positive and negative velocities
surrounding the droplet are due to the viscous effects. For further
processing, the droplet velocity and a drag coefficient $C_{D}$ are computed
as

$U_{c}={\int_{\Omega}\alpha^{(2)}u^{(2)}_{1}}\bigg{/}{\int_{\Omega}\alpha^{(2)}},\quad
C_{D}=\frac{\rho^{(2)}}{\rho^{(1)}}\frac{dU_{c}}{dt}\frac{\pi
r_{0}}{(U_{0}-U_{c})^{2}}.$ (4.7)

The computed $C_{D}$ is compared against a steady-state drag correlation for a
rigid cylinder $C_{D,S}(Re)$ [75] for $We=0.1$ in Fig. 15. The droplet with
$We=10$ starts deforming, and therefore, a comparison against the correlation
is not possible for that. These results show that establishing the viscous
flow surrounding the droplet takes approximately $5t_{0}$, and $C_{D}$ reaches
a steady value after that. The $C_{D}$ is overpredicted for the coarser grids.
However, a better match against the correlation $C_{D,S}(Re)$ is achieved with
grid-refinement. The $C_{D}$ for the finest grid is also $\approx$10% higher
than the correlation $C_{D,S}(Re)$, but this could be attributed to the
relative acceleration effect [81] that is not included in the correlation
$C_{D,S}(Re)$, but are present in the current simulations since the relative
velocity $U_{r}=U_{c}-U_{0}$ decreases with time. These results also show that
even though the $C_{D}$ with and without the SR do not match exactly, they are
very similar, and this further indicates the efficacy of the viscous DEM in
accurately modeling the non-conservative interface terms such that there is no
need to use the SR with the seven-equation model.

(a) $Re=81.3$, $We=0.1$, $t/t_{0}=15$

(b) $Re=813$, $We=10$, $t/t_{0}=39$

Figure 14: Streamlines and vorticity contours are shown surrounding a droplet
for the viscous computations. These results are without the SR with the
interface compression and for the grid with $N_{x}=100$ Figure 15: Computed
drag coefficient $C_{D}$ and its evolution with time is shown for $Re=81.3$,
$We=0.1$.

###  4.6 Shock-droplet interaction

This test simulates the early-stage deformation of a cylindrical water droplet
when hit by a shock wave. The setup follows the experiments by Igra et al.
[76] and other numerical works that used a five-equation DIM [77, 31]. A
circular droplet of diameter ($d_{0}$) 4.8 mm is placed in a 2D computational
domain of size 23$d_{0}$ $\times$ 12$d_{0}$. The background gas-phase is
initialized with a shock of strength $M=1.47$ that is located at a distance
$d_{0}$ ahead of the droplet center The pre-shock (background) pressure and
velocity are $p^{(1)}=p^{(2)}=101325$ Pa and $u^{(1)}=u^{(2)}=0$ m/s,
respectively. The corresponding post-shock conditions are subsonic, and they
are set as $p^{(1)}=p^{(2)}=238558$ Pa and $u^{(1)}=u^{(2)}=u_{0}=226$ m/s.
The water density is set as $\rho^{(2)}=1000$ m/s. The pre- and the post-shock
gas-phase densities are $\rho^{(1)}=1.204$ kg/m3 and $\rho^{(1)}=2.18$ kg/m3,
respectively. Note that these initial conditions are equivalent to the IC-2 as
described earlier in Section 4.3. The liquid water is modeled with a SG EOS
with $\gamma^{(2)}=6.12$ and $p_{0}^{(2)}=343.44\times 10^{6}$ Pa. The gas
phase is modeled with a CPG EOS with $\gamma^{(1)}=1.4$. The volume fraction
is initialized using Eq. (3.23) with $\delta=0.5$.

The computational domain is discretized using a uniform grid
$N_{x}=17,\,42,\,83,\,170,\,255,\,340$ cells across the droplet diameter. The
top and the bottom boundaries are set as slip walls. The left boundary behind
the shock is set as a non-reflective subsonic inflow at the post-shock
velocity $u_{0}$. The right boundary, ahead of the shock, is set as a
supersonic outflow to let the incoming shock wave pass through without any
reflections. The Reynolds number ($Re$) and Weber number ($We$) corresponding
to this condition are $Re=4860$ and $We=72320$, respectively. Under these
conditions, the droplet deformation is primarily driven by inertia during the
early stages. Therefore, the surface tension and viscous effects can be
neglected, as confirmed in the earlier numerical studies [77, 31]. Note that
these effects would still be relevant near the final stages of the breakup and
at much smaller length scales that are not resolved here. The simulations are
conducted and compared with and without the surface tension and with and
without the interface compression.

A non-dimensional time is defined as
$\tau=t/\left(u_{0}/d_{0}\sqrt{\rho^{(1)}/\rho^{(2)}}\right)$, where the
denominator characterizes time for breakup via Rayleigh–Taylor or
Kelvin–Helmholtz instabilities [82]. The simulation results at various times
are shown in Fig. 16 for the case with $N_{x}=83$. These are with the
interface compression and with/without the stiff relaxation. Pressure contours
and Numerical Schliren ($|\partial\rho^{(1)}/\partial x_{\textrm{i}}|$) are
used to identify the shock structure. The droplet is using a
$\alpha^{(2)}=0.5$ contour. As the initial shock hits the droplet, there are
transmitted and reflected waves. Unlike the previous 1D case, the shock
structures are no longer planar because of the circular droplet. The waves
that propagate inside the cylinder transmit and reflect through the material
interface. The droplet remains circular during the earlier stages (up to
$\tau\approx 0.17$). However, it starts to deform due to the vortical
structures created around it later. These features are qualitatively similar
to the experimental observations [76] and previous numerical results [77, 31].
The results with and without the SR match closely, reconfirming the ability of
the DEM to alleviate the need to use an SR.

The experiments measured the droplet motion and the droplet deformation in
four geometrical parameters. 1) Leading edge drift ($x_{L}$), that is, the
x-location of the droplet leading edge measured with respect to its initial
location, 2) droplet centerline width ($w$), which measures the width of the
cylinder along the $x$-axis, 3) droplet spanwise diameter ($d$), that is the
width of the droplet in the $y$-direction, and 4) droplet area ($A$). All
these are computed from the simulation flow-fields using a threshold volume
fraction $\alpha^{(2)}_{T}$, considering that $\alpha^{(2)}>\alpha^{(2)}_{T}$
is the droplet shape whose geometrical parameters have to be determined.

Because of the uncertainty in the experimental measurements, it is unclear
what value $\alpha^{(2)}_{T}$ should take to best match the data. Even though
the gas-liquid interface is supposed to remain sharp in reality, this region
may appear diffused depending on the measurement technique because of the
small droplets strip surrounding the parent droplet surface. In addition to
this, the interface also numerically diffuses for the simulations. The results
are analyzed for $\alpha^{(2)}_{T}=0.1,\,0.5,\,0.9$, but only those with
$\alpha^{(2)}_{T}=0.9$ are shown here. Sensitivity of five-equation DIM
results in the choice of $\alpha^{(2)}_{T}$ can be found elsewhere [31] and
our observations are similar for the cases without the interface compression.
The simulation results with the interface compression do not show a
sensitivity to $\alpha^{(2)}_{T}$ since a regular interface thickness is
maintained. However, those results are considered erroneous for this test, and
further details are discussed later.

The simulations with the SR and the interface compression take longer to run,
e.g., the simulations are 1.6 $\times$ slower with the interface compression
and 5.6-times slower with the SR. Because of the cost of these simulations,
the finer grids are not simulated until completion with the SR and the
interface compression. As discussed next, the results are still considered
sufficient to draw valuable conclusions.

The droplet leading edge ($x_{L}/d_{0}$) evolution is shown in Fig. 17. All
grid resolutions show a similar trend, but a grid convergence is achieved for
$N_{x}\geq 170$. Stiff relaxation or interface compression does not result in
a noticeable difference in leading-edge predictions. The results match the
experimental data for $\tau<0.7$, however, the difference eventually grows to
almost 100% by $\tau=1.3$. This difference could be due to the subsonic inflow
boundary conditions used on the left boundary for the simulations to maintain
the post-shock velocity. The experiments were conducted within a longer shock
tube where We did not need this.

The centerline droplet width ($w$) decreases with time due to the droplet
compression, and the spanwise droplet diameter ($d$) increases as the droplet
gets elongated. These are plotted in Fig. 18 and Fig. 19, respectively.
Without the SR and the interface compression, the predicted droplet width $w$
matches well against the data ($<$5% error) for all grid resolutions. On the
other hand, the droplet height $d$ is grossly underpredicted at coarser
resolutions, and a grid convergence is obtained only for $N_{x}\geq 255$. The
fine grid results capture the same trend as the experimental measurements
[76]. However, they show about 30-40% overprediction. Similar errors have also
been observed in earlier numerical studies using the five-equation DIM, but
the reasons are unclear. The droplet area is expected to reduce with time for
the experiments due to the stripping of small-sized droplets. The simulations
capture this without the SR and the interface compression as shown in Fig. 20.
The coarsest grid shows up to 25% underprediction. However, a grid convergence
is obtained for $N_{x}\geq 170$, and the finer grids’ results match the
measurements with ¡5% error.

Regarding the stiff relaxation solver usage, similar to the leading-edge
droplet predictions, the predictions for all three geometrical features ($d$,
$w$, $A$) are similar with/without the SR. The most significant differences
are observed later in $d$ and $A$ predictions, which are still $<$10-15%. This
further confirms the efficacy of the DEM in alleviating the need to use a
separate SR solver, which, as demonstrated for this case, could be as much as
5.6-times as costly.

The effect of the interface compression on the prediction is widely apparent
on $w$, $d$, and $A$. In reality, the droplet area would reduce with time due
to the stripping of small-sized droplets. One could estimate the size of these
stripped droplets as $\sim 6\ \mu$m based on $We\sim 1$. It is not possible to
resolve these small droplets with our computational grid, and in the DIM
simulations considered here, they would be represented as intermediate volume
fractions ($\epsilon<\alpha^{(2)}<1-\epsilon$). The interface compression
approach that is used here can only work with the resolved fluids, and it
wrongly interprets the regions with $\epsilon<\alpha^{(2)}<1-\epsilon$ as
numerically diffused regions, applying the compression there. As a result of
this, the results with the interface compression are insensitive to the
threshold volume fraction $\alpha^{(2)}_{T}$, but they do not match the data.
Transition to a dispersed phase Eulerian–Eulerian and Eulerian–we can consider
lagrangian modeling approaches such as [8] to model these regions with
$\epsilon<\alpha^{(2)}<1-\epsilon$, however, that is not the focus of this
work, and it would be considered in the future.

(a) $\tau=0.03$, w SR

(b) $\tau=0.03$, w/o SR

(c) $\tau=0.17$, w SR

(d) $\tau=0.17$, w/o SR

(e) $\tau=0.50$, w SR

(f) $\tau=0.50$, w/o SR

(g) $\tau=0.83$, w SR

(h) $\tau=0.83$, w/o SR

Figure 16: Pressure contours (bottom) and numerical Schlieren (top) image are
shown for the 2D shock-droplet interaction with $N_{x}=83$, with interface
compression, and with/without the stiff relaxations. The droplet is identified
using a $\alpha^{(2)}=0.5$ contour (thick red line). Figure 17: Time-evolution
of the droplet leading edge ($x_{L}$) compared against the experimental
measurements [76].  Figure 18: Time-evolution of the droplet width ($d$)
compared against the experimental measurements [76].  Figure 19: Time-
evolution of the droplet height ($d$) compared against the experimental
measurements [76].  Figure 20: Time-evolution of the droplet area ($A$)
compared against the experimental measurements [76].

###  4.7 Detonation-droplet interaction

One could encounter a detonation-droplet interaction in practical applications
such as a liquid-fueled or liquid-assisted rotating detonation engine (RDE)
[83, 84] or an energetic material packed with metal particles [85]. Some
recent efforts have focused on studying this phenomenon experimentally [53].
However, the simulations here serve as a qualitative demonstration without
detailed measurements for quantitative validation.

A 2D computational domain of size $(1500,100)$ mm is used. A liquid droplet
($\kappa=2$) of diameter $d_{0}$ mm is initialized at $(x,y)=(750,50)$ mm. A
uniform grid with $N_{x}=100$ cells across the droplet diameter is used. The
left and the right boundaries are set as supersonic outflows to let the
internal waves pass through without reflection. The top and the bottom
boundaries are set as slip walls. The entire domain is at rest initially
($u^{(1)}=u^{(2)}=0$) and the pressure are set as $p^{(1)}=p^{(2)}=10^{5}$ Pa.
The densities are $\rho^{(1)}=1.204$ kg/m3, and $\rho^{(2)}=1000$ kg/m3. The
gas-phase ($\kappa=1$) contains gaseous $\rm{H_{2}}$ and air at stoichiometry
to mimic conditions in a typical rotating detonation engine [63].

In order to initiate a 2D detonation wave, $x<10$ mm is initially filled with
detonation products at $T^{(1)}=3500$ K and $p^{(1)}=5\times 10^{6}$ Pa. A
circular disturbance of radius 1 mm is asymmetrically added in the domain to
trigger the transverse detonation waves [86].The droplet volume fraction
$\alpha^{(2)}$ is initialized using Eq. (3.23) with $\delta=0.5$. The droplet
diameter $d_{0}=40$ mm is relatively large compared to real applications.
However, it is chosen here considering the cost and the grid requirements for
resolving a droplet.

The liquid phase is modeled using an SG EOS with $\gamma^{(2)}=4.1$ and
$p^{0}=6\times 10^{8}$, and the gas phase is modeled using a TPG EOS that uses
NASA polynomial curve-fits for the species. The droplet would eventually
vaporize because of the high temperature behind the detonation. However, for
the current conditions, the ratio of the droplet vaporization time
($t_{drop}$) to the detonation passage time ($t_{deton}$) is estimated to be
substantial, $t_{drop}/t_{deton}\sim 10^{6}$, suggesting that one can simulate
the early stages of this detonation-droplet interaction without using any
heat/mass-transfer models. Similar to the previous test in Section 4.6, the
inertia primarily drives the droplet deformation at these high-speed
conditions. Therefore, the viscous and the surface tension effects are
neglected for this test.

The simulation results at various time instances are shown in Fig. 21. The
results shown here are with the interface compression scheme and the SR.
However, simulations without them have also been conducted. As a result of the
initiation, a 2D detonation wave propagates in the gas phase before reaching
the droplet. The structure of this detonation wave is confirmed by a pressure
spike (of 2.55 MPa), an increase in the temperature (to 3078 K) due to coupled
shock, and burning at the front. The detonation speed before reaching the
droplet is measured as 1968 m/s, which is extremely close to the
Chapman–-Jouguet velocity of 1964 m/s computed using a Chemical Equilibrium
Applications (CEA) code at the same conditions.

Similar to the shock-droplet interaction, as the detonation wave hits the
droplet, it results in the transmission and reflection of non-planar waves. A
circular reflected shock is observed propagating back into the detonation
products. As the wave propagates after crossing the droplet, it eventually
regains its reactive detonation structure. The droplet remains circular at
initial times. However, it deforms and elongates, as was also the case with
the shock-droplet interaction.

Figure 21: Pressure contours (first row), temperature contours (second row),
numerical Schlieren (third row), and flow-field quantities along the
centerline are shown for the 2D shock droplet interaction test. The droplet is
identified using a $\alpha^{(2)}=0.5$ contour (thick red line).

##  5 Conclusions

A seven-equation diffused interface method (DIM) framework is developed in
this work. Unlike the commonly used reduced DIM models, the seven-equation
model does not enforce a velocity/pressure equilibrium. Alleviating a strict
pressure relaxation allows this method for arbitrarily high pressure/density
ratio and arbitrary EOS of the phases, and alleviating the strict velocity
equilibrium would enable it to be used in a dispersed phase regime where the
MPE are unresolved. The seven-equation model has been used for simplistic 1D
tests [48, 28]. However, it is extended to be used in 2D/3D with arbitrary
EOS, surface tension, viscous effects, and multi-species reacting flows.
Various tests, i.e., multiphase shock tube, shock-material interaction,
oscillating droplet, drag in a viscous medium, shock-droplet interaction, and
detonation-droplet interaction, are conducted to validate different parts of
the method. Following are some novel contributions and conclusions of this
work

  * 1.

A discrete equations method (DEM) [48] is used in this work to compute the
non-conservative interphase exchange terms accurately. It is demonstrated via
1D and 2D tests that this relieves the need to use a stiff relaxation solver,
which is otherwise required to maintain pressure–velocity equilibrium between
the phases [28] and could make the simulations as much as 5$\times$ slower.
For the cases where the initial conditions are phase-specific, e.g., a shock-
material interaction problem with the shock initialized in one of the phases,
the use of the pressure/velocities equilibrium at the start results in errors
due to mismatching pre-/post-shock conditions between the phases.

  * 2.

A stiff relaxation solver is still consistent with the previous works that
used a seven-equation model, but certain modifications have been made.
Adaptive-step gradient descent is used for pressure relaxation to make the
scheme robust and useful for arbitrary EOS. The developed seven-equation
framework is robust for arbitrarily high density/pressure ratios
($\rho^{(2)}/\rho^{(1)}\sim 1-10^{3}$, $\rho^{(2)}/\rho^{(1)}\sim 1-10^{5}$),
and arbitrary EOS (SG, CPG, TPG, MG).

  * 3.

To include surface tension, a modification is proposed in the inviscid fluxes
of DEM and the corresponding HLLC Riemann solver. This proposed method can
handle a surface-tension-related pressure jump across a material interface in
an exact numerical manner. The pressure jump within a circular droplet is
captured with $<2$% error with numerically computed curvature. An oscillating
droplet where the surface potential energy converts back and forth into the
kinetic energy is simulated, and the oscillation period matches the exact
value with a $<15$% error. The subsequent peaks of kinetic energy reduce due
to numerical diffusion. However, the results improve with grid refinement.

  * 4.

Modifications for modeling 2D and 3D viscous effects within the DEM framework
have been included. Droplet evolution in a viscous medium is simulated. The
computed drag coefficient $C_{D}$ matches a steady-state drag correlation for
a low Weber number case where the droplet does not deform. Vortex shedding is
observed for a higher $Re$ case, but the droplet also deforms due to the
higher $We$.

  * 5.

An interface compression scheme initially used for the five-equation model is
modified for the seven-equation model. The interface compression was necessary
for the surface tension tests and the viscous drag computation. Otherwise, the
droplet interface gets irregular, leading to difficulties in curvature
computation and boundary and shear layer predictions surrounding the droplet
surface.

  * 6.

A shock-droplet interaction case is simulated at experimental conditions, and
the computed droplet deformation parameters are shown to match against
available data. This would create many small droplets due to the atomization
process at such high $Re$ and $We$ in reality. However, since the grid does
not resolve these, they manifest as regions of intermediate volume fractions
($\epsilon>\alpha^{(2)}>1-\epsilon$) with DIM. As a result, the interface
compression scheme was detrimental for this case since it would treat these
regions as numerically diffused regions and apply compression. We can consider
modeling these regions using a dispersed-phase approach in the future.

  * 7.

Finally, the scheme has been demonstrated to work for multi-species, reacting
flows via a detonation-droplet interaction problem. This test is a qualitative
demonstration without any validation/verification data.

This work does not consider heat and mass transfer at the droplet interface,
which is an obvious extension of this work. Infinite [42, 43] or finite-rate
[87] thermochemical relaxation solvers can be used. A key benefit of the
seven-equation model is its possible extension to dispersed-phase regions
where an equilibrium between the velocities is not maintained. It is shown in
this work that with the use of DEM, there is no need to enforce a stiff
relaxation for resolved multiphase simulations, and the next step could be to
use drag and heat and mass transfer laws for interphase exchange in the
regions where there are unresolved droplets. Such coupling of the seven-
equation DIM with the dense-to-dilute dispersed model [8] will be considered
in the future.

## Acknowledgment

This work was supported by NASA Glenn Research Center [grant number
NNX15AU91A].

## References

  * Lefebvre [1998] A. H. Lefebvre, Gas turbine combustion, CRC press, 1998.
  * Sutton [2006] G. P. Sutton, History of liquid propellant rocket engines, AIAA, 2006.
  * Tarver et al. [1996] C. M. Tarver, S. K. Chidester, A. L. Nichols,  Critical conditions for impact-and shock-induced hot spots in solid explosives,  The Journal of Physical Chemistry 100 (1996) 5794–5799.
  * Vehring [2008] R. Vehring,  Pharmaceutical particle engineering via spray drying,  Pharmaceutical research 25 (2008) 999–1022.
  * Gorokhovski and Herrmann [2008] M. Gorokhovski, M. Herrmann,  Modeling primary atomization,  Annu. Rev. Fluid Mech. 40 (2008) 343–366.
  * Saurel and Pantano [2018] R. Saurel, C. Pantano,  Diffuse-interface capturing methods for compressible two-phase flows,  Annual Review of Fluid Mechanics 50 (2018) 105–130.
  * Balachandar and Eaton [2010] S. Balachandar, J. K. Eaton,  Turbulent dispersed multiphase flow,  Annual Review of Fluid Mechanics 42 (2010) 111–133.
  * Panchal and Menon [2021] A. Panchal, S. Menon,  A hybrid Eulerian-Eulerian/Eulerian-Lagrangian method for dense-to-dilute dispersed phase flows,  Journal of Computational Physics 439 (2021) 110339.
  * Bryngelson et al. [2019] S. H. Bryngelson, K. Schmidmayer, T. Colonius,  A quantitative comparison of phase-averaged models for bubbly, cavitating flows,  International Journal of Multiphase Flow 115 (2019) 137–143.
  * Bryngelson et al. [2021] S. H. Bryngelson, R. O. Fox, T. Colonius,  Conditional moment methods for polydisperse cavitating flows,  arXiv preprint arXiv: 2112.14172 (2021).
  * Sethian and Smereka [2003] J. A. Sethian, P. Smereka,  Level set methods for fluid interfaces,  Annual Review of Fluid Mechanics 35 (2003) 341–372.
  * Sussman et al. [2007] M. Sussman, K. M. Smith, M. Y. Hussaini, M. Ohta, R. Zhi-Wei,  A sharp interface method for incompressible two-phase flows,  Journal of computational physics 221 (2007) 469–505.
  * Hirt and Nichols [1981] C. W. Hirt, B. D. Nichols,  Volume of fluid (VOF) method for the dynamics of free boundaries,  Journal of Computational Physics 39 (1981) 201–225.
  * Sussman and Puckett [2000] M. Sussman, E. G. Puckett,  A coupled level set and volume-of-fluid method for computing 3d and axisymmetric incompressible two-phase flows,  Journal of Computational Physics 162 (2000) 301–337.
  * Brackbill et al. [1992] J. Brackbill, D. B. Kothe, C. Zemach,  A continuum method for modeling surface tension,  Journal of Computational Physics 100 (1992) 335–354.
  * Fedkiw et al. [1999] R. P. Fedkiw, T. Aslam, B. Merriman, S. Osher,  A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method),  Journal of Computational Physics 152 (1999) 457–492.
  * Elghobashi [2019] S. Elghobashi,  Direct numerical simulation of turbulent flows laden with droplets or bubbles,  Annual Review of Fluid Mechanics 51 (2019) 217–244.
  * Das and Udaykumar [2020] P. Das, H. Udaykumar,  A sharp-interface method for the simulation of shock-induced vaporization of droplets,  Journal of Computational Physics 405 (2020) 109005.
  * Terashima and Tryggvason [2009] H. Terashima, G. Tryggvason,  A front-tracking/ghost-fluid method for fluid interfaces in compressible flows,  Journal of Computational Physics 228 (2009) 4012–4037.
  * Barton et al. [2011] P. T. Barton, B. Obadia, D. Drikakis,  A conservative level-set based method for compressible solid/fluid problems on fixed grids,  Journal of Computational Physics 230 (2011) 7867–7890.
  * Baer and Nunziato [1986] M. Baer, J. Nunziato,  A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials,  International journal of multiphase flow 12 (1986) 861–889.
  * Liu et al. [1994] X.-D. Liu, S. Osher, T. Chan,  Weighted essentially non-oscillatory schemes,  Journal of computational physics 115 (1994) 200–212.
  * Van Leer [1979] B. Van Leer,  Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method,  Journal of Computational Physics 32 (1979) 101–136.
  * Petitpas et al. [2009] F. Petitpas, J. Massoni, R. Saurel, E. Lapebie, L. Munier,  Diffuse interface model for high speed cavitating underwater systems,  International Journal of Multiphase Flow 35 (2009) 747–759.
  * Schmidmayer et al. [2020] K. Schmidmayer, F. Petitpas, S. Le Martelot, É. Daniel,  ECOGEN: An open-source tool for multiphase, compressible, multiphysics flows,  Computer Physics Communications 251 (2020) 107093.
  * Saurel et al. [2016] R. Saurel, P. Boivin, O. Le Métayer,  A general formulation for cavitating, boiling and evaporating flows,  Computers & Fluids 128 (2016) 53–64.
  * Rodio and Abgrall [2015] M. G. Rodio, R. Abgrall,  An innovative phase transition modeling for reproducing cavitation through a five-equation model and theoretical generalization to six and seven-equation models,  International Journal of Heat and Mass Transfer 89 (2015) 1386–1401.
  * Saurel and Abgrall [1999] R. Saurel, R. Abgrall,  A multiphase godunov method for compressible multifluid and multiphase flows,  Journal of Computational Physics 150 (1999) 425–467.
  * Kapila et al. [2001] A. Kapila, R. Menikoff, J. Bdzil, S. Son, D. S. Stewart,  Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations,  Physics of Fluids 13 (2001) 3002–3024.
  * Allaire et al. [2002] G. Allaire, S. Clerc, S. Kokh,  A five-equation model for the simulation of interfaces between compressible fluids,  Journal of Computational Physics 181 (2002) 577–616.
  * Meng and Colonius [2015] J. Meng, T. Colonius,  Numerical simulations of the early stages of high-speed droplet breakup,  Shock Waves 25 (2015) 399–414.
  * Liu et al. [2019] N. Liu, Z. Wang, M. Sun, R. Deiterding, H. Wang,  Simulation of liquid jet primary breakup in a supersonic crossflow under adaptive mesh refinement framework,  Aerospace Science and Technology 91 (2019) 456–473.
  * Murrone and Guillard [2005] A. Murrone, H. Guillard,  A five equation reduced model for compressible two phase flow problems,  Journal of Computational Physics 202 (2005) 664–698.
  * Petitpas et al. [2007] F. Petitpas, E. Franquet, R. Saurel, O. Le Metayer,  A relaxation-projection method for compressible flows. Part II: Artificial heat exchanges for multiphase shocks,  Journal of Computational Physics 225 (2007) 2214–2248.
  * Perigaud and Saurel [2005] G. Perigaud, R. Saurel,  A compressible flow model with capillary effects,  Journal of Computational Physics 209 (2005) 139–178.
  * Abgrall and Rodio [2014] R. Abgrall, M. G. Rodio,  Discrete equation method (DEM) for the simulation of viscous, compressible, two-phase flows,  Computers & Fluids 91 (2014) 164–181.
  * Saurel et al. [2008] R. Saurel, F. Petipas, R. Abgrall,  Modelling phase transition in metastable liquids: application to cavitating and flashing flows,  Journal of Fluid Mechanics 607 (2008) 313–350.
  * Bryngelson et al. [2021] S. H. Bryngelson, K. Schmidmayer, V. Coralic, J. C. Meng, K. Maeda, T. Colonius,  MFC: An open-source high-order multi-component, multi-phase, and multi-scale compressible flow solver,  Computer Physics Communications 266 (2021) 107396.
  * Schmidmayer et al. [2020] K. Schmidmayer, S. H. Bryngelson, T. Colonius,  An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics,  Journal of Computational Physics 402 (2020) 109080.
  * Wood [1930] A. B. Wood, A Textbook of Sound, G. Bell and Sons Ltd., 1930.
  * Saurel et al. [2009] R. Saurel, F. Petitpas, R. A. Berry,  Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures,  Journal of Computational Physics 228 (2009) 1678–1712.
  * Zein et al. [2010] A. Zein, M. Hantke, G. Warnecke,  Modeling phase transition for compressible two-phase flows applied to metastable liquids,  Journal of Computational Physics 229 (2010) 2964–2998.
  * Pelanti and Shyue [2019] M. Pelanti, K.-M. Shyue,  A numerical model for multiphase liquid–vapor–gas flows with interfaces and cavitation,  International journal of multiphase flow 113 (2019) 208–230.
  * Demou et al. [2022] A. D. Demou, N. Scapin, M. Pelanti, L. Brandt,  A pressure-based diffuse interface method for low-mach multiphase flows with mass transfer,  Journal of Computational Physics 448 (2022) 110730.
  * Rodriguez et al. [2021] M. Rodriguez, S. H. Bryngelson, S. Cao, T. Colonius,  Acoustically-induced bubble growth and phase change dynamics near compliant surfaces,  in: 11th International Symposium on Cavitation, 2021\.
  * Dorschner et al. [2020] B. Dorschner, L. Biasiori-Poulanges, K. Schmidmayer, H. El-Rabii, T. Colonius,  On the formation and recurrent shedding of ligaments in droplet aerobreakup,  Journal of Fluid Mechanics 904 (2020).
  * Coralic and Colonius [2014] V. Coralic, T. Colonius,  Finite-volume weno scheme for viscous compressible multicomponent flows,  Journal of Computational Physics 274 (2014) 95–121.
  * Abgrall and Saurel [2003] R. Abgrall, R. Saurel,  Discrete equations for physical and numerical compressible multiphase mixtures,  Journal of Computational Physics 186 (2003) 361–396.
  * Chang and Liou [2007] C.-H. Chang, M.-S. Liou,  A robust and accurate approach to computing compressible multiphase flow: Stratified flow model and AUSM+-up scheme,  Journal of Computational Physics 225 (2007) 840–873.
  * Tokareva and Toro [2010] S. Tokareva, E. F. Toro,  HLLC-type riemann solver for the Baer-Nunziato equations of compressible two-phase flow,  Journal of Computational Physics 229 (2010) 3573–3604.
  * Nguyen and Dumbser [2015] N. T. Nguyen, M. Dumbser,  A path-conservative finite volume scheme for compressible multi-phase flows with surface tension,  Applied Mathematics and Computation 271 (2015) 959–978.
  * Panchal [2022] A. Panchal, Modeling moderately dense to dilute multiphase flows, Ph.D. thesis, Georgia Institute of Technology, 2022\.
  * Dyson et al. [2022] D. Dyson, S. Vasu, A. Arakelyan, N. Berube, S. Briggs, J. Ramirez, E. M. Ninnemann, K. Thurmond, G. Kim, W. H. Green, et al.,  Detonation wave-induced breakup and combustion of RP-2 fuel droplets,  in: AIAA SciTech 2022 Forum, 2022, pp. AIAA–2022–1453.
  * Shukla et al. [2010] R. K. Shukla, C. Pantano, J. B. Freund,  An interface capturing method for the simulation of multi-phase compressible flows,  Journal of Computational Physics 229 (2010) 7411–7439.
  * Dodd and Ferrante [2016] M. S. Dodd, A. Ferrante,  On the interaction of Taylor length scale size droplets and isotropic turbulence,  Journal of Fluid Mechanics 806 (2016) 356–412.
  * Fechter et al. [2017] S. Fechter, C.-D. Munz, C. Rohde, C. Zeiler,  A sharp interface method for compressible liquid–vapor flow with phase transition and surface tension,  Journal of Computational Physics 336 (2017) 347–374.
  * Schmidmayer et al. [2017] K. Schmidmayer, F. Petitpas, E. Daniel, N. Favrie, S. Gavrilyuk,  A model and numerical method for compressible flows with capillary effects,  Journal of Computational Physics 334 (2017) 468–496.
  * Pelanti and Shyue [2014] M. Pelanti, K.-M. Shyue,  A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves,  Journal of Computational Physics 259 (2014) 331–357.
  * Tanguy et al. [2007] S. Tanguy, T. Ménard, A. Berlemont,  A level set method for vaporizing two-phase flows,  Journal of Computational Physics 221 (2007) 837–853.
  * Jain et al. [2020] S. S. Jain, A. Mani, P. Moin,  A conservative diffuse-interface method for compressible two-phase flows,  Journal of Computational Physics 418 (2020) 109606.
  * Chiapolino et al. [2017] A. Chiapolino, R. Saurel, B. Nkonga,  Sharpening diffuse interfaces with compressible fluids on unstructured meshes,  Journal of Computational Physics 340 (2017) 389–417.
  * Akiki et al. [2017] M. Akiki, T. P. Gallagher, S. Menon,  Mechanistic approach for simulating hot-spot formations and detonation in polymer-bonded explosives,  AIAA Journal 55 (2017) 585–598.
  * Salvadori et al. [2022] M. Salvadori, P. Tudisco, D. Ranjan, S. Menon,  Numerical investigation of mass flow rate effects on multiplicity of detonation waves within a H2/Air rotating detonation combustor,  International Journal of Hydrogen Energy 47 (2022) 4155–4170.
  * Baurle et al. [1994] R. Baurle, G. Alexopoulos, H. Hassan,  Assumed joint probability density function approach for supersonic turbulent combustion,  Journal of Propulsion and Power 10 (1994) 473–484.