\section{Introduction}
\label{Sec1}

The main difference
between single and
multi-phase flow is
in the topological
changes of the
boundaries separating
coexisting phases
and/or flow
domain(s).
The wide
spectrum of
physical
phenomena
at these
boundaries
governs
the multi-phase flow
and have
to be
accounted for
in all
multi-phase
flow
regimes.
The numerical solution
of the flow problems
with moving boundary
conditions is source
of complexity of the
numerical algorithms
and has led to the
development of the
abundance of
the computational
techniques
\cite{anderson1998,osher03,prosperetti06,trygg11}.

In fluid mechanics,
the most popular
approach
splits
the control
volume  containing
neighboring phases
with the sharp
interface
smooth enough
to compute
mean Gaussian
curvature
needed to determine
capillary terms
\cite{osher03,trygg11}.
While the division
of the given control
volume
is not arbitrary
(as one would expect
in the continuum mechanical model),
the sharp interface
model guarantees
the highest
possible numerical
resolution
of the interface
(at worst in the single
control volume and its
nearest neighbors).

In the phase field methods,
the control volume
is divided
by the diffusive interface
\cite{anderson1998}.
The main advantage
of the phase field
methods is in their
physical interpretation
based on theory of
the Ginzburg-Landau
free energy functional \cite{kim12}.
However,
the order
parameter
in the phase field
model does not
have clear
physical interpretation.
This is caused by
the fact
Allen-Cahn \cite{allen1979}
and Cahn-Hilliard \cite{cahn1958}
equations governing
the phase field
does not
have a mathematical
form of
the classical transport
equation.
To overcome
these difficulties
and preserve
conservation
of the order parameter
the Lagrange multiplayer
techniques has been
adopted
\cite{kim12,kim14}.

In the fluid
mechanics
community
it is assumed,
that the numerical
sharp and
diffusive
interface
models
follow the thermodynamic
dividing plane model of Gibbs \cite{gibbs1874}
and the smooth transition region
model of van der Waals
\cite{waals1979}.
We  note,
that if this
indeed is the case,
they
inherit
the known
shortcomings
of these models, too.
For instance,
as the position of the
Gibbs dividing plane $\mathbf{x}_\Gamma\,[m]$ is
arbitrary \cite{lang15,faust2018}
the (signed) distance
measured from
this position
is not defined
as well.
In the phase
field models
it is often
assumed the radius of
the curvature $R\,[m]$
has to be much
larger than the interface
thickness $\epsilon_h\,[m]$.
Otherwise a
numerical robustness
 and
physical interpretation
of the  phase
field model(s)
order
parameter are
impaired
\cite{kim12,lang15}.
As phase transitions
(e.g. boiling)
typically
start at the molecular
level from
the mixture of
fluid and its vapor
(no bubble exists at
this stage, yet)
the assumption
about $\epsilon_h/R\!\ll\!1$
ratio
seems to
restrict
the range
of the phase
field model(s)
application.

Whilst
the mathematical and
following physical
imperfections
of the aforementioned
model equations
can be compensated
by a careful choice
of the  sophisticated numerical
techniques
(e.g. higher-order advection schemes \cite{liu1994},
      the height function method \cite{cummins05},
	  adaptive mesh refinement \cite{popinet09},
      Lagrange multiplayer \cite{kim14},
      etc.)
leading, eventually, to plethora
of successful numerical predictions;
the most striking flaw of
the deterministic sharp/diffusive
interface models  seems
to be their detachment from
the experimental reality
\cite{vrij1973,aarts2004}.

Moreover,
one could also ask:
how likely is it ,
that solution of
the  five \cite{cahn1958,allen1979,hirt81,osher1988,olsson05},
different (sets of)
apparently unrelated
partial differential
equations,
describing
the same physical phenomenon(a),
using different
boundary conditions,
leads to the same
(or very similar) predictions?

The answers
to some of these  questions
are given in \cite{twacl17}.
Unlike described by the sharp
and diffusive interface models
the mesoscopic interface $\Gamma$
between two phases
is rough and/or turbulent.
It is disturbed
by thermal fluctuations
in the domain where
particles
of coexisting
phases violently
exchange energy
dependent on
the temperature
and cohesive forces of
the given system.
For this reason,
the width of
this domain is
estimated to be
$\epsilon_h \!\sim\! \sqrt{k_B T/\sigma} \! > \! 0$
where
$k_B\,[J/K]$ is
the Boltzman constant,
$T\,[K]$ is
the absolute
temperature and
$\sigma\,[J/m^2]$ is
the surface tension
coefficient
\cite{vrij1973,aarts2004}.
The discrepancy between
true nature of the interface
and its smooth, deterministic models
has
led the present author to
development of the description
that is alternative for the
sharp and diffusive interface
models and yet it shows
they can be considered as
complementary
components
of the same
statistical description
\cite{twacl15,twacl17,twacl21}.

In the present work,
recently
derived \cite{twacl15}
equation describing
evolution of
the intermittency region
is coupled with
the Navier-Stokes
incompressible
flow
solver
and used to
simulate
flow of
the rising,
two-dimensional
gas bubble.
In addition,
we use
non-equilibrium
solution
of the modified
Ginzburg-Landau
functional \cite{twacl17,twacl21}
derived by
the present author
to introduce
the new mechanism
of the
topological changes
in two-phase flow.
To this goal,
the two cases are considered:
first, when characteristic
length scale field $\epsilon_h \lr \bx,t \rr$
 governing
the local intermittency region
width is constant,
and second,
when $\epsilon_h \lr \bx,t \rr$
is given by
the prescribed
analytical
function.

This
paper is organized
as follows.
In the next section,
derivation of the
intermittency region
evolution equation
and its
relation with
the modified
Ginzburg-Landau
free energy
functional
is recalled
\cite{twacl17,twacl21}.
Next, we
present the
numerical
method and
set-up of
the two-dimensional
numerical
experiment.
Finally,
obtained results
are discussed
in context of
difference between
equilibrium and
non-equilibrium
effects in
two-phase flows
and differences
between the predictions
of the deterministic
and the statistical
intermittency
region models.

\section{Statistical model of the intermittency region}
\label{Sec2}

In this section we briefly recall
arguments already presented in
works \cite{twacl15,twacl17,twacl21}.
Therein the derivation,
physical interpretation
and convergence
of numerical solution
of
the  intermittency region
evolution equation
is investigated.
The new contribution in
these work is description
of the intermitency region
in terms of probability of finding
the mesoscopic sharp interface $\Gamma$.
Moreover,
establishing its relation to
the minimization procedure
of the modified Ginzburg-Landau
free energy functional.

\subsection{Derivation of the cumulative distribution function transport equation}
\label{ssec21}

We
assume,
the interface
between two-phases is
domain named
the microscopic intermittency region.
Therein, the
position of the sharp interface $\Gamma$
(herein defined on the molecular level)
can be found with
the non-zero probability
described by the cumulative
distribution function $\alpha\,[-]$.
It is important
to note,
the intermittency
region paradigm
was first introduced and
developed for modeling of
turbulence interface interactions
\cite{brocchini01a,mwaclawczyk11,waclawczyk2015}.
In \cite{twacl17}
we have drew  analogy
between the interface
interacting with turbulence
and sharp interface $\Gamma$
agitated by
thermal fluctuations.
In both cases,
these phenomena
can be described
as stochastic
processes.
This allowed us
to use the conditional
averaging procedure \cite{pope88}
and eddy viscosity
model \cite{waclawczyk2015}
to close unknown
correlations
as in
turbulence-interface
interaction model.
As a result,
equation
for evolution
of the probability
of finding the
mesoscopic sharp interface
$\Gamma$: $\alpha\,[-]$,
 is derived
\begin{equation}
\pd{\alpha}{t} \!+\! \mathbf{w} \nabla \alpha
\!=\! \nabla \cdot \left( D \nabla \alpha \right)
-  |\avG{\mathbf{n}_\Gamma}| \nabla{\left( D \Sigma \right)}  \!\cdot\! \mathbf{n}_\gamma
\label{eq1}
\end{equation}
where $\mathbf{w}\,[m/s]$
is  velocity of
the regularized
interface $\gamma$
defined by the expected
position of $\Gamma$: $\alpha \left( \psi\!=\!0 \right) \!=\!1/2$ ,
$\mathbf{n}_\gamma$ is vector
normal to $\gamma$,
$\mathbf{n}_\Gamma$ is vector
normal to $\Gamma$,
$\avG{\cdot}$ is
the conditional average operator
and
$D\!=\!C \left( \mathbf{x},t \right) \epsilon_h \lr \bx,t \rr\,[m^2/s]$
is the diffusivity
coefficient.
$C \left( \mathbf{x},t \right)$,
$\epsilon_h \lr \bx,t \rr$
are velocity and
length scales
characterizing
the intermittency
region, respectively.
The first RHS term in \Eq{eq1}
is responsible for spreading
of the intermittency
region width,
the second
one has been shown to be
the counter gradient
diffusion
responsible
for its contraction
\cite{mwaclawczyk11,waclawczyk2015}.

As \Eq{eq1} contains
an unclosed RHS term with
the unknown
surface density
$\Sigma\,[m^2/m^3]$,
we have used results
in \cite{cahn1958,olsson05,chiu11}
for its conservative closure
and its physical interpretation \cite{twacl21}.
Hence, $\alpha\,[-]$ evolution equation reads
\begin{equation}
\pd{\alpha}{t} \!+\! \nabla \! \cdot \! \left( \mathbf{w} \alpha \right) \!=\! \nabla \cdot \left[ D |\nabla \alpha|\mathbf{n}_\gamma \!-\! C \alpha \left( 1-\alpha \right) \mathbf{n}_\gamma \right].
\label{eq2}
\end{equation}
The coefficients $C\,[m/s],\, D\,[m^2/s]$
in \Eq{eq2} specify
the characteristic
length $\epsilon_h \!\sim\! D/C\,[m]$
and time $\tau_h \! \sim \! D/C^2\,[s]$
scales
governing its
solution
and in general
can be
functions of
space and time.
Thus,
in works \cite{kmtw14,twaclawczyketal14,waclawczyk2015}
it was shown \Eq{eq2} can
be used to predict
the intermittency region
evolution due to interaction
of turbulent eddies with
the macroscopic
interface $\gamma$
(the same interface as
in the standard volume of fluid (VOF)
and level-set (SLS)
numerical models).

The steady
state
solution
of \Eq{eq2}
with $\epsilon_h\!=\! const.$
and
$\mathbf{w}\!=\!\mathbf{u} \!=\!0$
is given by
the regularized
Heaviside function
\begin{equation}
\alpha \lr  \psi \rr = \frac{1}{1+\exp{\left( -\psi  \left( \bx,t \right)/\epsilon_h \right)}}=
\frac{1}{2} \left[ 1+\tanh{\left( \frac{\psi  \left( \bx,t \right)}{2\epsilon_h} \right)} \right]
\label{eq3}
\end{equation}
and its inverse
function that is
the signed distance
from the expected position of
the regularized interface $\gamma$
defined by the level-set $\psi \left( \alpha \!=\!1/2 \right) \!=\! 0$
\begin{equation}
\psi \left( \alpha \right) = \epsilon_h \ln{\left[ \frac{\alpha\left( \psi \right)}
	{1 - \alpha \left(  \psi \right)} \right]}.
\label{eq4}
\end{equation}
As noticed
by the
present author
\cite{twacl15},
\Eqs{eq3}{eq4}
are known
to characterize
the cumulative
distribution $\alpha \lr  \psi \rr$,
and quantile $\psi  \left( \alpha \right)$
functions
of the logistic
distribution.
Additionally,
the gradient
of $\alpha \lr  \psi \rr$
given by the formula
\begin{equation}
\nabla \alpha = \frac{\tilde{\delta} \lr  \alpha \rr}{\epsilon_h} \nabla \psi,
\label{eq5}
\end{equation}
where $\tilde{\delta} \lr  \alpha \rr/\epsilon_h \!=\! \alpha \left( 1\!-\!\alpha \right)/\epsilon_h$ is
the probability density function
of the logistic distribution.
Computation of $\nabla \alpha\,[1/m]$
with \Eq{eq5}
bridges
the numerical problems
with approximation
of steep functions
gradients (see $\alpha$ when $\epsilon_h \rightarrow 0$ and analysis in  \cite{twacl15}).
\Eq{eq5} can also be interpreted as
calculation of $\nabla \alpha$
in the $\mathbf{n}_\gamma$ direction
\begin{equation}
\nabla \alpha \cdot \mathbf{n}_\gamma \!=\! \pd{\alpha}{\psi} \!=\! \frac{\tilde{\delta} \lr  \alpha \rr}{\epsilon_h}.
\label{eq6}
\end{equation}
This observation was
recently used in \cite{holger2021}
to reformulate the
modified Ginzburg-Landau
free energy functional
for its solution
in the local
coordinate system
normal in every
point to
the interface
$\gamma$.
While \Eq{eq6} is exact,
the Laplacian
$\nabla \!\cdot\! \left\{ \nabla \left[ \alpha \lr  \psi \rr \!\cdot\! \mathbf{n}_\gamma \right] \right\}$
in the normal direction
can be only approximated
in similar way
when the interface
curvature $\kappa \!\ll\! 1$.

Curiously,
the definition
of the chemical potential
and the corresponding
modified Ginzburg-Landau
functional
in \cite{twacl17} and \cite{holger2021}
is the same (to details
in the coefficients and the fact
that different solution norms
are used \cite{mirja20}).
The advantage
of description
used herein
over the classical phase
field models
is  sound physical
interpretation
of the order
parameter $\alpha\,[-]$, its gradient $\nabla \alpha\,[1/m]$
and the order parameter inverse function $\psi\,[m]$,
see Eqs.~(\ref{eq3},\ref{eq5},\ref{eq4})
respectively.

This makes the difference
as for instance,
one can use the physical
meaning of the cumulative
distribution function given
by \Eq{eq3}
to pinpoint
the position
of the Gibbs dividing
surface $\mathbf{x}_\Gamma$.
The Gibbs
dividing surface
can be
defined
as the expected
position of
the sharp interface $\Gamma$
agitated by
the stochastic
thermal fluctuations.
This position
is given by the level-sets
$\alpha \left( \psi\!=\!0 \right)\!=\!1/2$
and describes
the two dimensional
smooth surface that
is called
the sharp interface
in the deterministic
models of
the intermittency
region.

In \cite{twacl15}
it was noticed that
substitution of \Eq{eq5}
into \Eq{eq2} results in
\begin{equation}
\pd{\alpha}{t} + \mathbf{w} \nabla \alpha =  \nabla \cdot \left[ C  \tilde{\delta} \lr  \alpha \rr \left( |\nabla\psi|-1 \right) \mathbf{n}_\gamma \right],
\label{eq7}
\end{equation}
where
$\mathbf{n}_\gamma \!=\! \nabla \alpha /|\nabla \alpha|\!=\! \nabla \psi /|\nabla \psi|$.
In the present work
we separate the advection
and re-initialization
steps in \Eq{eq7},
which leads to
\begin{equation}
\pd{\alpha}{t} \!+\! \mathbf{w} \nabla \alpha \!=\!
\pd{\alpha}{t} \!+\! \frac{\tilde{\delta} \lr  \alpha \rr}{\epsilon_h} \mathbf{w} \!\cdot\! \nabla \psi  \!=\! 0,
\label{eq8}
\end{equation}
\begin{equation}
\pd{\alpha}{\tau} \!=\!  \nabla\! \cdot \! \left[ C  \tilde{\delta} \lr  \alpha \rr \left( \left |\nabla\psi \right| \!-\! 1 \right) \mathbf{n}_\gamma \right],
\label{eq9}
\end{equation}
where $\tau\,[s]$
is time needed to
obtain the equilibrium
solution of \Eq{eq9}.
In our works,
this form
of \Eq{eq2}
is preferred
for numerical
solution
and theoretical
analysis.

At first,
the choice of
the partial
differential algebraic
equation (\ref{eq4}-\ref{eq7})
for the model
of two-phase
flow in fluid
mechanics appears
to be surprising.
However,
it is known \cite{ascher1998},
the equations of this type
describe systems with
the large separation of
scales (e.g. chemical kinetics in quasi steady state
        and partial equilibrium approximations,
        molecular dynamics, optimal control etc.).
This suits our goal
of derivation of the
intermittency region model
that is closer to
experimental reality
than known  in the literature
thermodynamic
descriptions of
the interface
between two-phases.

\subsection{Re-initialization as minimization of the modified Ginzburg-Landau functional}
\label{ssec22}

In  \cite{twacl21} it was shown,
when the
characteristic
length $\epsilon_h\,[m]$
and velocity
$C\,[m/s]$
scales
governing evolution
of the intermittency
region are constant,
this region is remaining
in the equilibrium
state according
to the modified
Ginzburg-Landau
functional describing
interfacial energy
of the intermittency
region \cite{twacl17}.
In this case,
the solution
of \Eq{eq9} in the Conservative Level-Set method (CLS) \cite{olsson05}
is called
re-initialization
of the function $\alpha$.
In this
numerical method,
the re-initialization
step is used to restore
deformed due to numerical
errors in the advection step
function $\alpha$
to its original
form given by \Eq{eq3}.

In the limit $\epsilon_h \! \rightarrow \! 0$,
\Eqs{eq5}{eq9}
can be used to  derive
the re-initialization equation
of the signed-distance function $\psi$
used in the Standard Level-Set (SLS) method
\cite{twacl15,twacl17}.
Here,
the goal of
the re-initialization step
is similar.
The goal function
in reconstruction
procedure is $\psi\,[m]$ disturbed
by the numerical
errors introduced
during
the advection
step
\cite{osher03}.
\Eqs{eq4}{eq7} explain
that these two functions
(and two models of
the intermittency
region)
are closely related.
Moreover,
the present
author
has proven \cite{twacl17},
the re-initialization
procedure
in the SLS \cite{osher03}
and CLS \cite{olsson05} level-set
methods
are equivalent
to minimization
of the
free
energy
functional
containing
the term
which accounts
for
the regularized
interface $\gamma$
deformation.

We note
that before
work \cite{twacl17}
the presence
of term
with $\nabla \alpha$
and the interface
curvature $\kappa$
in
the definition
of chemical potential
was postulated \cite{folch1999}.
Extension of
the original
free energy
Ginzburg-Landau functional
by the additional term containing
the interface curvature
was also proposed
in \cite{jamet2008}.
However,
as noticed
in  \cite{mirja20}
the derivations
in \cite{folch1999,jamet2008}
and \cite{twacl15,twacl17}
were carried out
in the different
norms.
Moreover,
as described
in \Sec{Sec1}
the present
work is
driven by the
attempt to paint
more
physical
and general
picture
of the interface.
In particular
when it is
agitated by
the thermal
or turbulent
fluctuations,
whereas
the  works \cite{folch1999,jamet2008}
use mainly
mathematical
arguments.

\subsection{Non-equilibrium of energy ratio in intermittency region}
\label{ssec23}

In most of
the analytical and numerical
discussions
presented in the literature
it is assumed,
the interface
or intermittency
region are passive
actors advected
by the fluid
velocity.
In our view,
the role
of intermittency
region in the multi-phase
flow is much
more complex.
The intermittency
region is open system
that can interplay
with neighboring
phases through exchange
of the mass and/or energy.
In particular,
in the case of phase changes
it mediates in the exchange
of mass and energy.
For these
reasons,
the intermittency
region between
two weakly miscible
phases
may not
be in
the equilibrium
state as it
is typically
assumed in
the known
sharp and
diffusive
interface
models.
Next
it is demonstrated,
that local
non-equilibrium
effects can
influence
distribution
of the material
properties
of the phases,
and hence,
affect
the two-phase
flow scenario.

Driven
by this
general picture
in the recent
work \cite{twacl21}
the non-equilibrium
solution of
\Eqs{eq7}{eq8}
was introduced.
Therein we have
shown, that statistical
model of the intermittency region
can also be
extended to the case
when $\epsilon_h \lr \bx,t \rr\!\ne\!const.$.
In particular
we have derived
the free energy
functional with
$\nabla \alpha \!\cdot\! \nabla \epsilon_h$
term  showing
the variation of
volume and
interfacial
forces work ratio
should affect the
cumulative distribution
function $\alpha \lr  \psi \rr$
shape.
Moreover,
we have demonstrated
the stationary solution
of \Eq{eq9}
can be used
as approximate
solution of the
non-equilibrium
problem
when used
in a spirit of
re-initialization
as in the CLS and
SLS numerical
methods.

In what follows we repeat
derivation of the approximate
solution already presented
in \cite{twacl21}.
For mathematical details
regarding the modified Ginzburg-Landau
free energy functional
the interested reader should
refer to \cite{twacl21}.

The equilibrium
condition
obtained
from
the stationary
solution
to \Eq{eq9}
reads
\begin{equation}
\nabla \alpha \!=\! |\nabla \alpha| \mathbf{n}_\gamma \!=\! \frac{\alpha \left( 1\!-\!\alpha \right)}{\epsilon_h \lr \bx,t \rr} \mathbf{n}_\gamma.
\label{eq10}
\end{equation}
\Eq{eq10} is
formulated
in the direction $\mathbf{n}_\gamma$
normal to the  regularized
interface $\gamma$,
hence,
it may
be rewritten as
\begin{equation}
\pd{\alpha}{\psi} \left| \pd{\psi}{\mathbf{x}} \right| \!=\! \frac{1}{\epsilon_h \lr \bx,t \rr}  \alpha \left( 1 \!-\! \alpha \right),
\label{eq11}
\end{equation}
where it
is assumed
$\pdl{\alpha}{\psi} \! > \! 0$
meaning
$\alpha \lr  \psi \rr$
is expected
to be
the cumulative
distribution
function
with infinite
support
analogously
to \Eq{eq5}.
Next, we assume
$|\nabla \psi| \!\equiv\! 1$
in \Eqs{eq10}{eq11}.
As a result,
substitution
of \Eq{eq10}
into \Eq{eq2} with
$D \left( \mathbf{x},t \right)\!=\!C \epsilon_h \lr \bx,t \rr$
let us derive
\Eq{eq7} with
the variable
characteristic
length scale
governing width
of the intermittency
region.
The assumption
$|\nabla \psi| \!\equiv\! 1$
means
the signed distance
function $ \psi  \left( \bx,t \right) $
spans
the space where
surface
averaged oscillations
of the sharp interface
$\Gamma$
take place.
On average,
these
oscillations occur
only in  the direction
$\mathbf{n}_\gamma$ normal
to the expected
position $\psi\!=\!0$
of the regularized
interface $\gamma$.

Further,
it is noticed
at each point
$\left( \mathbf{x},t \right)$
of the field
$\epsilon_h \lr \bx,t \rr$
the signed
distance
function
$\psi  \left( \bx,t \right)$ is given.
Hence,
the knowledge of
the field $\psi  \left( \bx,t \right)$
gives
$\left( \mathbf{x},t \right)$
and thus
$\epsilon_h \lr \bx,t \rr$.
Therefore,
we introduce
$\epsilon_h^{\psi} \left( \mathbf{x},t \right)$
denoting
$\epsilon_h \lr \bx,t \rr$
determined
using $\psi  \left( \bx,t \right)$.
This
let us
to integrate
\Eq{eq11}
in the local
coordinate
system attached
to the
regularized
interface $\gamma$.
As $\gamma$
is defined
by $\psi  \left( \bx,t \right)\!=\!0$,
$\psi  \left( \bx,t \right)$ is
the normal coordinate
with the origin at
$\psi  \left( \bx,t \right)\!=\!0$
of this local
system.
At each
fixed
point
of given
$\alpha \lr  \psi \rr,\,\psi  \left( \alpha \right),\,\epsilon_h \lr \bx,t \rr$
fields
this
integration
reads
\begin{equation}
\int^{1/2}_{\alpha \lr  \psi \rr} \frac{d\alpha'}{\alpha' \left( 1\!-\!\alpha' \right)} = \int^{0}_{\psi  \left( \alpha \right)} \frac{d \psi'}{\epsilon_h^{\psi'} \! \lr \bx,t \rr}.
\label{eq12}
\end{equation}
The integration
(\ref{eq12})
is performed
from the
arbitrary point
located at the signed-distance
from the regularized interface
$\alpha \lr  \psi \rr\!-\!\psi  \left( \alpha \right)$
to the expected
position of the
regularized
interface $\psi \left( \alpha \!=\!1/2 \right) \!=\! 0$.
One notes
the LHS integration
in \Eq{eq12}
does not
assume
or result in
any specific form/shape
of the function $\alpha \lr  \psi \rr$.

To reconstruct
the equilibrium
solution when
$\epsilon_h \lr \bx,t \rr \!=\! const.$
it is necessary
to preserve the mapping
between $\alpha \lr  \psi \rr\!-\!\psi  \left( \alpha \right)$,
see \Eq{eq4}.
For this reason,
it is more convenient
to reformulate
the RHS integral
in \Eq{eq12} using
a variable
substitution
as follows
\begin{equation}
\int^{0}_{\psi  \left( \alpha \right)} \frac{d \psi'}{\epsilon_h^{\psi'} \! \lr \bx,t \rr}
\!=\! \psi  \left( \alpha \right) \int^0_1 \frac{dt'}{\epsilon_h^{t' \psi} \! \lr \bx,t \rr} \!=\! \psi  \left( \alpha \right) \mathcal{I \lr \psi \rr}
\label{eq13}
\end{equation}
where $t' \! \in \!  [ 0,1 ]$
is the  parameter
such that $\psi' \!=\! t' \psi$
and $d \psi' \!=\! dt' \psi$,
furthermore
$\mathcal{I \lr \psi \rr}$ is used to denote
the integral on the RHS
of \Eq{eq13}.

After integration
of \Eq{eq11} with \Eq{eq13}
one obtains
\begin{equation}
\psi  \left( \alpha \right)  = \frac{1}{\mathcal{I \lr \psi \rr}} ln \left[ \frac{\alpha \lr  \psi \rr}{1-\alpha \lr  \psi \rr} \right].
\label{eq14}
\end{equation}
It is noted, \Eq{eq14}
can be the
source of
numerical
problems when $\alpha \rightarrow 0$
or $\alpha \rightarrow 1$.
In these limits,
due to the finite computer
arithmetic
computation of the
natural logarithm is
ill conditioned
(see discussion of this
problem in \cite{twacl15})
and
some measures must
be taken to account
for it.
In the code,
it is chosen to solve
\Eq{eq9} only in the
region where
the numerical
solution is still
possible assuming:
$|\nabla \psi|\!=\!1$
and hence $\pdl{\alpha}{\tau}\!=\!0$,
outside.

At
the given,
arbitrary
point $\left( \mathbf{x},t \right)$,
the
signed
distance $\psi  \left( \bx,t \right)$
has the known value.
For this reason,
at the point $(\mathbf{x},t)$
the integral
$\mathcal{I \lr \psi \rr} \!=\! const.$
and thus
the inverse
relation
is also true
\begin{equation}
\alpha \lr  \psi \rr \!=\! \frac{1}{1+\exp{ \left[ - \psi  \left( \alpha \right) \mathcal{I \lr \psi \rr} \right] }}.
\label{eq15}
\end{equation}
In the present work,
\Eq{eq15} is called
the approximate,
semi-analytical solution
of \Eq{eq9}.

When
the field
$\epsilon_h \lr \bx,t \rr \!=\! const.$,
\Eq{eq12} and \Eq{eq13}
reduce to
the equilibrium
solution,
which is guaranteed
by the
definition of $\mathcal{I \lr \psi \rr}$.
Thus,
the mapping given
by \Eq{eq14} or
the form of $\alpha \lr  \psi \rr$
given by \Eq{eq15}
can be employed during
numerical solution
of the system given
by \Eqs{eq8}{eq9}
to model
how
the
$\epsilon_h \lr \bx,t \rr$
field
is affecting
changes of
the
cumulative
distribution
function
$0 \!<\! \alpha \lr  \psi \rr \!<\! 1$
profile.

In what follows,
we will couple
described above
approximate solution
of \Eq{eq9} given by \Eq{eq15}
with the Navier-Stokes
equation solver
to predict how
non-equilibrium
affects affect
the flow
of
 the two-dimensional
rising gas bubble.
\begin{multicols}{1}
	\begin{figure*}[htbp]
		\begin{minipage}{.5\textwidth}
			\centering  \includegraphics[width=0.6\textwidth]{./pics/case.png}

			\caption{The schematic sketch of the flow case, all
				dimensions are in meters. The colored spot shows
				fixed position
				of variation in
				the characteristic length scale field
				$\epsilon_h \lr \bx,t \rr$ given by \Eq{eq19}.}
			\label{fig1}
		\end{minipage}

	\end{figure*}
\end{multicols}
\section{Numerical method}
\label{Sec3}

This section describes
the numerical method used
in the present paper
for coupling of Eqs.~(\ref{eq8}, \ref{eq9}, \ref{eq15})
with the incompressible
Navier-Stokes solver.
In addition,
details
of incorporation
of the non-equilibrium
effects into the
two-phase
flow solver
are described.

\subsection{Flow solver and one-fluid model implementation}
\label{ssec31}

In the present paper,
the second-order accurate finite
volume flow solver
in the collocated variable
arrangement is used \cite{schaefer06,peric02}.
For detailed
description of the algorithm(s)
and presentation of the results,
see
previous
works by the present author \cite{waclawczyk05,waclawczyk06,waclawczyk07,waclawczyk08,waclawczyk08_2,waclawczyk08_3,waclawczyk08_phd,waclawczyk2015}.
As in these
references, herein,
the incompressible
Navier-Stokes
equation
in the conservative
form is
solved.
The momentum
conservation
equation is
coupled with
the Poisson
equation to bind
 pressure
and velocity fields
in the SIMPLE
procedure.
The equilibrium between
the pressure gradient and mass forces
is established using
the same numerical
operators for their
discretization \cite{zun07}.
In this standard solver
of incompressible flow,
the two-phase
flow model
using the one-fluid
formulation \cite{trygg11}
is implemented;
in the present
work the presence
of the capillary
forces
is neglected.
Material properties of
one-fluid change according to
relations
\begin{eqnarray}
 \rho = \rho_1 \alpha + \rho_2 \left( 1-\alpha \right),\\
 \label{eq16}
 \mu  = \mu_1 \alpha + \mu_2 \left( 1-\alpha \right),
 \label{eq17}
\end{eqnarray}
where $\rho\,[kg/m^3],\,\mu\,[Pa \!\cdot\! s]$ denote
density and dynamic viscosity
of the mixture of
two phases; $\rho_k,\,\mu_k,\,\,k=1,2$
are   density and dynamic viscosity
of pure fluid and gas
phases, respectively.
\begin{figure*}[htbp]
	\begin{minipage}{.49\textwidth}
		\centering  \includegraphics[width=0.9\textwidth]{./pics/mass_conservation/cmass3}
	\end{minipage}%
	\begin{minipage}{.49\textwidth}
		\centering  \includegraphics[width=0.9\textwidth]{./pics/mass_conservation/cmass4}
	\end{minipage}
	\caption{The mass conservation of gas during
		coupled simulations with constant $\epsilon_h \lr \bx,t \rr\!=\!\epsilon_{h,b}$
		and variable  characteristic length scale field $\epsilon_h \lr \bx,t \rr$.
		Total mass is computed using the c.d.f-s $\alpha \lr  \psi \rr$  (left)
		and $\alpha_d \! \lr  \psi \rr$(right).
		$\epsilon_{h,max}$ denotes the maximum magnitude of $\epsilon_h \lr \bx,t \rr$ given by \Eq{eq12}.
		}
	\label{fig5}
\end{figure*}

The advection \Eq{eq8}
is discretized in time using
the second-order accurate implicit
Euler method.
For discretization of
the convective term in space,
the deferred-correction
method with
the second-order
accurate TVD MUSCL
ﬂux limiter \cite{xue98}
is employed.
The same
temporal and spatial
discretization is used
for velocity components
in the Navier-Stokes
equation.
Re-initialization equation
(\ref{eq9}) is integrated
in time
using the explicit third-order
accurate TVD Runge-Kutta
method \cite{gottlieb98}.
The constrained interpolation
\cite{twacl17,twacl21}
is used to discretize
$\tilde{\delta} \lr  \alpha \rr\!=\!\alpha \left( 1-\alpha \right)$,
in \Eq{eq9}.
For more details
about numerical schemes
and discretization
used to solve \Eqs{eq8}{eq9},
refer to
Appendixes in \cite{twacl15,twacl17,twacl21}.

It is noticed,
\Eq{eq7}
in the limit $\epsilon_h \rightarrow 0$
guarantees the statistical model
of the intermittency region
will converge to the solution obtained
using the sharp interface model(s).
We choose
base or minimal
width
of the intermittency
region as $\epsilon_{h,b}\!=\!\sqrt{2}\Delta x/4$;
in  \cite{twacl15}
it was checked this
selection guarantees the second
order convergence of \Eq{eq9}
and its higher-order derivatives.

Furthermore it is noted,
in the  statistical model
of the intermittency region,
unlike in
some phase
field solvers
\cite{kim12},
there
is no danger
the material
properties
of the fluids
will have
nonphysical
values.
This is due
to the fact
that $\alpha \lr  \psi \rr$
is the cumulative
distribution
function
(see \Eq{eq3})
and it
is incorporated
into the solution
procedure.
Therefore,
from its definition
it is always
bounded $0\!<\!\alpha \lr  \psi \rr\!<\!1$.
\begin{figure*}[t]

	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/alpha_bepsh10.0049.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/alpha_bepsh10.0099.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/alpha_bepsh10.0149.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/alpha_bepsh10.0199.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/alpha_bepsh10.0249.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/alpha_bepsh10.0299.png}
	\end{minipage}%
	\begin{minipage}{.15\textwidth}
		\includegraphics[width=0.2\textwidth,height=1.25\textwidth,angle=0]{./pics/alpha_bar.png}
	\end{minipage}
	%
	%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/dalpha_bepsh10.0049.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/dalpha_bepsh10.0099.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/dalpha_bepsh10.0149.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/dalpha_bepsh10.0199.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/dalpha_bepsh10.0249.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/dalpha_bepsh10.0299.png}
	\end{minipage}%
	\begin{minipage}{.15\textwidth}
		\includegraphics[width=0.2\textwidth,height=1.25\textwidth,angle=0]{./pics/alpha_bar.png}
	\end{minipage}

	\caption{Level-sets $\left[ 0.01,0.5,0.99 \right]$ and color maps
		of the cumulative distribution functions: $0\!<\!\alpha \lr  \psi \rr\!<\!1$
		(white contours, top row) and
		$0\!<\!\alpha_d \! \lr  \psi \rr\!<\!1$
		(black contours, bottom row).
		%
		Position of both functions is
		affected by the characteristic length
		scale field $\epsilon_h \lr \bx,t \rr\,[m]$ with $\epsilon_{h,max}\!=\!11\epsilon_{h,b}\,[m]$,
		see \Fig{fig3} (bottom row).
		%
		The snapshots from the left
		to the right are taken every $\Delta t \!=\! 0.05\,[s]$.}
	\begin{picture}(0,0)
	\put(4,271){{{\color{white}{\small{t=0.05\,[s]}}}}}
	\put(84,271){{{\color{white}{\small{t=0.10\,[s]}}}}}
	\put(164,271){{{\color{white}{\small{t=0.15\,[s]}}}}}
	\put(244,271){{{\color{white}{\small{t=0.20\,[s]}}}}}
	\put(324,271){{{\color{white}{\small{t=0.25\,[s]}}}}}
	\put(404,271){{{\color{white}{\small{t=0.30\,[s]}}}}}
	\end{picture}
	\label{fig2}
\end{figure*}

\begin{figure*}[t]

	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh1.0_pic/vofandDvof_bepsh1.0049.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh1.0_pic/vofandDvof_bepsh1.0099.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh1.0_pic/vofandDvof_bepsh1.0149.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh1.0_pic/vofandDvof_bepsh1.0199.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh1.0_pic/vofandDvof_bepsh1.0249.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh1.0_pic/vofandDvof_bepsh1.0299.png}
	\end{minipage}%
	\begin{minipage}{.15\textwidth}
		\includegraphics[width=0.2\textwidth,height=1.25\textwidth,angle=0]{./pics/ndiff_bar1.png}
	\end{minipage}
	%
	%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/vofandDvof_bepsh10.0049.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/vofandDvof_bepsh10.0099.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/vofandDvof_bepsh10.0149.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/vofandDvof_bepsh10.0199.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/vofandDvof_bepsh10.0249.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/vofandDvof_bepsh10.0299.png}
	\end{minipage}%
	\begin{minipage}{.15\textwidth}
		\includegraphics[width=0.2\textwidth,height=1.25\textwidth,angle=0]{./pics/ndiff_bar2.png}
	\end{minipage}

	\caption{Level-sets $\left[ 0.01,0.5,0.99 \right]$ of the probability $\alpha \lr  \psi \rr$ with $\epsilon_{h,b}\!=\!const.$
		(white contours) and function $\alpha_d \! \lr  \psi \rr$
		approximating the exact solution
		of \Eq{eq9} (black contour).
		$\alpha_d \! \lr  \psi \rr$ is
		affected by the characteristic length
		scale field $\epsilon_h \lr \bx,t \rr\,[m]$ (colour map).
		$\epsilon_{h,max}\!=\!\epsilon_{h,b}\!=\!const.$ (top row),
		$\epsilon_{h,max}\!=\!11\epsilon_{h,b}$ (bottom row).
		The snapshots from left
		to right are taken every $\Delta t \!=\! 0.05\,[s]$.}
	%
	\begin{picture}(0,0)
	\put(4,273){{{\color{white}{\small{t=0.05\,[s]}}}}}
	\put(84,273){{{\color{white}{\small{t=0.10\,[s]}}}}}
	\put(164,273){{{\color{white}{\small{t=0.15\,[s]}}}}}
	\put(244,273){{{\color{white}{\small{t=0.20\,[s]}}}}}
	\put(324,273){{{\color{white}{\small{t=0.25\,[s]}}}}}
	\put(404,273){{{\color{white}{\small{t=0.30\,[s]}}}}}
	\end{picture}
	\label{fig3}
\end{figure*}
\subsection{Non-equilibrium effects in the two-phase flow}
\label{ssec32}

In the case
when $\epsilon_h \lr \bx,t \rr\!=\!const.$,
and the regularized
interface velocity $\mathbf{w}$
in \Eq{eq8} is prescribed
by  known analytical formula,
the numerical solution
of the transport equations \Eqs{eq8}{eq9}
with constraint given by \Eq{eq4}
are described in
every detail
in \cite{twacl15}.
In the present work,
the non-equilibrium
solution of \Eq{eq9},
approximated
by \Eq{eq15}
is implemented in
described
in the previous
section
two-phase
flow solver.

To account for
non-equilibrium
due to variations
in the characteristic
length scale field $\epsilon_h \lr \bx,t \rr$,
first, \Eqs{eq8}{eq9} are
solved
with $\epsilon_{h,b}\!=\!const.$,
see results presented
in the top rows of
\Figs{fig2}{fig4}
and in the bottom
row of \Fig{fig3}
(white contours).
This solution provides
$\psi  \left( \alpha \right)$ and $\alpha \lr  \psi \rr$
fields at every
time step.
Further,
$\psi  \left( \alpha \right)$ is
considered
to define
the local
coordinate systems
(where normal in
 every point of
 the interface $\gamma:$
 $\psi \left( \alpha\!=\!1/2 \right) \!=\! 0$
 is given by $\mathbf{n}_\gamma$).
These local
systems are
used to compute
the integral
\Eq{eq18} in \Eq{eq15}.
The field created
by these local
coordinate systems
is visualized
in \Fig{fig4}
as the level-sets
of $\psi  \left( \alpha \right)$
function.

In \cite{twacl21},
 physical
interpretation of $\psi  \left( \alpha \right)$
field was
proposed.
It is interpreted
as the space where
the averaged oscillations
of the mesoscopic, sharp
interface $\Gamma$
take place.
The
oscillations
of $\Gamma$ (in the present work) defined
at the molecular level,
are governed by
thermal fluctuations
and for this reason
can not be
modeled
directly in
the continuum
mechanical
description.
\Eq{eq15}
mimics their
variation in
the averaged
sense,
see \cite{twacl21}
for full discussion.

Let us consider
the last
outer iteration of
the SIMPLE algorithm,
after the
single time
step $\Delta t\,[s]$
were
the momentum conservation
and Poisson
equations are
iterated to
enforce stronger
coupling between
the pressure
and velocity
fields.
At the last outer iteration,
$N_\tau\!=\!4$ re-initialization
steps $\Delta \tau \!=\! D/C^2\,[s]$
is carried out
advancing \Eq{eq9}
in time $\tau$.
As discussed
in \Sec{ssec22},
this is equivalent
to the minimization
of the modified
Ginzburg-Landau
free energy
functional.

Before the
solution of \Eq{eq9}
in time $\tau$,
the integral
\begin{equation}
\begin{split}
&\mathcal{I \lr \psi \rr} \!=\! \int^0_1 \frac{dt'}{\epsilon_h^{t' \psi} \! \lr \bx,t \rr}  \! \approx \! \\
 &\qquad\qquad \frac{1}{6} \left[ \frac{1}{ \epsilon_h \texttt{(inp)} }
+ \frac{4}{ \epsilon_h \texttt{(inm)} }
+ \frac{1}{ \epsilon_h \texttt{(int)}} \right].
\label{eq18}
\end{split}
\end{equation}
is computed
in each control
volume \texttt{(inp)} where the
signed distance
function $\psi  \left( \alpha \right)$
is reconstructed.
Other
Indexes in \Eq{eq18}
denote, respectively,
control volumes where:
the interface $\alpha \left( \psi \!=\!0 \right) \!=\! 1/2$ is located $\texttt{(int)}$
and in the distance in between $\texttt{(inm)}$.

It was checked
in \cite{twacl21},
approximation
of the  integral (\ref{eq18})
using the third-order
accurate
Simpson rule
is sufficient
to obtain
$\mathcal{I \lr \psi \rr}$ field.
Therein, more detailed
description of the integration
scheme used to evaluate $\mathcal{I \lr \psi \rr}$
field is given, too.
Next, obtained
in each control
volume values of $\mathcal{I \lr \psi \rr}$
are used to
compute $\alpha_d \! \lr  \psi \rr$
field given
by \Eq{eq15}
and than material
properties are
determined setting
$\alpha \lr  \psi \rr\!=\!\alpha_d \! \lr  \psi \rr$
in
\Eqs{eq16}{eq17}.
Afterwards,
the mass
and diffusive fluxes
are calculated
at the faces of all
control volumes.
Finally,
the segregated
Navier-Stokes equation
solver proceeds
to the solution
on the next time
step $t\!+\!\Delta t$
solving euqations
for velocity components
and repeating
described
procedure.

\section{Numerical experiment}
\label{Sec4}

In this section,
simple two-phase
flow experiment
is carried out
to demonstrate
the first application
of the numerical method \cite{twacl21}
used for the approximate solution
of Eqs.~(\ref{eq4}, \ref{eq8}, \ref{eq9})
with \Eq{eq15} which
can account for
the non-equilibrium
effects due to
$\epsilon_h \lr \bx,t \rr \! \ne  \! const.$
In particular,
during discussion
of the numerical results obtained
in this study,
the differences
 between
statistical and deterministic
description of the two-phase
flow separated by the
intermittency region
in the equilibrium
or/and non-equilibrium
states
are emphasized.

\subsection{Simulation set-up}
\label{sec41}

In what follows,
the two-phase
flow of
the rising
due to positive buoyancy
gas bubble is
studied.
For this
demonstration,
density and viscosity
contrasts
are selected to be
$\rho_1/\rho_2\!=\!6$ and
$\mu_1/\mu_2\!=\!20$,
respectively.
The gravitational
acceleration is chosen
to have its standard magnitude
$\mathbf{g}\!=[0,9.8]\,[m/s^2]$.
The two-phase
flow problem
is solved
in the domain
$\Omega = <0,0.1>\!\times\!<0,0.1>\,[m^2]$,
 discretized
using
uniform grid with
$2^7\times2^7$
control volumes.
Time
step size
is set to
$\Delta t \!=\! 2.5\cdot10^{-4}\,[s]$
what guarantees
the maximum Courant
number $Cu_{max} \!<\! 0.25\,[-]$
during the whole
simulation lasting
$T\!=\![0,0.5]\,[s]$.

To couple the pressure
and velocity fields,
$N_{out} \!=\! 4$
outer iterations
per time step
$\Delta t$ are used.
The same time step
$\Delta t\,[s]$ is used
to advance
\Eq{eq8} in time $t\,[s]$.
Re-initialization equation (\ref{eq8})
is solved in time $\tau\,[s]$
using the time
step size
$\Delta \tau \!=\! 0.5 \!\cdot\! D/C^2 \!=\! 0.5 \!\cdot\! \epsilon_{h,b}\,[s]$
and  $N_\tau\!=\!4$ iterations per $\Delta \tau$.
Model function $\alpha_d \! \lr  \psi \rr$ given by \Eq{eq15}
is computed after each
re-initialization cycle.
\begin{figure*}[t]

	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh1.0_pic/ndel_bepsh1.0049.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh1.0_pic/ndel_bepsh1.0099.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh1.0_pic/ndel_bepsh1.0149.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh1.0_pic/ndel_bepsh1.0199.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh1.0_pic/ndel_bepsh1.0249.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh1.0_pic/ndel_bepsh1.0299.png}
	\end{minipage}%
    \begin{minipage}{.15\textwidth}
    	\includegraphics[width=0.2\textwidth,height=1.25\textwidth,angle=0]{./pics/ndel_bar.png}
    \end{minipage}
	%
	%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/ndel_bepsh10.0049.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/ndel_bepsh10.0099.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/ndel_bepsh10.0149.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/ndel_bepsh10.0199.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/ndel_bepsh10.0249.png}
	\end{minipage}%
	\begin{minipage}{.16\textwidth}
		\centering\includegraphics[width=0.95\textwidth,height=1.5\textwidth,angle=0]{./pics/gauss_bepsh10.0_pic/ndel_bepsh10.0299.png}
	\end{minipage}%
	\begin{minipage}{.15\textwidth}
		\includegraphics[width=0.2\textwidth,height=1.25\textwidth,angle=0]{./pics/ndel_bar.png}
	\end{minipage}

	\caption{Level-sets $\left\{ - 3 \!\cdot\! n \!\cdot\! \epsilon_{h,b}, 0, 3 \!\cdot\! n \! \cdot \! \epsilon_{h,b}  \right\}$ where $n\!=\!1,\ldots,8$
		     of the signed distance function field $\psi  \left( \alpha \right)$.
		     The color map
			 depicts joint probability $\alpha_d \! \lr  \psi \rr \left( 1\!-\!\alpha_d \! \lr  \psi \rr \right)$
		     of finding the phase $\alpha_d \! \lr  \psi \rr$ used to define the macroscopic interface $\gamma$ position.
			 $\epsilon_{h,max}\!=\!\epsilon_{h,b}\!=\!const.$ (top row),
			 $\epsilon_{h,max}\!=\!11\epsilon_{h,b}$ (bottom row).
			 The snapshots from the left
			 to the right are taken every $\Delta t \!=\! 0.05\,[s]$.
		 	}
   	\begin{picture}(0,0)
   	\put(4,273){{{\color{white}{\small{t=0.05\,[s]}}}}}
   	\put(84,273){{{\color{white}{\small{t=0.10\,[s]}}}}}
   	\put(164,273){{{\color{white}{\small{t=0.15\,[s]}}}}}
   	\put(244,273){{{\color{white}{\small{t=0.20\,[s]}}}}}
   	\put(324,273){{{\color{white}{\small{t=0.25\,[s]}}}}}
   	\put(404,273){{{\color{white}{\small{t=0.30\,[s]}}}}}
   	\end{picture}
	\label{fig4}
\end{figure*}

In \Fig{fig1},
sketch of the
computational domain
with the initial position
of the gas bubble $[0.05,0.02]\,[m]$
having diameter
$D\!=\!0.02\,[m]$
is presented.
The same figure depicts
constant in time position $\mathbf{x}_r\!=\![0.06,0.04]\,[m]$
of the maximal value
of the characteristic length scale
field $\epsilon_h \lr \bx,t \rr$ (colored drop),
see also the color map in \Fig{fig3}.

The spatial variation
of the characteristic
length scale field
$\epsilon_h \lr \bx,t \rr$ is predefined
to be bell shaped
with the analytical formula
\begin{equation}
\epsilon_h \lr \bx,t \rr \!=\! \epsilon_{h,b} \left\{ 1 + A \!\cdot\! \exp{\left[ - \left( \frac{r}{20\epsilon_{h,b}} \right)^2 \right]} \right\}
\label{eq19}
\end{equation}
where $r\!=\!\left[ \left( x-x_r \right)^2+\left( y-y_r \right)^2 \right]^{1/2}$,
$\mathbf{x}_r$ denotes center
of $\epsilon_h \lr \bx,t \rr$,
and the base width of the intermittency
region is  $\epsilon_{h,b}=\sqrt{2}\Delta x/4$.

To analyze influence
of $\epsilon_h \lr \bx,t \rr$ magnitude
on the obtained results,
in particular
on the mass
conservation computed
using  $\alpha \lr  \psi \rr$ or $\alpha_d \! \lr  \psi \rr$
functions
several values
of $A\!=\!\left\{ 2.5,5.0,7.5,10 \right\}$
are used
affecting
the maximal
elevation of
$\epsilon_h \lr \bx,t \rr$ field.
These $A$ values, result
in  $\epsilon_{h,max} \!=\! max \left[ \epsilon_h \lr \bx,t \rr \right] \!=\!\left\{ 3.5,6,8.5,11 \right\}$
elevations, see
results in \Fig{fig5}.

It is noticed,
\Figs{fig2}{fig4}
show only the results
of numerical simulation
in the two cases
$\epsilon_{h,max}/\epsilon_{h,b} \!=\! \left\{ 1,\,11 \right\}$.
When $\epsilon_{h,max}/\epsilon_{h,b}\!=\!1$
than $\epsilon_h \lr \bx,t \rr=const.$;
 in the
case $\epsilon_{h,max}/\epsilon_{h,b}\!=\!11$,   $\epsilon_h \lr \bx,t \rr$
is defined by \Eq{eq19}
with $A=10$.

\section{Discusion of the results}
\label{sec5}

\subsection{Conservation of mass}
\label{ssec51}

At first the mass conservation
of the present numerical method
is discussed, see \Fig{fig5}.
As in the present
model we are using
the functions $\alpha \lr  \psi \rr$
(in the equilibrium case when $\epsilon_h \lr \bx,t \rr\!=\!\epsilon_{h,b}\!=\!const.$)
and $\alpha_d \! \lr  \psi \rr$ (to account for variations in $\epsilon_h \lr \bx,t \rr$ field),
they both can be used
to compute conservation
of mass or area in
incompressible
two-phase flow.
In the statistical interpretation,
the mass conservation
is equivalent of the probability
(of finding phase chosen
to identify the macroscopic
interface $\gamma$) conservation.

In \Fig{fig5} history of
the mass variation
in the present test case
is presented.
This diagram is obtained using
integration of $\alpha \lr  \psi \rr$, $\alpha_d \! \lr  \psi \rr$
functions in the computational
domain there where $\psi  \left( \alpha \right)$
is reconstructed.
It is noted, the total mass/area
is slightly higher then $S_{ext}=\pi D^2/4$
used to normalize data in \Fig{fig5}.
This is because $\alpha \lr  \psi \rr$, $\alpha_d \! \lr  \psi \rr$
are the cumulative distribution functions
with the infinite support and
they are resolved  that
way as accurate as it is possible
with the double precision computer
arithmetic.
In both cases presented
in \Fig{fig5} the law of mass
conservation is satisfied
during the whole simulation.
In this test, the
coarse mesh is used
($\sim 26$ grid nodes per buble
diameter $D$ at $t=0\,[s]$),
the relatively exact
mass conservation is due to
the conservative form
of the transport \Eq{eq7}.

In the case of
the mass recording
obtained
using $\alpha \lr  \psi \rr$
in \Fig{fig5} (left)
one observes
how coupling
through continuity
equation with
the non-equilibrium solution
affects its conservation.
Counter intuitively,
the larger $\epsilon_h \lr \bx,t \rr$
field magnitude $\epsilon_{h,max}$ is,
the lower mass error
is obtained.

In \Fig{fig5} (right),
conservation
of $\alpha_d \! \lr  \psi \rr$
is depicted.
Herein, one notes interaction
with the variable characteristic
length scale field $\epsilon_h \lr \bx,t \rr$
causes
fluctuations of
the mass magnitude
in $\pm3\%$ range.
However,
due to coupling
with the conserved
$\alpha \lr  \psi \rr-\psi  \left( \alpha \right)$ functions
the total mass
is conserved as well.

The time moments
of oscillations peaks
in \Fig{fig5} (right)
can be compared
with the two-phase
flow field(s)
in \Figs{fig2}{fig4}.
One can observe,
these peaks
occur when the non-symmetrical
modes of $\alpha_d \! \lr  \psi \rr$ are induced
by $\epsilon_h \lr \bx,t \rr$ field.
This can be obseved
at times
$t=0.015,0.02\,[s]$,
e.g. compare results
in \Fig{fig5}
and \Fig{fig2}.

\subsection{Impact of non-equilibrium effects on two-phase flow}
\label{ssec52}

Impact of the characteristic
length scale field $\epsilon_h \lr \bx,t \rr$ on
the obtained results is presented
in \Figs{fig3}{fig4},
compare results in the top
and bottom rows.
One observes,
in the case $\epsilon_h \lr \bx,t \rr=\epsilon_{h,b}$
investigated two-phase flow
is symmetric.
In spite they
are there, black
contours of $\alpha_d \! \lr  \psi \rr$
field can not be seen
in \Fig{fig3} (top row)
because in this
case $\alpha_d \! \lr  \psi \rr\!=\!\alpha \lr  \psi \rr$.

This is not the case
in the bottom row of
\Fig{fig3}.
Herein, the difference
between $\alpha \lr  \psi \rr$ (white contours)
and $\alpha_d \! \lr  \psi \rr$ (black contours)
functions is clearly
visible.
As in the non-equilibrium case,
the material properties of
gas and liquid phases
are computed using
\Eqs{eq15}{eq16} with
$\alpha \lr  \psi \rr$ replaced by $\alpha_d \! \lr  \psi \rr$.
This locally changes
the mass and diffusive
fluxes used to solve
the Navier-Stokes equations
and couple
the pressure
and velocity
fields.
Therefore now,
the two-phase flow
scenario is
affected by $\alpha_d \! \lr  \psi \rr$
variations
causing
the flow
symmetry
to be broken.
This can be
attributed to
differences in
topological changes
on the left and
right side of
the gas bubble,
compare results
in the top and bottom
rows in \Figs{fig3}{fig4}
for example
at times $t=0.25,\,0.3\,[s]$.

In our opinion,
the best illustration
of mechanisms in
the statistical two-phase
flow model is given
by the joint probabilities
$\alpha \lr  \psi \rr \left( 1-\alpha \lr  \psi \rr \right)$
or
$\alpha_d \! \lr  \psi \rr \left( 1-\alpha_d \! \lr  \psi \rr \right)$
depicted in the top
and bottom rows in \Fig{fig4},
respectively.
Variations
of the joint probability
used to identify macroscopic
interface $\gamma$ as $\psi \left( \alpha=1/2 \right)=0$,
illustrate
the probability distribution
of the mesoscopic
interface $\Gamma$
residence.
In other words,
in regions with higher
joint probability, it is
more likely to find the
mesoscopic interface $\Gamma$
(we recall in the present
work $\Gamma$ is defined on the
molecular level).

\subsection{Difference between statistical and deterministic two-phase flow models}
\label{ssec53}

Finally,
let us  discuss
the main
difference
between
deterministic
and statistical \cite{mwaclawczyk11,waclawczyk2015,twacl15,twacl17,twacl21}  two-phase flow
models used in the one-fluid model
framework.

In the  deterministic description
of the interface,
the physical interpretation
is attributed only to flow
domains where the order parameter or
the phase indicator function
has value zero or one \cite{trygg11,kim12}.
According to models of
Gibbs \cite{gibbs1874} and
van der Waals \cite{waals1979},
these domains are occupied
by the homogeneous gas or liquid
phases \cite{lang15}.
However,
in
physical/chemical reality,
gases or liquids will never
be homogeneous
(gas or liquid phases containing
100\% of one type molecules
are almost never observed in every
day experience).
Moreover,
this constraint imposed on
the deterministic sharp interface models
results in the presumption
that
 domains where
 the numerical
 solution of the order parameter
 or phase indicator transport
 equation is between zero and one
 (smeared, not
sharply resolved regions)
are non-physical.
Similarly,
diffusive interface
in the phase field models
is considered to be thin
transition region
between the  two-phases
without clear  physical
interpretation as well.
Nevertheless,
exceptionally,
the position
of the interface
in its sharp and/or
diffusive deterministic
models
is typically located
between zero and one.
As a consequence,
the only possibility
to increase
the accuracy
of the numerical
predictions
with the sharp or
diffusive
interface model(s)
is
addition of
the new
grid points using
the adaptive
refinement
what impairs
performance
of the solver
(steam is putted
to the wrong,
non-physical wheel).

In the statistical model
of  intermittency region,
the cumulative distribution
function describing
the probability
of finding
mesoscopic interface
$\Gamma$
(or  probability
of finding the phase chosen to
identify $\Gamma$ and $\gamma$)
is
bounded $0\!<\!\alpha \lr  \psi \rr\!<\!1$
and  conserved
from its definition,
see mathematical
form of \Eq{eq7}.
In this interpretation,
the expected position of the
interface $\Gamma$ defined as
$\alpha \left( \psi\!=\!0 \right) \!=\! 1/2$
explains where the macroscopic
interface $\gamma$ is
located (also
in the deterministic
interface models).
Moreover, values
of $\alpha \lr  \psi \rr$ have
the physical interpretation
over the entire
range of $\alpha \lr  \psi \rr$
field magnitude.
Hence,
regions
where
the two-phase flow
is not resolved
due to lack of
the spatial
resolution are
promoted to
physical
interpretation
as well.
This feature of
the statistical model
of the intermittency region
is general and applies
equally to the equilibrium
and non-equilibrium flow
cases.
In order not
to be too vague,
let us analyze
the results
obtained
in the
present numerical
experiment
in the context
of above
explanations.

It is recalled,
$\alpha \lr  \psi \rr$ and $\alpha_d \! \lr  \psi \rr$
are probabilities of finding
the liquid phase,
$1\!-\!\alpha \lr  \psi \rr$ and $1\!-\!\alpha_d \! \lr  \psi \rr$
are probabilities of finding
the gas phase.
We note, at $t=0.3\,[s]$ (see \Figs{fig2}{fig4}) some
features of the rising gas bubble are
not well resolved.
Contours of $\alpha \lr  \psi \rr$ (white)
and $\alpha_d \! \lr  \psi \rr$ (black) capture only values
$\left\{ 0.5,0.99 \right\}$, see \Figs{fig2}{fig3}.
Unlike in the deterministic
interface models,
the present numerical solution,
although not fully resolved,
still carries
a physical
information.
Namely it
means, that
finding of
the  gas phase
in such regions
is not certain
as
 $1\!-\!\alpha \lr  \psi \rr \!<\! 0.99$
and this is why
$\alpha \lr  \psi \rr\!=\!0.01$ or  $\alpha_d \! \lr  \psi \rr\!=\!0.01$
contours can not be seen.
\section{Conclusions}
\label{Conclusions}

In the present paper,
the statistical model of
the intermittency region
between two phases \cite{twacl17,twacl21}
was coupled with
the incompressible
Navier-Stokes equation
solver
to study
differences between
the equilibrium and
non-equilibrium solution
of \Eq{eq7} in
the case of
rising gas bubble.
Results
of the present
numerical simulations
show,
inclusion of
recently
proposed
approximate
solution to \Eq{eq7}
permits to account
for non-equilibrium effects
in the
two-phase flow model
rooted in the one-fluid
framework.

Moreover, we
have discussed
differences
between deterministic
sharp and diffusive
interface models
and statistical
interface model.
It is concluded,
the main difference
between
these two descriptions
of two-phase flow
is that the results of
the statistical
two-phase flow
model can be interpreted
over the entire range
of the cumulative
distribution function
$\alpha \lr  \psi \rr$ values
also when
two-phase
flow is
not fully
resolved.

\bibliographystyle{spmpsci}