Modeling and Simulation for Steady State and Transient Pipe Flow of Condensate Gas

1. Introduction

Condensate gas is mainly demonstrated by methane. However, it also contains a lot of heavier contents like C₅ or C₅+ and some non-hydrocarbon mixture as well (Mokhatab et al, 2006). After recovering from gas wells, condensate gas needs liquid separation, gas purification and condensate stabilization treatment in the processing plant to meet the quality requirements. Processing plants far away from the gas well with long distances of two-phase flow in one condensate gas pipeline will take less investment than adjacent process plant with two single phase pipelines which are dry gas pipeline and liquid phase pipeline (Li, 2008).

If the operation temperature somewhere in the condensate gas pipeline is lower than the gas dew point, liquid condensation would occur, subjecting the pipeline to two phase flow (Potocnik, 2010). While gas and its condensate flow simultaneously, mass transfer takes place continuously due to the change in pressure and temperature conditions. This leads to compositional changes and associated fluid property changes and also makes the hydraulic and thermal calculations of condensate gas more complex than normal gas. The condensate gas pipeline model which is established and solved based on the principle of fluid mechanics can simulate hydraulic and thermal parameters under various operation conditions. By means of technical support, this model is of great importance in the pipeline design and safety operation aspects (Mokhatab, 2009).

2. Thermodynamic model

The purpose of the thermodynamic model is three-fold. First, it defines the transition between single phase/two phase conditions (point of condensate inception in the pipeline or gas dew point). Second, it is used for the prediction of properties for the flowing fluids (gas and its condensate). And lastly, it derives the mass exchange between the flowing phases (Adewumi et al, 1990 ; Estela-Uribe et al, 2003). This work uses the BWRS equation of state (EOS) to implement the thermodynamic model as it has proven reliable for gas condensate system (McCain et al, 1990). Most property predictions are derived from the equation of state (i.e., densities values, densities values and their derivatives with respect to pressure and temperature, departure enthalpies, heat capacities, and Joule-Thompson coefficients). Additionally, phase equalibria are calculated on the basis of flash calculation method. Expressions for such parameters as fugacity are elaborated in standard textbooks, where the theory and relevant procedures for flash calculation are well documented (API, 2005).

3. General model of condensate gas pipeline

3.1. Basic assumptions

One of the most fundamental approaches used to model two phase flow is the two-fluid model (Ayala et al, 2003). It consists of separate mass, momentum and energy conservation equations written for each of the phases. This results in up to six differential equations. Furthermore, the model is more complicate as relating parameters in the equations vary with the fluid flowing. Thus, simple and reasonable assumptions could help to reduce the unknowns in the model (Hasan et al, 1992). The hydrodynamic model has three major inherent assumptions:

1. Gas and liquid average flow velocities are calculated according to the section area occupied by each phase respectively.

2. Mass transfer takes place between gas and liquid phase. However, two phases are assumed to be at thermodynamic equilibrium at every point within the pipe.

3. In the transient process, the pipeline assumed to be isothermal. The liquid hold up within pipeline, the wall friction force of each phase, the drag force between two phases are assumed to be the same as steady state.

3.2. Description of general model

The general model for two phase flow is built on the basis of mass, momentum and energy conservation equations of each phase (Li et al, 2009; Ayala et al, 2003).

Gas phase continuity equation:

$$\frac{\partial }{\partial t}\left({\rho }_g\phi A\right)+\frac{\partial }{\partial x}\left({\rho }_g\phi A{w}_g\right)=\Delta {\dot{m}}_{gl}$$E1

Liquid phase continuity equation:

$$\frac{\partial }{\partial t}\left({\rho }_l{H}_{L}A\right)+\frac{\partial }{\partial x}\left({\rho }_l{H}_{L}A{w}_l\right)=\Delta {\dot{m}}_{lg}$$E2

Where, $$\Delta {\dot{m}}_{gl}$$is mass rate of phase change from gas to liquid, (kg/s.m); $$\Delta {\dot{m}}_{lg}$$is mass rate of phase change from liquid to gas, (kg/s.m). H_L is liquid hold up; $\phi$is gas hold up; ..${\rho }_g$..is density of the gas phase, kg/m³; ${\rho }_l$is density of the liquid phase, kg/m³; A is pipe cross sectional area, m²; t is time, s; x is length along the pipe length, m; w_g is velocity of the gas phase, m/s; w_l is velocity of the liquid phase, m/s.

$$\Delta {\dot{m}}_{gl}$$is defined as (3).

$$\Delta {\dot{m}}_{gl}=\left[{\left(\frac{\partial Y_s}{\partial P}\right)}_T\frac{\partial P}{\partial x}+{\left(\frac{\partial Y_s}{\partial P}\right)}_T\frac{\partial P}{\partial t}\frac{\partial t}{\partial x}\right]\cdot \left(m_g+m_l\right)+\left[{\left(\frac{\partial Y_s}{\partial T}\right)}_P\frac{\partial T}{\partial x}+{\left(\frac{\partial Y_s}{\partial T}\right)}_P\frac{\partial T}{\partial t}\frac{\partial t}{\partial x}\right]\cdot \left(m_g+m_l\right)$$E3

$$Y_s=\frac{m_g}{m_g+m_l}$$E4

.

$$m_g={\rho }_g{w}_g\phi A$$E5

$$m_l={\rho }_l{w}_g{H}_{L}A$$E6

Where, P is pressure, Pa; T is temperature, K; m_g is gas mass fraction in two-phase fluid system; m_l is liquid mass fraction in two-phase fluid system; Y_s is mass fraction of gas.

Gas phase momentum equation.

$$\frac{\partial }{\partial t}\left({\rho }_g\phi w_{g}A\right)+\frac{\partial }{\partial x}\left({\rho }_g\phi w_g^{2}A\right)+A_g\frac{\partial P}{\partial x}=\Delta {\dot{m}}_{gl}{w}_a-F_{gw}-F_{gi}-{\rho }_{g}g\phi A\sin \theta$$E7

Where, A_g is pipe cross sectional area occupied by gas phase, m²;F_gw is wall shear force of gas phase, N/m³; F_gi is interfacial drag force on gas phase, N/m³; g is acceleration of gravity, m/s²; θ is pipeline slope, rad.

Liquid phase momentum equation.

$$\frac{\partial }{\partial t}\left({\rho }_l{H}_L{w}_{l}A\right)+\frac{\partial }{\partial x}\left({\rho }_l{H}_L{w}_l^{2}A\right)+A_l\frac{\partial P}{\partial x}=\Delta {\dot{m}}_{lg}{w}_a-F_{lw}-F_{li}-{\rho }_{l}g{H}_{L}A\sin \theta$$E8

Where, A_l is pipe cross sectional area occupied by liquid phase, m²; F_lw is wall shear force of liquid phase, N/m³; F_li is interfacial drag force on liquid phase, N/m³; w_a is transition velocity between the gas phase and liquid phase, m/s.

Gas-liquid phase mixture energy equation

$$\begin{array}{l}\frac{\partial }{\partial t}\left[A{\rho }_g\phi \left(h_g+\frac{w_g^2}{2}+gz\right)+A{\rho }_l{H}_L\left(h_l+\frac{w_l^2}{2}+gz\right)\right]\\ +\frac{\partial }{\partial x}\left[A{\rho }_g\phi w_g\left(h_g+\frac{w_g^2}{2}+gz\right)+A{\rho }_l{H}_L{w}_l\left(h_l+\frac{w_l^2}{2}+gz\right)\right]+K\pi D\left(T-T_0\right)=0\end{array}$$E9

Where, K is overall heat transfer co-efficiency, W/(m² K); D is external diameter, m; T₀ is environmental temperature, K; h_l is enthalpy of liquid phase in pipeline, kJ/kg; h_g is enthalpy of the gas phase in pipeline, kJ/kg; z is pipeline elevation, m.

Add (7) and (8), obtain:

$$\begin{array}{l}\frac{\partial }{\partial t}\left({\rho }_g\phi w_{g}A\right)+\frac{\partial }{\partial t}\left({\rho }_l{H}_L{w}_{l}A\right)+\frac{\partial }{\partial x}\left({\rho }_g\phi w_g^{2}A\right)+\frac{\partial }{\partial x}\left({\rho }_l{H}_L{w}_l^{2}A\right)+A\frac{\partial P}{\partial x}\\ =-F_{gw}-F_{lw}-\left({\rho }_g\phi +{\rho }_l{H}_L\right)gA\sin \theta \end{array}$$E10

Equation (10) is transformed into (11) through eliminating the pressure terms.

$$\begin{array}{l}\frac{\partial }{\partial t}\left({\rho }_g{w}_g\right)-\frac{\partial }{\partial t}\left({\rho }_l{w}_l\right)+\frac{\partial }{\partial x}\left({\rho }_g{w}_g^2\right)-\frac{\partial }{\partial x}\left({\rho }_l{w}_l^2\right)=\\ \Delta {\dot{m}}_{gl}{w}_a\left(\frac{1}{A_g}+\frac{1}{A_l}\right)-\frac{F_{gw}}{A_g}-\frac{F_{gi}}{A_g}+\frac{F_{lw}}{A_l}+\frac{F_{li}}{A_l}-\left({\rho }_l{H}_L-{\rho }_g\phi \right)g\sin \theta \end{array}$$E11

In view of the slow transient behaviour in condensate gas pipeline, we can obtain equation (12) by ignoring velocity variation terms in equation (11) (Li et al, 1998):

$$-\frac{F_{gw}}{A_g}-F_{gi}\left(\frac{1}{A_g}+\frac{1}{A_l}\right)+\frac{F_{lw}}{A_l}-\left({\rho }_l{H}_L-{\rho }_g\phi \right)g\sin \theta =0$$E12

Equations (1), (2), (9), (10), (12) construct the basic model for condensate gas pipe flow simulation.

4. Constitutive equations

The condensate gas flow model is one dimensional two-fluid multiphase hydrodynamic model which adapts to different flow patterns in pipeline. According to Cindric and Shoham, the flow patterns in horizontal pipeline are stratified flow, intermittent flow, annular flow, dispersed flow and these in vertical pipeline are bubble flow, slug flow, churn flow, annular flow (Mokhatab et al, 2006). Because of the constitutive equations is dependent on the flow pattern, one of the greatest difficulties in the analysis of two-phase flow in pipeline is defining appropriate constitutive equations for relating relevant forces-such as the steady drag force and interfacial force.

Considering low liquid hold up, the flow pattern in the condensate gas pipeline is stratified flow which has explicit interface between the liquid and gas phase, as depicted in Fig.1. Then, we can obtain the calculation methods of unknowns which dependent on the constitutive equations (Taitel et al, 1995 ; Chen et al, 1997 ).

Figure 1.

Stratified flow in condensate gas pipeline

Wall shear force of each phase is expressed as follow.

$$F_{kw}=-{\tau }_{kw}{S}_k$$E13

Where, k=g when the equation is applied for the gas phase; k=l when the equation is applied for the liquid phase; S_k is defined as follow:

$$S_k\text{=}\frac{\text{total}\text{wall}\text{area}\text{wetted}\text{by}\text{phase}k}{\text{Totoal volume}}$$E14

${\tau }_{kw}$is defined as:

$${\tau }_{kw}=\frac{1}{2}{\lambda }_k{\rho }_k\left|w_k\right|w_k$$E15

In which,$\lambda$ is Fanning factor which is calculated by Colebrook & White empirical correlation.

If Reynolds Number Re≤ 2000

$$\lambda =\frac{16}{Re}$$E16

If Re≥ 2000

$$\frac{1}{\sqrt{\lambda }}=3.48-4\lg \left(\frac{2\epsilon }{D}+\frac{9.35}{Re\sqrt{\lambda }}\right)$$E17

Where, $\epsilon$is absolute roughness of pipeline wall, m.

Interfacial force between phases is defined as follow:

$$F_{gi}=-F_{li}=-{\tau }_i{S}_i$$E18

Where

$$S_i\text{=}\frac{\text{Totalsurfaceareaofcontactbetweenphases}}{\text{Totoalvolume}}$$E19

$${\tau }_i=\frac{1}{2}{\lambda }_i{\rho }_g\left|w_g-w_l\right|\left(w_g-w_l\right)$$E20

The interfacial friction factor ${\lambda }_i$ is calculated with Hanrrity correlation.If $w_{sg}\le w_{sg\cdot t}$

$${\lambda }_i={\lambda }_g$$E21

If $w_{sg}\ge w_{sg\cdot t}$

$${\lambda }_i={\lambda }_g\left[1+15\sqrt{\frac{h_l}{D}}\left(\frac{w_{sg}}{w_{sg\cdot t}}-1\right)\right]$$E22

Where,

$$w_{sg}=\frac{Q_G}{A}$$E23

$$w_{sg\cdot t}=5\sqrt{\frac{101325}{P}}$$E24

Where, w_sg is reduced velocity of the gas phase, m/s; w_sg._t is reduced velocity for indentifying the transition from stratified flow pattern to smooth stratified flow pattern, m/s; Q_G is flow rate of the gas phase, m³/s.

5. Steady state analysis of condensate gas pipeline

5.1. Basic equations

While steady operation, the variation of each parameter in equations (1), (2), (9), (10) with time can be ignored. Expand the equations above and the following equations used for steady state simulation can be obtained (Li et al, 2009):

Gas phase continuity equation:

$$A{w}_g\phi {(\frac{\partial {\rho }_g}{\partial P})}_T\frac{dP}{dx}+A{w}_g\phi {(\frac{\partial {\rho }_g}{\partial T})}_P\frac{dT}{dx}+A{\rho }_g\phi \frac{d{w}_g}{dx}+A{\rho }_g{w}_g\frac{d\phi }{dx}=\Delta {\dot{m}}_{gl}$$E25

Liquid phase continuity equation:

$$A{w}_l{H}_L{(\frac{\partial {\rho }_l}{\partial P})}_T\frac{dP}{dx}+A{w}_l{H}_L{(\frac{\partial {\rho }_l}{\partial T})}_P\frac{dT}{dx}+A{\rho }_l{H}_L\frac{d{w}_l}{dx}+A{\rho }_l{w}_l\frac{d{H}_L}{dx}=\Delta {\dot{m}}_{lg}$$E26

Gas-liquid phase mixture momentum equation:

$$A\frac{dP}{dx}+{\rho }_g{w}_g\phi A\frac{d{w}_g}{dx}+{\rho }_l{w}_l{H}_{L}A\frac{d{w}_l}{dx}=-F_{gw}-F_{lw}-\left({\rho }_g\phi +{\rho }_l{H}_L\right)gA\sin \theta$$E27

Gas-liquid phase mixture energy equation:

$$\begin{array}{l}\left[{\rho }_g{w}_g\phi A{\left(\frac{\partial h_g}{\partial P}\right)}_T+{\rho }_l{w}_l{H}_{L}A{\left(\frac{\partial h_l}{\partial P}\right)}_T\right]\frac{dP}{dx}\\ +\left[{\rho }_g{w}_g\phi A{\left(\frac{\partial h_g}{\partial T}\right)}_P+{\rho }_l{w}_l{H}_{L}A{\left(\frac{\partial h_l}{\partial T}\right)}_P\right]\frac{dT}{dx}+{\rho }_g{w}_g^2\phi A\frac{d{w}_g}{dx}+{\rho }_l{w}_l^2{H}_{L}A\frac{d{w}_l}{dx}\\ =-\left({\rho }_{g}g{w}_g\phi A+{\rho }_{l}g{w}_l{H}_{L}A\right)\frac{dz}{dx}-K\pi D\left(T-T_0\right)-\Delta {\dot{m}}_{gl}\left(h_g-h_l+\frac{w_g^2-w_l^2}{2}\right)\end{array}$$E28

The system of simultaneous differential equations composed of (25)-(28) can be written in their non-conservative form.

$$A\frac{dU}{dx}=D$$E29

Where $$A=\left[\begin{array}{cccc}{a}_{11}& a_{12}& a_{13}& a_{14}\\ a_{21}& a_{22}& a_{23}& a_{24}\\ a_{31}& a_{32}& a_{33}& a_{34}\\ a_{41}& a_{42}& a_{43}& a_{44}\end{array}\right]$$,$$D=\left[\begin{array}{c}{D}_1\\ D_2\\ D_3\\ D_4\end{array}\right]$$ ,$$U=\left[\begin{array}{c}P\\ T\\ w_g\\ w_l\end{array}\right]$$ ,$$a_{11}=A{w}_g\phi {(\frac{\partial {\rho }_g}{\partial P})}_T$$,$$a_{12}=A{w}_g\phi {(\frac{\partial {\rho }_g}{\partial T})}_P$$ ,$$a_{13}=A{\rho }_g\varphi$$ ,$$a_{14}=0$$ ,$$a_{21}=A{w}_l{H}_L{(\frac{\partial {\rho }_l}{\partial P})}_T$$ ,$$a_{22}=A{w}_l{H}_L{(\frac{\partial {\rho }_l}{\partial T})}_P$$,$$a_{23}=0$$ ,$$a_{24}=A{\rho }_l{H}_L$$ ,$$a_{31}=A$$ ,$$a_{32}=0$$ ,$$a_{33}={\rho }_g{w}_g\phi A$$ ,$$a_{34}={\rho }_l{w}_l{H}_{L}A$$ ,$$a_{41}={\rho }_g{w}_g\phi A{(\frac{\partial h_g}{\partial P})}_T+{\rho }_l{w}_l{H}_{L}A{(\frac{\partial h_l}{\partial P})}_T$$$$a_{42}={\rho }_g{w}_g\phi A{(\frac{\partial h_g}{\partial T})}_P+{\rho }_l{w}_l{H}_{L}A{(\frac{\partial h_l}{\partial T})}_P$$,$$a_{43}={\rho }_g{w}_g^2\phi A$$ ,$$a_{44}={\rho }_l{w}_l^2{H}_{L}A$$ $$D_1=\Delta {\dot{m}}_{gl}-A{\rho }_g{w}_g\frac{d\phi }{dx}$$,$$D_2=\Delta {\dot{m}}_{lg}-A{\rho }_l{w}_l\frac{d{H}_L}{dx}$$ ,$$D_3=-F_{gw}-F_{lw}-\left({\rho }_g\phi +{\rho }_l{H}_L\right)gA\sin \theta$$ $$D_4=-({\rho }_{g}g{w}_g\phi A+{\rho }_{l}g{w}_l{H}_{L}A)\frac{dz}{dx}-K\pi D(T-T_0)-\Delta {\dot{m}}_{gl}(h_g-h_l+\frac{w_g^2-w_l^2}{2})$$

5.2. Model solving

Steady state condensate gas model is formed by 5 equations which are (29) and (12). There are five unknowns, liquid holdup (H_L), pressure (P), temperature (T), gas and liquid velocity (w_g and w_l), in the model. Thus, the closure of the model is satisfied.

To solve the model, the liquid hold up is obtained by solving (12) firstly. And then, pressure (P), temperature (T), velocity of the gas phase (w_g), and the velocity of the liquid phase (w_l) are obtained by solving (29). The procedures for solving (29) are presented in details as follow:

1. Suppose the pipeline is composed of a lot of pipes with different slope. Divide each pipe into small blocks with the step length of △x and input the start point data.

2. Establish steady equation (29) on each block section. Input the boundary conditions at the initial point of pipeline (pressure, temperature, gas velocity, and liquid velocity). According to the thermodynamic model, calculate the thermophysic parameters such as density of the gas and liquid phase, gas fraction. Because there is no slip between the two phases at initial point, the liquid hold up can be gained by its relationship with mass flow rate of the gas phase and liquid phase.

3. Set dU/dx as unknowns, and simplify (29) with Gaussian elimination method, then we can obtain more explicit form of (29).

4. Work out pressure (P_i), temperature (T_i), gas and liquid velocity (w_gi and w_li) by four-order Runge - Kutta Method.

5. Figure out liquid holdup (H_Li) by equation (12).

6. Resolve equations (29) by Adams predictor-corrector formula until the reasonable unknowns of this grid section are all gotten.

7. Repeat the second step to the sixth step until reach the last block section which is also the end of this pipeline.

In order to make the numerical calculation converges more quickly, the Adams predictor-corrector and Runge - Kutta Method should be applied simultaneous. As the two methods have four-order accuracy, the desired accuracy also can be improved. The flow chart of the whole solving procedures is depicted in Fig.3.

Figure 3.

Solving procedures of steady state model

6. Transient analysis of condensate gas pipeline

6.1. Basic equations of transient analysis

Opposite to the steady state simulation, the parameters in the general model are dependent on time. Expand (1), (2) and (10), and the following equations can be obtained (Masella et al, 1998).

Gas phase continuity equation:

$$A\phi {(\frac{\partial {\rho }_g}{\partial P})}_T\frac{\partial P}{\partial t}+A{\rho }_g\frac{\partial \phi }{\partial t}+A{w}_g\phi {(\frac{\partial {\rho }_g}{\partial P})}_T\frac{\partial P}{\partial x}+A{\rho }_g\phi \frac{\partial w_g}{\partial x}+A{\rho }_g{w}_g\frac{\partial \phi }{\partial x}=\Delta {\dot{m}}_{gl}$$E31

Liquid phase continuity equation:

$$A{H}_L{(\frac{\partial {\rho }_l}{\partial P})}_T\frac{\partial P}{\partial t}+A{\rho }_l\frac{\partial H_L}{\partial t}+A{w}_l{H}_L{(\frac{\partial {\rho }_l}{\partial P})}_T\frac{\partial P}{\partial x}+A{\rho }_l{H}_L\frac{\partial w_l}{\partial x}+A{\rho }_l{w}_l\frac{\partial H_L}{\partial x}=\Delta {\dot{m}}_{lg}$$E32

Momentum equation

$$\begin{array}{l}{\rho }_g\phi A\frac{\partial w_g}{\partial t}+{\rho }_l{H}_{L}A\frac{\partial w_l}{\partial t}+A\frac{\partial P}{\partial x}+{\rho }_g{w}_g\phi A\frac{\partial w_g}{\partial x}+{\rho }_l{w}_l{H}_{L}A\frac{\partial w_l}{\partial x}\\ =-F_{gw}-F_{lw}-\left({\rho }_g\phi +{\rho }_l{H}_L\right)gA\sin \theta \end{array}$$E33

The transient flow model can be represented by (31) ~ (33) and (12). The unknowns are the pressure P, flow velocity of the gas phase$$w_g$$, flow velocity of the liquid phase $$w_l$$and liquid holdup$$H_L$$. Notice that (31) ~ (33) are a set of partial differential equations so that they can be recast to the following matrix form.

$$B\frac{\partial U}{\partial t}+A\frac{\partial U}{\partial x}=D$$E34

Where, $$A=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]$$,$$B=\left[\begin{array}{ccc}{b}_{11}& b_{12}& b_{13}\\ b_{21}& b_{22}& b_{23}\\ b_{31}& b_{32}& b_{33}\end{array}\right]$$ ,$$D=\left[\begin{array}{c}{D}_1\\ D_2\\ D_3\end{array}\right]$$ ,$$U=\left[\begin{array}{c}P\\ w_g\\ w_l\end{array}\right]$$ $$a_{11}=A{w}_g\phi {(\frac{\partial {\rho }_g}{\partial P})}_T$$,$$a_{12}=A{\rho }_g\varphi$$$$$$,$$a_{13}=0$$ ,$$a_{14}=A{\rho }_g{w}_g$$ ,$$a_{21}=A{w}_l{H}_L{(\frac{\partial {\rho }_l}{\partial P})}_T$$ ,$$a_{22}=0$$ ,$$a_{23}=A{\rho }_l{H}_L$$ ,$$a_{24}=A{\rho }_l{w}_l$$ ,$$a_{31}=A$$ ,$$a_{32}={\rho }_g{w}_g\phi A$$ ,$$a_{33}={\rho }_l{w}_l{H}_{L}A$$ ,.$$a_{34}=0$$.,$$b_{11}=A\varphi {(\frac{\partial {\rho }_g}{\partial P})}_T,$$$$b_{12}=0$$,$$b_{13}=0$$ ,$$b_{14}=A{\rho }_g$$ ,$$b_{21}=A{H}_L{(\frac{\partial {\rho }_l}{\partial P})}_T$$ ,$$b_{22}=0$$ ,$$b_{23}=0$$ ,$$b_{24}=A{\rho }_l$$ ,$$b_{31}=0$$ ,$$b_{32}={\rho }_g\varphi A$$ ,$$b_{33}={\rho }_l{H}_{L}A$$ ,$$b_{34}=0$$$$D_1=\Delta {\dot{m}}_{gl}-b_{14}\frac{\partial \varphi }{\partial t}-a_{14}\frac{\partial \varphi }{\partial x}$$,$$D_2=\Delta {\dot{m}}_{lg}-{\text{b}}_{\text{24}}\frac{\partial \phi }{\partial t}-a_{24}\frac{\partial H_L}{\partial x}$$ ,$$D_3=-F_{gw}-F_{lw}-\left({\rho }_g\varphi +{\rho }_l{H}_L\right)gA\sin \theta$$ The characteristic determinant of (34) is as follow:

$$\left|\begin{array}{ccc}{\lambda }_0{w}_g\phi {\left(\frac{\partial {\rho }_g}{\partial P}\right)}_T+{\lambda }_1\phi {\left(\frac{\partial {\rho }_g}{\partial P}\right)}_T& {\lambda }_0& 0\\ {\lambda }_0{w}_l{H}_L{\left(\frac{\partial {\rho }_l}{\partial P}\right)}_T+{\lambda }_1{H}_L{\left(\frac{\partial {\rho }_l}{\partial P}\right)}_T& 0& {\lambda }_0\\ {\lambda }_0& {\lambda }_0{w}_g+{\lambda }_1& {\lambda }_0{w}_l+{\lambda }_1\end{array}\right|=0$$E35

Where, $${\lambda }_0$$and $${\lambda }_1$$ are the eigenvalues of A and B. If $$\lambda ={\lambda }_0/\lambda$$ established, one solution of $$\lambda$$is 0 and the other two non-zero solutions should be expressed as follow:

$$\begin{array}{l}\left[\phi {\left(\frac{\partial {\rho }_g}{\partial P}\right)}_T+H_L{\left(\frac{\partial {\rho }_L}{\partial P}\right)}_T\right]{\lambda }^2+2\left[\phi w_g{\left(\frac{\partial {\rho }_g}{\partial P}\right)}_T+H_L{w}_l{\left(\frac{\partial {\rho }_l}{\partial P}\right)}_T\right]\lambda \\ +w_l^2{H}_L{\left(\frac{\partial {\rho }_l}{\partial P}\right)}_T+w_g^2\phi {\left(\frac{\partial {\rho }_g}{\partial P}\right)}_T-1=0\end{array}$$E36

The condition for having real solutions of (36) is:

$$\left|w_g-w_l\right|<\sqrt{\frac{1}{H_L}{\left(\frac{\partial P}{\partial {\rho }_l}\right)}_T+\frac{1}{\phi }{\left(\frac{\partial P}{\partial {\rho }_g}\right)}_T}$$E37

In (37), $${\left(\frac{\partial P}{\partial {\rho }_l}\right)}_T$$and $${\left(\frac{\partial P}{\partial {\rho }_g}\right)}_T$$represent square of liquid and gas isothermal wave velocities separately. Thus the right side of (37) is greater than the wave velocity, which means (37) can be established for raw gas pipelines.

According to the analysis above, it can be concluded that under the giving range of operation conditions, all the eigenvaluse of (35) are real different number. Therefore, (34) is strictly hyperbolic and this type will not change as its well posedness.

6.2. Boundary conditions

The boundary conditions include the input and output conditions as well as conditions, such as valves shut off, compressors shut down and etc., which will lead to operating changes. For this transient model, the boundary conditions are considered as: the starting pressure and flow rate over time, ending pressure and flow rate over time.

$${P|}_{x=0}=f_{p1}\left(t\right)$$E38

$${Q|}_{x=0}=f_{q1}\left(t\right)$$E39

$${P|}_{x=L}=f_{p2}\left(t\right)$$E40

$${Q|}_{x=L}=f_{q2}\left(t\right)$$E41

6.3. Solutions of the transient model

Equation (34) is a set of non-liner partial differential equations and the method of characteristic (MOC) and implicit difference method are always adopted to gain numerical solutions. By both methods have advantages and disadvantages. There is strict restriction with the time length of each step (Fig.4), thus the stability of the implicit difference method is better. But at each time layer, differential equations of all the grids should be solved simultaneously, which causes the computing time is large. Different with the implicit difference method, for MOC, the time length of each step is restricted in a relatively short range by its stability conditions. The advantage of MOC is the unknowns of each grid at each time layer can be solved dependently, and the amount of equations solved simultaneously is reduced. Considering the transient feature of low-liquid loading multiphase flow in pipeline is slow transient behaviour which means the transient process lasts long time, the former method is adopted.

Mesh the length variable x and time variable t into grids as shown in Fig. 4. In each block, the differential equations of (31) ~ (33) can be obtained by the implicit different method (Li et al, 2011).

Figure 4.

Differential blocks of pipelines

Gas phase continuity equation:

$$\begin{array}{l}{b}_{11}^n\frac{P_i^{n+1}+P_{i-1}^{n+1}-P_i^n-P_{i-1}^n}{2\Delta t}-b_{14}^n\frac{H_{Li}^{n+1}+H_{Li-1}^{n+1}-H_{Li}^n-H_{Li-1}^n}{2\Delta t}+\\ +a_{11}^n\frac{\xi (P_i^{n+1}-P_{i-1}^{n+1})+(1-\xi )(P_i^n-P_{i-1}^n)}{\Delta x}\\ +a_{12}^n\frac{\xi (w_{gi}^{n+1}-w_{gi-1}^{n+1})+(1-\xi )(w_{gi}^n-w_{gi-1}^n)}{\Delta x}-\\ -a_{14}^n\frac{\xi (H_{Li}^{n+1}-H_{Li-1}^{n+1})+(1-\xi )(H_{Li}^n-H_{Li-1}^n)}{\Delta x}=\Delta {\dot{m}}_{gl}\end{array}$$E42

In which $$\xi$$is the weight coefficient. Set$$\lambda =\frac{2\Delta t}{\Delta x}$$, then above equation reduces to:

$$\begin{array}{l}(b_{11}^n-a_{11}^n\lambda \xi )P_{i-1}^{n+1}-a_{12}^n\lambda \xi w_{gi-1}^{n+1}-(b_{14}^n-a_{14}^n\lambda \xi )H_{Li-1}^{n+1}+(b_{11}^n+a_{11}^n\lambda \xi )P_i^{n+1}+a_{12}^n\lambda \xi w_{gi}^{n+1}\\ -(b_{14}^n+a_{14}^n\lambda \xi )H_{Li}^{n+1}=2\Delta t\Delta {\dot{m}}_{gl}+b_{11}^n(P_i^n+P_{i-1}^n)-b_{14}^n(H_{Li}^n+H_{Li-1}^n)-a_{11}^n\lambda (1-\xi )(P_i^n-P_{i-1}^n)\\ -a_{12}^n\lambda (1-\xi )(w_{gi}^n-w_{gi-1}^n)+a_{14}^n\lambda (1-\xi )(H_{Li}^n-H_{Li-1}^n)\end{array}$$E43

Similarly, liquid phase continuity equation is as follow:

$$\begin{array}{l}(b_{21}^n-a_{21}^n\lambda \xi )P_{i-1}^{n+1}-a_{23}^n\lambda \xi w_{li-1}^{n+1}+(b_{24}^n-a_{24}^n\lambda \xi )H_{Li-1}^{n+1}+(b_{21}^n+a_{21}^n\lambda \xi )P_i^{n+1}+a_{23}^n\lambda \xi w_{li}^{n+1}\\ +(b_{24}^n+a_{24}^n\lambda \xi )H_{Li}^{n+1}=2\Delta t\Delta {\dot{m}}_{\lg }+b_{21}^n(P_i^n+P_{i-1}^n)+b_{24}^n(H_{Li}^n+H_{Li-1}^n)\\ -a_{21}^n\lambda (1-\xi )(P_i^n-P_{i-1}^n)-a_{23}^n\lambda (1-\xi )({\upsilon }_{li}^n-{\upsilon }_{li-1}^n)-a_{24}^n\lambda (1-\xi )(H_{Li}^n-H_{Li-1}^n)\end{array}$$E44

Momentum equation:

$$\begin{array}{l}(b_{32}^n-a_{32}^n\lambda \xi )w_{gi-1}^{n+1}+(b_{33}^n-a_{33}^n\lambda \xi )w_{li-1}^{n+1}-a_{31}^n\lambda \xi P_{i-1}^{n+1}+(b_{32}^n+a_{32}^n\lambda \xi )w_{gi}^{n+1}+(b_{33}^n+a_{33}^n\lambda \xi )w_{li}^{n+1}\\ +a_{31}^n\lambda \xi P_i^{n+1}=2\Delta t{D}_3^n+b_{32}^n(w_{gi}^n+w_{gi-1}^n)+b_{33}^n(w_{li}^n+w_{li-1}^n)-\lambda a_{32}^n(1-\xi )(w_{gi}^n-w_{gi-1}^n)\\ -a_{33}^n\lambda (1-\xi )(w_{li}^n-w_{li}^n)-\lambda a_{31}^n(1-\xi )(P_i^n-P_{i-1}^n)\end{array}$$E45

Based on the above three differential equations, the pressure, flow velocity of the gas phase as well as that of the liquid phase can be obtained. The solution procedures are as follow:

1. Divide pipeline system into series of pipeline segments with different angles. Mesh each segment into blocks with the step length of △x;

2. Get the operation parameters on each grid node by steady state solutions;

3. Put the parameters on 0 time layer into (43) ~ (45) and get the solutions of operation parameters (the pressure P, gas flow rate w_g and liquid flow rate w_l) at 1 time layer;

4. Get the solution of liquid holdup$$H_L$$at 1 time layer by (12);

5. Repeat step 3 ~ step 4 until reach the giving time layer;

The flow chart of the solution procedures is shown in Fig. 5.

7. Application

7.1. Steady state simulation

Table 1 is the basic data of a condensate gas pipeline. The pressure at starting point is 11.0MPa and the temperature is 330K. The mass of gas flow rate is 29.0kg/s and liquid flow is 0kg/s. The components of the condensate gas are shown in Table 2. Calculate the operation parameters of the pipeline by steady state model and determine the location where phase change occurs.

Segments	Length; (km)	Diameter; (mm)	Absolute roughness; (mm)	Total diathermanous factor (W/m2.K)	Ambient; temperature; (K)	Step length; (m)
30	40.0	355.6	0.0457	6.05	330	1333.33

Table 1.

Basic data of the condensate gas pipeline

Figure 5.

Solving procedures of transient model

Component	Fraction	Component	Fraction
C1	0.841454	C6	0.002818
C2	0.075518	C7	0.001535
C3	0.039954	C8	0.001442
C4	0.009476	CO2	0.012166
C5	0.007121	N2	0.008517

Table 2.

Basic data of the condensate gas pipeline

The solutions can be illustrated as Fig. 6. According to the solutions above, the liquid begins to condensate from the section of 12.3km because the hydrocarbons enter the anti-condensate region.

Figure 6.

The steady state pressure variations of the condensate gas pipeline

Figure 7.

The steady state temperature variations of the condensate gas pipeline

Figure 8.

The steady state gas velocity variations of the condensate gas pipeline

Figure 9.

The steady state liquid velocity variations of the condensate gas pipeline

Figure 10.

The steady state liquid holdup variations of the condensate gas pipeline

The feature of condensate gas pipelines is phase change may occur during operating. This leads to a lot of new phenomena as follow:

1. It can be seen from Fig.6 that the pressure drop curve of two phase flow is significantly different from of gas flow even the liquid holdup is quite low. The pressure drop of gas flow is non-linear while the appearance of liquid causes a nearly linear curve of the pressure drop. This phenomenon is expressed that the relatively low pressure in the pipeline tends to increase of the gas volume flow; the appearance of condensate liquid and the temperature drop reduce the gas volume flow.

2. It can be seen from Fig. 7 that the temperature drop curve of two phase flow is similar to single phase flow. The temperature drop gradient of the first half is greater than the last half because of larger temperature difference between the fluid and ambient.

3. It can be seen from Fig. 8 and Fig.9 that the appearance of two phase flow lead to a reduction of gas flow velocity as well as an increase of liquid flow velocity. The phenomenon also contributes to the nearly linear drop of pressure along the pipeline.

4. The sharp change of liquid flow velocity as shown in Fig. 9 is caused by phase change. The initial flow velocity of liquid is obtained by flash calculation which makes no consideration of drag force between the phases. Therefore, an abrupt change of the flow rate before and after the phase change occurs as the error made by the flash calculation cannot be ignored. The two-fluid model which has fully considerate of the effect of time is adopted to solve the flow velocity after phase change and the solutions are closer to realistic. It is still a difficulty to improve the accuracy of the initial liquid flow rate at present. The multiple boundaries method is adapted to solve the steady state model. But the astringency and steady state need more improve while this method is applied to non-linear equations.

5. As shown in Fig.10, the liquid hold up increases behind the phase transition point (two-phase region). Due to the increasing of the liquid hold up is mainly constraint by the phase envelope of the fluid, increasing amount is limited.

The steady state model can simulate the variation of parameters at steady state operation. Actually, there is not absolute steady state condition of the pipeline. If more details of the parameters should be analyzed, following transient simulation method is adopted.

7.2. Transient simulation

Take the previous pipeline as an example, and take the steady state steady parameters as the initial condition of the transient simulation. The boundary condition is set as the pressure at the inlet of pipeline drops to 10.5MPa abruptly at the time of 300s after steady state. The simulation results are shown in Fig11-Fig.15.

Figure 11.

Pressure variation along the pipeline

Figure 12.

Temperature variation along the pipeline

Figure 13.

Velocity of the gas phase variation along the pipeline

Figure 14.

Velocity of the liquid phase variation along the pipeline

Figure 15.

Liquid hold up variation along the pipeline

Compared with steady state, the following features present.

1. Fig.11 depicts the pressure along the pipeline drops continuously with time elapsing after the inlet pressure drops to 10.5MPa at the time of 300s as the changing of boundary condition.

2. Fig.12 shows the temperature variation tendency is nearly the same as steady state. The phenomenon can be explained by the reason that the energy equation is ignored in order to simplify the transient model. The approximate method is reasonable because the temperature responses slower than the other parameters.

3. As depicted in Fig.13, there are abrupt changes of the gas phase velocity at the time of 300s. The opposite direction flow occurs because the pressure at the inlet is lower than the other sections in the pipeline. However, with the rebuilding of the new steady state, the velocity tends to reach a new steady state.

4. Fig.14 shows the velocity variation along the pipeline. Due to the loss of pressure energy at the inlet, the liquid velocity also drops simultaneously at the time of 300s. Similar to gas velocity, after 300s, the liquid velocity increases gradually and tends to reach new steady state with time elapsing.

5. Due to the same liquid hold up equation is adopted in the steady state and transient model, the liquid hold up simulated by the transient model and steady state mode has almost the same tendency (Fig.15). However, the liquid hold up increases because of the temperature along the pipeline after 300s is lower than that of initial condition.

Sum up, the more details of the results and transient process can be simulated by transient model. There are still some deficiencies in the model, which should be improved in further work.

8. Conclusions

In this work, a general model for condensate gas pipeline simulation is built on the basis of BWRS EOS, continuity equation, momentum equation, energy equation of the gas and liquid phase. The stratified flow pattern and corresponding constitutive equation are adopted to simplify the model.

By ignoring the parameters variation with time, the steady state simulation model is obtained. To solve the model, the four-order Runge - Kutta method and Gaussian elimination method are used simultaneously. Opposite to steady state model, the transient model is built with consideration of the parameters variation with time, and the model is solved by finite difference method. Solving procedures of steady-state and transient models are presented in detail.

Finally, this work simulated the steady-state and transient operation of a condensate gas pipeline. The pressures, temperatures, velocity of the gas and liquid phase, liquid hold up are calculated. The differences between the steady-state and transient state are discussed. The results show the model and solving method proposed in this work are feasible to simulate the steady state and transient flow in condensate gas pipeline. Nevertheless, in order to expand the adaptive range the models, more improvements should be implemented in future work (Pecenko et al, 2011).

Acknowledgments

This paper is a project supported by sub-project of National science and technology major project of China (No.2008ZX05054) and China National Petroleum Corporation (CNPC) tackling key subject: Research and Application of Ground Key Technical for CO₂ flooding, JW10-W18-J2-11-20.