# A kinetic energy–and entropy-preserving scheme for
compressible two-phase flows

Suhas S. Jain<EMAIL_ADDRESS>Parviz Moin<EMAIL_ADDRESS>Center for
Turbulence Research, Stanford University, California, USA 94305

###### Abstract

Accurate numerical modeling of compressible flows, particularly in the
turbulent regime, requires a method that is non-dissipative and stable at high
Reynolds ($Re$) numbers. For a compressible flow, it is known that discrete
conservation of kinetic energy is not a sufficient condition for numerical
stability, unlike in incompressible flows.

In this study, we adopt the recently developed conservative diffuse-interface
method (Jain, Mani & Moin, J. Comput. Phys., 2020) along with the five-
equation model for the simulation of compressible two-phase flows. This method
discretely conserves the mass of each phase, momentum, and total energy of the
system and consistently reduces to the single-phase Navier-Stokes system when
the properties of the two phases are identical. We here propose discrete
consistency conditions between the numerical fluxes, such that any set of
numerical fluxes that satisfy these conditions would not spuriously contribute
to the kinetic energy and entropy of the system. We also present a set of
numerical fluxes—which satisfies these consistency conditions—that results in
an exact conservation of kinetic energy and approximate conservation of
entropy in the absence of pressure work, viscosity, thermal diffusion effects,
and time-discretization errors. Since the model consistently reduces to the
single-phase Navier-Stokes system when the properties of the two phases are
identical, the proposed consistency conditions and numerical fluxes are also
applicable for a broader class of single-phase flows.

To this end, we present coarse-grid numerical simulations of compressible
single-phase and two-phase turbulent flows at infinite $Re$, to illustrate the
stability of the proposed method in canonical test cases, such as an isotropic
turbulence and Taylor-Green vortex flows. A higher-resolution simulation of a
droplet-laden compressible decaying isotropic turbulence is also presented,
and the effect of the presence of droplets on the flow is analyzed.

###### keywords:

compressible flows , turbulent flows , two-phase flows , phase-field method ,
split-flux forms , non-dissipative schemes

††journal: Journal of Computational Physics

##  1 Introduction

Compressible turbulent flows are encountered in a wide range of aerodynamic
flows and other high-speed flow applications. Particularly for two-phase
flows, the applications include atomization (Sallam et al., 2006), droplet
combustion (Law, 1982), bubble cavitation (Plesset and Prosperetti, 1977),
Rayleigh-Taylor instability (Youngs, 1984), and Richtmyer-Meshkov instability
(Brouillette, 2002) flows. Accurate numerical modeling of such flows requires
a method that is non-dissipative and is stable at high Reynolds numbers
($Re$). For simulations of turbulent flows, it is known that the correct
representation of evolution of kinetic energy is important for achieving
accurate numerical simulations. For incompressible simulations, Morinishi et
al. (1998) showed that the numerical instabilities can be suppressed in the
limit of zero viscosity without adding numerical dissipation if kinetic energy
is discretely conserved. For a compressible flow, however, it is known that
discrete conservation of kinetic energy is not a sufficient condition for
numerical stability, unlike in incompressible flows, and additional
constraints on the entropy of the system are needed to suppress the non-linear
instabilities (Honein, 2005).

The first class of methods that improve numerical stability consists of those
that try to minimize aliasing errors. Feiereisen (1983) proposed using skew-
symmetric splitting for convective terms to ensure that the method does not
spuriously contribute to the kinetic energy of the system. Blaisdell et al.
(1996), Kravchenko and Moin (1997), and Ducros et al. (2000) showed that the
use of skew-symmetric splitting for the convective terms reduces the aliasing
errors compared to the convective form of the discretization. Lee (1993)
showed that the use of the internal energy equation instead of the total
energy equation results in better numerical stability, and Nagarajan et al.
(2003) used a staggered scheme to show better numerical stability compared to
a collocated scheme. The skew-symmetric formulations have since been used
widely to achieve improved numerical stability for the simulations of
compressible turbulent flows (see, Subbareddy and Candler, 2009, Kok, 2009,
Morinishi, 2010, Pirozzoli, 2010, 2011). However, these methods are known to
become unstable at coarse resolutions for high $Re$, even at low Mach numbers
and in the absence of shocks, without a subgrid-scale model. For high enough
resolutions, in the limit of direct numerical simulation, all these methods
can yield stable numerical simulations.

Another class of methods consists of those that try to suppress numerical
instabilities by satisfying an entropy condition. Schemes that satisfy entropy
condition are known to exhibit stable density fluctuations. Honein and Moin
(2004) proposed a modification to the total energy equation that mimics
solving an entropy transport equation with skew-symmetric splitting, and
showed that the modified total energy equation results in higher conservation
of entropy and excellent numerical stability at high $Re$. Tadmor (1987, 2003)
proposed entropy-conservative numerical fluxes that can be combined with
additional dissipation to obtain a method that satisfies a mathematical
entropy condition for a general system of conservation laws. Following this
idea, and choosing the physical thermodynamic entropy as the mathematical
entropy function, Chandrashekar (2013) derived numerical fluxes that conserve
kinetic energy and entropy for the Euler and Navier-Stokes equations. More
recently, Coppola et al. (2019), Kuya et al. (2018) proposed schemes that have
good kinetic energy and entropy conservation properties. Hereafter, we refer
to these class of schemes as kinetic energy–and entropy-preserving (KEEP)
schemes.

It is worth mentioning that in the context of incompressible two-phase flows,
various methods have been proposed that conserve kinetic energy of the system
(Fuster, 2013, Valle et al., 2020, Mirjalili and Mani, 2021, Jain, 2022) and
therefore are stable for high-$Re$ flows. For compressible two-phase flows, on
the other hand, it was found that the use of existing KEEP schemes in the
literature, that were developed for single-phase flows, does not result in
stable numerical simulations (Jain and Moin, 2020). In addition to the
conservation of entropy and kinetic energy, the non-linear stability
properties, such as boundedness and total-variation diminishing (TVD)
properties of the volume-fraction field, are important to maintain realizable
values of the other field variables, such as density of the two-phase mixture.
We also need thermodynamic consistency at the interface to avoid spurious
numerical solutions. It is even more crucial to have all these non-linear
stability properties, if one is using a non-dissipative scheme, to obtain
stable numerical simulations. Therefore, unlike single-phase compressible
flows, entropy stability alone does not imply numerical stability for two-
phase compressible flows.

The focus of the present work is, therefore, on the development of a robust,
numerically stable scheme—a KEEP scheme—that works for both single-phase and
two-phase compressible flows. We adopt the recently developed conservative
diffuse-interface method (Jain et al., 2020) along with the five-equation
model (Allaire et al., 2002, Kapila et al., 2001) for modeling the system of
compressible two-phase flows. Some favorable properties of this method are:
(a) it discretely conserves the mass of each phase, momentum, and total energy
of the system, (b) it satisfies the boundedness and TVD properties for the
volume-fraction field, (c) it has consistency corrections (consistent
regularization) in the momentum and energy equations that maintain mechanical
equilibrium at the interface, and (d) it consistently reduces to the single-
phase Navier-Stokes system when the properties of the two phases are
identical. The proposed framework of KEEP scheme presented in this work,
however, can be adopted for any of the other interface-capturing methods, such
as a volume-of-fluid or a level-set method. We, here, propose a general
framework for the development of KEEP schemes for single-phase and two-phase
flows. We propose discrete consistency conditions for the numerical fluxes of
volume fraction, mass, momentum, kinetic energy, and internal energy, such
that an exact conservation of kinetic energy and an approximate conservation
of entropy are achieved in the absence of pressure work, viscosity, thermal
diffusion effects, and time-differencing errors. To this end, we present
numerical simulations of single-phase and two-phase turbulent flows at finite
and infinite $Re$ to illustrate the stability and accuracy of the method in
isotropic turbulence and Taylor-Green vortex flows.

The rest of this paper is organized as follows. The conservative diffuse-
interface method used in this work with a five-equation model is described in
Section 2 along with the derivation of analytical transport equations for
kinetic energy and entropy. The discrete consistency conditions between the
numerical fluxes are presented in Section 3 along with the proposed KEEP
scheme. The under-resolved simulations of single-phase and two-phase Taylor-
Green vortex and isotropic turbulence at infinite $Re$ are presented in
Section 4; and the resolved simulations of droplet-laden isotropic turbulence
at finite $Re$ are presented in Section 5, followed by a summary and
concluding remarks in Section 6.

##  2 Conservative diffuse-interface method

The conservative diffuse-interface method for compressible two-phase flows
(Jain et al., 2020) with the underlying five-equation model is given in Eqs.
(1)-(4) along with the mixture equation of state (EOS) in Eq. (5).

$\frac{\partial\phi_{1}}{\partial t}+\frac{\partial u_{j}\phi_{1}}{\partial
x_{j}}=(\phi_{1}+\zeta_{1})\frac{\partial u_{j}}{\partial
x_{j}}+\frac{\partial a_{1j}}{\partial x_{j}},$ (1) $\frac{\partial
m_{l}}{\partial t}+\frac{\partial u_{j}m_{l}}{\partial x_{j}}=\frac{\partial
R_{lj}}{\partial x_{j}},\hskip 14.22636ptl=1,2,$ (2) $\frac{\partial\rho
u_{i}}{\partial t}+\frac{\partial\rho u_{i}u_{j}}{\partial
x_{j}}+\frac{\partial p}{\partial x_{i}}=\frac{\partial u_{i}f_{j}}{\partial
x_{j}}+\frac{\partial\tau_{ij}}{\partial
x_{j}}+\sigma\kappa\frac{\partial\phi_{1}}{\partial x_{i}}+\rho g_{i},$ (3)
$\frac{\partial E}{\partial t}+\frac{\partial\left(E+p\right)u_{j}}{\partial
x_{j}}=\frac{\partial\tau_{ij}u_{i}}{\partial x_{j}}+\frac{\partial
kf_{j}}{\partial
x_{j}}+\sum_{l=1}^{2}\frac{\partial\rho_{l}h_{l}a_{lj}}{\partial
x_{j}}+\sigma\kappa u_{i}\frac{\partial\phi_{1}}{\partial x_{i}}+\rho
u_{i}g_{i},$ (4) $p=\frac{\rho
e+\left(\sum_{l=1}^{2}\frac{\phi_{l}\beta_{l}}{\alpha_{l}}\right)}{\left(\sum_{l=1}^{2}\frac{\phi_{l}}{\alpha_{l}}\right)}.$
(5)

This form of the model has a volume-fraction advection equation [Eq. (1)], a
mass balance equation for each of the phases $l$ [Eq. (2)], a momentum
equation [Eq. (3)], and a total energy equation [Eq. (4)]. If a general EOS
for phase $l$ is written as $p_{l}=\alpha_{l}\rho_{l}e_{l}+\beta_{l}$, where
$\alpha_{l}$ and $\beta_{l}$ are parameters in the EOS that are generally
determined using experimental measurements, then by invoking the isobaric
closure law for pressure in the mixture region

$p=p_{1}=p_{2},$ (6)

the generalized mixture EOS can be written as in Eq. (5). Here, a linear
stiffened-gas EOS has been assumed for simplicity. However, following the same
procedure, a generalized mixture EOS can be written for any complex cubic EOS.

We specifically use the five-equation model by Kapila et al. (2001) over the
Allaire et al. (2002) model, which includes an additional dilatational term in
the volume-fraction equation, due to its improved accuracy for flows that
involve high compression and expansion, as illustrated by Tiwari et al. (2013)
and Schmidmayer et al. (2020). Murrone and Guillard (2005) also showed that
the five-equation model by Kapila et al. (2001) admits conservative transport
equations for entropy of each phase. For weakly compressible flows, we found
that both five-equation models yield identical results.

In Eqs. (1)-(5), $\phi_{l}$ is the volume fraction of phase $l$ that satisfies
the condition $\sum_{l=1}^{2}\phi_{l}=1$; $\rho_{l}$ is the density of phase
$l$; $\rho$ is the total density, defined as
$\rho=\sum_{l=1}^{2}\rho_{l}\phi_{l}$; $u_{i}$ is the velocity; $p$ is the
pressure; $e$ is the specific mixture internal energy, which can be related to
the specific internal energy of phase $l$, $e_{l}$, as
$e=\sum_{l=1}^{2}\rho_{l}e_{l}$; $k=u_{i}u_{i}/2$ is the specific kinetic
energy; $E=\rho(e+k)$ is the total energy of the mixture per unit volume; and
the function $\zeta_{1}$ is given by

$\zeta_{1}=\frac{\rho_{2}c_{2}^{2}-\rho_{1}c_{1}^{2}}{\frac{\rho_{1}c_{1}^{2}}{\phi_{1}}+\frac{\rho_{2}c_{2}^{2}}{\phi_{2}}},$
(7)

where $c_{l}$ is the speed of sound for phase $l$. If each of the phases is
assumed to follow a stiffened-gas EOS, then the parameters in the EOS can be
written as $\alpha=\gamma-1$ and $\beta=-\gamma\pi$, where $\gamma$ is the
polytropic coefficient and $\pi$ is the reference pressure. Then, $c_{l}$ can
be defined as

$c_{l}=\sqrt{\gamma_{l}\Big{(}\frac{p+\pi_{l}}{\rho_{l}}\Big{)}}.$ (8)

In Eq. (4), $h_{l}=e_{l}+p/\rho_{l}$ represents the specific enthalpy of the
phase $l$ and can be expressed in terms of $\rho_{l}$ and $p$ using the
stiffened-gas EOS as

$h_{l}=\frac{(p+\pi_{l})\gamma_{l}}{\rho_{l}(\gamma_{l}-1)}.$ (9)

In Eqs. (1)-(4), $\sigma$ is the surface-tension coefficient,
$\kappa=-\partial n_{1j}/\partial x_{j}$ is the curvature of the interface,
$g_{i}$ is the gravitational acceleration, and
$a_{1i}=a_{i}(\phi_{1})=\Gamma\\{\epsilon\partial\phi_{1}/\partial
x_{i}-\phi_{1}(1-\phi_{1})n_{1i}\\}$ is the volumetric interface-
regularization flux for phase $1$, which satisfies the condition
$a_{i}(\phi_{1})=-a_{i}(\phi_{2})$. $n_{1i}=(\partial\phi_{1}/\partial
x_{i})/|(\partial\phi_{1}/\partial x_{i})|$ is the normal of the interface for
phase $1$; and $\Gamma$ and $\epsilon$ are the interface parameters, where
$\Gamma$ represents a regularization velocity scale and $\epsilon$ represents
an interface thickness scale. $R_{li}=\rho_{l}a_{li}$ is the interface-
regularization flux in the mass balance equation for phase $l$, and
$m_{l}=\rho_{l}\phi_{l}$ is the mass per unit total volume for phase $l$.
Summing up the mass balance equations for phases $1$ and $2$ in Eq. (2), the
total mixture mass balance equation (a modified continuity equation) can be
written as

$\frac{\partial\rho}{\partial t}+\frac{\partial\rho u_{j}}{\partial
x_{j}}=\frac{\partial f_{j}}{\partial x_{j}},$ (10)

where $f_{i}=\sum_{l=1}^{2}R_{li}=\sum_{l=1}^{2}\rho_{l}a_{li}$ is the net
interface-regularization flux for the mixture mass equation.

Invoking Stokes’s hypothesis, we write the Cauchy stress tensor as
$\tau_{ij}=2\mu D_{ij}+(\beta-2\mu/3)(\partial u_{k}/\partial
x_{k})\delta_{ij}$, where $\mu$ is the dynamic viscosity of the mixture
evaluated using the one-fluid mixture rule as
$\mu=\sum_{l=1}^{2}\phi_{l}\mu_{l}$; $D_{ij}=\\{\partial u_{i}/\partial
x_{j}+\partial u_{j}/\partial x_{i}\\}/2$ is the strain-rate tensor; and
$\beta$ is the bulk viscosity.

###  2.1 Transport equations for kinetic energy and internal energy

The transport equation for kinetic energy can be obtained by multiplying
$u_{i}$ to the momentum equation. In the absence of viscosity, surface
tension, and gravity effects, this can be written as

$u_{i}\left(\frac{\partial\rho u_{i}}{\partial t}+\frac{\partial\rho
u_{i}u_{j}}{\partial x_{j}}+\frac{\partial p}{\partial x_{i}}=\frac{\partial
u_{i}f_{j}}{\partial x_{j}}\right).$ (11)

Invoking the product rule and algebraically manipulating the equation, we
arrive at the form

$\displaystyle\frac{\partial}{\partial
t}\left(\rho\frac{u_{i}u_{i}}{2}\right)+\left(\frac{u_{i}u_{i}}{2}\right)\left\\{\frac{\partial\rho}{\partial
t}+\frac{\partial\rho u_{j}}{\partial x_{j}}-\frac{\partial f_{j}}{\partial
x_{j}}\right\\}$ (12) $\displaystyle+\frac{\partial}{\partial x_{j}}\left(\rho
u_{j}\frac{u_{i}u_{i}}{2}\right)-\frac{\partial}{\partial
x_{j}}\left(f_{j}\frac{u_{i}u_{i}}{2}\right)+\frac{\partial pu_{j}}{\partial
x_{j}}-p\frac{\partial u_{j}}{\partial x_{j}}=0.$

Now, subtracting out the continuity equation [Eq. (10)] multiplied by the term
$u_{i}u_{i}/2$, we obtain the compressible two-phase kinetic energy transport
equation:

$\displaystyle\frac{\partial}{\partial
t}\left(\rho\frac{u_{i}u_{i}}{2}\right)+\frac{\partial}{\partial
x_{j}}\left(\rho u_{j}\frac{u_{i}u_{i}}{2}\right)-\frac{\partial}{\partial
x_{j}}\left(f_{j}\frac{u_{i}u_{i}}{2}\right)+\frac{\partial pu_{j}}{\partial
x_{j}}-p\frac{\partial u_{j}}{\partial x_{j}}=0.$ (13)

The kinetic energy is not a conserved quantity in the absence of viscosity,
surface tension, and gravity force, as can be seen from the governing equation
due to the pressure term, which is not in a divergence form. The pressure term
represents the reversible exchange of energy between the kinetic and internal
energy. However, the convective and interface-regularization terms are in a
divergence form and their contribution to kinetic energy should be either only
through the boundary terms or zero, in a periodic domain. Maintaining this
property discretely is important for these terms to not contribute spuriously
to the kinetic energy.

The internal energy transport equation was derived by Jain et al. (2020), and
can be written as

$\frac{\partial\rho e}{\partial t}+\frac{\partial\rho eu_{j}}{\partial
x_{j}}+\frac{\partial pu_{j}}{\partial x_{j}}-u_{j}\frac{\partial p}{\partial
x_{j}}=\sum_{l=1}^{2}\frac{\partial\rho_{l}h_{l}a_{lj}}{\partial x_{j}},$ (14)

in the inviscid limit. Summing up the internal energy and the kinetic energy
transport equations, we can arrive at the total energy equation in Eq. (4).

###  2.2 Transport equations for individual phase and mixture entropy

Starting from the internal energy equation in Eq. (14) and rewriting this in
terms of material derivative, we obtain

$\frac{D\rho e}{Dt}+\rho h\frac{\partial u_{j}}{\partial
x_{j}}=\sum_{l=1}^{2}\frac{\partial\rho_{l}h_{l}a_{lj}}{\partial x_{j}},$

where $D/Dt=\partial/\partial t+u_{j}\partial/\partial x_{j}$ represents a
material derivative. Now, expressing the mixture internal energy in terms of
phase quantities as

$\frac{D\rho
e}{Dt}=\sum_{l=1}^{2}\frac{D(\phi_{l}\rho_{l}e_{l})}{Dt}=\sum_{l=1}^{2}\left\\{\phi_{l}\frac{D(\rho_{l}e_{l})}{Dt}+\rho_{l}e_{l}\frac{D\phi_{l}}{Dt}\right\\},$

and then utilizing Gibbs’s relation to express internal energy in terms of
specific entropy of each phase, $s_{l}$, as

$\mathrm{d}(\rho_{l}e_{l})=\rho_{l}\mathrm{d}e_{l}+e_{l}\mathrm{d}\rho_{l}=\rho_{l}T_{l}\mathrm{d}s_{l}+h_{l}\mathrm{d}\rho_{l},$

we arrive at the relation

$\sum_{l=1}^{2}\left\\{\rho_{l}\phi_{l}T_{l}\frac{Ds_{l}}{Dt}+\phi_{l}h_{l}\frac{D\rho_{l}}{Dt}+\phi_{l}h_{l}\rho_{l}\frac{\partial
u_{j}}{\partial
x_{j}}+\rho_{l}e_{l}\frac{D\phi_{l}}{Dt}\right\\}=\sum_{l=1}^{2}\frac{\partial\rho_{l}h_{l}a_{lj}}{\partial
x_{j}}.$ (15)

Substituting for $D\rho_{l}/Dt$ and $D\phi_{l}/Dt$ from Eqs. (1)-(2) in Eq.
(15), we obtain an equation that relates $Ds_{1}/Dt$ and $Ds_{2}/Dt$ as

$\rho_{1}\phi_{1}T_{1}\frac{Ds_{1}}{Dt}+\rho_{2}\phi_{2}T_{2}\frac{Ds_{2}}{Dt}=\sum_{l=1}^{2}\left\\{\frac{\partial\rho_{l}h_{l}a_{lj}}{\partial
x_{j}}-h_{l}\frac{\partial\rho_{l}a_{lj}}{\partial x_{j}}\right\\}.$ (16)

However, this provides one relation for two unknowns, $Ds_{1}/Dt$ and
$Ds_{2}/Dt$. Another relation can be obtained starting from the isobaric law
and by taking the material derivative as

$p_{1}(\rho_{1},s_{1})=p_{2}(\rho_{2},s_{2}),$ $\left.\frac{\partial
p_{1}}{\partial\rho_{1}}\right\rvert_{s_{1}}\frac{D\rho_{1}}{Dt}+\left.\frac{\partial
p_{1}}{\partial
s_{1}}\right\rvert_{\rho_{1}}\frac{Ds_{1}}{Dt}=\left.\frac{\partial
p_{2}}{\partial\rho_{2}}\right\rvert_{s_{2}}\frac{D\rho_{2}}{Dt}+\left.\frac{\partial
p_{2}}{\partial s_{2}}\right\rvert_{\rho_{2}}\frac{Ds_{2}}{Dt}.$

Now, substituting for $D\rho_{l}/Dt$, using the relations $\partial
p_{l}/\partial\rho_{l}\rvert_{s_{l}}=c_{l}$ and $\left.\partial p_{l}/\partial
s_{l}\right\rvert_{\rho_{l}}=\rho_{l}\Gamma_{l}T_{l}$—where $\Gamma_{l}$ is
the Gruneisen coefficient for phase $l$ defined as
$\Gamma_{l}=\left.(1/\rho_{l})\partial p_{l}/\partial
e_{l}\right\rvert_{\rho_{l}}$ which is equal to $\alpha_{l}$ for the EOS
considered in this work—we arrive at the equation

$\displaystyle\rho_{1}\Gamma_{1}T_{1}\frac{Ds_{1}}{Dt}-\rho_{2}\Gamma_{2}T_{2}\frac{Ds_{2}}{Dt}=\left(\frac{\partial\phi_{1}}{\partial
t}+u_{j}\frac{\partial\phi_{1}}{\partial x_{j}}-\frac{\partial
a_{1j}}{\partial
x_{j}}\right)\left(\frac{\rho_{1}c_{1}^{2}}{\phi_{1}}+\frac{\rho_{2}c_{2}^{2}}{\phi_{2}}\right)$
(17) $\displaystyle+\frac{\partial u_{j}}{\partial
x_{j}}\left(\rho_{1}c_{1}^{2}-\rho_{2}c_{2}^{2}\right)+a_{2j}\frac{\partial\rho_{2}}{\partial
x_{j}}\left(\frac{c_{2}^{2}}{\phi_{2}}\right)-a_{1j}\frac{\partial\rho_{1}}{\partial
x_{j}}\left(\frac{c_{1}^{2}}{\phi_{1}}\right).$

Now, subtracting the volume-fraction equation in Eq. (1) from this equation,
we arrive at the second relation for the two unknowns, $Ds_{1}/Dt$ and
$Ds_{2}/Dt$, as

$\rho_{1}\Gamma_{1}T_{1}\frac{Ds_{1}}{Dt}-\rho_{2}\Gamma_{2}T_{2}\frac{Ds_{2}}{Dt}=a_{2j}\frac{\partial\rho_{2}}{\partial
x_{j}}\left(\frac{c_{2}^{2}}{\phi_{2}}\right)-a_{1j}\frac{\partial\rho_{1}}{\partial
x_{j}}\left(\frac{c_{1}^{2}}{\phi_{1}}\right).$ (18)

Now, solving for $Ds_{1}/Dt$ and $Ds_{2}/Dt$ using Eqs. (16) and (18), we
arrive at the transport equation for the individual phase entropy as

$\displaystyle\rho_{l}\phi_{l}\frac{Ds_{l}}{Dt}=\frac{1}{T_{l}\left(\sum_{k=1}^{2}\frac{\Gamma_{k}}{\phi_{k}}\right)}\left[a_{l^{\prime}j}\frac{\partial\rho_{l^{\prime}}}{\partial
x_{j}}\left(\frac{c_{l^{\prime}}^{2}}{\phi_{l^{\prime}}}\right)-a_{lj}\frac{\partial\rho_{l}}{\partial
x_{j}}\left(\frac{c_{l}^{2}}{\phi_{l}}\right)\right.$ (19)
$\displaystyle+\left.\frac{\Gamma_{l^{\prime}}}{\phi_{l^{\prime}}}\sum_{k=1}^{2}\left\\{\frac{\partial\rho_{k}h_{k}a_{kj}}{\partial
x_{j}}-h_{k}\frac{\partial\rho_{k}a_{kj}}{\partial x_{j}}\right\\}\right],$

where the subscript $l^{\prime}=1-l$ denotes another phase, different from
phase $l$. This relation in Eq. (19) shows that the material derivative of
$s_{l}$ is not zero. The non-zero terms on the right-hand side is a function
of interface-regularization flux $a_{li}$, which would go to zero if the
internal structure of the interface is in a perfect equilibrium state—a
hyperbolic-tangent function. Therefore, $Ds_{l}/Dt=0$ would be satisfied when
the interface is in an equilibrium state. Now using the mass balance equation
in Eq. (2) and rewriting Eq. (19) in a conservative form, a transport equation
for the evolution of total entropy of phase $l$, $\rho_{l}\phi_{l}s_{l}$, can
be written as

$\displaystyle\frac{\partial\rho_{l}\phi_{l}s_{l}}{\partial
t}+\frac{\partial\rho_{l}u_{j}\phi_{l}s_{l}}{\partial
x_{j}}=s_{l}\frac{\partial\rho_{l}a_{lj}}{\partial x_{j}}$ (20)
$\displaystyle+\frac{1}{T_{l}\left(\sum_{k=1}^{2}\frac{\Gamma_{k}}{\phi_{k}}\right)}\left[a_{l^{\prime}j}\frac{\partial\rho_{l^{\prime}}}{\partial
x_{j}}\left(\frac{c_{l^{\prime}}^{2}}{\phi_{l^{\prime}}}\right)-a_{lj}\frac{\partial\rho_{l}}{\partial
x_{j}}\left(\frac{c_{l}^{2}}{\phi_{l}}\right)\right.$
$\displaystyle+\left.\frac{\Gamma_{l^{\prime}}}{\phi_{l^{\prime}}}\sum_{k=1}^{2}\left\\{\frac{\partial\rho_{k}h_{k}a_{kj}}{\partial
x_{j}}-h_{k}\frac{\partial\rho_{k}a_{kj}}{\partial x_{j}}\right\\}\right].$

Here again, the right-hand side of Eq. (20) is a function of interface-
regularization flux $a_{li}$, and hence, goes to zero when the interface is in
an equilibrium state. Therefore, the entropy of each phase is not exactly
conserved in a diffuse-interface method—even in the inviscid limit—due to the
interface-regularization process which is an irreversible process (see
Appendix A for the numerical quantification of entropy change associated with
the interface-regularization process).

However, an approximate conservation of entropy can be sought in the limit of
equilibrium interface state ($a_{li}\rightarrow 0$). In this limit, we recover

$\frac{\partial\rho_{l}\phi_{l}s_{l}}{\partial
t}+\frac{\partial\rho_{l}u_{j}\phi_{l}s_{l}}{\partial x_{j}}=0,\hskip
14.22636ptl=1,2,$

and

$\frac{\partial\rho s}{\partial t}+\frac{\partial\rho u_{j}s}{\partial
x_{j}}=0,$ (21)

where $\rho s=\sum_{l=1}^{2}\rho_{l}\phi_{l}s_{l}$ is the total mixture
entropy. With this, the total integrated entropy of each of the phases,
$S_{l}=\int_{\Omega}\left(\rho_{l}\phi_{l}s_{l}\right)dV$, and the mixture
entropy, $S=\int_{\Omega}\left(\rho s\right)dV$, are conserved, where $\Omega$
is the computational domain.

###  2.3 Reduced system for single-phase flows

If the material properties such as density ($\rho_{1}=\rho_{2}=\rho$),
viscosity ($\mu_{1}=\mu_{2}=\mu$), and parameters in the EOS
($\alpha_{1}=\alpha_{2}=\alpha;\beta_{1}=\beta_{2}=\beta$) are identical for
the two phases, the two-phase system in Eqs. (1)-(4) will consistently reduce
to the single-phase Navier-Stokes system, in the absence of surface tension.
This can be written as

$\frac{\partial\rho}{\partial t}+\frac{\partial\rho u_{j}}{\partial x_{j}}=0,$
(22) $\frac{\partial\rho u_{i}}{\partial t}+\frac{\partial\rho
u_{i}u_{j}}{\partial x_{j}}+\frac{\partial p}{\partial x_{i}}=0,$ (23)
$\frac{\partial E}{\partial t}+\frac{\partial\left(E+p\right)u_{j}}{\partial
x_{j}}=0,$ (24) $p=\alpha\rho e+\beta,$ (25)

in the absence of viscous and gravity effects. Now, following the steps taken
in Sections 2.1-2.2, the kinetic energy, internal energy, and entropy
transport equations for the single-phase Navier-Stokes system can be written
as

$\frac{\partial}{\partial
t}\left(\rho\frac{u_{i}u_{i}}{2}\right)+\frac{\partial}{\partial
x_{j}}\left(\rho u_{j}\frac{u_{i}u_{i}}{2}\right)+\frac{\partial
pu_{j}}{\partial x_{j}}-p\frac{\partial u_{j}}{\partial x_{j}}=0,$ (26)
$\frac{\partial\rho e}{\partial t}+\frac{\partial\rho eu_{j}}{\partial
x_{j}}+\frac{\partial pu_{j}}{\partial x_{j}}-u_{j}\frac{\partial p}{\partial
x_{j}}=0,$ (27)

and

$\frac{\partial\rho s}{\partial t}+\frac{\partial\rho u_{j}s}{\partial
x_{j}}=0,$ (28)

respectively, in the absence of viscous and gravity effects.

##  3 Numerical flux forms

In this section, we first derive the consistency conditions between the
numerical fluxes that needs to be satisfied to achieve conservation of kinetic
energy and entropy in Sections 3.2-3.3. The consistency conditions are then
summarized in Section 3.4. Finally, the proposed KEEP scheme, which consists
of a set of numerical fluxes that satisfy the proposed consistency conditions,
for single-phase and two-phase systems are presented in Section 3.5.

Starting from the mass, momentum, and energy equations in Eqs. (2)-(4), and
neglecting the viscous, surface tension, gravity, and thermal conduction
terms; rewriting the total energy as $E=\rho(e+k)$; and utilizing the mixture
mass equation in Eq. (10) instead of the individual mass equations in Eq. (2),
we can write the conservative system of equations in a general form as

$\frac{\partial U_{i}}{\partial t}+\frac{\partial F_{ij}}{\partial
x_{j}}+\frac{\partial P_{ij}}{\partial x_{j}}=\frac{\partial R_{ij}}{\partial
x_{j}},$ (29)

where $U_{i}$ is the state vector for the conserved variables and $F_{ij}$,
$P_{ij}$, and $R_{ij}$ are the flux vectors for the convective, pressure, and
interface-regularization terms, respectively, and are given by

$\begin{gathered}U_{i}=\begin{bmatrix}\rho\\\ \rho u_{i}\\\ \rho e+\rho
k\end{bmatrix},F_{ij}=\begin{bmatrix}\rho u_{j}\\\ \rho u_{i}u_{j}\\\ \rho
eu_{j}+\rho ku_{j}\end{bmatrix},\\\ P_{ij}=\begin{bmatrix}0\\\ p\delta_{ij}\\\
pu_{j}\end{bmatrix},\mathrm{and}\
R_{ij}=\begin{bmatrix}\sum_{l=1}^{2}\rho_{l}a_{lj}=f_{j}\\\ u_{i}f_{j}\\\
kf_{j}+\sum_{l=1}^{2}\rho_{l}h_{l}a_{lj}\end{bmatrix}.\end{gathered}$ (30)

Adopting the finite-volume approach, the system of equations can be
discretized in space using the conservative numerical fluxes of the form

$\left.\frac{\partial q_{j}}{\partial
x_{j}}\right|_{m}\approx\frac{\hat{q}_{j}\rvert_{(m+\frac{1}{2})}-\hat{q}_{j}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}},$ (31)

where $\hat{q}_{j}\rvert_{(m\pm 1/2)}$ represents the numerical flux for a
generic flux function $q_{j}$ evaluated on the cell face, subscript $m$
denotes the grid-cell index, and $\Delta x_{j}$ represents the grid-cell size.
Now, the semi-discrete representation for the system of equations can be
written as

$\left.\frac{\partial U_{i}}{\partial
t}\right|_{m}+\frac{\hat{F}_{ij}\rvert_{(m+\frac{1}{2})}-\hat{F}_{ij}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}+\frac{\hat{P}_{ij}\rvert_{(m+\frac{1}{2})}-\hat{P}_{ij}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}=\frac{\hat{R}_{ij}\rvert_{(m+\frac{1}{2})}-\hat{R}_{ij}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}},$ (32)

where $\hat{F}_{ij}$, $\hat{P}_{ij}$, and $\hat{R}_{ij}$ are the numerical
flux vectors for the convective, pressure, and interface-regularization terms,
respectively. The components of the numerical flux vectors can be written as

$\begin{gathered}\hat{F}_{ij}\rvert_{(m\pm\frac{1}{2})}=\begin{bmatrix}\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("\rho
u_{j}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{\rho
u_{j}}{\tmpbox}\rvert_{(m\pm\frac{1}{2})}=\hat{C}_{j}\rvert_{(m\pm\frac{1}{2})}\\\
\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("\rho
u_{i}u_{j}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{\rho
u_{i}u_{j}}{\tmpbox}\rvert_{(m\pm\frac{1}{2})}=\hat{M}_{ij}\rvert_{(m\pm\frac{1}{2})}\\\
\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("\rho
eu_{j}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{\rho
eu_{j}}{\tmpbox}\rvert_{(m\pm\frac{1}{2})}+\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("\rho
ku_{j}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{\rho
ku_{j}}{\tmpbox}\rvert_{(m\pm\frac{1}{2})}=\hat{I}_{j}\rvert_{(m\pm\frac{1}{2})}+\hat{K}_{j}\rvert_{(m\pm\frac{1}{2})}\end{bmatrix},\\\
\hat{P}_{ij}\rvert_{(m\pm\frac{1}{2})}=\begin{bmatrix}0\\\
\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("p\delta_{ij}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{p\delta_{ij}}{\tmpbox}\rvert_{(m\pm\frac{1}{2})}=\hat{P}_{ij}\rvert_{(m\pm\frac{1}{2})}\\\
\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("pu_{j}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{pu_{j}}{\tmpbox}\rvert_{(m\pm\frac{1}{2})}=\hat{W}_{j}\rvert_{(m\pm\frac{1}{2})}\end{bmatrix},\\\
\mathrm{and}\
\hat{R}_{ij}\rvert_{(m\pm\frac{1}{2})}=\begin{bmatrix}\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("f_{j}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{f_{j}}{\tmpbox}\rvert_{(m\pm\frac{1}{2})}=\hat{F}_{j}\rvert_{(m\pm\frac{1}{2})}\\\
\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("u_{i}f_{j}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{u_{i}f_{j}}{\tmpbox}\rvert_{(m\pm\frac{1}{2})}=\hat{R}_{ij}\rvert_{(m\pm\frac{1}{2})}\\\
\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("kf_{j}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{kf_{j}}{\tmpbox}\rvert_{(m\pm\frac{1}{2})}+\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("\sum_{l=1}^{2}\rho_{l}h_{l}a_{lj}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{\sum_{l=1}^{2}\rho_{l}h_{l}a_{lj}}{\tmpbox}\rvert_{(m\pm\frac{1}{2})}=\hat{T}_{j}\rvert_{(m\pm\frac{1}{2})}+\hat{H}_{j}\rvert_{(m\pm\frac{1}{2})}\end{bmatrix}.\end{gathered}$
(33)

Using any form of numerical flux for the components in Eq. (33) leads to an
implementation that discretely conserves the variables in the state vector
$U_{i}$ in Eq. (32) due to the telescoping property. However, to achieve
secondary conservation of the kinetic energy and entropy, the numerical fluxes
need to be approximated in a consistent manner, such that the conservation can
be achieved implicitly.

In addition to these fluxes, the convective numerical fluxes in the volume
fraction equation [Eq. (1)] and individual phase mass equation [Eq. (2)] can
be written as

$\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("\phi
u_{j}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{\phi
u_{j}}{\tmpbox}\rvert_{(m\pm\frac{1}{2})}=\hat{\Phi}_{j}\rvert_{(m\pm\frac{1}{2})},$
(34)

and

$\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("m_{l}u_{j}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{m_{l}u_{j}}{\tmpbox}\rvert_{(m\pm\frac{1}{2})}=\hat{Q}_{lj}\rvert_{(m\pm\frac{1}{2})},$
(35)

respectively; and the interface-regularization fluxes in the volume fraction
equation and individual phase mass equation can be written as

$\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("a_{lj}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{a_{lj}}{\tmpbox}\rvert_{(m\pm\frac{1}{2})}=\hat{A}_{lj}\rvert_{(m\pm\frac{1}{2})},$
(36)

and

$\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("\rho_{l}a_{lj}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{\rho_{l}a_{lj}}{\tmpbox}\rvert_{(m\pm\frac{1}{2})}=\hat{B}_{lj}\rvert_{(m\pm\frac{1}{2})},$
(37)

A consistent numerical approximation for these fluxes in Eqs. (34)-(37) is
also crucial in achieving secondary conservation of the kinetic energy and
entropy.

###  3.1 Skew-symmetric-like splitting

The derivatives of the non-linear products that typically arise in the
convective terms can be discretized in different ways. If the generic non-
linear flux in Eq. (31) can be expressed as $q_{1}=abc$, where $a$, $b$, and
$c$ are functions of $x$, then according to Eq. (31), the three possible
numerical flux forms for $\hat{q}_{1}\rvert_{(m\pm 1/2)}$ are

$\begin{gathered}\hat{q}_{1}\rvert_{(m\pm\frac{1}{2})}=\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$abc$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})},\\\
\hat{q}_{1}\rvert_{(m\pm\frac{1}{2})}=\hbox{\vbox{\hrule height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$ab$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$c$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})},\\\
\hat{q}_{1}\rvert_{(m\pm\frac{1}{2})}=\hbox{\vbox{\hrule height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$a$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$b$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$c$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})},\\\
\end{gathered}$ (38)

which are called the divergence form, quadratic-split form, and cubic-split
form, respectively, and the overbar denotes an arithmetic average of quantity
evaluated at $m$ and $m\pm 1$, e.g., $\hbox{\vbox{\hrule height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$ab$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}$ is

$\hbox{\vbox{\hrule height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$ab$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}=\frac{ab\rvert_{m}+ab\rvert_{(m\pm
1)}}{2}.$ (39)

However, the interface-regularization flux terms contain non-linear product of
four or five quantities, and therefore, new numerical split forms need to be
defined. If $q_{1}=abcd$, then a quartic-split form of the numerical flux can
be defined as

$\hat{q}_{1}\rvert_{(m\pm\frac{1}{2})}=\hbox{\vbox{\hrule height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$a$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$b$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$c$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$d$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})},$ (40)

Similarly, if $q_{1}=abcde$, then a quintic-split form of the numerical flux
can be defined as

$\hat{q}_{1}\rvert_{(m\pm\frac{1}{2})}=\hbox{\vbox{\hrule height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$a$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$b$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$c$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$d$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$e$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})},$ (41)

These skew-symmetric-like split schemes are used here, because these schemes
result in reduced aliasing errors (Blaisdell et al., 1996, Chow and Moin,
2003, Kennedy and Gruber, 2008) and are known to improve the conservation
properties of the quadratic quantities (Kravchenko and Moin, 1997).

The numerical fluxes, introduced here, can also be equivalently rewritten in a
pointwise manner, e.g., the divergence form, quadratic-split form, and cubic-
split forms of the numerical flux can be written as

$\begin{gathered}\frac{\partial q_{1}}{\partial
x_{1}}\approx\left.\frac{\delta abc}{\delta x}\right|_{m},\\\ \frac{\partial
q_{1}}{\partial x_{1}}\approx\frac{1}{2}\left(\left.\frac{\delta abc}{\delta
x}\right|_{m}+\left.c\right|_{m}\left.\frac{\delta ab}{\delta
x}\right|_{m}+\left.ab\right|_{m}\left.\frac{\delta c}{\delta
x}\right|_{m}\right),\\\ \frac{\partial q_{1}}{\partial
x_{1}}\approx\frac{1}{4}\left(\left.\frac{\delta abc}{\delta
x}\right|_{m}+\left.c\right|_{m}\left.\frac{\delta ab}{\delta
x}\right|_{m}+\left.a\right|_{m}\left.\frac{\delta bc}{\delta
x}\right|_{m}+\left.b\right|_{m}\left.\frac{\delta ac}{\delta
x}\right|_{m}\right.\\\ \left.+\left.ab\right|_{m}\left.\frac{\delta c}{\delta
x}\right|_{m}+\left.bc\right|_{m}\left.\frac{\delta a}{\delta
x}\right|_{m}+\left.ac\right|_{m}\left.\frac{\delta b}{\delta
x}\right|_{m}\right),\end{gathered}$ (42)

respectively, where the $\delta/\delta x$ operator is the standard finite-
difference derivative.

###  3.2 Consistency conditions for volume fraction, mass, momentum, and
kinetic-energy fluxes

The kinetic energy equations in Eqs. (13), (26) are not solved directly, and
hence, we seek an implicit conservation of kinetic energy here. One approach
to achieve this is by following the same steps taken in Section 2.1
discretely, and by seeking those forms of numerical fluxes in the mass,
momentum, and energy equations that result in the implied discrete kinetic
energy equation that has discrete divergence operators in all the terms. In
this case, the contribution of the terms in the discrete kinetic energy
equation will be only through the boundary terms or zero, in a periodic
domain, due to the telescoping property of the divergence operator. Note that
we are seeking kinetic energy conservation in the limit of zero viscosity,
surface tension, gravity force, and reversible pressure work.

Theorem 3.1 proposes consistency conditions between $\hat{C}_{j}$ and
$\hat{M}_{ij}$, and $\hat{F}_{j}$ and $\hat{R}_{ij}$, such that the discrete
convective and interface-regularization terms do not contribute spuriously to
the total kinetic energy.

######  Theorem 3.1.

The total kinetic energy of the system can be discretely conserved in the
absence of reversible pressure work, viscosity, surface tension, gravity, and
time-differencing errors if the following consistency conditions are
satisfied: those between the convective fluxes of mass $\hat{C}$ and momentum
$\hat{M}$,

$\hat{M}_{ij}\rvert_{(m\pm\frac{1}{2})}=\hat{C}_{j}\rvert_{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u_{i}$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}$
(43)

and between the interface-regularization fluxes of mass $\hat{F}$ and momentum
$\hat{R}$,

$\hat{R}_{ij}\rvert_{(m\pm\frac{1}{2})}=\hat{F}_{j}\rvert_{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u_{i}$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}$
(44)

where $u$ is the velocity, $m$ is the grid index, and $i$ and $j$ are the
indices in the Einstein notation.

###### Proof.

Ignoring the time-stepping errors, we can obtain the semi-discrete kinetic
energy equation by multiplying the semi-discrete momentum equation by $u_{i}$

$u_{i}\frac{\partial\rho u_{i}}{\partial
t}+u_{i}\rvert_{m}\left(\frac{\hat{M}_{ij}\rvert_{(m+\frac{1}{2})}-\hat{M}_{ij}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)-u_{i}\rvert_{m}\left(\frac{\hat{R}_{ij}\rvert_{(m+\frac{1}{2})}-\hat{R}_{ij}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)+u_{i}(\mathrm{O.T.})=0,$ (45)

where $\mathrm{O.T.}$ represents other terms such as pressure, viscosity,
surface tension, and gravity force. Here, Eq. (45) is a discrete analogue of
Eq. (11). Now, suppose steps taken in going from Eq. (11) to Eq. (12) are
valid discretely, then, discretely following the steps in going from Eq. (12)
to Eq. (13), i.e., subtracting out the continuity equation contribution
multiplied by $u_{i}u_{i}/2$ term from this equation, we will arrive at the
discrete analogue of the kinetic energy equation in Eq. 13 with the required
form of the discrete divergence operators for the convective and interface-
regularization terms. This is sufficient to show

$\sum_{m}\frac{\partial\rho k}{\partial t}=0.$

Therefore, the required condition to achieve discrete kinetic energy
conservation is that the steps taken in going from Eq. (11) to Eq. (12) be
valid discretely. The discrete analogue of Eq. (12) can be written as

$\displaystyle\frac{\partial\rho k}{\partial
t}+\frac{1}{2}u_{i}u_{i}\rvert_{m}\frac{\partial\rho}{\partial
t}+\frac{1}{2}u_{i}u_{i}\rvert_{m}\left(\frac{\hat{C}_{j}\rvert_{(m+\frac{1}{2})}-\hat{C}_{j}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)-\frac{1}{2}u_{i}u_{i}\rvert_{m}\left(\frac{\hat{F}_{j}\rvert_{(m+\frac{1}{2})}-\hat{F}_{j}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)$ (46)
$\displaystyle+\left(\frac{\hat{K}_{j}\rvert_{(m+\frac{1}{2})}-\hat{K}_{j}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)-\left(\frac{\hat{T}_{j}\rvert_{(m+\frac{1}{2})}-\hat{T}_{j}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)+u_{i}(\mathrm{O.T.})=0.$

Therefore, Eq. (46) should be equivalent to Eq. (45) to achieve discrete
kinetic energy conservation. Seeking global conservation of kinetic energy, we
equate Eq. (45) and Eq. (46) under the summation over all grid cells as

$\displaystyle\sum_{m}\left\\{u_{i}\rvert_{m}\left(\frac{\hat{M}_{ij}\rvert_{(m+\frac{1}{2})}-\hat{M}_{ij}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)-u_{i}\rvert_{m}\left(\frac{\hat{R}_{ij}\rvert_{(m+\frac{1}{2})}-\hat{R}_{ij}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)\right\\}$ (47)
$\displaystyle=\sum_{m}\Bigg{\\{}\frac{1}{2}u_{i}u_{i}\rvert_{m}\left(\frac{\hat{C}_{j}\rvert_{(m+\frac{1}{2})}-\hat{C}_{j}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)-\frac{1}{2}u_{i}u_{i}\rvert_{m}\left(\frac{\hat{F}_{j}\rvert_{(m+\frac{1}{2})}-\hat{F}_{j}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)$
$\displaystyle+\left(\frac{\hat{K}_{j}\rvert_{(m+\frac{1}{2})}-\hat{K}_{j}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)-\left(\frac{\hat{T}_{j}\rvert_{(m+\frac{1}{2})}-\hat{T}_{j}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)\Bigg{\\}}.$

$\hat{K}_{j}\rvert_{(m\pm\frac{1}{2})}$ and
$\hat{T}_{j}\rvert_{(m\pm\frac{1}{2})}$ terms form telescoping series and
therefore their contribution cancels out. Expanding the summation and
separating the contributions from the convective and interface-regularization
terms, we arrive at the conditions

$\hat{M}_{j}\rvert_{(m\pm\frac{1}{2})}=\hat{C}_{j}\rvert_{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u_{i}$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})},$
(48)

and

$\hat{R}_{j}\rvert_{(m\pm\frac{1}{2})}=\hat{F}_{j}\rvert_{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u_{i}$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})},$
(49)

which concludes the proof. ∎

The consistency condition in Eq. (43) is similar to the kinetic energy
preserving (KEP) flux by Jameson (2008), which is a sufficient condition for
numerical stability of incompressible flows because of the absence of exchange
of energy between kinetic energy and internal energy, unlike for compressible
flows. For compressible flows, an additional constraint on the entropy of the
system is needed to achieve stable numerical simulations.

There are many approaches to achieve (exact or approximate) discrete
conservation of entropy. One option would be to directly solve a transport
equation for entropy instead of an energy equation, which was explored in
Honein (2005). They also proposed a modification to the total energy equation,
in Honein and Moin (2004), which would result in higher conservation of
entropy ($\rho s$ and $\rho s^{2}$) implicitly. Another option would be to
solve an internal energy equation or a total energy equation with consistent
numerical discretization that results in conservation of entropy implicitly.
This approach has been previously used by Chandrashekar (2013), Kuya et al.
(2018), and Coppola et al. (2019). This approach that involves solving an
energy equation with implicit conservation of entropy is a more preferred
option since it also results in the primary conservation of energy, which is
crucial for compressible flows. Therefore, we choose to adopt the second
option in this work.

Although, the consistency conditions in Theorem 3.1 result in discretely
conserving the total kinetic energy, in the absence of reversible pressure
work, this is not a sufficient condition for numerical stability of
compressible flows. This is because an incorrect local exchange of energy
between kinetic energy and internal energy could spuriously contribute to
entropy of the system. Hence, the first step toward achieving conservation of
entropy implicitly is to propose additional constraints on the transport of
kinetic energy, such that there is no spurious local exchange of energy
between kinetic energy and internal energy (a “local” kinetic energy
preservation). To this end, Theorem 3.2 proposes consistency conditions
between $\hat{C}_{j}$ and $\hat{K}_{j}$, and $\hat{F}_{j}$ and $\hat{T}_{j}$,
such that the discrete convective and interface-regularization terms do not
contribute spuriously to the local transport of kinetic energy.

######  Theorem 3.2.

The discrete kinetic energy of the system can be locally conserved in the
absence of reversible pressure work, viscosity, surface tension, gravity, and
time-differencing errors if the following consistency conditions are
satisfied: those between the convective fluxes of mass $\hat{C}$ and kinetic
energy $\hat{K}$,

$\hat{K}_{j}\rvert_{(m\pm\frac{1}{2})}=\hat{C}_{j}\rvert_{(m\pm\frac{1}{2})}\frac{u_{i}\rvert_{(m\pm
1)}u_{i}\rvert_{m}}{2},$ (50)

and between the interface-regularization fluxes of mass $\hat{F}$ and kinetic
energy $\hat{T}$,

$\hat{T}_{j}\rvert_{(m\pm\frac{1}{2})}=\hat{F}_{j}\rvert_{(m\pm\frac{1}{2})}\frac{u_{i}\rvert_{(m\pm
1)}u_{i}\rvert_{m}}{2},$ (51)

where $u$ is the velocity, $m$ is the grid index, and $i$ and $j$ are the
indices in the Einstein notation.

###### Proof.

Following the proof of Theorem 3.1, and now instead seeking local conservation
of kinetic energy, we equate Eq. (45) and Eq. (46) in a pointwise manner as

$\displaystyle
u_{i}\rvert_{m}\left(\frac{\hat{M}_{ij}\rvert_{(m+\frac{1}{2})}-\hat{M}_{ij}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)-u_{i}\rvert_{m}\left(\frac{\hat{R}_{ij}\rvert_{(m+\frac{1}{2})}-\hat{R}_{ij}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)$ (52)
$\displaystyle=\frac{1}{2}u_{i}u_{i}\rvert_{m}\left(\frac{\hat{C}_{j}\rvert_{(m+\frac{1}{2})}-\hat{C}_{j}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)-\frac{1}{2}u_{i}u_{i}\rvert_{m}\left(\frac{\hat{F}_{j}\rvert_{(m+\frac{1}{2})}-\hat{F}_{j}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)$
$\displaystyle+\left(\frac{\hat{K}_{j}\rvert_{(m+\frac{1}{2})}-\hat{K}_{j}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)-\left(\frac{\hat{T}_{j}\rvert_{(m+\frac{1}{2})}-\hat{T}_{j}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right).$

Making use of the conditions in Eqs. (48) and (49), and separating the
contributions from the convective and interface-regularization terms, we
arrive at the conditions

$\hat{K}_{j}\rvert_{(m\pm\frac{1}{2})}=\hat{C}_{j}\rvert_{(m\pm\frac{1}{2})}\frac{u_{i}\rvert_{(m\pm
1)}u_{i}\rvert_{m}}{2},$ (53)

and

$\hat{T}_{j}\rvert_{(m\pm\frac{1}{2})}=\hat{F}_{j}\rvert_{(m\pm\frac{1}{2})}\frac{u_{i}\rvert_{(m\pm
1)}u_{i}\rvert_{m}}{2},$ (54)

which concludes the proof. ∎

When individual mass balance equations in Eq. (2) are solved for each of the
phases, they are supposed to be solved consistently with other equations.
Hence, additional consistency conditions can be derived for these. Note that
the mass balance equations for phases $1$ and $2$ were summed to obtain the
total mixture mass balance equation in Eq. (10), and therefore, the same
criterion should hold for numerical fluxes, e.g., as

$\sum_{l=1}^{2}\left.\frac{\delta u_{j}m_{l}}{\delta
x_{j}}\right|_{m}=\left.\frac{\delta u_{j}\rho}{\delta x_{j}}\right|_{m}.$
(55)

This results in the consistency condition between $\hat{C}_{j}$ and
$\hat{Q}_{lj}$ as

$\sum_{l=1}^{2}\hat{Q}_{lj}\rvert_{(m\pm 1/2)}=\hat{C}_{j}\rvert_{(m\pm
1/2)}.$ (56)

Similarly, the consistency condition between $\hat{F}_{j}$ and $\hat{B}_{lj}$
can be written as

$\sum_{l=1}^{2}\hat{B}_{lj}\rvert_{(m\pm 1/2)}=\hat{F}_{j}\rvert_{(m\pm
1/2)}.$ (57)

In the limit of incompressibility, mass balance equations in Eq. (2) reduce to
the volume fraction equation in Eq. (1) (Jain et al., 2020). This can be used
to obtain a consistency condition between $\hat{C}_{j}$ and $\hat{\Phi}_{j}$.
The consistency condition states that the split-numerical flux form used in
$\hat{C}_{j}\rvert_{(m\pm\frac{1}{2})}$ should be the same as the one used in
$\hat{\Phi}_{j}\rvert_{(m\pm\frac{1}{2})}$.

###  3.3 Consistency conditions for internal energy, pressure, and enthalpy
fluxes

In addition to satisfying the local kinetic energy preservation in Theorem
3.2, in order to achieve implicit conservation of entropy without having to
solve for Eqs. (20) or (21), all the steps taken in Section 2.2 in arriving at
Eq. (21) starting from the internal energy in Eq, (14) need to hold
discretely. The mass balance equation and the volume-fraction advection
equations used in Eqs. (15) and Eq. (17) are directly solved and hence are
satisfied discretely. The only remaining constraint is that the method needs
to satisfy the interface-equilibrium condition (IEC).

According to Abgrall (1996), if $u^{k}\rvert_{m}=u_{0}$ and
$p^{k}\rvert_{m}=p_{0}$, where $k$ is the time-step index, any model or a
numerical scheme that satisfies $u\rvert_{m}^{k+1}=u_{0}$ and
$p\rvert_{m}^{k+1}=p_{0}$, $\forall m$, is said to satisfy the IEC. It is
known that the IEC needs to be satisfied to avoid spurious solutions and to
achieve stable numerical solutions. Jain et al. (2020) derived the analytical
formulation of the interface-regularization terms in the conservative diffuse-
interface method such that the resulting internal energy equation in Eq. (14)
satisfies the IEC, and showed that this results in approximate conservation of
entropy (analytically). We here, seek fluxes that verify discretely the IEC
which results in improved discrete entropy preservation. To do so, we
formulate consistency conditions between $\hat{I}_{j}$ and $\hat{C}_{j}$, and
$\hat{H}_{j}$ and $\hat{F}_{j}$ in Lemma 3.1 and the proof is presented in
Appendix B.

######  Lemma 3.1.

The consistency conditions between $\hat{I}_{j}$ and $\hat{C}_{j}$, and
$\hat{H}_{j}$ and $\hat{F}_{j}$, required to satisfy IEC, can be stated as:
the split-numerical flux form used in $\hat{I}_{j}\rvert_{(m\pm\frac{1}{2})}$
should be the same as the one in $\hat{C}_{j}\rvert_{(m\pm\frac{1}{2})}$, and
the split-numerical flux form used in $\hat{H}_{j}\rvert_{(m\pm\frac{1}{2})}$
should be the same as the one in $\hat{F}_{j}\rvert_{(m\pm\frac{1}{2})}$.

Finally, for the pressure flux terms, a discrete consistency condition can be
derived as follows. The sum of pressure work in the transport equation for
kinetic energy in Eq. (13) and the pressure dilatation term in the transport
equation for the internal energy in Eq. (14) gives the conservative pressure
term in the total energy equation in Eq. (4), an idea that was also used by
Kok (2009), Kuya et al. (2018). This condition should be satisfied discretely,
and can be written as

$\left.u_{j}\frac{\delta p}{\delta x_{j}}\right|_{m}+\left.p\frac{\delta
u_{j}}{\delta x_{j}}\right|_{m}=\left.\frac{\delta pu_{j}}{\delta
x_{j}}\right|_{m}$

Writing this in terms of fluxes, a consistency condition between $\hat{W}_{j}$
and $\hat{P}_{ij}$ is given by

$\hat{P}_{ij}\rvert_{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u_{i}$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}+p\rvert_{m}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u_{j}$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}=\hat{W}_{j}\rvert_{(m\pm\frac{1}{2})}.$

###  3.4 Summary of the consistency conditions

The consistency conditions used in this work can be summarized as follows:

  1. 1.

The kinetic energy in the system is discretely conserved, globally and
locally, in the absence of reversible pressure work, dissipative mechanisms,
and time-differencing errors, if the following relations are satisfied

$\hat{M}_{ij}\rvert_{(m\pm\frac{1}{2})}=\hat{C}_{j}\rvert_{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u_{i}$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})},$
$\hat{K}_{j}\rvert_{(m\pm\frac{1}{2})}=\hat{C}_{j}\rvert_{(m\pm\frac{1}{2})}\frac{u_{i}\rvert_{(m\pm
1)}u_{i}\rvert_{m}}{2},$
$\hat{R}_{ij}\rvert_{(m\pm\frac{1}{2})}=\hat{F}_{j}\rvert_{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u_{i}$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})},$
$\hat{T}_{j}\rvert_{(m\pm\frac{1}{2})}=\hat{F}_{j}\rvert_{(m\pm\frac{1}{2})}\frac{u_{i}\rvert_{(m\pm
1)}u_{i}\rvert_{m}}{2}.$

  2. 2.

The mixture mass flux can be obtained by summing the mass fluxes of the
individual phases, which results in the relations

$\sum_{l=1}^{2}\hat{Q}_{lj}\rvert_{(m\pm 1/2)}=\hat{C}_{j}\rvert_{(m\pm
1/2)},$ $\sum_{l=1}^{2}\hat{B}_{lj}\rvert_{(m\pm
1/2)}=\hat{F}_{j}\rvert_{(m\pm 1/2)}.$

  3. 3.

The split-numerical flux form used in $\hat{C}_{j}\rvert_{(m\pm\frac{1}{2})}$
should be the same as the ones used in
$\hat{\Phi}_{j}\rvert_{(m\pm\frac{1}{2})}$ and
$\hat{I}_{j}\rvert_{(m\pm\frac{1}{2})}$ to satisfy IEC.

  4. 4.

The split-numerical flux form used in $\hat{F}_{j}\rvert_{(m\pm\frac{1}{2})}$
should be the same as the one used in $\hat{H}_{j}\rvert_{(m\pm\frac{1}{2})}$
to satisfy IEC.

  5. 5.

The sum of pressure work in the kinetic energy equation and the pressure
dilatation term in the internal energy equation should be equal to the
conservative pressure-diffusion term in the total energy equation. This
results in the relation

$\hat{P}_{ij}\rvert_{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u_{i}$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}+p\rvert_{m}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u_{j}$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}=\hat{W}_{j}\rvert_{(m\pm\frac{1}{2})}.$

The consistency conditions $1$–$5$ completely define the numerical scheme,
provided $\hat{C}_{j}$ and $\hat{F}_{j}$ are chosen. Note that for the single-
phase Navier-Stokes system, the consistency conditions in $1$, $3$, and $5$
are sufficient to completely define the numerical scheme, provided
$\hat{C}_{j}$ is chosen.

###  3.5 The proposed KEEP scheme

The convective mass flux is a product of $\rho$ and $u_{j}$. Hence, using the
quadratic-split form defined in Eq. (38) for $\hat{C}_{j}$, we can derive the
numerical fluxes that result in a consistent set of numerical discretization
for the single-phase Navier-Stokes system, using the consistency conditions in
Section 3.4, as

$\displaystyle\hat{C}_{j}\rvert_{(m\pm\frac{1}{2})}=\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\rho$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u$\kern-1.00006pt}}}_{j}^{(m\pm\frac{1}{2})},$
(58) $\displaystyle\hat{M}_{ij}\rvert_{(m\pm\frac{1}{2})}=\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\rho$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u$\kern-1.00006pt}}}_{j}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u$\kern-1.00006pt}}}_{i}^{(m\pm\frac{1}{2})},$
$\displaystyle\hat{K}_{j}\rvert_{(m\pm\frac{1}{2})}=\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\rho$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u$\kern-1.00006pt}}}_{j}^{(m\pm\frac{1}{2})}\frac{u_{i}\rvert_{(m\pm
1)}u_{i}\rvert_{m}}{2},$
$\displaystyle\hat{I}_{j}\rvert_{(m\pm\frac{1}{2})}=(\hbox{\vbox{\hrule
height=0.5pt\kern 2.15277pt\hbox{\kern-1.00006pt$\rho
e$\kern-1.00006pt}}})^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u$\kern-1.00006pt}}}_{j}^{(m\pm\frac{1}{2})},$
$\displaystyle\hat{P}_{ij}\rvert_{(m\pm\frac{1}{2})}=\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$p$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\delta_{ij},$
$\displaystyle\hat{W}_{j}\rvert_{(m\pm\frac{1}{2})}=\frac{u_{j}\rvert_{(m\pm
1)}p\rvert_{m}+u_{j}\rvert_{m}p\rvert_{(m\pm 1)}}{2}.$

The mass-regularization flux in Eq. (10) consists of a diffusive flux term and
a non-linear sharpening flux term, that can be written as

$\hat{F}_{j}\rvert_{(m\pm\frac{1}{2})}=\sum_{l=1}^{2}\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("\rho_{l}a_{lj}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{\rho_{l}a_{lj}}{\tmpbox}\rvert_{(m\pm\frac{1}{2})}=\Gamma\sum_{l=1}^{2}\left.\left\\{\epsilon\underbrace{\left(\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("\rho_{l}\frac{\partial\phi_{l}}{\partial
x_{j}}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{\rho_{l}\frac{\partial\phi_{l}}{\partial
x_{j}}}{\tmpbox}\right)}_{\begin{subarray}{c}\textrm{diffusive}\\\
\textrm{flux}\end{subarray}}-\underbrace{\left(\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("\rho_{l}\phi_{l}(1-\phi_{l})n_{lj}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{\rho_{l}\phi_{l}(1-\phi_{l})n_{lj}}{\tmpbox}\right)}_{\begin{subarray}{c}\textrm{sharpening}\\\
\textrm{flux}\end{subarray}}\right\\}\right|_{(m\pm\frac{1}{2})}.$

The sharpening flux is a product of four variables $\rho_{l}$, $\phi_{l}$,
$1-\phi_{l}$, and $n_{li}$. Now, approximating this term using the quartic-
split form defined in Eq. (40), we can derive the remaining numerical fluxes,
using the consistency conditions in Section 3.4, as

$\displaystyle\hat{\Phi}_{j}\rvert_{(m\pm\frac{1}{2})}=\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\phi$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u$\kern-1.00006pt}}}_{j}^{(m\pm\frac{1}{2})},$
(59) $\displaystyle\hat{F}_{j}\rvert_{(m\pm
1/2)}=\sum_{l=1}^{2}\left(\hbox{\vbox{\hrule height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\rho_{l}$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("a")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{a}{\tmpbox}_{lj}\rvert_{(m\pm\frac{1}{2})}\right),$
$\displaystyle\hat{R}_{ij}\rvert_{(m\pm
1/2)}=\left\\{\sum_{l=1}^{2}\left(\hbox{\vbox{\hrule height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\rho_{l}$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("a")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{a}{\tmpbox}_{lj}\rvert_{(m\pm\frac{1}{2})}\right)\right\\}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u$\kern-1.00006pt}}}_{i}^{(m\pm\frac{1}{2})},$
$\displaystyle\hat{T}_{j}\rvert_{(m\pm\frac{1}{2})}=\left\\{\sum_{l=1}^{2}\left(\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\rho_{l}$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("a")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{a}{\tmpbox}_{lj}\rvert_{(m\pm\frac{1}{2})}\right)\right\\}\frac{u_{i}\rvert_{(m\pm
1)}u_{i}\rvert_{m}}{2},$ $\displaystyle\hat{H}_{j}\rvert_{(m\pm
1/2)}=\sum_{l=1}^{2}\left\\{\left(\hbox{\vbox{\hrule height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\rho_{l}h_{l}$\kern-1.00006pt}}}\right)^{(m\pm\frac{1}{2})}\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("a")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{a}{\tmpbox}_{lj}\rvert_{(m\pm\frac{1}{2})}\right\\},$

where

$\hat{a}_{lj}\rvert_{(m\pm\frac{1}{2})}=\Gamma\left\\{\epsilon\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("\left(\frac{\partial\phi_{l}}{\partial
x_{j}}\right)")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{\left(\frac{\partial\phi_{l}}{\partial
x_{j}}\right)}{\tmpbox}\rvert_{(m\pm\frac{1}{2})}-\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\phi$\kern-1.00006pt}}}_{l}^{(m\pm\frac{1}{2})}(1-\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\phi$\kern-1.00006pt}}}_{l}^{(m\pm\frac{1}{2})})\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$n$\kern-1.00006pt}}}_{lj}^{(m\pm\frac{1}{2})}\right\\}.$

The numerical fluxes in Eq. (59), in addition to Eq. (58), provides a
consistent set of discretization for a two-phase system.

Here, a quadratic-split form is used for $\hat{C}_{j}$ as opposed to the
divergence form because this results in reduced aliasing errors (see Appendix
C). In Eq. (58), the numerical fluxes $\hat{M}_{ij}$ and $\hat{I}_{j}$ assume
cubic-split and quadratic-split forms, respectively; and in Eq. (59),
$\hat{R}_{ij}$ and $\hat{H}_{j}$ assume quintic-split and quartic-split forms,
respectively. In the context of single-phase flows, alternate discrete forms
have been proposed for $\hat{I}_{j}$ in the literature that do not satisfy the
consistency conditions proposed in this work. Chandrashekar (2013) proposed

$\hat{I}_{j}\rvert_{(m\pm\frac{1}{2})}=\hbox{\vbox{\hrule height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\rho$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\left(\frac{2e\rvert_{(m\pm
1)}e\rvert_{m}}{e\rvert_{(m\pm 1)}+e\rvert_{m}}\right)\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u$\kern-1.00006pt}}}_{j}^{(m\pm\frac{1}{2})},$
(60)

as the form of $\hat{I}_{j}$ in their approximate entropy conserving flux
formulation, whereas, Kennedy and Gruber (2008), and Kuya et al. (2018) used

$\hat{I}_{j}\rvert_{(m\pm\frac{1}{2})}=\hbox{\vbox{\hrule height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\rho$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$e$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u$\kern-1.00006pt}}}_{j}^{(m\pm\frac{1}{2})}.$
(61)

We, here, use the form

$\hat{I}_{j}\rvert_{(m\pm\frac{1}{2})}=(\hbox{\vbox{\hrule height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\rho
e$\kern-1.00006pt}}})^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u$\kern-1.00006pt}}}_{j}^{(m\pm\frac{1}{2})},$
(62)

because this is required to satisfy IEC (consistency condition 3). We denote
the $\hat{I}_{j}$ proposed here in Eq. (62) as $QS$, Eq. (60) as $CS-H$, and
Eq. (61) as $CS$, where $QS$ and $CS$ stands for cubic-split and quadratic-
split forms, respectively; and $H$ is used to denote that a harmonic average
is used instead of an arithmetic average in Eq. (60). The different split-
numerical flux forms of $\hat{I}_{j}$ are summarized in Table 1. Similarly,
for $\hat{H}_{j}$, other splitting forms can be derived, but the form in Eq.
(59) is used here because this is required to satisfy IEC (consistency
condition 4).

| CS-H | CS | QS
---|---|---|---
$\hat{I}_{j}\rvert_{(m\pm\frac{1}{2})}$ | $\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.00006pt$\rho$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\left(\frac{2e\rvert_{(m\pm 1)}e\rvert_{m}}{e\rvert_{(m\pm 1)}+e\rvert_{m}}\right)\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.00006pt$u$\kern-1.00006pt}}}_{j}^{(m\pm\frac{1}{2})}$ | $\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.00006pt$\rho$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.00006pt$e$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.00006pt$u$\kern-1.00006pt}}}_{j}^{(m\pm\frac{1}{2})}$ | $(\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.00006pt$\rho e$\kern-1.00006pt}}})^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.00006pt$u$\kern-1.00006pt}}}_{j}^{(m\pm\frac{1}{2})}$
TABLE 1: Split-numerical flux forms for $\hat{I}_{j}$.

We showed that all three $QS$, $CS$, and $CS-H$ forms of $\hat{I}_{j}$ result
in approximate conservation of entropy for single-phase flows, when used
consistently with the other fluxes in Eqs. (58) (Jain and Moin, 2020).
However, only the proposed $QS$ form of $\hat{I}_{j}$ results in conservation
of entropy for two-phase flows. The $CS$ and $CS-H$ forms do not satisfy the
consistency conditions proposed in Section 3.4 and therefore result in
spurious mixing of internal energy at the interface. This is illustrated using
the numerical simulations in Section 4. This alludes to the fact that the set
of consistent numerical fluxes that can be derived for a single-phase Navier-
Stokes system that results in exact conservation of kinetic energy and
approximate conservation of entropy is not unique.

##  4 Coarse-grid simulations

It is well known that the coarse-grid simulation of high-$Re$ turbulent flow
that does not resolve turbulence is a good test of the robustness of a
numerical method. A formulation that does not conserve kinetic energy and
entropy would not be stable because of the lack of grid resolution to support
dissipative scales. These types of under-resolved simulations at high-$Re$
were previously studied, for single-phase flows, by Honein and Moin (2004),
Honein (2005), Mahesh et al. (2004), Hou and Mahesh (2005), Subbareddy and
Candler (2009), Gassner et al. (2016), Kuya et al. (2018), and Coppola et al.
(2019). Here, we propose that the under-resolved simulations at high-$Re$ can
also be used as a test of robustness of a numerical method for two-phase
flows. In this section, we present coarse-grid single-phase and two-phase
simulations, at infinite $Re$, of Taylor-Green vortex and isotropic turbulence
to evaluate the numerical stability of the proposed KEEP scheme. Throughout
this work, an RK4 time-stepping scheme is used for time discretization, and
the temporal discretization errors are assumed to be small. The time-step size
is chosen such that the acoustic time scales of both phases are resolved as
$\Delta t=CFL\min{(\Delta x/c_{1},\Delta x/c_{2})}$, where $CFL$ is the
Courant–Friedrichs–Lewy number chosen to be about $0.65$ in this work.
Throughout this work, the interface-regularization parameters are chosen to be
$\Gamma=|\vec{u}|_{max}$ and $\epsilon=\Delta x$, where $|\cdot|_{max}$
represents the maximum value in the domain.

To evaluate the stability, we look at the time evolution of the total kinetic
energy, $K$, of the system, normalized by the initial kinetic energy $K_{o}$,
defined as

$\frac{K}{K_{o}}=\frac{\int_{\Omega}(\rho u_{i}u_{i}/2)dV}{\int_{\Omega}(\rho
u_{i}u_{i}/2)_{o}dV}=\frac{\sum_{m}(\rho u_{i}u_{i})\rvert_{m}}{\sum_{m}(\rho
u_{i}u_{i})_{o}\rvert_{m}},$ (63)

and the relative change in the total entropy of each of the phases $l$ defined
as

$\frac{\Delta
S_{l}}{S_{lo}}=\frac{\int_{\Omega}\left\\{(\rho_{l}\phi_{l}s_{l})-(\rho_{l}\phi_{l}s_{l})_{o}\right\\}dV}{\lvert\int_{\Omega}(\rho_{l}\phi_{l}s_{l})_{o}dV\rvert}=\frac{\sum_{m}\left\\{\left.(\rho_{l}\phi_{l}s_{l})\right\rvert_{m}-\left.(\rho_{l}\phi_{l}s_{l})_{o}\right\rvert_{m}\right\\}}{\lvert\sum_{m}\left.(\rho_{l}\phi_{l}s_{l})_{o}\right\rvert_{m}\rvert},$
(64)

respectively, where the subscript $o$ denotes that the quantity is evaluated
at the initial time of $t=0$; $dV=\Pi_{j=1}^{3}\Delta x_{j}$ is the cell
volume; and $s_{l}$, for a stiffened-gas EOS, can be defined as

$s_{l}=c_{vl}\ln{\left(\frac{p+\pi_{l}}{\rho_{l}^{\gamma_{l}}}\right)}+s_{lo},$
(65)

where $c_{vl}$ is the specific heat at constant volume for phase $l$.

###  4.1 Inviscid Taylor-Green vortex

Here, an under-resolved simulation of Taylor-Green vortex at infinite $Re$ is
presented. The Taylor-Green vortex is a useful test case to evaluate the
proposed KEEP scheme because it possesses important physical properties of a
turbulent flow, albeit requiring a simple initial setup.

FIGURE 1: A schematic of the domain used in this work for the coarse-grid
simulations of Taylor-Green vortex and isotropic turbulence simulations.

The setup consists of a three-dimensional triply periodic domain of size
$2\pi\times 2\pi\times 2\pi$, discretized into a uniform grid of size
$64\times 64\times 64$. No subgrid-scale model is used in the computation. A
slab of size $2L$ of fluid $1$ is introduced in to the domain surrounded by
fluid $2$, as shown in Figure 1. The initial densities of fluids $1$ and $2$
are chosen to be $\rho_{1o}$ and $\rho_{2o}$, respectively. The viscosities of
fluids $1$ and $2$ are zero: $\mu_{1}=\mu_{2}=0$. The surface tension between
the two fluids is set to zero. The material properties of the fluids $1$ and
$2$ are chosen to be: $\gamma_{1}=\gamma_{2}=\gamma=1.4$, and
$\pi_{1}=\pi_{2}=0$. The length scale $L$ is chosen to be $1$. The initial
conditions for the volume fraction, density, velocities, and pressure are
given by

$\begin{gathered}\phi=1-\frac{1}{2}\left[1+\tanh\left\\{\frac{\sqrt{(x-\pi)^{2}}-L}{2\epsilon}\right\\}\right],\\\
\rho=\rho_{1o}\phi+\rho_{2o}(1-\phi),\\\
u=M_{o}\sin\left(\frac{x}{L}\right)\cos\left(\frac{y}{L}\right)\cos\left(\frac{z}{L}\right),\\\
v=-M_{o}\cos\left(\frac{x}{L}\right)\sin\left(\frac{y}{L}\right)\cos\left(\frac{z}{L}\right),\\\
w=0,\\\ p=\frac{1}{\gamma}+\frac{\rho
M_{o}^{2}}{16}\left\\{\cos\left(\frac{2x}{L}\right)+\cos\left(\frac{2y}{L}\right)\right\\}\left\\{\cos\left(\frac{2z}{L}\right)+2\right\\},\end{gathered}$
(66)

respectively, where $x$, $y$, and $z$ are the spatial coordinates; $u$, $v$,
and $w$ are the velocity components along the $x$, $y$, and $z$ directions,
respectively; and, $M_{o}$ is the initial Mach number of the flow defined
based on the surrounding fluid (fluid $2$) properties. The $Re$ of this case
is, therefore, $Re=\rho UL/\mu=\infty$, where $U$ is the characteristic
velocity. Note that when the densities of the two fluids are equal, the setup
reduces to a single-phase system (or a pseudo-two-phase system).

####  4.1.1 Single-phase flow

Here, the initial densities of the two fluids are chosen to be $\rho_{1o}=1$
and $\rho_{2o}=1$. Since there is no density difference across the interface
between fluids $1$ and $2$ and all other material properties of the two fluids
are either zero or equal, the system of equations discretely reduce to the
single-phase Navier-Stokes system in Section 2.3. Using the KEEP scheme
proposed in this work in Eq. (58), we compare the three numerical flux forms
for $\hat{I}_{j}$, namely, the $CS,CS$-$H,$ and $QS$ forms, introduced in
Section 3.5, and keeping all other fluxes unchanged.

The time evolution of the total kinetic energy and the individual entropy of
fluids $1$ and $2$ are plotted in Figures 2 and 3 for $M_{o}=0.05$ and
$M_{o}=0.2$, respectively. Note that the properties of fluids $1$ and $2$ are
the same and hence the tags $1$ and $2$ represent the same fluid, but that are
in separate regions in the domain. All three methods conserve kinetic energy
and entropy of the system and therefore are stable for single-phase flows.
This further verifies that the consistent set of numerical fluxes that result
in conservation of kinetic energy and entropy is not unique for a single-phase
Navier-Stokes system.

FIGURE 2: Inviscid single-phase Taylor-Green case ($\rho_{2}/\rho_{1}=1$) with
$M_{o}=0.05$. (a) The evolution of the total kinetic energy in the domain $K$
normalized by the initial value of the kinetic energy in the domain $K_{o}$.
(b) The change in entropy of fluid $1$ with time, normalized by the initial
value. (c) The change in entropy of fluid $2$ with time, normalized by the
initial value. Note that, here, fluids $1$ and $2$ are the same and hence the
tags $1$ and $2$ represent the same fluid, but that are in separate regions in
the domain. FIGURE 3: Inviscid single-phase Taylor-Green case
($\rho_{2}/\rho_{1}=1$) with $M_{o}=0.2$. (a) The evolution of the total
kinetic energy in the domain $K$ normalized by the initial value of the
kinetic energy in the domain $K_{o}$. (b) The change in entropy of fluid $1$
with time, normalized by the initial value. (c) The change in entropy of fluid
$2$ with time, normalized by the initial value. Note that, here, fluids $1$
and $2$ are the same and hence the tags $1$ and $2$ represent the same fluid,
but that are in separate regions in the domain.

####  4.1.2 Two-phase flow

In this section, the initial densities of the two fluids are chosen to be
$\rho_{1o}=0.1$ and $\rho_{2o}=1$ and the setup is, therefore, a two-phase
system with non-unity density difference across the interface between fluids
$1$ and $2$. Now using the KEEP scheme proposed in this work in Eqs.
(58)-(59), we again compare the three numerical flux forms for $\hat{I}_{j}$,
keeping all other fluxes unchanged. Snapshots of the flow are shown in Figure
6 at various times, illustrating the breakup of the fluid slab into bubbles
due to the breakdown of the underlying Taylor-Green vortex.

The time evolution of the total kinetic energy and the individual entropy of
fluids $1$ and $2$ are plotted in Figures 4 and 5 for $M_{o}=0.05$ and
$M_{o}=0.2$, respectively. Unlike for the single-phase system, here only the
$QC$ form proposed in this work conserves kinetic energy and entropy, and is
therefore stable for two-phase flows. The $CS$ and $CS$-$H$ forms on the other
hand, do not satisfy the consistency conditions proposed in Section 3.4, and
therefore, do not conserve kinetic energy and entropy. The simulations with
these forms diverge. Note that, here, different numerical flux forms for
$\hat{I}_{j}$ are compared because these forms are used before in the
literature in the context of single-phase flows. Similarly, a cubic-split form
for $\hat{H}_{j}$ can be constructed instead of the quadratic-split form
proposed in Eq. (59) and can be shown that the resulting simulation would not
be stable since the cubic-split form for $\hat{H}_{j}$ does not satisfy the
consistency conditions. Appendix D discusses other possible inconsistent flux
formulations and the effect of these formulations on the numerical stability.

FIGURE 4: Inviscid two-phase Taylor-Green case with $\rho_{2}/\rho_{1}=10$ and
$M_{o}=0.05$. (a) The evolution of the total kinetic energy in the domain $K$
with time normalized by the initial value of the kinetic energy in the domain
$K_{o}$. (b) The change in entropy of fluid $1$ with time, normalized by the
initial value. (c) The change in entropy of fluid $2$ with time, normalized by
the initial value. FIGURE 5: Inviscid two-phase Taylor-Green case with
$\rho_{2}/\rho_{1}=10$ and $M_{o}=0.2$. (a) The evolution of the total kinetic
energy in the domain $K$ with time normalized by the initial value of the
kinetic energy in the domain $K_{o}$. (b) The change in entropy of fluid $1$
with time, normalized by the initial value. (c) The change in entropy of fluid
$2$ with time, normalized by the initial value. FIGURE 6: Snapshots of the
inviscid two-phase Taylor-Green vortex at various times with the proposed KEEP
scheme in this work. The solid surface represents the interface ($\phi=0.5$
contour) between fluid $1$ and fluid $2$. The translucent surface represents
an isosurface of the Q-criterion colored by the vorticity magnitude.

Finally, the proposed KEEP scheme is verified for different initial Mach
numbers and densities of the fluids in Figure 7. As $M_{o}$ increases, the
kinetic energy does not stay constant due to the exchange of energy between
the kinetic energy and internal energy through reversible pressure work.

FIGURE 7: Inviscid two-phase Taylor-Green case with $\rho_{2}/\rho_{1}=1000$
for different Mach numbers (a) The evolution of the total kinetic energy in
the domain $K$ with time normalized by the initial value of the kinetic energy
in the domain $K_{o}$. (b) The change in entropy of fluid $1$ with time,
normalized by the initial value. (c) The change in entropy of fluid $2$ with
time, normalized by the initial value.

###  4.2 Two-phase isotropic turbulence

In this section, an under-resolved simulation of two-phase isotropic
turbulence at infinite $Re$ is presented. The setup consists of a three-
dimensional triply periodic domain of size $2\pi\times 2\pi\times 2\pi$,
discretized into a uniform grid of size $32\times 32\times 32$. No subgrid-
scale model is used in the computation. A slab of size $2L$ of fluid $1$ is
introduced in to the domain surrounded by fluid $2$, as shown in Figure 1. The
initial densities of fluids $1$ and $2$ are chosen to be $\rho_{1o}=0.1$ and
$\rho_{2o}=1$, respectively. The viscosities of fluids $1$ and $2$ are zero:
$\mu_{1}=\mu_{2}=0$. The material properties of the fluids $1$ and $2$ are
chosen to be: $\gamma_{1}=\gamma_{2}=\gamma=1.4$, and $\pi_{1}=\pi_{2}=0$. The
surface tension between the two fluids is set to zero. The length scale $L$ is
chosen to be $1$. The initial conditions for the volume fraction and density
are the same as the ones in Section 4.1. An initial energy spectrum used by
Honein (2005) defines the initial velocity field.

FIGURE 8: Inviscid two-phase isotropic turbulence case with
$\rho_{2}/\rho_{1}=10$ and $M_{to}=0.07$. (a) The evolution of the total
kinetic energy in the domain $K$ with time normalized by the initial value of
the kinetic energy in the domain $K_{o}$. (b) The change in entropy of fluid
$1$ with time, normalized by the initial value. (c) The change in entropy of
fluid $2$ with time, normalized by the initial value.

Using the KEEP scheme proposed in this work in Eqs. (58)-(59), we again
compare the three numerical flux forms for $\hat{I}_{j}$, keeping all other
fluxes unchanged. The time evolution of the total kinetic energy and the
individual entropy of fluids $1$ and $2$ are plotted in Figure 8 for an
initial turbulent Mach number of $M_{to}=0.07$. Similar to the two-phase
Taylor-Green system in Section 4.1.2, only the $QC$ form proposed in this work
conserves kinetic energy and entropy, and is therefore stable. The $CS$ and
$CS$-$H$ forms, however, do not satisfy the consistency conditions and
therefore do not conserve kinetic energy and entropy. The simulations with
these forms diverge.

##  5 Higher-resolution simulations of droplet-laden isotropic turbulence

In this section, simulations of droplet-laden isotropic turbulence is
presented to evaluate the accuracy of the method at finite $Re$. The setup
consists of a three-dimensional triply periodic domain of size
$\pi/2\times\pi/2\times\pi/2$, and is discretized into a uniform grid of size
$256^{3}$. No subgrid-scale model is used in the computation. A 1000,
initially spherical, droplets of fluid $1$ is introduced to the domain
surrounded by fluid $2$. The initial Taylor-Reynolds number of the flow is
$Re_{\lambda}\approx 75$ and the initial turbulent Mach number is
$M_{to}=0.07$, based on the surrounding fluid properties. The initial diameter
of the seeded droplets is $d_{o}\approx 2\lambda$. Here, $Re_{\lambda}$ was
chosen to match the incompressible droplet-laden isotropic turbulence study by
Dodd and Ferrante (2016).

The initial density of fluid $2$ is $\rho_{2o}=1$ and the density of the fluid
$1$ is varied: $\rho_{1o}=1,10,100,1000$. A reference case of single-phase
homogeneous isotropic turbulence (HIT) is also simulated with only fluid $2$
in the domain. The kinematic viscosities of fluids $1$ and $2$ are assumed to
be the same. The material properties of the fluids $1$ and $2$ are chosen to
be: $\gamma_{1}=\gamma_{2}=\gamma=1.4$, and $\pi_{1},\pi_{2}$ are chosen such
that the values of $M_{to}$ are identical in two fluids. The surface tension,
$\sigma$, between the two fluids is chosen such that the initial turbulent
Weber number is $We_{rms}=\rho_{2o}U_{rms}^{2}d_{o}/\sigma=1$. An initial
energy spectrum used by Honein (2005) defines the initial velocity field.

Snapshots from the simulations of droplet-laden isotropic turbulence for all
four values of $\rho_{1}/\rho_{2}$, at $t=20$, is presented in Figure 9. The
only difference between the single-phase case and the two-phase case with
$\rho_{1}/\rho_{2}=1$ is the presence of the surface tension effects at the
droplet interface. This inhibits breakup of the droplet fluid, and hence, the
droplet size increases with time. With an increase in the density of the
droplet fluid, it is evident that the sizes of the drops are smaller. This is
due to the increase in the inertial effects over the surface tension effects
that counter the breakup process.

FIGURE 9: Snapshots of the droplet-laden isotropic turbulence simulations for
various values of $\rho_{1}/\rho_{2}$ at $t=20$. The solid surface represents
the interface ($\phi=0.5$ contour) between fluid $1$ and fluid $2$. The
translucent surface represents an isosurface of the Q-criterion colored by the
vorticity magnitude.

The time evolution of the total kinetic energy is plotted in Figure 10 for all
four values of $\rho_{1}/\rho_{2}$ along with the reference single-phase case.
The relative enhancement of the decay rate of total kinetic energy at the
early times, due to the presence of drops (see the differences in $K$ between
the two-phase case of $\rho_{1}/\rho_{2}=1$ and the single-phase case), and
with an increase in the density of the droplet fluid, is in very good
agreement with the observations made by Dodd and Ferrante (2016). However, for
the $\rho_{1}/\rho_{2}=1$ case, the decay of $K$ slows down compared to the
single-phase case at later times and a crossover of $K$ at $t\approx 10$ can
be seen in Figure 10. This behavior can be explained using the effect of
surface tension on the flow. Writing the evolution equation for $K$ as

$\frac{dK}{dt}=-E+\Psi$ (67)

where $E$ and $\Psi$ are the total dissipation rate and the rate of change of
surface energy, which scales as

$\Psi\sim-\frac{dA}{dt},$ (68)

where $A$ is the total surface area of the drops in the domain. As the droplet
sizes increase with time for the $\rho_{1}/\rho_{2}=1$ case at later time in
the simulation, the total surface area $A$ decreases. Therefore, $\Psi$ is
positive and acts as a source of $K$. This results in slower decay of $K$ for
the two-phase case with $\rho_{1}/\rho_{2}=1$ compared to the single-phase
case, and hence, a crossover in $K$ is seen at $t\approx 10$. A similar
behavior is expected for the case of $\rho_{1}/\rho_{2}=10$, but at a much
later time, due to a relative slower increase in the droplet sizes compared to
the case of $\rho_{1}/\rho_{2}=1$. Figure 10 confirms this behavior, and a
crossover in $K$ can be seen for this case at $t\approx 18$.

FIGURE 10: The time evolution of the total kinetic energy in the domain $K$,
normalized by the initial value of the kinetic energy in the domain $K_{o}$,
for the droplet-laden isotropic turbulence simulations for various values of
$\rho_{1}/\rho_{2}$ and a reference single-phase case.

##  6 Summary and remarks

In this work, we developed a general framework for the derivation of
consistent numerical fluxes for compressible single-phase and two-phase flows,
and proposed a set of consistency conditions between the numerical fluxes of
volume fraction, mass, momentum, kinetic energy, and internal energy. The
proposed consistency conditions between the fluxes are required for achieving
discrete conservation of global kinetic energy, preservation of local kinetic
energy, and for maintaining interface-equilibrium condition. For
incompressible flows, conservation of global kinetic energy is known to be a
sufficient condition to achieve stable numerical simulations. However, for
compressible flows, in addition to conserving global kinetic energy, more
constraints such as local kinetic energy preservation and interface-
equilibrium condition are required to maintain an approximate discrete
conservation of entropy, and to achieve stable numerical simulations. The
proposed consistency conditions for the numerical fluxes are general enough
that they can be adopted for other multiphysics flow problems, such as
multicomponent flows.

Using the proposed framework, a numerical scheme that satisfies the
consistency conditions was derived for both single-phase and two-phase flows
and was verified that it results in an exact conservation of the kinetic
energy and approximate conservation of entropy (a KEEP scheme) in the absence
of pressure work, viscosity, thermal diffusion effects, and time-differencing
errors. We also showed that a KEEP scheme is not unique for single-phase
flows, and various forms of numerical fluxes can be derived that satisfy the
consistency conditions and hence would result in exact conservation of kinetic
energy and approximate conservation of entropy.

A conservative diffuse-interface method with a five-equation model was used as
the interface-capturing method for modeling the system of compressible two-
phase flows in this work. However, the proposed consistency conditions are
general, and can be used with any other interface-capturing method to derive a
KEEP scheme.

We used coarse-grid simulations of single-phase and two-phase Taylor-Green
vortex and isotropic turbulence at infinite $Re$, to test the numerical
stability of the proposed scheme. The observations from the numerical
experiments verified that the proposed scheme results in conservation of
kinetic energy and entropy, and hence, the scheme is superior in terms of
maintaining stability for long-time numerical integrations. A higher-
resolution simulation of droplet-laden decaying isotropic turbulence case (at
finite $Re$) is also presented at the end, and the effect of presence of
droplets on the flow is studied. This test case illustrates the stability and
accuracy of the proposed method in a complex setting.

## Acknowledgments

S. S. J. was supported by the Franklin P. and Caroline M. Johnson Stanford
Graduate Fellowship, and the authors acknowledge the support from the
Predictive Science Academic Alliance Program (PSAAP) III at Stanford
University. The authors also acknowledge the Argonne Leadership Computing
Challenge award 2021-22 which provided access to the Theta supercomputer,
which was used for simulations in this work. A preliminary version of this
work has been published as a technical report in the annual publication of the
Center for Turbulence Research (Jain and Moin, 2020) and is available
online111http://web.stanford.edu/group/ctr/ResBriefs/2020/29_Jain.pdf. S. S.
J. is thankful for Mr. Kihiro Bando, for reading the manuscript and for his
comments that helped improve the manuscript.

## Appendix A: Entropy change associated with interface-regularization process

To quantify the entropy change associated with the interface-regularization
process and to illustrate the conservation of entropy in the limit of
equilibrium interface state, a simple numerical test is performed. In this
test case, the interface thickness is varied in a stationary setup such that
the only terms that are non-zero are the interface-regularization terms,
facilitating the quantification of entropy change associated with these terms.

The setup consists of an infinite slab of fluid $1$ of width $2L$, where
$L=0.25$, surrounded by fluid $2$ in a two-dimensional periodic square domain,
as shown in Figure 11(a). The domain has dimensions of
$[-0.5,0.5]\times[-0.5,0.5]$, discretized into a uniform grid of size
$64\times 64$. The initial conditions for the volume fraction, density,
velocities, and pressure are given by
$\phi=1-{1}/{2}\left[1+\tanh\left\\{{(|x|-L)}/{(2\epsilon_{o})}\right\\}\right]$,
$\rho=\rho_{1o}\phi+\rho_{2o}(1-\phi)$, $u=0$, $v=0$, and $p=1$, respectively;
$\rho_{1o}=1$ and $\rho_{2o}=1000$ are the initial densities of the phases $1$
(air) and $2$ (water), respectively. The viscosities are taken to be zero:
$\mu_{1}=\mu_{2}=0$. The surface tension is set to zero.

FIGURE 11: A stationary slab of fluid $1$ in fluid $2$. (a) initial setup, (b)
A line plot of $\phi$ along $y=0$ from $x=0$ to $x=0.5$, showing the initial
and final states for three different scenarios of $\epsilon$. FIGURE 12: The
relative change in total entropy $\Delta S_{l}/S_{lo}$ of (a) fluid $1$, and
(b) fluid $2$, for three different values of $\epsilon$.

Three different scenarios are considered with different interface equilibrium
states that are controlled using the interface thickness parameter,
$\epsilon$: (1) $\epsilon=\Delta x$, (2) $\epsilon=1.25\Delta x$, and (3)
$\epsilon=0.75\Delta x$. In all the three cases, volume fraction is
initialized using the hyperbolic tangent function with $\epsilon_{o}=\Delta
x$. Hence, the interface is already in an equilibrium state in case $1$. But
for cases $2$ and $3$, the interface thickness increases and decreases,
respectively, until it reaches a new equilibrium state. The relative change in
total entropy of each phase $l$ is computed as a function of time and are
shown in Figures 12(a,b) for all three cases. Since in case $1$, the interface
is already in an equilibrium state, the change in entropy of both fluids $1$
and $2$ are negligible. For cases $2$ and $3$, where the interface thickness
increases and decreases, respectively, there is a small change (increase) in
entropy. But once the interface reaches a new equilibrium state in cases $2$
and $3$, the entropy of both fluids does not change, as expected. The initial
and final states of the interface at times, $t=0$ and $t=1$, respectively, are
shown in Figure 11(b) for all three scenarios. This numerical demonstration
verifies that the entropy is conserved in the limit of equilibrium interface
state and the change in entropy due to the interface-regularization process is
small and positive. Although the change in entropy due to the interface-
regularization terms are small and essentially negligible, this does not
warrant the use of inconsistent numerical fluxes for these terms. Failing to
satisfy the consistency conditions (KEEP scheme) presented in this work, these
interface-regularization terms could spuriously contribute to the kinetic
energy and entropy of the system (see, Appendix D).

## Appendix B: KEEP scheme satisfies the interface-equilibrium condition

The IEC provides a consistency condition to check and eliminate the forms of
the numerical discretizations that contribute spuriously to the solution.
Satisfying the IEC is crucial for a method to conserve kinetic energy and
entropy.

######  Lemma 6.1.

The proposed KEEP scheme satisfies the IEC defined in Section 3.3.

###### Proof.

### Part (a). Uniform velocity

Consider a finite-volume discretization of the total mass balance equation in
Eq. (10)

$\rho\rvert_{m}^{k+1}-\rho\rvert_{m}^{k}=-\Delta
t\left(\frac{\hat{C}_{j}\rvert_{(m+\frac{1}{2})}-\hat{C}_{j}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)^{k}+\Delta
t\left(\frac{\hat{F}_{j}\rvert_{(m+\frac{1}{2})}-\hat{F}_{j}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)^{k},$ (69)

where $k$ is the time-step index, and $m$ is the grid index. Now, consider a
finite-volume discretization of the momentum equation in Eq. (3) and ignoring
viscous, surface tension, and gravity terms, and assuming
$u_{i}\rvert_{m}^{k}=u_{0}$ and $p\rvert_{m}^{k}=p_{0}$, $\forall m$

$(\rho u_{i})\rvert_{m}^{k+1}-\rho\rvert_{m}^{k}u_{0}=-\Delta
t\left(\frac{\hat{M}_{ij}\rvert_{(m+\frac{1}{2})}-\hat{M}_{ij}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)^{k}+\Delta
t\left(\frac{\hat{R}_{ij}\rvert_{(m+\frac{1}{2})}-\hat{R}_{ij}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)^{k}.$ (70)

Utilizing the consistency conditions in Eqs. (43), (44), and assuming
$u_{i}\rvert_{m}^{k}=u_{0}$, $\forall m$, the discrete momentum equation can
be written as

$(\rho u_{i})\rvert_{m}^{k+1}-\rho\rvert_{m}^{k}u_{0}=-\Delta
t\left(\frac{\hat{C}_{j}\rvert_{(m+\frac{1}{2})}-\hat{C}_{j}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)^{k}u_{0}+\Delta
t\left(\frac{\hat{F}_{j}\rvert_{(m+\frac{1}{2})}-\hat{F}_{j}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)^{k}u_{0}.$ (71)

Subtracting this from the discrete mass balance equation in Eq. (69), gives
$u_{i}\rvert_{m}^{k+1}=u_{0}$, $\forall m$.

### Part (b). Uniform pressure

Consider a finite-volume discretization of the total energy equation [Eq. (4)]
without the viscous, surface tension, and gravity terms. Subtracting the
discrete version of the kinetic energy equation in Eq. (13) from the discrete
total energy equation (note that, this step requires that the consistency
conditions in Eqs. (50), (51) are satisfied), we arrive at

$\left(\rho e\right)\rvert_{m}^{k+1}-\left(\rho e\right)\rvert_{m}^{k}=-\Delta
t\left(\frac{\hat{I}_{j}\rvert_{(m+\frac{1}{2})}-\hat{I}_{j}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)^{k}+\Delta
t\left(\frac{\hat{H}_{j}\rvert_{(m+\frac{1}{2})}-\hat{H}_{j}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)^{k}$ (72)

Now, substituting for $\hat{I}_{j}\rvert_{(m\pm\frac{1}{2})}$ and
$\hat{H}_{j}\rvert_{(m\pm\frac{1}{2})}$ from Eqs. (58)-(59); expressing
$h_{l}$ in terms of $p$ and $\rho_{l}$ using Eq. (9) as
$h_{l}=(p-\beta_{l})/\alpha_{l}+p$; expressing the mixture internal energy in
terms of the individual species energies as $\rho
e=\sum_{l=1}^{2}\rho_{l}e_{l}\phi_{l}$; and assuming
$u_{i}\rvert_{m}^{k}=u_{0}$ and $p\rvert_{m}^{k}=p_{0}$, we obtain

$\displaystyle\sum_{l=1}^{2}\left(\phi_{l}\rho_{l}e_{l}\right)\rvert_{m}^{k+1}-\sum_{l=1}^{2}\left(\phi_{l}\rho_{l}e_{l}\right)\rvert_{m}^{k}=-\Delta
t\sum_{l=1}^{2}\left\\{\frac{\left(\hbox{\vbox{\hrule height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\rho_{l}e_{l}\phi_{l}$\kern-1.00006pt}}}^{(m+\frac{1}{2})}\right)-\left(\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\rho_{l}e_{l}\phi_{l}$\kern-1.00006pt}}}^{(m-\frac{1}{2})}\right)}{\Delta
x_{j}}\right\\}^{k}u_{0}$ (73) $\displaystyle+\Delta
t\left\\{\sum_{l=1}^{2}\left(\frac{p_{0}-\beta_{l}}{\alpha_{l}}\right)\left(\frac{\hat{a}_{lj}\rvert_{(m+\frac{1}{2})}-\hat{a}_{lj}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)\right\\}^{k},$

and expressing $e_{l}$ in terms of $p$ using the EOS results in the
discretized equation for pressure

$\displaystyle\left(\sum_{l=1}^{2}\frac{\phi_{l}\rvert_{m}}{\alpha_{l}}\right)^{k+1}p\rvert_{m}^{k+1}-\left(\sum_{l=1}^{2}\frac{\phi_{l}\rvert_{m}\beta_{l}}{\alpha_{l}}\right)^{k+1}-\left(\sum_{l=1}^{2}\frac{\phi_{l}\rvert_{m}}{\alpha_{l}}\right)^{k}p_{0}+\left(\sum_{l=1}^{2}\frac{\phi_{l}\rvert_{m}\beta_{l}}{\alpha_{l}}\right)^{k}$
(74) $\displaystyle=-\Delta
t\left[\sum_{l=1}^{2}\left\\{\frac{\left(\frac{\hbox{\vbox{\hrule
height=0.5pt\kern
1.50694pt\hbox{\kern-0.70004pt$\phi$\kern-0.70004pt}}}_{l}^{(m+\frac{1}{2})}}{\alpha_{l}}\right)p_{0}-\left(\frac{\hbox{\vbox{\hrule
height=0.5pt\kern
1.50694pt\hbox{\kern-0.70004pt$\phi$\kern-0.70004pt}}}_{l}^{(m+\frac{1}{2})}\beta_{l}}{\alpha_{l}}\right)-\left(\
\frac{\hbox{\vbox{\hrule height=0.5pt\kern
1.50694pt\hbox{\kern-0.70004pt$\phi$\kern-0.70004pt}}}_{l}^{(m-\frac{1}{2})}}{\alpha_{l}}\right)p_{0}+\left(\frac{\hbox{\vbox{\hrule
height=0.5pt\kern
1.50694pt\hbox{\kern-0.70004pt$\phi$\kern-0.70004pt}}}_{l}^{(m-\frac{1}{2})}\beta_{l}}{\alpha_{l}}\right)}{\Delta
x_{j}}\right\\}^{k}\right]u_{0}$ $\displaystyle+\Delta
t\left\\{\sum_{l=1}^{2}\left(\frac{p_{0}-\beta_{l}}{\alpha_{l}}\right)\left(\frac{\hat{a}_{lj}\rvert_{(m+\frac{1}{2})}-\hat{a}_{lj}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)\right\\}^{k}.$

Let $L(\phi_{l})$ be a finite-volume discretization of the volume fraction
advection equation for phase $l$ [Eq. (1)], given as

$L(\phi_{l})\equiv\phi_{l}\rvert_{m}^{k+1}-\phi_{l}\rvert_{m}^{k}=-\Delta
t\left(\frac{\hbox{\vbox{\hrule height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\phi$\kern-1.00006pt}}}_{l}^{(m+\frac{1}{2})}-\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\phi$\kern-1.00006pt}}}_{l}^{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)^{k}u_{0}+\Delta
t\left(\frac{\hat{a}_{lj}\rvert_{(m+\frac{1}{2})}-\hat{a}_{lj}\rvert_{(m-\frac{1}{2})}}{\Delta
x_{j}}\right)^{k}.$ (75)

Subtracting Eq. (74) from the equation
$\Big{(}\sum_{l=1}^{2}{L(\phi_{l})}/{\alpha_{l}}\Big{)}p_{0}-\Big{(}\sum_{l=1}^{2}{L(\phi_{l})\beta_{l}}/{\alpha_{l}}\Big{)}$,
results in $p_{i}\rvert_{m}^{k+1}=p_{0}$, $\forall m$, which concludes the
proof. ∎

Here, a stiffened gas EOS has been used in Eq. (74); however, a more general
cubic EOS can also be used to prove IEC. Note that when Eq. (74) was
subtracted from
$\Big{(}\sum_{l=1}^{2}{L(\phi_{l})}/{\alpha_{l}}\Big{)}p_{0}-\Big{(}\sum_{l=1}^{2}{L(\phi_{l})\beta_{l}}/{\alpha_{l}}\Big{)}$,
it was assumed that the split-numerical flux forms used in $\hat{I}_{j}$ and
$\hat{\Phi}_{j}$ were identical. Without this, it is not possible to show that
$p_{i}\rvert_{m}^{k+1}=p_{0}$, and therefore, the IEC cannot be proved.
Section 3.2 describes that the split-numerical flux forms used in
$\hat{\Phi}_{j}$ and $\hat{C}_{j}$ should be identical. Therefore, from the
transitive property, the split-numerical flux forms used in $\hat{I}_{j}$ and
$\hat{C}_{j}$ should be identical for the IEC to be satisfied. Similarly, the
split-numerical flux forms used in $\hat{H}_{j}$ and $\hat{F}_{j}$ should be
identical for the IEC to be satisfied.

## Appendix C: Effect of flux splitting

In this section, the effect of skew-symmetric splitting of the numerical
fluxes on the non-linear stability of the method is evaluated. Particularly,
the effect of splitting of mass flux $\hat{C}_{j}$ is assessed for two-phase
flows. In Section 3.5, a quadratic split form

$\hat{C}_{j}\rvert_{(m\pm\frac{1}{2})}=\hbox{\vbox{\hrule height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\rho$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u$\kern-1.00006pt}}}_{j}^{(m\pm\frac{1}{2})},$
(76)

was used. An alternative would be to use the divergence form for the mass flux
$\hat{C}_{j}$ as

$\hat{C}_{j}\rvert_{(m\pm\frac{1}{2})}=\hbox{\vbox{\hrule height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\rho
u_{j}$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}.$ (77)

| divergence $\hat{C}_{j}$ | quadratic split $\hat{C}_{j}$
---|---|---
$\hat{C}_{j}\rvert_{(m\pm\frac{1}{2})}$ | $\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.00006pt$\rho u_{j}$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}$ | $\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.00006pt$\rho$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.00006pt$u$\kern-1.00006pt}}}_{j}^{(m\pm\frac{1}{2})}$
$\hat{M}_{ij}\rvert_{(m\pm\frac{1}{2})}$ | $\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.00006pt$\rho u_{j}$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.00006pt$u$\kern-1.00006pt}}}_{i}^{(m\pm\frac{1}{2})}$ | $\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.00006pt$\rho$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.00006pt$u$\kern-1.00006pt}}}_{j}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.00006pt$u$\kern-1.00006pt}}}_{i}^{(m\pm\frac{1}{2})}$
$\hat{K}_{j}\rvert_{(m\pm\frac{1}{2})}$ | $\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.00006pt$\rho u_{j}$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\frac{u_{i}\rvert_{(m\pm 1)}u_{i}\rvert_{m}}{2}$ | $\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.00006pt$\rho$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-1.00006pt$u$\kern-1.00006pt}}}_{j}^{(m\pm\frac{1}{2})}\frac{u_{i}\rvert_{(m\pm 1)}u_{i}\rvert_{m}}{2}$
TABLE 2: Consistent momentum and kinetic energy fluxes for the divergence and
quadratic-split form of the mass flux.

Using the consistency conditions in Section 3.4, the consistent momentum and
kinetic energy fluxes for the divergence and quadratic-split forms of
$\hat{C}_{j}$ can be written as shown in Table 2. To compare the two
formulations, under-resolved simulations of two-phase Taylor-Green vortex at
infinite $Re$ presented in Section 4.1.2 are repeated using these
formulations, while keeping all other fluxes unchanged. The time evolution of
the total kinetic energy and the individual entropy of fluids 1 and 2 are
plotted in Figure 13 for $M_{o}=0.05$. The initial densities of the two fluids
are $\rho_{1o}=0.1$ and $\rho_{2o}=1$. The simulation that uses the divergence
form of the flux for $\hat{C}_{j}$ is not stable and diverges, albeit
$\hat{M}_{ij}$ and $\hat{K}_{j}$ satisfy the consistency conditions; whereas
the simulation with the proposed quadratic-split form of the flux for
$\hat{C}_{j}$ is stable due to the reduced aliasing errors associated with
these schemes (Blaisdell et al., 1996, Chow and Moin, 2003, Kennedy and
Gruber, 2008).

FIGURE 13: Inviscid two-phase Taylor-Green case with $\rho_{2}/\rho_{1}=10$
and $M_{o}=0.05$. (a) The evolution of the total kinetic energy in the domain
$K$ with time normalized by the initial value of the kinetic energy in the
domain $K_{o}$. (b) The change in entropy of fluid $1$ with time, normalized
by the initial value. (c) The change in entropy of fluid $2$ with time,
normalized by the initial value.

## Appendix D: Inconsistent flux formulations

In Section 4, the effect of inconsistent flux formulations for $\hat{I}_{j}$,
that violate the IEC, on the non-linear stability was already examined in
detail. Here, a similar study is performed with three different inconsistent
flux formulations for $\hat{M}_{ij}$, $\hat{K}_{j}$, and $\hat{T}_{j}$, and
are compared against the proposed KEEP formulation in Section 3.5.

An inconsistent flux formulation for $\hat{M}_{ij}$ can be written as

$\hat{M}_{ij}\rvert_{(m\pm\frac{1}{2})}=\hbox{\vbox{\hrule height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\rho
u_{j}$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u$\kern-1.00006pt}}}_{i}^{(m\pm\frac{1}{2})}.$
(78)

This form of $\hat{M}_{ij}$ in Eq. (80) will violate the consistency condition
1 in Section 3.4, when all other fluxes are kept unchanged as in Eqs. (58),
(59). An inconsistent formulation for $\hat{K}_{j}$ can be written as

$\hat{K}_{j}\rvert_{(m\pm\frac{1}{2})}=\hbox{\vbox{\hrule height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\rho$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$u$\kern-1.00006pt}}}_{j}^{(m\pm\frac{1}{2})}\frac{u_{i}\rvert_{(m\pm\frac{1}{2})}u_{i}\rvert_{(m\pm\frac{1}{2})}}{2}.$
(79)

This form of $\hat{K}_{j}$ in Eq. (80) will also violate the consistency
condition 1 in Section 3.4, when all other fluxes are kept unchanged as in
Eqs. (58), (59). An inconsistent formulation for $\hat{T}_{j}$ can be written
as

$\hat{T}_{j}\rvert_{(m\pm\frac{1}{2})}=\left\\{\sum_{l=1}^{2}\left(\hbox{\vbox{\hrule
height=0.5pt\kern
2.15277pt\hbox{\kern-1.00006pt$\rho_{l}$\kern-1.00006pt}}}^{(m\pm\frac{1}{2})}\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("a")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{a}{\tmpbox}_{lj}\rvert_{(m\pm\frac{1}{2})}\right)\right\\}\frac{u_{i}\rvert_{(m\pm\frac{1}{2})}u_{i}\rvert_{(m\pm\frac{1}{2})}}{2}.$
(80)

This form of $\hat{T}_{j}$ in Eq. (80) will also violate the consistency
condition 1 in Section 3.4, when all other fluxes are kept unchanged as in
Eqs. (58), (59).

To compare the three inconsistent formulations for $\hat{M}_{ij}$,
$\hat{K}_{j}$, and $\hat{T}_{j}$ in Eqs. (78)-(80) with the proposed KEEP
formulation in this work, under-resolved simulations of two-phase Taylor-Green
vortex at infinite $Re$ presented in Section 4.1.2 are repeated using these
formulations. The time evolution of the total kinetic energy and the
individual entropy of fluids 1 and 2 are plotted in Figure 14 for
$M_{o}=0.05$. The initial densities of the two fluids are $\rho_{1o}=0.1$ and
$\rho_{2o}=1$. The simulations that use the inconsistent formulations in Eqs.
(78)-(80) are not stable because they don’t satisfy the consistency
conditions, and hence they diverge; whereas the simulation with the proposed
KEEP formulation in Section 3.5, that satisfies the consistency conditions, is
stable.

FIGURE 14: Inviscid two-phase Taylor-Green case with $\rho_{2}/\rho_{1}=10$
and $M_{o}=0.05$. (a) The evolution of the total kinetic energy in the domain
$K$ with time normalized by the initial value of the kinetic energy in the
domain $K_{o}$. (b) The change in entropy of fluid $1$ with time, normalized
by the initial value. (c) The change in entropy of fluid $2$ with time,
normalized by the initial value.

## References

  * Sallam et al. (2006) K. A. Sallam, C. Aalburg, G. Faeth, K.-C. Lin, C. Carter, T. A. Jackson,  Primary breakup of round aerated-liquid jets in supersonic crossflows,  Atomization and Sprays 16 (2006).
  * Law (1982) C. Law,  Recent advances in droplet vaporization and combustion,  Progress in energy and combustion science 8 (1982) 171–201.
  * Plesset and Prosperetti (1977) M. S. Plesset, A. Prosperetti,  Bubble dynamics and cavitation,  Annual review of fluid mechanics 9 (1977) 145–185.
  * Youngs (1984) D. L. Youngs,  Numerical simulation of turbulent mixing by rayleigh-taylor instability,  Physica D: Nonlinear Phenomena 12 (1984) 32–44.
  * Brouillette (2002) M. Brouillette,  The richtmyer-meshkov instability,  Annual Review of Fluid Mechanics 34 (2002) 445–468.
  * Morinishi et al. (1998) Y. Morinishi, T. S. Lund, O. V. Vasilyev, P. Moin,  Fully conservative higher order finite difference schemes for incompressible flow,  Journal of computational physics 143 (1998) 90–124.
  * Honein (2005) A. E. Honein, Numerical aspects of compressible turbulence simulations, Ph.D. thesis, Stanford University, 2005\.
  * Feiereisen (1983) W. J. Feiereisen, Numerical simulation of a compressible, homogeneous, turbulent shear flow, Ph.D. thesis, Stanford University, 1983.
  * Blaisdell et al. (1996) G. Blaisdell, E. Spyropoulos, J. Qin,  The effect of the formulation of nonlinear terms on aliasing errors in spectral methods,  Applied Numerical Mathematics 21 (1996) 207–219.
  * Kravchenko and Moin (1997) A. Kravchenko, P. Moin,  On the effect of numerical errors in large eddy simulations of turbulent flows,  Journal of computational physics 131 (1997) 310–322.
  * Ducros et al. (2000) F. Ducros, F. Laporte, T. Soulères, V. Guinot, P. Moinat, B. Caruelle,  High-order fluxes for conservative skew-symmetric-like schemes in structured meshes: application to compressible flows,  Journal of Computational Physics 161 (2000) 114–139.
  * Lee (1993) S. Lee,  Large eddy simulation of shock turbulence interaction,  Center for Turbulence Research, Annual Research Briefs (1993) 73–84.
  * Nagarajan et al. (2003) S. Nagarajan, S. K. Lele, J. H. Ferziger,  A robust high-order compact method for large eddy simulation,  Journal of Computational Physics 191 (2003) 392–419.
  * Subbareddy and Candler (2009) P. K. Subbareddy, G. V. Candler,  A fully discrete, kinetic energy consistent finite-volume scheme for compressible flows,  Journal of Computational Physics 228 (2009) 1347–1364.
  * Kok (2009) J. Kok,  A high-order low-dispersion symmetry-preserving finite-volume method for compressible flow on curvilinear grids,  Journal of Computational Physics 228 (2009) 6811–6832.
  * Morinishi (2010) Y. Morinishi,  Skew-symmetric form of convective terms and fully conservative finite difference schemes for variable density low-mach number flows,  Journal of Computational Physics 229 (2010) 276–300.
  * Pirozzoli (2010) S. Pirozzoli,  Generalized conservative approximations of split convective derivative operators,  Journal of Computational Physics 229 (2010) 7180–7190.
  * Pirozzoli (2011) S. Pirozzoli,  Stabilized non-dissipative approximations of euler equations in generalized curvilinear coordinates,  Journal of Computational Physics 230 (2011) 2997–3014.
  * Honein and Moin (2004) A. E. Honein, P. Moin,  Higher entropy conservation and numerical stability of compressible turbulence simulations,  Journal of Computational Physics 201 (2004) 531–545.
  * Tadmor (1987) E. Tadmor,  The numerical viscosity of entropy stable schemes for systems of conservation laws. i,  Mathematics of Computation 49 (1987) 91–103.
  * Tadmor (2003) E. Tadmor,  Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems,  Acta Numerica 12 (2003) 451–512.
  * Chandrashekar (2013) P. Chandrashekar,  Kinetic energy preserving and entropy stable finite volume schemes for compressible euler and navier-stokes equations,  Communications in Computational Physics 14 (2013) 1252–1286.
  * Coppola et al. (2019) G. Coppola, F. Capuano, S. Pirozzoli, L. de Luca,  Numerically stable formulations of convective terms for turbulent compressible flows,  Journal of Computational Physics 382 (2019) 86–104.
  * Kuya et al. (2018) Y. Kuya, K. Totani, S. Kawai,  Kinetic energy and entropy preserving schemes for compressible flows by split convective forms,  Journal of Computational Physics 375 (2018) 823–853.
  * Fuster (2013) D. Fuster,  An energy preserving formulation for the simulation of multiphase turbulent flows,  Journal of Computational Physics 235 (2013) 114–128.
  * Valle et al. (2020) N. Valle, F. X. Trias, J. Castro,  An energy-preserving level set method for multiphase flows,  Journal of Computational Physics 400 (2020) 108991.
  * Mirjalili and Mani (2021) S. Mirjalili, A. Mani,  Consistent, energy-conserving momentum transport for simulations of two-phase flows using the phase field equations,  Journal of Computational Physics (2021) 109918.
  * Jain (2022) S. S. Jain,  Accurate conservative phase-field method for simulation of two-phase flows,  arXiv preprint arXiv:2203.05802 (2022).
  * Jain and Moin (2020) S. S. Jain, P. Moin,  A kinetic energy and entropy preserving scheme for the simulation of compressible two-phase turbulent flows,  Center for Turbulence Research Annual Research Briefs (2020).
  * 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 (2020) 109606.
  * 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.
  * 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.
  * Tiwari et al. (2013) A. Tiwari, J. B. Freund, C. Pantano,  A diffuse interface model with immiscibility preservation,  Journal of Computational Physics 252 (2013) 290–309.
  * 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.
  * 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.
  * Chow and Moin (2003) F. K. Chow, P. Moin,  A further study of numerical errors in large-eddy simulations,  Journal of Computational Physics 184 (2003) 366–380.
  * Kennedy and Gruber (2008) C. A. Kennedy, A. Gruber,  Reduced aliasing formulations of the convective terms within the navier–stokes equations for a compressible fluid,  Journal of Computational Physics 227 (2008) 1676–1700.
  * Jameson (2008) A. Jameson,  Formulation of kinetic energy preserving conservative schemes for gas dynamics and direct numerical simulation of one-dimensional viscous compressible flow in a shock tube using entropy and kinetic energy preserving schemes,  Journal of Scientific Computing 34 (2008) 188–208.
  * Abgrall (1996) R. Abgrall,  How to prevent pressure oscillations in multicomponent flow calculations: a quasi conservative approach,  Journal of Computational Physics 125 (1996) 150–160.
  * Mahesh et al. (2004) K. Mahesh, G. Constantinescu, P. Moin,  A numerical method for large-eddy simulation in complex geometries,  Journal of Computational Physics 197 (2004) 215–240.
  * Hou and Mahesh (2005) Y. Hou, K. Mahesh,  A robust, colocated, implicit algorithm for direct numerical simulation of compressible, turbulent flows,  Journal of Computational Physics 205 (2005) 205–221.
  * Gassner et al. (2016) G. J. Gassner, A. R. Winters, D. A. Kopriva,  Split form nodal discontinuous galerkin schemes with summation-by-parts property for the compressible euler equations,  Journal of Computational Physics 327 (2016) 39–66.
  * 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.