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
where
and
are independent cylindrical space variables, and
is the independent time variable. The dependent fluid variables include pressure
, its time average
, velocities in
x-direction and
r-direction (
and
), temperature
, its time average
, density
, its time average
, fluid absolute viscosity
, fluid effective bulk modulus
, 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
where
is the Laplace transform variable, and
is the flow rate and pressure along the transmission line. The distributed lumped parameters in Equations (6) and (7) are
,
, and
defined as the resistance, inertance, capacitance, and gravity term per unit length of transmission line and are presented in
Table 1, where
is the time averaged speed of sound in the fluid,
is the gravity acceleration, and
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
where
and
Illustrated in
Figure 1 is the sensitivity of the matrix transfer functions
,
,
, and
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 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 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
Impedance Representation: Admittance Representation: 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
where
is the liquid bulk modulus,
is the gas bulk modulus, and
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
where
and
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
hence, the equivalent speed of sound in the fluid is obtained as
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
where
is the pipe diameter,
is the pipe cross-sectional area,
is the flow rate, and
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
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.