Skip to main content

The 1D Flow Computation


Single phase flow or multiphase flow is modeled, depending on the present fluids. Mechanistic models for multiphase flow are generally more accurate than the correlations, but the costs of implementing a full mechanistic model, e.g. Zhang et al. (2006) is substantial.


NOTE!
The notations used in this chapter are independent from those used in the other chapters of this documentation.


Single Phase Flow Model#


The Darcy friction factor \lambda is applied for the calculation of pressure gradient in single phase flow. For laminar flow, the exact solution of Hagen-Poiseuille applies:

=64/R鈥夆仯e\lambda = 64/\text{R\!e}

where ReRe is the Reynolds number.

The well-established Colebrook-White formula for the friction factor in turbulent pipe flow has the disadvantage of being an implicit expression for the friction factor, and is therefore rarely used directly. The most widespread explicit approximation to Colebrook鈥檚 formula is the one by Haaland [1]:

1=1.8[6.9R鈥夆仯e+(D)1.11]\frac{1}{\sqrt{\lambda}} = -1.8 \left[ \frac{6.9}{\text{R\!e}} + \left(\frac{\epsilon}{D} \right)^{1.11} \right]


For the transition zone, one normally simply uses the maximum of the laminar and turbulent formulas.

A useful correlation that can be applied both for laminar, transitional and turbulent flow is the one by Churchill [2]:

=8[(8R鈥夆仯e)12+1(A+B)1.5]1/12\lambda = 8 \left[ \left(\frac{8}{\text{R\!e}} \right)^{12} + \frac{1}{(A+B)^{1.5}} \right]^{1/12}


where A=[2.457ln((7R鈥夆仯e)0.9+0.27D)]16A = \left[ -2.457\ln \left( \left(\frac{7}{\text{R\!e}} \right)^{0.9} + 0.27\frac{\epsilon}{D} \right) \right]^{16}

where \epsilon is the pipe roughness and DD is the pipe diameter.


Beggs and Brill's Correlation for Multiphase Flow#


Beggs and Brill, in spite of all its limitations, is the only empirical correlation that can be meaningfully applied to all inclination angles, and it involves a flow pattern map [3].

Let CLC_L be the no-slip holdup or liquid loading, CL=USL/UMC_L = U_{SL}/U_M where UM=USG+USLU_M = U_{SG} + U_{SL} is the mixture velocity, and let FrF_r be the mixture Froude number, Fr=UM2/gDF_r = UM^2/gD


Flow Pattern Map - Beggs & Brill#


The flow pattern map is defined by means of the following transition lines:

L1=316CL0.302L2=0.0009252CL2.4684鈥呪鈥呪鈥呪鈥呪L_1^* = 316C_{\text{L}}^{0.302} \hspace{5mm} L_2^* = 0.0009252C_{\text{L}}^{-2.4684}\;\;\;\; L3=0.1CL1.4516L4=0.5CL6.738L_3^* = 0.1C_{\text{L}}^{-1.4516} \hspace{3mm} L_4^* = 0.5C_{\text{L}}^{-6.738}

The flow regimes are then defined as follows:

  • Segregated flow if CL<0.01andF鈥夆仯r<L1orCL0.01andF鈥夆仯r<L2C_{\text{L}}< 0.01 \hspace{1mm}\text{and}\hspace{1mm} \text{F\!r}<L_1^*\hspace{3mm}\text{or} \hspace{3mm} C_{\text{L}}\geq 0.01\hspace{1mm}\text{and}\hspace{1mm}\text{F\!r}<L_2^*
  • Intermittent flow if 0.01CL<0.4andL3<F鈥夆仯rL10.01 \leq C_{\text{L}} < 0.4 \hspace{1mm}\text{and}\hspace{1mm} L_3^*<\text{F\!r}\leq L_1^* orCL0.4andL3<F鈥夆仯rL4\text{or} \hspace{3mm} C_{\text{L}}\geq 0.4\hspace{1mm}\text{and}\hspace{1mm}L_3^*<\text{F\!r}\leq L_4^*
  • Distributed flow if CL<0.4andF鈥夆仯rL1orCL0.4andF鈥夆仯r>L4C_{\text{L}}< 0.4 \hspace{1mm}\text{and}\hspace{1mm} \text{F\!r}\geq L_1^*\hspace{3mm}\text{or} \hspace{3mm} C_{\text{L}}\geq 0.4\hspace{1mm}\text{and}\hspace{1mm}\text{F\!r}>L_4^*
  • Transitional flow if CL0.1andL2<F鈥夆仯r<L3C_{\text{L}} \geq 0.1 \hspace{1mm}\text{and}\hspace{1mm} L_2^* < \text{F\!r} < L_3^*\hspace{3mm}

Liquid Holdup and Hydrostatic Head - Beggs & Brill#


Beggs & Brill gives one liquid holdup correlation for horizontal flow and another correlation for the inclination dependence.

The holdup for horizontal flow is given by the following correlations:

  • Segregated flow: HL(0)=0.98CL0.4846F鈥夆仯r0.0868H_{\text{L}}(0) = 0.98\frac{C_{\text{L}}^{0.4846}}{\text{F\!r}^{0.0868} }
  • Intermittent flow: HL(0)=0.845CL0.5351F鈥夆仯r0.0173H_{\text{L}}(0) = 0.845\frac{C_{\text{L}}^{0.5351}}{\text{F\!r}^{0.0173} }
  • Distributed flow: HL(0)=1.065CL0.5824F鈥夆仯r0.0609H_{\text{L}}(0) = 1.065\frac{C_{\text{L}}^{0.5824}}{\text{F\!r}^{0.0609} }
  • Transitional flow: HL(0)=AHSegregated(0)+(1A)HIntermittent(0)H_{\text{L}}(0) = A H_{\text{Segregated}}(0) + (1-A)H_{\text{Intermittent}}(0)
    where the weight factor AA is given by A=L3F鈥夆仯rL3L2A = \frac{L_3^* - \text{F\!r}}{L_3^* - L_2^*}

The holdup for inclined flow is given by the following correction factor:

HL()=B()HL(0)H_{\text{L}}(\theta) = B(\theta)H_{\text{L}}(0)


where B()=1+[sin(1.8)13sin3(1.8)]B(\theta) = 1 + \beta[\sin(1.8\theta) - \frac{1}{3}\sin^3(1.8\theta)]

The inclination factor \beta is a function of the flow regime.

Defining a liquid velocity number NvlN_{vl} as: Nvl=USL(Lg)1/4N_{vl} = U_{\text{SL}} \left( \frac{\rho_{\text{L}}}{g\sigma} \right)^{1/4}

We then get the following inclination factors \beta:

  • Uphill, segregated flow: =(1CL)ln[0.011Nvl3.539CL3.768F鈥夆仯r1.614]\beta = (1 - C_{\text{L}}) \ln \left[ \frac{0.011 N_{vl}^{3.539}}{C_{\text{L}}^{3.768} \text{F\!r}^{1.614}} \right]
  • Uphill, intermittent flow: =(1CL)ln[2.96CL0.305F鈥夆仯r0.0978Nvl0.4473]\beta = (1 - C_{\text{L}}) \ln \left[ \frac{2.96 C_{\text{L}}^{0.305} \text{F\!r}^{0.0978}}{N_{vl}^{0.4473}} \right]
  • Uphill, distributed flow: =0\beta = 0
  • Downhill, all flow regimes: =(1CL)ln[4.70Nvl0.1244CL0.3692F鈥夆仯r0.5056]\beta = (1 - C_{\text{L}}) \ln \left[ \frac{4.70 N_{vl}^{0.1244}}{C_{\text{L}}^{0.3692} \text{F\!r}^{0.5056}} \right]

Frictional Pressure Drop - Beggs & Brill#


The frictional pressure drop is given by:

(p/z)fric=(/2)fUM2/D(\partial p / \partial z)_{\text{fric}} = -(\lambda/2) \rho_f U_{\text{M}}^2/D


where \lambda is the Darcy friction factor.

A two-phase multiplier is used to calculate the effect of slip on frictional pressure losses. The multiplier is defined in terms of a parameter SS , which is defined as follows:

f2fnoslip=eS\frac{f_{2\phi}}{f_{\text{noslip}}} = e^{S}


Where fnoslipf_{noslip} is absed on a mixture Reynolds number R鈥夆仯eM=UMDNS/NS\text{R\!e}_{\text{M}}=U_{\text{M}}D \rho_{\text{NS}}/\mu_{\text{NS}}, and SS is a function of a nondimensional parameter yy.

y=lnCLHL()2y = \ln \frac{C_{\text{L}}}{H_{\text{L}}(\theta)^2}


The index NS denotes average values at no-slip conditions, i.e.

NS=CLL+(1CL)G\rho_{\text{NS}} = C_{\text{L}} \rho_{\text{L}} + (1-C_{\text{L}})\rho_{\text{G}}

NS=CLL+(1CL)G\mu_{\text{NS}} = C_{\text{L}} \mu_{\text{L}} + (1-C_{\text{L}})\mu_{\text{G}}

If 1<y<1.21 < y <1.2, then S=ln(2.2y1.2)S = ln(2.2y-1.2)

Otherwise, S=y0.01853y40.872y2+3.18y0.0523S = \frac{y}{0.01853y^4 - 0.872y^2 +3.18y - 0.0523 }


Zhang's Multiphase Flow Model#


The multiphase flow model is implemented based on the Zhang (2003) two-phase flow model [4]. The complete descriptions of the different flow regimes are added in this part.


Computational Sequence#


The overall computational sequence is presented in this section for Zhang's model. The inputs of the flow model include well diameters, inclination angles, and fluid properties for each phase containing density, viscosity, surface tension, and superficial velocity.

  • The first step is to decide if it is single phase flow or multiphase flow with all the inputs mentioned above.
  • The flow regime is predicted using simple correlations or models for multiphase flow.
  • The multiphase flow behavior is calculated using flow regime specific momentum equations and required closure relationships.
  • At last, the pressure gradient, liquid holdup, and slug characteristics are calculated in gas-liquid flow at all inclination angles.

Computing the transition from slug flow to dispersed bubble flow requires solving the full equations for slug flow. This involves solving one nonlinear algebraic equation for the liquid holdup in the film zone of the Taylor bubble in Zhang et al (2003a) [4].

L(vTvF)(vSvF)C(vTvC)(vSvC)lFFSFHLFA+CSC(1HLF)A+ISI(1HLFA1(1HLF)A)(LC)gsin=0\frac{\rho_{\text{L}}(v_{\text{T}} - v_{\text{F}})(v_{\text{S}} - v_{\text{F}}) - \rho_{\text{C}}(v_{\text{T}} - v_{\text{C}})(v_{\text{S}} - v_{\text{C}}) }{l_{\text{F}}} - \frac{\tau_{\text{F}} S_{\text{F}}}{H_{\text{LF}}A} + \frac{\tau_{\text{C}} S_{\text{C}}}{(1-H_{\text{LF}})A} +\tau_{\text{I}}S_{\text{I}}\left(\frac{1}{H_{\text{LF}}A} - \frac{1}{(1-H_{\text{LF}})A} \right) - (\rho_{\text{L}} - \rho_{\text{C}}) g \sin\theta = 0

And the transition from slug flow to bubbly flow is simple and explicit.

The transition from slug flow to stratified/annular flow again requires solving the above equation, but now with the additional condition that the slug fraction FS=LS/LU=0F_S = L_S/L_U = 0, giving the equations to solve the liquid holdup of the film. If the flow regime then is stratified, annular or slug flow, the equation then needs to be solved for the appropriate conditions. The momentum equation for stratified flow needs to remove the momentum exchange term in the above equation,

FSFHLFA+CSC(1HLF)A+ISI(1HLFA1(1HLF)A)(LC)gsin=0- \frac{\tau_{\text{F}} S_{\text{F}}}{H_{\text{LF}}A} + \frac{\tau_{\text{C}} S_{\text{C}}}{(1-H_{\text{LF}})A} + \tau_{\text{I}}S_{\text{I}}\left(\frac{1}{H_{\text{LF}}A} - \frac{1}{(1-H_{\text{LF}})A} \right)- (\rho_{\text{L}} - \rho_{\text{C}}) g \sin\theta = 0

Based on the above equation, the gas core shear stress is ignored for annular flow,

FSFHLFA+ISI(1HLFA1(1HLF)A)(LC)gsin=0- \frac{\tau_{\text{F}} S_{\text{F}}}{H_{\text{LF}}A} + \tau_{\text{I}}S_{\text{I}}\left(\frac{1}{H_{\text{LF}}A} - \frac{1}{(1-H_{\text{LF}})A} \right)- (\rho_{\text{L}} - \rho_{\text{C}}) g \sin\theta = 0

The selection between stratified and annular flow is based on the model by Grolman (1994).

The above momentum balance equation for each flow regime, which is actually a function of liquid holdup in the film HLFH_{LF} can't be solved directly, but the initial interval is known as [0,1]. A solution in the interval will be obtained using an interpolation method.

The algorithm is improved in the present work for computational robustness and accuracy based on Zhang's work [4]. A unified algorithm for different flow regimes is applied in the computational procedure, where the balance of momentum equation is required for each flow regime. In each iteration, the following caculation steps from 2-4 repeat until the convergence criterion is satisfied. If we rewrite the momentum equation so the left-hand side is equal to \xi, the convergence criterion arrives when \xi is a small value.


Slug Flow Model#


Zhang's model used the Taylor bubble region, and therefore needs to include the momentum transfer between the liquid in the slug body and in the film under the Taylor bubble [4].

The calculation steps are as follows:

  1. Compute the translational velocity (Taylor bubble velocity) vTv_T using the mixture velocity vSv_S, the slug length lSl_S, and the entrained fraction FEF_E. The slug liquid holdup HLSH_{LS} can be estimated with the existing correlation. Guess a value of HLFH_{LF}.
  2. The following equations need to be solved for vFv_F, vCv_C, and lFl_F. vF=vTHLSHLF(vTvSL)v_{\text{F}} = v_{\text{T}} - \frac{H_{\text{LS}}}{H_{\text{LF}}} ( v_{\text{T}} - v_{\text{SL}} ), vC=vT1HLS1HLF(vTvS)v_{\text{C}} = v_{\text{T}} - \frac{1- H_{\text{LS}}}{1 - H_{\text{LF}}} (v_{\text{T}} - v_{\text{S}}), and lF=lSHLSvSvSLvSLHLFvFl_{\text{F}} = l_{\text{S}} \frac{H_{\text{LS}}v_{\text{S}} - v_{\text{SL}} }{v_{\text{SL}} - H_{\text{LF}}v_{\text{F}} }. If the film length lF0l_F \leq 0, there is no solution for slug or bubble flow and we will get stratified/annular flow.
  3. Compute geometry and hydraulic diameter \Theta, SFS_F, SCS_C, SIS_I, and friction factors fFf_F, fCf_C, fIf_I. These equations are explicit and straightforward.
  4. Compute shear stresses f\tau_f, C\tau_C and I\tau_I.
  5. Check residual \xi and iterate until convergence is achieved.

The Stratified/Annular Flow Model#


A double-circle model to estimate the wall wetting of the liquid layer is applied in Zhang's work [4]. Zhang also uses a different computational scheme that we would probably like to avoid, solving for the liquid film velocity instead of the liquid holdup. There is an advantage in terms of robustness in solving for the liquid holdup, which is bounded, instead of for the liquid film velocity, which is unbounded.

  1. Compute the translational velocity (Taylor bubble velocity) vTv_T using the mixture velocity vSv_S, the slug length lSl_S, and the entrained fraction FEF_E. Guess a value of HLFH_{LF}.
  2. Calculate vFv_F, vCv_C, and lFl_F. vF=(1FE)vSLHLFv_{\text{F}} = \frac{(1-F_{\text{E}}) v_{\text{SL}}}{H_{\text{LF}}}, vC=vSL+vSGHLFvF1HLFv_{\text{C}} = \frac{v_{\text{SL}} + v_{\text{SG}} - H_{\text{LF}} v_{\text{F}} }{ 1 - H_{\text{LF}} } and HLC=FEvSL/vCH_{\text{LC}} = F_{\text{E}} v_{\text{SL}}/v_{\text{C}}
  3. Compute geometry and hydraulic diameters \Theta, SFS_F, SCS_C, SIS_I, and friction factors fFf_F, fCf_C, fIf_I. These equations are explicit and straightforward.
  4. Compute shear stresses f\tau_f, C\tau_C and I\tau_I for stratified or annular flow.
  5. Check residual \xi and iterate until convergence is achieved.

Closure Relations#


There are a substantial number of closure relations, typically semi-empirical correlations, needed to close the model and run the solver. Most of them are simple and explicit and listed in Zhang et al (2003a).

  • Wall fraction factors \Theta - Standard hydraulic approximation is applied with Blasius' correlation for turbulent flow Re>3000Re > 3000 and Hagen-Poiseuille for laminar flow Re<2000Re < 2000. For 2000<Re<30002000 < Re < 3000, an interpolation is applied. The wetted wall fraction is taken from Grolmam (1994).
  • Interfacial friction factor fIf_I - Andritsos and Hanratty for stratified flow, Asali for annular flow. Interfacial perimeter from Grolman. All given in the paper.
  • Droplet entrainment FEF_E in the gas core proposed by Oliemans et al (1986).
  • Liquid holdup in slug HLSH_{LS}- this model is a bit more complicated. Since this model requires input that has to be computed by the slug flow model, e.g. the film velocity and the gas core velocity, a simpler correlation has to be used in the first iteration on the slug flow model. For this purpose, the correlation by Gregory et al (1978) is used.
  • The translational velocity vTv_T of the liquid slugs is expressed as a function of mixture velocity proposed by Nicklin (1962). The estimation for the slug length lSl_S is listed in the paper [4].

Flow Regime Transitions#


The calculation of the transition between two flow patterns is summarized in the following.


The transition from Slug to Dispersed Bubble Flow


The transition from slug to dispersed bubble flow is calculated using a principle proposed by Barnea et al (1985). The basic assumption is that the liquid slugs accommodate the same gas volume fraction as dispersed bubble flow with the same mixture velocity on the slug/dispersed bubble flow transition boundary. The governing equation for the transition is claimed to be (13) in Zhang et al (2003b) [5].

32[fs2sVm2+d4LHLF(VTVF)(VmVT)ls]HLS(2.5sin2)6db(1HLS)\frac{3}{2} \left[\frac{f_s}{2}\rho_s V^2_m + \frac{d}{4}\frac{\rho_{\text{L}}H_{\text{LF}}(V_{\text{T}} - V_{\text{F}})(V_m - V_{\text{T}})}{l_s}\right] H_{\text{LS}} \leq \left(\frac{2.5 - |sin \theta|}{2} \right) \frac{6\sigma}{d_b}(1-H_{\text{LS}})

where the bubble size dsd_s is given by db=20.4(LG)gd_b = 2 \sqrt{\frac{0.4 \sigma}{(\rho_{\text{L}} - \rho_{\text{G}})g}}

Once the slug flow model has been solved, the equation above can be solved for the mixture velocity Vm=VSG+VSLV_m = V_{SG}+V_{SL}. Here we have reformulated the equation as inequality, and dispersed bubble flow occurs when VmV_m exceeds the critical value given by this inequality. This criterion applies only to dispersed bubble flow and not to bubbly flow.


The Transition from Slug to Bubbly Flow


The transition from slug to bubbly flow (low velocity bubbly flow) is assumed to occur at a mean void fraction of 0.25. In other words, we will not have slug flow when the mean void fraction is equal to or less than 0.25. This criterion applies only to bubbly flow and not to dispersed bubble flow. We will not use a criterion based on comparing the translational velocity and the small bubble velocity.


The Transition from Slug to Stratified/Annular Flow


This transition criterion is described in detail by the following equations in Zhang (2003a) [4], or alternatively by computing a film length that is negative or zero from the slug flow model. We will start out assuming the latter criterion.

HLF=(HLS(vTvS)+vSL)(vSG)+vSLFEvTvSLFEvTvSGH_{\text{LF}} = \frac{ (H_{\text{LS}}(v_{\text{T}} - v_{\text{S}}) + v_{\text{SL}})(v_{\text{SG}}) + v_{\text{SL}}F_{\text{E}} - v_{\text{T}}v_{\text{SL}}F_{\text{E}} }{v_{\text{T}}v_{\text{SG}} }

where vFv_F, HLCH_{LC}, and vCv_C is expressed as,

vF=vSL(1FE)HLFv_F = \frac{v_{\text{SL}}(1-F_{\text{E}})}{H_{\text{LF}}}

HLC=vSLFE(1HLF)vSvSL(1FE)H_{LC} = \frac{v_{\text{SL}}F_{\text{E}}(1-H_{\text{LF}})}{v_{\text{S}} - v_{\text{SL}}(1-F_{\text{E}})}

vC=vSvSL(1FE)1HLFv_C = \frac{v_{\text{S}} - v_{\text{SL}}(1-F_{\text{E}})}{1-H_{\text{LF}}}

When we compute slug/dispersed bubble/bubbly flow in downward flow, we assume that the Taylor bubbles in slug flow do not move against the flow or that the small bubbles in dispersed bubble/bubbly flow do not have a net velocity against the flow direction. In these two cases, we assume that we have stratified/annular flow. Slug flow is the prevailing regime when slug flow has a solution with a positive film length and we haven鈥檛 predicted a transition to dispersed bubble flow or bubbly flow, even if both slug and stratified/annular flow have physically realizable solutions.


The Transition from Stratified to Annular Flow


This transition boundary is governed by the correlation for wetted wall fraction, taken from the work of Grolman (1994), given in Zhang et al (2003a) [4].


Pressure Drop#


Stratified Flow


The general expression for the pressure drop for stratified flow is obtained by adding the momentum equations for gas and liquid:

pz=GSGAGLSLAL[HLL+(1HL)G]gsin\frac{\partial p}{\partial z} = - \frac{\tau_{\text{G}} S_{\text{G}}}{A_{\text{G}}} - \frac{\tau_{\text{L}} S_{\text{L}}}{A_{\text{L}}} - [ H_{\text{L}} \rho_{\text{L}} + (1-H_{\text{L}})\rho_{\text{G}} ]g\sin\theta

Here GG and LL are the gas and liquid shear stresses as outlined in the paper, SGS_G and SLS_L are the wetted perimeters and AGA_G and ALA_L are the phase areas, given by the phase fractions and pipe cross sectional area.

The wetted perimeters are given by SL=D,SG=D()S_{\text{L}} = D \sigma, \hspace{5mm} S_{\text{G}} = D (\pi - \sigma)

where \sigma is the liquid wetted angle, given by Biberg's formula.


Annular Flow

The expressions for annular flow are similar as for stratified flow, except for two things:

  • The liquid is the only phase in contact with the wall, so SG=0S_G = 0 in the pressure drop expression.
  • The liquid wetted perimeter is then SL=DS_L = \pi D.

Dispersed Bubble Flow or Bubbly Flow


The expressions for dispersed bubble flow and bubbly flow are similar as for annular flow, except that the gas and liquid are treated as a homogeneous mixture filling the entire pipe cross section. To calculate the liquid wall shear stress, we use the volume averaged density averaged by the mean liquid holdup, but we use the liquid viscosity in the Reynolds number.

The wetted perimeter is then equal to the pipe perimeter D\pi D and the velocity equal to the mixture velocity vSG+vSLv_{SG} + v_{SL}.


Slug Flow


In slug flow, the pressure drop in the film/bubble region is computed as for stratified or annular flow as given by the stratified-annular transition criterion. The pressure drop in the slug region is computed as for dispersed bubble flow. The average pressure drop is computed as the mean of the pressure drops in the film/bubble region and the slug region, weighted with the slug fraction FSF_S.The wetted perimeters and friction factors for the film/bubble region are computed using the velocities vFv_F and vCv_C, and the phase fractions HFH_F and 1HF1-H_F for the liquid and gas core respectively.

The wetted perimeter for the slug region is thus equal to the pipe perimeter D\pi D and the velocity equal to vsv_s, which equals the mixture velocity vSG+vSLv_{SG}+v_{SL}. To calculate the liquid wall shear stress, we use the volume averaged density averaged using the liquid holdup in slug $H_{LS}, but again we use the liquid viscosity in the Reynolds number.


References#


[1] Haaland, S.E., Simple and Explicit Formulas for the Friction Factor in Turbulent Flow. Journal of Fluids Engineering, 1983, Vol. 105, pp. 89-90.

[2] Churchill, S.W., Friction factor equation spans all fluid-flow regimes. Chemical engineering, 1977, Vol. 84, pp. 91-92.

[3] Beggs, H. D., Brill, J.P., A Study of Two-Phase Flow in Inclined Pipes. Journal of Petroleum Technology, May 1973, pp. 607-617.

[4] Zhang, H.Q., Wang, Q., Sarica, C. and Brill, J. P., Unified Model for Gas-Liquid Pipe Flow via Slug Dynamics -- Part 1: Model Development. Journal of Energy Resources Technology, December 2003, Vol. 125, pp. 266-273.

[5] Zhang, H.Q., Wang, Q., Sarica, C. and Brill, J. P., Unified Model for Gas-Liquid Pipe Flow via Slug Dynamics -- Part 2: Model Validation. Journal of Energy Resources Technology, December 2003, Vol. 125, pp. 274-283.