Next Article in Journal
Refraction of Oblique Shock Wave on a Tangential Discontinuity
Previous Article in Journal
Online Coupled Generalized Multiscale Finite Element Method for the Poroelasticity Problem in Fractured and Heterogeneous Media
Previous Article in Special Issue
An Overview of the Lagrangian Dispersion Modeling of Heavy Particles in Homogeneous Isotropic Turbulence and Considerations on Related LES Simulations
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

An Explicit Analytical Solution for Transient Two-Phase Flow in Inclined Fluid Transmission Lines

Cullen College of Engineering, Mechanical Engineering, University of Houston, 4722 Calhoun Rd., Houston, TX 77204, USA
*
Author to whom correspondence should be addressed.
Fluids 2021, 6(9), 300; https://doi.org/10.3390/fluids6090300
Submission received: 7 June 2021 / Revised: 25 July 2021 / Accepted: 13 August 2021 / Published: 24 August 2021
(This article belongs to the Special Issue Numerical Methods and Physical Aspects of Multiphase Flow)

Abstract

:
Due to the complex nonlinearity characteristics, analytical modeling of compressible flow in inclined transmission lines remains a challenge. This paper proposes an analytical model for one-dimensional flow of a two-phase gas-liquid fluid in inclined transmission lines. The proposed model is comprised of a steady-state two-phase flow mechanistic model in-series with a dynamic single-phase flow model. The two-phase mechanistic model captures the steady-state pressure drop and liquid holdup properties of the gas-liquid fluid. The developed dynamic single-phase flow model is an analytical model comprised of rational polynomial transfer functions that are explicitly functions of fluid properties, line geometry, and inclination angle. The accuracy of the fluid resonant frequencies predicted by the transient flow model is precise and not a function of transmission line spatial discretization. Therefore, model complexity is solely a function of the number of desired modes. The dynamic single-phase model is applicable for under-damped and over-damped systems, laminar, and turbulent flow conditions. The accuracy of the overall two-phase flow model is investigated using the commercial multiphase flow dynamic code OLGA. The mean absolute error between the two models in step response overshoot and settling time is less than 8% and 2 s, respectively.

Graphical Abstract

1. Introduction

The oil and gas industry has been looking for a long time into the development of reliable technologies and practices to be able to properly assess and understand the dynamics of multiphase flows within transmission lines. In upstream-downstream production systems, an accurate evaluation of the flow transients allows the determination of the dynamic forces acting on the pipeline wall and the expected flow rate facilitating the assessment of the well architecture and type of material as well as the choice of the safety equipment placement and redundancy. Onshore and offshore drilling is another field in which flow transients have an impact on the progress and safety of drilling operation. For instance, influxes entering the wellbore result in a change of the circulating fluid flow rate and pressure. Accurately knowing the flow dynamics allows the adjustment of the injected mudflow rate and weight to avoid undesirable events, such as gas kicks.
Looking into the standard tools used within the oil and gas industry, a variety of commercial packages have been introduced, such as OLGA and LEADAFLOW. Although these codes have been developed using physics based models, they have been simplified and calibrated using empirical correlations and in-field data. Such hybrid approached have been seen, despite the accuracy of these codes, altering the ability to generalize these methods in order to cover the wide range of applications as well as the variety of investigations and analyses users in such an industry are looking to perform. Hence, having a simplified, analytical yet accurate model is crucial so that transmission line flow models can be integrated with other oil and gas equipment models to analyze the behavior of the whole complex system, achieve tasks that require real-time calculations as well as having a fast and reliable tool that can be deployed during the design and optimization phase of oil and gas development projects.
Several methodologies have been followed in order to develop steady-state multiphase flow models. The main driver that differentiates each of the proposed methods is mainly the way the different fluid phases are treated. Specifically, three different approaches have been introduced within the literature: a homogenous approach was proposed in [1,2], in which the two phases are considered traveling at the same velocity; hence, the total flow can be assumed to be a single phase flow. In [3], a separated flow model is presented, in which the authors assumed that the two phases’ velocities are different, which affects the fluid conservation equations. A multi-fluid model was presented in [4] based on the assumption that separate conservation equations are introduced for each of the two phases to highlight the interaction between both phases. Finally, drift flux models were proposed in [5,6], in which a distribution parameter and an average of the local fluid velocity difference between the two phases are used to describe the overall flow.
Another modeling approach was introduced via mechanistic models as an alternative to the methods described above. Mathematically speaking, the main effort when developing most of the mechanistic models is to identify the flow regime in question. Numerous efforts have been placed toward this direction. Taitel and Dukler [7] started by laying out a procedure to determine the transitions between the different flow regimes. This was the beginning of multiple investigations in which researchers were looking into developing analytical approaches in order to accurately predict the different flow patterns as well as the regimes’ transition boundaries. However, the proposed techniques have been accompanied with shortcomings, such as only considering flow pattern determination [8] or the applicability to a limited range of pipe inclinations [9]. Petalas and Aziz [10] proposed a model that was ameliorated for the case of specific flow regimes in [11], to surpass these limitations. Although mechanistic models are essentially based on fundamental laws and they offer more accurate predictions when considering the variations of both fluid and geometry properties, they still depend on closure relationships based on observations.
In addition to steady-state flow conditions, numerous studies have been conducted to model transient two-phase flow systems based on either a two-fluid model [12] or a drift-flux model [13]. In both cases, the models are developed mainly on the basis of introducing a full set of conservation equations (mass, momentum, and energy). To determine approximated solutions of the Navier-Stokes equations, a numerical modeling approach was employed. In this case, when we are handling a complex system, thousands of simulations are required, which can be very time consuming. To overcome the numerical modeling limitations, reduced-order modeling techniques have been suggested to simulate two-phase dynamics in a fluid transmission line. In [14], a transient single-phase flow model is integrated with a steady-state two-phase model. The single-phase flow model is derived from the well-known dissipative model, which is regarded as the most accurate analytical model for transient laminar flow in horizontal transmission lines. Since the dissipative model is described in Laplace domain and its transfer functions contain hyperbolic and Bessel functions, approximated rational polynomial transfer functions are developed to enable simulations in the time domain. Similarly, to extend the results of the dissipative model to lines with different inclinations, a reduced-order model is developed in [15]. The resulting single-phase model is used in [16] to simulate transient two-phase flow in inclined transmission lines. The results are compared with those simulated using OLGA software showing good agreement. A major drawback of these approximated rational transfer functions is that the physical properties are lost in the model coefficients since they are calculated numerically.
The present development focuses on the development of an analytical model whose coefficients are functions of fluid properties and transmission line properties and inclination for single-phase flow conditions. A comprehensive mechanistic model is then adopted to estimate steady-state liquid holdup and pressure gradient, which will be used to derive equivalent single-phase fluid parameters, namely, altering the fluid parameters of density, viscosity, and bulk modulus as a function of gas volume fraction (GVF). These steady-state parameters are finally integrated within the analytical dynamic single-phase model to predict transient flow behavior in inclined transmission lines.
The remainder of the paper is structured as follows. In Section 2, a reduced-order solution of Navier-Stokes equations written in Laplace domain is derived for inclined transmission lines by accounting only for steady state friction losses. In Section 3, the solution is approximated using rational polynomial transfer functions and extended to include unsteady-state friction effects. In Section 4, the proposed model is applied to estimate transient air-water flow dynamics in inclined transmission lines. Concluding remarks are presented in Section 5.

2. Reduced Order Modeling of Inclined Compressible Single Phase Flow

The solution for transient flow in transmission line is obtained from the continuity, momentum, energy, and state equations
ρ ¯ u t = P x + μ ( 2 u r 2 + 1 r u r ) ρ g   sin ( θ ) ,
  ρ t + ρ ¯ ( u x + v r + v r ) = 0 ,
d T d t + T ¯   ( γ 1 ) ρ t = α ( 2 T r 2 + 1 r T r ) ,
  d ρ ρ ¯ = d P β e ,
              d ρ ρ ¯ + d T T ¯ = d P P ¯   ,
where x and r are independent cylindrical space variables, and t is the independent time variable. The dependent fluid variables include pressure P , its time average P ¯ , velocities in x-direction and r-direction ( u and v ), temperature T , its time average T ¯ , density ρ , its time average ρ ¯ , fluid absolute viscosity μ , fluid effective bulk modulus β e , diffusion coefficient α , specific heat ratio γ , and angle of inclination θ . For a horizontal flow, the set of Equations (1)–(5) describes the so-called dissipative model and is regarded as the accurate model for laminar flow case since it includes viscosity and heat transfer effects [15].
Similar to the approach adopted in [17], the second-order effects, namely, the unsteady friction losses and the heat transfer effects, are temporarily neglected and only laminar flow conditions are considered in order to retain a one-dimensional model. The only difference is that the present development accounts for transmission line inclination. The resulting reduced-order system of equations that describes transient flow in inclined transmission line can be written in Laplace domain as
I x s Q + R x Q + G x P = d P d x ,
C x s P = d Q d x ,
where s is the Laplace transform variable, and Q is the flow rate and pressure along the transmission line. The distributed lumped parameters in Equations (6) and (7) are R x , I x , C x , and G x , defined as the resistance, inertance, capacitance, and gravity term per unit length of transmission line and are presented in Table 1, where c ¯ is the time averaged speed of sound in the fluid, g is the gravity acceleration, and A is the transmission line cross-sectional area [17]
Using Laplace Transform, the analytical solution for Equations (6) and (7) can be represented as a two-port model (four-terminal network) and written in a matrix form as
[ P i n ( s ) Q o u t ( s ) ] = [ T F 11 ( s ) T F 12 ( s ) T F 21 ( s ) T F 22 ( s ) ] [ P o u t ( s ) Q i n ( s ) ] ,
where
T F 11 ( s ) = e G x L 2 c o s h 1 ( Γ L 2 ) 1 G x Γ   t a n h ( Γ L 2 ) ,
T F 12 ( s ) = 2 ( I x s   +   R x ) Γ t a n h ( Γ L 2 ) 1 G x Γ   t a n h ( Γ L 2 ) ,
T F 21 ( s ) = 2 C x s Γ   t a n h ( Γ L 2 ) 1 G x Γ   t a n h ( Γ L 2 ) ,
T F 22 ( s ) = e G x L 2 c o s h 1 ( Γ L 2 ) 1 G x Γ   t a n h ( Γ L 2 ) ,
and
Γ = 4 ( I x C x s 2 + R x C x s ) + G x 2
Illustrated in Figure 1 is the sensitivity of the matrix transfer functions T F 11 ( s ) , T F 12 ( s ) , T F 21 ( s ) , and T F 22 ( s ) to the transmission line inclination angle for the laminar flow case, given the parameters shown in Table 2.
At low frequencies, the inclination angle θ has no effect on the output Q o u t ( s ) since the line inlet and outlet flow rates are equal at steady state conditions and independent of the outlet pressure. However, the angle θ has an effect on the output P i n ( s ) as expected from the momentum equation. As the inclination angle increases, the difference between the line inlet and outlet pressures increases. It is also concluded that the system described in Equation (8) is sensitive to the inclination angle at middle frequencies (first harmonics). As the transmission line inclination increases, a slight decrease in natural frequencies is observed along with an increase in damping ratios.
The four-terminal network representation obtained in Equation (8) is causal and known as the Hybrid Representation. It is also referred to as the PQ-Model and represents mixed boundary conditions system, where a pressure at one end of the line and a flow rate at the other end are used as inputs [18]. Based on all possible boundary conditions, there are three other causalities to fluid transmission line, namely the Hybrid Representation (PQ-Model), the Impedance Representation (also known as the Q-Model), and the Admittance Representation (also known as the P-Model). Using the transfer function matrix presented in Equation (8), the other three causal representations to the inclined fluid transmission line problem are defined as
Hybrid Representation:
[ P o u t ( s ) Q i n ( s ) ] = [ T F 22 ( s ) T F 12 ( s ) T F 21 ( s ) T F 11 ( s ) ] [ P i n ( s ) Q o u t ( s ) ]
Impedance Representation:
[ P i n ( s ) P o u t ( s ) ] = [ T F 12 ( s ) 1     T F 11 2 ( s ) T F 11 ( s ) T F 12 ( s ) 1     T F 11 2 ( s ) T F 11 ( s ) T F 12 ( s ) 1     T F 11 2 ( s ) T F 12 ( s ) 1     T F 11 ( s ) 2 ] [ Q i n ( s ) Q o u t ( s ) ]
Admittance Representation:
[ Q i n ( s ) Q o u t ( s ) ] = [ T F 21 ( s ) 1     T F 22 2 ( s ) T F 22 ( s ) T F 21 ( s ) 1     T F 22 2 ( s ) T F 22 ( s ) T F 21 ( s ) 1     T F 22 2 ( s ) T F 21 ( s ) 1     T F 22 2 ( s ) ] [ P i n ( s ) P o u t ( s ) ]

3. Analytical Fluid Transmission Line Model

Developed in this section is the explicit analytical model for transient compressible single-phase flow in inclined fluid transmission lines. First, a modal approximation is developed such that the matrix transfer functions of the solution (8) are approximated by rational polynomial transfer functions. Since these transfer functions are accounted only for linear friction effects, the approximated transfer functions are then corrected to incorporate second order effects [17].

3.1. Modal Approximation

The obtained model in matrix form can be used to perform frequency-domain analyses. However, it is impractical for time domain analyses and a modal approximation procedure is required [17]. The matrix transfer functions in Equation (8) can be expressed as functions of four transcendental transfer functions that can be approximated as
c o s h 1 ( Γ L 2 ) C 1 i = 1 m a c i s + b c i s 2 + 2   ξ i ω n i   s + ω n i 2 ,
2 ( I x s + R x ) Γ t a n h ( Γ L 2 ) C 2 i = 1 m a r t i s + b r t i s 2 + 2   ξ i ω n i   s + ω n i 2 ,
2 C x s Γ   t a n h ( Γ L 2 ) C 3   i = 1 m s ( a c t i s + b c t i ) s 2 + 2   ξ i ω n i   s + ω n i 2 ,
G x Γ   t a n h ( Γ L 2 ) C 4 i = 1 m a g t i s + b g t i s 2 + 2   ξ i ω n i   s + ω n i 2 ,
where m is the desired number of modes to be included in the approximation. The denominator coefficients for each ith mode are the natural frequency ω n i and the damping ratio ξ i . The numerators coefficients a c i , b c i , a r t i , b r t i , a c t i , b c t i , a g t i , and b g t i are determined using the residue theorem [17]. The DC gain factors C 1 , C 2 , C 3 , and C 4 are incorporated to adjust the steady-state conditions due to the finite number of modes used to approximate the transcendental transfer functions. All the transfer functions in Equation (8) have the same poles and can be computed using the following series expansion
c o s h ( Γ L 2 ) = i = 1 ( 1 + Γ 2 L 2 4 π 2 ( i 0.5 ) 2 ) = 0
Two complex conjugate poles exist for each second-order mode, provided
L π R x C x I x 1 1 + G x 2 L 2 π 2 < 1
The corresponding complex conjugate poles of each ith mode are
s ¯ i , 1 = R x C x + j 4 π 2 ( i 0.5 ) 2   I x C x L 2 R x 2 C x 2 + G x 2 I x C x 2 I x C x ,
s ¯ i , 2 = R x C x j 4 π 2 ( i 0.5 ) 2   I x C x L 2 R x 2 C x 2 + G x 2 I x C x 2 I x C x
which yields to analytical expressions for the natural frequencies and damping ratios
ω n i = π L 1 I x C x ( i 0.5 ) 2 + ( G x L 2 π ) 2 ,
ξ i = L 2 π R x C x I x   1 ( i 0.5 ) 2 + ( G x L 2 π ) 2  
The numerator coefficients can be computed using the residue calculation technique [17] and they are expressed as
a c i = Res [ c o s h 1 ( Γ L 2 ) ] s ¯ i , 1 + Res [ c o s h 1 ( Γ L 2 ) ] s ¯ i , 2 ,
b c i = s ¯ i , 2 Res [ c o s h 1 ( Γ L 2 ) ] s ¯ i , 1 s ¯ i , 1 Res [ c o s h 1 ( Γ L 2 ) ] s ¯ i , 2 ,
a r t i = Res [ 2 ( I x s + R x ) H t a n h ( H L 2 ) ] s ¯ i , 1 + Res [ 2 ( I x s + R x ) Γ t a n h ( Γ L 2 ) ] s ¯ i , 2 ,
b r t i = s ¯ i , 2 Res [ 2 ( I x s + R x ) Γ t a n h ( Γ L 2 ) ] s ¯ i , 1 s ¯ i , 1 Res [ 2 ( I x s + R x ) Γ t a n h ( Γ L 2 ) ] s ¯ i , 2 ,
a c t i = s ¯ i , 1 1 Res [ 2 C x s Γ   t a n h ( Γ L 2 ) ] s ¯ i , 1 + s ¯ i , 2 1 Res [ 2 C x s Γ   t a n h ( Γ L 2 ) ] s ¯ i , 2 ,
b c t i = s ¯ i , 1 1 s ¯ i , 2 Res [ 2 C x s Γ   t a n h ( Γ L 2 ) ] s ¯ i , 1 s ¯ i , 2 1 s ¯ i , 1 Res [ 2 C x s Γ   t a n h ( Γ L 2 ) ] s ¯ i , 2 ,
a g t i = Res [ G x Γ   t a n h ( Γ L 2 ) ] s ¯ i , 1 + Res [ G x Γ t a n h ( Γ L 2 ) ] s ¯ i , 2 ,
b g t i = s ¯ i , 2 Res [ G x Γ   t a n h ( Γ L 2 ) ] s ¯ i , 1 s ¯ i , 1 Res [ G x Γ t a n h ( Γ L 2 ) ] s ¯ i , 2 ,
where the corresponding residues are obtained as
Res [ c o s h 1 ( Γ L 2 ) ] s = s ¯ i , k = [ 2 Γ L sinh ( Γ L 2 ) ] s = s ¯ i , k = [ 2 π ( 1 ) i + 1 ( i 0.5 ) ( 2 I x C x s + R x C x ) L 2 ] s = s ¯ i , k ,
Res [ G x Γ   t a n h ( Γ L 2 ) ] s = s ¯ i , k = [ 2 G x Γ   Γ L ] s = s ¯ i , k = [ G x ( 2 I x C x s + R x C x ) L ] s = s ¯ i , k ,
Res [ 2 ( I x s + R x ) Γ t a n h ( Γ L 2 ) ] s = s ¯ i , k = [ 4 ( I x s + R x )   Γ   Γ L ] s = s ¯ i , k = [ 2 ( I x s + R x ) ( 2 I x C x s + R x C x ) L ] s = s ¯ i , k ,
Res [ 2 C x s Γ   tanh ( Γ L 2 ) ] s = s ¯ i , k = [ 4 C x s Γ   Γ L ] s = s ¯ i , k = [ 2 C x   s ( 2 I x C x s + R x C x ) L ] s = s ¯ i , k
Finally, the approximated rational transfer functions in Equations (17)–(20) are adjusted such that their DC gains match accurately those of the analytical transfer functions in Equations (9)–(12). The steady state correction factors are obtained as
C 1 = [ c o s h ( G x L 2 ) i = 1 m ( 1 ) i ( 1 2 i ) π π 2 ( i 0.5 ) 2 + ( G x L 2 ) 2 ] 1 ,
          C 2 = C 3 ,
          C 3 = C 4 ,
C 4 = { [ i = 1 m 2 π 2 ( i 0.5 ) 2 ] 1 i f   θ = 0 t a n h ( G x L 2 ) G X L [ i = 1 m 1 π 2 ( i 0.5 ) 2   +   ( G x L 2 ) 2 ] 1 i f   θ 0
Hence, the resulting approximated matrix transfer functions in Equations (9)–(12) for oscillating pressure and flow rate in a fluid inclined transmission line can be written as
T F 11 ( s ) e G x L 2 C 1 i = 1 m ( 1 ) i ( 1     2 i ) π L 2 I x C x s 2   +   R x I x   s   +   1 L 2 I x C x ( π 2 ( i     0.5 ) 2   +   ( G x L 2 ) 2 ) 1 C 4 i = 1 m G x I x C x L s 2   +   R x I x   s   +   1 L 2 I x C x ( π 2 ( i     0.5 ) 2   +   ( G x L 2 ) 2 ) ,
T F 22 ( s ) e G x L 2 C 1 i = 1 m ( 1 ) i ( 1     2 i ) π L 2 I x C x s 2   +   R x I x   s   +   1 L 2 I x C x ( π 2 ( i     0.5 ) 2   +   ( G x L 2 ) 2 ) 1 C 4 i = 1 m G x I x C x L s 2   +   R x I x   s   +   1 L 2 I x C x ( π 2 ( i   0.5   ) 2   +   ( G x L 2 ) 2 ) ,
T F 12 ( s ) C 2 i = 1 m 2 C x L s   +   2 R x I x C x L s 2   +   R x I x   s   +   1 L 2 I x C x ( π 2 ( i 0.5 ) 2   +   ( G x L 2 ) 2 ) 1 C 4 i = 1 m G x I x C x L s 2   +   R x I x   s   +   1 L 2 I x C x ( π 2 ( i     0.5 ) 2   +   ( G x L 2 ) 2 ) ,
T F 21 ( s ) C 3 i = 1 m 2 I x L s s 2   +   R x I x   s   +   1 L 2 I x C x ( π 2 ( i     0.5 ) 2   +   ( G x L 2 ) 2 ) 1 C 4 i = 1 m G x I x C x L s 2   +   R x I x   s   +   1 L 2 I x C x ( π 2 ( i     0.5 ) 2   +   ( G x L 2 ) 2 )
Illustrated in Figure 2 is the accuracy of the developed approximated transfer functions. A comparison of the frequency responses of the transcendental (analytical) transfer functions (9)–(12) and their rational polynomial approximations (43)–(46) using eight second-order modes is performed using the study case presented in Table 2.
Approximated transfer functions proposed in Equations (43)–(46) are obtained for under-damped transmission lines. The same approach is applied to the other cases in which the inequality condition in (22) is not met. It is found that the approximated transfer functions (43)–(46) remain the same; hence, they are valid approximations for over/well-damped systems. Shown in Figure 3 is a frequency response evaluation of the proposed approximations for two different cases. The first case represents the horizontal transmission line presented in Table 2, which corresponds to an under-damped system ( ξ 1 = 0.25 ). In the second case, a well-damped transmission line ( ξ 1 = 3 ) is considered by decreasing the line cross-sectional area. Agreement is seen between the analytical solution and the proposed approximation for both cases.

3.2. Proposed Analytical Model

The analytical solution in Equation (8) is obtained by solving Equations (1)–(5) assuming Hagan-Poiseuille flow (losses are proportional to mean velocity) and no heat transfer effects [17]. Proposed in this section is an explicit analytical model obtained by incorporating unsteady friction and heat transfer effects in the approximated transfer functions using modal approximations of the dissipative model, which is known as the most accurate analytical model for laminar flow in horizontal transmission lines. Similar to the work presented in [17,19], two unsteady modifications factors are introduced to adjust the numerators and denominators of the approximated transfer functions. Hence, the updated analytical model describing transient flow in an inclined fluid transmission line is proposed as
[ P i n ( s ) Q o u t ( s ) ] = [ T F ˜ 11 ( s )   T F ˜ 12 ( s ) T F ˜ 21 ( s ) T F ˜ 22 ( s ) ] [ P o u t ( s ) Q i n ( s ) ] ,
where
T F 11 ˜ ( s ) = e G x L 2 C 1 i = 1 m κ i ( 1 ) i ( 1     2 i ) π L 2 I x C x s 2   +   τ i R x I x   s   +   κ i 1 L 2 I x C x ( π 2 ( i     0.5 ) 2   +   ( G x L 2 ) 2 ) 1 C 4 i = 1 m κ i G x I x C x L s 2   +   τ i R x I x   s   +   κ i 1 L 2 I x C x ( π 2 ( i     0.5 ) 2   +   ( G x L 2 ) 2 ) ,
T F 22 ˜ ( s ) = e G x L 2 C 1 i = 1 m κ i ( 1 ) i ( 1     2 i ) π L 2 I x C x s 2   +   τ i R x I x   s   +   κ i 1 L 2 I x C x ( π 2 ( i     0.5 ) 2   +   ( G x L 2 ) 2 ) 1 C 4 i = 1 m κ i G x I x C x L s 2   +   τ i R x I x   s   +   κ i 1 L 2 I x C x ( π 2 ( i     0.5 ) 2   +   ( G x L 2 ) 2 ) ,
T F 12 ˜ ( s ) = C 2 i = 1 m 2 C x L s   +   τ i 2 R x I x C x L s 2   +   τ i R x I x   s   +   κ i 1 L 2 I x C x ( π 2 ( i     0.5 ) 2   +   ( G x L 2 ) 2 ) 1 C 4 i = 1 m κ i G x I x C x L s 2   +   τ i R x I x   s   +   κ i 1 L 2 I x C x ( π 2 ( i     0.5 ) 2   +   ( G x L 2 ) 2 ) ,
T F 21 ˜ ( s ) = C 3 i = 1 m κ i 2 I x L S s 2   +   τ i R x I x   s   +   κ i 1 L 2 I x C x ( π 2 ( i 0.5 ) 2   +   ( G x L 2 ) 2 ) 1 C 4 i = 1 m κ i G x I x C x L s 2   +   τ i R x I x   s + κ i 1 L 2 I x C x ( π 2 ( i     0.5 ) 2   +   ( G x L 2 ) 2 ) ,
The expressions for the low frequency correcting factors remain the same except for C 2 , whose expression becomes
C 2 = { [ i = 1 m 2 τ i κ i π 2 ( i     0.5 ) 2 ] 1 i f   θ = 0 t a n h ( G x L 2 ) G X L [ i = 1 m τ i κ i π 2 ( i     0.5 ) 2   +   ( G x L 2 ) 2 ] 1 i f   θ 0
The frequency dependent modifications factors τ i and κ i are computed by adjusting the approximated transfer functions of Equation (8) using the dissipative model, which is the exact solution of horizontal laminar flow dynamics problem [17]. Shown in Figure 4 are the resulting coefficients for liquid and air cases as functions of the dimensionless characteristic root λ c i that is defined as
λ c i = 8   ( i 0.5 ) L R x I x C x
The analytical expressions of these modification factors are obtained as
τ i = A 1 t a n h ( A 2 l o g 10 ( λ c i ) + A 3 ) + A 4 t a n h ( A 5 l o g 10 ( λ c i ) + A 6 ) + A 7 ,
l o g 10 ( κ i ) = B 1 t a n h ( B 2 l o g 10 ( λ c i ) + B 3 ) + B 4 t a n h ( B 5 l o g 10 ( λ c i ) + B 6 ) + B 7
where the parameters A i and B i are given in Table 3 for both liquid and air cases [17]. Similar expressions can be derived for other gases.
The seconds-order effects on the transmission line frequency response are illustrated in Figure 5 for the case of laminar flow in a vertical transmission line given the properties shown in Table 2. It is seen that second-order effects result in much greater attenuation of mid- and high frequency components of the proposed model in Equation (47), compared to the linear friction model in Equation (8).
The accuracy of the proposed model is investigated by comparing its frequency response with that of the modified dissipative model developed in [15]. The dissipative model is known in the literature as the exact model for laminar flow in horizontal transmission lines. The model is extended to account for compressibility effects in inclined fluid lines and an exact analytical solution represented as a four-terminal network is derived. A comparison between the frequency responses of the model developed in and the proposed model of Equation (47) is depicted in Figure 6, showing a good agreement for different line inclinations.
The proposed model and the model developed in [15] are solutions of linearized ordinary differential equations that are derived by assuming laminar flow conditions. To extend the model to account for turbulent flow conditions, a similar technique to that presented in [14] is adopted. An effective resistance term is introduced as
R e f f = δ ( R t u r R l a m ) R ,
where δ = 0 and δ = 1 for the case of laminar flow and turbulent flow, respectively. The resistance terms are defined as
R t a m = λ l a m π   μ   R e 8   A 2   L ,
R t u r = λ t u r π   μ   R e 8   A 2   L ,
such that λ l a m and λ t u r are the friction factors for laminar flow and turbulent flow, respectively. The proposed model that can be applied to both laminar and turbulent flow conditions in inclined transmission lines is given as
[ P i n ( s ) Q o u t ( s ) ] = [ T F ˜ 11 ( s ) 1   +   R e f f T F ˜ 21 ( s )   T F ˜ 12 ( s )   +   R e f f 1   +   R e f f T F ˜ 21 ( s ) T F ˜ 21 ( s ) 1   +   R e f f T F ˜ 21 ( s ) T F ˜ 22 ( s ) 1   +   R e f f T F ˜ 21 ( s ) ] [ P o u t ( s ) Q i n ( s ) ]

4. Dynamic Modeling of Transient Two Phase Flow

As described in Figure 7, to model transient two-phase flow behavior, a steady-state model is first used to determine the flow pattern and estimate the liquid holdup and pressure drop quantities. Then, the equivalent fluid properties (density, viscosity, and bulk modulus) of the mixture are estimated using the results presented in [14], and finally used as model parameter inputs for the proposed analytical dynamic and single-phase model.
A mechanistic model is used to estimate pressure gradient (drop) and liquid holdup, given gas and liquid configuration in horizontal and inclined transmission lines [14]. The model was derived from fundamental laws of physics. An experimental database was used to calibrate the model unknown parameters. The model distinguishes between several flow patterns (regimes) such as bubble flow, dispersed bubble flow, stratified flow, annular-mist flow, intermittent flow, and forth flow. For flow regime identification, the steady-state two-phase flow model developed in [15] assumes the existence of a particular flow pattern and then examines various criteria that establish stable boundaries of the flow regime. A new flow pattern is assumed, and the procedure is repeated if the regime is found to be unstable. Once the flow pattern is obtained, fluid dynamic relationships are then used to estimate the pressure gradient and liquid holdup variables based on fluid mixture, phase material properties, and transmission line inclination. The procedure and corresponding mathematical relationships for flow pattern determination and pressure drop/liquid holdup calculation are detailed in [7,14].

4.1. Equivalent Fluid Parameters

Based on the steady-state pressure drop and liquid holdup obtained by the mechanistic model, the equivalent fluid phase parameters, namely, the bulk modulus, density, speed of sound, and viscosity, are derived using the gas and liquid single-phase properties [14]. The two-phase flow is then assumed to be one homogeneous mixture fluid, and the proposed analytical dynamic and single-phase model in Equation (47) is applied.
The bulk modulus of a fluid describes the material elasticity as it experiences a volumetric deformation, which is defined as the ratio of an infinitesimal pressure stress to its resulting relative change of the volume strain. However, in the present research for the liquid-gas two-phase flow, the existence of air in the oil flow can considerably reduce its effective bulk modulus. Thus, considering such deviation in two-phase flow, this mixture bulk modulus was determined as a combination of the liquid and the gas bulk moduli in parallel. The equivalent bulk modulus of the liquid-gas two-phase mixture can be expressed as
1 β e q = E L β L + 1 E L β G ,
where β L is the liquid bulk modulus, β G is the gas bulk modulus, and E L is the predicted liquid holdup given by the steady-state model. While the flow GVF is increasing, the equivalent bulk modulus of the two-phase mixture will decrease accordingly, based on Equation (56), which causes lower frequencies of oscillation in flow applications.
The equivalent density of the two-phase fluid is calculated as follows
ρ e q = E L ρ L + ( 1 E L ) ρ G ,
where ρ L and ρ G are respectively the liquid and gas density. Assuming a constant temperature inside the pipeline, the density of the equivalent fluid will be a function of the pressure only
ρ e q ( P ) = ρ 0 1 + 1 β e q ( P 0 P )
hence, the equivalent speed of sound in the fluid is obtained as
c e q = β e q ρ e q
For a liquid-gas two-phase flow problem, higher GVF will cause lower mixture density and its equivalent bulk modulus. Considering that the fluid sound speed can be more sensitive to its bulk modulus compared to its density, such higher GVF may lead to a lower sound speed when an equivalent fluid is characterized [14]. For such equivalent mixture flow, the basic Darcy friction factor is calculated by
f e q = 2 D A 2 Δ P s s ρ e q Q 2 ,
where D is the pipe diameter, A is the pipe cross-sectional area, Q is the flow rate, and Δ P s s is the steady-state predicted pressure drop. Then, the equivalent viscosity of the mixture fluid is derived based on this Darcy friction factor. In the case of laminar flow, the equivalent dynamic viscosity will be given by
μ e q = 1 64 ρ e q V m D f e q
If the flow is turbulent, the Colebrook equation [20] is used to recover the equivalent viscosity. Other correlations such as Goudar and Sonnad [21] can also be used. Once all the equivalent single-phase fluid properties are determined, the dynamic model proposed in Equation (47) is used to simulate transient behavior of the two-phase flow in inclined transmission lines.

4.2. Model Comparison with OLGA

To investigate the model accuracy under transient multiphase flow conditions, a time-domain analysis is performed by simulating an air/water flow in inclined transmission line using OLGA under different GVF conditions as described given Table 4.
Shown in Figure 8 are comparisons between a step-time response from OLGA and that from the proposed model using one mode approximation for the case of 10% GVF and for different inclination angles. The inlet flow rate is stepped as shown in Table 4 and the outlet pressure is kept constant at 100 bar. An agreement is observed between the two responses, indicating a match between the estimated system fundamental frequencies and damping ratios. Indeed, a small error is found with a mean absolute error (MAE) of less than 4% between the estimated overshoots and a MAE of less than 0.5 s between the estimated settling times. The same flow conditions are simulated for a fixed inclination angle of 20° upward and various GVF values. The inlet pressure transient responses are depicted in Figure 9. Agreement between the two models is noted, especially for lower GVF values. As GVF increases, more discrepancies between the results of the proposed model and those obtained from OLGA are observed. The MAE between the estimated overshoots for the 80% GVF case is around 8% and that between the estimated settling times is around 2 s.
As described previously, the accuracy of the fluid resonant frequencies predicted by the proposed model is precise and not a function of number of modes used in the approximation. Increasing the model order (complexity) is solely a function of the desired frequency range. Illustrated in Figure 10 is a comparison of inlet pressure step responses between the proposed model and OLGA for the 10% GVF flow case in an inclined line with 20° downward angle. It can be concluded that, unlike OLGA, the proposed model enables the estimation of high-order dynamic effects by increasing the number of approximated second-order modes.

5. Conclusions

Presented is an analytical model for transient laminar and turbulent two-phase flow in inclined transmission lines. Four transcendent transfer functions are obtained by solving the linearized single-phase Navier-Stokes equations. To enable system analysis and synthesis in the time domain, these transfer functions are approximated by finite-order rational polynomial transfer functions using residue theorem. A major benefit of the proposed model over existing models is that its coefficients are explicit functions of the transmission line and fluid properties rather than numerically derived from tables or graphs. The proposed dynamic single-phase model is then coupled with a mechanistic model for steady-state pressure drop and liquid holdup estimation of liquid-gas flow conditions. Model ability to capture the essential dynamics of two-phase flow in transmission lines is investigated upon comparison with OLGA under different GVF conditions and inclination angles, showing a good agreement of steady-state characteristics and transmission line natural frequency response. Discrepancies have been seen in the step response overshoots and settling times. However, for the case studies described in this paper, the mean absolute error between the estimated overshoots and settling times is less than 8% and 2 s, respectively.

Author Contributions

Conceptualization, M.A.F. and T.W.; methodology, T.W.; software, H.M., T.W., and Y.T.; validation, T.W., H.M., and Y.T.; writing—original draft preparation, T.W.; writing—review and editing, T.W., M.A.F., H.M., and Y.T.; supervision, M.A.F. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Bansal, P.K.; Rupasinghe, A.S. A Homogeneous Model for Adiabatic Capillary Tubes. Appl. Therm. Eng. 1998, 18, 207–219. [Google Scholar] [CrossRef]
  2. Awad, M.M.; Muzychka, Y.S. Effective Property Models for Homogeneous Two-Phase Flows. Exp. Therm. Fluid Sci. 2008, 33, 106–113. [Google Scholar] [CrossRef]
  3. Richter, H.J. Separated Two-Phase Flow Model: Application to Critical Two-Phase Flow. Int. J. Multiph. Flow 1983, 9, 511–530. [Google Scholar] [CrossRef]
  4. Stevanovic, V.; Prica, S.; Maslovaric, B. Multi-Fluid Model Predictions of Gas-Liquid Two-Phase Flows in Vertical Tubes. FME Trans. 2007, 35, 173–181. [Google Scholar]
  5. Choi, J.; Pereyra, E.; Sarica, C.; Park, C.; Kang, J.M. An Efficient Drift-Flux Closure Relationship to Estimate Liquid Holdups of Gas-Liquid Two-Phase Flow in Pipes. Energies 2012, 5, 5294–5306. [Google Scholar] [CrossRef] [Green Version]
  6. Coddington, P.; Macian, R. A Study of the Performance of Void Fraction Correlations Used in the Context of Drift-Flux Two-Phase Flow Models. Nucl. Eng. Des. 2002, 215, 199–216. [Google Scholar] [CrossRef]
  7. Taitel, Y.; Dukler, A.E. A Model for Predicting Flow Regime Transitions in Horizontal and Near Horizontal Gas-Liquid Flow. AIChe J. 1976, 22, 22–24. [Google Scholar] [CrossRef]
  8. Barnea, D. A Unified Model for Predicting Flow Pattern Transitions for the Whole Range of Pipe Inclinations. Int. J. Multiph. Flow 1987, 13, 1–12. [Google Scholar] [CrossRef]
  9. Ansari, A.M.; Sylvester, A.D.; Sarica, C.; Shoham, O.; Brill, J.P. A Comprehensive Mechanistic Model for Upward Two-Phase Flow in Wellbores. SPE Prod. Facil. 1994, 9, 143–151. [Google Scholar] [CrossRef] [Green Version]
  10. Petalas, N.; Aziz, K. Development and Testing of a New Mechanistic Model for Multiphase Flow in Pipes. In Proceedings of the ASME Fluids Engineering Division Summer Meeting, Part 1(of 2), San Diego, CA, USA, 7 July 1996; pp. 153–159. [Google Scholar]
  11. Petalas, N.; Aziz, K. A Mechanistic Model for Multiphase Flow in Pipes. J. Can. Pet. Technol. 2000, 39, 43–55. [Google Scholar] [CrossRef]
  12. Bendiksen, K.; Malnes, D.; Moe, R.; Nuland, S. The Dynamic Two-Fluid Model OLGA: Theory and Application. SPE Prod. Eng. 1991, 6, 171–180. [Google Scholar] [CrossRef]
  13. Pauchon, C.; Dhulesia, H.; Lopez, D.; Fabre, J. TACITE: A Comprehensive Mechanistic Model for Two-Phase Flow. In Proceedings of the 6th International Conference on Multiphase Production, Cannes, France, 16–18 June 1993. [Google Scholar]
  14. Meziou, A.; Chaari, M.; Franchek, M.; Borji, R.; Grigoriadis, K.; Tafreshi, R. Low-Dimensional Modeling of Transient Two-Phase Flow in Pipelines. J. Dyn. Syst. Meas. Control 2016, 138, 101008. [Google Scholar] [CrossRef]
  15. Omrani, A.; Franchek, M.A.; Grigoriadis, K. Transmission Line Modeling of Inclined Compressible Fluid Flows. ASME J. Dyn. Syst. Meas. Control 2018, 140, 011001. [Google Scholar] [CrossRef]
  16. Chaari, M.; Fekih, A.; Seibi, A.C.; Ben Hmida, J. A Generalized Reduced-Order Dynamic Model for Two-Phase Flow in Pipes. ASME J. Fluids Eng. 2019, 141, 101303. [Google Scholar] [CrossRef]
  17. Wassar, T.; Franchek, M.; Gutierrez, J. Reduced-order Modelling of Transient Flow in Transmission Lines Using Distributed Lumped Parameters. Int. J. Fluid Power 2017, 18, 153–166. [Google Scholar] [CrossRef]
  18. Matko, D.; Geiger, G.; Werener, T. Modelling of The Pipeline as a Lumped Parameter System. Automatika 2001, 42, 177–188. [Google Scholar]
  19. Yang, W.C.; Tobler, W.E. Dissipative Modal Approximation of Fluid Transmission Line Using Linear Friction Model. J. Dyn. Syst. Meas. Control Trans. ASME 1991, 113, 152–162. [Google Scholar] [CrossRef]
  20. Colebrook, C.F. Turbulent Flow in Pipes, with Particular Reference to the Transition Region between Smooth and Rough Pipe Laws. J. Inst. Civ. Eng. 1939, 11, 133–156. [Google Scholar] [CrossRef]
  21. Goudar, C.T.; Sonnad, J.R. Comparison of the Iterative Approximations of the Colebrook-White Equation. Hidrocarbon Process. 2008, 87, 79–83. [Google Scholar]
Figure 1. Transfer functions sensitivity to inclination angle.
Figure 1. Transfer functions sensitivity to inclination angle.
Fluids 06 00300 g001
Figure 2. Transcendental transfer functions approximation.
Figure 2. Transcendental transfer functions approximation.
Fluids 06 00300 g002
Figure 3. Modal approximation for under-damped and well-damped cases.
Figure 3. Modal approximation for under-damped and well-damped cases.
Fluids 06 00300 g003
Figure 4. Frequency dependent modification factors (reproduced from [17]).
Figure 4. Frequency dependent modification factors (reproduced from [17]).
Fluids 06 00300 g004
Figure 5. Frequency dependent effects on line frequency response.
Figure 5. Frequency dependent effects on line frequency response.
Fluids 06 00300 g005
Figure 6. Comparison between proposed model and Omrani et al. model [15].
Figure 6. Comparison between proposed model and Omrani et al. model [15].
Fluids 06 00300 g006
Figure 7. Multiphase flow modeling procedure.
Figure 7. Multiphase flow modeling procedure.
Fluids 06 00300 g007
Figure 8. Comparison of dynamic inlet pressure for 10% GVF.
Figure 8. Comparison of dynamic inlet pressure for 10% GVF.
Fluids 06 00300 g008
Figure 9. Comparison of dynamic inlet pressure for 20° inclination.
Figure 9. Comparison of dynamic inlet pressure for 20° inclination.
Fluids 06 00300 g009
Figure 10. Number of modes effect on pressure transient response.
Figure 10. Number of modes effect on pressure transient response.
Fluids 06 00300 g010
Table 1. Transmission line distributed parameter.
Table 1. Transmission line distributed parameter.
R x I x C x G x
8   π μ A 2 ρ ¯ A A c ¯ 2 ρ ¯ g   sin ( θ ) c ¯ 2
Table 2. Study case parameters.
Table 2. Study case parameters.
ParameterValueUnit
L500m
A4.9 × 10−4m2
ρ ¯ 500kg/m3
μ 6.6 × 10−4kg/m/s
β e 106Pa
Table 3. Values of the parameters Ak and Bk.
Table 3. Values of the parameters Ak and Bk.
LiquidAir
A10.9202.050
A21.5001.300
A3−2.100−2.100
A40.8700.800
A52.0002.530
A6−4.320−5.780
A72.4703.610
B10.1000.275
B21.1430.665
B3−1.147−0.649
B4−0.127−0.100
B50.7160.702
B60.455−1.323
B70.000−0.213
Table 4. OLGA simulation parameters.
Table 4. OLGA simulation parameters.
Transmission Line
LengthDiameterInclinationRoughness
100 m0.038 m0°, 30°, 60°, 90°3.3 × 10−6 m
Liquid
DensitySpeed of soundViscosity
884 kg/m31372 m/s0.075 Pa s
Gas
DensitySpeed of soundViscosity
82 kg/m3400 m/s1.5 × 10−5 Pa s
GVFFlow rate
10%, 25%, 50%Stepped from 1.8 × 10−3 to 0.002 m3/s
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Wassar, T.; Franchek, M.A.; Mnasri, H.; Tang, Y. An Explicit Analytical Solution for Transient Two-Phase Flow in Inclined Fluid Transmission Lines. Fluids 2021, 6, 300. https://doi.org/10.3390/fluids6090300

AMA Style

Wassar T, Franchek MA, Mnasri H, Tang Y. An Explicit Analytical Solution for Transient Two-Phase Flow in Inclined Fluid Transmission Lines. Fluids. 2021; 6(9):300. https://doi.org/10.3390/fluids6090300

Chicago/Turabian Style

Wassar, Taoufik, Matthew A. Franchek, Hamdi Mnasri, and Yingjie Tang. 2021. "An Explicit Analytical Solution for Transient Two-Phase Flow in Inclined Fluid Transmission Lines" Fluids 6, no. 9: 300. https://doi.org/10.3390/fluids6090300

APA Style

Wassar, T., Franchek, M. A., Mnasri, H., & Tang, Y. (2021). An Explicit Analytical Solution for Transient Two-Phase Flow in Inclined Fluid Transmission Lines. Fluids, 6(9), 300. https://doi.org/10.3390/fluids6090300

Article Metrics

Back to TopTop