Method for modelling multiphase flows in pipelines

A unified hydraulic model has been developed by the method according to the invention which is applicable to any slope and diameter of pipeline and can handle most of the steady state as well as transient multiphase flow regimes encountered in practice. The new modelling method differentiates two types of flow patterns: separated flow patterns (stratified or annular) and dispersed flow patterns. Intermittent flow patterns (slug, churn flow) are a combination of these two patterns. The same concept has been successfully applied for transition criteria between different flow regimes, insuring continuity of the solutions across the transitions. This requirement is very important for simulating transient phenomena. The transient resolution is achieved by an explicit time advancing scheme. The advantages of the method are; its ability to follow wave front propagation, an easy implementation for the resolution of complex pipeline networks. The performance of the resulting unified hydraulic model is demonstrated using a large number of experimental data.

BACKGROUND OF THE INVENTION 
1. Field of the Invention 
The invention relates to a method for modelling steady state as well as 
transient multiphase flows such as hydrocarbon mixtures circulating in 
pipeline networks taking into account a set of variables defining fluid 
properties and flow patterns, as well as dimensions or slope angles of the 
pipelines. 
2. Description of the Prior Art 
Prior art relating to multi-phase 
Fabre, J., et al 1983. Intermittent gas-liquid flow in horizontal or 
slightly inclined pipes, Int. Conference on the Physical Modelling of 
Multi-Phase Flow, Coventry, England, pp 233, 254 
Fabre, J., et al 1989. Two fluid/two flow pattern model for transient gas 
liquid flow in pipes, Int. Conference on Multi-Phase Flow, Nice, France, 
pp 269, 284, Cranfield, BHRA. 
SUMMARY OF THE INVENTION 
The method according to the invention includes numerical resolution of a 
drift-flux type model applied to any situation encountered in multiphase 
production. The state of a multi-phase flow is determined by solving a set 
of four conservation or transport equations. Two equations relate to the 
mass of each phase in a mixture, a second one, to the mixture momentum and 
a third one, for the mixture energy. The missing information about the 
slip between phases is restored by a steady state closure model depending 
on the flow regime. 
In view of identifying the flow regimes, it is considered that each of them 
is a space vs time combination of two basic patterns: separated flow i.e 
flows stratified in a vertical direction or flows radially separated in a 
pipeline, and dispersed flows. Any intermediate or intermittent flow 
corresponding to liquid or gas slugs is then considered as a combination 
of the two basic patterns and is characterized by a parameter .beta. 
representing the fraction of separated flow which insures continuity of 
some closure relations. 
To represent closely the actual flowing of hydrocarbons in pipelines it is 
neccessary to build a model which is continuous across flow regime 
transitions and also to select a limited set of closure relations well 
qualified against experimental results and continuous with respect to pipe 
slope and fluid properties. 
The method according to the invention is used for modelling steady state 
and transient multiphase flows of fluids such as hydrocarbon mixtures 
circulating in pipelines taking into account a set of variables defining 
fluid properties and flow patterns and also dimensions or slope angles of 
the pipelines, by numerical resolution of a drift-flux type model, 
comprising solving a set of transport equations including mass 
conservation equations for each phase, a mixture momentum equation and a 
mixture energy equation, and closure relations suitable for dispersed flow 
regimes, a separated flow regime and any intermediate intermittent flow 
regime with formation of liquid or gas slugs, said method also comprising: 
characterizing flow regimes by a parameter .beta. representing the fraction 
of a flow in a separated state, said parameter .beta. continuously ranging 
from 0 for dispersed flow regimes and 1 for separated flow regimes; 
determining any current flow regime while determining said set of transport 
equations by comparing current values of a liquid fraction in slugs with 
respect to a liquid fraction in a dispersed region of the flow, as well as 
gas slug velocity with respect to a critical velocity; and 
imposing while solving said closure relations continuity constraints at the 
boundaries between said regimes to the respective gas volume fractions and 
to the slug velocities. 
The method can be carried out for example by imposing as an upper limit for 
the volume fraction of gas in the dispersed region of a flow, the gas 
volume fraction, and having the slug velocity continuously tend to the gas 
velocity in the slugs as said volume fraction of gas reaches said upper 
limit, thereby smoothing transitions at the dispersed/intermittent 
boundary. 
The method can include using a explicit time advancing scheme. 
The method according to the invention allows building a model provided with 
closure relations better accounting for transitions between intermediate 
intermittent flow regimes and the two basic patterns. It also offers a 
faster way to reach sought solutions by a better arrangment of the 
different intermediate steps. An early accounting for values of some 
parameters allows for determination of the true flow regime and thus 
avoids useless computing steps which are otherwise all achieved when tests 
relating to the values of such parameters are as usual effected afterwards 
.

DESCRIPTION OF THE PREFERRED EMBODIMENTS 
The method for modelling liquid/gas flow of fluids in pipelines includes 
solving a set of transport and closure equations while giving particular 
attention to providing realistic smooth solutions for flow regime 
transitions. The method is based upon: 
a selection of a set of transport equations which insures a good 
representation of the properties of the fluid progressing along pipes, 
the choice of a limited set of closure relations well qualified against 
experimental results and which has been made continuous with respect to 
slope and fluid properties, and 
the concept of flow pattern transitions based on the continuity of the 
calculated variables. 
The set of transport equations includes one mass conservation equation for 
each of the phases (Equations. 1 & 2), one mixture momentum equation 
(Equation 3) and one mixture energy equation (Equation.4). 
##EQU1## 
Subscript k stands for the phase (G for the gas phase and L for the liquid 
phase), the superscripts S, D for the two basic regimes (Separated flow 
and Dispersed flow). r, R, U and H are the density, the volume fraction, 
the velocity and the enthalpy respectively. q is the inclination angle 
with respect to horizontal. T.sup.w and Q.sup.w are the contributions of 
the wall friction and heat transfer across the wall. .beta. is a factor 
representing the fraction of a flow in a separated state. It ranges from 0 
for the flows totally in a dispersed state to 1 for the flows totally in 
the separated state. 
Five main variables are unknown: 
R.sub.G the gas volume fraction, 
P the pressure, 
T the temperature, 
U.sub.G the gas velocity, 
U.sub.L the liquid velocity. 
Thus, an additional equation is required to close this system of 4 partial 
differential equations: 
EQU Hydro[U.sub.G,U.sub.L,R.sub.G,P,T]=0 (5) 
Equation (5) reduces to a slip equation in the case of dispersed flow. It 
is a macroscopic momentum balance in the case of a stratified flow and a 
set of algebraic equations in the case of a intermittent flow. Thus, it is 
flow regime dependent, and as explained later special treatment is used to 
avoid any discontinuities when flow regime transitions occur. 
In addition, two complementary terms M.sub.c and E.sub.c appear in the 
momentum and energy equations. These terms vanish for dispersed flow 
(.beta.=0), and for separated flow (.beta.=1); in intermittent flow, they 
account for the non homogeneous distribution of a void (or gas) and 
velocities in the separated and dispersed parts. They are expressed versus 
the secondary variables U.sub.k.sup.S, U.sub.k.sup.D, R.sub.k.sup.D and 
.beta.. Moreover, the wall friction and heat fluxes may be written: 
EQU T.sup.w =.beta.T.sup.wS +(1-.beta.)T.sup.wD (6) 
EQU Q.sup.w =.beta.Q.sup.wS +(1-.beta.)Q.sup.wD (7) 
in which the contributions of separated and dispersed parts appear. 
A unified formulation is used for the gas and liquid wall friction and the 
interfacial friction and thereby the form of these closure relations is 
independent of the flow regime and has the same structure as a single 
phase flow relationship. Furthermore, it will be shown that it is possible 
to obtain continuous solutions of the hydrodynamic model across flow 
regime transitions, even though the analytical form of the hydrodynamic 
model for the two basic regimes is very different. 
Hereafter are discussed the closure relations and the constraints which are 
added to provide continuity at the boundaries between the general case 
(0&lt;.beta.&lt;1) corresponding to intermittency, and the two degenerated cases 
.beta.=0 (dispersed flow) and .beta.=1 (separated flow). 
Dispersed flow (.beta.=0) 
The closure equation (5) reduces to a drift flux relation expressing the 
gas velocity versus the mixture velocity U.sup.D and the bubble drift 
velocity V.sub.B. 
EQU U.sub.G.sup.D =C.sub.0 U.sup.D +V.sub.B (d.sub.B,R.sub.G.sup.D,.theta.)(8) 
The resolution algorithm for the dispersed flow regime is given in FIG. 2. 
Separated flow (.beta.=1) 
For separated flow, Eq. (5) is a combination of the gas and liquid momentum 
equations from which the pressure gradient has been eliminated: 
##EQU2## 
FIG. 3 shows the algorithm for obtaining a iterative solution for equation 
(10) from expressions of the wall and interfacial friction, 
T.sub.L.sup.wS, T.sub.G.sup.wS, T.sup.i. 
Intermittent flow (0&lt;.beta.&lt;1) 
In intermittent flowing, we have a periodic structure of dispersed and 
separated flows. 
The solution of the model depends on an accurate prediction of the slug 
velocity. So, we have chosen to express the slug velocity by the relation: 
##EQU3## 
where the coefficients C.sub.0 and C.sub..infin. are a function of the 
inclination angle, the Bond number Bo, Reynolds number Re and Froude 
number Fr defined respectively as: 
##EQU4## 
The solution algorithm for intermittent flow is given in FIG. 4. Flow 
regime transitions 
The algorithms as shown on FIGS. 2-4 correspond to a calculation with 
imposed flow regime. Hereafter we show how transition criteria lead to 
discontinuities in the calculated variables, why this is incompatible with 
the averaged nature of the equations of motion, and how the transition 
rules of the method according to the invention have been selected to avoid 
such discontinuities. 
A standard approach to determine flow regimes is to apply a set of 
transition criteria. However, if the criteria are not closely coupled with 
the models for the individual flow regimes, they lead to discontinuities 
in the calculated variables. This is best illustrated in FIG. 5 where the 
hydrodynamic models (Hydro 1 and Hydro 2) for flow regimes 1 and 2 are 
represented. This graph schematically represents the dependence of one of 
the variables (e.g. the gas volume fraction R.sub.G) as a function of the 
other dependent variables (e.g. the superficial liquid velocity U.sub.LS). 
FIG. 5 then shows that unless the transition criterion crosses the 
intersection between the two hydrodynamic models, a discontinuity in one 
of the dependent variables will occur at the transition. 
The averaged equations of motion are based on the assumption that the 
length scale of the flow is small with respect to the length of the 
control volume over which the equations are averaged. Similarly, the flow 
regime must be understood as an averaged property of the flow structure in 
the control volume. Thus, even if locally, the void fraction and slip 
velocity between phases may present a discontinuous behaviour, these 
variables cannot be discontinuous in an averaged sense. 
Proposed treatment of flow regime transitions 
The method according to the invention then provides special treatment of 
flow discontinuities. The constraint that the calculated variables be 
continuous across the transitions, implies that the hydrodynamic models 
representing the different flow regimes, though very different in form, 
should lead to continuous solutions across transition boundaries. To 
obtain this, the following rules are implemented: 
The intermittent flow regime is considered as the basic flow regime. The 
parameter .beta.: 
##EQU5## 
ranges from 0 to 1, and is a result of the calculation. 
In relation 12: 
R.sub.G is the gas volume fraction; 
R.sub.G U.sub.G, the gas superficial velocity; and 
superscript P and B refer respectively to the separated and dispersed 
regions. 
When .beta.=0., (R.sub.G.sup.B =R.sub.G) the hydrodynamic model H reduces 
to the dispersed flow model, that is a slip equation in the liquid slug. 
The transition between dispersed and intermittent flow occurs when the 
void fraction in the liquid slug is equal to the void fraction of pure 
dispersed flow. 
When .beta.=1, (R.sub.G.sup.P =R.sub.G) the hydrodynamic model reduces to 
the stratified flow model, i.e a momentum balance in the separated region. 
It must be insured however that: 
no solution can be found for .beta. outside of the interval [0,1], and 
the solution tends smoothly towards these limits. 
The dispersed/intermittent boundary (.beta.=0) is dealt with by imposing 
the constraint R.sub.G.sup.B .ltoreq.R.sub.G to the constitutive relation 
for the holdup in the slug. Furthermore a smooth transition is guaranteed 
by having the slug velocity tend continuously towards U.sub.G.sup.B i.e. 
the gas velocity in the slug as R.sub.G.sup.B .fwdarw.R.sub.G : 
For the stratified/intermittent boundary the following general relations 
are used: 
##EQU6## 
In intermittent flow, the gas bubbles in slugs coalesce with the front of 
the gas pocket, and bubbles are released at the gas pocket tail into the 
following liquid slug. Thus: 
EQU V.sub.P .gtoreq.U.sub.G. 
Let us consider now the limiting case when the intermittent configuration 
tends towards the separated configuration. We have R.sub.G 
.ltoreq.R.sub.G.sup.P and R.sub.G .fwdarw.R.sub.G.sup.P 
.fwdarw.R.sub.G.sup.S. It turns out that the limiting value of the slug 
velocity when .beta..fwdarw.1 is U.sub.G : By combining (12) and (14), we 
obtain: 
##EQU7## 
and as .beta..fwdarw.1, U.sub.G.sup.P .fwdarw.U.sub.G, and V.sub.P 
.fwdarw.U.sub.G. 
Though in principle this result is true, in practice, the slug velocity is 
represented by a closure law or constitutive relation and V.sub.P &gt;U.sub.G 
irrespective of the value of .beta.. To force this limit, (FIG. 7A) we put 
an upper bound for the slug velocity defined as: 
##EQU8## 
It is easy to verify that when .beta.=1, R.sub.G.sup.S =R.sub.G.sup.P and 
U.sub.G R.sub.G =U.sub.G.sup.P R.sub.G.sup.P, thus V.sub.p,strat =V.sub.P 
=U.sub.G. When .beta.&lt;1 it is assumed that the slip in separated flow 
(FIG. 7B) is greater than in intermittent flow, thus R.sub.G.sup.P 
&gt;R.sub.G &gt;R.sub.G.sup.S, and from relations (14) and (16), V.sub.P 
&lt;V.sub.p,strat. 
It is also important to note that the determination of the actual flow 
regime can take place in the method of the invention at an early stage of 
the general algorithm. As can be seen on FIG. 1, testing may be made on 
the value of the liquid proportion RL in the slugs with respect to the 
corresponding proportion in the dispersed regime. In the same way, testing 
may also be made on the gas slug velocity V.sub.p with respect to Vmax. As 
a result the proper flow regime can be determined and only the computation 
of this actual flow regime is achieved. 
Numerical scheme 
The model requires the solution of a hyperbolic set of partial differential 
equations of the form given in equations (1) to (4): 
##EQU9## 
together with the algebraic constraint: 
EQU Hydro(U)=0 (18) 
U=T(R.sub.G,U.sub.G,U.sub.L,P,T) is the set of dependent variables. At each 
time step, W is computed by a time advancing scheme, and the dependent 
variables are calculated with the help of the hydrodynamic model (Eq. 19). 
The time advancing scheme is explicit, so that the state of the system at 
time (n+1) is calculated on the basis of the flux term F and source terms 
Q evaluated at time n. 
Explicit schemes are stable provided that the Courant Friedriech Levy 
criterion well known in the art is satisfied. In practice, this condition 
imposes a limitation on the time step .DELTA.t depending on the speed of 
the fastest wave propagating in the system c.sub.max : 
##EQU10## 
The method according to the invention comprises use of a three point 
predictor corrector scheme which is an explicit time advancing scheme with 
a second order accuracy as described by: 
Lerat, A., 1981, "Sur le calcul des solutions faibles des systemes 
hyperboliques de lois de conservation a l'aide de schemas aux 
differences." Thesis, ONERA, France. 
A time advancing scheme was chosen because of the following arguments: 
It has the advantage of a clear separation between the hydrodynamic 
calculation module, namely the function Hydro H and the treatment of the 
numerical scheme. This allows an independent evolution of the physical 
model and the numerical scheme. Thus such explicit scheme is easy to 
implement. 
For transient simulation of complex networks, it allows for the solution of 
individual network components sequentially instead of requiring their 
simultaneous solution. Thus, the selected scheme gives more flexibility 
for the treatment of complex transient networks. 
A good accuracy for the description of transients is needed in order to 
predict precisely the volume of liquid slugs travelling in the pipeline 
system and reaching the pipeline outlet. The numerical scheme gives second 
order accuracy in time and space. It also allows good front tracking 
capabilities. This will be very important when estimating the size of 
large liquid slugs reaching the pipeline outlet. 
FIGS. 8A and 8B show a typical test case on of horizontal pipe 10 km long. 
The pressure at the outlet is decreased in 100 seconds from 10 bars to 5 
bars. FIG. 7A shows that the pressure wave propagates very fast towards 
the inlet to generate a void fraction front which propagates to the outlet 
in 2500 seconds. FIGS. 8A and 8B illustrate the importance of having a 
numerical scheme able to track accurately wave fronts propagating in the 
system. 
Special attention has been drawn to the treatment of boundary conditions. 
The hyperbolic nature of the problem makes it very important to treat 
precisely the transport of information through the pipe, and in particular 
the information going in and out of the system. Four types of waves are 
identified in the analytical model: 
Void fraction waves travelling downstream. 
Pressure waves travelling downstream. 
Pressure waves travelling upstream. 
Enthalpy waves travelling with the mean flow. 
The applied treatment of boundary conditions takes into account the 
information going out of the system through a set of compatibility 
relations, which must be solved together with the imposed boundary 
conditions. Thus, it depends on the direction of propagation of the 
different wave systems. Under such conditions, the solution of the 
following set of equations gives a rigorous treatment of the boundary 
conditions: 
At the inlet: 
Mass flowrate of gas 
Mass flowrate of liquid 
Temperature 
Compatibility condition for pressure waves travelling upstream 
Hydrodynamic function 
At the outlet: 
Pressure 
Compatibility condition for enthalpy propagation 
Compatibility condition for pressure waves travelling downstream 
Compatibility condition for void fraction waves 
Hydrodynamic function 
Validation of the hydrodynamic function 
The hydrodynamic module of the method according to the invention has been 
validated extensively and notably on experimental data available from 
measurements on a test pipe loop. The data cover a wide range of 
inclinations (-3.degree., -0.5.degree., 0.degree., 1.degree., 4.degree., 
15.degree., 45.degree., 75.degree. and 90.degree.) with 3" and 6" piping. 
diesel oil and condensate have been used alternatively together with 
natural gas as working fluids. The fluid properties, which were measured 
in the laboratory pressures up to 5 MPa, were investigated over a wide 
range of flowrates. the horizontal and slightly inclined data bank 
contains about 1750 points of which 900 are in stratified flow, 600 are in 
slug flow, and 250 are in dispersed flow. The data bank on vertical and 
highly inclined flow consists of about 700 measurement points (260 for 
bubbly flow, 400 for slug flow and 40 for annular flow). FIGS. 8A and 8B 
for example show comparisons of the model predictions for inclinations of 
15.degree., 45.degree., 75.degree. and 90.degree. of the pipes. 
The hydrodynamic function has also been tested successfully on a number of 
operating North Sea fields