# Modeling transport of scalars in two-phase flows with
a diffuse-interface method

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

###### Abstract

In this article, we propose a novel scalar-transport model for the simulation
of scalar quantities in two-phase flows with a phase-field method (diffuse-
interface method). In a two-phase flow, the scalar quantities typically have
disparate properties in two phases, which results in effective confinement of
the scalar quantities in one of the phases, in the time scales of interest.
This confinement of the scalars lead to the formation of sharp gradients of
the scalar concentration values at the interface, presenting a serious
challenge for its numerical simulations.

To overcome this challenge, we propose a model for the transport of scalars.
The model is discretized using a central-difference scheme, which leads to a
non-dissipative implementation that is crucial for the simulation of turbulent
flows. Furthermore, the provable strengths of the proposed model are: (a) the
model maintains the positivity property of the scalar concentration field, a
physical realizability requirement for the simulation of scalars, when the
proposed criterion is satisfied, (b) the proposed model is such that the
transport of the scalar concentration field is consistent with the transport
of the volume fraction field, which results in the enforcement of the
effective zero-flux boundary condition for the scalar at the interface; and
therefore, prevents the artificial numerical diffusion of the scalar across
the interface.

Finally, we present numerical simulations using the proposed model in a wide
range of two-phase flow regimes, spanning laminar to turbulent flows; and
assess: the accuracy and robustness of the model, the validity of the
positivity property of the scalar concentration field, and the enforcement of
the zero-flux boundary condition for the scalar at the interface.

###### keywords:

phase-field method , scalars , two-phase flows , conservative schemes ,
positivity , electrokinetics , turbulent flows

††journal: Journal of Computational Physics

##  1 Introduction

The transport of scalars in a two-phase flow is an important problem that
finds applications in wide range of natural phenomena and industrial
processes. A scalar quantity can represent: temperature field in modeling
boiling and evaporation phenomena (Villegas et al., 2016), dissolved gas
concentration in modeling oceanic carbon sequestration process (Lal, 2008),
salt concentration in modeling electrochemical systems (Fernandez et al.,
2014), surfactant concentration in modeling Marangoni effects (Takagi and
Matsumoto, 2011), etc.

The transport of scalars in a two-phase flow, typically, involves very
disparate length and time scales; and the scalar quantities often experience
very large and small diffusivities and mobilities in different phases. For
example, the diffusion coefficient of CO2 in air is $1.6\times 10^{-5}\
\mathrm{m^{2}/s}$; whereas, in water it is $1.6\times 10^{-9}\
\mathrm{m^{2}/s}$, which is $\approx 10,000$ times smaller than in air.
Similarly, the thermal diffusivity of air, at room temperature, is $1.9\times
10^{-5}\ \mathrm{m^{2}/s}$; whereas, that of water is $1.43\times 10^{-7}\
\mathrm{m^{2}/s}$, which is $\approx 100$ times smaller compared to that of
air. In a worst-case scenario, the diffusivities in the two phases could be so
different that the ratio of diffusivities in two phases could tend to
infinity. These disparate properties of the scalar quantities in different
phases result in scalars effectively being confined to one of the phases in
the time scales of interest. This poses a numerically challenging task to
resolve the gradient of the scalar concentration at the material interface,
and it usually leads to numerical leakage (artificial numerical diffusion) at
the interface and negative values of the scalar concentration field.

Several studies in the past have quantified, and tackled this issue of
numerical leakage in the context of a volume-of-fluid (VOF) method. Alke et
al. (2009) and Bothe and Fleckenstein (2013) used a piecewise linear interface
calculation (PLIC) algorithm for the transport of species at the interface to
prevent artificial mass transfer across the interface. Hassanvand and
Hashemabadi (2012) demonstrated that with no special treatment for the scalar-
transport equation, the scalar quantity artificially diffuses across the
interface into the other phase, even though the diffusivity of the scalar in
the other phase is set to zero. They proposed a fix to the issue by modifying
the flux of the scalar and by setting it to zero at the interface. Berry et
al. (2013) combined the zero-flux boundary condition for charged ions with the
ion-transport equation and derived a modified transport equation that
implicitly achieves the ion-impenetrable boundary condition at the interface.

In this study, however, we use a diffuse-interface method for accurate
modeling of the interface between two fluids. Advantages of a diffuse-
interface method over a geometric VOF method is the absence of geometric
reconstruction of the interface that results in a low-cost, load-balanced, and
a highly-scalable method, see Jain et al. (2020) for the parallel-scalability
performance of a diffuse-interface method on large scale computing machines.
However, a VOF method is known to be more accurate compared to a diffuse-
interface method on grids of same size. Therefore, on a per-cost basis, both
VOF and diffuse-interface methods are known to perform similarly for
incompressible two-phase flows (Mirjalili et al., 2019). Compared to a level-
set method, a diffuse-interface method also maintains discrete conservation of
volume of each of the phases. For a more detailed discussion and a comparison
of various interface-capturing methods, see the recent review by Mirjalili et
al. (2017).

Unlike the numerous advances in modeling the transport of scalars in the
context of a VOF method, there have been very few studies for scalar transport
in a diffuse-interface method. Scalar equation algorithm (SEA) is a method
that is similar to a diffuse-interface method, wherein the transport of volume
fraction is done using the finite-difference schemes with no geometric
reconstruction. Pericleous et al. (1995) used the SEA method for modeling
free-surface flows and adopted the Van Leer total-variation diminishing (TVD)
scheme for the discretization of the energy equation instead of an upwinding
scheme to avoid the artificial heat loss across the interface due to numerical
smearing. However, the mere use of a Van Leer TVD scheme for the
discretization of the scalar-transport equation does not completely eliminate
the artificial diffusion of the scalar across the interface (Mehdi-Nejad et
al., 2004). Most VOF methods tackle the issue of artificial numerical
diffusion of the scalar across the interface by geometrically advecting the
scalar concentration field in a manner consistent with the advection of the
volume fraction field (Davidson and Rudman, 2002, Alke et al., 2009, Bothe and
Fleckenstein, 2013, Berry et al., 2013). However, in a diffuse-interface
method, advection of the volume fraction field is performed non-geometrically
with the use of finite-difference schemes, which is a main contributing factor
in reducing the cost and improving the scalability of the method compared to a
VOF method. But the lack of geometric fluxing in diffuse-interface method
poses a challenge towards modeling scalars and to prevent the artificial
numerical diffusion (numerical leakage) of the scalar across the interface.

To the best of our knowledge, still lacking are numerical methods with
simulation capabilities that can accurately capture the behavior of the scalar
concentration fields without any numerical leakage from one phase to the other
while still maintaining the positivity of the scalar concentration values. For
example, in the work of Davidson and Rudman (2002), negative values of the
scalar were reset to zero to prevent the method from diverging. To address
these deficiencies, we have developed a general scalar-transport model for
two-phase flows, particularly for material interfaces modeled using a diffuse-
interface method. The proposed model is general enough that the scalar can
represent soluble surfactants, electrolytes, temperature, chemical species in
combustion modeling, etc. We adopt the two-scalar approach (Davidson and
Rudman, 2002, Alke et al., 2009, Bothe and Fleckenstein, 2013, Berry et al.,
2013, Ma and Bothe, 2013) where a separate scalar-transport equation is solved
for each of the phases but each of the scalar-transport equations are solved
in the entirety of the domain. Our newly developed scalar-transport model
prevents the artificial numerical diffusion of scalars from one phase to the
other, even in challenging flow environments such as in the presence of
electrokinetic effects and in turbulent flows. We discretize the proposed
model equation using a second-order central scheme in space due to its non-
dissipative nature, that is crucial for the simulation of turbulent flows
(Mittal and Moin, 1997). We have proved that the resulting discrete equation
maintains the positivity of the scalar concentration field throughout the
simulation, provided the grid resolution is sufficient to resolve the length
scales present in the flow. We have thus derived a “positivity criterion” that
determines this length scale that needs to be resolved to maintain the
positivity of the scalar concentration field for two-phase flows.

The rest of this paper is organized as follows. The conservative diffuse-
interface method, that is used for modeling two-phase flows in this work, is
presented in Section 2. The discrete representation of the scalar
concentration field adopted in this work is formally presented in Section 3
along with some of the previous scalar-transport models available in the
literature. The newly proposed scalar-transport model in this work is
presented in Section 4; followed by the presentation of the proof of
positivity for the scalar concentration field, and the definition of the
positivity criterion in terms of the grid resolution requirements, in Section
5. Finally, the adopted numerical strategy is summarized in Section 6; and
results from the numerical simulations are presented in Section 7,
illustrating the applicability of the proposed scalar-transport model for
modeling passive and active scalars in a wide range of two-phase flow settings
spanning laminar to turbulent flow regimes. The summary of results and
conclusions are presented in Section 8.

##  2 Conservative phase-field method

The first step towards simulating scalars in a two-phase flow is to choose an
interface-capturing method that can accurately simulate interfaces in
challenging flow environments. In this work, we choose a phase-field for the
reasons described in the previous section. Within the class of phase-field
methods, the popular approaches are based on the Cahn-Hilliard equation by
Cahn and Hilliard (1958) and the Allen-Cahn equation by Allen and Cahn (1979)
which were originally proposed to describe the phase separation and coarsening
phenomena in solids and the motion of antiphase boundaries in crystalline
solids, respectively. More recently, these equations have been successfully
used in fluid dynamics to model the transport of material interfaces between
two fluids (Anderson et al., 1998, Lowengrub and Truskinovsky, 1998, Chen et
al., 1998, Jacqmin, 1999, Liu and Shen, 2003, Badalassi et al., 2003, Yue et
al., 2004, Yang et al., 2006, Kim, 2012). However, there are some downsides to
the use of these equations. The Cahn-Hilliard equation conserves the volume of
each phase but contains a fourth-order derivative term which requires special
care for its discretization. The Allen-Cahn equation, on the other hand, does
not conserve the volume of each phase but involves a second-order derivative
term. To have the advantages of both equations, Chiu and Lin (2011) developed
the conservative phase-field method for incompressible two-phase flows by
combining the reinitialization step of the conservative level-set method by
Olsson and Kreiss (2005) with the phase-field equation.

In this work, we choose to use the conservative phase-field method for
incompressible two-phase flows. In this method the phase-field equation is
written as

$\frac{\partial\phi}{\partial
t}+\vec{\nabla}\cdot(\vec{u}\phi)=\vec{\nabla}\cdot\Big{[}\Gamma\Big{\\{}\epsilon\vec{\nabla}\phi-\phi(1-\phi)\vec{n}\Big{\\}}\Big{]},$
(1)

where $\phi$ is the phase-field variable which represents volume fraction of
the given phase,
$\Gamma\\{\epsilon\vec{\nabla}\phi-\phi(1-\phi)\vec{n}\\}=\vec{a}(\phi)$ term
on the right-hand side is the flux of the interface regularization (diffusion-
sharpening) term, $\vec{n}=\vec{\nabla}\phi/|\vec{\nabla}\phi|$ is the unit
normal vector to the interface, and $\Gamma$ and $\epsilon$ are the interface
parameters, where $\Gamma$ represents an artificial regularization velocity
scale and $\epsilon$ represents an interface thickness scale. This equation
satisfies both $\phi_{1}$ and $\phi_{2}$, where $\phi_{1}$ and
$\phi_{2}=1-\phi_{1}$ are the phase-field variables for phases $1$ and $2$,
respectively. We can show that the interface regularization term satisfies
$\vec{a}(\phi_{1})=-\vec{a}(\phi_{2})$ for phases $1$ and $2$. The phase-field
equation [Eq. (1)] can be solved in conjunction with the momentum balance
equation

$\frac{\partial\rho\vec{u}}{\partial
t}+\vec{\nabla}\cdot(\rho\vec{u}\otimes\vec{u}+p\mathds{1})=\vec{\nabla}\cdot\underline{\underline{\tau}}+\vec{\nabla}\cdot(\vec{S}\otimes\vec{u}),$
(2)

where

$\vec{S}=\Gamma\left[\epsilon\vec{\nabla}\rho-\frac{(\rho_{1}-\rho)(\rho-\rho_{2})}{\rho_{1}-\rho_{2}}\vec{n}\right]$
(3)

is the net mass regularization flux, $\rho_{l}$ is the density of the phase
$l$ and is a constant, and $\rho=\sum_{l=1}^{2}\rho_{l}\phi_{l}$ is the total
density of the mixture. Here, $\vec{S}$ is added to the incompressible
momentum equation so that the resulting momentum is transported in a manner
consistent with the transport of volume fraction using the phase-field
equation (Mirjalili and Mani, 2020). A similar consistency corrections for the
momentum equation have been proposed for the Cahn-Hilliard and Allen-Cahn
models in Huang et al. (2020a) and Huang et al. (2020b), respectively. The
Cauchy stress tensor is written as
$\underline{\underline{\tau}}=2\mu\mathbb{D}-2\mu(\vec{\nabla}\cdot\vec{u})\mathds{1}/3$,
where $\mu$ is the dynamic viscosity of the mixture evaluated using the one-
fluid mixture rule (Kataoka, 1986) as $\mu=\sum_{l=1}^{2}\phi_{l}\mu_{l}$,
$\mu_{l}$ is the dynamic viscosity of the phase $l$,
$\mathbb{D}=\\{(\vec{\nabla}\vec{u})+(\vec{\nabla}\vec{u})^{T}\\}/2$ is the
strain-rate tensor, and $p$ is the pressure. In this work, both the phases are
assumed to be incompressible. Therefore, a pressure field $p$ that satisfies
the criterion of zero divergence of the velocity field can be computed using
the fractional-step method of Kim and Moin (1985). However, this does not
limit the applicability of the newly proposed scalar-transport model in this
work for incompressible flows; it can also be used with the conservative
diffuse-interface method proposed by Jain et al. (2020) for compressible two-
phase flows because of the same form of the equilibrium kernel function—a
hyperbolic-tangent function—for the phase field $\phi$ (see, Section 4.1).

It can be shown that the phase field $\phi$ is bounded between $0$ and $1$
(Mirjalili et al., 2020), provided $\Gamma$ and $\epsilon$ are chosen such
that they satisfy the criterion

$\frac{\epsilon}{\Delta
x}\geq\frac{\left(\frac{|u|_{\mathrm{max}}}{\Gamma}+1\right)}{2},$ (4)

where $\Delta x$ is the grid size and $|u|_{\mathrm{max}}$ is the maximum
value of the magnitude of the velocity in the domain. This method also
inherently satisfies the TVD property for the phase field (Jain et al., 2020)
with the use of a central-difference scheme for the discretization of all the
operators and without having to add any flux limiters that destroy the non-
dissipative nature of the scheme.

##  3 Problem description

FIGURE 1: Schematic representing the domain filled with two phases and a
scalar quantity that is confined to phase $1$. Color represents the scalar
concentration.

Consider the schematic of a domain, $\Omega$, with two phases $1$ and $2$
shown in Figure 1. Let $c$ represent the concentration (amount of scalar per
unit volume) of a scalar quantity in the domain. The scalar is said to be
conserved, $\int_{\Omega}c\ dV=\mathrm{constant}$, if it is not being
generated or destroyed. As described in Section 1 the ratio of diffusivites of
the scalar in two phases is typically very large, hence the scalar will be
confined to one of the phases as illustrated in the Figure 1. Now
reformulating and generalizing the problem at hand, i.e., the diffusivity is,
say, a finite value $D$ in phase $1$ and practically zero in phase $2$, then
the scalar is confined only to phase $1$ in the domain.

Since the scalar is confined to phase $1$, we can then define another variable
$\tilde{c}$ that represents the local concentration of the scalar in phase $1$
(amount of scalar per unit volume of the phase $1$). Then, the relation
between $c$ and $\tilde{c}$ is

$c=\phi\tilde{c},$ (5)

where $\phi$ represents the volume fraction of the phase (volume of the phase
per unit total volume) where the scalar is present (phase $1$ in this case).
By definition, the value of $\tilde{c}$ is undefined in phase $2$. In the
continuum limit of an infinitely sharp interface, $\phi$ reduces to a
Heaviside function with a value of one inside phase $1$ and zero inside phase
$2$. Therefore, in this limit, the scalar concentration can be represented as

$c=\Bigg{\\{}\begin{aligned} &\tilde{c}\hskip 28.45274pt\mathrm{inside\ the\
phase\ where\ scalar\ is\ present}\\\ &0\hskip
28.45274pt\mathrm{elsewhere}.\end{aligned}$ (6)

Now, assuming that $\tilde{c}$ satisfies a generic transport (advection-
diffusion) equation within phase $1$, an evolution equation for $\tilde{c}$
can be written as

$\frac{\partial\tilde{c}}{\partial
t}+\vec{\nabla}\cdot(\vec{u}_{c}\tilde{c})=\vec{\nabla}\cdot(D\vec{\nabla}\tilde{c}),$
(7)

where $\vec{u}_{c}=\vec{u}+\vec{u}_{r}$ represents the total convective
velocity of the scalar, $\vec{u}$ represents the fluid velocity, $\vec{u}_{r}$
represents any effective velocity with which the scalar is being advected
relative to the fluid (e.g., electromigration velocity
$\vec{u}_{r}=\nu\vec{E}$, where $\nu$ is the electrical mobility and $\vec{E}$
is the electric field), and $D$ represents the diffusivity of the scalar.
Since the scalar is confined to phase $1$ with the interface acting as the
boundary, this represents a classical boundary value problem. The appropriate
no-flux boundary condition for this problem can be written as

$(D\nabla\tilde{c}-\vec{u}_{r}\tilde{c})\cdot\vec{n}=0.$ (8)

Combining Eq. (7) and Eq. (8), Berry et al. (2013) derived the modified
transport equation for charged ions that can implicitly achieve the ion-
impenetrable boundary condition as

$\frac{\partial(\phi\tilde{c})}{\partial
t}+\vec{\nabla}\cdot(\phi\vec{u}_{c}\tilde{c})=\vec{\nabla}\cdot(D\phi\vec{\nabla}\tilde{c}).$
(9)

They used this modified transport equation along with a sharp-interface method
(coupled level-set and volume-of-fluid method) to simulate the electrokinetics
of liquid-liquid systems.

In the discrete limit when the interface is not represented as a sharp
boundary—more so with the use of diffuse-interface methods—however, the
interface is not infinitely sharp; and the smallest interface thickness that
can be handled on an Eulerian grid is of the order of grid-cell size. In this
limit, Eq. (6) would not hold anymore, which is illustrated in the schematic
shown in Figure 2.

FIGURE 2: Schematic representing the discrete $c$ and $\tilde{c}$ quantities.
Opacity of the color is proportional to the field value in that particular
cell. Here, “conc.” represents the local scalar concentration, “vol. frac.” is
the local volume fraction of the phase and $cv$ represents a control volume.

In this discrete limit and with the use of finite-difference schemes, i.e.,
without any special geometric treatment, we can observe that $\tilde{c}$ is
not a conserved quantity, i.e., $\int_{\Omega}\tilde{c}\
dV\neq\mathrm{constant}$, because of the values that it takes in the cells
that contain the interface. This is the reason most VOF methods use geometric
fluxing for the transport of the scalar quantities. However, $\int_{\Omega}c\
dV=\mathrm{constant}$ still holds, and therefore, $c$ is still a conserved
quantity in this discrete limit. Hence, we could instead write a transport
equation for $c$ and look for modifications that mimic the no-flux boundary
condition at the interface. In other words, we seek a transport equation for
$c$ that results in transport of the scalar quantity consistent with the
transport of the phase field (see Figure 4; we define consistency more
rigorously in Section 4.1), such that there is no artificial numerical
diffusion of the scalar at the interface.

With this notion, we could start with a generic form of the transport equation
for $c$ as

$\frac{\partial c}{\partial
t}+\vec{\nabla}\cdot(\vec{u}_{c}c)=\vec{\nabla}\cdot(D\vec{\nabla}c).$ (10)

A straightforward modification (a naive approach) to this equation could be to
multiply $\phi$ to $\vec{u}_{r}$ and $D$, such that the flux
$D\vec{\nabla}c-\vec{u}_{r}c$ goes to zero as $\phi$ goes to zero outside
phase $1$ as

$\frac{\partial c}{\partial
t}+\vec{\nabla}\cdot(\vec{u}c+\phi\vec{u}_{r}c)=\vec{\nabla}\cdot(D\phi\vec{\nabla}c).$
(11)

Though this simple modification to achieve no-flux boundary condition sounds
promising, the scalar concentration field $c$ that we obtain by solving Eq.
(11) will not be consistent with the phase field$\phi$ and it results in the
artificial numerical diffusion across the interface, especially for diffuse-
interface methods. However, this approach has been used previously with a VOF
method in the context of modeling heat and mass transfer across interfaces
(Davidson and Rudman, 2002).

Figure 3 shows the application of Eq. 11 for the evolution of the scalar
quantity in a stationary one-dimensional drop. Since the initial concentration
$c$ of the scalar was uniform everywhere within the drop, the initial state of
$c$ should be maintained throughout the simulation. However, the model in Eq.
11 results in predicting a lower value of $c$ with an $\mathcal{O}(1)$ error
due to the numerical leakage of the scalar across the interface as illustrated
in Figure 3 (b), thus affecting the overall accuracy of the solution. For
problems with a two-way coupling between the flow and the transport of scalar,
the local concentration of the scalar field is crucial in predicting accurate
flow fields. For example, the transport of reacting species in combustion
modeling and ion transport in electrokinetics both modify the flow field
depending on the local concentration of the scalar being transported.
Therefore, the spurious leakage of the scalar and the resulting
$\mathcal{O}(1)$ error in these cases, could be detrimental to the overall
accuracy of the simulation. In the worst-case scenario, the numerical solution
for the scalar concentration field can admit unphysical negative values near
the interface and often within the phase with the lower diffusivity value,
which might result in unrealizable scalar concentration fields.

FIGURE 3: Simulation results obtained by solving Eq. (11), showing the one-
dimensional drop ($\phi$ field) and the scalar concentration ($c$ field) at
(a) initial time and (b) later times. Parameters used for this simulation are
$D=0.01$, $\Delta x=0.01$, $u=0$, and $u_{r}=0$. The arrow denotes the time
evolution of the scalar concentration field.

One trivial approach to address this issue of unphysical negative values is to
reset them to zero every time step (Davidson and Rudman, 2002). However, this
would destroy the conservative property of the scalar quantity. Other common
and more sophisticated ways to address this issue is to use flux limiters and
positivity-preserving limiters (Laney, 1998). However, the numerical
diffusivity associated with these schemes inevitably adds to the unphysical
leakage of the scalar across the interface into the impermeable phase. Hence,
in the current work, we present a new model equation for the transport of
scalars in two-phase flows along with a consistent numerical discretization
scheme that overcomes the challenges that were presented here.

##  4 Proposed model equation for the transport of scalars in two-phase flows

With the objective of developing a scalar-transport model that does not admit
negative values for the scalar concentration field, nor permit unphysical
leakage of the scalar across the interface, we propose a scalar-transport
model for the scalar quantity that is confined to phase $1$ as

$\frac{\partial c}{\partial
t}+\vec{\nabla}\cdot(\vec{u}c+\phi\vec{u}_{r}c)=\vec{\nabla}\cdot\left[D\left\\{\vec{\nabla}c-\frac{(1-\phi)\vec{n}c}{\epsilon}\right\\}\right],$
(12)

where $\vec{n}$ is the normal vector to the interface, $\epsilon$ is the same
interface parameter that is present in Eq. (1). The proposed model equation in
Eq. (12) is a modification to Eq. (11), which failed to prevent the numerical
leakage of the scalar and maintain the no-flux boundary condition for the
scalar at the interface (Figure 3). Away from the interface, in the bulk
region of phase $1$ where the scalar is present ($\phi\rightarrow 1$), the
proposed model in Eq. (12) reduces to

$\frac{\partial c}{\partial
t}+\vec{\nabla}\cdot(\vec{u}c+\vec{u}_{r}c)=\vec{\nabla}\cdot\big{(}D\vec{\nabla}c\big{)},$
(13)

which is the transport equation for scalars in a single-phase flow. Hence, the
model has no adverse effect on the transport of scalars away from the
interface.

Since we adopt the two-scalar approach, wherein a separate scalar equation is
written for each of the phases, an equation for the scalar in phase $2$,
represented by the phase field $1-\phi$, can be written by replacing $\phi$
with $1-\phi$ in Eq. 12 as

$\frac{\partial c^{\prime}}{\partial
t}+\vec{\nabla}\cdot\\{\vec{u}c^{\prime}+(1-\phi)\vec{u}_{r}c^{\prime}\\}=\vec{\nabla}\cdot\left\\{D\left(\vec{\nabla}c+\frac{\phi\vec{n}c^{\prime}}{\epsilon}\right)\right\\},$
(14)

where $c^{\prime}$ represents the concentration of the scalar quantity that is
confined to phase 2. Furthermore, the transfer of scalar across the interface
can also be straightforwardly accounted for in this two-scalar approach by
introducing sink/source terms (Davidson and Rudman, 2002) into the scalar-
transport model in Eqs. (12),(14) and will be explored in a future work.

###  4.1 Consistency and equilibrium solution

As described in Section 3, the consistency of the scalar concentration $c$ and
the phase field $\phi$ is crucial in preventing the unphysical numerical
leakage of the scalar across the interface. The $c$ and $\phi$ are formally
said to be consistent if they possess the same equilibrium kernel function—a
hyperbolic tangent function—such that the local concentration field
$\tilde{c}$, which is defined as the ratio $\tilde{c}=c/\phi$, stays constant
across the interface. The local concentration field $\tilde{c}$ then
represents the interfacial concentration of the scalar quantity. The proposed
scalar-transport model in Eq. (12) was derived such that $c$ satisfies this
consistency condition.

To realize the above defined consistency between $c$ and $\phi$, we can look
at the steady-state equilibrium solutions to the phase-field equation and the
proposed scalar-transport equation. At steady state, for $\vec{u}=0$ and
$\vec{u}_{r}=0$ and in one dimension, the phase-field equation in Eq. (1)
reduces to the form

$0=\vec{\nabla}\cdot\left[\Gamma\epsilon\left\\{\vec{\nabla}\phi-\frac{(1-\phi)\vec{n}\phi}{\epsilon}\right\\}\right]\Rightarrow\frac{d^{2}\phi}{dx^{2}}-\frac{1}{\epsilon}\frac{d\left\\{(1-\phi)\phi\right\\}}{dx}=0,$
(15)

for $\vec{n}=+1$, and the proposed scalar-transport equation in Eq. (12)
reduces to the form

$0=\vec{\nabla}\cdot\left[D\left\\{\vec{\nabla}c-\frac{(1-\phi)\vec{n}c}{\epsilon}\right\\}\right]\Rightarrow\frac{d^{2}c}{dx^{2}}-\frac{1}{\epsilon}\frac{d\left\\{(1-\phi)c\right\\}}{dx}=0.$
(16)

FIGURE 4: Equilibrium solutions for $\phi$ and $c$, showing consistency
between the scalar concentration field and the volume fraction field,
therefore, a constant value of the local concentration $\tilde{c}=c/\phi$.
Here, $c_{0}$ is chosen to be equal to $5$ for the sake of illustration.

Now, assuming that the interface is at the origin as shown in Figure 4, and
integrating Eq. 15 along with the boundary conditions

$\phi=\Bigg{\\{}\begin{aligned} &0\hskip 28.45274ptx\rightarrow-\infty\\\
&0.5\hskip 28.45274ptx=0,\end{aligned}$ (17)

we obtain

$\phi=\frac{e^{(x/\epsilon)}}{1+e^{(x/\epsilon)}}=\frac{1}{2}\left\\{1+\tanh{\left(\frac{x}{2\epsilon}\right)}\right\\}.$
(18)

Using the equilibrium solution for $\phi$ in Eq. (18) and solving for $c$ by
integrating the Eq. (16) and using the boundary conditions

$c=\Bigg{\\{}\begin{aligned} &0\hskip 28.45274ptx\rightarrow-\infty\\\
&c_{0}\hskip 28.45274ptx=0,\end{aligned}$ (19)

we obtain

$c=2c_{0}\frac{e^{(x/\epsilon)}}{1+e^{(x/\epsilon)}}=c_{0}\left\\{1+\tanh{\left(\frac{x}{2\epsilon}\right)}\right\\}.$
(20)

Hence, the equilibrium kernel functions of $c$ and $\phi$ are both hyperbolic
tangent functions of the spatial coordinate along the interface normal (here
$x$), and therefore, are consistent at the interface. This results in a
constant value of $\tilde{c}=c/\phi=2c_{0}$ across the interface (see Figure
2). If we choose $c_{0}$ to be equal to $0.5$, the equilibrium kernel
functions of $c$ and $\phi$ are identical, and the resulting interfacial
scalar concentration is $\tilde{c}=1$. For other values of $c_{0}$, $c$ is
only scaled by a constant factor, but is still a hyperbolic tangent function
as illustrated in Figure 4. Hence, the scalar concentration field obtained
from solving the proposed model in Eq. (12) is consistent with the phase field
from Eq. (1), which results in the transport of the scalar inside the confined
phase without any unphysical numerical leakage across the interface that would
be seen otherwise.

##  5 Positivity of scalars

Positivity of the scalar concentration field is a crucial realizability
requirement that needs to be satisfied at all times in the simulation. A
common approach to achieve this is to use flux limiters and positivity-
preserving limiters (Laney, 1998). However, these limiters add artificial
numerical dissipation to the scheme and also lead to unphysical numerical
leakage of the scalar across the interface. Hence, we use a central-difference
scheme to discretize the operators in our system of equations because of its
well-known non-dissipative property (Moin and Verzicco, 2016), and derive a
criterion [Eq. (21); Figure 5] for the choice of grid size to be used in the
simulation. We propose in Theorem 5.1 that the criterion in Eq. (21) is a
sufficient condition for the proposed model transport equation in Eq. (12) to
maintain the positivity of the scalar concentration field $c$ at all times
during the simulation.

######  Theorem 5.1.

If $\phi^{k}_{i}$ is bounded between $0$ and $1$, $\forall k\in\mathds{Z}^{+}$
and $\forall i$, on a uniform one-dimensional grid, then $c^{k}_{i}\geq 0$
holds $\forall k\in\mathds{Z}^{+}$, where $k$ is the time-step index and $i$
is the grid index, provided

$\Delta
x\leq\left(\frac{2D}{|u|_{\mathrm{max}}+|u_{r}|_{\mathrm{max}}+\frac{D}{\epsilon}}\right),$
(21)

and

$\Delta t\leq\frac{\Delta x^{2}}{2D}$ (22)

are satisfied, where $\Delta x$ is the grid-cell size, $\Delta t$ is the time-
step size, $|u|_{\mathrm{max}}$ and $|u_{r}|_{\mathrm{max}}$ are the maximum
fluid velocity and the maximum relative velocity of the scalar in the domain,
respectively.

###### Proof.

Consider the discretization of Eq. (12) on a one-dimensional uniform grid

$\displaystyle c^{k+1}_{i}=c^{k}_{i}+\Delta
t\left[-\left(\frac{u^{k}_{i+1}c^{k}_{i+1}-u^{k}_{i-1}c^{k}_{i-1}}{2\Delta
x}\right)-\left(\frac{{u^{k}_{r}}_{i+1}\phi^{k}_{i+1}c^{k}_{i+1}-{u^{k}_{r}}_{i-1}\phi^{k}_{i-1}c^{k}_{i-1}}{2\Delta
x}\right)\right]$ (23) $\displaystyle+\Delta
t\left[D\left(\frac{c_{i+1}^{k}-2c_{i}^{k}+c_{i-1}^{k}}{\Delta
x^{2}}\right)-\frac{D}{\epsilon}\left\\{\frac{(1-\phi_{i+1}^{k})n_{i+1}^{k}c_{i+1}^{k}-(1-\phi_{i-1}^{k})n_{i-1}^{k}c_{i-1}^{k}}{2\Delta
x}\right\\}\right].$

This can be rearranged as

$c_{i}^{k+1}=\tilde{C}_{i-1}^{k}c_{i-1}^{k}+\tilde{C}_{i}^{k}c_{i}^{k}+\tilde{C}_{i+1}^{k}c_{i+1}^{k},$
(24)

where $\tilde{C}$’s are coefficients given by

$\tilde{C}_{i-1}^{k}=\frac{\Delta tu^{k}_{i-1}}{2\Delta x}+\frac{\Delta
t{u^{k}_{r}}_{i-1}\phi^{k}_{i-1}}{2\Delta x}+\frac{\Delta tD}{\Delta
x^{2}}+\frac{\Delta tD}{2\epsilon\Delta x}(1-\phi_{i-1}^{k})n_{i-1}^{k},$ (25)
$\tilde{C}_{i+1}^{k}=-\frac{\Delta tu^{k}_{i+1}}{2\Delta x}-\frac{\Delta
t{u^{k}_{r}}_{i+1}\phi^{k}_{i+1}}{2\Delta x}+\frac{\Delta tD}{\Delta
x^{2}}-\frac{\Delta tD}{2\epsilon\Delta x}(1-\phi_{i+1}^{k})n_{i+1}^{k}$ (26)

and

$\tilde{C}^{k}_{i}=1-\frac{2\Delta tD}{\Delta x^{2}}.$ (27)

######  Lemma 5.1.1.

A scheme is said to maintain positivity (also called the “boundedness”
criterion in Patankar (1980), Versteeg and Malalasekera (2007)) if
$\tilde{C}$’s are all positive (Laney, 1998).

It is given that $0\leq\phi^{k}_{i}\leq 1$ holds $\forall k\in\mathds{Z}^{+}$,
which implies that $(1-\phi_{i-1}^{k})n_{i-1}^{k}\geq-1$. Using this in Eq.
(25), we obtain

$\displaystyle\tilde{C}_{i-1}^{k}\geq\frac{\Delta tu^{k}_{i-1}}{2\Delta
x}+\frac{\Delta t{u^{k}_{r}}_{i-1}\phi^{0}_{i-1}}{2\Delta x}+\frac{\Delta
tD}{\Delta x^{2}}-\frac{\Delta tD}{2\epsilon\Delta x}$ (28)
$\displaystyle\Rightarrow\tilde{C}_{i-1}^{k}\geq-\frac{\Delta t}{2\Delta
x}\bigg{(}|u|^{k}_{max}+|u_{r}|^{k}_{max}+\frac{D}{\epsilon}\bigg{)}+\frac{\Delta
tD}{\Delta x^{2}}$

Now, invoking the condition in Eq. (21), we can show that
$\tilde{C}_{i-1}^{k}\geq 0$ holds. Using similar arguments, we can show that
$\tilde{C}_{i+1}^{k}\geq 0$ holds. Invoking the condition in Eq. (22), we can
also show that $\tilde{C}_{i}^{k}\geq 0$ holds. Thus, Lemma 5.1.1 proves that
$c^{k}_{i}\geq 0$ is satisfied $\forall k\in\mathds{Z}^{+}$, which concludes
the proof.

∎

Now, generalizing Theorem 5.1 for three dimensions, the time-step restriction
required for the positivity of $c$—assuming an isotropic mesh—can be written
as

$\Delta t\leq\frac{\Delta x^{2}}{6D}.$ (29)

Note that $\phi$ is assumed to be bounded between $0$ and $1$ throughout the
simulation. This boundedness of $\phi$ can be achieved by choosing appropriate
values for the interface parameters $\Gamma$ and $\epsilon$ such that they
satisfy the criterion in Eq. (4). For more details on the boundedness of
$\phi$, see Mirjalili et al. (2020) for incompressible flows and Jain et al.
(2020) for compressible flows.

The criterion on the time-step size in Eq. (22) and Eq. (29) are the Courant-
Friedrich-Levy (CFL) conditions for the scalar diffusion process in one and
three dimensions, respectively, and are typically already satisfied in an
explicit time-marching scheme to achieve temporal stability. Hence, the only
additional criterion that needs to be satisfied to maintain the positivity of
the evolution of the scalar concentration field is the restriction on the grid
size given in Eq. (21). In the proof of Theorem 5.1, a first-order Euler time-
stepping scheme was used to arrive at the restrictions on the time-step size
in Eq. (22) and Eq. (29); however, these criteria are sufficient to maintain
the positivity of $c$ with most higher-order explicit time-stepping schemes
since the diffusive CFL conditions for the scalar in Eq. (22) and Eq. (29) are
less restrictive for higher-order time-stepping schemes.

It was shown in Section 4 that the proposed model transport equation in Eq.
(12) reduces to a generic scalar-transport equation [Eq. (13)] for single-
phase flows, in the bulk region away from the interface. Now, repeating the
analysis in Theorem 5.1 for the scalar-transport equation for the bulk region
[Eq. (13)], we can show that the criterion that needs to be satisfied to
maintain the positivity of evolution of the scalar concentration field is the
restriction on the grid size given by

$\Delta
x\leq\bigg{(}\frac{2D}{|u|_{\mathrm{max}}+|u_{r}|_{\mathrm{max}}}\bigg{)}.$
(30)

A similar analysis holds good if one is interested in the large-eddy
simulation (LES) of scalars in a two-phase turbulent flow, instead of a direct
numerical simulation (DNS). Here, the diffusivity $D$ of the scalar can be
replaced by the effective diffusivity (sum of resolved and subgrid
contributions) of the scalar in Eq. (21) and Eq. (30).

###  5.1 Spatial resolution requirements

FIGURE 5: Graphical representation of the positivity criteria. The solid line
is the positivity criterion in Eq. (31) for the proposed scalar-transport
model [Eq. (12)], and the dashed line is the positivity criterion in Eq. (32)
for the scalar-transport model in the bulk region away from the interface [Eq.
(13)]. The red dots represent the positivity criteria for the choice of
$\epsilon=\Delta x$ and the arrow represents the additional restriction
imposed on the grid size at the interface compared to the bulk region away
from the interface.

Now, rewriting the positivity criterion for the proposed scalar-transport
model in Eq. (21) in terms of the non-dimensional groups $[\epsilon/\Delta x]$
and $[D/\\{(u+u_{r})\epsilon\\}]$ as

$\Big{[}\frac{\epsilon}{\Delta
x}\Big{]}\geq\frac{1+\Big{[}\frac{D}{(u+u_{r})\epsilon}\Big{]}}{2\Big{[}\frac{D}{(u+u_{r})\epsilon}\Big{]}}$
(31)

and the positivity criterion in Eq. (30) for the bulk region away from the
interface as

$\Big{[}\frac{\epsilon}{\Delta
x}\Big{]}\geq\frac{1}{2\Big{[}\frac{D}{(u+u_{r})\epsilon}\Big{]}},$ (32)

where $u=|u|_{max}$ and $u_{r}=|u_{r}|_{max}$, symbol $[\cdot]$ represents a
non-dimensional group. By plotting $[\epsilon/\Delta x]$ vs
$[D/\\{(u+u_{r})\epsilon\\}]$ in Figure 5, we can see that any grid size
$\Delta x$ for a given $\epsilon$ and $D$ that lies above the lines maintains
the positivity of the scalar concentration field throughout the simulation.
Since the solid line is above the dashed line in Figure 5 for all values of
$[D/\\{(u+u_{r})\epsilon\\}]$, the proposed model transport equation in the
full form in Eq. (12) imposes a more restrictive condition on the grid size
$\Delta x$ for a given $\epsilon$ and $D$ than the reduced form of the
transport equation for the bulk region in Eq. (13). Therefore, the criterion
in Eq. (21) should be used to choose the grid size throughout the domain.

In the conservative phase-field method, $\epsilon$ is typically chosen to be
equal to $\Delta x$. With this choice of $\epsilon$, and recognizing that the
non-dimensional group $[D/\\{(u+u_{r})\epsilon\\}]=[D/\\{(u+u_{r})\Delta
x\\}]$ is the inverse of cell-Peclet number $Pe_{c}$, the positivity criterion
for the proposed model in Eq. (21) can be written in terms of cell-Peclet
number as

$1\geq\frac{1+\Big{[}\frac{1}{Pe_{c}}\Big{]}}{2\Big{[}\frac{1}{Pe_{c}}\Big{]}}\Rightarrow
Pe_{c}\leq 1,$ (33)

and the positivity criterion in Eq. (30) for the bulk region away from the
interface can be written as

$1\geq\frac{1}{2\Big{[}\frac{1}{Pe_{c}}\Big{]}}\Rightarrow Pe_{c}\leq 2.$ (34)

This criterion in Eq. (34) was first proposed by Patankar (1980) and Versteeg
and Malalasekera (2007) for a generic scalar-transport equation for single-
phase flows. Hence, the presence of a material interface and the effective no-
flux boundary condition for the scalar that prevents it from artificially
diffusing across the interface has introduced a more restrictive criterion on
the grid size in terms of the cell-Peclet number as $Pe_{c}\leq 1$, which
otherwise would have been $Pe_{c}\leq 2$ in the absence of any interface. This
is also graphically shown (arrow) in Figure 5. Therefore, assuming a uniform
grid throughout the domain, the grid size for the simulation of scalars with
material interfaces (two-phase flow) should be twice as small compared to the
grid size for the simulation of scalars in the absence of material interfaces
(single-phase flow). One could also use a twice-refined grid only around the
interface using an adaptive-mesh refinement (AMR) technique instead of using a
uniform grid throughout the domain, to reduce the cost of computation.

The restriction on the grid size in Eq. (33) can also be motivated by the
physical scales involved in the problem. A cell-Peclet number of $Pe_{c}\leq
1$, implies that the grid size is $\Delta x\leq D/(u+u_{r})$. Here,
$D/(u+u_{r})\sim l_{c}$ represents the characteristic length scale $l_{c}$ of
the problem. Therefore the cell-Peclet number restriction is alluding to the
fact that the grid size should be small enough to resolve the physical length
scales present in the problem, i.e., $\Delta x\lesssim l_{c}$.

##  6 Numerical strategy

The scalar-transport model proposed in this work is implemented in the CTR-
DIs2D and CTR-DIs3D solvers (Jain et al., 2018, 2020). These solvers can
handle both compressible and incompressible flows. For incompressible flows, a
finite-volume discretization strategy on a staggered grid has been employed
wherein the phase field, the pressure field, and the scalar concentration
fields are stored at the cell centers; and the components of the velocity
field vector are stored at the cell faces where all the fluxes are evaluated.
This choice of discretization is adopted, for incompressible flows, to avoid
the spurious checkerboarding of the pressure field (Patankar, 1980). The
pressure-Poisson equation is solved with a geometric-multigrid preconditioned
conjugate gradient method using the HYPRE package (Falgout and Yang, 2002).

We use the fourth-order Runge-Kutta (RK4) time-stepping scheme and the second-
order central-differencing scheme for the discretization of the spatial
operators. This choice of numerical scheme has some advantages, particularly
for the simulation of turbulent flows due to its (a) non-dissipative nature,
(b) low aliasing error, (c) easy boundary treatment, (d) low cost, and (e)
improved stability (Moin and Verzicco, 2016). With the appropriate choice of
$\Delta x$, $\Gamma$, and $\epsilon$, we can achieve the positivity for the
scalar concentration field (see, Section 5) and the boundedness and TVD
properties for the phase field (see, Section 2) even with the use of a
central-difference scheme for the spatial operators, which would otherwise
admit oscillatory solutions due to the associated dispersion errors.

##  7 Simulation results

The proposed scalar-transport model can be used with a wide range of two-phase
flows, from laminar to turbulent flow regimes. To illustrate this, multiple
test cases are presented in this section, starting from the simple one-
dimensional cases of a droplet advection with a scalar in Section 7.1. These
cases were used to assess the validity of the positivity criterion of the
proposed model by testing the model for various choices of parameters on the
positivity map in Figure 5. This is followed by the two-dimensional cases of a
bubble in a concentration-polarization region in Section 7.2 which involves
modeling the transport of a scalar quantity around the bubble, and charged
ions in a drop in Section 7.3. To simulate charged ions, the proposed scalar-
transport model was recast into a Nernst-Planck equation and was solved in
conjunction with the Gauss’s law. Finally, a three-dimensional case of a
scalar in a turbulent two-phase flow is presented in Section 7.4. In all the
above cases, the scalar is confined to one of the phases, and therefore, the
ability of the model to prevent the unphysical numerical leakage of the scalar
across the interface is also evaluated.

###  7.1 Advection of a droplet along with a scalar

In this section, a one-dimensional drop of radius $R=0.25$ is advected in a
periodic domain along with a scalar quantity that is confined to the drop, to
test the validity of the positivity criterion in Eq. (33). Four different
choices of the parameters are made as shown in Figure 6(a) in terms of the
cell-Peclet number ($Pe_{c}=(u+u_{r})\Delta x/D$), such that the two out of
the four choices ($Pe_{c}=1,0.8$) satisfy the positivity criterion and the
other two ($Pe_{c}=2,4$) violate the criterion. This test case is repeated
with a zero and a non-zero relative velocity of the scalar with respect to the
fluid, and with a non-uniform initial conditions for the scalar concentration
field as further described in the subsequent Sections 7.1.1, 7.1.2, and 7.1.3,
respectively. For all the simulations in this section, the drop is initially
placed in the domain centered at $x_{0}=0.5$ along with the scalar, and then
advected with a uniform prescribed velocity. The scalar has a uniform
concentration of $1$ inside the drop and $0$ outside, as shown in Figure 6(b),
unless specified otherwise. The domain is discretized into a uniform grid of
size $\Delta x=0.01$. The interface parameters are chosen as $\Gamma=100$ and
$\epsilon=\Delta x$.

FIGURE 6: (a) Four different choices of the parameters in terms of the cell-
Peclet number $Pe_{c}$, represented by ovals on the positivity plot, that are
selected to test the validity of the positivity criterion. (b) Initial
conditions for the volume fraction field $\phi$ and the scalar concentration
field $c$.

####  7.1.1 Scalar advection along with the drop ($\vec{u}_{r}=0$)

Here in this section, the uniform fluid velocity is chosen to be
$\vec{u}=100$, the relative velocity of the scalar with respect to the fluid
is $\vec{u}_{r}=0$, the domain range is $[0,1]$, and the total integration
time is $t_{\mathrm{end}}=10$. The diffusivity of the scalar is chosen based
on the required cell-Peclet number. The drop advects to the right due to the
imposed fluid velocity and returns to its original position at $t=0.01$ due to
the periodic boundary condition, and this process repeats $1000$ times until
the time $t=10$. The final state of the drop and the scalar concentration
field at time $t=10$ is shown in Figure 7 along with the minimum value of the
scalar concentration field for all four choices of the parameters shown in
Figure 6(a). For the cases with $Pe_{c}=1$ and $0.8$, that satisfy the
positivity criterion in Eq. (33), the scalar concentration $c$ is positive
throughout the domain; and for the cases with $Pe_{c}=2$ and $4$, that violate
the positivity criterion, the scalar concentration $c$ admits negative values
close to the interface, as expected.

FIGURE 7: Final state of the drop and the scalar concentration field at time
$t=10$ for the case of scalar advection along with the drop. Four plots
represent the four different parameters chosen to test the positivity of the
scalar: (a) $Pe_{c}=1$, (b) $Pe_{c}=0.8$, (c) $Pe_{c}=2$, and (d) $Pe_{c}=4$.

Comparing the results in Figure 7 with those in Figure 3, where the scalar was
found to artificially leak outside the drop resulting in an $\mathcal{O}(1)$
error for the scalar concentration values, it is easy to see the role and
importance of the proposed scalar-transport model in maintaining the
consistency between $c$ and $\phi$ and in the prevention of any artificial
numerical leakage of the scalar. To further quantify the error and
inconsistency between $c$ and $\phi$, $c-\phi$ and the local concentration
$\tilde{c}$ are plotted in Figure 8 for the same four choices of the
parameters shown in Figure 6(a). Since $c$ is expected to be same as $\phi$ in
this case, we could use $c-\phi$ as a metric to evaluate how far the solution
deviates compared to the expected one. As can be seen in Figure 8, the
solution is exact for $Pe_{c}=1$, but there is a small error for $Pe_{c}=0.8$
case. The error also seems to increase with the increasing value of $Pe_{c}$
for $Pe_{c}>1$. Small non-zero values of $c-\phi$ close to the interface imply
that the scalar concentration field is out-of-equilibrium by a small amount. A
similar behavior to $c-\phi$ can be seen in the plots of $\tilde{c}$. Here,
$\tilde{c}$ is computed as

$\tilde{c}=\Bigg{\\{}\begin{aligned} &\frac{c}{\phi}&\hskip
28.45274pt\phi>0.01\\\ &0&\hskip 28.45274pt\mathrm{else},\end{aligned}$ (35)

to avoid division by a small number.

FIGURE 8: $c-\phi$ and $\tilde{c}$ at the final time of $t=10$ for the case of
scalar advection along with the drop. Four plots represent the four different
parameters: (a) $Pe_{c}=1$, (b) $Pe_{c}=0.8$, (c) $Pe_{c}=2$, and (d)
$Pe_{c}=4$.

####  7.1.2 Scalar advection relative to the drop ($\vec{u}_{r}\neq 0$)

Here, the uniform fluid velocity is chosen to be $\vec{u}=50$, the relative
velocity of the scalar with respect to the fluid is $\vec{u}_{r}=50$, the
domain range is $[0,1]$, and the total integration time is
$t_{\mathrm{end}}=10$. The diffusivity of the scalar is chosen based on the
required cell-Peclet number. The drop advects to the right due to the imposed
fluid velocity and returns to its original position at $t=0.02$ due to the
periodic boundary condition, and this process repeats $500$ times until the
end of the simulation. The scalar also advects to the right, relative to the
drop, and accumulates on the right end of the drop since it is confined to the
drop. It reaches a steady state when the diffusion balances the advection due
to the relative velocity. The final state of the drop and the scalar
concentration field at time $t=10$ are shown in Figure 9 along with the
minimum values of the scalar concentration field for all four choices of the
parameters shown in Figure 6. For the cases with $Pe_{c}=1$ and $0.8$, that
satisfy the positivity criterion, the scalar concentration $c$ is positive
throughout the domain; and for the case with $Pe_{c}=4$, that violate the
positivity criterion, the scalar concentration $c$ admits negative values for
$Pe_{c}=4$, as expected. However the positivity is still maintained for
$Pe_{c}=2$ though the criterion is violated. This is because the criterion is
only a sufficient, not a necessary condition for the scalar concentration
field to remain positive and therefore presents the most restrictive condition
such that the positivity of the scalar is satisfied even in some situations
when the criterion is violated. This shows the robustness of the positivity
criterion that is crucial in maintaining the realizable values of the scalar
concentration field throughout the duration of the simulation.

FIGURE 9: Final state of the drop and the scalar concentration field at time
$t=10$ for the case of scalar advection relative to the drop. Four plots
represent the four different parameters chosen to test the positivity of the
scalar: (a) $Pe_{c}=1$, (b) $Pe_{c}=0.8$, (c) $Pe_{c}=2$, and (d) $Pe_{c}=4$
with $u_{r}\neq 0$.

####  7.1.3 Non-uniform scalar concentration within the drop

The simulations in the previous Sections 7.1.1 and 7.1.2 had uniform initial
scalar concentration fields within the drop. However, in reality, the scalar
concentration could be non-uniform inside the drop. To evaluate the model for
this scenario, we repeat the same test case in Section 7.1.2 with $Pe_{c}=1$,
except that the total integration time is $t_{\mathrm{end}}=0.05$ and a larger
domain of range $[0,5]$ is chosen with periodic boundary conditions. The drop
is placed at the same location of $x_{0}=0.5$ with a non-uniform scalar
concentration field given by the function with a compact support

$c=\Bigg{\\{}\begin{aligned}
&Ae^{\big{\\{}\frac{1}{(4x-2)^{2}-1}\big{\\}}}&\hskip
28.45274ptx\in[0.25,0.75]\\\ &0&\hskip 28.45274pt\mathrm{else},\end{aligned}$
(36)

as shown in Figure 10(a). Here, $A$ is a factor chosen to make the quantities
$\int_{\Omega}c\ dV$ and $\int_{\Omega}\phi\ dV$ discretely equal, such that
as the scalar diffuses and reaches a steady state within the drop, we should
expect to see the scalar concentration values reach a uniform value of $1$
within the drop. Initial conditions for $\phi$ and $c$, the time evolution of
$c$, and the final state are shown in Figure 10(a,b,c), respectively.

FIGURE 10: The advection of a drop along with a non-uniformly distributed
scalar quantity inside the drop. (a) The initial conditions for $\phi$ and
$c$. (b) The time evolution of $c$. (c) The final states of $\phi$ and $c$.

The scalar undergoes diffusion within the drop without any numerical leakage
or encountering negative values, and at the final time, the scalar
concentration is uniform within the drop, as expected, which illustrates the
robustness of the proposed model.

###  7.2 Bubble in a concentration polarization region

In this section, a two-dimensional test case of transport of a passive scalar
around the bubble in a channel is presented. A schematic of the domain is
shown in Figure 11. The domain is a square channel of size $L\times L$, where
$L=0.1$. The streamwise direction has periodic boundary conditions and the
other direction has no-slip walls on both ends. A gas bubble of radius
$R=d/2=0.02$ is initially placed at the center of the channel, surrounded by a
liquid. Because of the periodic boundary conditions, the single bubble in a
channel setup essentially represents a train of equally-spaced bubbles in a
channel.

FIGURE 11: A schematic of the domain for the case of scalar diffusion around a
bubble in a channel. Here, $c_{U}$ and $c_{L}$ represents the boundary
conditions for the scalar concentration on the upper and lower walls,
respectively.

Here, the scalar is confined to the liquid region. If the scalar quantity
represents the temperature field, then this corresponds to the physical
scenario of $Pr_{l}/Pr_{g}\rightarrow 0$, where $Pr$ is the Prandtl number,
and the subscripts $l$ and $g$ represents liquid and gas fluid quantities. If
the scalar instead represents dissolved salts or dye, then this corresponds to
the scenario of $Sc_{l}/Sc_{g}\rightarrow 0$, where $Sc$ is the Schmidt
number. In either case, the ratio of diffusivities is $D_{g}/D_{l}\rightarrow
0$, which essentially leads to the confinement of the scalar in the liquid
region.

In the current setup, the diffusivities are chosen to be $D_{l}=0.01$ and
$D_{g}=0$. The scalar concentration is initially zero everywhere in the
domain. The upper and lower walls have Dirichlet boundary conditions of
$c_{U}=1$ and $c_{L}=0$, respectively. Without the presence of the bubble in
the channel, the scalar concentration will remain uniform along the streamwise
direction. The final steady state scalar concentration profile along the wall-
normal direction can be obtained by solving a one-dimensional steady state
heat equation, which would give a linear profile with a slope of
$(c_{U}-c_{L})/L=10$. When the bubble is present, it experiences the varying
scalar concentration values around it which makes it a good test case to
evaluate the proposed model. Moreover, this test case was chosen because the
setup represents the concentration polarization region and the gas bubbles
that are formed at the electrode in a electrochemical system.

In the subsequent Sections 7.2.1 and 7.2.2, two situations that correspond to
(a) a stationary bubble, and (b) a moving bubble, have been studied. The
importance of the proposed scalar-transport model in Eq. 12 and the positivity
criterion in Eq. 33 is also highlighted by comparing the results obtained
using this model against the results obtained by solving the Eq. 11.
Hereafter, we denote the solution obtained from Eq. 11 as the “previous
approach”.

####  7.2.1 Stationary bubble

Here, the bubble is assumed to be stationary, and the fluid velocity is zero
everywhere in the domain. This setup, therefore, requires solving only the
scalar-transport equation. The time scale of diffusion of the scalar in the
domain is $\tau_{d}\sim L^{2}/D=1$, and we expect the scalar concentration to
reach a steady state after this time. Therefore, the total time of the
simulation is taken to be $t_{end}=2\tau_{d}=2$.

FIGURE 12: The final time steady state scalar concentration field obtained
using, (a) the previous approach on a grid of size $128\times 128$, (b) the
previous approach on a grid of size $512\times 512$, (c) the proposed model on
a grid of size $128\times 128$, and (d) the proposal model on a grid of size
$512\times 512$. The white solid line represents the interface; and the color
plot represents the scalar concentration field.

The final steady state scalar concentration fields obtained by solving the
proposed model and from the previous approach are shown in Figure 12, on two
different grids of sizes, $128\times 128$ and $512\times 512$. The previous
approach results in the leakage of the scalar into the bubble for both grids;
and increasing the grid resolution resulted in less leakage of the scalar.
However, the results from the proposed model show no signs of leakage of the
scalar into the bubble on either grids. The leakage of the scalar in the
previous approach is more evident in Figure 13, where $c$ and $\phi$ are
plotted along the wall-normal direction at $x=0$.

FIGURE 13: The scalar concentration field $c$ and the volume fraction field
$\phi$ at the final time of $t=2$ along the $x=0$ line for the case of scalar
diffusion around a stationary bubble in a channel, simulated on the grids of
size (a) $128\times 128$, and (b) $512\times 512$.

To further quantify the leakage seen in Figure 13, a quantity $c_{e}$ defined
as

$c_{e}=\Bigg{\\{}\begin{aligned} &|c-\phi|&\hskip 28.45274pt\phi<10^{-3}\\\
&0&\hskip 28.45274pt\mathrm{else},\end{aligned}$ (37)

is computed along the wall-normal direction at $x=0$ and is plotted in Figure
14. Since $c_{e}$ is the difference between $c$ and $\phi$ within the bubble
region, it approximately represents the effective amount of scalar that leaked
into the bubble region. An exact method would result in $c_{e}$ being zero,
which can be seen for the proposed method in Figure 14. Whereas, the previous
approach results in the leakage of the large amounts of scalar into the bubble
region.

FIGURE 14: The scalar leakage error $c_{e}$ at the final time of $t=2$ along
the $x=0$ line for the stationary-bubble case, simulated on the grids of size
(a) $128\times 128$, and (b) $512\times 512$.

Integrating $c_{e}$ along a line, an integrated leakage error metric
$I_{error}$ can be defined as

$I_{error}=\int_{s}c_{e}\ dxdy.$ (38)

Values of $I_{error}$ obtained by integrating $c_{e}$ in Figure 14 are listed
in Table 1. The integrated leakage error with the proposed method is four
orders of magnitude lower compared to the previous approach.

$I_{error}$ | $128\times 128$ | $512\times 512$
---|---|---
previous approach | $5.96\times 10^{-3}$ | $6.31\times 10^{-3}$
proposed method | $8.78\times 10^{-7}$ | $2.34\times 10^{-7}$
TABLE 1: Integrated leakage error for the stationary-bubble case.

####  7.2.2 Moving bubble

Here, in this section, two situations of moving bubble in a channel are
considered. In the first case, the bubble is assumed to be translating at a
fixed velocity without undergoing deformation. This setup, therefore, requires
solving the scalar-transport equation and the phase-field equation. Later, in
the second case, the bubble and the surrounding liquid is initialized with the
same velocity, but in this case the bubble is allowed to deform in an
accelerating flow. Therefore, this requires solving the proposed model coupled
with the hydrodynamics.

FIGURE 15: The scalar concentration field for the previous approach at (a)
$t=0.1$, (b) $t=2$, and for the proposed model at (c) $t=0.1$, (d) $t=2$. The
white solid line represents the interface; and the color plot represents the
scalar concentration field.

For the first case of a translating bubble, the fluid velocity is $\vec{u}=1$
everywhere in the domain. This case is more challenging compared to the
stationary bubble case in Section 7.2.1, because of the non-zero droplet-
Peclet number ($Pe=ud/D=4$) of this flow. The scalar concentration fields
obtained by solving the proposed model and from the previous approach are
shown in Figure 15 at two different times, $t=0.1$, and $t=2$, on a grid of
size $128\times 128$. The results from the previous approach not only show the
leakage of scalar into the bubble, but also exhibit the characteristic
dispersion errors associated with a non-dissipative scheme, at early times
($t=0.1$). At the final time of $t=2$, the scalar almost completely leaks into
the bubble. The dispersion errors generate negative values of the scalar
concentration and results in realizability issues of the scalar concentration
field. The negative values in the results from the previous approach is
evident in Figure 16, where $c$ and $\phi$ are plotted along the wall-normal
direction at $x=0$. Unlike the results from the previous approach, the
proposed model neither show signs of leakage of the scalar into the bubble,
nor the dispersion errors.

FIGURE 16: The scalar concentration field $c$ and the volume fraction field
$\phi$ along the $x=0$ line for the case of scalar diffusion around a moving
bubble in a channel, simulated on a grid of size $128\times 128$, at (a)
$t=0.1$, and (b) $t=2$.

The scalar leakage error $c_{e}$ defined in Eq. (37) and computed along the
wall-normal direction at $x=0$ is plotted in Figure 17. The non-zero values of
$c_{e}$ for the previous approach imply a significant leakage of the scalar
into the bubble region. The integrated leakage error defined in Eq. (38) and
evaluated along the wall-normal direction at $x=0$ is listed in Table 2.
Similar to the stationary-bubble case, the integrated leakage error with the
proposed method is significantly lower compared to the previous approach.

FIGURE 17: The scalar leakage error $c_{e}$ along the $x=0$ line for the moving-bubble case, simulated on a grid of size $128\times 128$, at (a) $t=0.1$, and (b) $t=2$. $I_{error}$ | $t=0.1$ | $t=2$
---|---|---
previous approach | $3.54\times 10^{-2}$ | $3.79\times 10^{-2}$
proposed method | $6.27\times 10^{-7}$ | $6.28\times 10^{-7}$
TABLE 2: Integrated leakage error for the moving-bubble case.

In the second case, for a deforming bubble, the initial fluid velocity is
taken to be $\vec{u}=1$ everywhere in the domain. The flow is sustained with a
unit acceleration, $g=1$ along the streamwise direction. The density of the
gas and the liquid are $\rho_{g}=1$ and $\rho_{l}=10$, respectively. The
viscosity of the gas and the liquid are $\mu_{g}=1.81\times 10^{-5}$ and
$\mu_{l}=8.9\times 10^{-4}$, respectively. This corresponds to a liquid
Reynolds number of about $Re_{l}=\rho_{l}ud/\mu_{l}=450$. The surface tension
is varied from $\sigma=0$ to $\sigma=5$, to study the effect of bubble breakup
on the evolution of the scalar concentration field. This corresponds to the
Eötvös numbers of $Eo=\infty$ to $Eo=3.2\times 10^{-3}$, respectively, where
the Eötvös number is defined as $Eo=\Delta\rho gd^{2}/\sigma$.

FIGURE 18: The evolution of the scalar concentration field around a deforming
bubble in a channel for a surface tension value of $\sigma=0$. The white solid
line represents the interface; and the color plot represents the scalar
concentration field.

Figure 18 shows the evolution of the bubble shape and the scalar concentration
field at various times up to $t=10$ for the surface tension value $\sigma=0$
($Eo=\infty$), simulated on a grid of size $128\times 128$. The bubble deforms
due to the relative difference in the body force between the bubble and the
surrounding liquid, which is effectively a buoyancy force, and develops skirt
on the two ends as can be seen at $t=1.5$ and $t=2$. Clift et al. (2005)
characterized the shapes of the bubble on a $Re-Eo$ plot, and the $Re=450$ and
$Eo=\infty$ corresponds to a skirted spherical-cap regime, which matches well
with our observation, although the current setup is a two dimensional one.
Later, the bubble breaks up leaving behind a trail of small bubbles, which
eventually fills the channel as can be seen at $t=10$. With an increase in the
surface tension values (decrease in $Eo$) the deformation of the bubble is
reduced. Figure 19 shows the time evolution of bubble shape and the scalar
concentration for the surface tension value of $\sigma=5\times 10^{-4}$
($Eo=32$).

FIGURE 19: The evolution of the scalar concentration field around a deforming
bubble in a channel for a surface tension value of $\sigma=5\times 10^{-4}$.
The white solid line represents the interface; and the color plot represents
the scalar concentration field.

One consequence of the bubble deformation and breakup is that the flux of the
scalar across the channel in the wall-normal direction is reduced because of
the reduced cross-sectional area. This can be quantified by calculating the
flux of the scalar at the wall defined as $f_{w}=D\lvert\vec{\nabla}c\rvert$.
The average flux of the scalar, $\langle f_{w}\rangle$, is plotted in Figure
20 versus time for various values of surface tension, where the average is
computed over top and bottom walls. For the case of $\sigma=0$, it was seen in
Figure 18 that the bubble would eventually span the whole channel blocking the
flow of scalar. As a result, the flux of the scalar can be seen to drop to
zero for time $t\gtrapprox 8$. With an increase in the surface tension values,
the deformation of the bubble decreases, and therefore the flux of the scalar
increases. At the surface tension value of $\sigma=5$, the bubble essentially
remains circular throughout simulation; and the average scalar flux reaches a
limiting value of $\langle f_{w}\rangle=0.075445$ that was obtained for the
stationary bubble case in Section 7.2.1.

FIGURE 20: The average scalar flux versus time for various values of surface
tension. The thin black dashed line represents the average flux for the
stationary bubble case in Section 7.2.1.

###  7.3 Charged ions in a drop

In this section, the simulation of reorganization of ions with unbalanced
charge within a drop will be presented. This case illustrates the
applicability of the proposed scalar-transport model in Eq. 12 for modeling
electrokinetics in two-phase flows. To model this phenomenon, the scalar-
transport equation is recast into the Nernst-Planck equation

$\frac{\partial c^{\pm}}{\partial
t}+\vec{\nabla}\cdot(\vec{u}c^{\pm}\pm\phi\lambda\vec{E}c^{\pm})=\vec{\nabla}\cdot\left[D^{\pm}\left\\{\vec{\nabla}c^{\pm}-\frac{(1-\phi)\vec{n}c^{\pm}}{\epsilon}\right\\}\right],$
(39)

where $c^{+}$ and $c^{-}$ represent the cationic and anionic concentration
fields, respectively; $D^{+}$ and $D^{-}$ represent the cationic and anionic
diffusion coefficients, respectively; $\lambda$ is the electrical mobility and
$\vec{E}$ is the electric field. Here, the relative velocity $\vec{u}_{r}$ is
replaced by the effective electromigration velocity experienced by the ions,
and is given by $\vec{u}_{r}=\lambda\vec{E}$. To close the system of
equations, the Nernst-Planck equation is combined with the Gauss’s law

$\vec{\nabla}\cdot(\varepsilon\vec{E})=\rho^{f},$ (40)

where $\rho^{f}=z^{+}ec^{+}-z^{-}ec^{-}$ is the free charge density; $z^{+}$
and $z^{-}$ represent the cationic and anionic valences, respectively; $e$ is
the elementary charge; and $\varepsilon$ is the electrical permittivity of the
electrolyte. Making an electrostatic approximation, i.e., assuming that the
time variation of the magnetic field is much slower compared to that of the
electric field, the electric field can be shown to be irrotational using
Faraday’s law (Griffiths, 2005). Therefore, the electric field can be written
as

$\vec{E}=-\vec{\nabla}\psi,$ (41)

where $\psi$ is the electrostatic potential. Now, invoking the Einstein-
Smoluchowski relation

$D^{\pm}=\frac{\lambda k_{B}T}{z^{\pm}e},$ (42)

where $k_{B}$ is the Boltzmann constant, and by re-scaling the electrostatic
potential $\psi$, with the thermal voltage $V_{T}$, as

$\tilde{\psi}=\frac{\psi}{V_{T}},$ (43)

where $V_{T}=k_{B}T/(z^{\pm}e)$, and $T$ is the absolute temperature, the
Nernst-Planck equation in Eq. 39 can be rewritten in the form

$\frac{\partial c^{\pm}}{\partial
t}+\vec{\nabla}\cdot\big{\\{}\vec{u}c^{\pm}\mp\phi
D^{\pm}(\vec{\nabla}\tilde{\psi})c^{\pm}\big{\\}}=\vec{\nabla}\cdot\left[D^{\pm}\left\\{\vec{\nabla}c^{\pm}-\frac{(1-\phi)\vec{n}c^{\pm}}{\epsilon}\right\\}\right],$
(44)

which will be used in the current setup. Electrostatic phenomena and
hydrodynamics can be coupled through the Maxwell stress tensor either as a
body force or as a divergence of the tensor (Saville, 1997, Berry et al.,
2013), and will be part of a separate study and is beyond the scope of this
work.

In this study, the setup consists of a stationary circular liquid drop of
radius $R=0.2$, in a quiescent surrounding medium of another fluid. The
density of both fluids are chosen to be $\rho_{d}=\rho_{s}=1$, where the
subscripts $d$ and $s$ represents the drop and surrounding fluid properties.
The viscosity of both fluids are chosen to be $\mu_{d}=\mu_{s}=10^{-3}$. The
ion concentrations are initially uniform within the drop and are normalized by
the reference value $c_{ref}=10^{9}$ as $\tilde{c}^{\pm}=c^{\pm}/c_{ref}$. The
normalized initial ion concentrations are chosen to be $\tilde{c}^{+}=1$, and
$\tilde{c}^{-}=0$. The diffusivity of the ions are chosen to be
$D^{\pm}_{d}=0.01$ and $D^{\pm}_{s}=0$.

FIGURE 21: (a) The normalized ion concentration field, and (b) the re-scaled
electrostatic potential field for the case of charged ions in a drop at time
$t=1$, simulated on a grid of size $512\times 512$.

Figure 21 shows the ion concentration field $\tilde{c}^{+}$ and the resulting
electrostatic potential field $\tilde{\psi}$ at time $t=1$, simulated on a
grid of size $512\times 512$. The ions in the drop experiences the Columbic
repulsive force because of the unbalanced net positive charge. This force
drives the ions apart, towards the interface. This is counteracted by the
diffusion in the opposite direction. The ions, however, cannot cross the
interface because of the non-conducting surrounding medium; and therefore,
they arrive at a state where the electromigration is balanced by the
diffusion. Figure 22 shows the ion concentration $\tilde{c}^{+}$ and the
resulting electrostatic potential $\tilde{\psi}$ at time $t=1$ along the line
$x=0$ in the domain. The results of the simulation from four different grid
sizes were chosen to show the convergence of the results.

FIGURE 22: (a) The normalized ion concentration $\tilde{c}^{+}$, and (b) the
re-scaled electrostatic potential $\tilde{\psi}$ along the $x=0$ line for the
case of charged ions in a drop at time $t=1$, simulated on grids of different
sizes.

###  7.4 Droplet-laden turbulent channel flow

Up until here in Section 7, simple tests cases in 1D and 2D were presented.
Here, a direct numerical simulation of a droplet-laden turbulent channel flow
with passive scalar quantity that is confined to the surrounding fluid
(carrier fluid) is presented. This test case illustrates: the validity of the
positivity criterion, the applicability and robustness of the proposed scalar-
transport model, and the prevention of unphysical numerical leakage by the
model for complex two-phase turbulent flow regimes.

The setup consists of a three-dimensional channel of size $8H\times 2H\times
4H$ in $x,y,z$ as shown in Figure 23, where $H=1$ is the channel half height.
The streamwise direction (along the $x$ coordinate) and the spanwise direction
(along the $z$ coordinate) have periodic boundary conditions, and the wall-
normal direction (along the $y$ coordinate) has no-slip walls.

FIGURE 23: A schematic of the domain for the case of droplet-laden channel
flow.

A unity density ratio $\rho*=\rho_{d}/\rho_{s}=1$ and a unity viscosity ratio
$\mu^{*}=\mu_{d}/\mu_{s}=1$ are considered in this case such that the
kinematic viscosities $\nu=\mu_{d}/\rho_{d}=\mu_{s}/\rho_{s}=0.01$ are same
for the droplet fluid and the surrounding fluid. The computation is performed
for a Reynolds number of $Re=UH/\nu=3000$, defined based on the mean
centerline velocity $U$. This corresponds to a shear Reynolds number of
$Re_{\tau}=u_{\tau}H/\nu=150$, which is similar to the droplet-laden turbulent
flow study by Scarbolo et al. (2015).

The domain is discretized into a uniform grid of size $256\times 64\times
128$, which results in a grid spacing of size $\Delta y^{w}=4.6875$ in wall
units, where the superscript $w$ represents a non-dimensional quantity scaled
by the wall variables; e.g. $\Delta y^{w}=(\Delta y)u_{\tau}/\nu$, where $\nu$
is the kinematic viscosity, $u_{\tau}=\sqrt{\tau_{w}/\rho}$ is the shear
velocity, and $\tau_{w}$ is the wall shear stress. No subgrid-scale model is
used in the computation, since the grid resolution is considered to be
sufficiently fine to resolve the essential turbulent scales (Kim et al.,
1987).

A single-phase channel flow with surrounding fluid properties is simulated
until it reaches a statistically stationary state. Then, droplets of diameter
$d=0.4$ are added to the channel at a distance of $y^{+}=54$ from the walls,
and are placed in the flow such that their centers are equally spaced in the
$x-z$ plane. A total of $100$ droplets are added and the resulting volume
fraction is $0.054$. A passive scalar quantity of diffusivities $D_{s}=0.5$
and $D_{d}=0$, in the surrounding fluid and the droplet fluid, respectively,
is added to the surrounding fluid with a uniform initial concentration of $1$
at the same time when the droplets are introduced to the flow. The scalar is
confined to the surrounding fluid region since the ratio of diffusivities in
two fluids is $D_{d}/D_{s}=0$. The corresponding Schmidt numbers are
$Sc_{s}=0.02$ and $Sc_{d}=\infty$, in the surrounding fluid and the droplet
fluid, respectively. Since the $Sc_{s}$ is less than $1$, we expect that the
scalar microscale, $\eta_{s}$, to be larger than the Kolmogorov scale, $\eta$,
in the flow and it scales as $\eta_{s}\sim(D_{s}/\nu)^{3/4}\eta$ (Davidson,
2015). Therefore, the grid requirement for the surrounding flow would be
sufficient to resolve the smallest scales of the scalar. However, if $Sc_{s}$
was larger than $1$, then a more refined grid for the scalar needs to be used
to resolve the Batchelor scale (see, Schwertfirm and Manhart, 2007).

FIGURE 24: The snapshots of the droplet-laden turbulent channel flow at
various times. The left column shows the three-dimensional view of the
droplets and the flow structures in the channel. The blue colored surfaces
represent the droplet interface. The flow structures are shown on the bottom
half of the channel by plotting the isosurfaces of the second-invariant of the
velocity gradient tensor, $Q=500$, colored by the local vorticity magnitude in
the flow. The velocity magnitude is also plotted on a $x-z$ plane at a
distance of $y=0.2H$ from the bottom wall. The right column shows the scalar
concentration field on a $x-z$ cross-section midplane in the channel. The
solid white lines represent the droplet interface.

The surface tension between the droplet fluid and the surrounding fluid is
chosen such that the droplet-shear Weber number is $We_{\tau,d}=\rho
u_{\tau}^{2}d/\sigma=1$. Since the droplets are of size $d^{w}=96$ in wall
units, they are larger than the Kolmogorov scale of $\eta^{w}\approx 2$ wall
units; and they undergo breakup and coalescence. The snapshots of the flow are
shown in Figure 24 at various times. It appears from the snapshots that the
droplet size is increasing, and therefore the number of droplets is decreasing
with time, which is consistent with the observation of Scarbolo et al. (2015).
For more details on the dynamics of the breakup and coalescence of the
droplets, refer to the study by Scarbolo et al. (2015). Figure 24 also shows
the scalar concentration field along the $x-z$ cross-section center plane of
the channel at various times. Since the droplets were seeded at a distance of
$y^{+}=54$ from the walls, the cross-section plane do not cut any of the
droplets at the early time of $t^{w}=13.5$, therefore no droplets can be seen
at this time. However, at later times, the droplets undergo coalescence and
breakup and therefore can be seen in this view. Evidently, there was no
leakage of the scalar quantity into the droplets or violation of the
positivity of the scalar throughout the simulation. This illustrates the
robustness of the proposed method and the validity of the positivity
criterion.

In this case, the maximum cell-Peclet number of the simulation is around
$Pe_{c}=1.875$; and the positivity of the scalar is maintained throughout the
simulation thought it violates the Positivity criterion by a small amount.
This is because the Positivity criterion is a sufficient condition, and
therefore cell-Peclet number values higher than $1$ might also result in
maintaining positivity for the scalar quantity, as already illustrated in
Section 7.1.2 for a 1D case. However, repeating the droplet-laden channel flow
simulation with a coarser grid (half the size of the original grid size in
every direction) with a maximum cell-Peclet number of $Pe_{c}=3.75$, resulted
in violation of the positivity of the scalar. Therefore, the grid size needs
to selected such that the cell-Peclet number is maintained close to or less
than $1$.

##  8 Conclusions

In this work, we proposed a novel transport model for the simulation of
scalars in two-phase flows. The scalars are usually confined to one of the
phases in a two-phase flow due to its disparate values of diffusivity and
mobility in the two phases; and this typically poses a challenge for any
numerical method in resolving the gradient of the scalar at the material
interface.

We therefore developed and verified a general scalar-transport model for two-
phase flows, particularly for interfaces modeled using a phase-field (diffuse-
interface) method. We showed that our newly proposed model equation prevents
the artificial numerical diffusion (unphysical leakage) of the scalar from one
phase to the other, while maintaining the positive values for the scalar
concentration field throughout the simulation, albeit the use of central-
difference schemes for the discretization of the operators. The use of
central-difference scheme is to achieve a non-dissipative implementation that
is crucial for the simulation of turbulent flows.

We proved that the model maintains the positivity of the scalar concentration
values—a crucial realizability requirement for the simulation of
scalars—provided the grid resolution was fine enough to satisfy the given
positivity criterion. The grid resolution required to satisfy the positivity
criterion is also based on the fact that the grid size should be small enough
to resolve the smallest physical scales present in the flow.

The prevention of unphysical numerical leakage of the scalar across the
interface was achieved by enforcing consistency between the transport of the
scalar concentration field and the phase field. It was shown that the
equilibrium solution for the proposed scalar-transport model also exhibits a
hyperbolic tangent function similar to the equilibrium solution for the phase-
field equation, thereby making the transport of the scalar and the phase field
consistent.

At the end, the proposed model was assessed of its accuracy, robustness,
effectiveness in maintaining the positivity, and in preventing unphysical
leakage across the interface, by simulating a wide range of two-phase flows,
starting from simple one-dimensional droplet advection flows to complex three-
dimensional droplet-laden turbulent flows. The droplet advection flows were
used to verify the positivity criterion by choosing parameters that both
violate and satisfy the positivity criterion, and by showing that the
positivity is maintained in all the cases that satisfied the criterion. The
proposed model was also recast into a Nernst-Planck equation and was used to
simulate electrokinetics of two-phase flows. Finally, the droplet-laden
turbulent channel flow showed the robustness of the model in simulating
scalars in complex high-Reynolds-number turbulent two-phase flows.

## Acknowledgments

This investigation was supported by the Office of Naval Research, Grant #
N00014-19-1-2425. S. S. Jain is also funded by the Franklin P. and Caroline M.
Johnson Fellowship. The computations in this paper were performed on the
Shepard, Armstrong and Yellowstone clusters at the Stanford HPC Center,
supported through awards from the National Science Foundation, DOD HPCMP, and
Office of Naval Research. A preliminary report on this work has been published
in Jain and Mani (2019) as a Center-for-Turbulence-Research Annual Research
Brief and the authors acknowledge Dr. Ronald Chan’s helpful comments on this
report.

## References

  * Villegas et al. (2016) L. R. Villegas, R. Alis, M. Lepilliez, S. Tanguy,  A ghost fluid/level set method for boiling flows and liquid evaporation: application to the leidenfrost effect,  Journal of Computational Physics 316 (2016) 789–813.
  * Lal (2008) R. Lal,  Carbon sequestration,  Philosophical Transactions of the Royal Society B: Biological Sciences 363 (2008) 815–830.
  * Fernandez et al. (2014) D. Fernandez, P. Maurer, M. Martine, J. Coey, M. E. Mobius,  Bubble formation at a gas-evolving microelectrode,  Langmuir 30 (2014) 13065–13074.
  * Takagi and Matsumoto (2011) S. Takagi, Y. Matsumoto,  Surfactant effects on bubble motion and bubbly flows,  Annual Review of Fluid Mechanics 43 (2011) 615–636.
  * Alke et al. (2009) A. Alke, D. Bothe, M. Kroeger, H. Warnecke,  Vof-based simulation of conjugate mass transfer from freely moving fluid particles,  Computational Methods in Multiphase Flow V (2009) 157–168.
  * Bothe and Fleckenstein (2013) D. Bothe, S. Fleckenstein,  A volume-of-fluid-based method for mass transfer processes at fluid particles,  Chemical Engineering Science 101 (2013) 283–302.
  * Hassanvand and Hashemabadi (2012) A. Hassanvand, S. H. Hashemabadi,  Direct numerical simulation of mass transfer from taylor bubble flow through a circular capillary,  International journal of heat and mass transfer 55 (2012) 5959–5971.
  * Berry et al. (2013) J. Berry, M. Davidson, D. J. Harvie,  A multiphase electrokinetic flow model for electrolytes with liquid/liquid interfaces,  Journal of Computational Physics 251 (2013) 209–222.
  * Jain et al. (2020) S. S. Jain, A. Mani, P. Moin,  A conservative diffuse-interface method for compressible two-phase flows,  Journal of Computational Physics 418 (2020) 109606.
  * Mirjalili et al. (2019) S. Mirjalili, C. B. Ivey, A. Mani,  Comparison between the diffuse interface and volume of fluid methods for simulating two-phase flows,  International Journal of Multiphase Flow 116 (2019) 221–238.
  * Mirjalili et al. (2017) S. Mirjalili, S. S. Jain, M. Dodd,  Interface-capturing methods for two-phase flows: An overview and recent developments,  Center for Turbulence Research Annual Research Briefs (2017) 117–135.
  * Pericleous et al. (1995) K. Pericleous, K. Chan, M. Cross,  Free surface flow and heat transfer in cavities: the sea algorithm,  Numerical Heat Transfer 27 (1995) 487–507.
  * Mehdi-Nejad et al. (2004) V. Mehdi-Nejad, J. Mostaghimi, S. Chandra,  Modelling heat transfer in two-fluid interfacial flows,  International journal for numerical methods in engineering 61 (2004) 1028–1048.
  * Davidson and Rudman (2002) M. R. Davidson, M. Rudman,  Volume-of-fluid calculation of heat or mass transfer across deforming interfaces in two-fluid flow,  Numerical Heat Transfer: Part B: Fundamentals 41 (2002) 291–308.
  * Ma and Bothe (2013) C. Ma, D. Bothe,  Numerical modeling of thermocapillary two-phase flows with evaporation using a two-scalar approach for heat transfer,  Journal of Computational Physics 233 (2013) 552–573.
  * Mittal and Moin (1997) R. Mittal, P. Moin,  Suitability of upwind-biased finite difference schemes for large-eddy simulation of turbulent flows,  AIAA journal 35 (1997) 1415–1417.
  * Cahn and Hilliard (1958) J. W. Cahn, J. E. Hilliard,  Free energy of a nonuniform system. i. interfacial free energy,  The Journal of chemical physics 28 (1958) 258–267.
  * Allen and Cahn (1979) S. M. Allen, J. W. Cahn,  A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening,  Acta metallurgica 27 (1979) 1085–1095.
  * Anderson et al. (1998) D. M. Anderson, G. B. McFadden, A. A. Wheeler,  Diffuse-interface methods in fluid mechanics,  Annual review of fluid mechanics 30 (1998) 139–165.
  * Lowengrub and Truskinovsky (1998) J. Lowengrub, L. Truskinovsky,  Quasi–incompressible cahn–hilliard fluids and topological transitions,  Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 454 (1998) 2617–2654.
  * Chen et al. (1998) L. Q. Chen, J. Shen, et al.,  Applications of semi-implicit fourier-spectral method to phase field equations,  Computer Physics Communications 108 (1998) 147–158.
  * Jacqmin (1999) D. Jacqmin,  Calculation of two-phase navier–stokes flows using phase-field modeling,  Journal of Computational Physics 155 (1999) 96–127.
  * Liu and Shen (2003) C. Liu, J. Shen,  A phase field model for the mixture of two incompressible fluids and its approximation by a fourier-spectral method,  Physica D: Nonlinear Phenomena 179 (2003) 211–228.
  * Badalassi et al. (2003) V. E. Badalassi, H. D. Ceniceros, S. Banerjee,  Computation of multiphase systems with phase field models,  Journal of computational physics 190 (2003) 371–397.
  * Yue et al. (2004) P. Yue, J. J. Feng, C. Liu, J. Shen,  A diffuse-interface method for simulating two-phase flows of complex fluids,  Journal of Fluid Mechanics 515 (2004) 293.
  * Yang et al. (2006) X. Yang, J. J. Feng, C. Liu, J. Shen,  Numerical simulations of jet pinching-off and drop formation using an energetic variational phase-field method,  Journal of Computational Physics 218 (2006) 417–428.
  * Kim (2012) J. Kim,  Phase-field models for multi-component fluid flows,  Communications in Computational Physics 12 (2012) 613–661.
  * Chiu and Lin (2011) P.-H. Chiu, Y.-T. Lin,  A conservative phase field method for solving incompressible two-phase flows,  Journal of Computational Physics 230 (2011) 185–204.
  * Olsson and Kreiss (2005) E. Olsson, G. Kreiss,  A conservative level set method for two phase flow,  Journal of computational physics 210 (2005) 225–246.
  * Mirjalili and Mani (2020) S. Mirjalili, A. Mani,  Consistent, energy-conserving momentum transport for simulations of two-phase flows using the phase field equations,  Journal of Computational Physics (2020) 109918.
  * Huang et al. (2020a) Z. Huang, G. Lin, A. M. Ardekani,  Consistent, essentially conservative and balanced-force phase-field method to model incompressible two-phase flows,  Journal of Computational Physics 406 (2020a) 109192\.
  * Huang et al. (2020b) Z. Huang, G. Lin, A. M. Ardekani,  Consistent and conservative scheme for incompressible two-phase flows using the conservative allen-cahn model,  Journal of Computational Physics 420 (2020b) 109718\.
  * Kataoka (1986) I. Kataoka,  Local instant formulation of two-phase flow,  International Journal of Multiphase Flow 12 (1986) 745–758.
  * Kim and Moin (1985) J. Kim, P. Moin,  Application of a fractional-step method to incompressible navier-stokes equations,  Journal of computational physics 59 (1985) 308–323.
  * Mirjalili et al. (2020) S. Mirjalili, C. B. Ivey, A. Mani,  A conservative diffuse interface method for two-phase flows with provable boundedness properties,  Journal of Computational Physics 401 (2020) 109006.
  * Laney (1998) C. B. Laney, Computational gasdynamics, Cambridge university press, 1998.
  * Moin and Verzicco (2016) P. Moin, R. Verzicco,  On the suitability of second-order accurate discretizations for turbulent flow simulations,  European Journal of Mechanics-B/Fluids 55 (2016) 242–245.
  * Patankar (1980) S. Patankar, Numerical heat transfer and fluid flow, CRC press, 1980.
  * Versteeg and Malalasekera (2007) H. K. Versteeg, W. Malalasekera, An introduction to computational fluid dynamics: the finite volume method, Pearson Education, 2007.
  * Jain et al. (2018) S. S. Jain, A. Mani, P. Moin,  A conservative diffuse-interface method for the simulation of compressible two-phase flows with turbulence and acoustics,  Center for Turbulence Research Annual Research Briefs (2018) 47–64.
  * Falgout and Yang (2002) R. D. Falgout, U. M. Yang,  hypre: A library of high performance preconditioners,  in: International Conference on Computational Science, Springer, pp. 632–641.
  * Clift et al. (2005) R. Clift, J. R. Grace, M. E. Weber, Bubbles, drops, and particles, Courier Corporation, 2005.
  * Griffiths (2005) D. J. Griffiths, Introduction to electrodynamics, 2005.
  * Saville (1997) D. Saville,  Electrohydrodynamics: the taylor-melcher leaky dielectric model,  Annual review of fluid mechanics 29 (1997) 27–64.
  * Scarbolo et al. (2015) L. Scarbolo, F. Bianco, A. Soldati,  Coalescence and breakup of large droplets in turbulent channel flow,  Physics of Fluids 27 (2015) 073302.
  * Kim et al. (1987) J. Kim, P. Moin, R. Moser,  Turbulence statistics in fully developed channel flow at low reynolds number,  Journal of fluid mechanics 177 (1987) 133–166.
  * Davidson (2015) P. A. Davidson, Turbulence: an introduction for scientists and engineers, Oxford university press, 2015\.
  * Schwertfirm and Manhart (2007) F. Schwertfirm, M. Manhart,  DNS of passive scalar transport in turbulent channel flow at high Schmidt numbers,  International Journal of Heat and Fluid Flow 28 (2007) 1204–1214.
  * Jain and Mani (2019) S. S. Jain, A. Mani,  Scalar transport modeling for two-phase flows with a diffuse-interface method,  Center for Turbulence Research Annual Research Briefs (2019) 137–152.