# A Vortex Damping Outflow Forcing for Multiphase Flows with Sharp Interfacial
Jumps

Akash Dhruv  Mathematics and Computer Science Division,
Argonne National Laboratory, Lemont, IL, USA

###### Abstract

Outflow boundaries play an important role in multiphase fluid dynamics
simulations that involve transition between liquid and vapor phases. These
flows are dominated by low Weber numbers and a sharp jump in pressure,
velocity, and temperature. Inadequate treatment of these jumps at the outlet
generates undesirable fluid disturbances that propagate upstream and lead to
instabilities within the computational domain. To mitigate these disturbances,
we introduce a forcing term that can be applied to incompressible Navier-
Stokes equations to enforce stability in the numerical solution. The forcing
term acts as a damping mechanism to control vortices that are generated by
droplet/bubbles in multiphase flows, and is designed to be a general
formulation that can be coupled with a fixed pressure outflow boundary
condition to simulate a variety of multiphase flow problems. We demonstrate
its applicability to simulate pool and flow boiling problems, where bubble-
induced vortices during evaporation and condensation present a challenge at
the outflow. Validation and verification cases are chosen to quantify accuracy
and stability of the proposed method in comparison to established benchmarks
and reference solutions, along with detailed performance analysis for three-
dimensional simulations on leadership supercomputing platforms. Computational
experiments are performed using Flash-X, which is a composable open-source
software instrument designed for multiscale fluid dynamics simulations on
heterogeneous architectures.

###### keywords:

Multiphase flows, Outflow boundary condition, Pool boiling, Flow boiling,
Level-set method, Ghost fluid method, Adaptive mesh refinement

††journal: Journal of Computational Physics

##  1 Introduction

Outflow boundary conditions are used to enforce far-field behavior of flow
variables in computational fluid dynamics (CFD) simulations, and are typically
implemented to mitigate numerical instabilities caused by pressure
fluctuations and turbulence at the exit of a computational domain. Their
design and application are largely determined by nature of the problem,
numerical formulation, limitations of the computational box size, and their
effect on the overall accuracy and reliability of the solution. Present
literature covers a variety of boundary conditions for the simulation of
single and multiphase incompressible flows [14, 13, 3, 25, 20, 2], each with
their advantages and drawbacks when coupling pressure and velocity to enforce
momentum and continuity equations.

Within the context of a fractional step, staggered-grid, predictor-corrector
(projection) method for the incompressible Navier-Stokes equations, the
simplest implementation involves setting a fixed-pressure outflow, and forcing
velocity-gradients normal to the boundary to be zero. This approach implicitly
satisfies mass conservation through coupling between predictor, pressure
Poisson, and corrector equations that are solved in a sequence to enforce
momentum and continuity constraints; and is suitable for flows without an
explicitly defined inlet or with uniform pressure distribution at far-field
boundaries, as in large-scale atmospheric simulations.

For turbulent flows with an explicitly defined inlet, such as flow over an
aircraft geometry, a different approach is used. This involves setting a zero
pressure-gradient at the outflow, accompanied by solving a convection equation
for velocity at the boundary during the corrector step. The convective
velocity of the equation is set equal to the maximum or mean velocity of the
flow at the outflow boundary, which forces the vortices to exit the domain
without generating a backflow. Due to the nature of the implementation, that
is, replacing the solution of incompressible Navier-Stokes equations with a
convective transport equation, explicit mass conservation is enforced by
scaling the velocities at the outflow based on the mean velocity at the inlet.
This popular treatment was implemented by Vanella et al. for single-phase
fluid-structure interaction simulations within FLASH [38], a multiphysics
simulation software framework.

Outflow treatments in general influence the flow behavior upstream and
introduce bias in numerical solution. However, they serve has as a more
favorable alternative to using larger simulations domains which increase the
computational cost, and help mitigate numerical instabilities that arise from
non-physical backflow propagations. As a result, this topic has garnered
significant attention from researchers in various domains to design
formulations that reduce outflow-induced errors in the numerical solution. The
most sophisticated approach involves imposing a traction boundary condition
that is derived from the weak form of Navier-Stokes equations [35, 3], and
applied to modify the formulation of the Pressure Poisson Equation (PPE) for
the fractional-step projection method for incompressible Navier-Stokes
equations [20, 25, 3].

Extension of outflow boundary condition to multiphase flows requires special
considerations to satisfy mass conservation for individual phases, and
mitigate pressure and velocity fluctuations due to interfacial jumps. Traction
boundary conditions are excellent candidates for designing a general solution
to model outflow in multiphase flow simulations, however, the complexity lies
in designing appropriate traction force that can balance effect of various
multiphase dynamics. Many have investigated this problem in detail. Dong et
al. [14, 13] have made contributions in this area for flows involving
multiphase jets with explicitly defined inlet and outlet. Their formulation is
based on a phase field approach to interface tracking with a mollified forcing
for multiphase boundary conditions, which is applicable to the simulation of
flows where inertial force dominates surface tension, resulting in high Weber
numbers ($\text{We}>10^{3}$). Recently, Bozonnet et al. [3] made advances in
this area by designing a formulation that allows for stability in multiphase
flow simulations where backflow is a natural solution.

For flows involving phase transition, where inertial and surface tension
forces are in a balance ($\text{We}=1$) new challenges related to mass
conservation and sharp treatment of interfacial jumps have emerged. These
challenges particularly arise due to large surface tension forces at the
triple contact point of liquid, vapor and outflow boundary, which lead to
instability in the pressure solver, eventually affecting the divergence of
velocity and overall solution of momentum equations. In addition, local effect
of evaporation poses difficulties with explicit mass balance within the
domain, presenting a new research avenue to develop outflow boundary condition
for these types of problems.

Numerical simulations of multiphase flows with sharp interfacial jumps have
been conducted by many researchers in the last decade [16, 32, 22, 29, 8, 9,
31, 36], with contributions ranging from fundamental mathematical formulations
to production simulations that mimic experimental setups and provide new
insights into the physics of bubble/droplet interactions. Considerable effort
has also been made to model phase-transition. These efforts range from
developing numerical methods [33, 21, 24] and benchmarks [34] to performing
complex simulations to gain physical insights into the process of boiling
across a range of parametric space [37, 30]. In our previous work [8], we
extended FLASH to focus on the problem of pool boiling and designed an
adaptive mesh refinement (AMR)-based solver for high-fidelity simulations, and
performed an extensive numerical study to validate and quantify effects of
gravity on boiling heat transfer mechanisms [9]. The outflow boundary
condition consisted of a fixed pressure treatment with normal velocity
gradients set to zero at the boundary. To avoid numerical instabilities at the
outflow due to multiphase interactions described above, we introduced
artificial condensation for bubbles based on the treatment proposed by Sato et
al. [28] which leads to complete condensation from vapor to liquid within a
buffer region near the outflow. This treatment mimicked a type of
sponge/nudging layer that has been widely used to damp flows widely in oceanic
simulations [2]. As we tried to exercise our solver for a range of different
cases for pool and flow boiling simulations, spanning a parametric space of
different degrees of subcooling, wall superheat, and vapor-liquid properties,
numerical challenges started to appear. The buffer region was characterized by
its length and the strength of artificial subcooling. To completely condense
vapor bubbles into liquid, the length and strength of the buffer region had to
be adjusted on the basis of the degree of subcooling, gravity, and the rest of
the parametric space. This became heuristic and often tedious process when
deploying production simulations. Additionally, artificial condensation
introduced non-physical vortices that necessitated development of a more
formal treatment.

Extending traction boundary conditions for phase transition problems in our
AMR based solver poses several challenges. Existing multiphase formulations
are implemented using a variable coefficient PPE [3]. Within our framework, we
use a multigrid Poisson solver that would require building the coefficient
matrix associated with this PPE every time step as the multiphase system
evolves. For computational efficiency, we solve a constant coefficient PPE for
the variable-density system based on the methodology described in [12].
Extending multiphase traction boundary condition is not trivial for this
approach. Additionally, using homogeneous boundary conditions for the pressure
solver and formalizing the sponge/nudging layer based approach can serve as a
foundation for developing more robust traction boundary conditions in the
future. Therefore, in this work, we focus on treating outflow in multiphase
flows using a forcing term for velocity and pressure gradients. The
formulation is implemented for a sharp interface ghost fluid method (GFM) but
can be easily extended to other interface capturing methods with minor
modifications. The forcing term is used with a fixed pressure outflow
condition to mitigate numerical instabilities that arise in problems like pool
and flow boiling. The numerical framework is implemented within Flash-X [15],
a multiphysics simulation software instrument derived from FLASH, and is
actively being developed for performance portability and heterogeneous
computing.

The paper is structured as follows, Section 2 provides an overview of the
mathematical formulation and governing equations, followed by Section 3 which
provides details of numerical implementation and performance metrics. Section
4 discusses the outflow forcing, and finally results and conclusions are
provided in Sections 5 and 6, respectively.

##  2 Governing Equations

Figure 1 provides a schematic of pool and flow boiling simulations that serve
as fundamental tools to understand design of thermal cooling systems for
industrial, automotive, and electronic applications. Boiling occurs when a
liquid coolant undergoes evaporation on the surface of a solid heater,
resulting in formation of gas (vapor) bubbles which induce turbulence and
improve heat transfer efficiency. The dynamics of the bubbles is governed by
the balance between gravity ($g$), surface tension ($\sigma$) and evaporative
heat flux. In the simulation, the liquid-gas interface, $\Gamma$, is tracked
using a level-set, $\phi$, a signed distance function that is positive inside
the gas and negative inside the liquid. $\phi=0$, represents the implicit
location of $\Gamma$.

During flow boiling, the liquid coolant flows over the heater surface with a
bulk temperature, $T_{b}$ ($T_{bulk}$) and a velocity, $U_{b}$. The heater
acts as a no-slip wall with constant temperature, $T_{w}$ ($T_{wall}$), and
promotes the evaporation and formation of bubbles. The corresponding boundary
conditions for velocity ($u$), pressure ($P$) and temperature ($T$) are shown
in Figure 1a using tensor notation to represent derivatives and direction.

The departure of bubbles from the heater surface and their subsequent movement
through the liquid produce vortices that induce turbulence and circulatory
flow within the computational domain. Vortices are the result of interfacial
dynamics related to evaporation, condensation, and surface tension, which also
appear during pool boiling. Pool boiling is similar to flow boiling without
any liquid inflow. See the boundary conditions in Figure 1 b for more details.

The dynamics of boiling are modeled using a coupled system of momentum and
energy equations. Assuming an incompressible flow, the equations can be
described for each phase using tensor notations for three dimensional space.
The equations for liquid phase are written as,

$\frac{\partial u_{i}}{\partial
t}+u_{j}\partial_{j}u_{i}=-\frac{1}{\rho_{L}^{\prime}}\partial_{i}P+\frac{1}{\rho_{L}^{\prime}\>\text{Re}}\partial_{j}\Big{(}\mu_{L}^{\prime}\partial_{j}u_{i}+\mu_{L}^{\prime}\partial_{i}u_{j}\Big{)}+\frac{g_{i}}{\text{Fr}^{2}}$
(1a) $\frac{\partial T}{\partial
t}+u_{j}\partial_{j}T=\frac{1}{\rho_{L}^{\prime}C_{p_{L}}^{\prime}\text{Pe}}\partial_{j}\Big{(}k_{L}^{\prime}\partial_{j}T\Big{)}$
(1b)

And for the gas phase as,

$\frac{\partial u_{i}}{\partial
t}+u_{j}\partial_{j}u_{i}=-\frac{1}{\rho_{G}^{\prime}}\partial_{i}P+\frac{1}{\rho_{G}^{\prime}\>\text{Re}}\partial_{j}\Big{(}\mu_{G}^{\prime}\partial_{j}u_{i}+\mu_{G}^{\prime}\partial_{i}u_{j}\Big{)}+\frac{g_{i}}{\text{Fr}^{2}}$
(2a) $\frac{\partial T}{\partial
t}+u_{j}\partial_{j}T=\frac{1}{\rho_{G}^{\prime}C_{p_{G}}^{\prime}\text{Pe}}\partial_{j}\Big{(}k_{G}^{\prime}\partial_{j}T\Big{)}$
(2b)

Figure 1: Schematic of multiphase heat transfer problems: (a) Flow boiling (b)
Pool boiling, showing vortices generated by gas bubbles along with boundary
conditions for velocity ($u_{i}$), pressure ($P$), and temperature ($T$). The
level-set function ($\phi$) implicitly tracks the liquid-gas interface
($\Gamma$). The inset (c) provides details of the contact line, the pressure
jump, and the evaporative mass flux.

.

$u$ is the velocity, ${P}$ is the pressure, and ${T}$ is the temperature
everywhere in the domain. The equations are non-dimensionalized with reference
quantities from the liquid phase resulting in Reynolds number (Re)
$=\rho_{L}u_{0}l_{0}/\mu_{L}$, Prandtl number (Pr) $=\mu_{L}C_{p_{L}}/k_{L}$,
Froude number (Fr) $=u_{0}/\sqrt{gl_{0}}$, and Peclet number (Pe)
$=\text{Re}\>\text{Pr}$. $u_{0}$ and $l_{0}$ are reference velocity and length
scales, and $g$ is acceleration due to gravity. The symbols $\rho_{L}$,
$\mu_{L}$, $C_{p_{L}}$ and $k_{L}$ represent density, viscosity, specific heat
and thermal conductivity of the liquid. For boiling problems, the reference
length scale is set to the capillary length
$l_{0}=\sqrt{\sigma/(\rho_{L}-\rho_{G})\>g}$ (where $\sigma$ is the surface
tension), and the velocity scale is the terminal velocity
$u_{0}=\sqrt{g\>l_{0}}$. The reference temperature scale is given by
$(T-T_{bulk})/\Delta T$, where $\Delta T=T_{wall}-T_{bulk}$. Additionally, the
non-dimensional fluid properties in each phase are given by,

$y_{L}^{\prime}=1,\quad y_{G}^{\prime}=\frac{y_{G}}{y_{L}}$ (3)

where $y$ corresponds to $\rho$ (density), $\mu$ (viscosity), $C_{p}$
(specific heat), and $k$ (thermal conductivity).

The momentum equations strictly follow incompressible formulation reducing the
viscous term in Equation 1a and 1b to,

$\partial_{j}\Big{(}\mu^{\prime}\partial_{j}u_{i}+\mu^{\prime}\partial_{i}u_{j}\Big{)}=\partial_{j}\Big{(}\mu^{\prime}\partial_{j}u_{i}\Big{)}+\mu^{\prime}\partial_{i}\Big{(}\partial_{j}u_{j}\Big{)}=\partial_{j}\Big{(}\mu^{\prime}\partial_{j}u_{i}\Big{)}$
(4)

where $\mu^{\prime}$ is $\mu_{L}^{\prime}$ or $\mu_{G}^{\prime}$ depending on
the phase.

This reduction is based on the constraint that divergence of velocity is zero
in both phases. The effect of evaporation, which leads to non-zero divergence
near the interface, is handled separately using interfacial jump conditions
discussed further in the text. Unless specified otherwise, physical quantities
are assumed to be non-dimensional in the following discussions. Additionally,
Einstein’s convention for summation is observed for repeating tensor indices.
In special cases where this convention is violated—specifically, when
repeating indices appear more than twice in a term—an explicit summation
operator is provided.

The continuity equation is given by,

$\partial_{j}u_{j}=-\dot{m}\>n_{\Gamma_{j}}\partial_{j}\frac{1}{\rho^{\prime}}\>\Biggr{|}_{\>\Gamma}$
(5)

where,
$n_{\Gamma_{i}}=\partial_{i}\phi/\sqrt{\partial_{j}\phi\partial_{j}\phi}$, is
a unit vector normal to the liquid-gas interface, and, $\dot{m}$, is the
evaporative mass flux. It is important to emphasize that the fluid flow is
fully incompressible within each phase, and the right-hand side of Equation 5
accounts for the mass balance due to the change in phase from liquid to gas
and vice versa. Furthermore,
$\dot{m}\>n_{\Gamma_{j}}\partial_{j}\frac{1}{\rho^{\prime}}$, is zero
everywhere except at the interface.

Figure 1c shows a schematic of the heat flux from liquid to gas,
$n_{\Gamma_{j}}\partial_{j}T_{L}$ and the heat flux from gas to liquid,
$n_{\Gamma_{j}}\partial_{j}T_{G}$, normal to the liquid-gas interface,
$\Gamma$. The difference between these heat fluxes determines the mass
transfer from one phase to the other. The value of ${\dot{m}}$ at the
interface is calculated as,

$\dot{m}=\frac{\text{St}}{\text{Re}\>\text{Pr}}\Bigg{(}{n_{\Gamma_{j}}k_{L}^{\prime}\>\partial_{j}T_{L}\>\bigr{|}}_{\>\Gamma}-n_{\Gamma_{j}}k_{G}^{\prime}\>{\partial_{j}T_{G}\>\bigr{|}}_{\>\Gamma}\Bigg{)}$
(6)

where St is the Stefan number given by, $\text{St}=C_{p_{L}}\Delta T/Q_{l}$,
with ${Q_{l}}$ as latent heat of evaporation. The heat flux from the liquid
region, $\partial_{j}T_{L}$, contributes to evaporation, and the heat flux
from the gas region, $\partial_{j}T_{G}$, accounts for condensation.

The level-set function representing the interface is computed using the
convection equation,

$\frac{\partial\phi}{\partial t}+u_{\Gamma_{j}}\partial_{j}\phi=0,$ (7)

where, $u_{\Gamma_{i}}=u_{i}+(\dot{m}/\rho^{\prime})n_{\Gamma_{i}}$, is the
interface velocity. The convection of level-set is accompanied by a selective
reinitialization technique to mitigate diffusive errors that are generated
from numerical descretization of the convection term. The detailed approach
for this is discussed in [8, 5]

Reorganizing equation for,
$u_{\Gamma_{i}}=u_{i}+(\dot{m}/\rho^{\prime})n_{\Gamma_{i}}$, results in the
following expression for jump in velocity normal to the liquid-gas interface:

$[u]_{\Gamma_{i}}=u_{G_{i}}-u_{L_{i}}=\dot{m}\>n_{\Gamma_{i}}\Bigg{[}\frac{1}{\rho_{G}^{\prime}}-\frac{1}{\rho_{L}^{\prime}}\Bigg{]}$
(8)

The mass flux, $\dot{m}$, surface tension, $\sigma$, and viscous stresses
contribute towards a similar jump in pressure: where, ${\kappa}$ is the
interface curvature. The non-dimensional form for pressure jump using,
${\text{We}}=\rho_{l}u_{0}^{2}l_{0}/\sigma$, is consequently given by:

$[P]_{\Gamma}=P_{g}-P_{l}=\frac{\kappa}{\text{We}}-\Bigg{[}\frac{1}{\rho_{G}^{\prime}}-\frac{1}{\rho_{L}^{\prime}}\Bigg{]}\dot{m}^{2}$
(9)

The effect of viscous jump is assumed to negligible in this formulation due to
the smeared treatment of viscosity near the interface described in the
following section. This assumption fits with the formulation in [19].

Finally, the boundary condition for temperature at the liquid-gas interface is
given by,

$T_{\Gamma}=T_{sat}$ (10)

where, $T_{sat}$ corresponds to saturation temperature.

The jump conditions in Eqs 8, 9, 10 are modeled using a GFM [16, 8, 5].

##  3 Numerical Formulation

For the purposes of numerical formulation Eqs 1, 2, and 8-10 are merged into a
single-fluid approach, which has variable properties in space (either liquid
or gas) depending on the location of the interface. These equations are
discretized on a block-structured Cartesian grid in the full domain containing
both phases. The merged equations are written as follows,

$\frac{\partial u_{i}}{\partial
t}+u_{j}\partial_{j}u_{i}=\frac{1}{\rho^{\prime}}\Big{(}-\partial_{i}P+\text{H}_{\Gamma}\text{S}^{\Gamma}_{\partial_{i}P}+\text{S}^{\text{O}}_{\partial_{i}P})+\frac{1}{\rho^{\prime}\>\text{Re}}\partial_{j}\Big{(}\mu^{\prime}\partial_{j}u_{i}+\mu^{\prime}\partial_{i}u_{j}\Big{)}+\frac{g_{i}}{\text{Fr}^{2}}+\text{H}_{\Gamma}\text{S}^{\Gamma}_{u_{i}}+\text{S}^{\text{O}}_{u_{i}}$
(11a) $\frac{\partial T}{\partial
t}+u_{j}\partial_{j}T=\frac{1}{{\text{Pe}}}\partial_{j}\Big{(}\alpha^{\prime}\partial_{j}T\Big{)}+\text{H}_{\Gamma}\text{S}^{\Gamma}_{T}$
(11b)

where thermal diffusivity
$\alpha^{\prime}=k^{\prime}/(\rho^{\prime}C_{p}^{\prime})$ is used to
consolidate the the diffusion term for the energy equation. The effects of
interfacial jumps at $\Gamma$ ($\phi=0$) in Eqs 8, 9, and 10 are given by the
source terms $\text{S}^{\Gamma}_{u_{i}}$, $\text{S}^{\Gamma}_{\partial_{i}P}$,
and $\text{S}^{\Gamma}_{T}$ respectively. These terms are only relevant when
$\text{H}_{\Gamma}=1$. Based on $\phi$, the value of $\text{H}_{\Gamma}$ is 1
only in the cells near $\phi=0$ and is zero otherwise. The source terms,
$\text{S}^{\text{O}}_{u_{i}}$ and $\text{S}^{\text{O}}_{\partial_{i}P}$,
correspond to the treatment of the outflow which is the main focus of this
work.

A finite difference spatial and temporal discretization is implemented for Eqs
11 \- 10, representing a full system for multiphase fluid dynamics equations
solution. The vectors along with $\rho^{\prime}$ are defined at face-centers,
and the scalars along with remaining properties are defined at cell-centers.
The non-dimensional fluid properties at cell-centers are given as,

$\mu^{\prime}=\mu_{L}^{\prime}\>(1-\text{H}_{\epsilon})+\mu_{G}^{\prime}\>\text{H}_{\epsilon},\quad\alpha^{\prime}=\alpha_{L}^{\prime}\>(1-\text{H}_{\epsilon})+\alpha_{G}^{\prime}\>\text{H}_{\epsilon}$
(12)

where $\text{H}_{\epsilon}$ is a smeared Heaviside function,

$\text{H}_{\epsilon}=\begin{cases}0,&\text{if}\>\phi<-\epsilon\\\
\frac{1}{2}[1+\frac{\phi}{\epsilon}+\frac{1}{\pi}\text{sin}(\frac{\pi\phi}{\epsilon})],&\text{if}\>|\phi|\leq\epsilon\\\
1,&\text{if}\>\phi>\epsilon\\\ \end{cases}$ (13)

and $\epsilon$ is a tolerance that is chosen as $1.5$ times the grid spacing.

The face-centered density is defined as,

$\rho^{\prime}=\rho_{L}^{\prime}\>(1-\text{H})+\rho_{G}^{\prime}\>\text{H},\quad$
(14)

where H is a sharp Heaviside function,

$\text{H}=\begin{cases}0,&\text{if}\>\phi<0\\\ 1,&\text{if}\>\phi\geq 0\\\
\end{cases}$ (15)

Additionally, a special smearing strategy is applied for face-centered
density. Consider Figure 2, which describes the computational stencil in
direction $x_{1}$ with cell-centers ($I_{C}$) and faces ($I_{F}$). Near the
liquid gas interface, $\Gamma$, the density is smeared using,

$\rho^{\prime}\Bigr{|}_{I_{F}+1}=\rho_{L}^{\prime}\>(1-\beta)+\rho_{G}^{\prime}\>\beta$
(16)

where,
$\beta={|\phi^{n+1}_{I_{c}}|}\Big{/}\Big{(}{|\phi^{n+1}_{I_{c}}|+|\phi^{n+1}_{I_{c}+1}|}\Big{)}$.

Figure 2: Staggered computational grid highlighting faces ($\times$) and cell-
centers ($\bullet$), $I_{C}$ denotes the cell-centers and $I_{F}$ if face-
center. The liquid-gas interface, $\Gamma$, separates the two phases.

The advection terms are calculated using a fifth-order Hamilton-Jacobi
weighted essentially non-oscillatory (WENO) method described in [18, 4], and
the diffusion terms are calculated using a second-order central difference
discretization. Temporal integration of the convective and diffusion terms is
implemented using a second-order Adams-Bashforth scheme.

The momentum Eq 11a follows a fractional-step predictor-corrector approach to
enforce continuity in Eq 5. The multistep method is given as follows,

$\frac{u^{*}_{i}-u^{n}_{i}}{\Delta
t}=\frac{3}{2}\Theta^{n}_{i}-\frac{1}{2}\Theta^{n-1}_{i}+\frac{g_{i}}{\text{Fr}^{2}}+\text{S}^{\Gamma}_{u_{i}}+\text{S}^{\text{O}}_{u_{i}}$
(17)

Eq 17 is the predictor step that does not include contribution from pressure
or surface tension forces. $\Theta_{i}$, represents the contribution of
advection and diffusion terms, and $n$ denotes the solution in time as the
simulation progresses with a time step, $\Delta t$. Assuming,
$\zeta^{n+1}_{i}=\partial_{i}P^{n+1}-\text{S}^{\Gamma}_{\partial_{i}P}-\text{S}^{\text{O}}_{\partial_{i}P}$,
the corrector equation becomes,

$\frac{u^{n+1}_{i}-u^{*}_{i}}{\Delta
t}=-\frac{1}{\rho_{G}^{\prime}}\zeta^{n+1}_{i}-\frac{1}{\rho^{\prime}}(2\zeta^{n}_{i}-\zeta^{n-1}_{i})+\frac{1}{\rho_{G}^{\prime}}(2\zeta^{n}_{i}-\zeta^{n-1}_{i})$
(18a) $\frac{u^{n+1}_{i}-u^{*}_{i}}{\Delta
t}=-\frac{1}{\rho_{G}^{\prime}}\partial_{i}P^{n+1}+\frac{1}{\rho_{G}^{\prime}}(\text{S}^{\Gamma}_{\partial_{i}P}+\text{S}^{\text{O}}_{\partial_{i}P})+\Bigg{(}\frac{1}{\rho_{G}^{\prime}}-\frac{1}{\rho^{\prime}}\Bigg{)}(2\zeta^{n}_{i}-\zeta^{n-1}_{i})$
(18b)

The splitting of the pressure term is implemented to convert the variable
coefficient Poisson equation into a configuration with constant coefficients.
This splitting is based on the formulation described in [12], which separates
the pressure term into a constant term (the unknown) and a variable term that
utilizes linear extrapolation from time steps $n$ and $n-1$. Consequently, the
divergence of the corrector step in Equation 18b leads to the following
pressure equation,

$\partial_{j}\partial_{j}P^{n+1}=\frac{\rho_{G}^{\prime}}{\Delta
t}\Bigg{(}\partial_{j}u^{*}_{j}-\partial_{j}u^{n+1}_{j}\Bigg{)}+\sum_{j=1}^{3}\partial_{j}\Bigg{(}1-\frac{\rho_{G}^{\prime}}{\rho^{\prime}}\Bigg{)}(2\zeta^{n}_{j}-\zeta^{n-1}_{j})\\\
+\partial_{j}\text{S}^{\Gamma}_{\partial_{j}P}+\partial_{j}\text{S}^{\text{O}}_{\partial_{j}P}$
(19)

$\partial_{j}u^{n+1}_{j}$ is calculated using Equation 5. In reference to
Figure 2, the source terms $\text{S}^{\Gamma}_{u_{1}}$ and
$\text{S}^{\Gamma}_{\partial_{1}P}$ are given as,

$\text{S}^{\Gamma}_{u_{1}}\Bigr{|}_{I_{F}+1}=\frac{1}{\rho^{\prime}\>\text{Re}}\partial_{1}\Big{[}\mu^{\prime}\partial_{1}(\hat{u}_{1}-u_{1})\Big{]}\>\Biggr{|}_{I_{F}+1}$
(20a)
$\text{S}^{\Gamma}_{\partial_{1}P}\Bigr{|}_{I_{F}+1}=\frac{\kappa^{n+1}_{I_{c}}(1-\theta)+\kappa^{n+1}_{I_{c}+1}\theta}{\text{We}\>\Delta
x_{1}}-\Bigg{(}\frac{1}{\rho_{G}^{\prime}}-\frac{1}{\rho_{L}^{\prime}}\Bigg{)}\frac{\dot{m^{2}}^{n+1}_{I_{c}}(1-\theta)+\dot{m^{2}}^{n+1}_{I_{c}+1}\theta}{\Delta
x_{1}}$ (20b)

where curvtature, $\kappa=\partial_{j}n_{\Gamma_{j}}$, and $\hat{u}_{1}$ is
the ghost fluid velocity that is applied to the viscous term in Equation 11a.
Using the definition of interface velocity in Equation 8,
$\hat{u}_{1}=u_{1}+(\dot{m}/\rho^{\prime})n_{\Gamma_{1}}$, modifies Equation
20a to,

$\text{S}^{\Gamma}_{u_{1}}\Bigr{|}_{I_{F}+1}=\frac{1}{\rho^{\prime}\>\text{Re}}\partial_{1}\Bigg{(}\frac{\mu^{\prime}\dot{m}}{\rho^{\prime}}\partial_{1}n_{\Gamma_{1}}\Bigg{)}\>\Biggr{|}_{I_{F}+1}$
(21)

The corresponding discretization of the energy and level-set equations is
discussed at length in [8, 5], along with the verification of the overall
second-order accuracy of the multiphase solver. Calculation of the forcing
term $\text{S}^{\Gamma}_{T}$ is similar to formulation of
$\text{S}^{\Gamma}_{u_{1}}$ for the stencil in Figure 2,

$\text{S}^{\Gamma}_{T}\Bigr{|}_{I_{C+1}}=\frac{1}{\text{Pe}}\partial_{1}\Big{[}{\alpha^{\prime}}\partial_{1}(\hat{T}-T)\Big{]}\>\Biggr{|}_{I_{C+1}}$
(22)

In this case, $\hat{T}$ is the ghost fluid temperature which is applied to the
diffusion term in Equation 11b and calculated as,

$\hat{T}_{I_{C+1}}=T_{I_{C}}+\Delta
x_{1}\frac{T_{\Gamma}-T_{I_{C}}}{\theta\Delta x_{1}}$ (23)

where $T_{\Gamma}=T_{sat}$ from Equation 10. Ghost fluid boundary conditions
are only applied to diffusion terms, while the WENO advection scheme
effectively preserves the sharp jump in properties and does not require
special treatment.

Block-structured, octree-based adaptive mesh refinement (AMR) is used to
selectively refine the grid in the regions of interest. The default AMR
backend within FLASH was implemented using Paramesh along with a customized
multigrid Poisson solver that leverages HYPRE to calculate the solution at the
coarsest refinement level. Flash-X, on the other hand has AMReX as an
additional option for the management of the AMR network. AMReX [39], which is
natively a patch-based library, can also mimic the octree mode when needed. We
use the octree mode of AMReX for a direct comparison of the performance of the
two AMR backends and their multigrid solvers in multiphase fluid dynamics
applications within Flash-X. Figure 3a provides a schematic of AMR around the
liquid-gas interface located at $\phi=0$. The refinement criteria is set to
satisfy the condition $\lvert\phi\rvert\leq\epsilon$, where $\epsilon$ is the
tolerance that sets a range for $\phi$ to provide thresholds for refinement
and derefinement.

Figure 3: (a) Schematic of a multi-level block structured AMR grid around the
liquid-gas interface, $\phi=0$, (b) Distribution of scaled pressure,
$\frac{P-P_{low}}{P_{high}-P_{low}}$, highlighting the jump in value between
phases, and (c) Weak scaling of pressure Poisson solver on Summit [1]

Figure 3b shows an example of a scaled pressure distribution for a pool
boiling problem that highlights the sharp jump in values that are effectively
captured by our numerical implementation. The corresponding weak scaling study
of CPU-only performance of the Poisson solver is shown in Figure 3c, which
clearly demonstrates an order of magnitude improvement in performance with
AMReX and its native multigrid solver in comparison to implementation with
Paramesh + HYPRE.

##  4 Outflow Forcing

The goal of proposed outflow forcing is to improve on previous implementations
[8, 28] and provide a more concrete approach that avoids the need for
introducing artificial condensation within the computational domain.
Artificial condensation, as discussed previously, is a non-physical artifact
which was implemented to ensure existence of only liquid phase at the outflow
boundary. Figure 4 describes this behavior in time for a saturated nucleate
boiling problem. Even though this approach helps with the multiphase
instabilities at the boundary, introduction of artificial condensation affects
the overall quality of the numerical simulation over time. In the present
approach, we focus on relaxing the flow to a steady state instead by deriving
forcing terms that have physical meaning within the context of multiphase
flows.

Figure 4: Effect of artificial condensation on saturated nucleate boiling.
Physical condensation does no exist since both the bubble and bulk liquid are
at the saturation temperature, however, outflow forcing in [8, 28] introduces
artificial condensation which is non-physical ($t_{1}<t_{2}<t_{3})$.

Within our numerical formulation, vectors and scalars are spatially
represented on a staggered grid. Vectors such as velocity are located on faces
($\times$), and scalars such as pressure are located at the cell centers
($\bullet$) as shown in the stencil in Figure 5. The implementation of a
staggered grid enables a tight coupling between the pressure gradient and the
velocity during the corrector step in Equation 18 which translates into a more
robust implementation of the Poisson equation for pressure in Equation 19. The
predictor-corrector Equation 17 and 18 at the outflow have the following form,

$u^{*}_{1}\Bigr{|}_{B}=u^{n}_{1}\Bigr{|}_{B}+\Delta
t\Bigg{(}\frac{3}{2}\Theta^{n}_{1}-\frac{1}{2}\Theta^{n-1}_{1}+\text{S}^{\Gamma}_{u_{1}}+\text{S}^{\text{O}}_{u_{1}}\Bigg{)}\>\Biggr{|}_{B}+\Delta
t\frac{g_{1}}{\text{Fr}^{2}}$ (24)
$\begin{split}u^{n+1}_{1}\Bigr{|}_{B}=u^{*}_{1}\Bigr{|}_{B}+\mathchoice{\leavevmode\hbox
to488.15pt{\vbox to54pt{\pgfpicture\makeatletter\hbox{\hskip
244.07611pt\lower-34.5pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-244.07611pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\displaystyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{blue!17}$\displaystyle\mathstrut\Delta
t\Bigg{(}\frac{1}{\rho_{G}^{\prime}}(\text{S}^{\Gamma}_{\partial_{1}P}+\text{S}^{\text{O}}_{\partial_{1}P})+\Bigg{(}\frac{1}{\rho_{G}^{\prime}}-\frac{1}{\rho^{\prime}}\Bigg{)}(2\zeta^{n}_{1}-\zeta^{n-1}_{1})\Bigg{)}\>\Biggr{|}_{B}$}}{\hbox{\pagecolor{blue!17}$\textstyle\mathstrut\Delta
t\Bigg{(}\frac{1}{\rho_{G}^{\prime}}(\text{S}^{\Gamma}_{\partial_{1}P}+\text{S}^{\text{O}}_{\partial_{1}P})+\Bigg{(}\frac{1}{\rho_{G}^{\prime}}-\frac{1}{\rho^{\prime}}\Bigg{)}(2\zeta^{n}_{1}-\zeta^{n-1}_{1})\Bigg{)}\>\Biggr{|}_{B}$}}{\hbox{\pagecolor{blue!17}$\scriptstyle\mathstrut\Delta
t\Bigg{(}\frac{1}{\rho_{G}^{\prime}}(\text{S}^{\Gamma}_{\partial_{1}P}+\text{S}^{\text{O}}_{\partial_{1}P})+\Bigg{(}\frac{1}{\rho_{G}^{\prime}}-\frac{1}{\rho^{\prime}}\Bigg{)}(2\zeta^{n}_{1}-\zeta^{n-1}_{1})\Bigg{)}\>\Biggr{|}_{B}$}}{\hbox{\pagecolor{blue!17}$\scriptscriptstyle\mathstrut\Delta
t\Bigg{(}\frac{1}{\rho_{G}^{\prime}}(\text{S}^{\Gamma}_{\partial_{1}P}+\text{S}^{\text{O}}_{\partial_{1}P})+\Bigg{(}\frac{1}{\rho_{G}^{\prime}}-\frac{1}{\rho^{\prime}}\Bigg{)}(2\zeta^{n}_{1}-\zeta^{n-1}_{1})\Bigg{)}\>\Biggr{|}_{B}$}}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to488.15pt{\vbox to54pt{\pgfpicture\makeatletter\hbox{\hskip
244.07611pt\lower-34.5pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-244.07611pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\textstyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{blue!17}$\displaystyle\mathstrut\Delta
t\Bigg{(}\frac{1}{\rho_{G}^{\prime}}(\text{S}^{\Gamma}_{\partial_{1}P}+\text{S}^{\text{O}}_{\partial_{1}P})+\Bigg{(}\frac{1}{\rho_{G}^{\prime}}-\frac{1}{\rho^{\prime}}\Bigg{)}(2\zeta^{n}_{1}-\zeta^{n-1}_{1})\Bigg{)}\>\Biggr{|}_{B}$}}{\hbox{\pagecolor{blue!17}$\textstyle\mathstrut\Delta
t\Bigg{(}\frac{1}{\rho_{G}^{\prime}}(\text{S}^{\Gamma}_{\partial_{1}P}+\text{S}^{\text{O}}_{\partial_{1}P})+\Bigg{(}\frac{1}{\rho_{G}^{\prime}}-\frac{1}{\rho^{\prime}}\Bigg{)}(2\zeta^{n}_{1}-\zeta^{n-1}_{1})\Bigg{)}\>\Biggr{|}_{B}$}}{\hbox{\pagecolor{blue!17}$\scriptstyle\mathstrut\Delta
t\Bigg{(}\frac{1}{\rho_{G}^{\prime}}(\text{S}^{\Gamma}_{\partial_{1}P}+\text{S}^{\text{O}}_{\partial_{1}P})+\Bigg{(}\frac{1}{\rho_{G}^{\prime}}-\frac{1}{\rho^{\prime}}\Bigg{)}(2\zeta^{n}_{1}-\zeta^{n-1}_{1})\Bigg{)}\>\Biggr{|}_{B}$}}{\hbox{\pagecolor{blue!17}$\scriptscriptstyle\mathstrut\Delta
t\Bigg{(}\frac{1}{\rho_{G}^{\prime}}(\text{S}^{\Gamma}_{\partial_{1}P}+\text{S}^{\text{O}}_{\partial_{1}P})+\Bigg{(}\frac{1}{\rho_{G}^{\prime}}-\frac{1}{\rho^{\prime}}\Bigg{)}(2\zeta^{n}_{1}-\zeta^{n-1}_{1})\Bigg{)}\>\Biggr{|}_{B}$}}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to493.53pt{\vbox to27pt{\pgfpicture\makeatletter\hbox{\hskip
246.76604pt\lower-7.5pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-246.76604pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\scriptstyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{blue!17}$\displaystyle\mathstrut\Delta
t\Bigg{(}\frac{1}{\rho_{G}^{\prime}}(\text{S}^{\Gamma}_{\partial_{1}P}+\text{S}^{\text{O}}_{\partial_{1}P})+\Bigg{(}\frac{1}{\rho_{G}^{\prime}}-\frac{1}{\rho^{\prime}}\Bigg{)}(2\zeta^{n}_{1}-\zeta^{n-1}_{1})\Bigg{)}\>\Biggr{|}_{B}$}}{\hbox{\pagecolor{blue!17}$\textstyle\mathstrut\Delta
t\Bigg{(}\frac{1}{\rho_{G}^{\prime}}(\text{S}^{\Gamma}_{\partial_{1}P}+\text{S}^{\text{O}}_{\partial_{1}P})+\Bigg{(}\frac{1}{\rho_{G}^{\prime}}-\frac{1}{\rho^{\prime}}\Bigg{)}(2\zeta^{n}_{1}-\zeta^{n-1}_{1})\Bigg{)}\>\Biggr{|}_{B}$}}{\hbox{\pagecolor{blue!17}$\scriptstyle\mathstrut\Delta
t\Bigg{(}\frac{1}{\rho_{G}^{\prime}}(\text{S}^{\Gamma}_{\partial_{1}P}+\text{S}^{\text{O}}_{\partial_{1}P})+\Bigg{(}\frac{1}{\rho_{G}^{\prime}}-\frac{1}{\rho^{\prime}}\Bigg{)}(2\zeta^{n}_{1}-\zeta^{n-1}_{1})\Bigg{)}\>\Biggr{|}_{B}$}}{\hbox{\pagecolor{blue!17}$\scriptscriptstyle\mathstrut\Delta
t\Bigg{(}\frac{1}{\rho_{G}^{\prime}}(\text{S}^{\Gamma}_{\partial_{1}P}+\text{S}^{\text{O}}_{\partial_{1}P})+\Bigg{(}\frac{1}{\rho_{G}^{\prime}}-\frac{1}{\rho^{\prime}}\Bigg{)}(2\zeta^{n}_{1}-\zeta^{n-1}_{1})\Bigg{)}\>\Biggr{|}_{B}$}}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to414.62pt{\vbox to26pt{\pgfpicture\makeatletter\hbox{\hskip
207.30896pt\lower-6.5pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-207.30896pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\scriptscriptstyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{blue!17}$\displaystyle\mathstrut\Delta
t\Bigg{(}\frac{1}{\rho_{G}^{\prime}}(\text{S}^{\Gamma}_{\partial_{1}P}+\text{S}^{\text{O}}_{\partial_{1}P})+\Bigg{(}\frac{1}{\rho_{G}^{\prime}}-\frac{1}{\rho^{\prime}}\Bigg{)}(2\zeta^{n}_{1}-\zeta^{n-1}_{1})\Bigg{)}\>\Biggr{|}_{B}$}}{\hbox{\pagecolor{blue!17}$\textstyle\mathstrut\Delta
t\Bigg{(}\frac{1}{\rho_{G}^{\prime}}(\text{S}^{\Gamma}_{\partial_{1}P}+\text{S}^{\text{O}}_{\partial_{1}P})+\Bigg{(}\frac{1}{\rho_{G}^{\prime}}-\frac{1}{\rho^{\prime}}\Bigg{)}(2\zeta^{n}_{1}-\zeta^{n-1}_{1})\Bigg{)}\>\Biggr{|}_{B}$}}{\hbox{\pagecolor{blue!17}$\scriptstyle\mathstrut\Delta
t\Bigg{(}\frac{1}{\rho_{G}^{\prime}}(\text{S}^{\Gamma}_{\partial_{1}P}+\text{S}^{\text{O}}_{\partial_{1}P})+\Bigg{(}\frac{1}{\rho_{G}^{\prime}}-\frac{1}{\rho^{\prime}}\Bigg{)}(2\zeta^{n}_{1}-\zeta^{n-1}_{1})\Bigg{)}\>\Biggr{|}_{B}$}}{\hbox{\pagecolor{blue!17}$\scriptscriptstyle\mathstrut\Delta
t\Bigg{(}\frac{1}{\rho_{G}^{\prime}}(\text{S}^{\Gamma}_{\partial_{1}P}+\text{S}^{\text{O}}_{\partial_{1}P})+\Bigg{(}\frac{1}{\rho_{G}^{\prime}}-\frac{1}{\rho^{\prime}}\Bigg{)}(2\zeta^{n}_{1}-\zeta^{n-1}_{1})\Bigg{)}\>\Biggr{|}_{B}$}}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}-\mathchoice{\leavevmode\hbox
to241.38pt{\vbox to20.28pt{\pgfpicture\makeatletter\hbox{\hskip
120.68904pt\lower-4.31888pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-120.68904pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\displaystyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{teal!17}$\displaystyle\mathstrut\Delta
t\frac{1}{\rho_{G}^{\prime}}\frac{P^{n+1}|_{G_{C}}-P^{n+1}|_{I_{C}}}{\Delta
x_{1}}$}}{\hbox{\pagecolor{teal!17}$\textstyle\mathstrut\Delta
t\frac{1}{\rho_{G}^{\prime}}\frac{P^{n+1}|_{G_{C}}-P^{n+1}|_{I_{C}}}{\Delta
x_{1}}$}}{\hbox{\pagecolor{teal!17}$\scriptstyle\mathstrut\Delta
t\frac{1}{\rho_{G}^{\prime}}\frac{P^{n+1}|_{G_{C}}-P^{n+1}|_{I_{C}}}{\Delta
x_{1}}$}}{\hbox{\pagecolor{teal!17}$\scriptscriptstyle\mathstrut\Delta
t\frac{1}{\rho_{G}^{\prime}}\frac{P^{n+1}|_{G_{C}}-P^{n+1}|_{I_{C}}}{\Delta
x_{1}}$}}$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to241.38pt{\vbox to20.28pt{\pgfpicture\makeatletter\hbox{\hskip
120.68904pt\lower-4.31888pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-120.68904pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\textstyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{teal!17}$\displaystyle\mathstrut\Delta
t\frac{1}{\rho_{G}^{\prime}}\frac{P^{n+1}|_{G_{C}}-P^{n+1}|_{I_{C}}}{\Delta
x_{1}}$}}{\hbox{\pagecolor{teal!17}$\textstyle\mathstrut\Delta
t\frac{1}{\rho_{G}^{\prime}}\frac{P^{n+1}|_{G_{C}}-P^{n+1}|_{I_{C}}}{\Delta
x_{1}}$}}{\hbox{\pagecolor{teal!17}$\scriptstyle\mathstrut\Delta
t\frac{1}{\rho_{G}^{\prime}}\frac{P^{n+1}|_{G_{C}}-P^{n+1}|_{I_{C}}}{\Delta
x_{1}}$}}{\hbox{\pagecolor{teal!17}$\scriptscriptstyle\mathstrut\Delta
t\frac{1}{\rho_{G}^{\prime}}\frac{P^{n+1}|_{G_{C}}-P^{n+1}|_{I_{C}}}{\Delta
x_{1}}$}}$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to168.96pt{\vbox to14.2pt{\pgfpicture\makeatletter\hbox{\hskip
84.48225pt\lower-3.02321pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-84.48225pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\scriptstyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{teal!17}$\displaystyle\mathstrut\Delta
t\frac{1}{\rho_{G}^{\prime}}\frac{P^{n+1}|_{G_{C}}-P^{n+1}|_{I_{C}}}{\Delta
x_{1}}$}}{\hbox{\pagecolor{teal!17}$\textstyle\mathstrut\Delta
t\frac{1}{\rho_{G}^{\prime}}\frac{P^{n+1}|_{G_{C}}-P^{n+1}|_{I_{C}}}{\Delta
x_{1}}$}}{\hbox{\pagecolor{teal!17}$\scriptstyle\mathstrut\Delta
t\frac{1}{\rho_{G}^{\prime}}\frac{P^{n+1}|_{G_{C}}-P^{n+1}|_{I_{C}}}{\Delta
x_{1}}$}}{\hbox{\pagecolor{teal!17}$\scriptscriptstyle\mathstrut\Delta
t\frac{1}{\rho_{G}^{\prime}}\frac{P^{n+1}|_{G_{C}}-P^{n+1}|_{I_{C}}}{\Delta
x_{1}}$}}$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to120.69pt{\vbox to10.14pt{\pgfpicture\makeatletter\hbox{\hskip
60.34436pt\lower-2.15944pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-60.34436pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\scriptscriptstyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{teal!17}$\displaystyle\mathstrut\Delta
t\frac{1}{\rho_{G}^{\prime}}\frac{P^{n+1}|_{G_{C}}-P^{n+1}|_{I_{C}}}{\Delta
x_{1}}$}}{\hbox{\pagecolor{teal!17}$\textstyle\mathstrut\Delta
t\frac{1}{\rho_{G}^{\prime}}\frac{P^{n+1}|_{G_{C}}-P^{n+1}|_{I_{C}}}{\Delta
x_{1}}$}}{\hbox{\pagecolor{teal!17}$\scriptstyle\mathstrut\Delta
t\frac{1}{\rho_{G}^{\prime}}\frac{P^{n+1}|_{G_{C}}-P^{n+1}|_{I_{C}}}{\Delta
x_{1}}$}}{\hbox{\pagecolor{teal!17}$\scriptscriptstyle\mathstrut\Delta
t\frac{1}{\rho_{G}^{\prime}}\frac{P^{n+1}|_{G_{C}}-P^{n+1}|_{I_{C}}}{\Delta
x_{1}}$}}$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}\end{split}$ (25)
Solution from PPEPressure correction and surface tension forces

During the solution of Equation 24 and 25, the pressure ($P^{n+1}$) at the
outlet is set equal to $0$ by extrapolating values from $I_{C}$ to $G_{c}$,

$P^{n+1}|_{G_{C}}=-P^{n+1}|_{I_{C}}$ (26)

Consequently, $u^{n+1}_{1}$ at the domain boundary $B$ is computed using the
pressure gradient during the corrector step in Equation 25. The velocity
values at the face-center guard-cells, $G_{F}$, are then extrapolated from
$I_{F}$ to satisfy a zero gradient in the direction normal to the boundary,
$u^{n+1}_{1}|_{G_{F}}=u^{n+1}_{1}|_{I_{F}}$. The solution of PPE corrects the
predicted velocity, $u^{*}_{1}$ in Equation 25 and implicitly satisfies mass
conservation, therefore the value of $u^{n+1}_{1}$ at domain boundary should
remain untouched after the solution. For inlet, slip and no-slip boundary
conditions, the pressure gradient at the boundary is set to zero which also
enforces continuity for $u^{n+1}_{1}$ at the boundary $B$.

The overall algorithm for the fractional-step method is described as follows,

  1. 1.

Solve level-set advection Equation 7, and compute variable properties
described in Equation 12 and 16.

  2. 2.

Compute source terms $\text{S}^{\Gamma}_{u_{i}}$,
$\text{S}^{\Gamma}_{\partial_{i}P}$, and $\text{S}^{\Gamma}_{T}$.

  3. 3.

Solve advection-diffusion Equation 11b for temperature.

  4. 4.

Compute $\dot{m}$ using Equation 6

  5. 5.

Compute outflow forcing terms $\text{S}^{O}_{u_{i}}$,
$\text{S}^{O}_{\partial_{i}P}$ described in Equation 31 and 32.

  6. 6.

Solve Equation 17 for predicted velocity, $u^{*}_{i}$.

  7. 7.

Enforce boundary conditions for $u^{*}_{i}$, including at face-center
boundary, $B$.

  8. 8.

Solve PPE Equation 19 using multigrid methods supplied by AMReX. Boundary
conditions for $P^{n+1}$ are built into the linear system.

  9. 9.

Solve Equation 18 for corrected velocity, $u^{n+1}_{i}$.

  10. 10.

Enforce boundary conditions only $u^{n+1}_{i}$, except at face-center
boundary, $B$.

The forcing terms for outflow $\text{S}^{\text{O}}_{u_{i}}$ and
$\text{S}^{\text{O}}_{\partial_{i}P}$ are applied to internal cells within a
buffer region of length, $l_{b}$, near the outflow boundary. To define the
outflow region we use a weighting expression,

$h_{i}=\Bigg{(}\frac{2}{1+e^{4({x^{max}_{i}-x_{i}})/l_{b}}}\Bigg{)}$ (27)

where, $x$ is the spatial location of the grid point, $x^{max}$ is the upper
bound of the domain, $i=\\{1,2,3\\}$ is the tensor notation for the spatial
directions. Fig 5b shows the profile of this function for different values of
$l_{b}$. The choice of this profile was motivated by the need to concentrate
the effect of the forcing function only near the boundary while also avoid
conditional statements within stenciled computation to enable loop unrolling
and optimization. The exponential function in Equation 27 served as a suitable
candidate for this requirement. Choice for $l_{b}$ is governed by the size of
disturbances that exit out of the computational domain. For multiphase flows
like boiling described in Fig 1 these disturbances are proportional to the
capillary length, $l_{c}=\sqrt{\sigma/(\rho_{L}-\rho_{G})\>g}$, of the liquid-
gas pair and the width of the channel in direction parallel to boundary. As we
demonstrate in Section 5, choosing $l_{b}$ within this range adequately damps
fluctuations that cause instabilities in the flow and preserves statistical
behavior of the solution in the upstream region.

The phase averaged flow within the buffer region is given as the integral over
control volume of the domain, $V$,

$\overline{u}_{L_{i}}=\sum_{j=1}^{3}\delta_{ij}\frac{\int_{V}(1-\text{H})(u_{j}^{n}h_{j}\hat{o}_{j})dV}{\int_{V}(1-\text{H})(h_{j}\hat{o}_{j})dV},\>\>\overline{u}_{G_{i}}=\sum_{j=1}^{3}\delta_{ij}\frac{\int_{V}\text{H}(u_{j}^{n}h_{j}\hat{o}_{j})dV}{\int_{V}\text{H}(h_{j}\hat{o}_{j})dV}$
(28)

where, $\hat{o}_{j}=1$, if the flow is defined in the direction $j$, or $0$
otherwise. $\delta_{ij}$ is the Kronecker delta, and H is the sharp Heaviside
function from Equation 15a,

$\delta_{ij}=\begin{cases}1,&\text{if}\>i=j\\\ 0,&\text{if}\>i\neq j\\\
\end{cases}$ (29)

The effective outflow velocity, $u_{E}$, is given by,

$u_{E_{i}}=\text{max}\big{(}{u_{c}},\overline{u}_{L_{i}}\big{)}(1-\text{H})+\text{max}\big{(}{u_{c}},\overline{u}_{G_{i}}\big{)}(\text{H})$
(30)

setting the lower bound to terminal velocity, ${u_{c}=\sqrt{gl_{c}}}$. The
outflow forcing terms are subsequently given by,

$\begin{split}\text{S}^{\text{O}}_{u_{i}}=\mathchoice{\leavevmode\hbox
to489.41pt{\vbox to30.42pt{\pgfpicture\makeatletter\hbox{\hskip
244.7036pt\lower-11.00012pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-244.7036pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\displaystyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{red!17}$\displaystyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{j}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{red!17}$\textstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{j}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{red!17}$\scriptstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{j}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{red!17}$\scriptscriptstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{j}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to489.41pt{\vbox to30.42pt{\pgfpicture\makeatletter\hbox{\hskip
244.7036pt\lower-11.00012pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-244.7036pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\textstyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{red!17}$\displaystyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{j}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{red!17}$\textstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{j}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{red!17}$\scriptstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{j}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{red!17}$\scriptscriptstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{j}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to370.85pt{\vbox to25.39pt{\pgfpicture\makeatletter\hbox{\hskip
185.4257pt\lower-7.00009pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-185.4257pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\scriptstyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{red!17}$\displaystyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{j}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{red!17}$\textstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{j}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{red!17}$\scriptstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{j}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{red!17}$\scriptscriptstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{j}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to291.81pt{\vbox to22.71pt{\pgfpicture\makeatletter\hbox{\hskip
145.90707pt\lower-5.00006pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-145.90707pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\scriptscriptstyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{red!17}$\displaystyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{j}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{red!17}$\textstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{j}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{red!17}$\scriptstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{j}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{red!17}$\scriptscriptstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{j}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}+\mathchoice{\leavevmode\hbox
to537.71pt{\vbox to54.13pt{\pgfpicture\makeatletter\hbox{\hskip
268.85341pt\lower-34.7085pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-268.85341pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\displaystyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{teal!17}$\displaystyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{k}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{teal!17}$\textstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{k}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{teal!17}$\scriptstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{k}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{teal!17}$\scriptscriptstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{k}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to537.71pt{\vbox to54.13pt{\pgfpicture\makeatletter\hbox{\hskip
268.85341pt\lower-34.7085pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-268.85341pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\textstyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{teal!17}$\displaystyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{k}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{teal!17}$\textstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{k}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{teal!17}$\scriptstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{k}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{teal!17}$\scriptscriptstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{k}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to499.61pt{\vbox to26.39pt{\pgfpicture\makeatletter\hbox{\hskip
249.80545pt\lower-8.00009pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-249.80545pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\scriptstyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{teal!17}$\displaystyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{k}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{teal!17}$\textstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{k}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{teal!17}$\scriptstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{k}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{teal!17}$\scriptscriptstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{k}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to388.86pt{\vbox to22.71pt{\pgfpicture\makeatletter\hbox{\hskip
194.43217pt\lower-5.00006pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-194.43217pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\scriptscriptstyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{teal!17}$\displaystyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{k}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{teal!17}$\textstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{k}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{teal!17}$\scriptstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{k}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}{\hbox{\pagecolor{teal!17}$\scriptscriptstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\bBigg@{4}[\frac{\text{min}\Big{(}1,\big{\lvert}\frac{u_{E_{k}}}{u_{j}^{n}}\big{\rvert}\Big{)}u_{j}^{n}-u_{j}^{n}}{\Delta
t}\bBigg@{4}]$}}$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}\\\
-\mathchoice{\leavevmode\hbox to238.74pt{\vbox
to17.5pt{\pgfpicture\makeatletter\hbox{\hskip
119.37021pt\lower-10.00012pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-119.37021pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\displaystyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{blue!17}$\displaystyle\mathstrut\sum_{j=1}^{3}h_{j}\hat{o}_{j}\lvert
u_{E_{j}}\rvert\partial_{j}u_{i}^{n}$}}{\hbox{\pagecolor{blue!17}$\textstyle\mathstrut\sum_{j=1}^{3}h_{j}\hat{o}_{j}\lvert
u_{E_{j}}\rvert\partial_{j}u_{i}^{n}$}}{\hbox{\pagecolor{blue!17}$\scriptstyle\mathstrut\sum_{j=1}^{3}h_{j}\hat{o}_{j}\lvert
u_{E_{j}}\rvert\partial_{j}u_{i}^{n}$}}{\hbox{\pagecolor{blue!17}$\scriptscriptstyle\mathstrut\sum_{j=1}^{3}h_{j}\hat{o}_{j}\lvert
u_{E_{j}}\rvert\partial_{j}u_{i}^{n}$}}$}} }}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to238.74pt{\vbox to17.5pt{\pgfpicture\makeatletter\hbox{\hskip
119.37021pt\lower-10.00012pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-119.37021pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\textstyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{blue!17}$\displaystyle\mathstrut\sum_{j=1}^{3}h_{j}\hat{o}_{j}\lvert
u_{E_{j}}\rvert\partial_{j}u_{i}^{n}$}}{\hbox{\pagecolor{blue!17}$\textstyle\mathstrut\sum_{j=1}^{3}h_{j}\hat{o}_{j}\lvert
u_{E_{j}}\rvert\partial_{j}u_{i}^{n}$}}{\hbox{\pagecolor{blue!17}$\scriptstyle\mathstrut\sum_{j=1}^{3}h_{j}\hat{o}_{j}\lvert
u_{E_{j}}\rvert\partial_{j}u_{i}^{n}$}}{\hbox{\pagecolor{blue!17}$\scriptscriptstyle\mathstrut\sum_{j=1}^{3}h_{j}\hat{o}_{j}\lvert
u_{E_{j}}\rvert\partial_{j}u_{i}^{n}$}}$}} }}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to180.45pt{\vbox to14.5pt{\pgfpicture\makeatletter\hbox{\hskip
90.22574pt\lower-7.00009pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-90.22574pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\scriptstyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{blue!17}$\displaystyle\mathstrut\sum_{j=1}^{3}h_{j}\hat{o}_{j}\lvert
u_{E_{j}}\rvert\partial_{j}u_{i}^{n}$}}{\hbox{\pagecolor{blue!17}$\textstyle\mathstrut\sum_{j=1}^{3}h_{j}\hat{o}_{j}\lvert
u_{E_{j}}\rvert\partial_{j}u_{i}^{n}$}}{\hbox{\pagecolor{blue!17}$\scriptstyle\mathstrut\sum_{j=1}^{3}h_{j}\hat{o}_{j}\lvert
u_{E_{j}}\rvert\partial_{j}u_{i}^{n}$}}{\hbox{\pagecolor{blue!17}$\scriptscriptstyle\mathstrut\sum_{j=1}^{3}h_{j}\hat{o}_{j}\lvert
u_{E_{j}}\rvert\partial_{j}u_{i}^{n}$}}$}} }}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to141.59pt{\vbox to12.5pt{\pgfpicture\makeatletter\hbox{\hskip
70.79605pt\lower-5.00006pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-70.79605pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\scriptscriptstyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{blue!17}$\displaystyle\mathstrut\sum_{j=1}^{3}h_{j}\hat{o}_{j}\lvert
u_{E_{j}}\rvert\partial_{j}u_{i}^{n}$}}{\hbox{\pagecolor{blue!17}$\textstyle\mathstrut\sum_{j=1}^{3}h_{j}\hat{o}_{j}\lvert
u_{E_{j}}\rvert\partial_{j}u_{i}^{n}$}}{\hbox{\pagecolor{blue!17}$\scriptstyle\mathstrut\sum_{j=1}^{3}h_{j}\hat{o}_{j}\lvert
u_{E_{j}}\rvert\partial_{j}u_{i}^{n}$}}{\hbox{\pagecolor{blue!17}$\scriptscriptstyle\mathstrut\sum_{j=1}^{3}h_{j}\hat{o}_{j}\lvert
u_{E_{j}}\rvert\partial_{j}u_{i}^{n}$}}$}} }}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}\end{split}$ (31)
Scaling in streamwise directionScaling in normal directionconvection

$\text{S}^{\text{O}}_{\partial_{i}P}=-\mathchoice{\leavevmode\hbox
to185.61pt{\vbox to18.75pt{\pgfpicture\makeatletter\hbox{\hskip
92.80408pt\lower-10.00012pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-92.80408pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\displaystyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{red!17}$\displaystyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\text{S}^{\Gamma}_{\partial_{i}P}$}}{\hbox{\pagecolor{red!17}$\textstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\text{S}^{\Gamma}_{\partial_{i}P}$}}{\hbox{\pagecolor{red!17}$\scriptstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\text{S}^{\Gamma}_{\partial_{i}P}$}}{\hbox{\pagecolor{red!17}$\scriptscriptstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\text{S}^{\Gamma}_{\partial_{i}P}$}}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to185.61pt{\vbox to18.75pt{\pgfpicture\makeatletter\hbox{\hskip
92.80408pt\lower-10.00012pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-92.80408pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\textstyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{red!17}$\displaystyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\text{S}^{\Gamma}_{\partial_{i}P}$}}{\hbox{\pagecolor{red!17}$\textstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\text{S}^{\Gamma}_{\partial_{i}P}$}}{\hbox{\pagecolor{red!17}$\scriptstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\text{S}^{\Gamma}_{\partial_{i}P}$}}{\hbox{\pagecolor{red!17}$\scriptscriptstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\text{S}^{\Gamma}_{\partial_{i}P}$}}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to136.59pt{\vbox to14.22pt{\pgfpicture\makeatletter\hbox{\hskip
68.29614pt\lower-7.00009pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-68.29614pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\scriptstyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{red!17}$\displaystyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\text{S}^{\Gamma}_{\partial_{i}P}$}}{\hbox{\pagecolor{red!17}$\textstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\text{S}^{\Gamma}_{\partial_{i}P}$}}{\hbox{\pagecolor{red!17}$\scriptstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\text{S}^{\Gamma}_{\partial_{i}P}$}}{\hbox{\pagecolor{red!17}$\scriptscriptstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\text{S}^{\Gamma}_{\partial_{i}P}$}}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to103.91pt{\vbox to12.22pt{\pgfpicture\makeatletter\hbox{\hskip
51.95746pt\lower-5.00006pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-51.95746pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\scriptscriptstyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{red!17}$\displaystyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\text{S}^{\Gamma}_{\partial_{i}P}$}}{\hbox{\pagecolor{red!17}$\textstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\text{S}^{\Gamma}_{\partial_{i}P}$}}{\hbox{\pagecolor{red!17}$\scriptstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\text{S}^{\Gamma}_{\partial_{i}P}$}}{\hbox{\pagecolor{red!17}$\scriptscriptstyle\mathstrut\sum_{j=1}^{3}\delta_{ij}h_{j}\hat{o}_{j}\text{S}^{\Gamma}_{\partial_{i}P}$}}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}-\mathchoice{\leavevmode\hbox
to360.87pt{\vbox to18.75pt{\pgfpicture\makeatletter\hbox{\hskip
180.43434pt\lower-10.00012pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-180.43434pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\displaystyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{teal!17}$\displaystyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\text{S}^{\Gamma}_{\partial_{j}P}$}}{\hbox{\pagecolor{teal!17}$\textstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\text{S}^{\Gamma}_{\partial_{j}P}$}}{\hbox{\pagecolor{teal!17}$\scriptstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\text{S}^{\Gamma}_{\partial_{j}P}$}}{\hbox{\pagecolor{teal!17}$\scriptscriptstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\text{S}^{\Gamma}_{\partial_{j}P}$}}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to360.87pt{\vbox to18.75pt{\pgfpicture\makeatletter\hbox{\hskip
180.43434pt\lower-10.00012pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-180.43434pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\textstyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{teal!17}$\displaystyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\text{S}^{\Gamma}_{\partial_{j}P}$}}{\hbox{\pagecolor{teal!17}$\textstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\text{S}^{\Gamma}_{\partial_{j}P}$}}{\hbox{\pagecolor{teal!17}$\scriptstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\text{S}^{\Gamma}_{\partial_{j}P}$}}{\hbox{\pagecolor{teal!17}$\scriptscriptstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\text{S}^{\Gamma}_{\partial_{j}P}$}}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to264.61pt{\vbox to14.22pt{\pgfpicture\makeatletter\hbox{\hskip
132.30386pt\lower-7.00009pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-132.30386pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\scriptstyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{teal!17}$\displaystyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\text{S}^{\Gamma}_{\partial_{j}P}$}}{\hbox{\pagecolor{teal!17}$\textstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\text{S}^{\Gamma}_{\partial_{j}P}$}}{\hbox{\pagecolor{teal!17}$\scriptstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\text{S}^{\Gamma}_{\partial_{j}P}$}}{\hbox{\pagecolor{teal!17}$\scriptscriptstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\text{S}^{\Gamma}_{\partial_{j}P}$}}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to200.43pt{\vbox to12.22pt{\pgfpicture\makeatletter\hbox{\hskip
100.21684pt\lower-5.00006pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-100.21684pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\scriptscriptstyle\definecolor{currentcolor}{rgb}{0,0,0}\mathchoice{\hbox{\pagecolor{teal!17}$\displaystyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\text{S}^{\Gamma}_{\partial_{j}P}$}}{\hbox{\pagecolor{teal!17}$\textstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\text{S}^{\Gamma}_{\partial_{j}P}$}}{\hbox{\pagecolor{teal!17}$\scriptstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\text{S}^{\Gamma}_{\partial_{j}P}$}}{\hbox{\pagecolor{teal!17}$\scriptscriptstyle\mathstrut\sum_{j=1}^{3}\sum_{k=1}^{3}\delta_{ij}(1-\delta_{jk})h_{k}\hat{o}_{k}\text{S}^{\Gamma}_{\partial_{j}P}$}}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{}{{
{}{}{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}$ (32) Forcing in
streamwise directionForcing in normal direction

Figure 5: (a) Staggered computational grid near the outflow with faces
($\times$) and cell-centers ($\bullet$), $I$ denotes the internal cells, $G$
denotes the guard cells, and $B$ is the domain boundary. (b) Shape of outflow
profiles, $h_{2}$, for different buffer region lengths, $l_{b}$.

Each term in Equations 31 and 32 is highlighted to show its significance. The
buffer region represents a far-field with a mean velocity, $u_{e}$. Local
velocities with magnitudes greater than $u_{e}$ are scaled in both streamwise
and normal direction. It is important to note that scaling does not alter the
direction of the velocity and only ensures that magnitude of local velocities
cannot be greater than $u_{e}$. Natural backflow can still exist with $P=0$
outflow boundary condition. Enforcing the condition $u_{e}\geq u_{c}$ ensures,
at the very least, that capillary flow conditions are satisfied within the
buffer region. The convection term adds another forcing to allow turbulent
vortices to exit the domain. $\text{S}^{\text{O}}_{\partial_{i}P}$ counteracts
pressure jumps, $\text{S}^{\Gamma}_{\partial_{i}P}$, and completely nullifies
its contribution at the outflow, which aligns with the far-field assumption
within the buffer region. In the following section we present results from
some numerical experiments that demonstrate efficacy of this outflow forcing.

##  5 Results

This section presents a set of verification problems to demonstrate the
stability of outflow forcing described in the section above. Problems are
presented with increasing complexity, starting from benchmark multiphase flow
problems to homogeneous evaporation and heterogeneous pool/flow boiling. We
use laboratory notebooks that organizes each computational experiment using a
set of configuration files to enable data curation and reproducibility[10].
The laboratory notebook and Flash-X source code are open source to allow
community development and contribution [11, 7, 6].

###  5.1 Comparison with multiphase flow benchmarks

We start with a benchmark problem of a two-dimensional rising bubble without
phase transition. We chose this problem because of availability of data from
multiple solvers [17]. Figure 6(a) describes the computational setup used in
developing these benchmarks, and Figure 6(b) describes the setup that replaces
the no-slip boundary condition at the top with a pressure outflow and buffer
regions to allow for testing our formulation. We perform simulations using
both configurations to allow for comparing results from different Flash-X
implementations and other solvers.

Figure 6: Schematic of the benchmark rising bubble problem (a) With a no-slip
boundary condition at the top of the domain, and (b) With outflow boundary
condition along with buffer regions for testing

Flow configurations provided in [17] allowed us to calculate non-dimensional
parameters described in Table 1, along with the reference length scale
($l_{0}=0.5m$), and time scale ($t_{0}=0.7s$). $l_{0}$ corresponds to the
diameter of the bubble. Simulations are initialized by placing the bubble with
$D=1$ at the center $(0,0)$ of the computational domain, which is also scaled
using $l_{0}$.

Parameter | $\rho_{G}/\rho_{L}$ | $\mu_{G}/\mu_{L}$ | Re | Fr | We
---|---|---|---|---|---
Values | $0.001$ | $0.01$ | $35$ | $1$ | $125$
Table 1: Parameters for the rising bubble benchmark problem

Starting with a resolution of $D/80$, we simulate all the configurations
described in Figure 6. To compare simulation results, we track evolution of
circularity ($C$), bubble center of mass ($X_{2}$), and bubble rise velocity
($U_{2}$) with time. Equations for these quantities are described in [17] and
elaborated below,

  1. 1.

$C$: Ratio of area-equivalent circle to perimeter of the bubble. For a perfect
circle, $C=1$, which decreases as the bubble deforms.

  2. 2.

$X_{2}$: Normal component of the bubble center of mass, computed using
integral of area occupied by the bubble.

  3. 3.

$U_{2}$: Normal component of the mean velocity of the bubble.

Figure 7: Comparison of benchmark quantities for different simulation
configurations. The no-slip case refers to the setup reported in [17], which
is compared with the setup that uses the outflow forcing.

From the results we can see that the outflow forcing closes matches the no-
slip solution for $l_{b}=\\{0.5,1\\}$, and starts to deviate when $l_{b}>D$,
the size of the disturbance in this case. The deviation is particularly
prominent in the evolution of rise velocity, $U_{2}$ (Figure 7c), which is a
direct consequence of discrepancy in the bubble break-up pattern at $t=3.5s$
(Figure 7d). Sensitivity of bubble break-up event for this benchmark problem
has been reported in previous studies [17] as well, with different numerical
formulations exhibiting different behaviors.

Figure 8: Comparison of benchmark quantities for different spatial resolutions. Coarser grids, $D/20$ and $D/40$ demonstrate noisy behavior while $D/80$ and $D/160$ show good convergence, although minor discrepancy can still be observed in the bubble break-up behavior. Due to the cost associated with simulating this problem we could perform a study for $D/320$. $\Delta x_{i}$ | $\Delta t$ | $||e_{1}||$ | $\text{ROC}_{1}$ | $||e_{2}||$ | $\text{ROC}_{2}$ | $||e_{\infty}||$ | $\text{ROC}_{\infty}$
---|---|---|---|---|---|---|---
D/20 | $5.0E-05$ | $2.9E-02$ | - | $3.8E-02$ |  | $6.4E-02$ | -
D/40 | $1.25E-05$ | $1.4E-02$ | 1.1 | $1.9E-02$ | 1.0 | $5.8E-02$ | -
D/80 | $3.125E-06$ | $5.8E-03$ | 1.3 | $7.2E-03$ | 1.4 | $1.3E-02$ | 2.2
Table 2: Relative norms and convergence rate of bubble circularity, $C$.  $\Delta x_{i}$ | $\Delta t$ | $||e_{1}||$ | $\text{ROC}_{1}$ | $||e_{2}||$ | $\text{ROC}_{2}$ | $||e_{\infty}||$ | $\text{ROC}_{\infty}$
---|---|---|---|---|---|---|---
D/20 | $5.0E-05$ | $9.8E-03$ | - | $1.9E-02$ |  | $3.5E-02$ | -
D/40 | $1.25E-05$ | $5.7E-03$ | 0.8 | $1.0E-02$ | 0.9 | $1.9E-02$ | 0.88
D/80 | $3.125E-06$ | $2.1E-03$ | 1.44 | $3.8E-03$ | 1.39 | $8.5E-03$ | 1.16
Table 3: Relative norms and convergence rate of bubble center of mass, $X_{2}$.  $\Delta x_{i}$ | $\Delta t$ | $||e_{1}||$ | $\text{ROC}_{1}$ | $||e_{2}||$ | $\text{ROC}_{2}$ | $||e_{\infty}||$ | $\text{ROC}_{\infty}$
---|---|---|---|---|---|---|---
D/20 | $5.0E-05$ | $6.5E-02$ | - | $8.8E-02$ | - | $2.5E-01$ | -
D/40 | $1.25E-05$ | $4.0E-02$ | 0.7 | $6.3E-02$ | - | $1.6E-01$ | 0.64
D/80 | $3.125E-06$ | $1.5E-02$ | 1.4 | $2.5E-02$ | 1.33 | $6.0E-02$ | 1.4
Table 4: Relative norms and convergence rate of bubble velocity, $U_{2}$.

Next we perform a grid refinement study using the outflow forcing
configuration with $l_{b}=1$ to estimate the convergence rate of our numerical
formulation. The results for this study are reported in Figure 8. Coarser
grids, $D/20$ and $D/40$, demonstrate a lot of noise in the evolution of
circularity ($C$) and rise velocity ($U_{2}$). Even though good convergence is
achieved between $D/80$ and $D/160$, the solutions are far from grid
independent, which is clearly evident in the evolution of the rise velocity
($U_{2}$) which is again a result of discrepancy in the bubble break-up
behavior. We could not perform a simulation at resolution of $D/320$ due to
the time-stepping constraints imposed by the explicit nature of our
formulation.

Tables 2, 3, and 4, provide error norms and convergence rates for the spatial
resolution study in Figure 7. The norms are computed using $D/160$ as the
reference solution. Convergence rates are computed using,

$ROC=\frac{log(||e||_{coarse}/||e||_{fine})}{log(\Delta x_{coarse}/\Delta
x_{fine})}$ (33)

Due to the noisy behavior on coarser grids, the convergence rates between
$D/20$ and $D/40$ the convergence rates between them are less than linear.
However rates between $D/40$ and $D/80$ are greater than linear nearing a
second-order convergence.

Finally we compare our simulations results from the finest resolution to the
studies reported previously [17] (Figure 9. Perfect agreement between results
is not obtained due to the sensitivity of the study. However, general trend in
circularity, $C$, and velocity ($U_{2}$) is captured. Notice the timing of two
peaks for $U_{2}$. The first is peak observed when bubble reaches terminal
velocity, and the second is observed near the breakup.

Figure 9: Comparison of Flash-X simulation results for rising bubble benchmark
with previously reported studies in [17]. Due to the sensitivity of the
problem the agreement is not perfect between any of the dataset, however the
general trend in circularity ($C$) and velocity ($U_{2}$) is captured.

###  5.2 Bubble evaporation in superheated liquid

Two dimensional homogeneous evaporation of a gas (vapor) bubble in superheated
liquid is another case to compare the effect of the outflow forcing with a
solution where no forcing is applied
($\text{S}^{\text{O}}_{u_{i}}=\text{S}^{\text{O}}_{\partial_{i}P}=0$). The
schematic of the problem is provided in Figure 10 with dimensions of the
domain in directions $x_{1}$ and $x_{2}$. The simulation is initialized by
placing a spherical bubble of radius $r=0.25$ and center at the coordinates
$(0,0)$. The temperature inside the bubble is set to $T_{sat}=0$, and the
temperature inside the liquid is set to $T_{sup}=1$. Table 5 provides values
for the remaining non-dimensional parameters for the simulation. A resolution
of $D/16$ along with a time step of $2.5E-05$ is used for this study.

Figure 10: Schematic of the evaporating bubble problem with corresponding outflow profiles, $h_{2}$, for different buffer region lengths, $l_{b}$. Parameter | $\rho_{G}/\rho_{L}$ | $\mu_{G}/\mu_{L}$ | $C_{p_{G}}/C_{p_{L}}$ | $k_{G}/k_{L}$ | $\alpha_{G}/\alpha_{L}$ | Re | Pr | St | Fr | We
---|---|---|---|---|---|---|---|---|---|---
Values | $0.0012$ | $0.00025$ | $0.7$ | $0.2$ | $238$ | $200$ | $7$ | $0.004$ | $1$ | $50$
Table 5: Parameters for evaporating bubble case

No-slip boundary conditions are applied everywhere, except for the outflow
where $P=0$. No slip boundary conditions set the velocities and normal
derivative of pressure at the boundaries to zero. The effect of the outflow
forcing is tested for different buffer lengths, $l_{b}=\\{0.5,1.0,1.5\\}$, and
compared with the no forcing solution.

Figure 11: Comparison of temperature (top), and evolution of liquid-gas
interface, $\phi=0$ (bottom), between no forcing solution and solution
obtained using outflow forcing for varying $l_{b}$. Figure 12:  Evolution of
evaporative mass flux, $\dot{m}$ (left), and divergence of velocity,
$\sum_{j=1}^{3}\partial_{j}u_{j}+\sum_{j=1}^{3}\dot{m}\>n_{\Gamma_{j}}\partial_{j}\frac{1}{\rho^{\prime}}\>\bigr{|}_{\>\Gamma}$
(right) demonstrating the divergence preserving nature of the outflow forcing
($l_{b}=1.0$).

Figure 11 provides the temporal evolution for different cases, which show that
the effects of the outflow forcing only appear in the buffer region and
maintain the accuracy of the flow in the downstream areas. The temperature
contours and bubble shape are consistent for $x^{max}_{i}-x_{i}>l_{b}$. For
$x^{max}_{i}-x_{i}\leq l_{b}$, the bubble undergoes acceleration due to the
presence of the convective term in Equation 31, which changes the shape of the
liquid-gas interface and the corresponding temperature distribution. This test
problem provides good verification of the sharp nature of the outflow forcing.

Divergence preserving nature of our formulation is evident in Figure 12 where
temporal evolution of divergence is of order ${10E-13}$. The figure also show
evolution of evaporative mass flux, $\dot{m}$, which exists only near the
liquid-gas interface.

###  5.3 Two dimensional pool boiling simulation with a single nucleation site

The next problem we consider is the simplified version of the pool boiling
problem described in Figure 1b, where we look at dynamics associated with a
single nucleation site. This simulation not only adds an additional layer
complexity due to the heterogeneous nature of the boiling process (boiling
over a heater), but also provides an opportunity to understand the nature of
vortex-induced instability at the outflow by considering effects associated
with a single bubble. The size of the computational domain is shown in Figure
13, with the heater completely occupying the lower boundary in the direction
$x_{2}$.

Figure 13: (a) Comparison of pressure and velocity during the intial cycles of
bubble growth and departure for simulation without outflow forcing (top) and
simulation with $l_{b}=1.0$ (bottom). (b) Evolution of divergence with time.

We follow the nucleation site model discussed in our previous work [8, 9, 5],
which provides site location and nucleation wait time, $t_{wait}$, as input to
our simulation. The nucleation site model seeds a bubble of radius $0.1$ at a
given nucleation site if it has been vacant (covered by liquid) for
$t=t_{wait}$. For the test problem with a single nucleation site, the bubble
is seeded at location $(0,0)$.

The dynamics of the contact line is modeled by satisfying a contact angle,
$\psi$, for the level-set function, $\partial_{2}\phi=-cos(\psi)$, at the
heater boundary. $\psi$ is calculated using,

$\psi=\begin{cases}\psi_{r},&\text{if }u_{base}<0\\\
\frac{\psi_{a}-\psi_{r}}{u_{lim}}u_{base}+\psi_{r},&\text{if }0\leq
u_{base}\leq u_{lim}\\\ \psi_{a},&\text{if }u_{base}>u_{lim}\\\ \end{cases}$
(34)

where, $\psi_{r}$ is the receding and $\psi_{a}$ is the advancing contact
angle, $u_{base}$ is the near wall radial velocity around the interface in the
plane parallel to heater and $u_{lim}$ is the limiting velocity. This model is
a variation of that proposed by Mukherjee et al. [23]. They assume inertial
effects to be relevant both in the advancing and receding conditions, whereas
our model accounts for inertial effects only in the advancing state. From our
numerical experiments and based on results from Mukherjee et al. we found that
an optimum value for $u_{lim}$ is $20-25\%$ of the characteristic velocity,
$u_{0}=\sqrt{g\>l_{0}}$, where $g$ is acceleration due to gravity, and $l_{0}$
is the characteristic length scale given by
$l_{0}=\sqrt{\sigma/(\rho_{l}-\rho_{v})\>g}$. Table 6 provides the values of
the remaining input parameters for this simulation.

Parameter | $\rho_{G}/\rho_{L}$ | $\mu_{G}/\mu_{L}$ | $C_{p_{G}}/C_{p_{L}}$ | $k_{G}/k_{L}$ | $\alpha_{G}/\alpha_{L}$ | Re | Pr | St | Fr | We | $\psi_{r}$ | $\psi_{a}$ | $t_{wait}$ | $T_{sat}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Values | $0.0049$ | $0.02$ | $0.735$ | $0.143$ | $39$ | $303$ | $7$ | $0.16$ | $1$ | $1$ | $45\degree$ | $45\degree$ | 0.4 | 0
Table 6: Parameters for pool boiling with single nucleation site

Five different sets of studies were performed to compare the reference
solution ($L_{x_{2}}=24$), no forcing solution ($L_{x_{2}}=6$), and three
outflow forcing solution ($L_{x_{2}}=6$) with different values of $l_{b}$. All
simulations were performed using a fixed time step, $\Delta t=10^{-4}$, and
resolution, $\Delta x_{i}=0.02$. The comparison in Figure 13a for the initial
cycles of bubble growth and departure shows the instabilities created at the
exit when forcing is not applied. Each cycle consists of a single bubble
nucleating, growing, and departing from the heater surface. The disturbances
are the result of sharp interfacial jumps in pressure associated with
$\text{We}=1$. The simulation between the two cases starts to deviate at
$t=53$, when the disturbances in the reference solution travel downstream and
create unrealistic flow patterns that deform the nucleation of the bubble in
the heater for the third cycle of bubble growth and departure. Evolution of
divergence in Figure 15b, shows effect of this instability on growth of
divergence in the domain which eventually causes the simulation to blow-up. By
setting the outflow forcing within, $l_{b}=1.0$, quasi-steady state can be
reached.

Figure 14: (a) Quasi-steady behavior of bubble growth and departure cycle for
different values of $l_{b}$, (b) Temporal variation in the heat flux
$\tilde{Q}$, and (c) Growth of equivalent bubble diameter $D_{g}$.

The temporal and quasi-steady behavior shown in Figure14 compares results
obtained from different values of $l_{b}$ against each other and the reference
solution ($L_{x_{2}}=24$). See the progression of the bubble growth and
departure cycle from $t_{c}=0-12$ in Figure14a. A more quantified comparison
is found in Figures 14b and c, which compares the behavior of the mean wall
heat flux, $\tilde{Q}$, and the equivalent diameter during bubble growth,
$D_{g}$ as a function of time. $\tilde{Q}$ is an integral value computed
using,

$\tilde{Q}=\frac{\int_{S_{heater}}I_{L}\partial_{2}TdS}{\int_{S_{heater}}dS}$
(35)

where, $S_{heater}$, is the surface area of the heater, $I_{L}$ is the
indicator function for the liquid, and $\partial_{2}T$ is the temperature
gradient in the cell closest to the heater wall.

The equivalent diameter, $D_{g}$, is computed using,

$D_{g}=2\sqrt{\frac{A_{bubble}}{\pi}}$ (36)

where, $A_{bubble}$, is the bubble area during the growth period. The temporal
variation for both $\tilde{Q}$ and $D_{g}$ is consistent between the two runs
with negligible differences.

It is evident that the outflow forcing introduces a phase difference in
computation of these quantities, which is an expected effect of this
formulation. However, the temporal trend is statistically consistent which is
an important result. The value for $\overline{Q}$,

$\overline{Q}=\frac{\sum_{n=1}^{Nt}\tilde{Q}}{Nt}$ (37)

for different outflow forcing shows a variation of $3\%$ from the reference
solution. However, there is certainly scope to improve upon this forcing and
design a more robust method using traction boundary conditions.

Figure 15: Effect of convection term in Equation 31 term on the solution for
pool boiling problem. Figure 16:  Spatio-temporal statistics for the single
bubble pool boiling problem showing differences in mean and fluctuations for
velocity normal to heater surface ($u_{2}$) and temperature ($T$), along with
turbulent heat flux ($\overline{\overline{u_{2}^{\prime}T^{\prime}}}$).

The effect of convection term in Equation 31 is described in Figure 15 for the
initial cycle of bubble growth and departure, which demonstrates that the
convection term improves the quality of the outflow forcing in comparison to
the solution where the contribution of this term is zero.

Figure 16 provides a description of spatio-temporal statistics for this
problem which gives an overview of the upstream effects of the proposed
outflow forcing. The statistics represent the mean and fluctuations which are
described using following equation for a physical variable $f(x_{1},x_{2},t)$,

$\overline{f}=\frac{\sum_{n=1}^{Nt}f}{Nt}\quad\quad\quad\overline{\overline{f}}=\frac{\sum_{n=1}^{N_{x_{1}}}\sum_{n=1}^{Nt}f}{Nt\>{N_{x_{1}}}}$
(38)

Here $\overline{f}$ is mean over time and $\overline{\overline{f}}$ represents
mean over time and spatial direction $x_{1}$. Trends for velocity
($\overline{\overline{u_{2}}}$) and temperature ($\overline{\overline{T}}$)
show very negligible differences when compared to the reference solution from
the longer domain. Additionally, the turbulent heat flux
($\overline{\overline{u_{2}^{\prime}T^{\prime}}}$), which is an important
metric also collapses on the same trend line. The significance of turbulent
heat flux is discussed at length in our previous work [9]. The velocity and
temperature fluctuations that make up the turbulent heat flux are given by,

$f^{\prime}=f-\overline{f}\quad\quad\quad
f^{\prime\prime}=\sqrt{\frac{\sum_{n=1}^{N_{x_{1}}}\sum_{n=1}^{Nt}(f-\overline{f})^{2}}{Nt\>{N_{x_{1}}}}}$
(39)

The trends for $u_{2}^{\prime\prime}$ and $T^{\prime\prime}$ which represent
the fluctuations in velocity and temperature over space and time show the most
discrepancy. The velocity fluctuations ($u_{2}^{\prime\prime}$) quantify the
effect of the outflow forcing on the rising bubble vortices which causes
significant deviation from the reference solution. This trend is expected due
to the damping nature of the forcing. The upstream solution near the heater
wall remains unaffected and the overall trend for turbulent heat flux make
this an effective strategy for treating outflows.

###  5.4 Two dimensional pool boiling simulations with multiple nucleation
sites

For this study, we performed a gravity scaling analysis and compared it
against the heat flux model proposed by Raj et al. [27, 26]. Multiple
nucleation sites are calculated using a Halton sequence based on the approach
in our previous work [8, 9]. A total of $30$ nucleation sites are seeded using
the nucleation site model discussed above. The site distribution is obtained
from our previous work and is scaled from 3D to 2D for the present study. The
study parameters are provided in Table 7 and the resolution, $\Delta
x_{i}=0.03$.

Figure 17: Vorticity contours during quasi-steady for different gravity levels, $t_{1}<t_{2}<t_{3}<t_{4}<t_{5}$ (top), comparison of variation in $\tilde{Q}$ on earth gravity for two different simulations using only artificial condensation versus outflow forcing described in this paper (bottom-left), and the gravity scaling of quasi-steady state heat flux, $\overline{Q}$ (bottom-right). Parameter | $\rho_{G}/\rho_{L}$ | $\mu_{G}/\mu_{L}$ | $C_{p_{G}}/C_{p_{L}}$ | $k_{G}/k_{L}$ | $\alpha_{G}/\alpha_{L}$ | Re | Pr | St | We | $\psi_{r}$ | $\psi_{a}$ | $t_{wait}$ | $T_{sat}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---
Values | $0.0083$ | $1$ | $0.83$ | $0.25$ | $36$ | $238$ | $8.4$ | $0.5298$ | $1$ | $45\degree$ | $90\degree$ | 0.4 | 0.15
Table 7: Parameters for pool boiling with multiple nucleation sites

Figure 17 provides details of the computational domain described in the
schematic in Figure 1 with the heater located between $x_{1}=-5$ and
$x_{1}=5$. To demonstrate the effectiveness of the proposed formulation we
show the evolution of $\tilde{Q}$ for two different simulations using only
artificial condensation versus the outflow forcing described in Section 4.
Simulation using artificial condensation near the outflow deviates over time
from the quasi-steady behavior, resulting in an error in the heat flux
calculation. This error is caused by non-physical gusts that appear in the
numerical solution due to limitations of the artificial condensation
formulation. Figure 18 provides a schematic of this gust event from time time
$t_{1}$ to $t_{4}$,

Figure 18: Schematic of non-physical gusts introduced by the artificial
condensation method when applied to production problem of pool boiling
involving multiple bubbles ($t_{1}<t_{2}<t_{3}<t_{4}$). This numerical
artifact leads to errors in measurement of wall heat flux.

To compare results with previously reported trends, we performed a range of
simulations for different gravity values, $Fr=1-100$. The scaling,
$a/g=\frac{1}{Fr^{2}}$, scales gravity values with respect to the value on
Earth gravity, $g$. Gravity separates pool boiling in two separate regimes
dominated by buoyancy or surface tension. The regimes are termed buoyancy-
dominated boiling (BDB) and surface tension-dominated boiling (SDB), which are
separated by the transitional acceleration, $a_{trans}$, given by

$a_{trans}=\frac{4.41}{Bo\>L^{2}_{h}}$ (40)

where, $L_{h}=10$, is the length of the heater and $Bo=1$, is the Bond number
for the simulation. Equation 40 was identified by Raj et al. [26, 27] from the
experimental data they collected in varying gravity environments. The choice
of $l_{b}$ is governed by capillary length associated with $a_{trans}$ which
determines the size of disturbances that will exit of the computational
domain.

In BDB, the decrease in gravitational acceleration increases the size of
departing bubbles and reduces the strength of turbulent vortical structures
associated with movement and condensation of smaller bubbles within the
domain, see comparison between vorticity magnitudes for $a/g=1$ and $a/g=0.1$
in Figure 17. The scaling of quasi-steady heat flux, $\overline{Q}$ (Equation
37) with respect to gravity, $\overline{Q}_{a}/\overline{Q}_{g}$, follows the
slope $m_{BDB}$.

In SDB, the dynamics are dominated by the presence of a central bubble that
does not depart the heater surface, but rather acts as a vapor sink for
smaller satellite bubbles, as shown in Figure 17 for $a/g=0.001$. Heat flux
drops sharply by a value $K_{jump}$, which is dependent on the size of the
central bubble [9], and the slope, $m_{SDB}=0$.

The gravity scaling of heat flux computed from the simulations accurately
matches the expected trend of the model proposed by Raj et al. [26, 27], and
serves as validation of the outflow forcing.

###  5.5 Performance evaluation for three-dimensional flow boiling simulations

The final case we present is a demonstration of the capability of our
formulation to deal with flow boiling problems that involve transport of bulk
liquid through the domain, and provide detailed performance analysis of
various components of the formulation. The schematic is shown in Figure 1 a,
and the corresponding input values are shown in Table 8, with resolution,
$\Delta x_{i}=0.03$.

Parameter | $\rho_{G}/\rho_{L}$ | $\mu_{G}/\mu_{L}$ | $C_{p_{G}}/C_{p_{L}}$ | $k_{G}/k_{L}$ | $\alpha_{G}/\alpha_{L}$ | Re | Pr | St | Fr | We | $\psi_{r}$ | $\psi_{a}$ | $t_{wait}$ | $T_{sat}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Values | $0.0083$ | $1$ | $0.83$ | $0.25$ | $36$ | $238$ | $8.4$ | $0.5298$ | $1$ | $1$ | $45\degree$ | $90\degree$ | 0.4 | 0.15
Table 8: Parameters for flow boiling

The simulations are run with the bulk liquid velocity $U_{b}=u_{0}$. Figure 19
shows the dynamics of the flow in earth gravity that is highly bubbly and
turbulent. We performed a weak scaling study for this problem by increasing
the size of the computational domain in proportion to the number of processes
(ranks), to maintain consistent loading of 10 AMR blocks per rank. The
experiments were performed on Summit[1], which is a leadership computing
facility located at Oak Ridge National Laboratory. The reproducibility capsule
and data for this problem is hosted at [6].

Figure 19: Three dimensional simulation of flow boiling on earth gravity,
$a/g=1$, showing bubble dynamics (top), temperature distribution (middle), and
weak scaling performance of different components of the solver (bottom) [6].

Each simulation was run for approximately $7000s$ of wall time and the
performance was logged using native implementation of timers in Flash-X.
Figure 19 shows cost per iteration for each run. As is the case with many
incompressible Navier-Stokes simulations, the solution of Pressure Poisson
Equation (PPE) is the highest contributor to the cost, followed by
communication that is triggered by halo exchanges during guard-cell filling
events. Computation of WENO advection terms for Equation 11 and 7 comes in
third as a highest contributor on lower ranks. WENO advection stencils
typically involve conditional statements that interfere with loop unrolling
and optimization which contributes towards its poor computational performance.
On higher ranks, Input/Output (IO) events dominate over computation of WENO
advection. This is attributed to the increase in the size of domain which in-
turn increases the size of data that is written to the parallel file system.

Outflow forcing formulation proposed in this work shows good performance in
comparison to other components of the algorithm, and costs slighlty less than
computation of interfacial jumps. Overall, the multiphase solver in Flash-X
can compute approximately $10ms$ of physical time in $3600s$ of wall time for
three-dimensional flow boiling problems. This is a significant improvement
over its predecessor FLASH [8]. Table 9 provides the metrics for each run.

Ranks | SimTime (Non-Dimensional) | SimTime ($ms$) | Wall Time ($s$)
---|---|---|---
210 | 1.6 | 16 | 3600
504 | 1.2 | 12 | 3600
1050 | 1.1 | 11 | 3600
2100 | 1.0 | 10 | 3600
Table 9: Weak scaling of physical time per wall time for three-dimensional
flow boiling simulation

##  6 Conclusion

We present an outflow forcing treatment for multiphase flow simulations
involving phase transition and sharp interfacial jumps. The forcing term is
designed for a fractional step predictor-corrector formulation of
incompressible Navier-Stokes equations, that is coupled with a fixed pressure
outflow boundary condition to enable stable exit of multiphase disturbances
from the domain boundary. The forcing is applied within a buffer region near
the outflow to velocity during the predictor step and to pressure jumps during
the solution of pressure Poisson equation.

The concept of an effective outflow velocity is introduced which is used as a
parameter for damping and convection of bubble induced vortices, and is
calculated based on the characteristic flow velocity and mean velocity for
each phase within the outflow region. This strategy can be easily extended to
design outflow treatment for different numerical implementations of multiphase
flow solution.

We demonstrate the applicability and accuracy of this formulation by providing
results from various homogeneous and heterogeneous evaporation problems with
low Weber numbers that result in large instabilities at the outflow if not
handled properly. We conduct a systematic study on sensitivity of the
numerical solution based on length of the buffer region, and show that the
outflow forcing maintains numerical stability whilst preserving statistical
behavior of the solution. We also perform a gravity scaling study for
subcooled pool boiling and compare the results with heat flux models developed
from experimental data to show the ability of our solver to predict different
bubble regimes and their transition accurately. Finally, we show how the new
outflow forcing performs during simulation of high-fidelity flow boiling
problems.

Since the proposed solution introduces damping terms in the solution of
multiphase system they do not perform very well on severely truncated domains.
This presents an opportunity to design more robust traction boundary
conditions for these class of simulations. We hope to build on existing
methods and address these challenges in our future work.

Acknowledgements: The author would like to thank Dr. Anshu Dubey (Argonne
National Laboratory), Dr. Klaus Weide (University of Chicago), and the
reviewers for providing constructive feedback while developing the manuscript.

The material is based upon work supported by Laboratory Directed Research and
Development (LDRD) funding from Argonne National Laboratory, provided by the
Director, Office of Science, of the U.S. Department of Energy under contract
DE-AC02-06CH11357 and the Exascale Computing Project (17-SC-20-SC), a
collaborative effort of the US Department of Energy Office of Science and the
National Nuclear Security Administration.

We also acknowledge that The City of Chicago is located on land that is and
has long been a center for Native peoples. The area is the traditional
homelands of the Anishinaabe, or the Council of the Three Fires: the Ojibwe,
Odawa, and Potawatomi Nations. Many other Nations consider this area their
traditional homeland, including the Myaamia, Ho-Chunk, Menominee, Sac and Fox,
Peoria, Kaskaskia, Wea, Kickapoo, and Mascouten.

The submitted manuscript has been created by UChicago Argonne, LLC, operator
of Argonne National Laboratory (“Argonne”). Argonne, a U.S. Department of
Energy Office of Science laboratory, is operated under Contract No. DE-
AC02-06CH11357. The U.S. Government retains for itself, and others acting on
its behalf, a paid-up nonexclusive, irrevocable worldwide license in said
article to reproduce, prepare derivative works, distribute copies to the
public, and perform publicly and display publicly, by or on behalf of the
Government. The Department of Energy will provide public access to these
results of federally sponsored research in accordance with the DOE Public
Access Plan. http://energy.gov/downloads/doe-public-access-plan.

## References

  * [1] SUMMIT, Oak Ridge National Laboratory.  https://www.olcf.ornl.gov/summit/.
  * [2] Blayo, E., and Debreu, L. Revisiting open boundary conditions from the point of view of characteristic variables.  Ocean Modelling 9, 3 (2005), 231–252.
  * [3] Bozonnet, C., Desjardins, O., and Balarac, G. Traction open boundary condition for incompressible, turbulent, single- or multi-phase flows, and surface wave simulations.  Journal of Computational Physics 443 (2021), 110528.
  * [4] Delaney, K. An Adaptive Mesh Refinement Solver for Multiphase Incompressible Flows with Large Density Ratios.  PhD thesis, George Washington University, 2014.
  * [5] Dhruv, A. A Multiphase Solver for High-Fidelity Phase-Change Simulations over Complex Geometries.  PhD thesis, George Washington University, 2021.
  * [6] Dhruv, A. Github: Lab-Notebooks/Flow-Boiling-Performance, Nov. 2023.
  * [7] Dhruv, A. Github: Lab-Notebooks/Outflow-Forcing-Revised, Nov. 2023.
  * [8] Dhruv, A., Balaras, E., Riaz, A., and Kim, J. A formulation for high-fidelity simulations of pool boiling in low gravity.  International Journal of Multiphase Flow 120 (2019), 103099.
  * [9] Dhruv, A., Balaras, E., Riaz, A., and Kim, J. An investigation of the gravity effects on pool boiling heat transfer via high-fidelity simulations.  International Journal of Heat and Mass Transfer 180 (2021), 121826\.
  * [10] Dhruv, A., and Dubey, A. Managing software provenance to enhance reproducibility in computational research.  Computing in Science & Engineering 25, 3 (2023), 60–65.
  * [11] Dhruv, A., and Hassan, S. M. S. Github: Lab-Notebooks/Outflow-Forcing-BubbleML, Nov. 2023.
  * [12] Dodd, M. S., and Ferrante, A. A fast pressure-correction method for incompressible two-fluid flows.  Journal of Computational Physics 273 (2014), 416–434.
  * [13] Dong, S. An outflow boundary condition and algorithm for incompressible two-phase flows with phase field approach.  Journal of Computational Physics 266 (2014), 47–73.
  * [14] Dong, S., Karniadakis, G. E., and Chryssostomidis, C. A robust and accurate outflow boundary condition for incompressible flow simulations on severely-truncated unbounded domains.  Journal of Computational Physics 261 (2014), 83–105.
  * [15] Dubey, A., Weide, K., O’Neal, J., Dhruv, A., Couch, S., Harris, J. A., Klosterman, T., Jain, R., Rudi, J., Messer, B., Pajkos, M., Carlson, J., Chu, R., Wahib, M., Chawdhary, S., Ricker, P. M., Lee, D., Antypas, K., Riley, K. M., Daley, C., Ganapathy, M., Timmes, F. X., Townsley, D. M., Vanella, M., Bachan, J., Rich, P. M., Kumar, S., Endeve, E., Hix, W. R., Mezzacappa, A., and Papatheodore, T. Flash-x: A multiphysics simulation software instrument.  SoftwareX 19 (2022), 101168.
  * [16] Gibou, F., Chen, L., Nguyen, D., and Banerjee, S. A level set based sharp interface method for the multiphase incompressible Navier-Stokes equations with phase change.  Journal of Computational Physics 222, 2 (2007), 536–555.