Skip to main content

A new drift-flux model for all gas volume fractions

This page contains a description of a revised drift-flux model for all gas volume fractions#

Printable version#

Oliasoft Technical Documentation - Drift-flux Model

1 Introduction#

The classical drift-flux model for two-phase gas-liquid pipe flow describes the slip between the gas and liquid phases as the combined effect of non-uniform distribution of gas and liquid across the pipe cross section and the additional effect of gas buoyancy as well as local hydraulic gradients near the tip of a long bubble (Taylor bubble).

The main assumption in the model is that these two effects are additive:

UG=C0UM+UdU_G = C_0U_M + U_d

Where UGU_G is the gas velocity and UMU_M the mixture velocity given by the sum of the gas and liquid superficial velocities


The distribution coefficient C0C_0 and the drift velocity UdU_d are often taken to be flow regime dependent, since the physical phenomena they reflect vary qidely with the spatial distribution of gas and liquid.

Note that continuity implies that

USG=Ξ±GUGU_{SG} = \alpha_GU_G

USL=(1βˆ’Ξ±G)ULU_{SL} = (1-\alpha_G)U_L

UM=Ξ±GUG+(1βˆ’Ξ±G)ULU_M = \alpha_GU_G + (1-\alpha_G)U_L

From the above we can infer

UL=UMβˆ’Ξ±GUG1βˆ’Ξ±G=1βˆ’Ξ±GC01βˆ’Ξ±GUMβˆ’Ξ±1βˆ’Ξ±GUdU_L = \frac{U_M-\alpha_GU_G}{1-\alpha_G} = \frac{1-\alpha_GC_0}{1-\alpha_G}U_M - \frac{\alpha}{1-\alpha_G}U_d

We also have

Ξ±G=USGC0UM+Ud\alpha_G = \frac{U_{SG}}{C_0U_M+U_d}

The classical-type drift-flux model of Bhagwat and Ghajar (2014) has been employed for two-phase gas-liquid flow, since this model has been developed to cover all pipe or well inclinations and thus seems to be the most general and up to data drift-flux model available. We have then added the three-phase corrections from Shi et al (2003) on top of the basic two-phase gas-liquid drift-flux model.

2 The standard drift-flux model for two-phase flow#

Here we will first go through the two-phase drift-flux model of Bhagwat and Ghajar (2014), which is valid for all inclinations, and then add the three-phase corrections from Shi et al (2003) in the following section.

Distribution coefficient for gas-liquid slip#

Bhagwat and Ghajar correlated the distribution parameter as follows:

C0=2βˆ’(ρGρL)21+(Re1000)2+[(1+(ρGρL)2cosΞΈ)/(1+cosΞΈ)](1βˆ’Ξ±5)+C0,11+(1000Re)2C_0 = \frac{2-(\frac{\rho_G}{\rho_L})^2}{1+(\frac{Re}{1000})^2} + \frac{[(1+(\frac{\rho_G}{\rho_L})^2cos\theta)/(1+cos\theta)]^{(\frac{1-\alpha}{5})}+C_{0,1}}{1+(\frac{1000}{Re})^2}

Where ΞΈ\theta is the pipe inclination from the horizontal, and the mixture Reynolds number is given by

Re=UMρLDH/μLRe = U_M\rho_LD_H/\mu_L

The hydraulic diamter DHD_H is taken as the pipe diameter in pipe flow. In annulus flow, the hydraulic diameter is taken as the difference between inner and outer diameter, DH=Douterβˆ’DinnerD_H = D_{outer} - D_{inner}

The distribution parameter for horizontal flow is given by

C0,1=C1(1βˆ’ΟGρL)[(2.6βˆ’Ξ²)0.15βˆ’ftp](1βˆ’x)3/2C_{0,1} = C_1(1-\sqrt{\frac{\rho_G}{\rho_L}})[(2.6-\beta)^{0.15}-\sqrt{f_{tp}}](1-x)^{3/2}

except for gravity dominated flow (βˆ’50°≀θ≀0Β°-50\degree \leq \theta \leq 0\degree and FrSG≀0.1Fr_{SG} \leq 0.1), where C0,1C_{0,1} is set to zero. The constant C1C_1 is 0.2 for circular and annular cross sections and 0.4 for rectangular cross sections.

x is the so-called two-phase flow quality, which is equal to the following gas mass fraction:

x=ρGUSGρGUSG+ρLUSLx = \frac{\rho_GU_{SG}}{\rho_GU_{SG} + \rho_LU_{SL}}

Ξ²\beta is the following gas volume fraction (GVF)

Ξ²=USGUSG+USL\beta = \frac{U_{SG}}{U_{SG}+U_{SL}}

The two-phase Fanning friction factor can then be computed using Churchill’s friction formula with the two-phase Reynolds number and hydraulic diameter defined above:

ftp=fChurchill(Re,Ο΅/DH)f_{tp} = f_{Churchill}(Re,\epsilon/D_H)

Remember that the Fanning friction factor f is 1/4 of the Darcy friction factor Ξ»\lambda.

Drift velocity for gas-liquid slip#

The drift velocity UdU_d (denoted U_{gm} in the paper of Bhagwat and Ghajar) is given by the following equations:

Ud=(0.35sinΞΈ+0.45cosΞΈ)(1βˆ’Ξ±G)gDH(1βˆ’ΟGρL)C2C3C4U_d = (0.35sin\theta+0.45cos\theta)\sqrt{(1-\alpha_G)gD_H(1-\frac{\rho_G}{\rho_L})}C_2C_3C_4

Where the coefficients C2C_2, C3C_3 and C4C_4 are given by

C2={(0.434log⁑10(ΞΌL0.001))0.15Β forΒ ΞΌL0.001>101Β forΒ ΞΌL0.001≀10C_{2}=\left\{\begin{array}{cc}\left(\frac{0.434}{\log _{10}\left(\frac{\mu_{L}}{0.001}\right)}\right)^{0.15} & \text { for } \frac{\mu_{L}}{0.001}>10 \\ 1 & \text { for } \frac{\mu_{L}}{0.001} \leq 10\end{array}\right.

C3={(La0.025)0.9Β forΒ La<0.0251Β forΒ Laβ‰₯0.025C_{3}=\left\{\begin{array}{cc}\left(\frac{L a}{0.025}\right)^{0.9} & \text { for } L a<0.025 \\ 1 & \text { for } \quad L a \geq 0.025\end{array}\right.

C4={βˆ’1Β forΒ βˆ’50βˆ˜β‰€ΞΈβ‰€0∘ andΒ FrSG≀0.11Β otherwiseΒ C_{4}= \begin{cases}-1 \text { for } & -50^{\circ} \leq \theta & \leq 0^{\circ} \text { and } F r_{S G} \leq 0.1 \\ 1 & \text { otherwise }\end{cases}

The Laplace number La is given by

La=Οƒg(ρLβˆ’ΟG)DH2L a=\sqrt{\frac{\sigma}{g\left(\rho_{\mathrm{L}}-\rho_{G}\right) D_{H}^{2}}}

Here FrSGFr_{SG} is the Froude number based on the gas phase

FrSG=Ug(1βˆ’ΟGρL)gDHF r_{S G}=\frac{U_{g}}{\sqrt{\left(1-\frac{\rho_{G}}{\rho_{L}}\right) g D_{H}}}

Drift flux model for three-phase flow#

In the original design, the three-phase drift-flux model of Shi et al (2003) was used. There are some limitations in this model regarding pipe inclination and other parameters that we want to avoid. We have therefore opted to combine the two-phase drift-flux model of Bhagwat and Ghajar with the oil/water slip model of Hasan and Kabir (1999) which is used by Shi et al. for their three-phase corrections. We take care to ensure that no inconsistencies arise from combining these two models.

From two-phase to three-phase#

For the three-phase case, we first use the two-phase model to compute the gas and liquid velocities assuming known gas and liquid volume fractions and mixture velocity (transient case). Since we only will be using the three-phase model for dynamic kill simulations, we only need to consider the transient case here. We can thus assume that the gas, oil and water holdups are known from the previous time step or iteration.

As a logical analogy to the drift-flux relation for gas and liquid, we assume a quasilinear relation between the oil velocity UOU_O and the volume average liquid velocity ULU_L:

UO=C0β€²UL+Udβ€²U_{O}=C_{0}^{\prime} U_{L}+U_{d}^{\prime}

Where we let the prime (') denote the distribution coefficient and drift velocity for the oil-water slip.

Once the oil veolcity has been computed, the water velocity is then given by the continuity.
This gives

UW=USLβˆ’Ξ±OU0Ξ±W=USLβˆ’Ξ±OUO1βˆ’Ξ±Gβˆ’Ξ±OU_{W}=\frac{U_{S L}-\alpha_{O} U_{0}}{\alpha_{W}}=\frac{U_{S L}-\alpha_{O} U_{O}}{1-\alpha_{G}-\alpha_{O}}

Distribution coefficient for oil-water slip#

The oil-water distribution coefficient is then given by

C0β€²={Aβ€²Β forΒ Ξ±o≀B1β€²1Β forΒ Ξ±0β‰₯B2β€²Aβ€²βˆ’(Aβ€²βˆ’1)Ξ±0βˆ’B1β€²B2β€²βˆ’B1β€²Β forΒ B1′≀α0≀B2β€²C_{0}^{\prime}=\left\{\begin{array}{cl}A^{\prime} & \text { for } \quad \alpha_{o} \leq B_{1}^{\prime} \\ 1 & \text { for } \quad \alpha_{0} \geq B_{2}^{\prime} \\ A^{\prime}-\left(A^{\prime}-1\right) \frac{\alpha_{0}-B_{1}^{\prime}}{B_{2}^{\prime}-B_{1}^{\prime}} & \text { for } \quad B_{1}^{\prime} \leq \alpha_{0} \leq B_{2}^{\prime}\end{array}\right.

Where the parameters Aβ€²A^{\prime}, B1β€²B_1^{\prime} and B2β€²B_2^{\prime} have the values 1.2, 0.4 and 0.7 respectively.

Drift velocity for oil-water slip#

The drift velocity for the oil-water slip is given by

Udβ€²=1.53VCβ€²(1βˆ’Ξ±O)2(cos⁑θ)0.5(1+sin⁑θ)2U_{d}^{\prime}=1.53 V_{C}^{\prime}\left(1-\alpha_{O}\right)^{2}(\cos \theta)^{0.5}(1+\sin \theta)^{2}

Where the characteristic velocity VCβ€²V_C^{\prime} is given by

VCβ€²=[(ρWβˆ’Οo)gΟƒOWρW2]1/4V_{C}^{\prime}=\left[\frac{\left(\rho_{W}-\rho_{o}\right) g \sigma_{O W}}{\rho_{W}^{2}}\right]^{1 / 4}

3 Improved/revised drift-flux model#

THe classical drift-flux model is invalid as the gas volume fraction (GVF) approaches unity, and as a result, the liquid velocity may approach ±∞\pm \infin , and when USLβ†’0U_{SL} \rightarrow 0 , it can be easily shown that the gas velocity approaches infinity already as Ξ±Gβ†’1/C0\alpha_G \rightarrow 1/C_0.

The deeper reason for this is that the drift-flux model is based on the physics of slug flow and/or bubbly flow, which can generally be classified as liquid dominated flow regimes (flow patterns). The form of the drift-flux equation is the same as the equation for the velocity of Taylor bubbles (long bubbles) in slug flow, and hence the behaviours of the two equations are similar under a wide range of conditions.

Since bubbly flow (or dispersed bubble flow) can be regarded as a special case of slug flow without long bubbles, the drift-flux equation can be adapted of bubbly flow by imposing a distribution parameter C0C_0 close to unity. Any dependence of the slip on pipe inclination and fluid properties will then enter through the drift velocity UdU_d.

To account for this limitation the classical drift flux model will be replaced by an alternative model for cases with a high GVF. In the classical drift-flux model, the gas velocity is expressed as a linear function of the mixture velocity. In the alternative model, the gas velocity is expressed as a linear function of the liquid velocity, allowing for greater flexibility of the gas and liquid velocities to adjust to each other under extreme conditions like for example very high gas volume fractions.

The normal situation in horizontal or upward flow, or in high velocity downward flow, is that the gas velocity is higher than the mixture velocity, which in turn is higher than the liquid velocity. Thus, the coefficients in the new, alternative (hybrid) drift-flux model will be different from those in the classical model, and/or the distribution coefficient C0C_0 must be higher in the new alternative (hybrid) model, since UG/ULU_G/U_L should be higher than UG/UMU_G/U_M in most cases.

It is conceivable that both the distribution coefficient C0C_0 and the drift velocity UdU_d might be higher in the new model than in the classical model; however, from functional analysis it may be argued that the constant term should be the same, since the velocity dependence is contained in the linear C0C_0 term. For the time being, we will assume that the drift velocity term will be carried over from the classical model to the new model.


The classical drift-flux model is very well-established and tested for liquid dominated flows, i.e. slug and bubbly flow. The classical drift-flux model was implemented in the kill module in Oliasoft’s software, and should be kept more or less unchanged for liquid dominated flow. The new approach is to employ the new drift-flux model for gas dominated flow, and then use an interpolation between the two models for intermediate gas volume fractions. From experience we have defined gas dominated flow as multiphase flow with an input gas volume fraction of more than 90%. This may sound like a very strict definition, but dynamically it makes sense since a GVF of 90% corresponds to a gas mass fraction of 50% at a gas density of about 90 kg/m3m^3 (depending on the liquid density).

4 The new (hybrid) drift-flux model#

Two-phase gas liquid flow#

The model equations for the new drift-flux models are written as follows:

UG=C1UL+UdU_{G}=C_{1} U_{L}+U_{d}

The distribution coefficient C1C_1 has been tuned against an experimental data base consisting of approximately 4800 data points for horizontal and moderately inclined flow, 2400 points for vertical flow, and 660 points for steeply inclined flow. C1C_1 is given by

C1=CHcos⁑2θ+CVsin⁑2θC_{1}=C_{H} \cos ^{2} \theta+C_{V} \sin ^{2} \theta


CH={1+kΞ²5(ρLρG)0.85max⁑(1,FrcFrd)cos⁑θ,ΞΈβ‰₯01,ΞΈ<0C_{H}=\left\{\begin{array}{cc}1+k \beta^{5}\left(\frac{\rho_{L}}{\rho_{G}}\right)^{0.85} \max \left(1, \frac{F r_{c}}{F r_{d}}\right) \cos \theta, & \theta \geq 0 \\ 1, & \theta<0\end{array}\right.


CV=1C_V = 1

Again, Ξ²\beta is the gas volume fraction (GVF) and ΞΈ\theta the inclination relative to horizontal. THe tuning constants are given by the coefficient k = 0.11 and critical Forude number FrcFr_c = 1.

We impose the upper limit CH≀5C_H \leq 5, otherwise CHC_H could become very large for high density ratios ρL/ρG\rho_L/\rho_G. The upper limit of 5 has been determined by tuning, but does not apply to atmospheric air/water data, which are frequently cited in the literature but irrelevant for our application.

The "densitometric# Froude number is given by

Frd=ρGUSG2(ρLβˆ’ΟG)gDcos⁑θF r_{d}=\frac{\rho_{G} U_{S G}^{2}}{\left(\rho_{L}-\rho_{G}\right) g D \cos \theta}

The justification for the distinction between upflow (ΞΈ>0\theta > 0) and downflow (ΞΈ<0\theta < 0) is that the slip in near horizontal gas-liquid flow tends to be much higher for small to moderate upward inclinations than for small to moderate downward inclinations, except for very high mixture velocities, where the flow can be viewed as independent of gravity.

The drift velocity is kept the same as for the classical drift-flux model, since this coefficient is independent of the gas and liquid superficial velocities in the most common cases. The drift velocity is consequently computed from the classical Bendiksen correlation:

Ud=(0.54cos⁑θ+0.35sin⁑θ)(1βˆ’ΟG/ρL)gDU_{d}=(0.54 \cos \theta+0.35 \sin \theta) \sqrt{\left(1-\rho_{G} / \rho_{L}\right) g D}

For downward inclined flow we have to make a correction, since the tip of the Taylor bubble will switch direction for low velocities:

ΞΈ<0\theta<0 and Fr<1:{C1=0.9Ud=(βˆ’0.54cos⁑θ+0.35sin⁑θ)(1βˆ’ΟG/ρL)gDF r<1:\left\{\begin{array}{c}C_{1}=0.9 \\ U_{d}=(-0.54 \cos \theta+0.35 \sin \theta) \sqrt{\left(1-\rho_{G} / \rho_{L}\right) g D}\end{array}\right.

Where Fr is a suitable Froude number, e.g.

Fr=UM/(1βˆ’ΟG/ρL)gDF r=U_{M} / \sqrt{\left(1-\rho_{G} / \rho_{L}\right) g D}

Now for the new drift flux equation, we can compute the void fraction (gas volume fraction) in the following way:

Define V=USG+C1USL+UdV=U_{S G}+C_{1} U_{S L}+U_{d}

Ξ±G={USGV,Ud=0Vβˆ’V2βˆ’4USGUd2Ud,Udβ‰ 0\alpha_{G}=\left\{\begin{array}{c}\frac{U_{S G}}{V}, U_{d}=0 \\ \frac{V-\sqrt{V^{2}-4 U_{S G} U_{d}}}{2 U_{d}}, U_{d} \neq 0\end{array}\right.

We can trivially compute Ξ±L=1βˆ’Ξ±G,UG=USGΞ±G,UL=USLΞ±L\alpha_{L}=1-\alpha_{G}, U_{G}=\frac{U_{S G}}{\alpha_{G}}, U_{L}=\frac{U_{S L}}{\alpha_{L}}.

Interpolation method#

For a smooth but fast transition between the liquid dominated (low GVF) and gas dominated (high GVF) regions, a weighting based on a hyperbolic tangent function is applied:

Ο‰=12[1βˆ’tanh⁑(15(Ξ²βˆ’0.85))]\omega=\frac{1}{2}[1-\tanh (15(\beta-0.85))]

C1=Ο‰C1,liq+(1βˆ’Ο‰)C1,gasC_{1}=\omega C_{1, l i q}+(1-\omega) C_{1, g a s}

Ud=Ο‰Ud,liq+(1βˆ’Ο‰)Ud,gasU_{d}=\omega U_{d, l i q}+(1-\omega) U_{d, g a s}

Where Ξ²\beta is the gas volume fraction and liq and gas denote the values for low GVF and high GVF respectively.

5 Implementation#

In this section we discuss some implementation issues.

The steady state case#

For the steady state case, we normally assume that the gas and liquid superficial velocities USGU_{SG} and USLU_{SL} are known and that the void fraction Ξ±G\alpha_G is unknown. In this case we can compute the unknown void fraction and phase velocities in the standard drift-flux model as described previously.

The dynamic case#

For the dynamic or transient case, we normally assume that the mixture velocity UMU_{M} and the void fraction Ξ±G\alpha_G or liquid holdup (1βˆ’Ξ±G)(1-\alpha_G) are known. In this case, when using the standard drift-flux model we have simply

UG=C0UM+UdU_{G}=C_{0} U_{M}+U_{d}

UL=UMβˆ’Ξ±GUG1βˆ’Ξ±GU_{L}=\frac{U_{M}-\alpha_{G} U_{G}}{1-\alpha_{G}}

While for the hybrid drift-flux model for high GVF, we have to apply the equations from the previous section in a slightly modified manner.


[1] Shi, H., Holmes, J. A., Durlofsky, L. J., Aziz, K., Diaz, L. R., Alkaya, B. and Oddie, G. (2003): Drift-Flux Modeling of Multiphase Flow in Wellbores. SPE 84228.

[2] Bhagwat, S. M. and Ghajar, A. J. (2014): A flow pattern independent drift flux model based on void fraction correlation for a wide range of gas-liquid two-phase flow. Int. J. of Multiphase Flow 59, pp. 186-205.