Fourth International Conference on CFD in the Oil and Gas, Metallurgical & Process Industries

SINTEF / NTNU Trondheim, Norway

68 June 2005
PREDICTION OF TWO-PHASE PIPE FLOWS USING
SIMPLE CLOSURE RELATIONS IN A 2D TWO-FLUID MODEL
Svend Tollak Munkejord
1
Mona J. Mlnvik
2
Jens A. Melheim
1
Inge R. Gran
2
Robert Olsen
1
1
Norwegian University of Science and Technology (NTNU),
Department of Energy and Process Engineering, Kolbjrn Hejes veg 1A, NO-7491 Trondheim, Norway
2
SINTEF Energy Research, Energy Processes, Kolbjrn Hejes veg 1A, NO-7465 Trondheim, Norway
ABSTRACT
This paper presents a comparison of pressure drop and
liquid hold-up as calculated using different modelling
strategies; a two-dimensional two-uid model, a one-
dimensional simulator, and engineering correlations. The
numerical results are compared with experimental data. In-
terestingly, the two-dimensional two-uid model performed
well, using simple constitutive relations being functions of
the local ow conditions.
In the present work, a two-uid model has been imple-
mented in the framework of a two-dimensional multiphase
Computational Fluid Dynamics (CFD) code. The govern-
ing equations were spatially discretized using the nite-
volume technique, and the time-integrator was an explicit
low-storage ve-step fourth-order Runge-Kutta scheme.
Numerical results from the present program, the OLGAS
simulator, the Beggs and Brill correlation, and the Friedel
and the Premoli correlations have been compared to 52 data
points from the TILDA two-phase pipe ow database at the
SINTEF multiphase ow laboratory. The data were taken in
a 0.189m inner-diameter pipe. The liquid volumetric ux
j

varied between 0.1 and 1.0m/s, the gas volumetric ux

j
g
was in the range from 0.5 to 12m/s, and the pressure
p ranged from 20 to 90bar. Inclination angles of 0 and 1
were used, with the majority of the data coming from the
horizontal pipe.
We nd it interesting to note the relatively good corres-
pondence between experimental data and the computed res-
ults of the CFD program, particularly when considering the
simple constitutive relations employed.
Keywords: multiphase ow, pipe ow, computational
uid dynamics (CFD), closure relations, pressure drop, li-
quid hold-up, explicit scheme
NOMENCLATURE
Latin letters
A Area m
2
A Control surface m
2
B
k
Displacement factor
C Friction parameter, see equation (14) m
1
C

, C
1
, C
2
Constants in the k- turbulence model
c Speed of sound m/s
d Diameter m

Corresponding author. E-mail: stm (a) pvv.ntnu.no

E Constant in wall function
F Force per unit volume N/m
3
f Body force eld m/s
2
h Channel height m
I Unitary tensor
j Volumetric ux (supercial velocity) m/s
k Turbulence energy m
2
/s
2
N Number of data points
n Normal distance from wall m
n Unit normal vector m
P
k
Production rate of turbulence energy m
2
/s
3
p Pressure Pa
p Local pressure variation, see equation (6) Pa
Re Reynolds number
t Time s
T
k,w
Wall-function variable, see equation (29) Pa s/m
u Velocity vector m/s
V Volume m
3
V Control volume m
3
x Length coordinate m
y
+
Non-dimensional distance from wall in turbulent
shear layer
y
+
0
Constant in wall function
Greek letters
Volume fraction

Liquid hold-up
Dissipation rate of turbulence energy m
2
/s
3
Friction factor between phases k and l, see equa-
tion (13)
Mass source kg/(m
3
s)
Constant in wall function
Molecular viscosity Pa s

T
Eddy (turbulent) viscosity Pa s
Density kg/m
3
Turbulent Prandtl/Schmidt number

w
Wall shear stress Pa
Stress tensor Pa
General variable
Subscripts
d Drag
1
g Gas
h Homogeneous (no-slip)
i Direction i in Cartesian coordinates
j Direction j in Cartesian coordinates
k Pertaining to phase k
l Pertaining to phase l
Liquid
m Mixture
P Near-wall point, see equation (28)
par Parallel to wall, see equation (28)
w Wall
Supercripts
Reference
T Turbulent
1 INTRODUCTION
1.1 Background
A major challenge faced by the oil and gas industry stems
from the tendency towards elds on deeper water situated at
greater distance from the shore. This increases the drive to-
wards eld developments based on sub-sea processing and
multiphase ow transportation. The older two-phase ow
installations have gone off plateau and many of them now
produce signicant amounts of water. Water introduces new
challenges related to ow assurance, exemplied by gas hy-
drates and scale deposits. Sub-sea processing, increased
transport distance, and higher pipeline inclination in hilly
terrain require qualied multiphase ow simulation tools.
Currently, one-dimensional simulation tools are routinely
employed by engineers for designing and operating mul-
tiphase ow of oil and gas in long pipelines. OLGA is a
dynamic modied two-uid model developed primarily to
simulate slow transients in two- and three-phase pipeline
ow. The rst version of OLGA was developed in 1979 and
the model has been further developed by the Institute for
Energy Technology (IFE) and Scandpower since then. The
model uses two momentum equations, one representing the
gas and possible liquid droplets and one representing the re-
maining liquid ow. The equations are solved numerically
by a semi-implicit time-marching scheme.
One-dimensional models such as OLGA have enjoyed
considerable success in predicting the intrinsically one-
dimensional situation encountered in long pipelines. How-
ever, the above-mentioned problems exhibit two- and three-
dimensional ows, e.g. in short pipes, bends and junctions,
that are not suitable for one-dimensional approaches.
In the present work, a two-uid model has been imple-
mented in the framework of a two-dimensional multiphase
Computational Fluid Dynamics (CFD) code. It differs from
most other CFD simulators in that an explicit Runge-Kutta
scheme is employed for advancing the solution in time. The
presented results show that the multiphase equations are
solved robustly and accurately. The method is suitable for
efcient execution on massively parallel computers.
In this paper, we study two-phase high-pressure pipe
ow. This is an important case since benchmark experi-
mental data exist and because the performance of the simu-
lator can be compared directly with dedicated pipeline mod-
els. Nevertheless, the motivation for the development of
a CFD simulator is to predict the ow in multidimensional
geometries, such as valves, bends, separators, etc.
1.2 Previous work
Early numerical models of multiphase ows were most of-
ten tailored to a particular geometry, and specic ow re-
gimes. Examples of such models are OLGAS (Bendiksen
et al., 1991) and various ones used in the nuclear industry.
Improved and specialized models continue to be presen-
ted. Lin and Lopez (1997) made a two-uid model of
wavy separated two-phase ow, in which the shape of the
gasliquid interface was assumed to be known a priori, and
where the momentum transfer was calculated using wave
theory. Another approach was used by Newton and Behnia
(2001), who solved the steady axial (one-dimensional) mo-
mentum equation together with a low-Reynolds-number k-
turbulence model. A bipolar grid conveniently tted the
geometry of the stratied ow.
The use of two and three-dimensional CFD for multiphase
ows can facilitate the solution of ows in complex geo-
metries, as well as ow phenomena which otherwise could
not be calculated. Early two-uid ow codes primarily
originate from single-uid ow codes based on the Impli-
cit, Continuous-uid, Eulerian (ICE) solution scheme de-
veloped by Harlow and Amsden (1971) or the Semi-Implicit
Method for Pressure-Linked Equations (SIMPLE) by Pa-
tankar and Spalding (1972). Examples of early two-uid
solution schemes are variants of the Implicit MultiField
method (IMF) (Amsden and Harlow, 1974, 1978; Rivard
and Torrey, 1977), the Inter-Phase Slip-Algorithm (IPSA)
(Pun et al., 1979; Spalding, 1980), Gas And Liquid Ana-
lyzer (GALA) (Spalding, 1977), and TRAC (Liles and Reed,
1978). Such semi-implicit methods are in use in several
commercial CFD codes.
Moura and Rezkallah (1993) studied the two-phase ow
distribution in a T-junction, and reported good agreement
between calculated results and experimental data for the
phase separation. Alajbegovi c et al. (1999) calculated
particle/liquid pipe ows, and the radial position of the
particles, and the axial and uctuating velocities, compared
fairly well with experiments.
Saurel and LeMetayer (2001) developed a multiphase
model for compressible ows with interfaces, shocks, det-
onation waves and cavitation. This shows the potential of
CFD models to solve a range of different problems. Al-
though commercially available CFD codes are able to suc-
cessfully predict some multiphase ows (Brown, 2002),
several difcult problems persist, including such topics as
accuracy and convergence.
Little work concerning CFD calculations of pressure drop
and liquid hold-up of multiphase ow in horizontal pipes
has been published. However, since relevant high-pressure
experimental data for such setups exist, this is an interesting
benchmark case.
1.3 Outline of paper
In the following section, our mathematical model will be
presented in detail. Next, the numerical solution method is
briey described. Then the setup of the channel ow case
study is explained, and the results of the present model are
compared to those of other models and to experimental data.
Finally, conclusions are drawn.
2
2 MULTI-FLUID MODEL
The multi-uid equations can be deduced from the basic
equations of single-phase ow. Several approaches have
been used, for instance temporal averaging, volume aver-
aging, and ensemble averaging (see e.g. Ishii, 1975; Soo,
1989; Drewand Passman, 1999). It is also possible to derive
the multi-uid equations without averaging, but by intro-
ducing distributions, i.e., generalized functions (Kataoka,
1986). In the present work, the volume-averaging ap-
proach presented by Soo (1989, Chapter 2) and Soo (1990,
Chapter 6) has been employed. The different deductions
lead to equations having the same form. However, the form
of the unclosed terms and the interpretation of the variables
in the equations may differ.
2.1 Basic equations
The volume-averaged continuity equation is

t

k
+
k
u
k
=
1
V
_
A
k

k
(u
k
u
s
) n
k
dA =
k
,
(1)
and the volume-averaged momentum equation may be writ-
ten as

t

k
u
k
+
k
u
k
u
k
= p
k
+
k
+
k
f
+
1
V
_
A
k
(p
k
I +
k
) n
k
dA

1
V
_
A
k

k
u
k
(u
k
u
s
) n
k
dA. (2)
Herein subscript k refers to phase k; hence Einsteins sum-
mation rule does not apply. Further, the volume average of
a variable is dened by

k
=
1
V
_
V
k

k
dV, (3)
and the intrinsic volume average is
i

k
=
1
V
k
_
V
k

k
dV. (4)
That is,
k
is averaged over the whole control volume V,
whereas
i

k
is averaged only over the part of the control
volume where phase k is present, V
k
.
The integrals in Equations (1) and (2) appear due to the
averaging. The average of time derivatives is rewritten us-
ing the Reynolds transport theorem, whereas the average of
divergence and gradient terms is rewritten using the Slattery
averaging theorem (Slattery, 1967; Whitaker, 1969). The
resulting transfer integrals describe phase interactions due
to pressure, viscous stresses, and inertia forces.
In the present work, the energy equation has not been
considered. Further, the experimental data used in the com-
parison were adiabatic. It has been assumed that no mass
transfer between the phases takes place:

k
=0, (5)
which implies that the integrand of the second transfer in-
tegral in Equation (2) is zero. However, it is straightforward
to include other expressions for mass transfer in our numer-
ical model.
The transfer integral due to pressure and viscous stresses
needs to be modelled. However, the pressure p
k
on the in-
terface is generally unknown. It is therefore convenient to
assume that it may be written as the sum of the mean pres-
sure in the control volume, and a local variation:
p
k
=
i
p
k
+p
k
. (6)
A similar splitting was made by Drew and Passman (1999,
Section 11.3.2).
The following relation between the gradient of the
volume fraction
k
= V
k
/V and the surface integral of the
unit normal vector can be found, using the Slattery aver-
aging theorem:

k
=
1
V
_
A
k
n
k
dA. (7)
Insert Equations (6) and (7) into the pressure-part of the
transfer integral of Equation (2):
1
V
_
A
k
p
k
I n
k
dA =
i
p
k

k

1
V
_
A
k
p
k
n
k
dA. (8)
For the pressure-gradient term of Equation (2), we have:
p
k
=
k

i
p
k

i
p
k

k
, (9)
and the
i
p
k

k
term from the pressure gradient term and
from the transfer integral are seen to cancel. Therefore,
the volume-averaged momentum equation without phase
change may be written as

t

k
u
k
+
k
u
k
u
k
=
k

i
p
k
+
k

+
k
f +
1
V
_
A
k
(p
k
I +
k
) n
k
dA. (10)
2.2 Equation of state
The following equation of state gave the relation between
pressure and density:
p
k
=c
2
k
(
k

k
), (11)
where the speed of sound c
k
and the reference density

k
were constants.
2.3 Interface forces
The integral in Equation (10) was modelled the following
way:
1
V
_
A
k
(p
k
I +
k
) n
k
dA =(1B
k
)
i
p
k

k
+F
d
,
(12)
where B
k
is a displacement factor close to unity, which
can be regarded as a simplied model for forces causing
dispersion of the volume fraction prole, e.g. in intermittent
ow. See Appendix A.
F
d
is the drag force per unit volume, for which the fol-
lowing model was assumed (see e.g. Moura and Rezkallah,
1993):
F
d
=

l=k

kl
u
l
u
k
(u
l
u
k
), (13)
where

kl
=C
m

l
, k =l and [C] =m
1
, (14)
3
and where
u
2
u u. (15)
Herein,
m
is the mixture density:

m
=

k
. (16)
The drag force F
d
on phase k is due to the presence
of phase l, and it may be observed that the drag force is
proportional to the square of the relative velocity between
the phases. Further, the factor
kl
is proportional to the
mixture density and the area density (interfacial area per
unit volume) between the phases, and it approaches zero
when one of the volume fractions vanishes. C is a model
constant with dimension 1/L, where the appropriate length
scale depends upon the ow regime, among other things.
Due to surface tension and other effects, the pressure may
in general differ between the phases. In the present work,
however, the pressure was taken to be the same in all phases.
2.4 Turbulence model
In the present work, turbulence was modelled by introdu-
cing a spatial uctuation in the velocity (see e.g. Nigmat-
ulin, 1991, page 35):
u
k
=u
k
+u

k
. (17)
u
k
is the velocity of phase k at a point inside the control
volume V at time t , whileu
k
is its mean, and u

k
is its spatial
uctuation. The mean velocity is dened by:
u
k


k
u
k

, (18)
which is mass-weighted or Favre averaged. Hence, the
volume average of the spatial uctuation is zero:
_

k
u

k
_
=0. (19)
Volume averaging can only be applied to quantities per
volume or area. Therefore, when for instance velocity is
averaged, it is necessary to include the density:
i
u
k
=
1

k
V
_
V
k

k
u
k
dV =

k
u
k

. (20)
It follows from Equation (18) and Equation (20) that u
k
=
i
u
k
when the Favre averaging is based on volume aver-
aging. Therefore, the rules for Favre averaging can be ap-
plied to the volume-averaged equations.
The resulting Favre-averaged continuity equation is,
when Equation (5) has been taken into account:

t

k
+
k

i
u
k
=0, (21)
and for the momentum Equation (10) we obtain

t

k

i
u
k
+
k

i
u
k

i
u
k
=
k

i
p
k

(1B
k
)
i
p
k

k
+
_

k
+
T
k
_
+
k
f +F
d
,
(22)
after using Equation (12). Here,
T
k
is the uctuating stress
tensor:

T
k
=
_

k
u

k
u

k
_
, (23)
Table 1: Constants in the standard k- model.
C

C
1
C
2

k

0.09 1.44 1.92 1.0 1.3

which in the present work was modelled by the Reynolds
stress tensor.
In the present work, a two-equation turbulence model
was used. The single-phase k- model (Launder and Spald-
ing, 1974) was extended to multiphase ows by appropri-
ately accounting for the volume fraction of the phases. The
equation for the turbulence energy k was as follows, written
in Cartesian tensor notation:

k
k
k

t
+

k
k
k

i
_
u
k, j
_
x
j
=

x
j
__

k
+

T
k

k
_

i
k
k

x
j
_
+P
k
k

k

k
(24)
and the equation for the dissipation rate of turbulence en-
ergy was:

k

t
+

k

i
_
u
k, j
_
x
j
=

x
j
__

k
+

T
k

x
j
_
+C
1
i

i
k
k

P
k
k
C
2
i

i
k
k

k
, (25)
where the production P
k
k
of turbulence energy is given by:
P
k
k
=
_

T
k
_

i
_
u
k,i
_
x
j
+

i
_
u
k, j
_
x
i
_

2
3

i j
_

T
k

i
_
u
k,l
_
x
l
+
k
k
k

__

i
_
u
k,i
_
x
j
. (26)
The model for the eddy viscosity was:

T
k
=C

i
k
k

2
i

. (27)
The constants are equal to Launder and Spalding (1974)s
constants, given in Table 1. Multiphase effects and buoy-
ancy effects that might affect the production and dissipation
of turbulence energy, have not been considered in the above
equations. This is a simplication. Multiphase effects are,
however, accounted for to the extent that when the velocity
gradients change, so does the production term (26) of the
turbulence energy transport Equation (24).
2.5 Wall boundary
For the wall boundaries, a standard wall function for smooth
walls was used (see Schlichting, 1979, page 602). It was
extended to multiphase ow by accounting for the volume
fraction of the phases, but disregarding other possible mul-
tiphase effects.
The wall shear stress was represented by the following
expression:

k,w
=
k
T
k,w
_
i
u
k

P,par
u
wall
_
, (28)
4
Table 2: The constants in the wall function.
E y
+
0
0.4187 9.0 11.63
where u
wall
is the velocity of the wall, zero in the present
study, and
i
u
k

P,par
is the velocity vector parallel to the
wall. T
k,w
is given by the wall function:
T
k,w
=
_

k
/n
P
y
+
k,P
y
+
0
,
(
i

k

k,1
)
1/2
/ln
_
Ey
+
k,P
_
y
+
k,P
> y
+
0
,
(29)
where

k,1
=C
1/2

i
_
k
k,P
_
, (30)
and
y
+
k,P
=
_
i

k

k,1
_
1/2
n
P

k
. (31)
In Equations (28)(31), subscript P refers to the near-
wall point in the computational domain. n
P
is the normal
distance from point P to the wall.
k,1
is the shear stress
in the internal sub-layer and is assumed to be equal to
k,w
,
and y
+
k,P
is a dimensionless wall coordinate. The wall func-
tion constant for smooth walls are given in Table 2.
2.6 Numerics
Here we will briey describe the numerical solution
method. The equations were discretized using the nite-
volume method on a curvilinear non-orthogonal colocated
grid. The simulation program could use the power-law as
well as the second-order upwind discretisation scheme (Me-
laaen, 1992). In the present case, the differences between
the two were small.
The result of the application of discrete approximations
for the spatial terms is that the governing equations consti-
tute a set of coupled ordinary differential equations (ODEs).
The remaining derivatives are then discretized in time by
time steps.
Families of low-storage, explicit Runge-Kutta schemes
have been derived and successfully applied for integrating
the compressible Navier-Stokes equations (Kennedy et al.,
2000; Williamson, 1980). In the present study, the temporal
integration of the semi-discretized governing equations was
performed using the ve-step fourth-order scheme by Car-
penter and Kennedy (1994).
It is the Courant-Friedrichs-Lewy (CFL) number that de-
termines the maximum time-step length, and here, the fol-
lowing estimate was used:
CFL =(|u| +c)
t
x
, (32)
where c is the speed of sound, t is the time-step length
and x is the grid spacing in the relevant direction. All of
the present computations were run with CFL = 1, and the
solutions were veried to be independent of the time-step
length.
To avoid checkerboard oscillations, a fourth-order ar-
ticial dissipation lter was applied to the continuity and
momentum equations. It used the Rhie and Chow (1983)
Table 3: The range of experimental data employed in the
present study.
Quantity Range Unit
Liquid volumetric ux, j

0.11.0 m/s
Gas volumetric ux, j
g
0.512 m/s
Pressure, p 2090 bar
Pipe inclination angle, 01
momentum interpolation approach as a starting point. See
also for instance Ferziger and Peri c (1999, Section 7.5.2).
Instead of using a pressure-based lter, which is sufcient
in incompressible ow, the present lter is based on for
the continuity equation and u for the momentum equa-
tions.
The present numerical method has been tested and the
solutions have been found to converge by grid renement.
3 CASE STUDY: CHANNEL FLOW
3.1 Introduction
One of the main points of CFD is the calculation of ows in
complex geometries. However, experimental data for such
geometries are scarce. The objective of the present study
has therefore been to compare CFD calculations with exper-
imental data from full-scale tests with real conditions. This
was regarded as an important step in the validation of the
CFD code.
The calculated results were compared to data from the
SINTEF two-phase pipe ow database TILDA, which is used
in the oil and gas industry. The data were taken in a pipe
with an inner diameter d
i
of 0.189m (8

pipe) with two-

phase ow of naphtha and nitrogen. The range of exper-
imental data employed in the present study is shown in
Table 3. Most of the data were horizontal. The ow re-
gime was reported to be outside the slug range. The data
are for high pressures, where there is a strong gas-liquid in-
teraction. Hence, the slip between the phases is low.
The CFD calculations were performed two-
dimensionally, with the hydraulic diameter corresponding
to the inner diameter of the pipe, that is, with a channel
height of h = 1/2d
i
. Idelchik (1994) reports that for
turbulent ow, the pressure drop in a two-dimensional
channel should be expected to be 10% higher than that
of a channel with a square cross-section. However, for
calculations of single-phase ow, this could neither be
conrmed using the present computer program, nor with
the commercial code Fluent 6.0. Therefore, in the present
work, the calculated pressure drops were directly compared
to the experimental data.
The experimental test loop was modelled as a 200m long
straight channel. With such a length, the ow was fully de-
veloped at the outlet. At the inlet, gas and liquid mass uxes
were specied, and plug velocity proles were assumed. At
the outlet, pressure was specied. Further, the liquid hold-
up (mean liquid volume fraction) was calculated as

=
1
h
_
h
0

dy, (33)
while the gas and liquid volumetric uxes (also referred to
as supercial velocities in the literature) were calculated
5
Table 4: Model parameters in the interface force models.
Quantity Value Unit
Displacement factor, B
k
0.99992
Friction parameter, C 30 m
1
as
j
k
=
1
h
_
h
0

k
u
k
dy. (34)
The inlet gas and liquid mass uxes were varied in such
a way that the outlet gas and liquid volumetric uxes co-
incided with the values used in the experiments. Thus the
calculated pressure gradient and liquid hold-up at the outlet
could be compared to the experimental values.
For the present calculations, an orthogonal, Cartesian
equidistant computational grid of 22 points (20 control
volumes) in each coordinate direction was used. Typically,
when the number of control volumes was doubled, the pre-
dicted pressure gradient and liquid hold-up changed in the
order of 1%. This was considered to be less than the exper-
imental uncertainty, and also less than the uncertainty in the
currently employed constitutive relations.
It was necessary to assume values for the model para-
meters introduced in the interface force models. The dis-
placement factor B
k
(see Equation (12)) affects the volume
fraction proles. For all the calculations, it was set equal
to B
k
= 0.99992, see Table 4, which gave a relatively
smooth transition between the liquid-dominated and the
gas-dominated zones. To satisfy Newtons third law, B
k
was
set equal in both phases. The interface friction parameter C
(see Equation (14)) was assigned a value of C = 30m
1
.
This implied a length scale in the order of centimetres,
which seemed reasonable in the present case.
C and B
k
were kept constant in all the calculations, i.e.,
different ow regimes were not explicitly accounted for.
Importantly, no attempt was made to t the calculated val-
ues to the experimental data by adjusting the parameters.
3.2 Calculated results and comparison with
experimental data
In this section, the results obtained with the present program
are compared with experimental data, and also with results
from the 1D simulator OLGAS, and from engineering cor-
relations.
In the present study, results from a total of N = 52 data
points were evaluated. In addition to comparing our calcu-
lations with experimental data, we compared them to results
obtained by using the stationary oil and gas pipeline simu-
lator OLGAS 2000 v. 2.00. It is an engineering tool being
in use in the oil and gas industry. OLGAS 2000 calculates
the two-phase ow one-dimensionally, and its models have
been t to experimental data. OLGAS is a stationary version
of the OLGA model (Bendiksen et al., 1991).
Furthermore, our calculations were compared to two en-
gineering correlations; the Beggs and Brill (1973) correl-
ation as presented in Brill and Mukherjee (1999) (without
the Payne et al. modications), and a combination of the
Friedel (1979) correlation for the frictional pressure drop
and the Premoli et al. (1971) correlation for the volume
fraction (to determine the gravitational pressure drop and
the liquid hold-up). The former correlation was developed
specically for two-phase ow in oil and gas pipelines,
gas volumetric flux (m/s)
l
i
q
u
i
d
h
o
l
d
-
u
p
(

)
0 5 10
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8 present
OLGAS
Beggs & Brill
Premoli
Homogeneous
Figure 1: Calculated liquid hold-up

as a function of gas
volumetric ux j
g
, for a pressure p = 45bar and
a liquid volumetric ux j

=1m/s.
whereas the latter was made for air-water ow. These cor-
relations are relatively well-known, and will not be repeated
here.
3.2.1 Liquid hold-up vs. gas volumetric ux
Figure 1 shows liquid hold-up as a function of gas volumet-
ric ux for a pressure of p =45bar and a liquid volumetric
ux of j

=1m/s in a horizontal pipe. The results from the

present calculations are plotted together with data obtained
with the OLGAS program, the Premoli correlation, and the
Beggs and Brill correlation. The liquid hold-up calculated
for homogeneous ow (equal gas and liquid velocities, or
no-slip) is also shown in the gure. In the usual situation
that the liquid does not ow faster than the gas, the liquid
hold-up cannot attain lower values than those of homogen-
eous ow.
All the calculation methods captured the generally fall-
ing trend of the liquid hold-up

for increasing gas volu-

metric uxes j
g
. The present program performed well, al-
though with a slight underprediction for low j
g
and a slight
overprediction for high j
g
. The Premoli correlation repro-
duced the experimental data well throughout, except, curi-
ously, for j
g
= 4m/s. The OLGAS program predicted the
liquid hold-up slightly too high for j
g
= 11m/s, but with
an increasing overprediction for decreasing gas volumetric
uxes. The Beggs and Brill correlation performed oppos-
itely, with a slight underprediction for low gas volumetric
uxes, but with a substantial overprediction for high ones.
A similar plot, but for a pressure of 20bar and a li-
quid volumetric ux of 0.1m/s, is given in Figure 2. For
this case, the present model still predicts the high-gas-
volumetric-ux data well, but it signicantly underpredicts
the liquid hold-up for low gas volumetric uxes. The same
can be said for the Premoli correlation, although with a
slightly less pronounced underprediction for the low gas
uxes. The OLGAS program and the Beggs and Brill cor-
relation performed similarly to what they did for the p =
45bar, j

=1m/s case.
3.2.2 Pressure drop vs. gas volumetric ux
Figure 3 compares calculated pressure drops
p
x
for a pres-
sure of p =45bar and a liquid volumetric ux of j

=1m/s
6
gas volumetric flux (m/s)
l
i
q
u
i
d
h
o
l
d
-
u
p
(

)
0 5 10
0
0.1
0.2
0.3
0.4
0.5 present
OLGAS
Beggs & Brill
Premoli
Homogeneous
Figure 2: Calculated liquid hold-up

as a function of gas
volumetric ux j
g
, for a pressure p = 20bar and
a liquid volumetric ux j

=0.1m/s.
gas volumetric flux (m/s)

p
r
e
s
s
u
r
e
g
r
a
d
i
e
n
t
(
P
a
/
m
)
0 5 10
0
100
200
300
400
500
600
700
present
OLGAS
Beggs & Brill
Friedel
Figure 3: Calculated pressure drop
p
x
as a function of gas
volumetric ux j
g
, for a pressure p = 45bar and
a liquid volumetric ux j

=1m/s.
in a horizontal pipe. In this case, all calculation methods re-
produced the experimental values reasonably well for the
low gas volumetric uxes. For increasing gas volumet-
ric uxes, the present program slightly overpredicted the
pressure drop, whereas the OLGAS program slightly under-
predicted it. This tendency of underprediction was more
distinct for the Beggs and Brill correlation, and even more
so for the Friedel correlation.
The case of a pressure of p = 20bar and a liquid volu-
metric ux of j

= 0.1m/s is given in Figure 4. For the

lowest gas volumetric uxes j
g
, the present program and the
Friedel correlation overpredicted the pressure drop, whereas
the OLGAS program and the Beggs and Brill correlation
matched the experimental data quite well. For the highest
gas volumetric ux, the present program and the Beggs and
Brill correlation signicantly underpredicted the pressure
drop. The Friedel correlation underpredicted slightly less,
while the pressure drop calculated by OLGAS only was a bit
too high.
We hypothesize that the tendency of the present model
increasingly to underpredict the pressure drop for increas-
ing gas volumetric uxes at low pressures and low liquid
volumetric uxes, may come as a result of different ow
regimes not being explicitly accounted for. These are con-
ditions where the slip factor is high, and the ow is expec-
ted to be inhomogeneous. Another contributing effect to the
gas volumetric flux (m/s)

p
r
e
s
s
u
r
e
g
r
a
d
i
e
n
t
(
P
a
/
m
)
0 5 10
0
25
50
75
100
125
150
175
200
present
OLGAS
Beggs & Brill
Friedel
Figure 4: Calculated pressure drop
p
x
as a function of gas
volumetric ux j
g
, for a pressure p = 20bar and
a liquid volumetric ux j

=0.1m/s.
gas volumetric flux (m/s)
l
i
q
u
i
d
h
o
l
d
-
u
p
(

)
0 5 10
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8 p=90 bar, j
l
=1 m/s
p=45 bar, j
l
=1 m/s
p=20 bar, j
l
=1 m/s
p=90 bar, j
l
=1 m/s, =1
0
p=45 bar, j
l
=1 m/s, =1
0
p=90 bar, j
l
=0.2 m/s
p=45 bar, j
l
=0.2 m/s
p=20 bar, j
l
=0.2 m/s
p=90 bar, j
l
=0.1 m/s
p=45 bar, j
l
=0.1 m/s
p=20 bar, j
l
=0.1 m/s
Figure 5: Computed liquid hold-up as a function of gas
volumetric ux.
underprediction may be the presently used wall function for
smooth walls.
3.3 Liquid hold-up and pressure drop by the present
program
Figure 5 shows liquid hold-up calculated by the present pro-
gram as a function of gas volumetric ux, for varying pres-
sure and liquid volumetric ux. Where the inclination angle
is not indicated, the channel was horizontal. It can be ob-
served that the liquid hold-ups fall on three lines; one for
each liquid volumetric ux, virtually independent of pres-
sure. (The curves have been t to the calculated points to
illustrate this.) Naturally, there is a falling tendency as the
liquid volumetric ux decreases. Further, the liquid hold-
up strongly increases when the gas volumetric ux is de-
creased.
Figure 6 gives the calculated pressure drop as a function
of gas volumetric uxes for different pressures and liquid
volumetric uxes. Line segments have been added between
the calculated points to facilitate the reading of the gure.
As expected, pressure drop increased with increasing gas
and liquid volumetric uxes. Further, it decreased for lower
pressures, mainly due to a lower gas density.
The pressure drop was higher for the inclined channel
( =1) because of the gravitational effect.
7
gas volumetric flux (m/s)

p
r
e
s
s
u
r
e
g
r
a
d
i
e
n
t
(
P
a
/
m
)
0 5 10
0
100
200
300
400
500
600
700
p=90 bar, j
l
=1 m/s
p=45 bar, j
l
=1 m/s
p=20 bar, j
l
=1 m/s
p=90 bar, j
l
=1 m/s, =1
0
p=45 bar, j
l
=1 m/s, =1
0
p=90 bar, j
l
=0.2 m/s
p=45 bar, j
l
=0.2 m/s
p=20 bar, j
l
=0.2 m/s
p=90 bar, j
l
=0.1 m/s
p=45 bar, j
l
=0.1 m/s
p=20 bar, j
l
=0.1 m/s
Figure 6: Computed pressure drop as a function of gas volu-
metric ux.
measured neg. pressure gradient (Pa/m)
c
a
l
c
u
l
a
t
e
d
n
e
g
.
p
r
e
s
s
u
r
e
g
r
a
d
i
e
n
t
(
P
a
/
m
)
0 100 200 300 400 500 600 700
0
100
200
300
400
500
600
700
j
l
=1 m/s
j
l
=1 m/s, =1 deg
j
l
=0.2 m/s
j
l
=0.1 m/s
90 bar
45 bar
20 bar
y=x + 30 %
30 %
Figure 7: Comparison between experimental and computed
pressure gradient.
3.4 Comparison of all calculations to experimental data
An overview of the results of the present model for the pres-
sure drop as compared to experimental data is given in Fig-
ure 7. All the evaluated data points can be seen in the graph.
The measured values are plotted along the abscissa, and the
calculated values are plotted along the ordinate. Therefore,
ideally, all the points should fall on the line y = x. This is,
however, not the case, due to the simplications and uncer-
tainties in the model, and to some extent also due to experi-
mental uncertainty.
Figure 8 focuses on the low-velocity data. Figures 7
8 show that most of the data points t between the 30%
lines. However, for the pressure drops below about 10Pa/m
(with liquid volumetric uxes j

below 0.2m/s and gas

volumetric uxes below j
g
= 1m/s), our model generally
overpredicted the pressure drop. The pressure drop was
more overpredicted for the lower pressures, where the dens-
ity difference between gas and liquid was larger.
The data with a high liquid volumetric ux were rep-
resented fairly well by the present model. For the low li-
quid volumetric uxes, the high-pressure-drop data (high
gas volumetric uxes) were mostly underpredicted, while
there was a tendency to overpredict the low-pressure-drop
data.
Figure 9 gives an overview for the results for the liquid
hold-up, while Figure 10 focuses on low liquid hold-up.
measured neg. pressure gradient (Pa/m)
c
a
l
c
u
l
a
t
e
d
n
e
g
.
p
r
e
s
s
u
r
e
g
r
a
d
i
e
n
t
(
P
a
/
m
)
0 10 20 30 40 50 60
0
10
20
30
40
50
60 j
l
=1 m/s
j
l
=1 m/s, =1 deg
j
l
=0.2 m/s
j
l
=0.1 m/s
90 bar
45 bar
20 bar
y=x + 30 %
30 %
Figure 8: Comparison between experimental and computed
pressure gradient (low-velocity data).
measured liquid hold-up ()
c
a
l
c
u
l
a
t
e
d
l
i
q
u
i
d
h
o
l
d
-
u
p
(

)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8 j
l
=1 m/s
j
l
=1 m/s, =1 deg
j
l
=0.2 m/s
j
l
=0.1 m/s
90 bar
45 bar
20 bar
y=x + 30 %
30 %
Figure 9: Comparison between experimental and computed
liquid hold-up.
measured liquid hold-up ()
c
a
l
c
u
l
a
t
e
d
l
i
q
u
i
d
h
o
l
d
-
u
p
(

)
0 0.05 0.1 0.15 0.2
0
0.05
0.1
0.15
0.2 j
l
=1 m/s
j
l
=1 m/s, =1 deg
j
l
=0.2 m/s
j
l
=0.1 m/s
90 bar
45 bar
20 bar
y=x + 30 %
30 %
Figure 10: Comparison between experimental and com-
puted liquid hold-up (low-liquid-hold-up data).
8
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1 2 3 4 5 6 7 8 9 10 11
Data selection
R
M
S

d
e
v
i
a
t
i
o
n

(

)
Present OLGAS B & B Friedel
Figure 11: RMS deviations for pressure drop. The x-axis
entries are explained in Table 5.
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1 2 3 4 5 6 7 8 9 10 11
Data selection
R
M
S

d
e
v
i
a
t
i
o
n

(

)
Present OLGAS B & B Premoli
Figure 12: RMS deviations for liquid hold-up. The x-axis
entries are explained in Table 5.
The calculated values compare reasonably well to the high-
pressure and high-liquid-volumetric-ux data. However,
there is a tendency towards underprediction for low pres-
sures and low liquid volumetric uxes.
3.5 Deviations between calculated values and
experimental data
The discrepancies between calculated values and experi-
mental data were evaluated by the root-mean-square devi-
ation, dened by:

rms
() =
_
_
1
N
N

i=1
_

calc,i

exp,i

exp,i
_
2
_
_
1
2
. (35)
The absolute deviation was also calculated, and gave the
same trends.
The results for the pressure gradient and the liquid hold-
up are shown in Table 5, and given graphically in Fig-
ures 11 and 12. The column labelled N lists the number
of data points. Note that a data point for, say, a pressure
of p = 90bar and a liquid volumetric ux of j

= 1m/s in
a horizontal pipe will be counted in the categories 1, 2, 5,
9, and possibly 10 and 11. The table column labelled F/P
refers to the Friedel correlation for the pressure drop and
the Premoli correlation for the liquid hold-up. The far-right
column contains data where the calculations of the present
program are compared with those of OLGAS instead of with
the experimental data.
The categories 10 and 11 are somewhat special, and are
further discussed in Section 3.5.3.
3.5.1 Pressure drop
All the tested calculation methods, except OLGAS, gener-
ally overpredicted the pressure drop. As can be observed,
OLGAS had the smallest overall deviation (0.15), whereas
the Beggs and Brill overall deviation was 0.24. This is not
surprising, as those methods have been specically tuned
to match the pressure drop. The present CFD program and
the Friedel correlation showed similar overall results, with
overall deviations of 0.43 and 0.42, respectively.
OLGAS had the smallest deviations for all the data ranges,
except for the inclined pipe (Case 8). (Case 9 consists of the
corresponding data points for a horizontal pipe.) In Case
11, i.e., for for high gas volumetric uxes, the OLGAS pre-
dictions were very close to the experimental data. For the
inclined pipe, the Friedel correlation performed best. The
present model predicted the pressure drop more success-
fully for the highest liquid volumetric uxes ( j

=1m/s).
It should be noted that in OLGAS, a pipe roughness of
2.9 10
5
m was employed. For the present program and
the two engineering correlations, the roughness was set to
zero, as that gave better correspondence to the experimental
data for the pressure drop.
3.5.2 Liquid hold-up
As for the liquid hold-up, the deviations given in Table 5,
and illustrated by Figure 12, show that the best prediction
was given by the Premoli correlation, with an overall RMS
deviation of 0.35. It had a slight trend towards underpredic-
tion. The present model gave a deviation of 0.47, and the
OLGAS code had virtually the same result.
It is interesting to note that the Premoli correlation gave
the best results for the liquid hold-up for the studied data
points, and that the present CFD program performed well,
even though neither of these two models were developed
specically for oil and gas pipelines. The present program
gave an overall deviation equal to that of OLGAS. The li-
quid hold-up predictions of the Beggs and Brill correlation
gave an overall deviation one decade higher than those of
the other models (3.6). (Note that as a result of this, in Fig-
ure 12, the bars representing the Beggs and Brill correlation
signicantly exceed the y axis scale for most of the cases).
This was due to the Beggs and Brill correlations failing to
predict the liquid hold-up for some of the data points with
high gas volumetric uxes j
g
, particularly for high pres-
sures. Moreover, even when these data points were disreg-
arded, the Beggs and Brill correlation showed the largest
deviations.
The Premoli correlation gave the smallest deviations for
all the cases studied, except for the inclined pipe case (Case
8), where OLGAS was best, and Case 5 ( j

=1m/s), where
the present program performed equally well. The present
program performed relatively well for the highest liquid
volumetric uxes and for medium and high pressures.
3.5.3 Liquid hold-up for large Reynolds numbers
This section provides a discussion of the categories 10 and
11 in Table 5.
In category 10, the data points have been discarded where
the measured liquid hold-up
,exp
was smaller than the li-
quid hold-up for homogeneous ow,
,h
. Here, the term
homogeneous ow denotes a ow where the mean velo-
cities of gas and liquid are equal, so that the liquid hold-up
9
Table 5: RMS deviations from experimental data.
Qty. Range
*
N Present OLGAS B & B F/P Present

p
x
(1) All 52 0.434 0.145 0.243 0.420 0.676
(2) p =90bar 18 0.214 0.117 0.179 0.209 0.315
(3) p =45bar 20 0.372 0.135 0.245 0.390 0.522
(4) p =20bar 14 0.665 0.184 0.304 0.619 1.087
(5) j

=1m/s 22 0.141 0.129 0.265 0.256 0.137

(6) j

=0.2m/s 15 0.461 0.142 0.253 0.505 0.667

(7) j

=0.1m/s 15 0.641 0.168 0.196 0.511 1.055

(8) =1 8 0.132 0.127 0.109 0.080 0.113
j

=1m/s, =1
(9) p =90/45bar, 8 0.113 0.110 0.328 0.299 0.132
j

=1m/s, =0
(10)
,exp
>
,h
42 0.467 0.158 0.249 0.458 0.743
(11) Re
g
> 10
7
12 0.230 0.044 0.143 0.131 0.216

(1) All 52 0.468 0.470 3.614 0.346 0.336

(2) p =90bar 18 0.673 0.679 5.540 0.486 0.241
(3) p =45bar 20 0.261 0.337 2.362 0.196 0.306
(4) p =20bar 14 0.368 0.254 1.049 0.293 0.460
(5) j

=1m/s 22 0.136 0.181 0.790 0.132 0.119

(6) j

=0.2m/s 15 0.388 0.404 2.564 0.278 0.392

(7) j

=0.1m/s 15 0.763 0.744 6.148 0.558 0.467

(8) p =90/45bar, 8 0.107 0.072 0.756 0.170 0.100
j

=1m/s, =1
(9) p =90/45bar, 8 0.172 0.212 0.873 0.077 0.109
(10)
,exp
>
,h
42 0.272 0.261 1.074 0.237 0.372
j

=1m/s, =0
(11) Re
g
> 10
7
12 0.832 0.847 7.273 0.570 0.079
*
Every data point matching the condition(s) listed, is counted in each range.

The data in this column is for the present program with OLGAS, and not the experimental data, as reference.
could be calculated according to the formula

,h
=
j

j
g
+ j

. (36)
A liquid hold-up smaller than that of homogeneous ow
would imply that the mean liquid velocity was larger than
the mean gas velocity. In the cases studied here, we con-
sider such a behaviour unphysical. This shows that although
the experiments have been conducted very carefully, exper-
imental uncertainties will be present due to the demanding
nature of these measurements.
Category 11 includes the data points where the gas Reyn-
olds number Re
g
> 10
7
. It is dened by
Re
g
=

g
j
g
d

g
. (37)
Consider Figure 13. It shows the absolute relative devi-
ation (i.e., N = 1 in Equation (35)) between liquid hold-up

calculated by the present program and experimental data,

as a function of gas Reynolds number Re
g
. It can be seen
that for low Reynolds numbers, the deviation has a mostly
falling trend as the Reynolds number increases. However, at
a Reynolds number of about Re
g
10
7
, the deviations start
to increase, and for some of the data points very much so.
Figure 14 shows a similar plot to that of Figure 13, but here,
the reference is the OLGAS calculations. For low Reynolds
numbers, the deviations are relatively high, and higher than
those of Figure 13. Hence, the liquid hold-up calculated
by the present program is closer to the experimental data
than to the OLGAS calculations in that range. However, as
the gas Reynolds number increases, the deviations decrease,
and as opposed to the experimental data-based deviations,
they continue to be low for the highest Reynolds numbers.
This effect is illustrated by Case 11 (Re
g
> 10
7
) in
Table 5. All the calculation methods gave high deviations
gas Reynolds number ()
a
b
s
o
l
u
t
e
r
e
l
a
t
i
v
e
d
e
v
i
a
t
i
o
n
(

)
0.0x10
+00
5.0x10
+06
1.0x10
+07
1.5x10
+07
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
p=90 bar, j
l
=1 m/s
p=45 bar, j
l
=1 m/s
p=20 bar, j
l
=1 m/s
p=90 bar, j
l
=1 m/s, =1
0
p=45 bar, j
l
=1 m/s, =1
0
p=90 bar, j
l
=0.2 m/s
p=45 bar, j
l
=0.2 m/s
p=20 bar, j
l
=0.2 m/s
p=90 bar, j
l
=0.1 m/s
p=45 bar, j
l
=0.1 m/s
p=20 bar, j
l
=0.1 m/s
Figure 13: Absolute relative deviation between liquid hold-
up calculated by the present program and exper-
imental data, as a function of gas Reynolds num-
ber.
10
gas Reynolds number ()
a
b
s
o
l
u
t
e
r
e
l
a
t
i
v
e
d
e
v
i
a
t
i
o
n
(

)
0.0x10
+00
5.0x10
+06
1.0x10
+07
1.5x10
+07
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
p=90 bar, j
l
=1 m/s
p=45 bar, j
l
=1 m/s
p=20 bar, j
l
=1 m/s
p=90 bar, j
l
=1 m/s, =1
0
p=45 bar, j
l
=1 m/s, =1
0
p=90 bar, j
l
=0.2 m/s
p=45 bar, j
l
=0.2 m/s
p=20 bar, j
l
=0.2 m/s
p=90 bar, j
l
=0.1 m/s
p=45 bar, j
l
=0.1 m/s
p=20 bar, j
l
=0.1 m/s
Figure 14: Absolute relative deviation between liquid hold-
up calculated by the present program and by OL-
GAS, as a function of gas Reynolds number.
in this case. This is owing to the fact that all the data
points
1
with a measured liquid hold-up smaller than the
homogeneous-ow liquid hold-up t into this category. It
can be observed that the RMS deviations for the calcula-
tions of the present program and for those of OLGAS are
nearly equal, with a value of 0.83 and 0.85, respectively,
whereas the RMS deviation between the calculated results
of the present program and those of OLGAS (last column in
Table 5) is signicantly smaller, with a value of 0.08. This
shows that neither calculation method captured the meas-
ured data in this case. A further illustration can be given
by comparing the Cases 1 and 10: All calculation methods
have a signicantly smaller deviation fromthe experimental
data when the
,exp
<
,h
data points are not considered.
4 CONCLUSIONS
We have presented a comparison of pressure drop and liquid
hold-up calculated using different modelling strategies; a
two-dimensional two-uid model, a one-dimensional sim-
ulator, and engineering correlations. The numerical results
have been compared with experimental data.
In the present work, a two-uid model has been imple-
mented in the framework of a two-dimensional multiphase
Computational Fluid Dynamics (CFD) code. The govern-
ing equations were spatially discretized using the nite-
volume technique, and the time-integrator was an explicit
low-storage ve-step fourth-order Runge-Kutta scheme.
Although the main purpose of a multiphase CFD model is
to facilitate the design of process equipment, an important
step in its validation is the comparison with good experi-
mental data in simple geometries. The test case of high-
pressure two-phase ow in pipes has been selected since it
is of importance to the oil and gas industry.
Numerical results from the present program, the OLGAS
simulator, the Beggs and Brill correlation, and the Friedel
and the Premoli correlations have been compared to 52 data
points from the TILDA two-phase pipe ow database. The
data were taken in a 0.189m inner-diameter pipe. The liquid
volumetric ux j

varied between 0.1 and 1.0m/s, the gas

volumetric ux j
g
was in the range from 0.5 to 12m/s, and
pressures p from 20 to 90bar were employed. Inclination
1
Except one
angles of 0 and 1 were used, with the majority of the data
coming from the horizontal pipe.
The present program performed well at the higher pres-
sures and at high liquid volumetric uxes.
The liquid hold-up was slightly overpredicted for high
gas volumetric uxes and slightly underpredicted for low
volumetric gas uxes. For low pressures, the underpredic-
tion at low volumetric uxes was signicant.
At high pressure and high liquid volumetric uxes, the
calculated pressure drop compared well to the experimental
data.
For low pressures and low liquid volumetric uxes, the
pressure drop was increasingly underpredicted for increas-
ing gas volumetric uxes. These are conditions where the
slip factor is high, and the ow is expected to be inhomo-
geneous, so the underpredicted pressure drop may in part be
due to different ow regimes not being explicitly accounted
for.
The deviations between experimental data and model
predictions are given in Table 5 on the preceding page. For
the liquid hold-up, the predictions of the Friedel correlation
had the lowest overall deviation. Interestingly, the present
program had the same overall deviation as that of OLGAS.
The Beggs and Brill correlation was not successful in pre-
dicting the liquid hold-up.
The smallest deviations between experimental and calcu-
lated pressure drops were obtained by the OLGAS program
and the Beggs and Brill correlation. This is unsurprising,
since the parameters in these models have been tuned to
match experimentally observed pressure drops. The present
program performed equally to OLGAS for high pressures
and high liquid volumetric uxes.
In some of the data points, the measured liquid hold-up
was lower than that of homogeneous ow, and might be er-
roneous. The overall deviation for the liquid hold-up was
signicantly reduced for all the calculation methods when
these data points were disregarded.
We nd it interesting to note the relatively good corres-
pondence between experimental data and the computed res-
ults of the CFD program, particularly when considering the
simple constitutive relations employed.
ACKNOWLEDGEMENTS
The authors especially thank Kjell Arne Jacobsen of SIN-
TEF Petroleum Research for providing data from the TILDA
database. We would like to express our gratitude to Lars
Hovden of Scandpower AS for performing the OLGAS cal-
culations. We would also like to thank Professor Morten
Chr. Melaaen at the Process Equipment Research Group at
Telemark University College for his co-operation and useful
suggestions for this article.
This work has been sponsored by the Research Council
of Norway and by SINTEF Energy Research.
REFERENCES
Alajbegovi c, A., Drew, D. A., and Lahey, R. T. An ana-
lysis of phase distribution and turbulence in dispersed
particle/liquid ows. Chemical Engineering Communic-
ations, volume 174: pages 85133, 1999.
Amsden, A. A. and Harlow, F. H. KACHINA: An Eulerian
computer program for multield uid ows. Technical
11
Report LA-NUREG-5680, Los Alamos Scientic Labor-
atory, December 1974.
Amsden, A. A. and Harlow, F. H. K-TIF: A two-uid com-
puter program for downcomer ow dynamics. Technical
Report LA-NUREG-6994, Los Alamos Scientic Labor-
atory, January 1978.
Beggs, H. D. and Brill, J. P. A study of two-phase ow in
inclined pipes. Journal of Petroleum Technology, pages
607617, May 1973.
Bendiksen, K. H., Malnes, D., Moe, R., and Nuland, S. The
dynamic two-uid model OLGA: Theory and applica-
tion. SPE Production Engineering, pages 171179, May
1991.
Brill, J. P. and Mukherjee, H. Multiphase ow in wells,
volume 17 of Henry L. Doherty series Monograph /
SPE, chapter Multiphase ow pressure gradient predic-
tion. Henry L. Doherty Memorial Fund of AIME, Society
of Petroleum Engineers, Richardson, Tex, USA, 1999.
Brown, G. J. Erosion prediction in slurry pipeline tee-
junctions. Applied Mathematical Modelling, volume 26,
no. 2: pages 155170, February 2002.
Carpenter, M. H. and Kennedy, C. A. Fourth-order 2N-
storage Runge-Kutta schemes. Technical Report TM
109111, NASA, April 1994.
Drew, D. A. and Passman, S. L. Theory of Multicomponent
Fluids. Number 135 in Applied Mathematical Sciences.
Springer-Verlag, New York, 1999. ISBN 0-387-98380-5.
Ferziger, J. H. and Peri c, M. Computational methods for
uid dynamics. Springer-Verlag, Berlin, second edition,
1999. ISBN 3-540-65373-2.
Friedel, L. Improved friction pressure drop correlations for
horizontal and vertical two phase pipe ow. In: Proceed-
ings, European Two Phase Flow Group Meeting. Ispra,
Italy, June 1979. Paper E2.
Harlow, F. H. and Amsden, A. A. Numerical uid dynamics
calculation method for all ow speeds. Journal of Com-
putational Physics, volume 8, no. 2: page 197, 1971.
Idelchik, I. E. Handbook of hydraulic resistance, chapter
2.2 Diagrams of friction coefcients. CRC Press, Boca
Raton, Florida, USA, third edition, 1994. ISBN 0-8493-
9908-4.
Ishii, M. Thermo-uid dynamic theory of two-phase ow.
Collection de la Direction des Etudes et Recherches
dElectricit de France. Eyrolles, Paris, 1975.
Kataoka, I. Local instant formulation of two-phase ow.
International Journal of Multiphase Flow, volume 12,
no. 5: pages 745758, 1986.
Kennedy, C. A., Carpenter, M. H., and Lewis, R. M. Low-
storage, explicit Runge-Kutta schemes for the compress-
ible Navier-Stokes equations. Applied Numerical Math-
ematics, volume 35, no. 3: pages 177219, November
2000.
Launder, B. E. and Spalding, D. B. The numerical compu-
tation of turbulent ows. Computer Methods in Applied
Mechanics and Engineering, volume 3: pages 269289,
1974.
Liles, D. R. and Reed, W. H. Semi-implicit method for two-
phase uid-dynamics. Journal of Computational Physics,
volume 26, no. 3: pages 390407, 1978.
Lin, A. and Lopez, D. Two-uid model of wavy separ-
ated two-phase ow. International Journal of Multiphase
Flow, volume 23, no. 6: pages 11311146, 1997.
Melaaen, M. C. Calculation of uid-ows with staggered
and nonstaggered curvilinear nonorthogonal grids The
theory. Numerical Heat Transfer Part B Fundamentals,
volume 21, no. 1: pages 119, JanuaryFebruary 1992.
Moura, L. F. M. and Rezkallah, K. S. Numerical study
on the two-phase ow distribution in a T-junction. In-
ternational Journal for Numerical Methods in Fluids,
volume 17: pages 257270, 1993.
Newton, C. H. and Behnia, M. A numerical model of strat-
ied wavy gas-liquid pipe ow. Chemical Engineering
Science, volume 56, no. 24: pages 68516861, Decem-
ber 2001.
Nigmatulin, R. I. Dynamics of Multiphase Media, volume 1.
Hemisphere Publishing Corporation, New York, 1991.
ISBN 0-89116-316-6.
Patankar, S. V. and Spalding, D. B. Calculation procedure
for heat, mass and momentum-transfer in 3-dimensional
parabolic ows. International Journal of Heat and Mass
Transfer, volume 15, no. 10: page 1797, 1972.
Premoli, A., Francesco, D., and Prina, A. A dimensionless
correlation for determining the density of two-phase mix-
tures. La Termotecnica, volume 25: pages 1726, 1971.
In Italian.
Pun, W. M., Spalding, D. B., Rosten, H., and Svensson,
H. Calculation of two-dimensional steady two-phase
ows. In: Two-Phase Momentum, Heat and Mass Trans-
fer in Chemical, Process, and Energy Engineering Sys-
tems, F. Durst, G. V. Tsiklauri, and N. H. Afgan, edit-
ors, volume 1 of Thermal and Fluids Engineering, pages
461470. Hemisphere, Washington, USA, 1979. ISBN
0-89116-104-X.
Rhie, C.-M. and Chow, W. L. Numerical study of the tur-
bulent ow past an isolated airfoil with trailing edge sep-
aration. AIAA Journal, volume 21, no. 11: pages 1525
1532, November 1983.
Rivard, W. C. and Torrey, M. D. K-FIX: A computer
program for transient, two-dimensional, two-uid ow.
Technical Report LA-NUREG-6623, Los Alamos Sci-
entic Laboratory, April 1977.
Saurel, R. and LeMetayer, O. A multiphase model for com-
pressible ows with interfaces, shocks, detonation waves
and cavitation. Journal of Fluid Mechanics, volume 431:
pages 239271, March 2001.
Schlichting, H. Boundary-Layer Theory. McGraw-Hill,
New York, 7 edition, 1979.
12
Sha, W. T. and Soo, S. L. Brief communication: On the
effect of P term in multiphase mechanics. Inter-
national Journal of Multiphase Flow, volume 5, no. 2:
pages 153158, April 1979.
Slattery, J. C. Flow of viscoelastic uids through porous
media. AIChE Journal, volume 13, no. 6: pages 1066
1071, November 1967.
Soo, S. L. Particulates and continuum: Multiphase uid
dynamics. Hemisphere, New York, 1989. ISBN0-89116-
918-0.
Soo, S. L. Multiphase uid dynamics. Science Press,
Beijing, 1990. ISBN 0-291-39781-6, 70-300-0102.
Spalding, D. B. The calculation of free-convection phenom-
ena in gas-liquid mixtures. In: Heat Transfer and Turbu-
lent Buoyant Convection, D. B. Spalding and N. Afgan,
editors, volume 2 of Thermal and Fluids Engineering,
pages 569586. Hemisphere, Washington, USA, 1977.
ISBN 0-07-059926-2.
Spalding, D. B. Numerical computation of multi-phase uid
ow and heat transfer. In: Recent Advances in Numer-
ical Methods in Fluids, C. Taylor and K. Morgan, edit-
ors, volume 1, pages 139168. Pineridge Press, Swansea,
UK, 1980. ISBN 0-906674-07-7.
Whitaker, S. Advances in theory of uid motion in
porous media. Industrial and Engineering Chemistry,
volume 61, no. 4: pages 1428, December 1969.
Williamson, J. H. Low-storage Runge-Kutta schemes.
Journal of Computational Physics, volume 35, no. 1:
pages 4856, 1980.
A THE RANGE OF THE DISPLACEMENT FACTOR
B
k
Sha and Soo (1979) and Soo (1990, pages 319321) main-
tained that the displacement factor B
k
can attain values in
the range 01. It can, however, be shown that B
k
will be
close to unity. Consider Equation (12). By rearranging, and
using Equation (7), we get:
B
k
=1
1
V
_
A
k
(p
k
I +
k
) n
k
dAF
d
1
V
_
A
k
i
p
k
I n
k
dA
(38)
In the numerator on the right-hand side of Equation (38),
the terms will tend to cancel. Further, each of them will be
small compared to the denominator. Hence B
k
will be close
to unity.
13