Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: 5.1.1 - The Effect of Surface Tension and Diffusion on One-dimensional Modelling of Slug Flow
ALL VOLUMES CITATION THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00102023/00118
 Material Information
Title: 5.1.1 - The Effect of Surface Tension and Diffusion on One-dimensional Modelling of Slug Flow Instabilities
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Issa, R.I.
Montini, M.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
 Subjects
Subject: slug flow
two-fluid model
surface tension
diffusion
well-posedness
 Notes
Abstract: This paper investigates the use of additional closure relationships, such as surface tension and diffusion, with the aim of obtaining a well-posed 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 well-posed. However, more work is needed to determine a universal relationship for the diffusion coefficients necessary to render the system well-posed for all the cases in the slug flow regime.
General Note: The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF-2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: Bio-Fluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and Nano-Scale Multiphase Flows; Microgravity in Two-Phase Flow; Multiphase Flows with Heat and Mass Transfer; Non-Newtonian Multiphase Flows; Particle-Laden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows
 Record Information
Bibliographic ID: UF00102023
Volume ID: VID00118
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: 511-Issa-ICMF2010.pdf

Full Text



7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30- June 4, 2010


The effect of Surface Tension and Diffusion on One-dimensional
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, Two-fluid model, Surface tension, Diffusion, Well-posedness




Abstract

This paper investigates the use of additional closure relationships, such as surface tension and diffusion, with the
aim of obtaining a well-posed 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 well-posed. However, more work is needed to determine a universal
relationship for the diffusion coefficients necessary to render the system well-posed 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 (kgs-lim1)
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 one-dimensional two-fluid 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 three-dimensional nature of two-phase
flow, and in particular of phenomena such as turbulence
and the presence of slugs, the one-dimensional 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 one-dimensional two-fluid model is ill-posed 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 well-posedness 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 well-posedness 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 cross-sectional 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 gas-liquid 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 Blasius-like 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 liquid-wall 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 Well-posedness
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 gas-wall 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 well-posedness 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 ill-posed 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 well-posedness 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
gas-liquid 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 well-posed only if the
four roots are real and therefore it is possible to find
the following criterion for the well-posedness 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 Well-posedness limit ..
....... ... ................... 0 Ill-posed point
[ Well-posed 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 well-posedness 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 well-posed. Conversely, in
the area above the line the system is characterized by its
non-hyperbolicity and case II symbolised by a circle on
the map is, therefore, ill-posed.
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 well-posedness 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 well-posedness 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 well-posedness 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 Kelvin-Helmholtz
(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 wave-number 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(ak-ib+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)


S-where
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: ill-posed
- case I: well-posed







1 1 I """I 1 1
---------------


le-06
le-06


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 well-posed 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 well-posedness of the problem. On the
other hand for an ill-posed 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 well-posedness 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 well-posed 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 cut-off 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 cut-off
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 cut-off wavelength beyond which
the perturbations are stable: the higher the value of p
the longer the cut-off wavelength.
Once the stability analysis has provided an indication
regarding the behaviour of the system, a confirmation of
the well-posedness 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 ill-posed 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 well-posed 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 -
le-06


- 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 time-averaged slug frequency is plotted against
the size of the grid spacing: for a well-posed 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 ill-posed, 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 ill-posed case there


0.6

0.5

0.4

S0.3
C-

0.2

0.1

0
0


o case II: ill-posed
0 case I: well-posed
- 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 two-fluid 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 well-posed 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 well-posedness
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 well-posed 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 ill-posed 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 cut-off 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 well-posedness of
the system and a demonstration that the introduction of
axial diffusion is able to render the system of equations
well-posed while retaining the ability to capture natural
hydrodynamic instabilities leading to slug flow.

Conclusions

The one-dimensional two-fluid model for two-phase
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 two-fluid 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 Kelvin-Helmholtz stability analysis. A necessary
condition for a system to be well-posed 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 two-fluid model it
is possible to obtain a well-posed 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
well-posedness 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; CD-adapco;
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. 77-83 (1980)

Barnea, D. and Taitel, Y. Kelvin-Helmholtz stability
criteria for stratified flow: viscous versus non-viscous
(inviscid) approaches. Int. J. Multiphase Flow, Vol.
19(4), pp. 639-649 (1993)

Bonizzi, M. & Issa, R.I. A model for simulating gas
bubble entrainment in two-phase horizontal slug flow.
Int. J. Multiphase Flow, Vol. 29(11), pp. 1685-1717
(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. 1-17 (2008)

Ishii, M. & Hibiki, T. Thermo-Fluid 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 two-fluid
model. Int. J. Multiphase Flow, Vol. 29(1), pp. 69-95
(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. 79-98 (1986)

Prosperetti, A. & Jones, A.V. The linear stability of
general two-phase flow models. Int. J. Multiphase Flow,
Vol. 13(2), pp. 161-171 (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 two-phase equation
system. Nucl. Sci. Eng., Vol. 66(1), pp. 93-102 (1978)

Spedding, P.L. & Hand N.P. Prediction in stratified gas-
liquid co-current flow in horizontal pipelines. Int. J. Heat
Mass Transfer, Vol. 40(8), pp. 1923-1935 (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. 47-55 (1976)

Wallis, G.B. One-dimensional two-phase flow.
McGraw-Hill (1969)




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - Version 2.9.7 - mvs