7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
The effect of Surface Tension and Diffusion on Onedimensional
Modelling of Slug Flow
R.I. Issa* and M. Montini*
Department of Mechanical Engineering, Imperial College London, Exhibition Road, SW7 2AZ, London, UK
r.issa @imperial.ac.uk and m.montini07@0imperial.ac.uk
Keywords: Slug flow, Twofluid model, Surface tension, Diffusion, Wellposedness
Abstract
This paper investigates the use of additional closure relationships, such as surface tension and diffusion, with the
aim of obtaining a wellposed system of equations for the transient one dimensional two fluid model in the slug flow
regime. It is shown that these forces are able to damp short wavelength perturbations and bound the growth rate of
the instabilities. It was found that only with the introduction of diffusive terms in all the governing equations of the
system it is possible to obtain results converging to a unique solution regardless of the mesh used in their numerical
solution and hence the system is deemed to be wellposed. However, more work is needed to determine a universal
relationship for the diffusion coefficients necessary to render the system wellposed for all the cases in the slug flow
regime.
Nomenclature
Roman symbols
A cross sectional area (mn2)
A jacobian matrix of coefficients
B jacobian matrix of coefficients
C vector of algebraic terms
C friction factor coefficient
D diameter (m)
F friction force (Nm 3)
F algebraic terms (Nm 3)
f friction factor
g gravitational constant (ms 2)
h height of the phase (m)
i imaginary unit
K stability constant
k wavenumber (m 1)
L length of the pipe (m)
m friction factor power
P cross sectional averaged pressure (Nm 2)
p interfacial pressure (Nm 2)
R radius (m)
Re Reynolds number
S "wetted" perimeter (m)
u velocity (ms 1)
v vector of independent variables
V wave velocity (ms 1)
Greek symbols
a volume fraction
13 pipe inclination ( )
r mass diffusivity (in2s 1)
7 stratification angle ( )
A wavelength (m)
An characteristic (m/s)
p dynamic viscosity (kgslim1)
v kinematic viscosity (n2s 1)
p density (kgm 3)
a surface tension coefficient (Nm 1)
T shear stress (Nm 2)
p derivative of the volume fraction
I' dynamic mass diffusion coefficient (kgs lm 1)
Uc angular frequency (s 1)
Superscripts
s superficial
T transposed
Subscripts
0 initial
9 gas
i interface
IV inviscid
k generic phase
1 liquid
n number of variables
rel relative
V viscous
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Introduction
The onedimensional twofluid model represents a sys
tem of partial differential equations (PDEs), which at
tempts to mathematically describe and simulate the
fundamental properties of the different regimes in multi
phase flow.
Despite the threedimensional nature of twophase
flow, and in particular of phenomena such as turbulence
and the presence of slugs, the onedimensional approach
is still the most practicable for modelling the flow in
long pipelines in view of the current computational
capabilities.
The variations for the flow properties of interest
are considered to occur mainly in the axial direction;
therefore averaging the governing equations over the
cross section of the pipe is considered a reasonable
approximation. However, the averaging process leads to
loss of information about the flow variations in the radial
direction as well as of interfacial forces. This can be
compensated by introducing closure terms which model
forces or effects that are subsumed in the averaging
process, such as surface tension and friction.
It is well known that the system of equations for
the onedimensional twofluid model is illposed unless
appropriate closures are introduced and numerous inves
tigations have focused on this fundamental issue (among
others Arai (1980), Holmas et al. (2008), Issa & Kempf
(2003), Ramshaw & Trapp (1978)). The hyperbolicity
of the system, hence the wellposedness of the problem,
is still an unsolved and challenging issue, which plays a
key role in the reliability of the solution of the system
and, therefore, of the validity of results computed by
numerical codes.
The aim of the present work is to explore the possi
bility of extending the wellposedness of the problem by
introducing new closures into the system that relate to
the effects of surface tension and axial diffusion. These
closures are first examined from a mathematical point
of view; they are then implemented in a numerical code
with the aim of investigating the applicability of the new
models in the slug flow regime.
Mathematical model
The one dimensional two fluid model considered in
this paper is based on the following equations for the
conservation of mass and momentum
a09p, 9 .. ..
Ot ox
a9 lp1 aolplul
Ot Oax
0, (1)
0, (2)
9x { Ig 9x )
09X 0 )
9 g
D

T1
ul
Figure 1: Pipe crosssectional area and lateral view with
relevant properties
9 .. .. 0 .. ?
dt x
X ([a
09olplul + 9a,?
Ot a'x
0i ( ai
Pig i i 11
ag cagpYgcos 13x
,,I,.,, ii F Fi, (3)
OPil
09X:
i '1
alpigcos13
09x
alplgsin F + Fi, (4)
pi = a 2, (5)
where the subscripts g and I represent the gas and liquid
phases respectively, while i defines quantities at the
interface between the two phases.
Referring to figure 1, in these equations ak (either
k = g or k = 1) represents the generic phase volume
fraction with the condition ag + a = 1 and Ak Aak
is the area occupied by the phase, where A is the cross
sectional area of the pipe; pk the phase density, Uk the
velocity and Pik the interfacial pressure for each phase.
Furthermore, g is the gravitational acceleration and 3 is
the angle between the pipe and the horizontal line.
The last term in each continuity equation represents
the axial diffusion term (Holmas et al. (2008)), where
k = PkFk is the dynamic mass diffusion coefficient
and Fe is the mass diffusivity. The second term on the
right hand side of the momentum equations represents
T9
the hydrostatic pressure variation (Wallis (1969)), wl
hi is the height of the liquid interface. The third t
on the right hand side of each momentum equa
represents axial diffusion (Arai (1980), Holmas e
(2008)), where Pk is the dynamic viscosity.
If the interface between two phases is curved,
interfacial pressure in each phase is different, the
ference being due to the surface tension. This ef
is represented by equation (5), where a is the sur
tension coefficient (R.lin.ls & Trapp (1978)).
Finally, the terms F,, F, and F, represent the
tional forces between each phase and the wall an(
the gasliquid interface respectively. Since the f
under consideration is either in the stratified or slug f
regimes, these forces may be modelled in genera
follows
9
Si
F, 
Al
^^Q~r)71
\^; ^gi)
1
 2 'f .., I.,
1
Ti 2f fipiul ul
2.
1
2
ua;I.. U1 ,
where Sk are the "wetted" perimeters, Tk represent
shear stresses and fk are friction factors.
Although the system is almost entirely mechanisti
requires closure models for these three friction fac
which are based on empirical correlations. Exten
studies on different correlations can be found in
literature (Lin & Hanratty (1986), Spedding & H
(1997), Issa & Kempf (2003)). However, these te
are all based on Blasiuslike equations that can be
in a generic form as
lere
erm
tion
t al.
the
dif
fect
Face
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Table 1: Geometry and flow properties
Pipe length Pipe diameter
L 36.0m D = 0.078m
Air density Water density
pg 1.253kg/m3 pi 998.2kg/m3
Air viscosity Water viscosity
Pg 1.79 10 5kg/ms /i1 1.14 10 3kg/ms
inc while the coefficients ,.. and m, are equal to 1 for
d at laminar flow or 0.25 if the flow becomes turbulent. In
low the case of the liquidwall friction factor the correlation
low used is the one suggested by Spedding & Hand (1997)
1 as where the friction factor fi is based on the superficial
liquid velocity u = ul al and the actual pipe diameter
D (rather then the hydraulic diameter) and the two
(6) coefficients are Ci 0.0262 and ml 0.139 for the
turbulent case (Ref > 2100) and C1 24 and ml 1
(7) for the laminar case (Ref > 2100).
All the data presented in this paper relate to ex
periments performed on the WASP facility in the De
partment of Chemical Engineering at Imperial College
London for flows of water and air at atmospheric tempe
(8)
rature and pressure. Relevant physical properties of the
the flows are summarised in table 1.
c, it Wellposedness
tors
sive The entire mathematical model (1)(5) is a system of
the partial differential equations and can be rewritten in a
and more compact matrix form as follows
rms
cast
A(v) v + B(v) v + C(v)
8t 09x
f ( D ,, 4A9
f, ( D, A (9)
(Di, (Dia 4A)
\ S, + Si
f C Dl Di (10)
(D I., ull ( m
where D1 and D, are the hydraulic diameters and the
coefficients Ck and mk take different values for the
different correlations as well as whether the regime is
laminar or turbulent.
In all the analyses presented here, the gaswall and
interfacial correlations used are those suggested by Tai
tel & Dukler (1976), where the coefficients C, and C,
are set to 16 in the case of laminar flow (Re, < 1180
and Re, < 1180 respectively) and are equal to 0.046 if
the flow is turbulent (Re, > 1180 and Re, > 1180),
where vT = [a,, u. l,/' pil] is a column vector of
the independent variables and A and B are the Jacobian
matrices of dimension n x n (where n 5 is the number
of independent variables) containing the coefficients of
the differential equations while C is a column vector of
dimension n which contains the algebraic terms.
For the purpose of examining the mathematical prop
erties of the system as presented herein, the two phases
are treated as incompressible, immiscible and isothermal
if not stated otherwise. This is a usual assumption for
such analyses and does not invalidate the findings in
any way as was found for example by Bonizzi & Issa
(2003). In any event, the simulations themselves that
are presented later, are carried out with compressibility
of the gas phase taken fully into account.
In order to verify the wellposedness of the system,
two approaches may be taken: a characteristics analysis
and a linear stability analysis.
The characteristics A, of the system are defined by
det(A AB) = 0
(13)
This equation has n roots A, for the n equations of
the system and on this basis it is possible to define a
classification for the system of PDEs. If all these roots
are real and distinct, the set of differential equations is
hyperbolic or if they are real and equal then the system
is parabolic, in both of which cases the problem is well
posed (Arai (1980), Holmas et al. (2008), Issa & Kempf
(2003), Ramshaw & Trapp (1978)). In the case of any
root being complex, the system of equations is elliptic,
in which case the equation set is illposed for a time
dependent problem.
Equation (13) is usually called the characteristics
dx B
equation and the roots A= ,B of the system (12)
dt A
determine the characteristic velocities along which the
physical information propagates through the flow fields.
Consequently the information travels at a certain finite
speed in a number of directions equivalent to the number
of real characteristics of the system.
Furthermore, it is worth noting that in equation (12)
the algebraic terms in the vector C do not feature in
the characteristics analysis, therefore the choice of the
friction factors does not affect the wellposedness of the
system.
To begin with, the effect of diffusion and surface
tension is not taken into account (Fk = k a 0).
Moreover, since the two interfacial pressures become
equal, the pressure must be considered the same at the
gasliquid interface pil = = P.
The system of governing equations (13) is, therefore,
reduced to four equations in four unknowns vT
[al,. ul, P]. By solving it, it is now possible to obtain
four roots, two of which are null A1,2 0. These
latter two are zero because of the assumption of in
compressible flow, while the remaining two correspond
to the continuity waves which are of the same order
of magnitude of the phasic velocities (A3 and
A4 .,).
The system of equations is wellposed only if the
four roots are real and therefore it is possible to find
the following criterion for the wellposedness of the
problem
( ul) < lp + Pl (Pl
Pip9
A
pg)g cos f .
dh1
(14)
For a stratified flow in equilibrium, the well
posedness condition (14) is graphically represented on
the flow regime map of Taitel & Dukler (1976) in figure
2 by the continuous thick line. In the area below the
line the system of equations is hyperbolic, hence case
C
t
C,
a
05
_~o.
0 0
C^
v?
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Bubbly flow Wellposedness limit ..
....... ... ................... 0 Illposed point
[ Wellposed point
Slug flow
case II
        
Stratified wavy
Stratified smooth
Annular flow
10
Superficial gas velocity, U [m/s]
Figure 2: Flow pattern map for horizontal pipes with
the wellposedness limit
Table 2: Flow conditions for cases I and II
I7f 01
case I 6.000 0.400 0.566
case II 6.532 0.532 0.606
I indicated by a square is wellposed. Conversely, in
the area above the line the system is characterized by its
nonhyperbolicity and case II symbolised by a circle on
the map is, therefore, illposed.
The flow conditions of these two cases which will be
considered from now on as test cases are summarised in
table 2.
Next we consider the model to include the surface
tension force but we continue to leave out diffusion. To
determine the characteristics of the system with surface
tension, the second order derivative must be reduced to
a first order form. This can be achieved by introducing
a new variable = a and replacing equation (5)
with the following two expressions (Rainish,.l & Trapp
(1978))
(A 00
P1 r a (15)
Pil  dAl 9X(
dhi
and
= 0. (16)
ax
Therefore, the system (12) has six equations and six
unknowns, where the column vector of independent
variables becomes
v [al, .. ul Pi, ]T. (17)
The characteristics analysis of the governing equations
including these terms gives six real roots, all of which
I I I I I I I I
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
are zero: such a result is unfortunately misleading since
the mixture of first and higher order derivative terms and
the way in which these terms are commonly reduced to
first order (for example in Arai (1980), Holmas et al.
(2008), Ramshaw & Trapp (1978)) seems to invalidate
the characteristics analysis. More specifically, in the
limit for ca 0 the characteristics analysis is not able
to recover the criterion for the wellposedness of the
system without surface tension (given in equation (14)).
If the surface tension is suppressed (a 0) and the
effects of diffusion for both mass and momentum are
reinstated as suggested by Holmas et al. (2008), a similar
treatment to that described in equations (15)(16) can be
applied to the diffusion terms in equations (1)(4). The
second order derivative has again been reduced to a first
order form and the characteristics analysis once more
gives six null roots (Holmas et al. (2008)).
Once again the result is clouded by the presence of a
mixture of first and second order derivatives. Indeed
also in this case, the wellposedness condition in the
limit for the diffusion coefficients which tend to zero
(pg,i, Fr,1 0) is not the same as the one found without
considering the diffusion terms. Hence the analysis is
inconclusive because of the way with which the second
order derivative terms have been reduced to first order.
The conclusions drawn above with respect to the
inclusion of either surface tension or diffusion terms
in the characteristics analysis hold also in the case of
inclusion of both effects. Hence, in these cases, the
analysis of the characteristics does not shed any light on
the issue of the wellposedness of the system.
Stability
The stability of stratified flow has been extensively
analysed in the past in connection with both wave
generation and transition to slug flow. This is of
particular interest to the study of slug flow as the main
mechanisms responsible for the transition from stratified
to slug flow is the growth of natural hydrodynamic
instabilities.
The present approach to the stability analysis follows
the work by Bamea & Taitel (1993) and by Lin &
Hanratty (1986), who used a linear KelvinHelmholtz
(KH) analysis to study the onset of slugging for flows in
a pipe.
In this linear analysis, a small disturbance 6v is added
to the unperturbed solution v. Hence the solution vector
v in equation (12) is replaced by v + Sv and the result
can be linearised with respect to the perturbation 6v, as
follows
ayv 98v ( A\ av
A(v). + B(v). 0 v 0 0v
iB iB ivv 9C
+ V v + v =0. (18)
K av a) x S v
The perturbation is then cast in the form of a Fourier
series given by
6v = 6voc(x t) (19)
where k is the wavenumber and c is the angular
frequency. Equation (18) then becomes
icA(v) vo + ikB(v) 6vo + v a v
09v at
09B\ 0v 0C
+ V. + v v a
i~v < }3r i~v
0, (20)
which has a non trivial solution if the determinant of the
coefficient matrix vanishes. The dispersion relation is,
therefore, determined by
det(iwA + ikB + D) = 0,
where
(d A v\ 0 fdB Ov) (OC)
D v at v ) v (22) v
Therefore, the dispersion equation for the complete
system of governing equations becomes
W2 2(akib+irk2) dk4 +ck2 + i(sk3 k) = 0,
(23)
where
a= P
1 p1lul( +, .
1 aiF g j F
b 2p = ,9uU .\u
Sr v 49u ,ai_
c 
P al
9
1\ 4 aft
S= 1 a + 
P dAl
1 C0F u'
Ag cos 3
 (p P) dA
dh/
1r PL +, + p+ir + 's
2p a g ya Oag
1 p U '/ ... lUl_ ) ,
r \ + + _
P \ al ay al ag
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
p= +
O1 O1
F + i 9 + TS +)
f A1 A Ai A9
(pi p) g sin 3.
(31)
(32)
4
The system is unstable whenever the imaginary part
of the roots w(k) of the dispersion equation (23) is
negative, which implies an exponential growth in time
of the perturbation 6v in equation (19). On the other
hand, if the imaginary part of w(k) is positive, the
system is stable and every perturbation is damped and
therefore decays. The condition for neutral stability can
be found from the dispersion equation for the special
case in which Im[w(k)] = 0, hence the stability
criterion for the viscous case becomes
(Vv Viv)2 + (c a2) dk2 < 0, (33)
Swhere
where V 
2(b
sk3
and Viv a.
 rk2)
Here Vv and Viv are the critical wave velocities at
the onset of the instability for the viscous and inviscid
analysis respectively.
Moreover, it is important to highlight that the stabi
lity of a system of equations closely depends on the
specific wavenumber k considered. This means that if
Im[w(k)] < 0 for at least one specific disturbance of
wavenumber k then the whole flow is unstable. On the
other hand, if Im[w(k)] > 0 for each and all the values
of k then the flow is stable.
By substituting expressions (24)(32) in (33) in the
simplest case without surface tension and diffusion,
the KH criterion for neutral stability can be eventually
expressed as follows (Barnea & Taitel (1993))
U < KPlag + Pgal (P1
Urel < K P (pPG
V PIP 9
A
pg)g cos3 dA (34)
dhi
where ue = ., ul is the slip velocity and for the
inviscid case the constant K = 1, while for the viscous
case
(Vv VwV)2
K= 1 v v (35)
P1 P s A
SpdA
\ P dhl
As a further result of the stability analysis, it is
possible to obtain from the dispersion equation (23) the
growth rate of the instabilities as a function of the wave
length. As shown in figure 3, for different liquid and gas
phase velocities, one observes that the rate at which an
instability characterized by a certain wavelength grows
can vary considerably, hence generating varying wave
and slug structures since the time necessary for the
1uuuu
C
ao
n nni
S case II: illposed
 case I: wellposed
1 1 I """I 1 1

le06
le06
0.001 1
dimensionless wavelength, WD []
Figure 3: Growth rate for the instabilities
instability to reach the top of the pipe can also vary
considerably.
In the case of a wellposed system the growth rate
of instabilities must be bounded for all possible wave
lengths; furthermore, infinitesimally short waves must
also be damped. This is a necessary, but not sufficient,
condition for the wellposedness of the problem. On the
other hand for an illposed case the growth rate of the
instabilities is unbounded and it grows indefinitely as the
wavelength becomes smaller. A small perturbation in
the initial condition can lead to very different solutions,
therefore failing to depend continuously on the initial
data (Prosperetti & Tryggvason (2007)).
Given the conclusion reached earlier that the charac
teristics analysis in the presence of high order derivative
terms is dubious, the results of the stability analysis can
be used instead to investigate the wellposedness of the
system in the presence of surface tension and diffusion.
Since the two test cases I and II are very close in the
flow pattern map (see figure 2) and their flow conditions
are similar, their slug characteristics are also similar.
Therefore, the actual behaviour of the instabilities which
are at the origin of the slugs (hence the growth rates from
the stability analysis) must be comparable. Bearing in
mind this reasoning, the ideal condition to render the
system wellposed is that in which a new force or effect
is able to bound the growth rate of the instabilities of
case II in an analogous way as it happens for case I.
The introduction of the surface tension term modifies
the neutral stability limit as follows (see for example
Barnea & Taitel (1993), Prosperetti & Jones (1987),
Ramshaw & Trapp (1978))
Ure < K p/la + py91l [(p1
V PIPg
pg)g cos 3 + ak2] dA.
(36)
(36)
The contribution of the new term depends on the
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0.01 0.1 1 10 100
dimensionless wavelength, X/D []
Figure 4: Growth rate for the instabilities with surface
tension
wavenumber k, hence on the wavelength of the
instabilities and is specifically important for short
wavelengths disturbances only. On the other hand, in
the approximation of long waves (k 0) it becomes
irrelevant (Barnea & Taitel (1993)).
Figure 4 shows the behaviour of the growth rate of
the instabilities as a function of the wavelength for
case II with surface tension. The experimental value
of this coefficient for water is a 0.072. The
surface tension suppresses short wavelength perturba
tions (negative growth rate) and it limits their growth
rate to a certain value, while it does not affect longer
wavelengths as expected. However, in this case the
peak value of the growth rate of the instabilities is
approximately two order of magnitude higher than the
maximum value computed for case I, which is well
posed with the standard model without surface tension.
Moreover, the growth rate of the instabilities keeps
growing until it reaches a certain cutoff wavelength
without stabilising to a certain plateau value as for case
I.
Figure 5 shows the behaviour of the growth rate
curves as a function of the dimensionless wavelength of
the instabilities for case II including viscous terms in the
gas and liquid momentum equations only, where both
surface tension and mass diffusion are suppressed. For
a value of the viscosity coefficient which corresponds
to the molecular viscosity (p/ H 10 3kg(ms) and
p, z 10 5kg(ms) 1) and also for turbulent viscosity
(approximately two order of magnitude bigger than the
molecular one), the shape of the curve is similar to that
of case I with a long plateau at peak value and a cutoff
wavelength beyond which the instabilities are damped.
However, also in these two cases the bounding value of
the growth rate is still extremely high in comparison to
that of case I.
In the last part of the stability analysis, for the sake
Figure 5: Growth rate for the instabilities with viscous
terms
of simplicity, the effect of surface tension in expression
(27) is neglected (a 0) and both the gas and liquid
coefficients are considered equal for both the mass and
momentum diffusion (tl /= = and F1 F, F).
Figure 6 shows the effect of the diffusion terms on
the stability of the system and a combination of mass
and momentum diffusion seems to be able to shape
the stability curve of case II to a behaviour similar to
that of case I. The coefficient for the mass diffusion
(figure 6(a)) has a strong influence on the peak value
of the growth rate of the instabilities: the higher the
value of F the lower the value of the plateau which
bounds the growth rate. On the other hand momentum
diffusion acts on the cutoff wavelength beyond which
the perturbations are stable: the higher the value of p
the longer the cutoff wavelength.
Once the stability analysis has provided an indication
regarding the behaviour of the system, a confirmation of
the wellposedness of the problem can only come from
the results of numerical simulations.
Numerical results
Issa & Kempf (2003) demonstrated that the results for
wave growth given by a numerical solution of the model
equation for an illposed problem are invalid because
the results of the computations do not converge to a
unique solution regardless of mesh density. This can be
demonstrated here for slug flow simulations.
The mathematical model in equations (1)(5) has
been implemented in the research code TRIOMPH by
Issa & Abrishami (1986), which has been described
in Issa & Kempf (2003), and has been tested starting
from wellposed stratified flow conditions. Here the
compressibility of the gas phase is taken into account
by means of the ideal gas law.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
10000
1000
100
10
1
0.1
0.01
0.001 
le06
 i '"' case II
**** case II, F=0.01 p=10 5
S case II,=0.1 P=105
 case II, F=0.5 =105
 case I
  I  

0.0001 0.01 1
dimensionless wavelength, X/D []
(a) Effect of F
0.0001 0.01 1
dimensionless wavelength, /D []
(b) Effect of p
Figure 6: Growth rate for the instabilities with
artificial diffusion for different values of the diffusion
coefficients
Figure 7 presents a clear example with some results
of numerical simulations of slug flow for different mesh
sizes for the two aforementioned test cases when both
surface tension and diffusion are absent. The com
puted timeaveraged slug frequency is plotted against
the size of the grid spacing: for a wellposed system
like in case I, the average slug frequency converges to
a unique value, which is in excellent agreement with
the experimental value. On the other hand, for case II,
which is illposed, a unique solution no longer exists and
consequently computed values of the frequency diverge
as the mesh is refined.
The behaviour of the slug frequency presented in
figure 7 can also be explained looking at the results of
the stability analysis presented in figure 3 for the same
cases. In numerical terms, as the mesh is refined, it
is possible to capture shorter instabilities. In the well
posed case, once the instabilities with the highest value
of the growth rate are captured the solution stabilises
at a certain value, while for an illposed case there
0.6
0.5
0.4
S0.3
C
0.2
0.1
0
0
o case II: illposed
0 case I: wellposed
 experimental value
, I I I ,
0.2 0.4 0.6 i
dimensionless mesh spacing, Ax/D []
Figure 7: Computed slug frequency against mesh
density
will always be shorter and faster instabilities which can
develop into slugs, thus the slug frequency increases
accordingly.
For case II, short wavelength instabilities still play a
dominant role in the simulations, leading to unphysical
and erroneous results. This is due to the fact that the
instabilities predicted by the twofluid model are not
always appropriate solutions reflecting the real flow
conditions, but could also be a manifestation of the
mathematical or numerical instability of the model (Issa
& Kempf (2003)).
Surface tension effects are negligible compared to
the hydrostatic pressure (ok2 < g(pi p,)) and
therefore they only affect very small wavelength instabi
lities. Even though the short wavelength instabilities are
damped and the growth rate is bounded, the numerical
simulations for the slug flow case II including surface
tension did not give promising results. For actual values
of the surface tension coefficients the computed slug
characteristics still showed a strong dependence on the
grid used. This is a consequence of the unrealistically
high growth rates displayed in figure 4 as well as the
absence of a plateau in the growth rate curve.
As suggested by Ramshaw & Trapp (1978) and
Holmas et al. (2008) the system of equations with
surface tension is unquestionably more physical than
the standard model, given that surface tension is a real
phenomenon. However, even though the surface tension
has a stabilising effect it does not dissipate energy and
hence it does not damp instabilities. Once the amplitude
of small instabilities grows beyond the limit of linearity,
their energy can be transferred to other wavelengths.
Even longer wavelength perturbations, which would be
naturally stable but not damped, will grow indefinitely,
not because of their natural instability but due to the
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
additional energy received from shorter wavelengths,
leading to an obviously unphysical instability of the
whole system. Instead, in the case of viscosity the short
wavelength instabilities are actually damped and their
energy is dissipated. Therefore, this seems to be a more
appropriate solution to the problem.
The idea behind the mass and momentum diffusion
terms is that they may not represent any specific phy
sical phenomenon, but are meant to suppress the short
wavelength instabilities which are captured by the model
despite not being physical.
On one hand, the behaviour of the growth rate of the
perturbations as a function of the wavelength must re
flect that of the wellposed case I, where the short wave
length instabilities are damped (negative growth rate)
and the growth rate is bounded by a certain value. On
the other hand, a clear evidence of the wellposedness
of the system emerges from grid independent results as
already shown in figure 7. The ideal set of coefficients
is able to render the system wellposed suppressing
the unphysical short wavelength instabilities without
affecting the stability of long wavelength perturbations
and hence the natural stability of the flow.
In case the amount of diffusion introduced in the sys
tem is not sufficient to damp these artificial instabilities,
the system will behave like an illposed system, in which
the computed results show a strong dependence on the
grid used. In case too much diffusion is introduced
in the equations, the system would behave like a well
posed case converging to a unique solution regardless
of the mesh. However, the slug characteristics (such
as frequency) predicted by the numerical simulations
would be in bad agreement with experiments.
Finally, if the results of the numerical simulations
converge to the same solution for different mesh re
finements and this solution is in good agreement with
experiments, the perfect amount of diffusion has been
introduced in the system. The stability curve for the
aforementioned optimal combination of coefficients is
graphically represented in both figure 6(a) and 6(b) by
the thick black line. Here the length of the plateau at
peak value does not play an important role as long as
the short wavelength instabilities are damped. Indeed
the cutoff wavelength in all cases is far beyond the
numerical capabilities of the grid used.
A sensitivity analysis aimed at finding the best combi
nation for the diffusion coefficients for case II, indicated
that values of F 0.1 and p 105 are optimal to
yield good solutions. This was based on both theoretical
findings from the linear stability analysis and results
obtained from numerical simulations.
The computed slug frequency for case II with and
without additional diffusion is plotted in figure 8. The
results for the computations with the standard model
0.5
0.4
>
0.3
0.2
on
O case II: standard model
O case II: axial diffusion
 experimental value
0 0.2 0.4 0.6
dimensionless mesh spacing, Ax/D []
Figure 8: Computed slug frequency against mesh
density with artificial diffusion for case II
are the same shown earlier in figure 7, where the slug
frequency increases as the mesh is refined. Introducing
the most favourable amount of diffusion, the numerical
computations using different grids converge to the same
value, which is consistent with experimental evidence.
This is a clear manifestation of the wellposedness of
the system and a demonstration that the introduction of
axial diffusion is able to render the system of equations
wellposed while retaining the ability to capture natural
hydrodynamic instabilities leading to slug flow.
Conclusions
The onedimensional twofluid model for twophase
flow in pipelines is based on area averaging of the
equations that results in a loss of information. This
must be compensated by the introduction of closure
relationships, which aim to better represent the physics
of the problem.
In this paper a mathematical and numerical analysis
of the twofluid model including the effects of surface
tension and axial diffusion is presented.
It is shown that for such closures the characteristics
analysis results are unreliable, because of the presence
of a mixture of first and high order derivative terms.
However, it is still possible to rely on the results of a
linear KelvinHelmholtz stability analysis. A necessary
condition for a system to be wellposed is that very
short wavelength perturbations must be damped and the
growth rate of the instabilities must be bounded by a
certain value.
The introduction of the model for the surface tension
is able to satisfy these conditions. However, surface
tension effects are negligible in comparison to gravity
in the range of wavelength which can be captured by a
numerical simulation.
On the other hand, artificial viscosity seems to be a
more appropriate effect to dissipate and damp unphys
ical short wavelength perturbation. Indeed by adding
artificial diffusion to the standard twofluid model it
is possible to obtain a wellposed system of equations.
This has been shown by looking at both theoretical
predictions from a linear stability analysis and results
of numerical simulations. The former shows that in
presence of diffusion short wavelength perturbations are
damped and the growth rate of the entire range of growth
rates is bounded. The latter show that the computed slug
characteristics converge to a unique solution regardless
of the mesh used which is a clear manifestation of the
wellposedness of the system. Moreover, the numerical
results are also in good agreement with experimental
data.
Unfortunately, these remarkable results have been
achieved for a single case through a sensitivity analysis
on the diffusion coefficients and knowing a priori the
target value. Therefore, on the basis of these very
promising results additional effort is required in order
to achieve a more general formulation of the diffusion
coefficients as a function of the flow conditions.
Acknowledgements
This work has been undertaken within the Joint Project
on Transient Multiphase Flows and Flow Assurance
TMF4. The authors wish to acknowledge the contri
butions made to this project by the UK Engineering
and Physical Sciences Research Council (EPSRC) and
the following: Advantica; BP Exploration; CDadapco;
Chevron; ConocoPhillips; ENI; ExxonMobil; FEESA;
IFP; Institutt for Energiteknikk; PDVSA (INTEVEP);
Petrobras; PETRONAS; ScandpowerPT; Shell; SIN
TEF; StatoilHydro and TOTAL. The authors wish to
express their sincere gratitude for this support.
References
Arai M. Characteristics and stability analysis for two
phase flow equation system with viscous term. Nucl. Sci.
Eng., Vol. 74(2), pp. 7783 (1980)
Barnea, D. and Taitel, Y. KelvinHelmholtz stability
criteria for stratified flow: viscous versus nonviscous
(inviscid) approaches. Int. J. Multiphase Flow, Vol.
19(4), pp. 639649 (1993)
Bonizzi, M. & Issa, R.I. A model for simulating gas
bubble entrainment in twophase horizontal slug flow.
Int. J. Multiphase Flow, Vol. 29(11), pp. 16851717
(2003)
Holmas, H., Sira, T., Nordsveen, M., Langtangen, H.P.
& Schulkes, R. Analysis of a 1D incompressible two
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
fluid model including artificial diffusion. IMA J Appl
Math, Vol. 1, pp. 117 (2008)
Ishii, M. & Hibiki, T. ThermoFluid Dynamics of Two
Phase Flow. Springer Science, New York (2006)
Issa, R.I. & Abrishami, Y. Computer modelling of
slugging flow. Technical report, Petroleum Eng. Dept.,
Imperial College London, July (1986)
Issa, R.I. & Kempf, M.H.W. Simulation of slug flow in
horizontal and nearly horizontal pipes with the twofluid
model. Int. J. Multiphase Flow, Vol. 29(1), pp. 6995
(2003)
Lin, P.Y. & Hanratty, T.J. Prediction of the initiation of
slugs with linear stability theory. Int. J. Multiphase Flow,
Vol. 12(1), pp. 7998 (1986)
Prosperetti, A. & Jones, A.V. The linear stability of
general twophase flow models. Int. J. Multiphase Flow,
Vol. 13(2), pp. 161171 (1987)
Prosperetti, A. & Tryggvason G. Computational meth
ods for multiphase flow, Cambridge University Press,
London (2007)
Ramshaw, J.D & Trapp, J.A. Characteristics, stability
and short wavelength phenomena in twophase equation
system. Nucl. Sci. Eng., Vol. 66(1), pp. 93102 (1978)
Spedding, P.L. & Hand N.P. Prediction in stratified gas
liquid cocurrent flow in horizontal pipelines. Int. J. Heat
Mass Transfer, Vol. 40(8), pp. 19231935 (1997)
Taitel, Y. & Dukler A.E. A model for prediction flow
regime transition in horizontal and near horizontal gas
liquid flow. AIChE J., Vol. 22(1), pp. 4755 (1976)
Wallis, G.B. Onedimensional twophase flow.
McGrawHill (1969)
