\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}
