Computers & Fluids 87 (2013) 115–131

Contents lists available at SciVerse ScienceDirect

Computers & Fluids

j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

A parallel adaptive mesh method for the numerical simulation

of multiphase flows
Joseph M. Rodriguez a,b, Onkar Sahni b, Richard T. Lahey Jr. c, Kenneth E. Jansen d,⇑
a
Lockheed Martin – KAPL, Inc., P.O. Box 1072, Schenectady, NY 12309, USA
b
Scientific Computation Research Center, Rensselaer Polytechnic Institute, 110 8th St., Troy, NY 12180, USA
c
Center for Multiphase Research, Rensselaer Polytechnic Institute, 110 8th St., Troy, NY 12180, USA
d
University of Colorado Boulder, UCB 492, Boulder, CO 80309, USA

a r t i c l e i n f o a b s t r a c t

Article history: Methodology for performing large scale parallel numerical simulations of two-phase flow is presented
Received 19 April 2012 that uses a three-dimensional (3-D) finite element method (FEM) solver, PHASTA, to perform the simu-
Received in revised form 30 March 2013 lation and a parallel mesh adaptation code, phParAdapt, to iteratively refine and coarsen regions of the
Accepted 1 April 2013
solution domain. This version of PHASTA uses a stabilized FEM to discretize the 3-D Incompressible
Available online 18 April 2013
Navier–Stokes (INS) equations plus a level set method to capture the two-phase interface. The phPar-
Adapt code takes the current parallel mesh along with the PHASTA flow solution and, using various
Keywords:
user-specified correction indicators as input, performs a parallel mesh refinement/coarsening procedure.
Multiphase flow
Adaptive mesh
Utilizing local mesh modification operators, element sizes are adjusted to equally distribute the errors
FEM across the distributed mesh during the mesh refinement/coarsening procedure. Numerical simulations
DNS ranging from simple canonical test problems up to an experimental annular steam/water flow condition
Massive parallel processing (MPP) were performed to illustrate this methodology. The annular flow simulation captured the instabilities
that generate both the ripple and roll wave like structures on the steam/water interface and the droplet
entrainment mechanisms expected in annular two-phase flows.
Ó 2013 Published by Elsevier Ltd.

1. Introduction ‘‘data’’ to support the development of physically-based closure

laws to be used in advanced-generation 3-D, two-fluid computa-
Multiphase flow occurs in many important industrial applica- tional multiphase fluid dynamic (CMFD) models [24].
tions such as phase-change heat exchangers, oil well pipelines, fos- This paper presents an adaptive mesh methodology for per-
sil-fired boilers and nuclear reactors. Many flow regimes are forming large scale parallel numerical simulations of multiphase
possible, but the most common one in power production and uti- flows. PHASTA [17,46] is a 3-D, FEM code which uses a level
lization technology is annular flow. This flow regime is quite com- set algorithm for interface tracking. The phParAdapt mesh-adapta-
plex and, as shown schematically in Fig. 1, it is characterized by a tion code [8,2,32] is used in conjunction with PHASTA to iteratively
droplet laden gas core and a liquid film on the conduit wall. The refine and/or coarsen the computational mesh based on various
dynamics occurring on the wavy interface between the liquid film correction indicators. Results are presented in this paper which
and gas core are important to the process of liquid droplet entrain- demonstrates the ability of these numerical methods to capture
ment into the gas core, the lateral force responsible for keeping the the physical mechanisms of annular two-phase flow.
liquid film on the conduit walls, and the turbulence structures in
the liquid and gas phases which effect the pressure drop. Due to
the difficulty in obtaining detailed local data in annular flow exper- 2. Discussion
iments, researchers currently do not have a fundamental under-
standing of these interfacial processes. With the advancement in As a result of the importance of annular two-phase flow,
modern computing capabilities, Direct Numerical Simulation numerous researchers [16,15,40,18,44,9] have developed phenom-
(DNS) can be used to help us understand the physical mechanisms enological models describing the physical processes of annular
in annular flow (and other flow regimes) and generate numerical flow based more on empiricism than detailed physical understand-
ing. In mechanistically based two-fluid computational multiphase
fluid dynamics (CMFD) model formulations [35,20,22], which at-
⇑ Corresponding author. tempt to capture the dynamics of two-phase flow, one uses closure
E-mail address: [email protected] (K.E. Jansen). laws based on the physical transport process at the interface.

0045-7930/$ - see front matter Ó 2013 Published by Elsevier Ltd.

http://dx.doi.org/10.1016/j.compfluid.2013.04.004
116 J.M. Rodriguez et al. / Computers & Fluids 87 (2013) 115–131

code also uses parallel mesh tools developed at RPI [8] and the par-
allel MeshSim libraries of Simmetrix Inc. [36].
Flow High Velocity Sections 3 and 4 present the governing equations and the FEM
Gas Core implementation in PHASTA, respectively. Canonical problems
demonstrating the ability of PHASTA to track an interface are dis-
Liquid film cussed in Section 5, the parallel mesh adaptation code is described
in Section 6 and an annular two-phase flow simulation is given in
Section 7.
Entrained gas Liquid droplets
bubbles in film in gas core
3. Governing equations
Fig. 1. Annular two-phase flow.

3.1. Incompressible Navier–Stokes (INS) equations

Significantly, having accurate closure laws is important for both
The spatial and temporal discretization of the INS equations
steady and, especially, transient calculations.
within PHASTA has been given previously by Whiting and Jansen
Given that detailed experimental methods are not available to
[46] and Nagrath et al. [31]. This section follows these formula-
fully characterize the wavy interface, one may use DNS ‘‘data’’ to
tions. The so-called strong form of the INS equations is given by:
support fundamental two-fluid model development. That is,
although DNS is not feasible as a practical tool for the design and Continuity : ui;i ¼ 0 ð1Þ
analysis of multiphase systems (due to the large computational ex-
penses), Lakehal [25] and Lahey [23] have shown how DNS may be Momentum : qui;t þ quj ui;j ¼ p ^ij;j þ ^f i
^ ;i þ s ð2Þ
used for closure model development. This approach is consistent
with a ‘‘roadmap’’ for basic research in multiphase flow which where q is density, ui is ith component of velocity, p ^ is static pres-
^ij is the stress tensor, and ^f i represents the component of the
sure, s
was assembled by leading experts in the field [14], which supports
the use of DNS (and improved experimental methods) to help body force along the ith coordinate. For incompressible flow, the
understand the microscopic mechanisms that affect the macro- stress tensor is related to the strain rate tensor, Sij, as:
scopic behavior of multiphase flows. The vast majority of DNS cal-
s^ij ¼ 2lSij ¼ l½ui;j þ uj;i  ð3Þ
culations done to date for two-phase flows have focused on bubbly
flow [6,42,31]. In contrast, very little has been done for annular Using the Continuum Surface Force (CSF) model of Brackbill
flows due to the computational costs associated with accurately et al. [5], the surface tension force is computed as a local interfacial
capturing the wave structure at the interface. Stability limits for force density, which is included in ^f i .
stratified two-phase flows have been investigated by Cao et al.
[7], Lakehal et al. [26] and Fulgosi et al. [12] but these were done 3.2. Level set method
for relatively low Reynolds numbers such that wave breaking
and entrainment mechanisms were avoided. Nevertheless, some The level set method of Sussman [39,38,37,33] involves model-
interesting annular flow simulations have been performed previ- ing the interface as the zero level set of a smooth function, /,
ously. For example, Li and Renardy [27] simulated unsteady, axi- where / represents the signed distance from the interface. Hence,
symmetric oil/water ‘‘annular flows’’ using a Volume of Fluid the interface is defined by / = 0. The scalar /, the so-called first
(VOF) method. Their computation was made at Reynolds numbers scalar, is convected within a moving fluid according to,
less than 10 and with a density ratio near unity. Consequently,
wave breaking and droplet entrainment processes did not occur. D/ @/
¼ þ u  r/ ¼ 0 ð4Þ
Also, using a level set method, Fukano and Inatomi [11] have sim- Dt @t
ulated the formation of the liquid film around the tube wall in hor-
izontal air/water annular flow. They specifically investigated the where u is the flow velocity vector. Phase-1, typically the liquid
source of the levitation force that causes the liquid film to spread phase, is indicated by a positive level set, / > 0, and phase-2 by a
to the upper surface of the tube. Although their simulation was negative level set, / < 0.
able to capture a physical mechanism which could explain the ob- Evaluating the jump in physical properties using a step change
served liquid film mobility, it was performed on a fairly coarse grid across the interface leads to poor computational results. Therefore,
which was incapable of resolving the droplets entrained in the gas the properties near an interface are defined using a smoothed
core or any bubbles ingested into the liquid film. Heaviside kernel function, He, given by:
The current investigation uses the PHASTA code with a level 8
>
<  0; / < e
set algorithm to simulate various two-phase flows and the phPar-  
Adapt code to perform iterative mesh adaptation based on suitable He ð/Þ ¼ 12 1 þ /e þ p1 sin pe/ ; j/j < e ð5Þ
>
:
correction indicators. PHASTA solves the Incompressible Navier– 1; />e
Stokes (INS) equations in three dimensions using a stabilized
FEM. A second-order accurate and stable generalized-a time inte- where the ‘‘interface’’ thickness was defined by 2e, as shown in
grator [17] was applied to the INS equations and used to march Fig. 2.
the solution in time [46]. PHASTA incorporates the level set meth- The properties around the interface are thus defined as,
od of Sussman [39,37,38,33] to resolve the interfaces in two-phase q2
flows by modeling an interface as the zeroth level set of a smooth qð/Þ ¼ q1 ðHe ð/Þ þ ð1  He ð/ÞÞÞ
q1
function. The Continuum Surface Force (CSF) model of Brackbill l ð6Þ
et al. [5] was used to represent the local surface tension force as lð/Þ ¼ l1 ðHe ð/Þ þ 2 ð1  He ð/ÞÞÞ
l1
a volume force applied over an interface of finite thickness. The
phParAdapt code performs parallel mesh adaptation in which the Although the solution may be reasonably good in the immedi-
local mesh size is constructed using selected correction indicators ate vicinity of the interface, the distance field may not be correct
like residual error, shear, and/or proximity to the interface. The throughout the domain due to varying fluid velocities throughout
 J.M. Rodriguez et al. / Computers & Fluids 87 (2013) 115–131 117

where He is the Heaviside function given in Eq. (5) and dk is the va-
Fluid 1 lue of the distance field, d, for iteration-k. Small changes in the vol-
ume from pseudo time step s0 to sk may be approximated by,
Fsv R 0
V kij  V 0ij  ðsk  s0 Þ Xe @H@e ðd
s dX
Þ

ρ1, p1 R 0 k 0
ð12Þ
 Xe H0e ðd Þðd  d ÞdX

where H0e is given by:

(
2ε 0 0; jdj > e
He ð/Þ ¼ rHe ð/Þ ¼    ð13Þ
1 1
2 e
þ 1e cos ped ; jdj 6 e
n
k
interface Hence, to minimize the deviation in volume, each new dij is con-
strained to satisfy the following:
Z
Fsv 0 k
H0e ðd Þðd  d ÞdX ¼ 0
0
ð14Þ
Xe
Fluid 2
This is accomplished by updating the value of the current level set
. ~k ; to the new level set field, dk , where dk is computed from,
scalar, dij ij ij

k ~k þ k ðsk  s0 ÞH0 ðd0 Þ

Fig. 2. Interface’s transition region of finite thickness 2e between fluids 1 and 2. dij ¼ d ij ij e ð15Þ

The scaling function kij is assumed constant within each cell,

the flow field which distort the level set contours (i.e., see Eq. (4)). and it is computed from Eqs. (14) and (15) to be:
Thus the level set is corrected with a redistancing operation by R  
0 ~k 0
solving the following PDE:  X e H0e ðd Þ dsk d
s0 dX
kij ¼ R 0
ð16Þ
@d 0
Xe ðH e ðd ÞÞ2 dX
¼ Sð/Þ½1  jrdj ð7Þ
@s k ~k , and
Note that when the distance field has converged, dij ¼ d ij
where d is a scalar (i.e., the second scalar) that represents the cor- kij = 0. In PHASTA, the volume constraint correction occurs at the
rected distance field and s is a pseudo time over which the PDE is global level after each redistancing iteration. The calculation of kij
solved to steady-state for each physical time step, Dt. This may be occurs at the element level and is projected back to the global level
more clearly expressed as the following transport equation: to be used by Eq. (15).
@d
þ w  rd ¼ Sð/Þ ð8Þ 3.4. Continuum Surface Force (CSF) model
@s
where
Surface tension is the result of sharp changes in the molecular
rd forces of attraction at an interface due to the discontinuous
w ¼ Sð/Þ ð9Þ
jrdj changes in fluid properties. It is difficult to accurately model this
local force for complex geometries. Instead, the CSF model of
The second scalar, d, is originally assigned the level set field, /, Brackbill et al. [5] was utilized to model the local surface tension
and is convected with a pseudo velocity w, where S(/) is defined force as a continuous, three-dimensional volume force which
as: smoothly varies across an ‘‘interface’’ of finite thickness.
8
> 1; / < e The pressure change across an interface, defined by fluids 1 and
< p/
Sð/Þ ¼ /
þ 1
sin ; e </<e ð10Þ 2, is given by the momentum jump condition at the interface:
> e p e
:
1; />e ðp2  p1 þ rj þ n ^ i  s1  n ^ i  s2  n
^i þ n ^ i Þn^i

Note that, in a continuous solution of Eq. (8), the zeroth level set þ ðn ^
^ i  s1  t  n ^ ^
^ i  s2  t  jrs rjÞt ¼ m _ 00 ðv 2  v 1 Þ
2 ð17Þ
or interface, / = 0, does not move since its convecting velocity, w, is
where pj, vj, and sj are the pressure, velocity vector, and viscous shear
zero. Solving the second scalar to steady-state restores the distance
stresses at the interface for fluid-j, r is the surface tension coefficient,
field to rd = ±1 but does not alter the location of the interface de-
j is the local curvature of the interface, m_ 002 is the mass flux due to
fined by Eq. (4). The scalar level set field, /, is then updated using
phase change for fluid 2, n ^ i is the unit normal vector to the interface
the steady solution of d.
and ^t is the unit tangent vector along the interface. For the special case
of inviscid, adiabatic conditions Eq. (17) reduces to:
3.3. Volume constraint on level set redistancing step
^ i  ðjrs rjÞ^t ¼ 0
ðp2  p1 þ rjÞn ð18Þ
Sussman et al. [37] and Sussman and Fatemi [38] have proposed
an additional constraint to be applied during the redistancing step The normal component for constant surface tension further re-
to help ensure the interface (/ = 0) does not move when Eq. (8) is duces to the well-known Laplace expression:
solved in a discretization. They found imposing this constraint im-
proved the convergence of the redistancing step while maintaining p2  p1 ¼ rj ð19Þ
the original zero level set. The essence of the constraint is to pre-
where the local curvature of the interface is defined as the diver-
serve the original volume (i.e., mass) of each phase during the
gence of the normal along the interface:
redistancing step. The volume of each phase within a cell, Vij, is de-
fined as: j ¼ rs  n^ ¼ ½I  n^ n^   r  n^ ð20Þ
Z
k
V kij ¼ He ðd ÞdX ð11Þ The surface force per interfacial area at a given position along
Xe the interface is thus given by:
118 J.M. Rodriguez et al. / Computers & Fluids 87 (2013) 115–131

f s ¼ rjn
^ dðx  xi Þ ð21Þ Skh ¼ fujuð; tÞ 2 H1 ðXÞN ; t 2 ½0; T; ujx2Xe 2 Pk ðXÞN ; uð; tÞ ¼ g on Cg g ð27Þ

where d(x  xi) is the Dirac delta function which restricts the sur- Wkh ¼ fwjwð; tÞ 2 H1 ðXÞN ; t 2 ½0; T; wjx2Xe 2 Pk ðXÞN ; wð; tÞ ¼ 0 on Cg g ð28Þ
face tension force to the interface, which is at x = xi. Consider the Pkh ¼ fpjpð; tÞ 2 H1 ðXÞN ; t 2 ½0; T; pjx2Xe 2 P k ðXÞN g ð29Þ
interface shown in Fig. 2, which is defined by the level set method,
where the interface is given a finite thickness of 2e. The local surface where h represents the characteristic length of the parameterized
tension force given in Eq. (21) may be expressed as a force density, domain, g is an approximation of the prescribed boundary condition
Fsv, such that, in the discretized domain, and P k ðXe ) is the piecewise polynomial
Z Z space defined on the element Xe for order k P 1. Note the velocity
lim Fsv dV ¼ f s dA ð22Þ and pressure variables use the same polynomial space, which is en-
e!0 DV DA abled through the stabilized formulation used in PHASTA.
A smoothing function was chosen that smoothly varies the The weak form is generated by dotting Eqs. (1) and (2) on the
properties across the interface. In particular, the smooth Heavi- left by the weight functions defined above and integrating over
side function, He, as defined in Eq. (5), is used in PHASTA for this the domain. Following the stabilized formulation of Taylor et al.
purpose. The foundation of the CSF is to use the gradient of the [41] the residuals of the continuity and momentum equations
smoothing function, DHe, defined in Eq. (13), to vary the surface are added together. The Galerkin finite element formulation is gi-
tension force across the finite thickness of the interface. This ven by finding u 2 Skh and p 2 P kh , such that:
function must zero out the force at the boundaries of the finite
BG ðwi ; q; ui ; pÞ ¼ 0
interface but provide a net force over the volume equivalent to R
the local area force, fs. Additional choices for smoothing functions BG ðwi ; q; ui ; pÞ ¼ X ½q;i ui þ wi fui;t þ uj ui;j  fi g
R ð30Þ
have been explored by Brackbill et al. [5] and Williams et al. [47]. þ wi;j fpdij þ sij gdX þ C qui n ^ i  wi ðsij  pdij Þn
h
^ j dC
In the limit, as the interpolation region around the interface
shrinks, rHe for any of these functions becomes a Dirac delta for all w 2 Wkh and q 2 P kh .
function: The pressure, viscous stress, and continuity terms are inte-
grated by parts eliminating all second order derivatives. The
limrHe ð/Þ ¼ dðx  xi Þ ð23Þ boundary integral term in Eq. (30) is a by-product of this and only
e!0
exists on Ch. Stabilization terms are added to the standard Galerkin
Since the interface is defined in terms of contours of the smoothing method to correct known stabilization issues for advection domi-
function, the unit normal at the interface must also be defined in nated problems. The resulting stabilized formulation used in PHAS-
terms of this smoothing function: TA is given by finding u 2 Skh and p 2 P kh , such that:
re H r/
^¼
n ¼ ð24Þ Bðwi ; q; ui ; pÞ ¼ BG ðwi ; q; ui ; pÞþ < stabilization terms >¼ 0 ð31Þ
jre Hj jr/j
R
Due to the differentiability of the smoothing function, the inte- Bðwi ; q; ui ; pÞ ¼ X ½q;i ui þ wi fui;t þ uj ui;j  fi g þ wi;j fpdij þ sij gdX
R
gral of the surface tension force in Eq. (22) may be expressed as,
Ch qui ni  wi ðsij  pdij Þnj dC
þ ^ ^
Z Z
rHe ð/Þ X
nel
R
f s dA ¼ lim rjð/Þ dV ð25Þ s þ q;i gLi þ sC wi;i ui;i dXe
DA e!0 DV ½He ð/Þ Xe ½ M fuj wi;j
e¼1
We note that according to Eq. (5), the jump in He across the X
nel
R D D D
interface, [He], is unity. Comparing Eq. (25) to Eq. (22), we see that þ Xe ½wi uj ui;j
uj wi;j uk ui;k dXe
þs
the surface tension force density is given by: e¼1

ð32Þ
Fsv ¼ rjð/ÞrHe ð/Þ ð26Þ
for all w 2 Wkh
and q 2 P kh ,
where Li represents the residual of the
The choice of the smoothing function kernel can impact the
accuracy and stability of the numerical computations. A discussion strong form of the ith momentum equation, Eq. (2);
of the desirable and required features of smoothing kernels has Li ¼ ui;t þ uj ui;j þ p;i  sij;j  fi ð33Þ
been given by Williams et al. [47].
The third line of Eq. (32) is the standard Streamwise Upwind Pet-
rov–Galerkin (SUPG) stabilization terms for incompressible flows
4. Finite element formulation [10]. The first term in the last line was added by Taylor et al. [41]
to overcome the lack of mass conservation introduced as a conse-
4.1. Incompressible Navier–Stokes (INS) equations – the numerical quence of the momentum stabilization in the continuity equation
D
method and is given by uj ¼ sM Li . The second term on this line was intro-
duced to stabilize this new advective term. The stabilization param-
The INS equations apply over a spatial domain defined by eters have been described in detail by Whiting and Jansen [46].
X  RN , where N equals three for our case. This domain is com- The weight functions and solution variables in Eq. (32) are ex-
posed of the interior, X, and the boundary, C. Note the boundary panded in terms of the finite element basis and the spatial integrals
is further subdivided into regions with natural boundary condi- are evaluated using Gauss quadrature. Since the weight functions
tions, Ch, and regions with essential boundary conditions, Cg, such are arbitrary, this results in a system of first-order, non-linear or-
that C = Ch [ Cg. The finite element formulation is based on the dinary differential equations. Finally, this system of non-linear or-
so-called weak form of Eqs. (1) and (2). To derive the weak form, dinary differential equations is discretized in time via a
continuous function spaces of trial solutions, S, and weighting generalized-a time integrator [17], resulting in a non-linear system
functions, W, are defined with square-integrable values and deriv- of algebraic equations. This system is then linearized using New-
atives (H1 functions) on X. The domain is discretized into nel finite ton’s method, which yields a linear algebraic system of equations
elements defining the element domain, Xe . Consequently, finite- that is solved (at each time step) and the solution is updated for
dimensional approximations, or subspaces, of S and W, specified each of the Newton iterations. The linear algebra solver of ACUSIM
as Skh and Wkh , respectively, are defined as: [1] is used to solve the linear system of algebraic equations.
 J.M. Rodriguez et al. / Computers & Fluids 87 (2013) 115–131 119

Fig. 3. Schematic of two-dimensional dam break problem. Water was initially confined in a square box. The computational domain was 5a  1.25a, and a 400  100 uniform
mesh was used.

Time = 0.0
T* = 0.0

Time = 0.10
T* = 1.31

Time = 0.20
T* = 2.621

Time = 0.30
T* = 3.931

Fig. 4. Interface and velocity vectors at various times for the 2-D dam break using a 400  100 grid.

4.2. Level set method

The level set scalar equations given by Eqs. (4) and (8) may be 2D Dam Break: Surge Front -vs- Time
cast into the following form: a=2.25 in
6.0
/;t þ Ai /;i þ S ¼ 0 ð34Þ
5.0 Martin & Moyce, 1952
where Ai is ui for scalar-1, /, and Sð/Þ Dd=jDdj for scalar-2, d, and S is
zero for scalar-1 while S(/) for scalar-2. The level set equations 4.0 PHASTA 400x100 dt=0.0001
Z* = s/a

represented by Eq. (34) are solved in the same manner as for the
3.0
INS equations given in Section 4.1. Temporal discretization is con-
trolled via the generalized-a time integrator while the same stabi- 2.0
lized finite element method formulation is used for the spatial
discretization. The weak form of Eq. (34), with the stabilization 1.0
terms added, is given by:
0.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0
Z nel Z
X
ðw/;t þ wAi /;i þ wSÞdX þ T e
L wsð/;t þ Ai /;i þ SÞdX ¼ 0 T* = nt*(g/a)0.5
X e¼1 Xe
Fig. 5. 2-D dam break: non-dimensional surge front position versus non-dimen-
ð35Þ sional time.
120 J.M. Rodriguez et al. / Computers & Fluids 87 (2013) 115–131

2D Dam Break: Height -vs- Time weight function and solution variable are expanded in terms of lin-
a=2.25 in ear basis functions:
1.2
X
nnp X
nnp
1.0 Martin & Moyce, 1952 w¼ wB NB ; /¼ /A N A ;
B¼1 A¼1
PHASTA 400x100 dt=0.0001
0.8 X
nnp X
nnp
H* = h/a

/;i ¼ /A NA;i ; /;t ¼ /A;t NA : ð36Þ

0.6 A¼1 A¼1

0.4 After evaluating the integrals using Gauss quadrature, Eq. (35)
may be expressed as a non-linear ordinary differential equation
0.2 (ODE) given by:

0.0 @R nþ1
0.0 1.0 2.0 3.0 4.0 5.0
nþ1
Du;t ¼ R ð37Þ
τ = t*(g/a)0.5 @ u;
|ffl{zffl}t
M
Fig. 6. 2-D dam break: non-dimensional water column height versus non-dimen-
sional time.
where R is the residual error of Eq. (35). PHASTA can solve Eq. (37)
implicitly or explicitly for either scalar. When solving implicitly, the
generalized-a method described previously is used for the first level
where LT is the advective operator, ui @=@xi , s is a stabilization param- set scalar, /, while a Backward Euler approach is used for the second
eter, and w is the weight function belonging to the weight space level set scalar, d (i.e., the redistancing scalar). A Forward Euler
w 2 Wh. The solution u belongs to the solution space / 2 Vh. The method is used as the explicit solver for both scalars. Once the

0.10
0.08
0.06
0.04
0.02
0.00
-0.02
-0.04
-0.06
-0.08
-0.10
0.0 0.5 1.0 1.5 2.0

Fig. 7. Schematic of 2-D traveling solitary wave. Water was started at rest with a Bousinesq profile, of height Ao, with a hydrostatic pressure balance. The computational
domain of 20h  2h was discretized with a 200  120 non-uniform mesh.

Fig. 8. Snap shots of the solitary wave calculation. The solitary wave is considered generated 0.6 s after the free fall of the initial profile. It is designated with the t0 = 0.0
marker. Grid: 200  120, wave speed = 1 m/s.
 J.M. Rodriguez et al. / Computers & Fluids 87 (2013) 115–131 121

Fig. 9. Snap shot of the velocity vectors of solitary wave at t0 = 0.4 s. The wave is traveling to the right at a wave speed = 1 m/s.

PHASTA and the results are given herein. The problem setup, along
2D Solitary Wave: Wave Amplitude with a schematic indicating the free surface before and after the
0.30 dam break, is given in Fig. 3. The problem consists of a column of
water being held in place by a membrane in a box. At t = 0, the
0.25
membrane was broken allowing the water column to collapse.
Amplitude (Amax/h)

The location of the surge front, s, and the height of the water col-
0.20
umn, h, was measured as a function of time. The numerical prob-
lem has no-slip wall boundary conditions in the x and y
0.15
Analytical Perturbation Sol'n by Mei
dimensions and periodicity in the z dimension (i.e., it was a 3-D
(1989) simulation with only one node in the z direction). The air and water
0.10 PHASTA - Surf. Tension
were at atmospheric pressure and 68 °F and their velocities were
0.05
PHASTA - No Surf. Tension initially zero, with a hydrostatic pressure profile in the water col-
umn. The domain was modeled using a uniform 400  100 mesh
0.00 (Dx = Dy = 0.000714 m). The ‘‘interface’’ thickness, for advection
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 of the interface using the level set method, was 2e = 0.004 or about
Time (s), t’ 5.5Dx. The ‘‘interface’’ thickness for the redistancing operation was
kept at 2e = 0.002 or about 3Dx. The level set equations were
Fig. 10. Solitary wave amplitude. (Grid: 200  120).
solved explicitly with the CFL number for all equations not exceed-
ing 0.3.
_ is solved the
change in the time derivative of the second scalar, Dd, Snapshots of the free surface position and velocity vectors in the
second scalar field is corrected using: computational domain are provided at various times in Fig. 4.
Computations were made to the point where the water climbs
¼ d þ Dd_ Ds
nþ1 n
d ð38Þ the adjacent wall. At the start of the dam break, the water acceler-
n+1 n
ates from its initial static condition due to the relatively large pres-
where d is the current solution, d is the old solution, and Ds is sure difference existing on the right side of the water column once
the pseudo time step for the second scalar’s partial differential the membrane is removed. The velocity of the surge front of the
equation (PDE), Eq. (8). The user defines the number of iterations, water continues to increase until the opposite wall is impacted.
or pseudo time steps, to solve the redistancing equation to stea- At each time, the maximum velocity was in the air just ahead of
dy-state. the surge front. The non-dimensional surge front location, Z, and
water column height, H, are plotted with the experimental data
5. Canonical problems of Martin and Joyce [29] in Figs. 5 and 6, respectively. PHASTA
did a good job predicting the data and also agrees with the numer-
Two canonical test problems were studied to assess the ability ical predictions presented by Kelecy and Pletcher [19] and Yue
of PHASTA to accurately capture the movement of interfaces. The et al. [49]. As can be seen, the predictions begin to deviate a little
first problem numerically tracks the front propagation of a 2-D from the data after T > 2.0. Yue et al. [49], who also used a level
dam break problem that has been experimentally and analytically set method, also found this deviation and attributed it to using
studied. The second problem tracks the propagation of a traveling too large an interface thickness, and they found the best numerical
solitary wave. The results are compared to data and analytical solu- prediction was when using 2e = 2Dx (note, the PHASTA computa-
tions. PHASTA was able to match the theoretical and experimental tion presented herein used 2e = 5Dx). In addition, PHASTA did a
results for both problems, without the need for mesh adaptation, good job conserving mass (i.e. within ±0.6%) even for a relatively
while incurring minimal mass error. coarse spatial mesh (i.e. 400  100).

5.1. Analysis of a 2-D dam break

5.2. Analysis of a solitary wave
The 2-D dam break problem was studied experimentally by
Martin and Joyce [29] and numerically investigated by Kelecy A 2-D traveling solitary wave, as analyzed by Yue et al. [49], was
and Pletcher [19] and Yue et al. [49]. It was also simulated using numerically simulated and compared to the analytical theory of
122 J.M. Rodriguez et al. / Computers & Fluids 87 (2013) 115–131

Fig. 11. Schematic of three-dimensional convecting sphere. A ‘‘red water/green water’’ problem with surface tension. The computational domain was 0.5  0.2  0.2 m.

3D View of Sphere Planar Cut of Sphere Time

Step

405

405 after
adapt

868

Fig. 12. Transport of sphere through the domain using adaptive-mesh approach: element size = 0.005 m near interface and 0.02 m elsewhere, mesh size 410,000
tetrahedrons. The black line indicates the interface. Mesh displayed in upper half of planar cuts while the color spectrum in the lower half indicates pressure (Pa).

pffiffiffiffiffiffi
Mei [30]. The solitary wave was generated by releasing water, ini- Cw ¼ gh 1:0 m=s ð40Þ
tially at rest, having a Boussinesq profile, Whitham [45], given by,
For this problem the initial amplitude, Ao, was selected to be
, pffiffiffiffiffiffiffiffiffiffiffiffiffi !!2
3Ao =h x 0.4 h. The domain was modeled with a 200  120 non-uniform grid
AðxÞ ¼ Ao cosh ð39Þ where the x-dimension had 200 elements of uniform width
2 h
(Dx = 0.1 h), while the y dimension had two regions of uniform
This profile is shown in Fig. 7, where h is the depth of the water width, and a periodic boundary condition was used across a single
and Ao is the initial height of the water on the left wall. A hydro- element in the z direction. The region containing the solitary wave
static pressure profile was applied across the domain. As seen in had 100 elements of width Dy = 0.01 h while the region outside the
Fig. 7, the water will fall under gravity to create a solitary wave. solitary wave had 20 elements of width Dy = 0.05 h.
For the case of zero bulk velocity for air and water, where Ao/ The solitary wave shown in Fig. 7 was fully formed at t = 0.6 s
h 1, the dispersion relation for shallow water waves [45] gives and this is marked as the t0 = 0.0 waveform in Fig. 8, which had
the wave speed, Cw, as: an amplitude of 0.185 h and a wave speed of 1.0 m/s, which agreed
 J.M. Rodriguez et al. / Computers & Fluids 87 (2013) 115–131 123

refinement/coarsening procedure to produce a computational

mesh which creates the local element length scales required to
adequately resolve the physics length scales. The process of mesh
adaptation utilizes local mesh modification operators [28] to im-
prove the mesh with respect to the desired mesh size field, and
executes over a distributed mesh on a parallel computer [8,2]. In
this process, the flow solution on the input mesh is mapped onto
the adapted mesh as the mesh is being modified locally, which pro-
vides updated PHASTA input to continue the simulation. Such a
solution strategy allows all the necessary steps of the simulation
to be carried out in a parallel manner without creating the bottle-
necks associated with a serial mesh adaptation process.
There are several options to compute the desired mesh size
Fig. 13. Schematic of the simulation domain. Length = 0.025 m = 1/8kw for RISO Run
field. In general, phParAdapt computes the new mesh size field
# 602. by comparing the local value of the correction indicators with
the target ones, determined by the user controlled parameters,
with Eq. (40). The velocity vectors at t0 = 0.4 given in Fig. 9 shows with the aim of equi-distributing the errors over the domain [4].
the typical vortex generated at the crest of the solitary wave. In locations where the local value of the correction indicators ex-
According to the perturbation theory of Mei [30], the amplitude ceeds the target value, the size field is refined, whereas in regions
of the solitary wave decreases with time due to viscous dissipation with these values lower than the target, the size field is coarsened.
according to, After adaptation, parallel load balancing is used to restore roughly
!1=2 the same number of elements to each processor.
1=4 mw Cwt The correction indicators used to determine the desired mesh
Amax ¼ A1=4
o max þ 0:08356 1=2 3=2
ð41Þ
Cw h h size field for the domain are either imported from the error indica-
tors output by PHASTA or computed within phParAdapt. Typical
where Ao_max was taken to be the amplitude of the solitary wave methods of generating the size field are:
once fully established at t = 0.6 s; hence, Ao_max = 0.185 h. The com-
puted wave amplitudes are compared with Eq. (41) in Fig. 10 and (1) PHASTA’s error indicators which provide the residual error
found to agree rather well, even if surface tension was not consid- in the momentum equations or the shear stress divergence
ered. The PHASTA predictions for times greater than t0 = 0.6 s began (which can be related to shear stress jumps across elements)
to be influenced by the opposite wall which is not in the analytical in the three coordinate directions. The analysis of a plunging
model. As for the dam break problem, the mass error associated liquid jet has indicated good success using these type error
with the calculation was low. indicators [13].
(2) Interpolation-based error indicators are determined by con-
6. Parallel mesh adaptation sidering the lowest-order error term which dominates when
interpolating a function. When using linear basis functions
6.1. Code description in the calculations, the interpolation errors are dictated by
the second derivatives that can be reconstructed to extract
The phParAdapt code takes the current parallel mesh along with error information [32]. Furthermore, directional error indi-
the PHASTA flow solution and error estimates as inputs and com- cators can be derived by using second derivatives in the form
putes the desired mesh size field required by the parallel mesh of a Hessian matrix [32]. Such a scheme, referred to as a

Fig. 14. Axial planar slice (yz plane) of initial mesh(including boundary layer mesh). Interface shown as dark bold line.
124 J.M. Rodriguez et al. / Computers & Fluids 87 (2013) 115–131

Time Film Interface Notes

Initial condition with

t = 0.0 perturbation on
sec
interface

Surface instabilities
t = 0.003 quickly form over
sec interface (on crests and
troughs of perturbation)

Ligaments form on wave

t = 0.006 crests and shear off
sec droplets into steam core

Two main crests remain

from initial perturbation
t = 0.009 that continue to break.
sec In between, the surface
resembles ripple wave
structures.

Fig. 15. 3D views of the interface from time steps 0 through 21,000 with uniform mesh. The steam/water interface is colored by the velocity magnitude.

Hessian strategy [21] is often applied in flow problems where / is the level set field indicating the signed distance to the
involving anisotropic solution features like shock waves interface.
and/or boundary layers. (1) A ramp proximity indicator that allows the user to specify
(3) A proximity indicator, E(/), which indicates proximity to the the desired mesh size field continuously. The general logic
interface, such as: sets the size near the interface to a specified low value,
h  i min_size, and set the size away from the interface to a spec-
8
pj/j ified
< 1  12 1 þ j/j 1
e þ p sin e ; j/j < e 8 max value, max_size:
Eð/Þ ¼ ð42Þ >
< min size; j/j < e
: ðmin sizeþmax sizeÞ
0; j/j > e Size ¼ ; e < j/j < 3e ð43Þ
>
:
2
max size; j/j > 3e
 J.M. Rodriguez et al. / Computers & Fluids 87 (2013) 115–131 125

Time Film Interface Notes

t = 0.012 Wave crests continue to

sec break.

One wave crest has

t = 0.015 collapsed, leaving only
sec
one left.

t = 0.018 Remaining wave crest is

sec attenuating.

No remaining wave
t = 0.021
sec crests from initial
perturbation.

Fig. 15. (continued)

where the thickness of the region of small elements around the capability of the phParAdapt code. This was a so-called ‘‘red water/
interface is given by e. green water’’ problem where there was surface tension, but the
It is also possible to use combinations of the above strategies properties of the two fluids were the same. The domain, a
(e.g., use of a proximity detector like (3) as a conditional sample 0.5 m  0.2 m  0.2 m box, was periodic in the axial direction
for error indication based on strategies (1) and/or (2) as was done (i.e., z direction) and had slip walls elsewhere. The properties were
by Galimov [13]). for water at standard temperature and pressure and were constant
throughout; however, since surface tension was applied at the
6.2. Example problem – a convecting sphere interface there was a pressure rise consistent with the Laplace
equation within the sphere. Two approaches were used to track
A simple convecting sphere in plug flow problem, shown in the sphere for one evolution through the domain. The first ap-
Fig. 11, was simulated to assess the grid refinement and coarsening proach, termed the ‘‘uniform’’ approach, used a uniformly sized
126 J.M. Rodriguez et al. / Computers & Fluids 87 (2013) 115–131

Time Center Plane (x=0)

t = 0.0
sec

t = 0.003
sec

t = 0.006
sec

t = 0.009
sec

Fig. 16. Velocity fringe plots over the center plane (x = 0) through from time steps 0–21,000 with uniform mesh. The bold dark line indicates the steam/water interface.

mesh (size = 0.005 m) where no mesh adaptation was performed. interface were refined and fine regions away from the interface
The second approach, termed the ‘‘adaptive-mesh’’ approach, used were coarsened. At step 868, the solution has again progressed to
phParAdapt to keep a mesh having a size of 0.005 m around the the point where it must be halted and adapted before it can again
sphere interface but used a mesh of size 0.02 m elsewhere. For be advanced further. There is naturally a tradeoff between the
the adaptive-mesh approach, a PERL script controlled the iterations thickness of the adaptivity band and adaptation frequency; if it is
between PHASTA and phParAdapt with PHASTA automatically very thick there the mesh will be much larger but the frequency
halting when the interface (i.e., zeroth level set) reached the limits of adaptation can be small while if it very thin then the number
of the fine mesh region. The uniform mesh had a mesh size of of elements remains smaller but adaptation must be more fre-
1,776,240 tetrahedrons which is more than four times the average quent. The choice of with e = 0.04 m in (43) represents a balance
size of the adapted mesh of 410,000 tetrahedrons. we have observed to be efficient.
The time step, Dt, was set such that CFL = 1.0, the pseudo time
step Ds was set such that CFLs = 0.25, els = 2D, elsd = 2.5D, where D
represents the local mesh size. The specified size field method gi- 7. Simulation of annular two-phase flow
ven in (43), with e = 0.04 m, min_size = 0.005 m, and max_si-
ze = 0.02 m, was used in the mesh adaptation process. Fig. 12 7.1. Problem setup
shows 3-D contours and planar cuts of the sphere at the initial,
middle, and final times during the initial transport of the sphere A physically interesting numerical simulation, similar to the
through the domain using the adaptive-mesh approach. At time experimental steam/water annular flow tests of Wurtz [48], was
step 405, as the sphere was crossing through the periodic plane, performed. This simulation was largely based on the flow condi-
the interface reached the coarse elements within the specified tol- tions of RISO run #602 but modeled a tube diameter of 0.019 m
erance of 0.010 m (i.e. 2Dmin). Fig. 12 shows the mesh before and which is 1.0 mm smaller than that used in the RISO test. RISO
after adaptation at time step 405, where coarse regions near the run #602 was conducted at 70 bar and the associated saturated
 J.M. Rodriguez et al. / Computers & Fluids 87 (2013) 115–131 127

Time Center Plane (x=0)

t = 0.012
sec

t = 0.015
sec

t = 0.018
sec

t = 0.021
sec

Fig. 16. (continued)

2
mixture temperature, had an inlet mass flux of 500 kg/s m , and an simulation, the average liquid and steam velocities of RISO run
exit quality of 30%. The flow had a liquid film turbulent Reynolds #602 were used. Note the average velocity can only be determined
number, Res = qusdm/l, of 792, where us is the friction veloc- once the void fraction is known. Using the experimentally
ity = (swall/q)0.5, swall is the experimentally measured wall shear measured mean film thickness, dm, from RISO run #602 and the
(= Pa), q is the density (=741 kg/m3), dm is the experimentally mea- assumption of negligible liquid droplet entrainment, the mean li-
sured mean film thickness (=0.94 mm), and l is dynamic viscosity quid and steam velocities were 2.64 m/s and 4.99 m/s,
(=9.13E5 Pa s). The computational domain, shown in Fig. 13, was respectively.
a 30° wedge of the 0.019 m diameter tube. Ideally, the domain The initial perturbation, Ddm, applied to the interface to expe-
length would be at least the experimentally measured mean wave- dite problem start-up was:
length, kw, of the disturbance waves, which for run #602 was
approximately 0.2 m; however, the domain length was limited 2p
Ddm ¼ b cos z ð44Þ
here to kw/8 to keep the problem size tractable during the code kp
assessment. Periodicity was enforced in the axial and circumferen-
tial directions. The centerline was not included in the domain due where b is the assumed amplitude, kp is the wavelength of the per-
to boundary condition limitations that arise in the flow solver. That turbation, and z is the axial coordinate. The perturbation’s wave-
is, the model was clipped at a radius of 0.0005 m, which is the size length was set equal to the critical wavelength for Helmholtz
of the coarsest elements expected in the model. instability:

1=2
2p ðql  qg Þg
7.2. Initial conditions
kc ¼ ; jc ¼ ð45Þ
jc r
Each phase was assigned a flat (i.e., plug) velocity profile with a where jc is the critical wave number, g is gravity and r is the sur-
magnitude equal to the average velocity of the phase. For this face tension.
128 J.M. Rodriguez et al. / Computers & Fluids 87 (2013) 115–131

(a) 2D slice of mesh

(b) Entrained Droplet (c) Liquid Film

Uniform
Mesh

Adapted
Mesh

Fig. 17. Uniform and adapted meshes at time step 6000. (a) Planar slice of adapted mesh (y–z plane) (b) expanded view of liquid droplet colored by pressure rise due to
surface tension (c) expanded view of mesh used for local liquid film.

8
7.3. The mesh >
> min size; j/j < e
>
<  
ðj/jeÞ
Size ¼ Min½2  min size; min size 1 þ e ; />e ð46Þ
The meshing strategy used was to retain a sufficiently refined >  
>
>
mesh around the steam/water interface, properly resolve the li- : Min½max size; min size 1 þ ðj/jeÞ ; / < e
e
quid film, and allow for a suitably coarse mesh in the steam core.
A boundary layer mesh was placed on the wall where the first where min_size = 0.000025 m, max_size = 0.0002 m, and
layer was within the viscous sublayer, with a height of y+ = 5. A e = 0.00015 m. This places a mesh of size min_size around the inter-
wall function was used, u+ = y+, to reduce the size of the near wall face over a thickness of 2e. The mesh size then linearly grows away
mesh. Note that Vassallo [43] and Azzopardi [3] found that the from the interface to a maximum size of 2  min_size in the liquid
turbulence in the liquid film, although enhanced by the effect of (/ > 0) and max_size in the steam (/ < 0). In addition, the mesh
waves on the interface, still obeys the general form of the law within 1 mm of the wall is kept at a size of min_size. A planar slice
of the wall, yet the coefficients for the log law region are different of the initial mesh, showing the initially perturbed interface, is seen
than for single phase flows. Placing the off-the-wall node at the in Fig. 14. The initial mesh was comprised of 57 million
outer edge of the viscous sublayer eliminates the need to know tetrahedrons.
the values of these wave-structure-dependent log-law
coefficients. 7.4. Simulation methodology
The size of the elements outside the boundary layer were con-
trolled via the user-defined size field option in phParAdapt using As discussed previously, a Perl script controls the simulation by
the size field criteria, similar to Eq. (43), given by: iteratively running PHASTA and phParAdapt. PHASTA runs until
 J.M. Rodriguez et al. / Computers & Fluids 87 (2013) 115–131 129

(a) droplet @ time

step=6000

(b) droplet @ time

step=6050

(c) droplet @ time step=6050

with adapted mesh

Fig. 18. Mesh adaptation around droplet (solid black line): (a) time step 6000: planar slice of adapted mesh around droplet (b) time step 6050: droplet interface moves to
edge of refined region (c) time step 6050: mesh adapted to place refined mesh around new droplet location.

the code detects that the interface (i.e., the zeroth level set) has tetrahedron mesh was used having an element
reached the edge of the refined mesh region placed around the size = 0.00005 m = 2  min_size. Three-dimensional views of the
interface. PHASTA then halts and the Perl script launches phPar- flow are provided in Fig. 15 for every 3000 time steps. Two-dimen-
Adapt to adapt the mesh based on the current location of the inter- sional velocity fringe plots of the center plane (x = 0) are provided
face and transfers the solution to the new adapted mesh. New for these same instances in Fig. 16, where the interface is depicted
PHASTA inputs are generated (i.e., geometry and solution data) by the solid bold line. Starting from the initial perturbation, surface
and PHASTA is restarted. The Perl script runs for a pre-determined instabilities quickly form along the interface. Those in between the
number of PHASTA-phParAdapt iterations, and PHASTA and phPar- crests of the initial perturbation are eventually damped out and
Adapt are run in parallel using the same number of processors. leave only ripple-wave-like structures. The instabilities on the per-
turbation wave crests tend to form ligaments that eventually shear
7.5. Results off and form entrained liquid droplets. Significantly, the predicted
phenomena are very similar to seminal experimental observations
The annular flow simulation presented herein was run using taken at the AERE Harwell [16].
2048 processors on a 70 TeraFLOP, 32,768 node IBM Blue-Gene The ability, and benefit, of the latest software to maintain a re-
supercomputer of the Computational Center for Nanotechnology fined mesh about the steam/water interface is indicated in Figs. 17
Innovations (CCNI) at RPI. PHASTA has demonstrated excellent and 18. From the expanded views in Fig. 17, which shows both the
scaling up to 32,000 processors on CCNI [34]. A uniform isotropic initial uniform and adapted meshes at time step 6000, we see that
130 J.M. Rodriguez et al. / Computers & Fluids 87 (2013) 115–131

in the adapted mesh there is a finer spatial mesh around the inter- [9] Fore LB, Beus SG, Bauer RC. Interfacial friction in gas-liquid annular flow:
analogies to full and transition roughness. Int J Multiphase Flow
face and within the liquid phase but it is coarsened elsewhere
2000;26(11):1755–69.
when comparing the adapted mesh to the original uniform mesh. [10] Franca LP, Frey SL. Stabilized finite-element methods. 2. The incompressible
The adapted mesh consists of 56 million elements, which is of sim- Navier–Stokes equations. Comput Methods Appl Mech Eng 1992;99(2–
ilar size to the original uniform mesh, which consisted of 53 mil- 3):209–33.
[11] Fukano T, Inatomi T. Analysis of liquid film formation in a horizontal annular
lion elements, but the adapted mesh significantly improves flow by DNS. Int J Multiphase Flow 2003;29(9):1413–30.
resolution of the interfacial and liquid film features. The mesh res- [12] Fulgosi M, Lakehal D, Banerjee S, De Angelis V. Direct numerical simulation of
olution around the interface has been increased by a factor of two turbulence in a sheared air-water flow with a deformable interface. J Fluid
Mech 2003;482:319–45.
which, if accomplished using an isotropic mesh, would have in- [13] Galimov A. An analysis of interfacial waves and air ingestion mechanisms.
creased the mesh size by a factor of eight (i.e., the initial 53 million Ph.D. thesis – engineering physics, Rensselaer Polytechnic Institute, Troy, NY;
elements would have increased to 424 million elements). Note the 2007.
[14] Hanratty TJ, Theofanous T, Delhaye J-M, Eaton J, McLaughlin J, Prosperetti A,
resolution of the boundary layer in the liquid film has also been et al. Workshop findings. Int J Multiphase Flow 2003;29(7):1047–59.
improved by the introduction of a boundary layer mesh on the [15] Henstock WH, Hanratty TJ. Interfacial drag and height of wall layer in annular
wall. To achieve this same resolution with isotropic tetrahedrons flows. AICHE J 1976;22(6):990–1000.
[16] Hewitt GF, Hall-Taylor NS. Annular two-phase flow. Oxford, New
would have increased the size of the computational mesh by about York: Pergamon Press; 1970.
45 million. Fig. 18 shows the ability of phParAdapt to maintain a [17] Jansen KE, Whiting CH, Hulbert GM. A generalized-alpha method for
refined mesh around the steam/water interface of a droplet. integrating the filtered Navier–Stokes equations with a stabilized finite
element method. Comput Methods Appl Mech Eng 2000;190(3–4):305–19.
Fig. 18a shows the droplet interface at time step 6000, after the
[18] Kataoka I, Ishii M. Entrainment and deposition rates for droplets in annular
mesh was adapted, and Fig. 18b shows the droplet near the edge two-phase flow. In: ASME-JSME thermal engineering joint conference,
of the refined mesh region 50 time steps later. The mesh was then Honolulu, Hawaii; 1983. p. 69–80.
adapted to keep fine mesh around the droplet in its new location. [19] Kelecy FJ, Pletcher RH. The development of a free surface capturing approach
for multidimensional free surface flows in closed containers. J Comput Phys
This new mesh is seen in Fig. 18c. Clearly, accurate simulations of 1997;138(2):939–80.
this type require both anisotropic boundary layer elements and [20] Kumar R, Maneri CC, Strayer TD. Modeling and numerical prediction of flow
adaptivity. Furthermore, flow solvers and adaptation capability boiling in a thin geometry. J Heat Transfer-Trans ASME 2004;126(1):22–33.
[21] Kunert G. Toward anisotropic mesh construction and error estimation in the
are required that function on massive parallel processing (MPP) finite element method. Numer Methods Partial Differen Equations
machines with highly scalable performance. 2002;18(5):625–48.
[22] Lahey RT. The simulation of multidimensional multiphase flows. Nucl Eng Des
2005;235(10–12):1043–60.
8. Conclusions [23] Lahey RT. The simulation of multidimensional multiphase flows. Nucl Eng Des
2005;235(10–12):1043–60.
[24] Lahey RTJ. On the direct numerical simulation of two-phase flows. Nucl Eng
An adaptive mesh methodology for performing large scale par- Des 2009;239(5):867–79.
allel numerical simulations of multiphase flows was presented that [25] Lakehal D. DNS and LES of turbulent multifluid flows. In: 3rd International
couples a FEM solver, an implemented level set method, with a symposium on two-phase flow modelling and experimentation, Pisa, Italy;
2004.
parallelized mesh adaptation algorithm that performs mesh adap- [26] Lakehal D, Fulgosi M, Yadigaroglu G, Banerjee S. Direct numerical simulation of
tation while equally distributing the size field ‘‘errors’’ to yield a turbulent heat transfer across a mobile, sheared gas–liquid interface. J Heat
well balanced distributed mesh. The benefit from this adaptive Transfer-Trans ASME 2003;125(6):1129–39.
[27] Li J, Renardy Y. Direct simulation of unsteady axisymmetric core-annular flow
mesh approach is shown in terms of a reduced overall mesh size with high viscosity ratio. J Fluid Mech 1999;391:123–49.
and/or improved resolution. [28] Li X, Shephard MS, Beall MW. 3D anisotropic mesh adaptation by mesh
modification. Comput Methods Appl Mech Eng 2005;194(48–49):4915–50.
[29] Martin JC, Joyce WJ. Part IV. An experimental study of the collapse of liquid
Acknowledgements columns on a rigid horizontal plane. Philos Trans Roy Soc Lond Ser A, Math
Phys Sci 1952;244:312–24.
[30] Mei CC. The applied dynamics of ocean surface waves. Singapore: River Edge,
The authors would like to acknowledge the financial support gi-
NJ; 1989.
ven this study by Lockheed Martin – KAPL and the computational [31] Nagrath S, Jansen KE, Lahey RT. Computation of incompressible bubble
resources provided by the Department of Defense (DOD) and dynamics with a stabilized finite element level set method. Comput
Rensselaer Polytechnic Institute (RPI) through the Computational Methods Appl Mech Eng 2005;194(42–44):4565–87.
[32] Sahni O, Müller J, Jansen KE, Shephard MS, Taylor CA. Efficient anisotropic
Center for Nanotechnology Innovation (CCNI). We further adaptive discretization of the cardiovascular system. Comput Methods Appl
acknowledge support from NSF 749152 and the US DOE DE- Mech Eng, John H. Argyris Memorial Issue. Part II 2006;195(41–43):5634–55.
FC02-06ER25769 for much of the flow solver and mesh adaptation [33] Sethian JA. Level set methods and fast marching methods. Cambridge
University Press; 1999.
development. Finally we acknowledge the support given us by [34] Shephard MS, Jansen KE, Sahni O, Diachin LA. Parallel adaptive simulations on
Simmetrix (mesh adaptation libraries) and ACUSIM (linear algebra unstructured meshes. J Phys: Conf Ser 2007 [012053].
libraries). [35] Siebert BW, Maneri CC, Kunz RF, Edwards DP. A four-field model and CFD
implementation for multi-dimensional, heated two-phase flows. Multiphase
flow ‘95, Kyoto, Japan; 1995. p. P2-19–28.
References [36] Simmetrix, http://www.simmetrix.com, Simmetrix Inc., 10 Halfmoon
Executive Park Drive, Clifton Park, NY 12065, USA; 2008.
[37] Sussman M, Almgren AS, Bell JB, Colella P, Howell LH, Welcome ML. An
[1] ACUSIM, http://www.acusim.com, 2685 Marine Way, Suite 1421, Mountain
adaptive level set approach for incompressible two-phase flows. J Comput
View, CA 94043, USA; 2008.
Phys 1999;148(1):81–124.
[2] Alauzet F, Li X, Seol E, Shephard M. Parallel anisotropic 3D mesh adaptation by
[38] Sussman M, Fatemi E. An efficient, interface-preserving level set redistancing
mesh modification. Eng Comput 2006;21(3):247–58.
algorithm and its application to interfacial incompressible fluid flow. Siam J Sci
[3] Azzopardi BJ. Turbulence modification in annular gas/liquid flow. Int J
Comput 1999;20(4):1165–91.
Multiphase Flow 1999;25(6–7):945–55.
[39] Sussman M, Fatemi E, Smereka P, Osher S. An improved level set method for
[4] Babuvska I, Rheinboldt WC. Error estimates for adaptive finite element
incompressible two-phase flows. J Comput Fluids 1998;27(5–6):663–80.
computations. SIAM J Numer Anal 1978;15(4):736–54.
[40] Tatterson DF, Dallman JC, Hanratty TJ. Drop sizes in annular gas–liquid flows.
[5] Brackbill JU, Kothe DB, Zemach C. A continuum method for modeling surface
AICHE J 1977;23(1):68–76.
tension. J Comput Phys 1992;100(2):335–54.
[41] Taylor CA, Hughes TJR, Zarins CK. Finite element modeling of blood flow in
[6] Bunner B, Tryggvason G. Direct numerical simulations of three-dimensional
arteries. Comput Methods Appl Mech Eng 1998;158(1–2):155–96.
bubbly flows. Phys Fluids 1999;11(8):1967–9.
[42] Tryggvason G, Bunner B, Esmaeeli A, Juric D, Al-Rawahi N, Tauber W, et al. A
[7] Cao Q, Sarkar K, Prasad AK. Direct numerical simulations of two-layer
front-tracking method for the computations of multiphase flow. J Comput
viscosity-stratified flow. Int J Multiphase Flow 2004;30(12):1485–508.
Phys 2001;169(2):708–59.
[8] de Cougny HL, Shephard MS. Parallel refinement and coarsening of tetrahedral
meshes. Int J Numer Methods Eng 1999;46(7):1101–25.
 J.M. Rodriguez et al. / Computers & Fluids 87 (2013) 115–131 131

[43] Vassallo P. Near wall structure in vertical air-water annular flows. Int J interfaces. Cambridge UK; New York: Cambridge University Press; 1999. xv,
Multiphase Flow 1999;25(3):459–76. 461 p.
[44] Whalley PB. Boiling, condensation, and gas–liquid flow. Oxford University [48] Wurtz J. An experimental and theoretical investigation of annular steam/water
Press; 1987. flow in tubes and annuli at 30 to 90 bar. RISO Report No. 372, RISO National
[45] Whitham GB. Linear and nonlinear waves. New York: Wiley; 1974. Laboratory; 1978.
[46] Whiting CH, Jansen KE. A stabilized finite element method for the [49] Yue WS, Lin CL, Patel VC. Numerical simulation of unsteady multidimensional
incompressible Navier–Stokes equations using a hierarchical basis. Int J free surface motions by level set method. Int J Numer Methods Fluids
Numer Methods Fluids 2001;35(1):93–116. 2003;42(8):853–84.
[47] Williams MW, Kothe DB, Puckett EG. Accuracy and convergence of continuum
surface-tension models. In: Shyy WW, Narayanan R, editors. Fluid dynamics at