# An interface capturing method for liquid-gas flows at low-Mach number

Federico Dalla Barba<EMAIL_ADDRESS>Nicoló Scapin
<EMAIL_ADDRESS>Marco E. Rosti<EMAIL_ADDRESS>Francesco Picano
<EMAIL_ADDRESS>Luca Brandt<EMAIL_ADDRESS>Department of
Industrial Engineering & CISAS, University of Padova, Padova, Italy  Complex
Fluids and Flows Unit, Okinawa Institute of Science and Technology Graduate
University (OIST), 1919-1 Tancha, Onna-son, Okinawa 904-0495, Japan
Department of Engineering Mechanics, Royal Institute of Technology (KTH),
Stockholm, Sweden  Department of Energy and Process Engineering, Norwegian
University of Science and Technology (NTNU), Trondheim, Norway

###### Abstract

Multiphase, compressible and viscous flows are of crucial importance in a wide
range of scientific and engineering problems. Despite the large effort paid in
the last decades to develop accurate and efficient numerical techniques to
address this kind of problems, current models need to be further improved to
address realistic applications. In this context, we propose a numerical
approach to the simulation of multiphase, viscous flows where a compressible
and an incompressible phase interact in the low-Mach number regime. In this
frame, acoustics is neglected but large density variations of the compressible
phase can be accounted for as well as heat transfer, convection and diffusion
processes. The problem is addressed in a fully Eulerian framework exploiting a
low-Mach number asymptotic expansion of the Navier-Stokes equations. A Volume
of Fluid approach (VOF) is used to capture the liquid-gas interface, built on
top of a massive parallel solver, second order accurate both in time and
space. The second-order-pressure term is treated implicitly and the resulting
pressure equation is solved with the eigenexpansion method employing a robust
and novel formulation. We provide a detailed and complete description of the
theoretical approach together with information about the numerical technique
and implementation details. Results of benchmarking tests are provided for
five different test cases.

###### keywords:

Compressible multi-phase flows, Volume-of-Fluid method, low-Mach number
asymptotic expansions, pressure-correction methods.

##  1 Introduction

Multiphase flows of two or more immiscible and viscous fluids are common in a
large variety of engineering applications and fundamental scientific problems.
Each phase is segregated and gives origin to complex and time-evolving free
boundaries where discontinuities in the flow fields exist [1]. The phases
mutually interact exchanging mass, momentum and energy across the free
boundaries, the latter undergoing large and complex deformations. It is
therefore clear how the description of the problem is extremely challenging
both from a theoretical and numerical point of view. Some of the major issues
affecting the modeling of multiphase flows arise from the discontinuities in
the flow variables and properties across the free boundaries, from the
necessity of numerically tracking and reconstructing the interfaces as well as
from the need to account for the effect of the surface tension and jump
conditions at the interfaces. These aspects are critical in the simulations of
incompressible and isothermal multiphase flows, but additional complexity is
add when the compressibility needs to be taken into account. In the latter
case, heat transfer between the fluid phases and the boundaries must be
considered in addition to mutual heat transfer between the phases and density
variations. In this paper, the emphasis is on the multiphase flows of two
immiscible viscous fluids, one of them being compressible, the other being
incompressible. The attention is paid in particular to the low-Mach number
flows, where the effect of compressibility is significant and large density
variations occur in the compressible phase while acoustics is negligible. This
is of great interest for several applications. To mention some, the simulation
of bubble-laden flows and boiling flows [2, 3, 4], as well as the simulation
of the fuel jet atomization processes in the combustion chambers of internal
combustion engines [5, 6].

A great amount of literature dealing with the numerical simulations of
multiphase, viscous fluids has been produced [7, 8, 9, 10, 11, 12, 13, 14].
Many approaches have been proposed both in the Eulerian and Lagrangian
frameworks as well as hybrid methods. Among the latter, the arbitrary
Lagrangian Eulerian (ALE) approach [15, 16, 17]. In this frame, interface-
conforming grids are used where boundary conditions can be accurately
prescribed on the free boundaries of the flow. The main advantage of the ALE
methods is the accurate treatment of the interfaces. Nevertheless, the
computational cost of this kind of simulations is large due to the adaptive
mesh adjustment needed to preserve the conformity of the grid to the time-
evolving interfaces. The need for re-meshing is removed in the frame of the
fixed-Eulerian-grid methods. These are hybrid Eulerian-Lagrangian approaches
also referred to as front-tracking methods [18]. The Immersed Boundary Method
(IBM) belongs to this class of numerical techniques [19, 20, 21]. This
approach consists in solving the governing equations for the flow on a fixed
Eulerian grid while tracking the free boundaries separating the different
phases of the flow by means of Lagrangian markers distributed over the
interfaces. An additional forcing is imposed to the fluid, within a
neighborhood of the interfaces, such that the boundary conditions are
satisfied within a certain degree of accuracy. This class of methods have been
successfully applied to the simulation of multi-phase flows [22, 23]. Even if
fixed-Eulerian-grid methods are computationally more efficient than
conforming-grid methods, they suffer from low accuracy in the reconstruction
and tracking of the interface. A popular alternative is the so-called front-
capturing methods that are based on a fully Eulerian treatment of the
interface tracking and reconstruction. These include essentially the Volume Of
Fluid (VOF) method and the Level-Set Method (LSM). The LSM uses a continuous
level-set function, usually the signed distance to the interface, to
distinguish between the different phases of the flow [24, 25, 26, 27, 28].
Interfaces are accurately defined by an assigned level of the level-set
function while the advection of the level-set function itself allows for an
accurate tracking of the interfaces. In the LSM framework the interface
curvature can be computed easily and accurately, nevertheless, these are not
mass-preserving methods. Indeed, the advection of the level-set function may
result in a mass loss or gain. On the contrary, the VOF approach, uses a
discontinuous colour function to represent each different phase [1, 29, 30,
31]. Interface reconstruction is achieved by a piece-wise representation of
the interface in each discrete cell of the computational grid. An accurate
computation of the numerical fluxes needed to advect the VOF function is
possible and numerical diffusion is avoided making VOF methods extremely
accurate from a mass-conservation point of view. A number of algorithms have
been developed to compute the exact interface surfaces from the discrete VOF
function. Among these are the Simple Line Interface Calculation (SLIC) method
[32], the Piecewise Linear Interface Calculation (PLIC) [33, 34, 35], the
Tangent of Hyperbola for Interface Capturing (THINC) [36] and its
multidimensional variant MTHINC [37].

A major part of the numerical approaches to the simulation of multiphase flows
reported by archival literature and referenced above were originally developed
for incompressible flows. A great effort has been spent in the last decades to
extend these methods to the simulation of compressible multiphase flows [38,
39, 40, 41, 42]; nonetheless this is still a very active area of research.
Among compressible flows, low-Mach number flows are of great interest for many
applications where large density variations occur at low speeds, low subsonic
regimes. When addressing the simulation of this flow regime, a big issue
arises from the limitation imposed on the time step by the fastest dynamics of
the flow. Indeed, in a compressible flow the speed of propagation of pressure
waves scales as $1/Ma$, $Ma$ being the Mach number. Many solutions to this
problem have been proposed, such as an implicit treatment of the acoustic
pressure [43, 44]. Nevertheless, if the case under examination is dominated by
free or forced convection where the amount of energy carried by the acoustics
is only a small fraction of the overall energy of the flow, a low-Mach number
asymptotic formulation of the Navier-Stokes equations can be used to
numerically address the problem. Large density variations can be accounted
for, completely neglecting acoustics, but still describing entropy and
vorticity modes as well as taking into account compressibility [45, 46]. In
this frame, the pressure is split into two different terms: a zero-order,
thermodynamic pressure, $p_{0}$, and a second-order pressure, $p_{2}$. The
former is governed by the thermodynamic properties of the flow while the
latter enters the computation in a similar fashion to that of pressure in
incompressible flows [45].

In this context, we propose a one-fluid fully Eulerian approach to the
numerical simulation of multiphase low-Mach number flows, based on the
solution of a low-Mach number asymptotic formulation of the compressible
Navier-Stokes equations. For the reconstruction and subsequent advection of
the interface between the compressible and incompressible phases, we adopt an
algebraic Volume-of-Fluid method (MTHINC [37]). However, the mathematical and
numerical framework can be extended in a straightforward manner to any kind of
interface capturing and tracking technique based on the sharp interface
approach. The proposed method is implemented in the frame of the pressure-
correction methods, taking advantage of a Fast-Fourier-Transform (FFT) based
solver for the Poisson equation governing the second-order-pressure of the
flow. The effect of the surface tension is accounted for by using the
continuum surface force (CSF) model by Brackbill [47]. The implementation is
built upon an extensively validated code for the simulation of incompressible
flows. The solver uses second order finite difference schemes for space
discretization on a fixed Eulerian grid and a second order of accuracy Adams-
Bashforth time-marching algorithm [48, 49, 12, 50]. We provide a detailed and
complete description of the theoretical approach together with information
about the numerical technique and implementation details. Because of its
numerical efficiency, we believe that this approach is one of the most
promising in order to efficiently address the simulation of multiphase, low-
Mach number flows, in particular when one of the two phases can be assumed to
be incompressible.

##  2 Governing equations

We consider the flow of two immiscible viscous fluids and let one of them be
compressible, the other being incompressible (e.g. a gas-liquid system).
Moreover, the compressible phase is assumed to evolve in the low-Mach number
regime. A phase indicator function $H$ is defined at position $\mathbf{x}$ and
time $t$ in order to distinguish between the two phases:

$H(\mathbf{x},t)=\begin{cases}1\hskip 14.22636pt\text{if $\mathbf{x}\in
V_{g}$}\mathrm{,}\\\ 0\hskip 14.22636pt\text{if $\mathbf{x}\in
V_{l}$}\mathrm{;}\end{cases}$ (1)

where $V_{g}$ and $V_{l}$ are the domains pertaining to the gas and liquid
phases, divided by a zero-thickness interface $S(t)=V_{g}\bigcap V_{l}$. The
dynamics of the liquid phase is governed by the standard incompressible
Navier-Stokes equations and details are not given here as the framework is
well established [29]. On the other hand, a more detailed analysis is
presented in the next section for the gas phase.

###  2.1 Asymptotic expansions for the compressible phase

Let now consider the compressible phase only and therefore, unless otherwise
stated, all the quantities defined in this section refer to it. In general, if
compressibility is taken into account, the gas phase can be described by the
following form of the Navier-Stokes and energy equations:

$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\mathbf{u})=0,$
(2) $\displaystyle\frac{\partial(\rho\mathbf{u})}{\partial
t}+\nabla\cdot(\rho\mathbf{u}\otimes\mathbf{u})=\nabla\cdot\boldsymbol{\tau}-\nabla{p}+\mathbf{f}_{\sigma}+\rho\mathbf{g},$
(3) $\displaystyle\frac{\partial(\rho e_{t})}{\partial
t}+\nabla\cdot(\rho\mathbf{u}e_{t})=\nabla\cdot(\boldsymbol{\tau}\cdot\mathbf{u})+\nabla\cdot(k\nabla
T)-\nabla\cdot({p\mathbf{u}})+\left(\mathbf{f}_{\sigma}+\rho\mathbf{g}\right)\cdot\mathbf{u},$
(4) $\displaystyle p=\rho\mathcal{R}T,$ (5)

where $\mathbf{u}=(u,v,w)$, $\rho$ and $p$ are the fluid velocity, density and
pressure, $\mathbf{g}$ is the gravitational acceleration and $k$ the thermal
conductivity. The specific total energy of the flow is denoted as $e_{t}$ and
includes the specific internal energy of the flow, $e$ and the specific
kinetic energy: $e_{t}=e+\mathbf{u}\cdot\mathbf{u}/2$. The Newton-Stokes
constitutive relation is assumed such that the viscous stress tensor is

$\boldsymbol{\tau}=\mu\left[(\nabla\mathbf{u}+\nabla\mathbf{u}^{T})-\frac{2}{3}(\nabla\cdot\mathbf{u})\mathbf{I}\right],$
(6)

with $\mathbf{I}$ the identity tensor and $\mu$ the dynamic viscosity. The
effect of the surface tension on the interfaces separating the incompressible
and compressible phases is modelled by a continuum surface force (CSF) [47]:

$\mathbf{f}_{\sigma}=\sigma\kappa\delta(\mathbf{x}-\mathbf{x}_{s})\mathbf{n},$
(7)

where $\sigma$ is the surface tension coefficient, $\kappa$ the curvature of
the interface and $\mathbf{n}$ the unit normal on the interface pointing
towards the compressible phase. A delta function,
$\delta(\mathbf{x}-\mathbf{x}_{s})$, is used in Eq. (7) to impose the force
density, ${\mathbf{f}}_{\sigma}$, only at the interface position
$\mathbf{x}_{s}$. The hypothesis of ideal gas is considered, such that the
constitutive relation (5) holds, where $\mathcal{R}=R/\mathcal{M}$, with
$\mathcal{M}$ the molar mass and $R=8.314$ $\mathrm{J/(mol\cdot K)}$ the
universal gas constant. The model could be easily improved by using a
different constitutive relation, such as Van der Waals equation or more
complex cubic equations of state [51]. Nonetheless, for the sake of simplicity
the discussion is limited to the use of Eq. (5). The hypothesis of calorically
perfect gas is assumed too, such that the constant-pressure and constant-
volume heat capacities, $c_{p}$ and $c_{v}$, do not depend on the
thermodynamic pressure and temperature. Under this hypothesis the energy Eq.
(4) can be formulated in terms of the sensible internal energy and enthalpy
only:

$\displaystyle e=\Delta h_{T_{ref}}^{0}+c_{v}(T-T_{ref})=\Delta
h_{T_{ref}}^{0}+c_{v}T,$ (8) $\displaystyle h=\Delta
h_{T_{ref}}^{0}+c_{p}(T-T_{ref})=\Delta h_{T_{ref}}^{0}+c_{p}T,$ (9)

with $h=e+p/\rho$ the fluid enthalpy, $T$ the temperature and $\Delta
h_{T_{ref}}^{0}$ the fluid enthalpy of formation evaluated at $T_{ref}=0$
$\mathrm{K}$. Eq. (2)-(5) can now be recast in non-dimensional form setting as
independent reference scales the density, $\tilde{\rho}$, the pressure,
$\tilde{p}$, the length, $\tilde{L}$, the velocity, $\tilde{U}$ as well as the
following derived quantities:

$\tilde{T}=\tilde{p}/(\mathcal{R}\tilde{\rho}),\
\tilde{t}=\tilde{L}/\tilde{U},\ \tilde{e}=\tilde{p}/\tilde{\rho},\
\tilde{f}=\tilde{\rho}\tilde{U}^{2}/\tilde{L}.$

The scales $\tilde{T}$, $\tilde{t}$, $\tilde{e}$ and $\tilde{f}$ are the
reference temperature, time, specific energy and force per unit volume,
respectively. In addition, reference values for the thermal diffusion
coefficient, heat capacities, dynamic viscosity, surface tension coefficient
of the fluid and gravitational acceleration are $\tilde{k},\ \tilde{c}_{p},\
\tilde{c}_{v},\ \tilde{\mu}$, $\tilde{\sigma}$ and $\tilde{g}$. After some
manipulation, the non-dimensional form of Eq. (2) - (5) reads:

$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\mathbf{u})=0,$
(10) $\displaystyle\frac{\partial(\rho\mathbf{u})}{\partial
t}+\nabla\cdot(\rho\mathbf{u}\otimes\mathbf{u})=\frac{1}{Re}\nabla\cdot\boldsymbol{\tau}-\frac{1}{Ma^{2}}\nabla{p}+\frac{\mathbf{f}_{\sigma}}{We}+\frac{\rho\mathbf{g}}{Fr^{2}},$
(11) $\displaystyle\frac{\partial(\rho e)}{\partial
t}+\nabla\cdot(\rho\mathbf{u}h)+Ma^{2}\left[\frac{\partial}{\partial
t}\left(\rho\frac{\mathbf{u}\cdot\mathbf{u}}{2}\right)+\nabla\cdot\left(\rho\frac{\mathbf{u}\cdot\mathbf{u}}{2}\right)\right]$
(12) $\displaystyle=\dfrac{\gamma}{\gamma-1}\frac{1}{RePr}\nabla\cdot(k\nabla
T)+Ma^{2}\left[\dfrac{1}{Re}\nabla\cdot(\boldsymbol{\tau}\cdot\mathbf{u})+\dfrac{\mathbf{f}_{\sigma}}{We}+\frac{\rho\mathbf{g}}{Fr^{2}}\right]\cdot\mathbf{u},$
$\displaystyle p=\rho T,$ (13)

where $Ma=\tilde{U}/\sqrt{\tilde{p}/\tilde{\rho}}$ is a pseudo-Mach number,
$Re=\tilde{\rho}\tilde{U}\tilde{L}/\tilde{\mu}$ the Reynolds number,
$Pr=\tilde{c}_{p}\tilde{\mu}/\tilde{k}$ the Prandtl number,
$We=\tilde{\rho}\tilde{U}^{2}\tilde{L}/\tilde{\sigma}$ the Weber number and
$Fr=\tilde{U}/\sqrt{g\tilde{L}}$ the Froude number, while
$\gamma=\tilde{c}_{p}/\tilde{c}_{v}$. The low-Mach number limit of Eq. (10) -
(13) can be derived from a single-scale asymptotic perturbation for small Mach
numbers [46]. The pseudo Mach number, $Ma$, is present in all the equations
only with the power of two, hence each quantity can be expanded in the
following way:

$\displaystyle\mathbf{f}(\mathbf{x},t)=\mathbf{f}_{0}(\mathbf{x},t)+\mathbf{f}_{2}(\mathbf{x},t)Ma^{2}+O(Ma^{3}),$
(14)

$f$ being a generic vectorial or scalar quantity. It is also possible to proof
that the following hold for the product of two scalar quantities:

$\displaystyle\left[f(\mathbf{x},t)g(\mathbf{x},t)\right]_{0}=f_{0}(\mathbf{x},t)g_{0}(\mathbf{x},t),$
(15)
$\displaystyle\left[f(\mathbf{x},t)g(\mathbf{x},t)\right]_{2}=f_{2}(\mathbf{x},t)g_{0}(\mathbf{x},t)+f_{0}(\mathbf{x},t)g_{2}(\mathbf{x},t).$
(16)

In the following, we report the procedure to obtain the low-Mach number limit
of the momentum Eq. (11), the algebraic manipulation being similar for the
remaining ones. Applying (14) to (11) leads to:

$\displaystyle\frac{\partial}{\partial
t}\left[(\rho\mathbf{u})_{0}+(\rho\mathbf{u})_{2}Ma^{2}+O(Ma^{3})\right]+$
(17)
$\displaystyle+\nabla\cdot\left[(\rho\mathbf{u}\otimes\mathbf{u})_{0}+(\rho\mathbf{u}\otimes\mathbf{u})_{2}Ma^{2}+O(Ma^{3})\right]=$
$\displaystyle=\frac{1}{Re}\nabla\cdot\left[{\boldsymbol{\tau}}_{0}+{\boldsymbol{\tau}}_{2}Ma^{2}+O(Ma^{3})\right]-\frac{1}{Ma^{2}}\nabla\left[p_{0}+p_{2}Ma^{2}+O(Ma^{3})\right]+$
$\displaystyle+\frac{1}{We}\left[\mathbf{f}_{\sigma 0}+\mathbf{f}_{\sigma
0}Ma^{2}+O(Ma^{3})\right]+\frac{1}{Fr}\left[(\rho\mathbf{g})_{0}+(\rho\mathbf{g})_{2}Ma^{2}+O(Ma^{3})\right].$

The terms in Eq. (17) with the same order in $Ma$ can be grouped. After some
manipulation, we obtain for the zeroth-order pressure:

$\displaystyle\nabla p_{0}=0,$ (18)

and for the second order correction:

$\displaystyle\frac{\partial\mathbf{u}_{0}}{\partial
t}+\mathbf{u}_{0}\cdot\nabla\mathbf{u}_{0}=\frac{1}{\rho_{0}}\left[\frac{1}{Re}\nabla\cdot\boldsymbol{\tau}_{0}-\nabla
p_{2}+\frac{\mathbf{f}_{\sigma 0}}{We}\right]+\frac{\mathbf{g}}{Fr^{2}}.$ (19)

Algebraic manipulation is completely omitted for the continuity and energy
equations. The reader is referred to existing references for additional
details [45, 46, 52]. The final low-Mach number relations are directly given
below:

$\displaystyle\frac{\partial\rho_{0}}{\partial
t}+\nabla\cdot(\rho_{0}\mathbf{u}_{0})=0,$ (20)
$\displaystyle\frac{\partial\mathbf{u}_{0}}{\partial
t}+\mathbf{u}_{0}\cdot\nabla\mathbf{u}_{0}=\frac{1}{\rho_{0}}\left[\frac{1}{Re}\nabla\cdot\boldsymbol{\tau}_{0}-\nabla
p_{2}+\frac{\mathbf{f}_{\sigma 0}}{We}\right]+\frac{\mathbf{g}}{Fr^{2}},$ (21)
$\displaystyle\frac{\partial(\rho_{0}e_{0})}{\partial
t}+\nabla\cdot(\rho_{0}\mathbf{u}_{0}h_{0})=\frac{\gamma}{\gamma-1}\frac{1}{RePr}\nabla\cdot(k\nabla
T_{0}),$ (22) $\displaystyle p_{0}=\rho_{0}T_{0}.$ (23)

As mentioned above, two different pressure terms appear in these equations:
the zeroth-order pressure $p_{0}$ and the second-order pressure $p_{2}$. The
former, which can be referred to as thermodynamic pressure, is determined by
the thermodynamic state of the flow, it is uniform across the spatial field
and it is a function of the time only. The latter, conversely, enters the
computation in a similar fashion to that of pressure in incompressible flows
(e.g. by imposing a prescribed value of the velocity divergence) and it is
obtained as solution of a Poisson equation that will be discussed in the
following. It is useful to remind that, for the chosen set of reference
scales, the non-dimensional sensible enthalpy reads
$h_{0}=\gamma/(\gamma-1)T_{0}$ while the sensible internal energy
$e_{0}=1/(\gamma-1)T_{0}$. From here on, the subscript referring to the order
of quantities will be omitted except for the pressure terms. After some
additional manipulation Eq. (20) - (23) can be recast as:

$\displaystyle\frac{\partial\mathbf{u}}{\partial
t}+\mathbf{u}\cdot\nabla\mathbf{u}=\frac{1}{\rho}\left[\frac{1}{Re}\nabla\cdot\boldsymbol{\tau}-\nabla
p_{2}+\frac{\mathbf{f}_{\sigma}}{We}\right]+\frac{\mathbf{g}}{Fr^{2}},$ (24)
$\displaystyle\frac{\partial T}{\partial t}+\mathbf{u}\cdot\nabla
T=\frac{T}{p_{0}}\left[\frac{1}{RePr}\nabla\cdot(k\nabla
T)+\frac{\gamma-1}{\gamma}\frac{dp_{0}}{dt}\right],$ (25)
$\displaystyle\nabla\cdot\mathbf{u}=\frac{1}{p_{0}}\left[\frac{1}{RePr}\nabla\cdot(k\nabla
T)+\frac{1}{\gamma}\frac{dp_{0}}{dt}\right],$ (26)
$\displaystyle\frac{dp_{0}}{dt}=\frac{\gamma}{V}\left[\frac{1}{RePr}\int_{S}k\nabla
T\cdot\mathbf{n}\ dS-p_{0}\int_{S}\mathbf{u}\cdot\mathbf{n}\ dS\right],$ (27)
$\displaystyle p_{0}=\rho T.$ (28)

###  2.2 Final form of the governing equations

Eq. (24), (25), (26), (27) and (28) are valid only for the compressible gas
phase and can be coupled with those for the incompressible phase using the
phase indicator function defined by Eq. (1). By doing so, a monolithic system
of equations can be obtained describing at the same time the dynamics of both
the compressible and incompressible phases. For convenience the system is
written in dimensional form,

$\displaystyle\frac{\partial\mathbf{u}}{\partial
t}+\mathbf{u}\cdot\nabla\mathbf{u}=\frac{1}{\rho}\left(\nabla\cdot\boldsymbol{\tau}-\nabla
p_{2}+\mathbf{f}_{\sigma}\right)+\mathbf{g},$ (29)
$\displaystyle\frac{\partial T}{\partial t}+\mathbf{u}\cdot\nabla
T=\frac{1}{\rho c_{p}}\left[\nabla\cdot(k\nabla T)+\frac{dp_{0}}{dt}H\right],$
(30)
$\displaystyle\nabla\cdot\mathbf{u}=\frac{1}{p_{0}}\left[\frac{\gamma_{g}-1}{\gamma_{g}}\nabla\cdot(k\nabla
T)-\frac{1}{\gamma_{g}}\frac{dp_{0}}{dt}\right]H,$ (31)
$\displaystyle\frac{dp_{0}}{dt}=\frac{\gamma_{g}}{V_{g}}\left(\frac{\gamma_{g}-1}{\gamma_{g}}\int_{S}k\nabla
T\cdot\mathbf{n}\ dS-p_{0}\int_{S}\mathbf{u}\cdot\mathbf{n}\ dS\right),$ (32)
$\displaystyle\rho_{g}=\dfrac{p_{0}}{\mathcal{R}T}.$ (33)

The variables $c_{p}$, $k$ and $\mu$ are the heat capacity at constant
pressure, thermal conductivity and dynamic viscosity of the mixture,
respectively. These quantities are averaged, together with the mixture
density, using the phase indicator:

$\displaystyle\rho=\rho_{g}H+\rho_{l}(1-H),$ (34) $\displaystyle
c_{p}=c_{p,g}H+c_{p,l}(1-H),$ (35) $\displaystyle k=k_{g}H+k_{l}(1-H),$ (36)
$\displaystyle\mu=\mu_{g}H+\mu_{l}(1-H),$ (37)

where $c_{p}$ refers to the specific heat capacity at constant pressure and
the subscripts “$g$” and “$l$” to the physical parameter of the gas and liquid
phases. The ratio $\gamma_{g}$ depends only on the local thermodynamic
variables, but in the present work is taken as a constant and equal to
$c_{p,g}/c_{v,g}$.

Finally, it is worth remarking that while Eq. (29) - (31) hold for both the
phases, Eq. (32) - (33) are meaningful for the compressible phase only. In
particular, Eq. (31) reduces to $\nabla\cdot\mathbf{u}=0$ in the liquid
regions, where $H=0$. In the monolithic system of equations above, the zeroth-
order pressure, $p_{0}$, results to be constant over each of the closed
regions in the computational domain filled by the compressible phase. Its
value is set by the constitutive law for ideal gases (33) while its temporal
rate of change is set by the energy balance provided by Eq. (32). In
particular, the surface integrals appearing in Eq. (32) have to be intended as
integral computed over both the frontier defined by the boundaries of the
computational domain and the interface $S$ separating the compressible and
incompressible regions. Nevertheless, for its numerical integration, it is
more convenient to reformulate Eq. (32) in terms of volume integral by means
of the divergence theorem:

$\displaystyle\frac{dp_{0}}{dt}=\frac{\gamma_{g}}{V_{g}}\left(\frac{\gamma_{g}-1}{\gamma_{g}}\int_{V}H\nabla\cdot(k\nabla
T)\ dV-p_{0}\int_{V}H\nabla\cdot\mathbf{u}\ dV\right),$ (38)

where the volume over which the integrals are computed is the total volume
$V=V_{g}\bigcup V_{l}$. Eq. (38) is derived by employing the change of
variable $V_{g}=HV$, with the differential $dV_{g}=HdV$, being $dH=0$ by
definition of the phase indicator function [29].

##  3 Numerical methodology

The numerical solution of Eq. (29)-(33) is addressed on a fixed regular
Cartesian grid (e.g. using a uniform and equal spacing, $\Delta x=\Delta
y=\Delta z$), with a marker-and-cell arrangement of velocity and pressure
points, whereas all scalar fields are defined at the cell centers. Hereafter,
we present the numerical discretization of the governing equations following
the same order in which they are solved.

###  3.1 Interface representation and advection

The first step of each iteration of the time-marching algorithm consists in
the reconstruction of the interface between the two phases and its subsequent
advection. As mentioned in the introductory section, we address both the
aspects in a fully Eulerian framework using the VOF method to distinguish
between each of the flow phases [1]. By a numerical point of view the
indicator function $H$, defined in Eq. (1), is updated on the computational
grid by the following advection equation:

$\frac{\partial\Phi}{\partial
t}+\nabla\cdot(\mathbf{u}H)=\Phi\nabla\cdot\mathbf{u},$ (39)

where the volume fraction, $\Phi$, is defined as the average value of the
color function over a discrete computational cell of volume $\Delta V=\Delta
x\Delta y\Delta z$:

$\Phi=\int_{\Delta V}H(\mathbf{x},t)dV.$ (40)

Coherently with the given definition of $H$ in Eq. (1), the volume fraction
satisfies $\Phi=1$ in cells occupied by the gas phase only, $\Phi=0$ in that
filled only by the liquid and $0<\Phi<1$ in the cells containing the liquid-
gas interface.

In the present work, we employ as VOF method the multi-dimensional tangent of
hyperbola for interface capturing method (MTHINC), originally developed by Ii
et al. [37] and more recently applied to complex flows cases both in laminar
[12, 53, 54] and in turbulent conditions [50]. The description of the
procedure is not reported here as our procedure is not limited to a specific
interface capturing/tracking method; additional details are provided in the
references [37, 12]. Once the interface is reconstructed, the advection step
is performed using the standard directional splitting approach, originally
developed by Puckett et al. [55] and Aulisa et al. [35]. An additional
correction is used in order to account for the non-zero divergence of the
interface velocity [56], that in the present framework coincides with the one-
fluid velocity $\mathbf{u}$. Finally, the thermodynamic properties independent
of temperature ($\mu$, $k$ and $c_{p}$) are updated using the relations (35),
(36) and (37).

###  3.2 Thermodynamic pressure and temperature equation

The next step consists in the computation of the updated thermodynamic
pressure $p_{0}^{n+1}$ and temperature $T^{n+1}$ using a second order Adams-
Bashforth time-marching algorithm:

$\displaystyle p_{0}^{n+1}=p_{0}^{n}+\Delta
t^{n+1}\left[\left(1+\dfrac{1}{2}\dfrac{\Delta t^{n+1}}{\Delta
t^{n}}\right)RP^{n}-\left(\dfrac{1}{2}\dfrac{\Delta t^{n+1}}{\Delta
t^{n}}\right)RP^{n-1}\right],$ (41) $\displaystyle T^{n+1}=T^{n}+\Delta
t^{n+1}\left[\left(1+\dfrac{1}{2}\dfrac{\Delta t^{n+1}}{\Delta
t^{n}}\right)RT^{n}-\left(\dfrac{1}{2}\dfrac{\Delta t^{n+1}}{\Delta
t^{n}}\right)RT^{n-1}\right].$ (42)

The variables $\Delta t^{n+1}$ and $\Delta t^{n}$ represent the time step
evaluated at the time levels $n+1$ and $n$, respectively. The time step is
chosen to fulfill the temporal stability requirements as explained in section
3.4. The terms $RP$ and $RT$ are the right-hand side of the thermodynamic
pressure and temperature, provided below in a semi-discrete notation:

$\displaystyle
RP^{n}=\frac{\gamma_{g}}{V_{g}^{n}}\left[\frac{\gamma_{g}-1}{\gamma_{g}}\int_{V^{n}}\Phi^{n}\nabla\cdot(k^{n}\nabla
T^{n})\ dV-p_{0}^{n}\int_{V^{n}}\Phi^{n}\nabla\cdot\mathbf{u}^{n}\ dV\right],$
(43) $\displaystyle RT^{n}=-\mathbf{u}^{n}\cdot\nabla
T^{n}+\frac{1}{\rho^{n}c_{p}^{n+1}}\left[\nabla\cdot(k^{n+1}\nabla
T^{n})+\left(\dfrac{dp_{0}}{dt}\right)^{n}\Phi^{n}\right],$ (44)

where the rate of change of the thermodynamic pressure is computed as
$(dp_{0}/dt)^{n}=RP^{n}$. The gas volume $V_{g}^{n}$ over which Eq. (41) is
integrated can be approximated as:

$V_{g}^{n}=\sum_{i=1}^{N_{x}}\sum_{j=1}^{N_{y}}\sum_{k=1}^{N_{z}}\Phi_{i,j,k}^{n}\Delta
x\Delta y\Delta z,$ (45)

$N_{x}$, $N_{y}$ and $N_{z}$ being the number of grid points along the $x$,
$y$ and $z$ directions, respectively. All the spatial terms in Eq. (43) and
(44) are discretized by second order central schemes, except for the
temperature convection. The discretization of this latter is based on the
$5^{th}$-order WENO5 scheme described in reference [57]. Once $T^{n+1}$ and
$p_{0}^{n+1}$ are known, the volumetric density field is updated using Eq.
(37) while the gas density is computed using the equation of state, e.g.
$\rho_{g}^{n+1}=p_{0}^{n+1}/\mathcal{R}T^{n+1}$. It should be noted that,
following this approach, the right-hand side the energy equation, $RT^{n+1}$,
contains the specific heat capacity, $c_{p}$, evaluated at the new time level,
$n+1$, while the density, $\rho$, is only available at the previous time
level, $n$. Despite this inconsistency, dictated by the coupling between the
density and temperature field, the approach adopted in the present work
provides convergent and stable results in all the validation cases discussed
in the following. This inconsistency could be fixed by computing $\rho^{n+1}$
via $\Phi^{n+1}$, but employing the old value of $p_{0}^{n}$ and $T^{n}$. Even
if we do not provide the numerical results, we tested this alternative
solution observing that the numerical results of the two approaches were
indistinguishable.

###  3.3 Flow solver

In order to impose that the velocity field $\mathbf{u}^{n+1}$ satisfies the
divergence constraint given by Eq. (50), a pressure-correction scheme based on
the Adams-Bashforth method is employed and summarized below in semi-discrete
notation:

$\displaystyle\mathbf{u}^{*}=\mathbf{u}^{n}+\Delta
t^{n+1}\left[\left(1+\dfrac{1}{2}\dfrac{\Delta t^{n+1}}{\Delta
t^{n}}\right)\mathbf{RU}^{n}-\left(\dfrac{1}{2}\dfrac{\Delta t^{n+1}}{\Delta
t^{n}}\right)\mathbf{RU}^{n-1}\right],$ (46)
$\displaystyle\nabla\cdot\left(\dfrac{1}{\rho^{n+1}}\nabla
p_{2}^{n+1}\right)=\dfrac{1}{\Delta
t^{n+1}}\left[\nabla\cdot\mathbf{u}^{*}-\nabla\cdot\mathbf{u}^{n+1}\right],$
(47) $\displaystyle\mathbf{u}^{n+1}=\mathbf{u}^{*}-\dfrac{\Delta
t^{n+1}}{\rho^{n+1}}\nabla p_{2}^{n+1},$ (48)

where $\mathbf{u}^{*}$ is the predicted velocity. The right-hand side
$\mathbf{RU}^{n}$ is computed as:

$\displaystyle\mathbf{RU}^{n}=-\mathbf{u}^{n}\cdot\nabla\mathbf{u}^{n}+\dfrac{1}{\rho^{n+1}}\left[\nabla\cdot\boldsymbol{\tau}(\mu^{n+1},\mathbf{u}^{n})+\mathbf{f}_{\sigma}^{n+1}+\rho^{n+1}\mathbf{g}\right]\mathrm{,}$
(49)

where both the convection and diffusion terms are discretized by central
schemes. More specifically, the former is discretized in divergence form as
$\nabla\cdot(\mathbf{u}\mathbf{u})-\mathbf{u}\nabla\cdot\mathbf{u}$ whereas
the latter is treated in a fully conservative form. The color function
gradients, $\nabla\Phi^{n+1}$, in the expression for the surface tension
force, $\mathbf{f}_{\sigma}$, are estimated with the Youngs’ method [33, 34].
Finally, the divergence constraint in Eq. (47) is computed directly as

$\displaystyle\nabla\cdot\mathbf{u}^{n+1}=\frac{1}{p_{0}^{n+1}}\left[\frac{\gamma_{g}-1}{\gamma_{g}}\nabla\cdot(k^{n+1}\nabla
T^{n+1})-\frac{1}{\gamma_{g}}\left(\dfrac{dp_{0}}{dt}\right)^{n+1}\right]\Phi^{n+1}\mathrm{.}$
(50)

It should be remarked, that the evaluation of the term $(dp_{0}/dt)^{n+1}$ in
Eq. (50) implies the use of $\mathbf{u}^{n+1}$, which is unknown. A possible
strategy is to introduce inside each time-step a further inner loop to couple
$p_{0}^{n+1}$ and $\mathbf{u}^{n+1}$ [58]. In the present work, we adopt a
simpler but equally robust solution based on an extrapolation in time of the
term $(dp_{0}/dt)^{n+1}$:

$\left(\dfrac{dp_{0}}{dt}\right)^{n+1}\simeq\left(1+\dfrac{\Delta
t^{n+1}}{\Delta
t^{n}}\right)\left(\dfrac{dp_{0}}{dt}\right)^{n}-\left(\dfrac{\Delta
t^{n+1}}{\Delta t^{n}}\right)\left(\dfrac{dp_{0}}{dt}\right)^{n-1},$ (51)

where $(dp_{0}/dt)^{n}=RP^{n}$ and $(dp_{0}/dt)^{n-1}=RP^{n-1}$ computed using
Eq. (43). The resulting system of equations (46), (47) and (48), is formally
identical to that used by any incompressible two-fluid solver where the one-
fluid velocity field is not divergence-free. In general, this may occur in
presence of phase change at the interface (see reference [3]), but in our
mathematical framework it can occur only in the regions occupied by the
compressible phase.

####  3.3.1 Pressure equation

A key feature of any two-fluid solver is the ability to impose accurately and
efficiently the divergence constraint on the velocity field, this task being
directly related to the numerical procedure used to solve the pressure
equation, Eq. (47). A possible approach is based on the use of iterative
multigrid solvers. Despite the success and the widespread use of these
solvers, the solution is not exact, but satisfied up to a controlled
tolerance, usually of the order $\varepsilon=10^{-7}$ \- $10^{-8}$. Moreover,
since the coefficients of the Poisson equation vary in space, the system
matrix must be recomputed at each time-step. Alternatively, when the pressure
boundary conditions are homogeneous [59], a possible solution is to transform
Eq. (47) into a constant coefficient equation and apply the method of the
eigenexpansion [60, 61] to solve the pressure equation exactly with spectral
accuracy. The resulting pressure equation is still to be solved in an
iterative manner starting with an initial guess. Different methods based on
the latter approach are available in literature, the main difference being how
the variable coefficient pressure equation is recast into a constant
coefficient problem. In the following, we review two of these methods, that
have been recently proposed and designed to efficiently solve in an iterative
manner the Poisson equation with the method of eigenexpansion. It should be
noted that, these methods have been successfully implemented in numerical
codes that share with our one a similar parallelization strategy based on the
2DECOMP&FFT [62] library. Finally, we will introduce a new methodology that
proves to be more efficient and suitable for two-phase flows with capillary
effects and sharp gradients between the two phases.

  * 1.

Method I: this method has been proposed by Motheau and Abraham [63]. It is
designed for low-Mach number reactive flows, aiming at decreasing the number
of iterations of the previously developed FFT-based solvers for combustion
applications. The methodology consists in a semi-implicit approach that first
requires an iterative procedure to solve the following constant coefficient
Poisson equation:

$\nabla^{2}p_{2}^{s+1}=\nabla\cdot\left[\left(1-\dfrac{\tilde{\rho}_{0}^{n+1}}{\rho^{n+1}}\right)\nabla
p_{2}^{s}\right]+\dfrac{\tilde{\rho}_{0}^{n+1}}{\Delta
t^{n+1}}\left(\nabla\cdot\mathbf{u}^{*}-\nabla\cdot\mathbf{u}^{n+1}\right)\mathrm{,}$
(52)

where, $p_{2}^{s+1}$ and $p_{2}^{s}$ are the hydrodynamic pressure at two
subsequent iterations and $\tilde{\rho}_{0}^{n+1}$ is the minimum value of
$\rho^{n+1}$ over the computational domain. After the iterative loop to
compute $p_{2}^{n+1}$, a modified correction step is applied:

$\mathbf{u}^{n+1}=\mathbf{u}^{*}-\Delta
t^{n+1}\left[\dfrac{1}{\tilde{\rho}_{0}^{n+1}}\nabla
p_{2}^{n+1}+\left(\dfrac{1}{\rho^{n+1}}-\dfrac{1}{\tilde{\rho}_{0}^{n+1}}\right)\nabla
p_{2}^{n+1,q}\right]\mathrm{,}$ (53)

where $p_{2}^{n+1,q}$ is the second-to-last hydrodynamic pressure of the
iterative procedure at the new time-level. The main advantage of this method
is the ability to effectively impose the velocity divergence up to machine
accuracy, by setting a residual threshold to solve Eq. (52) to
$\varepsilon_{t}=10^{-6}-10^{-8}$, being the residual
$\varepsilon=||p_{2}^{s+1}-p_{2}^{s}||$). Nevertheless, we find that using
this approach in the case of a rising bubble (see section 4.2), the number of
iterations required to achieve convergence is of the order of one hundred.

  * 2.

Method II: this approach, proposed by Bartholomew and Laizet [64] and designed
for non-Boussinesq gravity currents, is based on the rearrangement of Eq. (47)
as a constant coefficient Poisson equation:

$\nabla^{2}p_{2}^{s+1}=\nabla^{2}p_{2}^{s}+\tilde{\rho}\left[\dfrac{1}{\Delta
t^{n+1}}\left(\nabla\cdot\mathbf{u}^{*}-\nabla\cdot\mathbf{u}^{n+1}-\nabla\cdot\dfrac{1}{\rho^{n+1}}\nabla
p_{2}^{s}\right)\right]\mathrm{,}$ (54)

where $p_{2}^{s+1}$ and $p_{2}^{s}$ are the hydrodynamic pressure at two
subsequent iterations. Eq. (54) is solved in an iterative manner until
convergence. After that, the correction step (48) is performed to obtain the
new velocity field $\mathbf{u}^{n+1}$. The modified density $\tilde{\rho}$ is
taken as the harmonic mean between $\rho_{l}$ and $\rho_{g}$, as suggested by
the authors to improve convergence. To satisfy the divergence constraint up to
machine accuracy, the threshold residual should be set to
$\varepsilon_{t}=10^{-12}$. In section (4.2) we will show that this approach
requires a lower number of iterations than the previous one, but still higher
than the one we are going to present next.

  * 3.

Method III: the basic idea behind this third approach, proposed here, is to
rearrange the Poisson equation into a constant-coefficient form by employing
the correction step (48):

$\begin{split}\dfrac{1}{\rho^{n+1}}\nabla^{2}p_{2}^{n+1}+\nabla\left(\dfrac{1}{\rho^{n+1}}\right)\cdot\nabla
p_{2}^{n+1}&=\dfrac{1}{\Delta
t^{n+1}}\left(\nabla\cdot\mathbf{u}^{*}-\nabla\cdot\mathbf{u}^{n+1}\right)\mathrm{,}\\\
\nabla^{2}p_{2}^{n+1}-\dfrac{1}{\rho^{n+1}}\nabla\rho^{n+1}\cdot\nabla
p_{2}^{n+1}&=\dfrac{\rho^{n+1}}{\Delta
t^{n+1}}\left(\nabla\cdot\mathbf{u}^{*}-\nabla\cdot\mathbf{u}^{n+1}\right)\mathrm{.}\end{split}$
(55)

Using the vector calculus identity
$\rho\nabla\cdot\mathbf{u}=\nabla\cdot(\rho\mathbf{u})-\mathbf{u}\cdot\nabla\rho$,
we finally rewrite Eq. (55) as:

$\displaystyle\nabla^{2}p_{2}^{n+1}=\dfrac{1}{\Delta
t^{n+1}}\left[\nabla\cdot(\rho^{n+1}\mathbf{u}^{*})-\rho^{n+1}\nabla\cdot\mathbf{u}^{n+1}-\mathbf{u}^{n+1}\cdot\nabla\rho^{n+1}\right]$
(56) $\displaystyle+\dfrac{1}{\Delta
t^{n+1}}\underbrace{\left[\left(\mathbf{u}^{n+1}-\mathbf{u}^{*}+\dfrac{\Delta
t^{n+1}}{\rho^{n+1}}\nabla
p_{2}^{n+1}\right)\cdot\nabla\rho^{n+1}\right]}_{\text{$=0$ due to
Eq.~{}\eqref{eqn:corr}}}\mathrm{.}$

As $p_{2}^{n+1}$ and $\mathbf{u}^{n+1}$ are both unknown, we solve Eq. (56)
together with Eq. (48) by an iterative loop, as reported in the pseudocode 1.
Two interesting features emerge when using this method. First, the constant-
coefficient Poisson Eq. (56) is an equivalent and exact formulation of its
variable counterpart (47), derived using the correction step (48). This
represents a major difference with respect to the previous two methods, Eq.
(52) and (54), which are only a consistent but not exact recast of Eq. (47).
Second, this method allows us to define and control the residual of the
iterative procedure on the basis of the velocity divergence, which represents
the constraint to be imposed on the flow field. These two advantages come at
the cost of performing a correction step for each iteration of the loop.
Nevertheless, the additional computational cost is more than compensated by
the lower numbers of iterations required to achieve convergence as the
solution of the Poisson equation is often the most expensive part in standard
two-fluid solvers. As we will show in the result section, this approach
requires a significantly lower number of iterations.

Algorithm 1 Solution of the pressure equation with Method III

1: $s=0,$

2: $\mathbf{u}^{s}=\mathbf{u}^{*},$

3: $\varepsilon=a\varepsilon_{t}$ with $a>1.$

4: while $\varepsilon>\varepsilon_{t}$ do

5: $s=s+1\mathrm{,}$

6: $\nabla^{2}p_{2}^{s+1}=\frac{1}{\Delta
t^{n+1}}\left[\nabla\cdot\left(\rho^{n+1}\mathbf{u}^{*}\right)-\rho^{n+1}\nabla\cdot\mathbf{u}^{n+1}-\mathbf{u}^{s}\cdot\nabla\rho^{n+1}\right]\mathrm{,}$

7: $\mathbf{u}^{s+1}=\mathbf{u}^{*}-\frac{\Delta t^{n+1}}{\rho^{n+1}}\nabla
p_{2}^{s+1}\mathrm{,}$

8:
$\varepsilon=\displaystyle{\sum_{i=1}^{N_{x}}\sum_{j=1}^{N_{y}}\sum_{k=1}^{N_{z}}}||\nabla\cdot\mathbf{u}_{i,j,k}^{s+1}-\nabla\cdot\mathbf{u}_{i,j,k}^{n+1}||\mathrm{.}$

9: end while

10: $p_{2}^{n+1}=p_{2}^{s+1},$

11: $\mathbf{u}^{n+1}=\mathbf{u}^{s+1}.$

The proposed methodology is outlined in the pseudocode 1 where the iterations
are performed until $\varepsilon\leq\varepsilon_{t}$, is satisfied. The
residual $\varepsilon$ is computed with a summation over the whole
computational domain whereas the threshold value $\varepsilon_{t}$ is chosen
considering the trade-off between the number of iterations required to achieve
convergence and the minimization of the residual error. Unless otherwise
stated, in the present work, $\varepsilon_{t}$ equal to $10^{-14}$ has been
set in order to impose the divergence constraint on the final velocity field,
$\mathbf{u}^{n+1}$, with machine accuracy. Finally, it should be noted that
the first two terms of Eq. (56) do not vary over the solution cycle and can be
pre-computed just before the iterative procedure to reduce the execution time,
whereas the last term, $\mathbf{u}^{s}\cdot\nabla\rho^{n+1}$ needs to be
updated at every iteration of the pressure-correction loop.

Before proceeding to the discussion of the validation cases, it is worth
remarking that the proposed mathematical and numerical framework can be
naturally extended to any interface capturing and tracking method consistent
with the sharp-interface definition of the phase indicator function as in Eq.
(1), e.g. Volume-of-Fluid, Level-set and Front-Tracking methods. Furthermore,
as also in different implementations of the diffuse interface approach (e.g.
phase field models based on the Cahn-Hilliard and Cahn-Allen equations) the
numerical procedure is typically based on the solution of a variable
coefficient Poisson equation, we believe that the methodology proposed here
can be helpful also to generalize the phase field theory in a low-Mach number
framework.

###  3.4 Time step restriction

The time step $\Delta t^{n+1}$ is estimated from the stability constraints of
the overall system:

$\Delta t^{n+1}=C_{\Delta t}\min(\Delta t_{c},\Delta t_{\sigma},\Delta
t_{\mu},\Delta t_{e})^{n+1}\mathrm{,}$ (57)

where $\Delta t_{c}$, $\Delta t_{\sigma}$, $\Delta t_{\mu}$ and $\Delta t_{e}$
are the maximum allowable time steps due to convection, surface tension,
momentum diffusion and thermal energy diffusion, respectively. These can be
determined as suggested in reference [65]:

$\displaystyle\Delta t_{c}$ $\displaystyle=\left(\dfrac{|u_{x,\max}|}{\Delta
x}+\dfrac{|u_{y,\max}|}{\Delta y}+\dfrac{|u_{z,\max}|}{\Delta
z}\right)^{-1}\mathrm{,}$ (58) $\displaystyle\Delta t_{\mu}$
$\displaystyle=\left[\max\left(\dfrac{\mu_{g}}{\tilde{\rho}_{g}^{m}},\dfrac{\mu_{l}}{\rho_{l}}\right)\left(\dfrac{2}{\Delta
x^{2}}+\dfrac{2}{\Delta y^{2}}+\dfrac{2}{\Delta
z^{2}}\right)\right]^{-1}\mathrm{,}$ $\displaystyle\Delta t_{\sigma}$
$\displaystyle=\sqrt{\dfrac{(\tilde{\rho}_{g}^{m}+\rho_{l})\min(\Delta
x,\Delta y,\Delta z)^{3}}{4\pi\sigma}}\mathrm{,}$ $\displaystyle\Delta t_{e}$
$\displaystyle=\left[\max\left(\dfrac{k_{g}}{\tilde{\rho}_{g}^{m}},\dfrac{k_{l}}{\rho_{l}c_{p,l}}\right)\left(\dfrac{2}{\Delta
x^{2}}+\dfrac{2}{\Delta y^{2}}+\dfrac{2}{\Delta
z^{2}}\right)\right]^{-1}\mathrm{,}$

where $|u_{i,\max}|$ is an estimate of the maximum value of the $i^{th}$
component of the flow velocity and $\tilde{\rho}_{g}^{m}$ is the minimum gas
density computed over the computational domain to account for the compressible
effects. For the cases presented here, $C_{\Delta t}=0.25$ was found
sufficient for a stable and accurate time integration and, unless otherwise
stated, this value has been employed for the validation cases.

##  4 Validation and testing

In order to validate and test the proposed numerical approach, five different
flow configurations are considered, denoted C1a and C1b C2, C3 and C4. The
first two simulations, C1a and C1b, reproduce the two-dimensional flow
originating in a fluid system made of alternating gaseous and liquid bands at
different initial densities and temperatures confined in a periodic, free-slip
channel. We believe that these test cases are particularly significant to
highlight the capabilities of the proposed numerical methodology. The third
simulation, C2, reproduces a two-dimensional gas bubble rising in an
incompressible liquid medium and is used as a quantitative validation against
a reference case from archival literature. We select this test case in order
to perform the comparison among the different methods analyzed in section 3.3.
The previous setup is also used to study the flow in the presence of three
rising bubbles, case C3. The last simulation, C4, reproduces a time-evolving,
plane mixing layer originating between two streams at different temperatures
and opposite velocities. One of the streams is assumed to be compressible, the
other being incompressible. The effect of the temperature gradients on the
temporal evolution of the mixing layer is fully described by means of the low-
Mach number asymptotic approach, taking into consideration thermal diffusion
as well as density gradients in the flow.

| $\mu_{l}/\mu_{g}$ | $k_{l}/k_{g}$ | $c_{v}/c_{v,g}$ | $c_{p}/c_{p,g}$ | $\tilde{\sigma}$ | $\tilde{\mathcal{R}}$ | $\tilde{p}$ | $\tilde{g}$
---|---|---|---|---|---|---|---|---
| - | - | - | - | $\mathrm{N/m}$ | $\mathrm{J/(kg\cdot K)}$ | $\mathrm{Pa}$ | $\mathrm{m/s^{2}}$
4.1 | $1$ | $1$ | $1$ | $1$ | $24.5$ | $300$ | $9\cdot 10^{4}$ | 0.0
4.2 | $10$ | $20$ | $5.76$ | $4.18$ | $24.5$ | $300$ | $9\cdot 10^{4}$ | 0.98
4.3 | $10$ | $20$ | $5.76$ | $4.18$ | $24.5$ | $300$ | $9\cdot 10^{4}$ | 0.98
4.4 | $1$ | $1$ | $1$ | $1$ | $0$ | $300$ | $9\cdot 10^{4}$ | 0
Table 1: Physical parameters of the fluids for cases C1a, C1b, C2, C3 and C4:
viscosity ratio, $\mu_{l}/\mu_{g}$, thermal conductivity ratio, $k_{l}/k_{g}$,
specific heat capacity ratios at constant pressure, $c_{p}/c_{p,g}$ and
constant volume $c_{v}/c_{v,g}$, reference surface tension coefficient,
$\tilde{\sigma}$, reference specific constant of the gas phase,
$\tilde{\mathcal{R}}$, reference thermodynamic pressure, $\tilde{p}$, and
reference gravitational acceleration, $\tilde{g}$. The subscript, “$g$”,
refers to the gaseous phase and the subscript “$l$” to the liquid.

###  4.1 Expansion of gas bands enclosed by an incompressible medium

The test case C1 reproduces the two-dimensional, isochoric (C1a) and isobaric
(C1b) transformation of a compressible gas band enclosed within an
incompressible liquid medium. All the quantities are provided in the non-
dimensional frame, the reference values and the simulation parameter being
reported in table 1. The domain is rectangular and extends, in non-dimensional
units, for $L_{x}/{\tilde{L}}=4$ and $L_{y}/{\tilde{L}}=0.5$, ${\tilde{L}}$
being the reference length scale. The domain is discretized using $N_{x}\times
N_{y}=128\times 16$ nodes. Free-slip boundary conditions are applied to the
lower and upper edges of the computational domain while an adiabatic, zero-
gradient boundary condition is prescribed to the temperature equation. A
periodic boundary condition is applied along the $x$ direction. The isochoric
case, C1a, considers a rectangular gas band of width $b/{\tilde{L}}=1$ that
splits the domain into two parts filled by an incompressible liquid. The band
is initially centred around the axial position $x_{c}/{\tilde{L}}=3$ and
extends over the whole domain in the $y$ direction.

Figure 1: a) Schematic of the computational domain and initial conditions for
the test case C1a. b) Non-dimensional thermodynamic pressure,
$p_{0}/{\tilde{p}}$, in the gas regions as a function of the non-dimensional
time, $t/{\tilde{t}}$. c) Non-dimensional temperature, $T/T_{l,i}$, computed
on the middle-line of the domain as a function of non-dimensional coordinate
$x/{\tilde{L}}$. d) Non-dimensional density, $\rho/\rho_{l,i}$, computed on
the middle-line as a function of $x/{\tilde{L}}$. e) Volume fraction, $\Phi$,
computed on the middle-line as a function of $x/{\tilde{L}}$. The temperature,
density and volume fraction curves are provided for four different time
instants, $t/{\tilde{t}}=0.001$, $t/{\tilde{t}}=0.01$, $t/{\tilde{t}}=0.1$ and
$t/{\tilde{t}}=1$, ${\tilde{t}}$ being the reference time scale.

The geometrical configuration of the problem is provided in Fig. 1. The ratio
of the gas to the liquid temperature is initially fixed to
$(T_{g}/T_{l})_{i}=5/6$. The initial temperature and density fields are
uniform over each band, the only discontinuities being located on the liquid-
gas interface. Fig. 1 provides also the temporal evolution of the
thermodynamic pressure for the gas phase together with a plot of the density,
temperature and volume fraction as a function of the axial position,
$x/{\tilde{L}}$, at four different time instants. As a result of the initial
temperature gradient, a heat flux develops from the liquid region towards the
gas band. The temperature of the latter increases as shown in Fig. 1, while
its density decreases as can be observed in Fig. 1. It should be noted that,
the prescribed boundary conditions do not allow any volume change of the gas
region. Hence, the transformation is isochoric, the volume fraction field
remains unchanged and the gas band does not change the position of its
centroid neither its boundaries during the transient as can be seen in Fig. 1.
In these conditions, the energy transfer to the gas band enforces the
thermodynamic pressure to progressively increase until a uniform temperature
field is established over the entire domain. At this point, the thermodynamic
pressure settles to a constant value.

Figure 2: a) Schematic of the computational domain and initial conditions for
the test case C1b. b) Thermodynamic pressure, $p_{0}/{\tilde{p}}$, in the gas
regions as a function of the non-dimensional time, $t/{\tilde{t}}$. c)
Temperature, $T/T_{l,i}$, as a function of $x/{\tilde{L}}$. d) Density,
$\rho/\rho_{l,i}$, as a function of $x/{\tilde{L}}$. e) Volume fraction,
$\Phi$, as a function of $x/{\tilde{L}}$. The temperature, density and volume
fraction are provided for four different time instants, $t/{\tilde{t}}=5$,
$t/{\tilde{t}}=10$, $t/{\tilde{t}}=15$ and $t/{\tilde{t}}=20$.

In the second test case, C1b, we address the simulation of the isobaric
contraction and expansion of two separated gaseous bands enclosed within an
incompressible liquid medium. The domain size, discretization and boundary
conditions are unchanged and similarly for the fluid parameters which are
provided in table 1. The initial configuration of the fluid system is provided
in Fig. 2. Initially, the two rectangular gas bands extend over a length
$b/{\tilde{L}}=1$ along the $x$ direction, separated by an incompressible
liquid. The left-side band centroid is located at the axial position
$x_{c1}/{\tilde{L}}=1$ while the right-side band is centered around
$x_{c2}/{\tilde{L}}=3$. The ratio between the initial gas and liquid
temperatures is set to $(T_{g,1}/T_{l})_{i}=4/3$ for the left-side band and to
$(T_{g,2}/T_{l})_{i}=2/3$ for the right-side region. The initial temperature
and density fields are uniform over each of the five different regions
composing the fluid system. The initial thermodynamic pressure is the same in
the two gaseous regions. Fig. 2 and 2 provide the temperature and the density
fields as a function of $x/{\tilde{L}}$ at four different time instants. The
colder, right-side, band absorbs energy from the surrounding liquid medium
while the hotter band on the left-side releases energy to the surrounding
fluid. Hence, we observe the expansion of the colder gas band, simultaneously
with the equivalent compression of the hotter gas, as shown in Fig. 2
providing the volume fraction as a function of $x/{\tilde{L}}$ at four
different time instants. The volume of the liquid region included between the
two gaseous bands cannot change; however, due to the periodic boundary
condition along the $x$ direction, the liquid fluid can move from the right to
left side of the domain. The two bands do not change the position of their
center of mass during the expansion and contraction. Since in this case the
volume of each band can vary freely, the transformation is isobaric as can be
seen in Fig. 2. Even if we cannot provide an analytical solution for the cases
C1a and C1b, we believe that their numerical outcomes clearly show the
capability of the method to account for heat transfer, density and temperature
gradients as well as for compressibility effects in both isobaric and
isochoric conditions in the low-Mach number regime.

###  4.2 Rising bubble

The test case C2 addresses the simulation of a two-dimensional rising bubble
flow. A circular gaseous bubble of initial density $\rho_{g}$ and temperature
$T_{g}$ is immersed in a liquid fluid with a higher, constant density,
$\rho_{l}$, and temperature, $T_{l}$. Both the temperature and the density
fields are initially uniform within the bubble and the liquid phase while a
discontinuity exists across the interface. The initial configuration is
displayed in Fig. 3.

Figure 3: Schematic of the computational domain and initial configuration used
for the rising bubble simulation.

The rectangular computational domain extends, in non-dimensional units, for
$L_{x}/{\tilde{L}}\times L_{y}/{\tilde{L}}=1\times 2$. The domain is
discretized using $N_{x}\times N_{y}=128\times 256$ nodes. The initial bubble
diameter is $d_{i}/{\tilde{L}}=0.5$ while the bubble center is initially
located at $\mathbf{X}_{c,i}/\tilde{L}=(0.5,0.5)$. A no-slip and no-
penetration boundary condition is prescribed to the momentum equation along
the lower and the upper edges of the domain while a zero-gradient boundary
condition is applied to the temperature equation. A periodic boundary
condition is prescribed along the $x$ direction. The physical parameters of
the fluids are provided in table 1. We report the results of five different
test cases with different initial temperature ratios, $(T_{l}/T_{g})_{i}=1$,
$1.2$, $1.5$, $2$ and $3$, and corresponding density ratios,
$(\rho_{l}/\rho_{g})_{i}=10$, $8.33$, $6.67$, $5$ and $3.33$. We use as output
quantities the center of mass of the bubble and the bubble rising velocity.
The bubble centroid is defined as

$\mathbf{X}_{c}=(x_{c},y_{c})=\frac{\int_{V_{g}}\mathbf{x}dV}{\int_{V_{g}}dV}.$
(59)

In a similar fashion, the bubble rising velocity is defined as the mean
velocity with which the gas phase is moving,

$\mathbf{U}_{c}=(u_{c},v_{c})=\frac{\int_{V_{g}}\mathbf{u}dV}{\int_{V_{g}}dV}.$
(60)

Figure 4: Rising bubble: vertical position of the bubble centroid,
$y_{c}(t)/{\tilde{L}}$, and vertical component of the bubble rising velocity,
$v_{c}(t)/{\tilde{U}}$, versus time, $t/{\tilde{t}}$. The continuous black
line represents the outcomes of the present isothermal case whereas the dotted
line provides the results by Hysing et al. [66].

Fig. 4 displays the vertical position of the bubble centroid,
$y_{c}(t)/{\tilde{L}}$, and the vertical component of the bubble rising
velocity, $v_{c}(t)/{\tilde{U}}$, versus time, $t/{\tilde{t}}$, for each of
the initial temperature ratios given above. In the isothermal case the initial
temperature field is uniform over the entire domain, $(T_{l}/T_{g})_{i}=1$,
and the density ratio is set to $(\rho_{l}/\rho_{g})_{i}=10$. The results of
the present simulation are compared with that obtained by Hysing et al. [66].
As the density ratio, $(\rho_{l}/\rho_{g})_{i}$, is decreased, the rising
velocity of the bubble is initially lower than in the isothermal case due to
the lower buoyancy force exerted by the liquid on the gas bubble. Nonetheless,
the thermal diffusion reduces progressively the temperature gradient between
the two phases. The bubble heats-up and the density ratio,
$\rho_{l}/\rho_{g}$, increases. As a result, after an initial transient, the
terminal bubble rising velocity tends to settle to the same regime velocity as
that of the isothermal case, independently of the initial density ratio,
$(\rho_{l}/\rho_{g})_{i}$. Clearly, the initial differences in the rising
velocities lead to an offset in the position of the bubble centroid.

To conclude the analysis of this test case, we compare the three methods
presented in section 3.3 in terms of number of iterations needed to solve the
Poisson equation. Since the three methods require different tolerances to
satisfy the divergence constraint with the same accuracy, we set
$\varepsilon_{t}=10^{-8}$ for Method 1 and $\varepsilon=10^{-11}$ for Method 2
and 3 for a fair comparison. Using these different thresholds,
$\varepsilon_{t}$, leads to similar values of the residual (below $10^{-14}$),
computed as the difference between the velocity divergence and its constraint
according to Eq. (50). For what concerns Method 3, and only for this case, we
compute the residual $\varepsilon$ on the pressure between two consecutive
iterations and not on the velocity divergence. Once again, this choice is
motivated by a fair comparison with the other two methods. Indeed, for Method
1 and Method 2, the residual based on the pressure $p_{2}$ is the only
possible choice being the correction step performed only at the end. The
results provided in Fig. 5 clearly show that the current approach, Method 3,
requires a number of iterations between $1.5$ and $3$ times lower than that of
Method 1 and Method 2 to achieve a full convergence. As mentioned above, we
attribute this faster convergence to the exact way we recast the variable-
coefficient Poisson equation into a constant-coefficient problem using
directly the correction step.

Figure 5: Number of iterations as a function of time $t/\tilde{t}$ required to
solve the pressure Poisson equation for the rising bubble test case using the
method by Bartholomew and Laizet [64], Motheau and Abraham [63] and the
present method. The reference time scale $\tilde{t}$ is set equal to
$\sqrt{d_{0}/\mathbf{g}}$.

###  4.3 Multiple rising bubbles

In this section, we consider the same configuration adopted in the previous
test case to study the flow in the presence of three compressible bubbles,
rising in an incompressible liquid. The bubbles have the same initial
diameter, $d_{i}$ and are initially at rest in a rectangular domain of
dimensions $L_{x}/\tilde{L}\times L_{y}/\tilde{L}=7.2\times 12.8$, being the
reference length $\tilde{L}=d_{i}$. The initial position of the bubble
centroids are set to $(\mathbf{X}_{c,1}/\tilde{L})_{i}=(1.8,1.25)$ for Bubble
n.1, $(\mathbf{X}_{c,2}/\tilde{L})_{i}=(3.6,1.25)$ for Bubble n.2 and
$(\mathbf{X}_{c,3}/\tilde{L})_{i}=(5.4,1.25)$ for Bubble n.3. The initial
temperature and density of the liquid phase are $T_{l}$ and $\rho_{l}$,
respectively. In order to highlight the compressibility effects, the three
bubbles are initialized at three different temperatures,
$(T_{g,1}/T_{l})_{i}=1.5$, $(T_{g,2}/T_{l})_{i}=1.0$ and
$(T_{g,3}/T_{l})_{i}=0.75$ as reported on the left panel of Fig. 6. In order
to avoid bubble coalescence and merging, we consider a limited Weber number
$We=\rho_{g,0}^{m}\tilde{U}^{2}d_{0}/\tilde{\sigma}=0.125$ and we set the
Reynolds number $Re=\rho_{g,0}^{m}\tilde{U}d_{0}/\mu_{g}=125$, where
$\tilde{U}=\sqrt{|\mathbf{g}|d_{0}}$ and $\rho_{g,0}^{m}$ is the minimum
initial gas density. The Prandtl number $Pr=\mu_{g}k_{g}/c_{p,g}$ is set to
$0.7$. All the other dimensionless parameters are kept the same as in the
previous case and are reported in table 1.

Figure 6: Position of the three rising bubbles at $t/\tilde{t}=0$ (left panel)
and for $t/\tilde{t}>0$ (right panel). The interface position is taken from
the grid points where $\Phi=0.5$ and the bubble contour are plotted with at
the dimensionless physical time
$t/\tilde{t}=\\{0,\,1.5,\,3.0,\,4.5,\,6.0,\,7.5,\,9.0,\,10.5,\,12.0,\,13.5\\}$,
with the reference time scale being $\tilde{t}=\sqrt{d_{0}/|\mathbf{g}|}$.

Figure 7: Averaged gas temperature (left panel) and averaged gas density of
the three bubbles (right panel) versus non-dimensional time. The initial
temperature and density of Bubble n.2 are taken as a reference temperature,
$T_{iso}$, and density, $\rho_{iso}$. The reference time scale is
$\tilde{t}=\sqrt{d_{0}/|\mathbf{g}|}$.

Since the system is closed and thermally isolated, once the bubbles start to
rise, the heat transfer exchanged among each other and with the liquid medium
drives them towards the thermal equilibrium. In detail, Bubble n.1 starts to
cool down, Bubble n.2 maintains an almost constant average temperature,
whereas Bubble n.3 is heated up. As a result and owing to the variation of the
thermodynamic pressure, the first bubble contracts increasing its mean
density, the third bubble expands decreasing its mean density, whereas the
second one slightly expands mainly due to the variation of the thermodynamic
pressure. As shown in Fig. 7, after $t/\tilde{t}\approx 6$, being
$\tilde{t}=\sqrt{d_{0}/|\mathbf{g}|}$, the thermal equilibrium is globally
reached and the mean temperature and density remain approximately constant for
the three bubbles.

Figure 8: Normalized vertical velocity of the bubbles versus time. The
reference velocity scale is $\tilde{U}=\sqrt{|\mathbf{g}|d_{0}}$ while the
reference time scale is $\tilde{t}=\sqrt{d_{0}/|\mathbf{g}|}$.

Fig. 8 provides the mean vertical velocity of the centroid of each bubble. The
initial expansion and contraction of the bubbles affects the vertical
component of their rising velocities computed as in Eq. (60). In particular,
until $t/\tilde{t}\approx 1.8$, all the three bubbles move with a comparable
vertical velocity. After this initial stage, the third bubble starts to
accelerate and arrives first at the top wall, whereas the first one starts to
decelerate and moves along the vertical direction at an almost constant speed.
On the other hand, the second bubble accelerates towards the top wall, but at
a lower rate than the third one. The physical explanation for this behavior
relies in the modification induced by the initial expansion and contraction
stage of the bubbles, which determines an increase of the buoyancy forces for
Bubble n.3 and a reduction for Bubble n.1.

###  4.4 Mixing layer

As a final test case, C4, the numerical simulation of a two-dimensional,
temporal mixing layer is addressed. This considered flow configuration
develops in the region between two counter-directional flows, one of them
being compressible, the other incompressible. The streams move with opposite
velocities, $U_{g}$ and $U_{l}$. In these conditions, a Kelvin-Helmholtz
instability promotes the formation of well-defined coherent vortices in the
region separating the two streams. The latter enhance micro-mixing and
molecular diffusion promoting the exchange of momentum and energy between the
opposite streams. The computational domain consists of a square box of unit
size, $L_{x}/\tilde{L}\times L_{y}/\tilde{L}=1\times 1$, discretized using
$N_{x}\times N_{y}=512\times 512$ nodes. In the lower part of the
computational domain, $0<y/{\tilde{L}}\leq 0.5$, the incompressible flow moves
from the right to the left while in the upper part of the domain,
$0.5<y/{\tilde{L}}\leq 1$, the compressible stream moves in the opposite
direction. A no-slip boundary condition is prescribed to the momentum equation
along the top and bottom sides of the domain while a zero-gradient, adiabatic
boundary condition is assigned to the temperature equations along the same
boundaries. Periodic boundary conditions are applied to all quantities along
the flow direction, $x$.

To better characterize the mixing-layer flow, it is worth introducing a
length-scale based on the initial vorticity thickness in the mixing layer,
$\delta$, and a corresponding Reynolds number
$Re_{\delta}=U_{c}\delta/\nu_{g,i}$ with $\nu_{g,i}$ the kinematic viscosity
of the gas phase (evaluated at the initial condition) and $U_{c}$ a prescribed
convective velocity defined as $U_{c}=1/2(U_{g}-U_{l})$. The initial velocity
field is prescribed imposing a pseudo-perturbation on a mean profile according
to the following relations [67]:

$\displaystyle\dfrac{u(x,y,0)}{U_{c}}=\tanh\left(\frac{2{\tilde{L}}}{\delta}y\right)+\xi_{noise}\frac{\partial\psi}{\partial
y},$ (61)
$\displaystyle\dfrac{v(x,y,0)}{U_{c}}=-\xi_{noise}\frac{\partial\psi}{\partial
x},$ (62)
$\displaystyle\psi(x,y)=\exp\left(-\frac{{\tilde{L}}^{2}}{\delta^{2}}y^{2}\right)\left[\cos(4\pi
x)+0.03\sin(10\pi x)\right],$ (63)

where $u(x,y,0)$ and $v(x,y,0)$ are the horizontal and vertical components of
the initial velocity field. Moreover, the factor $\xi_{noise}=10^{-3}$ is
chosen to ensure that the velocity perturbations remain a small percentage of
the mean velocity, as suggested by the authors in Zayernouri et al. [67].
Prescribing the hyperbolic tangent velocity profile given by Eq. (61) - (62),
the wave-length associated with the initial vortex distribution results to be
approximately $\lambda\simeq 7\delta$ [68, 67]. Hence, given the domain size,
${\tilde{L}}$ and the desired number of vortexes in the periodic domain, $N$,
the initial vortex thickness is $\delta/{\tilde{L}}=1/(7N)$. In the present
case, the initial vorticity thickness is fixed to $\delta/{\tilde{L}}=1/28$
and $Re_{\delta}=200$. The non-dimensional viscosity, thermal conductivity and
the specific heat capacity ratios are kept equal to unity, while the density
ratio based on the initial gas density $(\rho_{l}/\rho_{g})_{i}$ is taken
equal to $5$. Finally, the Prandtl number is set to
$Pr=\mu_{g}c_{p,g}/k_{g}=8.92$ with $c_{p,g}$ and $k_{g}$ being the specific
heat capacity and thermal conductivity of the gas phase. The initial
temperature field is initialised according to the step-function,

$T(x,y)=\begin{cases}T_{l,i},&\mbox{if}\ 0\leq y/{\tilde{L}}\leq 0.5,\\\
T_{g,i},&\mbox{if}\ 0.5<y/{\tilde{L}}\leq 1.\end{cases}$ (64)

Figure 9: Sketch of the domain for the temporal mixing-layer simulation
showing the initial velocity and temperature fields.

While keeping fixed the initial liquid temperature $T_{l,i}$ and the initial
density ratio $(\rho_{l}/\rho_{g})_{i}$ between the two phases, different
initial gas temperatures, $T_{g,i}$, gas densities, $\rho_{g,i}$ and liquid
densities $\rho_{l,i}$, are prescribed to the compressible and incompressible
fluids, as sketched in Fig. 9. A first test case considers the isothermal flow
where $(T_{g}/T_{l})_{i}=1$, whereas three other cases address a temperature
ratio equal to $15/16$, $5/6$ and $3/4$, respectively.

Figure 10: Temporal evolution of the thermodynamic pressure, $p_{0}$, mean gas
temperature, $\overline{T}_{g}$, mean gas density, $\overline{\rho}_{g}$ and
mean kinetic energy, $\overline{E}_{k}$ for four different temperature ratios,
$(T_{g}/T_{l})_{i}\in[1.0,15/16,5/6,3/4]$. All the quantities are non-
dimensional using as reference values the quantities of the isothermal case,
except for the mean kinetic energy where we employ the initial value of
$\overline{E}_{k}$ of each case, $\overline{E}_{k,i}$. The reference time-
scale $\tilde{t}$ is set equal to $\tilde{L}/\tilde{U}_{c}$.

Fig. 10 provides the temporal evolution of the thermodynamic pressure (uniform
over the computational domain), the mean gas and liquid temperature, the mean
gas density and the mean kinetic energy for the two phases. The gas and liquid
temperatures are computed as integrals over the corresponding domains, whereas
the mean kinetic energy is estimated over the compressible and incompressible
regions:

$\displaystyle{\overline{T}}_{g}(t)=\frac{1}{V_{g}}\int_{V_{g}}T(x,y,t)\Phi(x,y,t)\
dx\ dy,$ (65)
$\displaystyle{\overline{T}}_{l}(t)=\frac{1}{V_{l}}\int_{V_{l}}T(x,y,t)(1-\Phi(x,y,t))\
dx\ dy,$ (66)
$\displaystyle{\overline{E}_{k}}(t)=\frac{1}{2L_{x}L_{y}}\int_{-L_{x}/2}^{L_{x}/2}\int_{-L_{y}/2}^{L_{y}/2}\rho(x,y,t){\mathbf{u}(x,y,t)}\cdot{\mathbf{u}}(x,y,t)\
dx\ dy,$ (67)

where $V_{g}$ is computed using Eq. (45) and $V_{l}=L_{x}L_{y}-V_{g}$, being
the domain closed. Once the mean gas temperature is known, the mean gas
density is computed directly from the equation of state whereas the liquid
density is constant and equal to $\rho_{l}$,

${\bar{\rho}}_{g}(t)=\dfrac{p_{0}(t)}{\mathcal{R}\overline{T}_{g}(t)}\mathrm{,}\hskip
14.22636pt\text{and}\hskip 14.22636pt{\bar{\rho}}_{l}(t)=\rho_{l}\mathrm{.}$
(68)

Figure 11: Contour plots of the dimensionless vorticity field
$\omega_{z}/\tilde{\omega}_{c}$ in the mixing layer for the isothermal case
(left panels) and for the $(T_{g}/T_{l})_{i}=3/4$ case (right panel).The
contour plots refer to the dimensionless physical time $t/\tilde{t}=0.01$,
$0.02$ and $0.03$, in order. The reference time scale is
$\tilde{t}=\tilde{L}/U_{c}$ whereas the vorticity one is
$\tilde{\omega}_{c}=U_{c}/\tilde{L}$.

In the isothermal case, the thermodynamic pressure, the mean liquid and gas
temperature and the mean gas density do not change over time as shown in Fig.
10. Moreover with the prescribed boundary conditions and in absence of
external forces, the mean kinetic energy in the incompressible case ( e.g.
$\left(T_{g}/T_{l}\right)=1$) monotonically decreases due to the internal
dissipation in the flow. On the other hand, as the temperature ratio is
reduced below unity, the turbulent mixing enhances the thermal diffusion
between the two fluids thus rapidly reducing the temperature gradients. As a
result, the temperature tends to rapidly increase in the colder, compressible
stream until a stationary condition is established. As the compressible phase
heats up, the thermodynamic pressure, $p_{0}$, and the mean gas density
increase. These effects modify the mean kinetic energy balance, that in the
compressible cases contains not only the viscous dissipation but also a
pressure work term proportional to the gas expansion. This last term modifies
the variation of $E_{k}/E_{k,i}$ for all the compressible cases and is
responsible for the initial increase in the mean kinetic energy observed for
the case $\left(T_{g}/T_{l}\right)_{i}=0.75$ up to $t/\tilde{t}\approx 1$,
$\tilde{t}=\tilde{L}/U_{c}$ being the reference time scale. However, once the
temperature gradients become negligible and the thermodynamic pressure has
reached a constant value, the compressible effects expire and the mean kinetic
energy variation is mainly governed by the viscous dissipation.

Finally, Fig. 11 displays the contour plots of the instantaneous vorticity
field at three different physical times for isothermal case and for the
initial temperature ratio $\left(T_{g}/T_{l}\right)_{i}=3/4$. Despite we limit
the analysis at the initial times in order to avoid the loss in resolution
induced by the formation of smaller and smaller scales, we observe that, in
general, the presence of a temperature gradient enhances the mixing and the
growth-rate of the vorticity thickness with respect to the reference,
isothermal case.

##  5 Final remarks

Multiphase, compressible flows are of great interest in a wide range of
scientific fields and engineering problems. In this context, we propose a
novel approach to the numerical simulation of multiphase, viscous flows where
a compressible gas phase and an incompressible liquid mutually interact in the
low-Mach number regime. The problem is addressed in the framework of a low-
Mach number asymptotic expansion of the compressible formulation of the
Navier-Stokes equations. In this limit, acoustics is neglected but large
density variations of the gas phase can be accounted for as well as heat
transfer between the phases and with the domain boundaries. A Volume of Fluid
approach is used to deal with the presence of different phases in the flow as
well as for interface tracking. The interface reconstruction is based on the
MTHINC method [37] while the effect of the surface tension is accounted for
using the continuum surface force (CSF) model [47]. The same set of equations
is used for both the gas and the liquid phase, the zero-divergence condition
being exactly imposed to the latter. To numerically solve this set of
equations, we have developed a massive parallel solver, second order accurate
both in time and space. The Poisson pressure equation is managed by a FFT-
based solver that allows for a numerically efficient and very fast solution
procedure. In addition, this choice is suited for code optimization and
adaptation of incompressible GPU codes that benefits of FFT-based solvers
(e.g. see [69]). The proposed iterative procedure shows to be more efficient
in terms of number of iterations than the two approaches available in
literature in the context of low-Mach number flows [63, 64]. The solver has
been build upon a code for incompressible flows which has undergone an
extensive validation campaign [7, 48, 49]. We provide a detailed and complete
description of the theoretical approach together with information about the
numerical technique and implementation details. In addition, we apply the
described numerical approach to the simulation of five different flow
configurations. The outcomes of two simulations reproducing the two-
dimensional expansion and contraction of rectangular gaseous bands enclosed in
an incompressible fluid and confined in a free-slip, periodic channel are
provided. Next, we address the simulation of two-dimensional rising bubbles.
First, we consider a single bubble and compare the results of our simulation
with the reference data by Hysing et al. [66] using as benchmark quantities
the bubble centroid and the bubble rising velocity. Second, we simulate the
evolution of three bubbles of the same size but with different initial
temperatures. Finally we discuss the outcome of a numerical simulation
reproducing a plane, temporal mixing layer and show how the compressibility of
the gas phase alters the development of the instability.

As the proposed mathematical and numerical framework is independent of the
capturing/tracking technique used to describe the interface topology, the
proposed methodology can be directly extended to other existing numerical
codes. We believe that the results presented here demonstrate that it is
possible to accurately address the numerical simulation of multiphase, viscous
flows in the low-Mach number regime, also when one of the phases can be
treated as incompressible. Further extensions of the present methodology may
concern the addition of more complex physical phenomena like phase change and
complex interfacial thermodynamics, as absorption-desorption processes.

## Acknowledgements

N.S., M.E.R. and L.B. acknowledge the support from the Swedish Research
Council via the multidisciplinary research environment INTERFACE, Hybrid
multiscale modelling of transport phenomena for energy efficient processes and
the Grant No. 2014-5001. The computational resources were provided by SNIC
(Swedish National Infrastructure for Computing).

## References

  * Hirt and Nichols [1981] C. W. Hirt, B. D. Nichols,  Volume of fluid (vof) method for the dynamics of free boundaries,  Journal of computational physics 39 (1981) 201–225.
  * Elghobashi [2019] S. Elghobashi,  Direct numerical simulation of turbulent flows laden with droplets or bubbles,  Annual Review of Fluid Mechanics 51 (2019) 217–244.
  * Tanguy et al. [2014] S. Tanguy, M. Sagan, B. Lalanne, F. Couderc, C. Colin,  Benchmarks and numerical methods for the simulation of boiling flows,  Journal of Computational Physics 264 (2014) 1–22.
  * Le Martelot et al. [2014] S. Le Martelot, R. Saurel, B. Nkonga,  Towards the direct numerical simulation of nucleate boiling flows,  International Journal of Multiphase Flow 66 (2014) 62–78.
  * Salvador et al. [2016] F. Salvador, J. Romero, M. Roselló, D. Jaramillo,  Numerical simulation of primary atomization in diesel spray at low injection pressure,  Journal of Computational and Applied Mathematics 291 (2016) 94–102.
  * Duret et al. [2013] B. Duret, J. Reveillon, T. Menard, F. Demoulin,  Improving primary atomization modeling through dns of two-phase flows,  International Journal of Multiphase Flow 55 (2013) 130–137.
  * Picano et al. [2015] F. Picano, W. P. Breugem, L. Brandt,  Turbulent channel flow of dense suspensions of neutrally buoyant spheres,  Journal of Fluid Mechanics 764 (2015) 463–487.
  * Balcázar et al. [2016] N. Balcázar, O. Lehmkuhl, L. Jofre, J. Rigola, A. Oliva,  A coupled volume-of-fluid/level-set method for simulation of two-phase flows on unstructured meshes,  Computers & Fluids 124 (2016) 12–29.
  * Cottet and Maitre [2016] G. H. Cottet, E. Maitre,  A semi-implicit level set method for multiphase flows and fluid-structure interaction problems,  Journal of Computational Physics 314 (2016) 80–92.
  * Ardekani et al. [2019] M. Ardekani, M. Rosti, L. Brandt,  Turbulent flow of finite-size spherical particles in channels with viscous hyper-elastic walls,  Journal of Fluid Mechanics 873 (2019) 410–440.
  * Costa et al. [2018] P. Costa, F. Picano, L. Brandt, W. P. Breugem,  Effects of the finite particle size in turbulent wall-bounded flows of dense suspensions,  Journal of fluid mechanics 843 (2018) 450–478.
  * Rosti et al. [2019] M. E. Rosti, F. De Vita, L. Brandt,  Numerical simulations of emulsions in shear flows,  Acta Mechanica 230 (2019) 667–682.
  * Dalla Barba and Picano [2018] F. Dalla Barba, F. Picano,  Clustering and entrainment effects on the evaporation of dilute droplets in a turbulent jet,  Physical Review Fluids 3 (2018) 034304.
  * Rosti et al. [2019] M. E. Rosti, F. Picano, L. Brandt,  Numerical approaches to complex fluids,  in: Flowing Matter, Springer, 2019, pp. 1–34.
  * Hirt et al. [1974] C. W. Hirt, A. A. Amsden, J. Cook,  An arbitrary lagrangian-eulerian computing method for all flow speeds,  Journal of computational physics 14 (1974) 227–253.
  * Hughes et al. [1981] T. J. Hughes, W. K. Liu, T. K. Zimmermann,  Lagrangian-eulerian finite element formulation for incompressible viscous flows,  Computer methods in applied mechanics and engineering 29 (1981) 329–349.
  * Ganesan and Tobiska [2012] S. Ganesan, L. Tobiska,  Arbitrary lagrangian–eulerian finite-element method for computation of two-phase flows with soluble surfactants,  Journal of Computational Physics 231 (2012) 3685–3702.
  * Unverdi and Tryggvason [1992] S. O. Unverdi, G. Tryggvason,  A front-tracking method for viscous, incompressible, multi-fluid flows,  Journal of computational physics 100 (1992) 25–37.
  * Peskin [1972] C. S. Peskin,  Flow patterns around heart valves: a numerical method,  Journal of computational physics 10 (1972) 252–271.
  * Breugem [2012] W. P. Breugem,  A second-order accurate immersed boundary method for fully resolved simulations of particle-laden flows,  Journal of Computational Physics 231 (2012) 4469–4498.
  * Mittal and Iaccarino [2005] R. Mittal, G. Iaccarino,  Immersed boundary methods,  Annual Review of Fluid Mechanics 37 (2005) 239–261.
  * Haeri and Shrimpton [2012] S. Haeri, J. Shrimpton,  On the application of immersed boundary, fictitious domain and body-conformal mesh methods to many particle multiphase flows,  International Journal of Multiphase Flow 40 (2012) 38–55.
  * Kempe and Fröhlich [2012] T. Kempe, J. Fröhlich,  An improved immersed boundary method with direct forcing for the simulation of particle laden flows,  Journal of Computational Physics 231 (2012) 3663–3684.
  * Sussman and Ohta [2009] M. Sussman, M. Ohta,  A stable and efficient method for treating surface tension in incompressible two-phase flow,  SIAM Journal on Scientific Computing 31 (2009) 2447–2471.
  * Sussman et al. [1994] M. Sussman, P. Smereka, S. Osher,  A level set approach for computing solutions to incompressible two-phase flow,  Journal of Computational physics 114 (1994) 146–159.
  * Osher and Sethian [1988] S. Osher, J. A. Sethian,  Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations,  Journal of computational physics 79 (1988) 12–49.
  * Gibou et al. [2018] F. Gibou, R. Fedkiw, S. Osher,  A review of level-set methods and some recent applications,  Journal of Computational Physics 353 (2018) 82–109.
  * Sethian and Smereka [2003] J. A. Sethian, P. Smereka,  Level set methods for fluid interfaces,  Annual review of fluid mechanics 35 (2003) 341–372.
  * Scardovelli and Zaleski [1999] R. Scardovelli, S. Zaleski,  Direct numerical simulation of free-surface and interfacial flow,  Annual review of fluid mechanics 31 (1999) 567–603.
  * Gueyffier et al. [1999] D. Gueyffier, J. Li, A. Nadim, R. Scardovelli, S. Zaleski,  Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows,  Journal of Computational physics 152 (1999) 423–456.
  * O’Brien and Bussmann [2018] A. O’Brien, M. Bussmann,  A volume-of-fluid ghost-cell immersed boundary method for multiphase flows with contact line dynamics,  Computers & Fluids 165 (2018) 43–53.
  * Noh and Woodward [1976] W. F. Noh, P. Woodward,  Slic (simple line interface calculation),  in: Proceedings of the fifth international conference on numerical methods in fluid dynamics June 28–July 2, 1976 Twente University, Enschede, Springer, 1976, pp. 330–340.
  * Youngs [1982] D. L. Youngs,  Time-dependent multi-material flow with large fluid distortion,  Numerical methods for fluid dynamics (1982).
  * Youngs [1984] D. L. Youngs,  An interface tracking method for a 3d eulerian hydrodynamics code,  Atomic Weapons Research Establishment (AWRE) Technical Report 44 (1984) 35\.
  * Aulisa et al. [2003] E. Aulisa, S. Manservisi, R. Scardovelli, S. Zaleski,  A geometrical area-preserving volume-of-fluid advection method,  Journal of Computational Physics 192 (2003) 355–364.
  * Xiao et al. [2005] F. Xiao, Y. Honma, T. Kono,  A simple algebraic interface capturing scheme using hyperbolic tangent function,  International Journal for Numerical Methods in Fluids 48 (2005) 1023–1040.
  * Ii et al. [2012] S. Ii, K. Sugiyama, S. Takeuchi, S. Takagi, Y. Matsumoto, F. Xiao,  An interface capturing method with a continuous function: The thinc method with multi-dimensional reconstruction,  Journal of Computational Physics 231 (2012) 2328–2358.
  * Pelanti [2017] M. Pelanti,  Low mach number preconditioning techniques for roe-type and hllc-type methods for a two-phase compressible flow model,  Applied Mathematics and Computation 310 (2017) 112–133.
  * Haimovich and Frankel [2017] O. Haimovich, S. H. Frankel,  Numerical simulations of compressible multicomponent and multiphase flow using a high-order targeted eno (teno) finite-volume method,  Computers & Fluids 146 (2017) 105–116.
  * Lunati and Jenny [2006] I. Lunati, P. Jenny,  Multiscale finite-volume method for compressible multiphase flow in porous media,  Journal of Computational Physics 216 (2006) 616–636.
  * Saurel and Lemetayer [2001] R. Saurel, O. Lemetayer,  A multiphase model for compressible flows with interfaces, shocks, detonation waves and cavitation,  Journal of Fluid Mechanics 431 (2001) 239–271.
  * Saurel and Abgrall [1999] R. Saurel, R. Abgrall,  A multiphase godunov method for compressible multifluid and multiphase flows,  Journal of Computational Physics 150 (1999) 425–467.
  * Colella and Pao [1999] P. Colella, K. Pao,  A projection method for low speed flows,  Journal of Computational Physics 149 (1999) 245–269.
  * Wall et al. [2002] C. Wall, C. D. Pierce, P. Moin,  A semi-implicit method for resolution of acoustic waves in low mach number flows,  Journal of Computational Physics 181 (2002) 545–563.
  * Müller [1998] B. Müller,  Low-mach-number asymptotics of the navier-stokes equations,  Journal of Engineering Mathematics 34 (1998) 97–109.
  * Majda and Sethian [1985] A. Majda, J. Sethian,  The derivation and numerical solution of the equations for zero mach number combustion,  Combustion science and technology 42 (1985) 185–205.
  * Brackbill et al. [1992] J. U. Brackbill, D. B. Kothe, C. Zemach,  A continuum method for modeling surface tension,  Journal of computational physics 100 (1992) 335–354.
  * Rosti and Brandt [2017] M. E. Rosti, L. Brandt,  Numerical simulation of turbulent channel flow over a viscous hyper-elastic wall,  Journal of Fluid Mechanics 830 (2017) 708–735.
  * Rosti et al. [2018] M. E. Rosti, L. Brandt, D. Mitra,  Rheology of suspensions of viscoelastic spheres: deformability as an effective volume fraction,  Physical Review Fluids 3 (2018) 012301.
  * Rosti et al. [2019] M. E. Rosti, Z. Ge, S. S. Jain, M. S. Dodd, L. Brandt,  Droplets in homogeneous shear turbulence,  Journal of Fluid Mechanics 876 (2019) 962–984.
  * Battista et al. [2014] F. Battista, F. Picano, C. M. Casciola,  Turbulent mixing of a slightly supercritical van der waals fluid at low-mach number,  Physics of Fluids 26 (2014) 055101.
  * Meister [1999] A. Meister,  Asymptotic single and multiple scale expansions in the low mach number limit,  SIAM Journal on Applied Mathematics 60 (1999) 256–271.
  * De Vita et al. [2020] F. De Vita, M. E. Rosti, S. Caserta, L. Brandt,  Numerical simulations of vorticity banding of emulsions in shear flows,  Soft Matter (2020).
  * De Vita et al. [2019] F. De Vita, M. E. Rosti, S. Caserta, L. Brandt,  On the effect of coalescence on the rheology of emulsions,  Journal of Fluid Mechanics 880 (2019) 969–991.
  * Puckett et al. [1997] E. G. Puckett, A. S. Almgren, J. B. Bell, D. L. Marcus, W. J. Rider,  A high-order projection method for tracking fluid interfaces in variable density incompressible flows,  Journal of computational physics 130 (1997) 269–282.
  * Scapin et al. [2020] N. Scapin, P. Costa, L. Brandt,  A volume-of-fluid method for interface-resolved simulations of phase-changing two-fluid flows,  Journal of Computational Physics (2020) 109251.
  * Castro et al. [2011] M. Castro, B. Costa, W. S. Don,  High order weighted essentially non-oscillatory weno-z schemes for hyperbolic conservation laws,  Journal of Computational Physics 230 (2011) 1766–1792.
  * Daru et al. [2010] V. Daru, P. Le Quéré, M.-C. Duluc, O. Le Maitre,  A numerical method for the simulation of low mach number liquid–gas flows,  Journal of Computational Physics 229 (2010) 8844–8867.
  * Costa [2018] P. Costa,  A fft-based finite-difference solver for massively-parallel direct numerical simulations of turbulent flows,  Computers & Mathematics with Applications 76 (2018) 1853–1862.
  * Wilhelmson and Ericksen [1977] R. B. Wilhelmson, J. H. Ericksen,  Direct solutions for poisson’s equation in three dimensions,  Journal of Computational Physics 25 (1977) 319–331.
  * Schumann and Sweet [1988] U. Schumann, R. A. Sweet,  Fast fourier transforms for direct solution of poisson’s equation with staggered boundary conditions,  Journal of Computational Physics 75 (1988) 123–137.
  * Li and Laizet [2010] N. Li, S. Laizet,  2decomp & fft-a highly scalable 2d decomposition library and fft interface,  in: Cray user group 2010 conference, 2010, pp. 1–13.
  * Motheau and Abraham [2016] E. Motheau, J. Abraham,  A high-order numerical algorithm for dns of low-mach-number reactive flows with detailed chemistry and quasi-spectral accuracy,  Journal of Computational Physics 313 (2016) 430–454.
  * Bartholomew and Laizet [2019] P. Bartholomew, S. Laizet,  A new highly scalable, high-order accurate framework for variable-density flows: Application to non-boussinesq gravity currents,  Computer Physics Communications 242 (2019) 83–94.
  * Kang et al. [2000] M. Kang, R. P. Fedkiw, X.-D. Liu,  A boundary condition capturing method for multiphase incompressible flow,  Journal of Scientific Computing 15 (2000) 323–360.
  * Hysing et al. [2009] S.-R. Hysing, S. Turek, D. Kuzmin, N. Parolini, E. Burman, S. Ganesan, L. Tobiska,  Quantitative benchmark computations of two-dimensional bubble dynamics,  International Journal for Numerical Methods in Fluids 60 (2009) 1259–1288.
  * Zayernouri and Metzger [2011] M. Zayernouri, M. Metzger,  Coherent features in the sensitivity field of a planar mixing layer,  Physics of Fluids 23 (2011) 025105.
  * Betchov and Szewczyk [1963] R. Betchov, A. Szewczyk,  Stability of a shear layer between parallel streams,  The physics of fluids 6 (1963) 1391–1396.
  * Costa et al. [2020] P. Costa, E. Phillips, L. Brandt, M. Fatica,  Gpu acceleration of cans for massively-parallel direct numerical simulations of canonical fluid flows,  Computers & Mathematics with Applications (2020).