7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Understanding and Prediction of Hydrodynamics of Multiphase Flow Using CFD
K. Hendrickson*, B. K. Campbell*, Y. Liu*, and R. Robertst
Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
t Chevron Corporation, ETC Advanced Production Systems, 1400 Smith Street, Houston, TX 77002
khendrk@mit.edu, bcampb@mit.edu, yuming@mit.edu, and roberrml@chevron.com
Keywords: direct numerical simulation, nonlinear interfacial instability, slug flow initiation and evolution
Abstract
This work describes the development of two computational capabilities designed to investigate the hydrodynamics
of multiphase (gas/liquid) flow in support of the design of deepwater pipelines and the development of advanced
flowmonitoring devices. Of particular interest is understanding the transition and nonlinear evolution of stratified
flows into liquid slugs in horizontal or incline pipes.
Two complementary computational capabilities are developed. The first scheme is a twodimensional method for the
computation of largescale flow in pipelines with L/D = 0(1 10,000), which is based on weakly nonlinear inviscid
multiphase flow. This capability (referred to as "LongPipe") is useful for understanding and predicting the initial
development and nonlinear evolution of interfacial disturbances over long distances. The other numerical method
carries out direct numerical simulations (DNS) of the NavierStokes equations. This capability (called "ShortPipe")
is applicable to short pipes with lengthdiameter ratio L/D = O(110), and is used for understanding detailed turbulent
flow dynamics of roll waves and fully developed slugs and for developing effective turbulent/viscous models for
largerscale applications. The combined framework will be useful for developing effective physicsbased modeling
and strategies for improving fullscale viscous applications such as RANS and phenomenological models for
multiphase roll wave and slug flow prediction in long pipes.
LongPipe simulations are performed to understand the nonlinear effects upon the instabilities of twophase stratified
flows and to study the mechanisms for the development of largeamplitude long interfacial waves by resonant
nonlinear wavewave interactions. Comparisons have been made between LongPipe predictions, existing theoretical
analysis, and laboratory observations on the initial growth and subsequent nonlinear evolution of wavy flows with
good agreement being observed in Campbell (2009). ShortPipe simulations are conducted to understand the detailed
flow dynamics in the later stage of development of breaking waves and slug formation. Of critical interest is the
development and validation of physicsbased turbulence and interface closure models appropriate for use in the
simulations of violent interfacial flow in pipelines.
Nomenclature Re = pjUID Reynolds number
t Time (s)
Roman symbols U Fluid current speed (ms 1)
Fr p Froude number We p U D Weber number
g Gravitational acceleration (ms 2) x Horizontal spatial coordinate (m)
h Equilibrium fluid depth (m) y Vertical spatial coordinate (m)
D Channel depth (m)
1 _Greek symbols
k Wavenumber (m 1) f st
c a Void fraction
N = L Viscosity ratio
N iscs it Pressure forcing coefficient
P Pressure (Nm 2) 6 Dirac delta function
R = Density ratio
P1
ir Free surface elevation (m)
o/0 Initial free surface elevation (m)
7 Surface tension (Nm 1)
SInterface curvature ()
p Viscosity (Nsm 2)
w Frequency (s 1)
T Viscous stress tensor
p Potential function (ms 2)
p Level set function
p Density (kgm3)
Subscripts
u Upper fluid
1 Lower fluid
Introduction
The transition from stratified flow to more complex flow
regimes represents a significant theoretical and engineer
ing challenge. The offshore oil industry's reliance on
long pipelines has raised interest in understanding the
complex flow physics which can arise in the pipeline
during the transport of multiphase mixtures from off
shore to land based processing facilities. Variations in
parameters such as flow rates, pressure drop, pipe incli
nation (due to subsea terrain) can lead to the formation
of liquid slugs within the pipe. The presence of these
slugs create igmiitk.iiI design and flow assurance chal
lenges. Being able to consistently and accurately predict
the conditions which cause these transitions is necessary
for robust design of engineering systems and remains an
active area of research.
The theoretical prediction of the transition to slug flow
has traditionally been addressed through linear stability
analysis of a stratified state. Numerous stability mod
els have been generated such as Taitel & Dukler (1976),
Lin & H.iII.ili\ (1986), Barnea & Taitel (1993), and Fu
nada & Joseph (2001). These models include a wide
range of physical features such as experimental correla
tions, surface tension, normal viscous stresses, and in
terfacial friction factors. A survey by Mata et al. (1993)
found that comparisons between the different stability
models resulted in a wide range of stability predictions
with generally poor comparisons to experimental mea
surements. While linear instability techniques may give
insight into initial instabilities, these results are not valid
once the waves reach finite amplitude states. Being able
to predict the transition from finite amplitude waves to
slug flow remains a difficult challenge which cannot be
solved without robust nonlinear analysis.
Later stages of the transition process are complicated
by the effects of turbulence. Once the waves gen
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
rated by the initial instability reach large amplitude
states, strong viscous and turbulent mechanisms gener
ate significantly more complex flow physics such as bub
ble/spray formation and breaking waves. As the waves
eventually bridge the pipe diameter, the slug body elon
gates and becomes turbulent and bubbly. The presence
of such effects makes it unlikely that accurate predic
tions of slug formation can be made without a detailed
understanding of the turbulent flow physics.
Recently, a study by Valluri et al. (2008) carried
out direct numerical simulations (DNS) of the Navier
Stokes equations to examine the initial evolution of the
interface into slugtype structures in a pressure driven
channel flow. Good agreement was found with the linear
OrrSommerfeld instability predictions and it was found
that large amplitude structures evolve through a non
linear interfacial mechanism. While these simulations
were stopped prior to the wave touching the channel lid,
these results show that DNS will be able to address is
sues which could not be easily examined theoretically or
experimentally.
In addition to slug formation conditions, predicting
the length of the slug body, the pressure drop along the
pipeline, and frequency of slugging is important for en
gineering design. Numerical methods, such as OLGA,
have been developed that can carry out the large scale
simulations needed for modeling flow through pipelines.
These methods are based on phenomenological models,
area averaged techniques or rely on RANS codes. Due
to the length scales involved in the problem, these ap
proaches are necessary in order to be computationally
tractable; however, as a consequence, they are not ca
pable of resolving the important small scaled turbulent
flow physics. More accurate predictions could be made
from these industrial codes if these phenomenological
models could be supplemented with detailed physics
based, multiphase closure models obtained from DNS,
which capture all relevant length and time scales of the
problem or large eddy simulations (LES), which uses
physicsbased closure models to approximate the sub
grid scale physics.
This paper describes two numerical methods which
have been developed for the purpose of determining
some of the underlying physical mechanisms which
cause the transition from a stratified flow to large am
plitude wavy states and the eventual formation of liquid
slugs. One method captures the nonlinear evolution of
instabilities grown from a flat interface using a potential
flow based highorder pseudospectral method. The sec
ond method carries out direct numerical simulations of
the twofluid NavierStokes equations for the purpose of
growing a large amplitude wave into a slug. This code
can illuminate the dominant physical processes involved
in the later stages of slug development, examine the con
tact problem as the wave touches the pipe wall, and pro
vide physical models for the turbulent behavior inside of
the slug body.
Numerical Methods
In order to accurately simulate the development of large
amplitude disturbances and the transition to slug flow, it
is necessary to capture the appropriate physics inherent
to the problem. In order to do this, while still remain
ing computationally tractable, two separate numerical
schemes are described. The first method, referred to as
"LongPipe", examines the evolution of the flow using a
potential flow formulation. The efficiency and accuracy
of this method makes it computationally inexpensive to
examine the initial growth and nonlinear evolution of a
flat interface into large amplitude waves. Once the am
plitude of the disturbances becomes large, viscous (tur
bulent) effects become dominant and are examined by a
secondary method referred to as "ShortPipe". Through
direct numerical simulations, more complicated phe
nomena such as spray formation, wave breaking, and the
early stages of slug flow development can be examined.
While the solution generated by LongPipe does not
satisfy the same conditions as a viscous solution (conti
nuity of velocity and stress at the interface and no slip
at the walls), it is successful in identifying key physics
which should be incorporated into the ShortPipe initial
condition. Following this procedure allows for the ex
amination of realistic large amplitude disturbances with
out the 'igiiii.iiill computational expense of using DNS
to simulate the growth of nonlinear waves from a flat in
terface. In the following sections, the formulation and
capabilities of both of these numerical schemes shall be
discussed.
Together these codes can aide in understanding the
physical processes governing multiphase flow regime
transitions and develop physicsbased models used for
industrial applications.
LongPipe
The formulation of LongPipe assumes that the flow is
composed of two incompressible immiscible fluids in a
horizontal channel. The coordinate system is located at
the equilibrium position between the two fluids with y
directed upwards. The vertical displacement of the in
terface away from the equilibrium position is described
by the function y = n(x, t). The upper and lower fluids
have densities p, and p, and flow with constant uniform
currents Uu and U1 respectively. The initial depths of
each of the fluid layers are hu and hi. Both fluids are
assumed to be irrotational and can be described by po
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
tential functions cp and p which are governed by
V2 0 < y < h
V2p 0
h < y < TK.
No flux conditions are enforced at the channel walls re
quiring
y h.
0 y = h.
Kinematic constraints are enforced between the two flu
ids requiring that the interface remain material for all
time
y = (5)
y = (6)
at ( 5x) ax
The dynamic boundary condition represents a balance
between the fluid pressures, surface tension, and can in
corporate a slope coherent pressure forcing term which
models the effects of shear stress at the interface
R [i + 1P, 2 + U. 0,4+ 1 1
at 2 0x FF,2 9
2 Nj~ 2 09 2 1 r/x
Re [ Y2 y2 We (1i +, ,)3/2
Y = (7)
where B is a dimensionless pressure forcing coefficient,
Re is the Reynolds number, and We is the Weber num
ber.
To solve this system of equations, an efficient high
order spectral method developed by Campbell (2009)
was utilized. This method efficiently tracks the evolution
of a large number of wave components in a broadband
spectrum and accounts for their nonlinear interactions up
to an arbitrary high order of nonlinearity using a pseudo
spectral approach. Validation tests have shown that the
method converges exponentially with interaction order
and the number of wave modes (or grid refinement).
Due to the efficiency of this method, LongPipe is ca
pable of the simulation of flow physics over large scale
(L/D ~ 1000). It can be utilized to generate further de
veloped nonlinear initial conditions for direct numerical
simulations which would significantly reduce the com
putational expense of ShortPipe by not requiring that it
simulate disturbances from a flat interface.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
ShortPipe
ShortPipe carries out direct numerical simulations of
the incompressible, twofluid NaiverStokes equations
for a domain with variable fluid properties. To account
for these effects, this Cartesian Grid method captures
the interface using the level set method, which advects
a signed distance function p (Y) for which p (7) 0
identifies the interface between the two fluids p (7) > 0
in the lower fluid and p (7) < 0 in the upper fluid. This
level set function is governed by
+ u V = 0
where i is the fluid velocity. Knowing implicitly the lo
cation of = 0 then allows for all of the fluid properties
within the domain to be identified by:
p (0) = p./pi + (1 p/lp)H (0) (9)
S( /) + //i + (1 1p/1)H () (10)
where H (y) is a Heaviside function. Applying this
technique then allows for the NaiverStokes equation to
be cast in the form:
D [P()]
Dt
1V V
VP + V (2p (0) 7)
Re
p() 1
p()k+ 1K6 () V (11)
where is the material derivative and 6(0) is a Dirac
delta function.
As in Valluri et al. (2008), the level set implementa
tion is smoothed over a few grid points c and the Heavi
side function H (&) becomes
(1 7 e lower fluid
H(4; e) f(; e) f E interface region .
0 7 e upper fluid
(12)
and f(0; c) is a smoothing function evaluating to 1, 1/2
and 0 for p ec, 0, e, respectively.
f( ; ) = (+ + sin( )) for < .
(13)
For this work, these equations are solved in two di
mensions using a secondorder finite volume method on
a structured grid over a rectangular domain. A third or
der QUICK scheme is used to calculate the advection
terms in the NavierStokes equations. The pressure is
calculated using a Chorin projection method where the
continuity equation is used to project the velocity onto a
divergencefree field. A Vcycle multigrid solver is used
to solve the resulting Poisson equation.
The level set equation is advected using a conservative
form and central differences. While typical implemen
tations use highorder WENO schemes for this equation
(Osher & Fedkiew 2003), experience has shown that this
type of implementation is not necessary if the interface
does not form sharp corners and the additional numer
ical dissipation and volume loss of the WENO scheme
is avoided. Through the course of advection, the level
set function loses its signed distance property and re
quires periodic reinitialization to a signed distance func
tion. The reinitialization stage is performed through a
psuedotime integration of a partial differential equation
as in Sussman & Fatemi (1999). Finally, time integra
tion of the Navier Stokes and level set equations is im
plemented using an explicit second order total variation
dimension (TVD) RungeKutta scheme as proposed in
Trygvasson et al. (2007).
Periodic boundary conditions were imposed along the
length of the domain and both no flux and no slip con
ditions were imposed along the channel walls. A fixed
pressure gradient, with amplitude P(t), along the length
of the channel is applied to enforce desired upper and
lower velocities. P(t) may be constant or a tanh func
tion. For all of the cases presented here, the flow is initi
ated from a smallamplitude two phase Airy wave solu
tion (Hendrickson 2004). If P(t) is constant, a mapped
linear Poisuielle flow solution is added to the Airy wave
solution.
Due to the fact that this code carries out direct numeri
cal simulations of the NavierStokes equations and is re
quired to resolve all small scale physics, this code is lim
ited to relatively small domain lengths (L/D ~ 1 10).
LongPipe Numerical Results
LongPipe was validated through several rigorous tests.
The first validation trial demonstrated convergence of
the method against the analytic solution of the linear
KelvinHelmholtz problem for both stable and unsta
ble waves over a large range of wavenumbers. It was
found that the numerical solutions converged exponen
tially to the analytic solutions with the refinement of the
time step or the grid spacing. This process validated the
numerical scheme for first order expansions.
As a secondary validation of LongPipe, a comparison
was made to a second order nonlinear solution of Eqns.
17 obtained by Campbell (2009). It was found that for
second order nonlinearity, an interface composed of two
specially selected linear modes with wavenumbers and
frequencies (ki, wc) and (k2, w2) could satisfy a triad
resonance condition of the form
k2 k = k3
W2 W1 = 3.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
102
10.
i10
10'
S10'
10' Fo k
200 300 400 500
Nondimensional Time
Figure 1: Comparisons of wave amplitude variations in
a resonant triad between the LongPipe sim
ulation (circle) and theoretical multiplescale
solution (solid line).
with (k3, W3) satisfying the linear dispersion relation
ship. When these conditions are satisfied, a set of non
linear ordinary differential equations of the form
da B23a2a (16)
dt
B13saa3
B12a a2
(where the represents the complex conjugate) are ob
tained by use of the method of multiple scales which de
scribe the resonant interactions. Comparisons between
the solution of these nonlinear interaction equations and
the numerical method validate the accuracy of Long
Pipe for calculating nonlinear solutions. A case with
D 2.54 cm, U, 3.5m, U, (p,.I iI" a 0.8,
R = 1.23 x 103, g 9.81 , ki 27m and
ki2 m1 is shown in Figure 1. In this trial, surface
tension, viscosity, and the slope coherent pressure forc
ing were not included. For this trial, a (0) a 2(0) 
10 5m and a3(0) Om. Good agreement is observed
between the theoretical and numerical solutions.
Traditionally, linear stability analysis is carried out
to determine if it is possible for long waves to be
come unstable; however, it is commonly observed ex
perimentally that flows predicted to be stable to Kelvin
Helmholtz experience wave growth. Through LongPipe
simulations, it was found that long waves could form
due to the transfer of energy from unstable short waves
to stable long waves through nonlinear resonant wave
wave interactions. Consider the case where D 0.5 m,
U, = 10m, U 1m, a 0.3, R 1.23 x 103,
g = 9.81m, ki = 327m k2 = 347m1, al(0)
Figure 2: Time variation of the modal amplitudes of
waves while undergoing resonant interactions
with k2 and 2i _. being unstable to Kelvin
Helmholtz instability.
a2(0) 10 m, a3(0) 0 and a 7.34 x 10 2N
For these flow conditions, all wavenumbers less than
k 347m 1 are stable. From linear theory, long waves
(A ~ O(D)) are predicted to be stable.
If nonlinear resonant interactions are considered, it is
found that k2 satisfies a second harmonic resonance (the
phase velocity of k2 and its second harmonic are equal).
This is unique because 2i _. is linearly unstable and re
sults in the transfer of energy between the wave modes
k2 and 2i _.. Simultaneously, k2 and kl form a triad res
onance with the difference component k = k2 k1.
This leads to the resonant exchange of energy amongst
the modes involved in the triad. These two resonant
mechanisms occur simultaneously and rapid growth of
the longest wave component k3 is found to occur as a
result of this cascaded of energy from the unstable short
waves to the stable long waves. When the interactions
begin, k2 generates its second harmonic 2li due to the
inherent nonlinear interactions of the interface. Since
2i _. falls in the unstable wavenumber range, it begins to
grow exponentially with time. Due to the second har
monic resonance, energy from the linear instability of
the l _2 mode is transferred back to k2. However, since
k2 is simultaneously involved in a triad resonance with
ki and k3, the energy gained by k2 from the second
harmonic resonance is distributed among the resonant
triad. This cascade of energy through the coupled res
onant mechanisms results in rapid growth of the long
mode k3 as can be seen in Figure 2. It can be seen that
all wave modes can grow simultaneously and can reach
amplitudes which are several orders of magnitude larger
than the initial values.
It should be noted that these resonance conditions are
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
1 3 7
x (M)
Figure 3: Large amplitude long wave generated by
spectra of nonlinear resonant interactions with
maximum disturbance amplitude ~ 0.25h,.
not a rare occurrence. For a given spectrum, it is possi
ble for there to be many wave modes which satisfy the
triad resonance condition. When these additional wave
modes are considered, ignikii.lli amounts of energy are
transferred to the long wave components. Consider the
case where the initial interface is composed of all modes
between k = 327rm landk2 347m 1 and all of the
higher harmonics as well as all sum and difference inter
actions are permitted. For this band of wavenumbers,
nearly all components form a triad resonance (or near
resonance) with other components in the spectrum and
participate in the transfer of energy from short to long
waves through a similar mechanism as was previously
shown. Simulations produced by LongPipe show that
the resulting nonlinear evolution of the interface gen
erates large amplitude long waves as shown in Figure
3. This long wave has an amplitude O(0.25h,) and a
wavelength O(2D).
In the last simulation, the numerical method had is
sues due to the growth of very short waves. Due to
the linear instability, the short waves grew quickly and
reached large steepness values. Since the interface was
represented through a Fourier series, it was possible for
the short waves to reach the breaking point and stop
the simulation. To overcome this, a wave breaking fil
ter based on the wave steepness was applied. This fil
ter acted on the short waves by monitoring their modal
steepness (ak). Once a short wave component reached
a particular threshold value, the amplitude of that com
ponent was reset to a small value ~ O(0.la). Several
tests were conducted to determine the effect that this
threshold value had on the global solution and the dif
ferences were found to be iniiigilii,.i Removing the
short waves was not acceptable because the short wave
components generated the unstable energy which was
transferred to the long wave components. Removing the
short waves completely cripples the efficiency of the res
a 20 30 4 0 W M EU W
Figure 4: Standard decay test of twophase Airy wave
using ShortPipe. (green) r(xp, t) (blue) the
oretical decay envelope.
onant mechanism. It is desirable to extend these simu
lations to the case of a broadband spectrum because as
the number of unstable modes increases, so will the am
plitude and growth rate; however, a more robust filtering
technique is needed for this.
ShortPipe Numerical Results
Of particular interest for ShortPipe simulations is en
suring that the dissipation associated with the numerical
scheme is not excessive compared to the viscous dissi
pation. A standard test for single phase viscous flow
solvers is to simulate an Airy wave and compare the
amplitude decay with theoretical values. In a similar
test simulating a twophase Airy wave, ShortPipe was
shown to predict the appropriate viscous dissipation over
a simulation of 36 wave periods. Figure 4 is a wave
probe record with Re=1000, We=oo, Fr=l, pu/pl=O.Ol,
pt,/bi=0.02 and wave amplitude of a 0.1/(27). The
theoretical decay rate is also shown (Hendrickson 2004).
The simulation was performed on a 1282 grid with dt ~
0.006 on 1 node of a dualprocessor, dualcore Athlon
based cluster, utilizing all four cores through a parallel
MPI implementation, in ~ 15 minutes for 15,000 time
steps.
As the purpose of ShortPipe is to understand the de
tailed flow dynamics in the later stages to develop and
validate physicsbased turbulence models, the simula
tions of multiphase pipe flows over a range of operating
parameters is critical. The DNS datasets generated as a
part of this effort will be used for a priori analysis to de
velop and validate physicsbased turbulence models for
use in larger scale simulations such as Large Eddy Sim
ulation (LES) and RANS. To date, ShortPipe has been
utilized to generate over a dozen datasets within two
classes of runs. The first class of runs, called Case A, are
I 1 1 I I I I I
0l01
o 005
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
100
80
60
40
20
20
40
60
80
100
(a)
200
160
120
80
40
40
80
120
160
200
Figure 5: Case A simulation of a twophase airwater
Airy wave with a pressure gradient by Short
Pipe. U,,=0.62 m/s, U,1=0.015 m/s, and
a=0.5. Re=2000, Fr=l, We=oo. In (a), black
line is interface between two fluids and chan
nel wall.
for air and water at a Reynolds number of 2000, which
equates to a channel height of only 1.25 cm. While
ShortPipe is capable of including the effects of surface
tension, none are currently included in these classes.
For Case A, Us,/Us8 ranged from 20100 resulting in
mainly entraining events. Figure 5 is a simulation of
a twophase Airy wave with wavelength equal to the
pipe diameter that has the pipe pressure increased from
zero over a period of time resulting in Us,=0.62 m/s and
U,1=0.015 m/s. Figure 5a shows the transverse vortic
ity contours while Figure 5b shows the velocity profile,
where select vertical slices are shown for clarity. Flow is
from left to right. Separation off of the crests impacts the
back face of the downstream wave. There is evidence of
an entrainment event as the surface begins to pinch off
at the crest of the wave. The resolution for this case
is 512x256 and is not sufficient to resolve the interaction
of the crest separation vortex impacting the interface and
thus the simulation stops.
The second class of runs, called Case B, are for
highviscosity oil (0.181 PAs) in a 5 cm channel with
Reynolds number 0(200) (Gokcal et al. 2006). These
cases have Us,/U,1 of 0.110 and result in mainly elon
gated bubble events, where the interface attaches to the
channel wall and produces little to no entrainment. Fig
ure 6 is a simulation of a small amplitude twophase Airy
S11100
60
40
S20
20
640
100
(b)
Figure 6: Case B simulation of a twophase Airy wave
of air with highviscosity oil with a pres
sure gradient by ShortPipe. U,,=0.09 m/s,
Us,=0.8 m/s, and a=0.9. Re=200, Fr=l,
We=oo. In (a), black line is interface between
two fluids and channel wall.
wave (ka 0.01) with a mapped linearized Poisuielle
flow added to the velocity profile. The wave is forced
via slopecoherent pressure forcing until it touches the
top of the channel, forming an air cavity. Shown in Fig
ure 6a is the transverse vorticity contours and in Figure
6b is the velocity field. Again, select vertical slices are
shown for clarity. The boundary layers for the walls and
interface are clearly visible as is the formation of a jet
inside the air cavity.
Conclusions
Two complementary computational capabilities, Long
Pipe and ShortPipe, form a framework that is used
to understand the transition and nonlinear evolution of
multiphase pipe flows in support of the design of deep
water pipelines. LongPipe is used to understand the
initial growth and nonlinear evolution of a flat interface
into large amplitude waves and ShortPipe is used to un
derstand complex interfacial phenomena such as wave
breaking, slug flow and entrainment. These two capa
bilities will be used in tandem to understand the physi
cal processes governing the transition of multiphase flow
regimes and develop physicsbased models for industrial
applications.
Preliminary studies using LongPipe have demon
200
160
120
80
40
40
80
120
160
200
mD
strated it's abilities to develop large amplitude long
waves from a "nearly" flat interface using the funda
mental physics of the system. A nonlinear resonant in
teraction theory was found which transfers the energy
generated by linearly unstable waves across the wave
spectrum. This cascade of energy focuses energy among
the long wave modes and generates large amplitude long
waves which could possibly lead to the development of
large amplitude roll waves or slugs. The nearterm goal
of LongPipe is to improve flow transitions criteria by
including nonlinear resonant wave interaction effects.
A stateoftheart multiphase CFD solver, ShortPipe,
has been developed and validated for use in DNS of
late stage flow dynamics in pipes. Its ability to sim
ulate a range of operating parameters (fluid constitu
tive properties, Reynolds numbers, U,,/Us, and chan
nel height) has been established. Future applications of
ShortPipe will be to further study simulations within the
two classes of waves presented at improved resolutions:
(i) larger Reynolds number airwater flows with large
U,g and (ii) lower Reynolds number airhighviscosity
oil flows where slugs and elongated bubbles are known
to form. A third class of simulations are currently be
ing investigated for SF6water flows at Reynolds num
ber 0(2500), as experimental results report that simi
lar operating conditions form roll waves with igmniik.ii.i
gas entrainment in the water phase. These three classes
of simulations will be used in a priori development and
validation of physicsbased turbulence closure modeling
for largescale applications.
Acknowledgements
This research was sponsored by Chevron Corporation.
References
Barnea, D. and Taitel, Y. 1993. KelvinHelmholtz sta
bility criteria for stratified flow: viscous versis non
viscous(inviscid) approaches. Int. J. Multiphase Flow
19, 639649.
Campbell, B. 2009. Nonlinear Effects on Interfacial
Wave Growth into Slug Flow. M.S. Thesis. Mas
sachusetts Institute of Technology.
Funada, T. and Joseph, D.D. 2001. Viscous potential
flow analysis of KelvinHelmholtz instability in a chan
nel. J. Fluid Mech. 445, 263283.
Gokcal, B., Wang, Q., Zhang, H.Q., and Sarica, C.
2006. Effects of High Oil Viscosity on Oil/Gas Flow Be
havior in Horizontal Pipes. Presented at the 2006 SPE
Annual Technical Conference and Exhibition, San An
tonio, TX, USA.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Hendrickson, K. 2005. NavierStokes Simulation of
Steep Breaking Water Waves with a Coupled AirWater
Interface. Ph. D. Thesis. Massachusetts Institute of
Technology.
Lin, P.Y. and H.iiii.ii\ T.J. 1986. Prediction of the initi
ation of slugs with linear stability criterion. Int. J. Mul
tiphase Flow 12, 7998.
Mata, C., Pereyra, E., Trallero, J.L., and Joseph, D. D.
2002. Stability of stratified gasliquid flows. Int. J. Mul
tiphase Flow 28, 12491268.
Osher, S. and Fedkiew, R. 2003. Level Set Methods and
Dynamic Implicit Surfaces. Applied Maethematical Sci
ences vol. 153. Springer.
Sussman, M. and Fatemi, E. 1999. An efficient,
interfacepreserving level set redistancing algorithm and
its application to interfacial incompressible fluid flow.
SIAM J. Sci. Comput. 20(4), 11651191.
Taitel,Y. and Dukler, A.E. 1976. A model for predicting
flow regime transitions in horizontal and near horizontal
gasliquid flow. AIChE J 22, 4755.
Tryggvason, G., Bunner, B., Esmaeeli, A., Juric, D., Al
Rawahi, N., Tauber, W., Han, J., Nas, S., and Jan, Y. J.
2001. A FrontTracking Method for the Computations
of Multiphase Flow. J. Comp. Phys. 169:2, 708 759.
Valluri, P., Spelt, P.D.M., Lawrence, C.J., and Hewitt,
G.F. 2008. Numerical simulation of the onset of slug ini
tiation in laminar horizontal channel flow. Int. J. Multi
phase Flow 34, 206225.
