A Hydraulic Model for Multiphase Flow Based on the Drift Flux Model in Managed Pressure Drilling

Abstract

Managed pressure drilling (MPD) is a drilling technique used to address the narrow density window under complex geological environments. It has widespread applications in the exploration and exploitation of oil and gas, both onshore and offshore. In this study, to achieve effective control of the downhole pressure to ensure safety, a gas–liquid two-phase flow model based on the drift flux model is developed to describe the characteristics of transient multiphase flow in the wellbore. The advection upwind splitting method (AUSM) numerical scheme is used to assist with calculation and analysis, and the monotonic upwind scheme for conservation laws (MUSCLs) technique with second-order precision is adopted in combination with the Van Leer slope limiter to improve precision. Relevant data sourced from prior literature are used to validate the suggested model, the results of which reveal an excellent statistical consistency. Further, the influences of various parameters in a field application, including backpressure, density, and mass flow, are analyzed. Over the course of later-stage drilling, a combination of wellhead backpressure and displacement is recommended to exercise control.

1. Introduction

In the oil and gas industry, the continuously increasing demand for underground resources has led to a shift in the focus on exploration and exploitation to deep or ultra-deep reservoirs. Due to the complicated geological conditions of these reservoirs, drilling technology aimed at creating a wellbore to reach the gas or oil reservoir encounters a multitude of challenges, particularly within narrow pore-fracture pressure windows [1,2]. Over the course of drilling, incidents occur frequently and cause various disadvantages, such as longer drilling times, higher costs, and vast risks. Conventional drilling is incapable of addressing these issues, and advanced technology-managed pressure drilling (MPD) is proposed as an alternative, which is useful for resolving complicated situations such as leakage, gas-kick, collapse, and trapping within a narrow density window through reducing the non-productive time and avoiding incidents [3,4]. Regarding offshore drilling, MPD shows remarkable advantages in coping with complex situations and has been widely used in drilling engineering, considering the continued exploration of challenging reservoirs [5,6]. Moreover, for the drilling and production of natural gas hydrate, MPD also creates incomparable advantages, including precise monitoring and control of pressure, avoidable formation damage and unwanted fractures, and rare gas kick and lost circulation [7,8].

Managing annular pressure in wellbores is considered a challenging task, as flow materials are likely to contain mud, oil, gas, and cuttings [9]. Typically, a gas–liquid two-phase flow is included in MPD with the homogeneous flow model (HFM), two-fluid model (TFM), and drift flow model (DFM). The HFM is incapable of ensuring accuracy in reflecting the flow characteristics of each phase, and its calculation accuracy is less than satisfactory; thus, it is rarely applied nowadays. A four-dimensional hypercube has been suggested to describe gas–liquid two-phase flow, and a set of hyperbolic partial differential equations (PDEs) can be derived [10,11]. Of these, the most complicated is the high-fidelity model, which has the capability of describing how physical variables (momentum, mass, and energy) are transmitted and meet the requirements of the thermodynamic equilibrium equation [12,13,14,15]. Despite this, these models are highly complicated to solve, and it is possible that the calculation process does not converge enough to use high-order numerical schemes. Therefore, the simplified DFM was first suggested by Zuber and Findlay [16] and was subsequently improved by Ishii [17]. In the DFM, mass conservation of the gas–liquid two-phase is taken into consideration separately, and the interactions occurring between the two phases are discounted. Relative to the TFM, the DFM demonstrates various advantages such as having a simple form and less computation; thus, it has been commonly applied in the theoretical research process [18,19,20]. In order to enhance its accuracy, the majority of the following studies concentrated on the empirical parameters of this constitutive relationship. In fact, the calculation process remains complicated to some extent. In addition, there are many more simplified equations that have been derived for application under different scenarios, including early gas-kick detection and parameter estimation [21,22,23,24]. However, these models give no consideration to how the characteristics of various parameters, which are primarily used for automation and are not suited to fine MPD, are distributed in the wellbore. Figure 1 shows the relation between the models [11].

Concerning the DFM, accuracy and efficiency are crucial for the mathematical algorithm. The DFM is a typical hyperbolic equation that can be solved using the finite volume method. The Godunov and Roe schemes have been suggested to solve the DFM, but they depend on algebraic Riemann solvers or analytic Jacobian matrices of conservation equations [25,26]. A new pressure-based method was developed to address the convergence problem of the drift flux model by Wang [27]. Among the different schemes, it has been confirmed that the advection upwind splitting method (AUSM) schemes are effective solutions to the DFM in terms of their simplicity, accuracy, and robustness [28,29,30]. Not only do these schemes display accuracy in flux difference splitting (FDS) in the boundary layer, they also exhibit robustness in flux vector splitting (FVS) when capturing strong discontinuity. AUSM schemes are relatively late in performing the calculation of gas–liquid two-phase flow, and nearly all of them were studied after the year 2000. Niu [31] applied the AUSM based on flux difference (AUSMD) scheme to solve the flow model of seven equations. Nevertheless, it was discovered that Niu’s method was susceptible to numerical oscillation. Subsequently, Paillère et al. [32] achieved success in applying AUSM+ to solve the TFM of six equations. Evje et al. [28,33] also used the AUSM+ and AUSM based on flux vector (AUSMV) methods on a flow model of four equations and a drift flux model of three equations. Moreover, this method is applicable to engineering computing in drilling [34,35].

According to a general drift flux model of compressed gas–liquid two-phase flow, a model to describe the transient flow dynamics in the wellbore for the MPD is established in this paper. Pressure variations of the throttle valve at the wellhead and flow change resulting from the pressure difference at the bottom hole are considered as well. The governing equation of the DFM is numerically solved by the AUSM scheme with a second-order spatial accuracy and Van Leer slope limiter. The data obtained from published papers are used to validate the established model. In addition, a field application case is researched to support data analysis to make operation decisions.

2. Mathematical Model

2.1. Governing Equations

The typical hydraulic cycle used in the process of MPD is shown in Figure 2. The mud from the rig pump flows into the drill pipe through the surface equipment and riser, and flows through the bit nozzle at the bottom of the drill pipe. Then it enters the annular space and flows upward. After reaching the surface, it is processed by the solids control equipment and returns to the mud pit, then enters the rig pump again. The hydraulic model is the basis of the design and calculation of pressure control. MPD requires a real-time estimate based on the measured data, and it inputs the calculated results into the control system to achieve automatic control. It is costly to use pressure while drilling (PWD) tools to measure the downhole pressure, so the hydraulic model is of great significance. Additionally, even in the case of PWD, the hydraulic model is an essential auxiliary analytical tool.

As discussed in the previous section, the DFM is primarily a simplification of high-fidelity models using a complex conservation equation. In DFM, a mixture momentum equation is adopted rather than using an isolated equation for each phase. To achieve a high level of precision, some closure relations are applied into the equations. The basic DFM includes two conservation equations for the gas and liquid and a mixture momentum equation [28]:

$$\frac{\partial }{\partial t}\left(\begin{array}{c}{\rho }_l{\alpha }_l\\ {\rho }_g{\alpha }_g\\ {\rho }_g{\alpha }_g{v}_g+{\rho }_l{\alpha }_l{v}_l\end{array}\right)+\frac{\partial }{\partial x}\left(\begin{array}{c}{\rho }_l{\alpha }_l{v}_l\\ {\rho }_g{\alpha }_g{v}_g\\ {\rho }_g{\alpha }_g{v}_g^2+{\rho }_l{\alpha }_l{v}_l^2+p)\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ -F_G-F_f\end{array}\right)$$

(1)

Equation (1) can be written in vector form:

$$\frac{\partial \mathbf{u}}{\partial t}+\frac{\partial \mathbf{F}(\mathbf{u})}{\partial x}=\mathbf{Q}(\mathbf{u})$$

(2)

where:

$$\mathbf{u}=\left(\begin{array}{c}{\rho }_l{\alpha }_l\\ {\rho }_g{\alpha }_g\\ {\rho }_g{\alpha }_g{v}_g+{\rho }_l{\alpha }_l{v}_l\end{array}\right),\mathbf{F}(\mathbf{u})=\left(\begin{array}{c}{\rho }_l{\alpha }_l{v}_l\\ {\rho }_g{\alpha }_g{v}_g\\ {\rho }_g{\alpha }_g{v}_g^2+{\rho }_l{\alpha }_l{v}_l^2+p\end{array}\right),\mathbf{Q}(\mathbf{u})=\left(\begin{array}{c}0\\ 0\\ -F_G-F_f\end{array}\right)$$

(3)

where ρ, p, and υ are the density, pressure, and velocity. α is the volume fraction. The subscripts l and g denote the liquid and gas phases, respectively.

F_G represents the gravity term:

$$F_G=g\left({\alpha }_{\mathrm{l}}{\rho }_{\mathrm{l}}+{\alpha }_{\mathrm{g}}{\rho }_{\mathrm{g}}\right)\cos \theta$$

(4)

F_f is the friction term:

$$F_f=\frac{2f_m{\rho }_m{v}_m\left|v_m\right|}{D-d}$$

(5)

where:

$${\rho }_m={\alpha }_l{\rho }_l+{\alpha }_g{\rho }_g$$

(6)

$$v_m={\alpha }_l{v}_l+{\alpha }_g{v}_g$$

(7)

f_m can be calculated by [26]:

$$f_m=\{\begin{array}{ll}\frac{16}{{\mathrm{Re}}_m}\hfill & \mathrm{if}{\mathrm{Re}}_m\le 2400\hfill \\ \frac{1}{{\left[4\log \left(\frac{\Delta /(D-d)}{3.7065}-\frac{5.0452}{{\mathrm{Re}}_m}\log \Lambda \right)\right]}^2}\hfill & {\mathrm{if}\mathrm{Re}}_m\ge 2400\hfill \end{array}$$

(8)

$$\Lambda =\frac{{[\Delta /(D-d)]}^{1.1098}}{2.8257}+{\left(\frac{7.149}{{\mathrm{Re}}_m}\right)}^{0.8981}$$

(9)

$${\mathrm{Re}}_m=\frac{{\rho }_m{v}_m(D-d)}{{\mu }_m}$$

(10)

where D and d denote the outer diameter and inner diameter of the annulus, respectively. Re is the Reynolds number; and μ_m is the viscosity.

2.2. Closure Equations

For a constant mass flow, the following equations must be satisfied [20]:

$$\rho a=(\rho +d\rho )(a+da)$$

(11)

$$\rho a=\rho a+ad\rho +\rho da+d\rho da$$

(12)

a is the speed of sound in the fluid.

Neglect the last term:

$$a=-\rho \frac{da}{d\rho }$$

(13)

Similarly, for the momentum:

$$p+\rho a^2=(p+dp)+(\rho +d\rho ){(a+da)}^2$$

(14)

$$da=-\frac{dp+a^{2}d\rho }{2a\rho }$$

(15)

$$a=\rho \frac{dp+a^{2}d\rho }{2a\rho d\rho }$$

(16)

$$a=\frac{dp}{2ad\rho }+\frac{a}{2}$$

(17)

$$a^2=\left(\frac{dp}{d\rho }\right)$$

(18)

Therefore, the density of the liquid is described as:

$${\rho }_l={\rho }_{l,0}+\frac{p-p_{l,0}}{a_l^2}$$

(19)

Similarly,

$${\rho }_g=\frac{p}{a_g^2}$$

(20)

The following relationship is satisfied:

$${\alpha }_g+{\alpha }_l=1$$

(21)

In addition, the slip velocity equation proposed by Zuber and Findlay [16] is represented by

$$v_g{=C}_0{v}_m{+v}_t$$

(22)

where C₀ is the distribution parameter and υ_t is the gas drift velocity.

2.3. Shi Slip Relation

In this research, a Shi slip relation was adopted, which can be found in [36]:

$$v_t=\frac{\left(1-C_0{\alpha }_g\right)V_{c}K\left({\alpha }_g,K_u,C_0\right)m(\varphi )}{C_0{\alpha }_g\sqrt{\frac{{\rho }_g}{\rho }}+1-C_0{\alpha }_g}$$

(23)

where:

$$K=\{\begin{array}{ll}1.53/C_0\hfill & \mathrm{if}\hspace{1em}\hspace{1em}{\alpha }_g\le a_1\hfill \\ \frac{1.53}{C_0}+\frac{C_0{K}_u-1.53}{2C_0}\left[1-\cos \left(\pi \frac{{\alpha }_g-a_1}{a_2-a_1}\right)\right]\hfill & \mathrm{if}\hspace{1em}\hspace{1em}{a}_1<{\alpha }_g\le a_2\hfill \\ K_u\hfill & \mathrm{if}\hspace{1em}\hspace{1em}{\alpha }_g>a_2\hfill \end{array}$$

(24)

with a₁ = 0.2 and a₂ = 0.4.

K_u is the critical Kutateladze number:

$$K_u={\left[\frac{C_{K_u}}{\sqrt{N_B}}\left(\sqrt{1+\frac{N_B}{C_{K_u}^2{C}_w}}-1\right)\right]}^{1/2}$$

(25)

where $C_{K_u}$ and C_w denote the fitted parameters, which are set to 142 and 0.008, respectively.

N_B is the Bond number:

$$N_B=\left[\frac{g\left({\rho }_l-{\rho }_g\right)}{{\sigma }_{gl}}\right]{(D-d)}^2$$

(26)

where σ_gl is the gas–liquid interfacial tension.

V_c is the characteristic velocity:

$$V_c={\left(\frac{{\sigma }_{gl}g\left({\rho }_l-{\rho }_g\right)}{{\rho }_l^2}\right)}^{1/4}$$

(27)

The distribution parameter C₀ can be determined by:

$$C_0=\frac{C_{\max }}{1+\left(C_{\max }-1\right){\gamma }^2}$$

(28)

C_max represents the value of the gas distribution parameter, which is usually set to 1.2.

γ is calculated by:

$$\gamma =\frac{\beta -B}{1-B}$$

(29)

where:

$$\beta =\max \left({\alpha }_g,F_v\frac{{\alpha }_g\left|v_m\right|}{v_{sgf}}\right)$$

(30)

where F_υ is a multiplier of the flooding velocity fraction that defaults to 1.0.

$$v_{sgf}=V_c{K}_u{\left(\frac{{\rho }_l}{{\rho }_g}\right)}^{1/2}$$

(31)

The above relation is based on the derivation of the vertical flow. In an inclined pipe, the terminal rise velocity can be indicated as follows:

$$V_{\infty \theta }=V_{\infty }m\left(\theta \right)$$

(32)

where θ is the angle of deviation and m is the drift velocity multiplier.

m(θ) represents the influence of the wellborn inclination:

$$m\left(\theta \right)=m_0{\left(\cos \theta \right)}^{n_1}{\left(1+\sin \theta \right)}^{n_2}$$

(33)

For the water–gas flow, m₀ = 1.28, n₁ = 0.5, n₂ = 1.7 (θ < 88°). When 88°≤ θ ≤ 90°, the relation has a large deviation for horizontal wells.

2.4. Primitive Variables

Equation (1) can be expressed as:

$${\partial }_t\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\end{array}\right]+{\partial }_x\left[\begin{array}{l}{v}_l{u}_1\hfill \\ v_g{u}_2\hfill \\ v_l^2{u}_1+v_g^2{u}_2+p\left(u_1,u_2\right)\hfill \end{array}\right]=\left[\begin{array}{c}0\\ 0\\ -F_G-F_f\end{array}\right]$$

(34)

where $u_1={\alpha }_l{\rho }_l$, $u_2={\alpha }_g{\rho }_g$ and $u_3={\alpha }_l{\rho }_l{v}_l+{\alpha }_g{\rho }_g{v}_g$.

Equation (21) can be expressed in terms of u₁ and u₂:

$$\frac{u_1}{{\rho }_l}+\frac{u_2}{{\rho }_g}=1$$

(35)

Substitute Equations (19) and (20) into Equation (35):

$$\frac{u_1}{{\rho }_{l,0}+\frac{p-p_{l,0}}{a_l^2}}+\frac{u_2{a}_g^2}{p}=1$$

(36)

It can be obtained by simple transformation:

$$p^2+Bp+C=0$$

(37)

where:

$$B=a_l^2\left({\rho }_{l,0}-\frac{p_{l,0}}{a_l^2}-u_1-{\left(\frac{a_g}{a_l}\right)}^2{u}_2\right)$$

(38)

$$C=u_2{\left(a_g{a}_l\right)}^2\left(\frac{p_{l,0}}{a_l^2}-{\rho }_{l,0}\right)$$

(39)

The pressure is then given by:

$$p\left(u_1,u_2\right)=\frac{-B+\sqrt{B^2-4C}}{2}$$

(40)

The volume fractions are determined as:

$${\alpha }_l=\frac{u_1}{{\rho }_l}$$

(41)

$${\alpha }_g=\frac{u_2}{{\rho }_g}$$

(42)

One can get the following equation by Equation (34):

$$u_1{v}_l{+u}_2{v}_g{=u}_3$$

(43)

According to Equation (22):

$$-C_o{\alpha }_l{v}_l+\left(1-C_o{\alpha }_g\right)v_g=v_t$$

(44)

A system of linear equations from Equations (43) and (44) can be expressed in matrix form:

$$\left(\begin{array}{cc}{u}_1& u_2\\ w_1& w_2\end{array}\right)\left(\begin{array}{c}{v}_l\\ v_g\end{array}\right)=\left(\begin{array}{c}{u}_3\\ v_t\end{array}\right)$$

(45)

Hence, we can get expressions for υ_l and υ_g:

$$v_l=\frac{u_3{w}_2-u_2{v}_d}{u_1{w}_2-u_2{w}_1}\Rightarrow \frac{u_3\left(1-C_o{\alpha }_g\right)-u_2{v}_d}{u_1\left(1-C_o{\alpha }_g\right)+u_2{C}_o{\alpha }_l}$$

(46)

$$v_g=\frac{u_1{v}_d-u_3{w}_1}{u_1{w}_2-u_2{w}_1}\Rightarrow \frac{u_1{v}_d+u_3{C}_o{\alpha }_l}{u_1\left(1-C_o{\alpha }_g\right)+u_2{C}_o{\alpha }_l}$$

(47)

2.5. Boundary Conditions

At the inlet, the convective flux is determined by the flow of the gas and liquid. The pressure flux can be obtained by extrapolation:

$$P_{inlet}=p(1)-0.5(p(1)-p(2))$$

(48)

It is necessary to analyze the two operating conditions of a shut-in and open well at the outlet boundary. During shut-in, the flow rates of gas and liquids at the wellhead are zero, and then the convective flux is zero. The pressure flux can be calculated as:

$$P_{outlet}=p(N)+0.5(p(N)-p(N-1))$$

(49)

While the well is open, the convective flux can be dealt with using the original variables at the outlet (velocity, density, volume fraction). The pressure flux is equal to the ground backpressure.

The two-phase flow in the choke valve can be represented as [20]:

$$\frac{{\rho }_l{v}_l{\alpha }_l}{\sqrt{{\rho }_l}}+\frac{{\rho }_g{v}_g{\alpha }_g}{\gamma \sqrt{{\rho }_g}}=K_c{z}_c\sqrt{\frac{2}{{\rho }_a}(P_{choke}-P_{atm}})$$

(50)

where K_c is the valve discharge coefficient, z_c is the choke valve opening index, and γ is the gas expansion factor. P_choke and P_atm represent the pressure of the choke and the atmosphere, respectively.

The inflow rate resulting from the pressure difference between the reservoir and the bottom hole can be expressed as:

$$p_r^2-p_{wf}^2=Aq+B{q}^2$$

(51)

P_r and P_wf are the reservoir pressure and bottom-hole pressure, respectively. q is the inflow rate. A and B are regression coefficients, which are mainly determined by the reservoir.

3. Numerical Scheme

3.1. Advection Upwind Splitting Method Scheme

The flux at the cell interface is composed of the convective flux and pressure flux, which can be represented as [26,31]:

$$\mathbf{F}{(\mathbf{u})}_{j+1/2}={\left({\mathbf{F}}_l^c+{\mathbf{F}}_g^c+{\mathbf{F}}^p\right)}_{j+1/2}$$

(52)

${\mathbf{F}}_l^c$ and ${\mathbf{F}}_g^c$ represent the liquid and gas phase convective flux, respectively, and ${\mathbf{F}}^p$ represents the pressure flux.

For the convective terms:

$${\mathbf{F}}_g^c=\left({\alpha }_g{\rho }_g\right){\Psi }_g^{\pm }={\left({\alpha }_g{\rho }_g\right)}_L{\Psi }_{gL}^{+}+{\left({\alpha }_g{\rho }_g\right)}_R{\Psi }_{gR}^{-}$$

(53)

$${\mathbf{F}}_l^c=\left({\alpha }_l{\rho }_l\right){\Psi }_l^{\pm }={\left({\alpha }_l{\rho }_k\right)}_L{\Psi }_{lL}^{+}+{\left({\alpha }_l{\rho }_l\right)}_R{\Psi }_{lR}^{-}$$

(54)

where Ψ is the speed function, l and g are liquid and gas phases, respectively, and L and R represent left and right cells, respectively.

$${\Psi }_{gL}^{+}={\overline{V}}^{+}{\left(\begin{array}{ccc}0& 1& v\end{array}\right)}^T,{\Psi }_{lL}^{+}={\overline{V}}^{+}{\left(\begin{array}{ccc}1& 0& v\end{array}\right)}^T$$

(55)

$${\Psi }_{gR}^{-}={\overline{V}}^{-}{\left(\begin{array}{ccc}0& 1& v\end{array}\right)}^T,{\Psi }_{lR}^{-}={\overline{V}}^{-}{\left(\begin{array}{ccc}1& 0& v\end{array}\right)}^T$$

(56)

where:

$${\overline{V}}^{\pm }(v,c,\chi )=\{\begin{array}{ll}\chi {\mathrm{V}}^{\pm }(v,c)+(1-\chi )\frac{v\pm \left|v\right|}{2}\hfill & \mathrm{if}\left|v\right|\le c\hfill \\ \frac{v\pm \left|v\right|}{2}\hfill & \mathrm{otherwise}\hfill \end{array}$$

(57)

χ is the weight coefficient and can be expressed as:

$${\chi }_L=\frac{2{(\rho /\alpha )}_L}{{(\rho /\alpha )}_L+{(\rho /\alpha )}_R}$$

(58)

$${\chi }_R=\frac{2{(\rho /\alpha )}_R}{{(\rho /\alpha )}_L+{(\rho /\alpha )}_R}$$

(59)

$$V^{\pm }(v,c)=\{\begin{array}{ll}\pm \frac{1}{4c}{(v\pm c)}^2& \left|v\right|\le c\\ \frac{1}{2}(v\pm \left|v\right|)& \left|v\right|>c\end{array}$$

(60)

where v is the mixture velocity and c is the sound velocity.

$$c=\{\begin{array}{cc}{a}_l\hfill & a_g<\epsilon \\ w\hfill & \epsilon \le a_g\le 1-\epsilon \\ a_g\hfill & a_g>1-\epsilon \end{array}$$

(61)

$$w=\sqrt{\frac{p}{{\alpha }_g{\rho }_l\left(1-C_0{\alpha }_g\right)}}$$

(62)

ε is a small parameter to ensure smooth transition.

The pressure flux is updated by:

$$p_{j+1/2}=P^{+}\left(v_{m,L},c_{j+1/2}\right)p_L+P^{-}\left(v_{m,R},c_{j+1/2}\right)p_R$$

(63)

$$P^{\pm }(v,c)={\overline{V}}^{\pm }(v,c)\{\begin{array}{ll}\frac{1}{c}\left(\pm 2-\frac{v}{c}\right)\hfill & \mathrm{if}\left|v\right|\le c\hfill \\ \frac{1}{v}\hfill & \mathrm{otherwise}\hfill \end{array}$$

(64)

3.2. Second-Order Accuracy

In this paper, we use a classical second-order monotonic upwind scheme for conservation laws (MUSCLs) technique to improve the accuracy:

$$w_{j+1/2,L}=w_j+\frac{1}{2}S\left({\varphi }_j\right)\left(w_{j+1}-w_j\right)$$

(65)

$$w_{j+1/2,R}=w_{j+1}-\frac{1}{2}S\left({\varphi }_{j+1}\right)\left(w_{j+2}-w_{j+1}\right)$$

(66)

where $w=\left(p,{\rho }_g,{\rho }_l,{\alpha }_g,{\alpha }_l,v_g,v_l\right)$.

S is the Van Leer slope limiter:

$$S(\varphi )=\frac{\left|\varphi \right|+\varphi }{1+\left|\varphi \right|}$$

(67)

3.3. Solution Method

When Equation (2) is discretized explicitly:

$${\mathbf{u}}_j^{n+1}={\mathbf{u}}_j^n-\frac{\Delta t}{\Delta x}\left({\mathbf{F}}_{j+1/2}^n-{\mathbf{F}}_{j-1/2}^n\right)+\Delta t{\mathbf{Q}}_j^n$$

(68)

Whereas Equation (2) is discretized implicitly:

$${\mathbf{u}}_j^{n+1}-{\mathbf{u}}_j^n+\frac{\Delta t}{\Delta x}\left({\mathbf{F}}_{j+1/2}^{n+1}-{\mathbf{F}}_{j-1/2}^{n+1}\right)-\Delta t{\mathbf{Q}}_j^{n+1}=0$$

(69)

The Newton–Raphson method is used for solving the implicit equation, so that the process of the calculation is stable.

3.4. Workflow

The primary calculation steps are shown in Figure 3, as follows. (1) Input the known parameters (well depth, bottom hole assembly (BHA), fluid properties, etc.) and complete the initialization according to the calculation requirements. (2) Calculate the raw data with the MUSCLs of second-order accuracy and obtain temporary data. (3) Use the AUSM method to analyze the elements from the bottom up, considering the Shi slip relationship, and then convert the conservative variables to primitive variables. (4) After one time step, store the result of this moment as the initial condition for calculating the next time step. (5) Repeat the above steps until the preset time period is over, and then store all of the data.

4. Experimental Validation

4.1. Laboratory Test

To complete the verification of the established model, we designed a simple experimental device, as shown in Figure 4. The device mainly includes a wellbore simulator, liquid circulation system, and air supply unit, as shown in Figure 5, Figure 6, Figure 7 and Figure 8. The wellbore is a 6 m long transparent glass tube with a diameter of 180 mm, inside of which is an 80 mm blue plastic tube. The fluid is injected into the blue string, and the air is injected at the bottom. The air enters the annulus, and a gas–liquid two-phase flow is formed. Multiple pressure sensors are set in the wellbore to record the changes in pressure.

During the experiment (Figure 9), the annulus was filled with liquid at the beginning with a flow rate of 0.0028 m³/s and then injected with gas at 25 s, of which the flow rate was 0.01 kg/s. Figure 8 shows the annulus flow dynamics. The flow state is mainly divided into three stages: the liquid phase flow stage (0–25 s), the transient gas–liquid two-phase flow stage (25–32 s), and the stable flow stage (32–60 s). Generally, the pressure in the wellbore increases first and then decreases, and finally, tends to be stable.

Figure 10 shows how the bottom hole pressure changes over time, and the simulated data are in good agreement with the experimental results. The error range is within ±10%.

4.2. Full-Scale Experiment

Lage and Fjelde [37] conducted their injection experiments in a real well. The vertical test well was 1275 m deep with a parasitic pipe used for gas injection at a depth of 760 m and a pressure sensor fitted at 605 m. At the start, the wellbore was filled with water at a flow rate of 10.11 m³/s prior to the air being injected into the wellbore. The change in gas flow is presented in Figure 11. Over the course of gas injection, there was a significant change in the gas volume. The total time taken to conduct this experiment amounted to 1500 s, and the injection depth was relatively deep, which was used to exhibit the characteristics of the gas invasion process in practice.

Figure 12 presents a diagram that shows the changes in pressure over time at a latitude of 605 m. In this paper, the original measurement data as well as the calculation results obtained by Meng et al. [35] and Lage and Fjelde [37] are compared. As revealed by this figure, there is a decent consistency between the simulation results and the on-site testing data. The curve obtained in this study is smoother than that obtained by Lage and Fjelde. Meanwhile, the measurement data were contrasted with the calculation results, which indicates that the resulting error is within ±5%.

In order to further understand the flow dynamics in the wellbore, we analyzed the variation of other parameters according to the simulated data. Figure 13 reflects the changes in pressure with time and depth. As can be seen from the figure, in both the measured and simulated data, the pressure shows a gradual decline, which is followed by a sharp fall, before an incremental increase again. The drop in pressure is due to the original fluid within the ambient environment being dispelled with the ingress of gas. It is known that gas has a lower density than fluid, causing the downhole pressure to decline on a continued basis. As the gas moves upwards, the pressure it is subjected to decreases, which causes the gas to expand further. As a result, more fluid is dispelled, and the drop in pressure occurs faster. Up to the point when the gas reaches the wellhead, the downhole pressure is at its minimum. With the continued ingress of gas, there will be an increase in downhole pressure to some extent until the gas–fluid two-phase flow stabilizes. Figure 14 shows the changes in the gas holdup. As revealed by this figure, there is a constant change in gas holdup at varying depths over time, which is consistent with the changes in the gas mass flow. This is most obvious at the wellhead. The gas holdup reaches its maximum at the wellhead when the maximum gas mass flow reaches the wellhead.

Figure 15 presents an analysis of the gas mass flow. As revealed in the figure, there is a significant change in gas mass flow, with the maximum at 0.265 kg/s. The variation in gas mass flow accounts for the changes to all other characterizing parameters for the wellbore. In practice, due to the difference between the geological pressure and that at the downhole, there will be a constant change in the gas mass flow, as a result of which the characteristics exhibited by the transient multiphase flow will change inside the wellbore. Figure 16 reflects the changes in fluid mass flow. The maximum fluid mass at the wellhead is 63 kg/s, and the minimum is 4 kg/s, which is because as the amount of filled gas is on the rise and the fluid inside the wellbore is dispelled, the fluid mass flow is significant when there is less gas ingress. When the gas mass increases, the fluid mass flow will drop noticeably.

5. Sensitivity Analysis

In the previous sections, the model proposed in this paper was validated and shown to be suited to industrial application in the oil and gas industry. One example is a well located in the north of Sichuan Province, where the MPD technique is applied. To complete the calculation, it was supplemented with the following parameters (

Table 1. Geometric parameters.

Well Depth (m)	Outer Tubing Diameter (m)	Casing Diameter (m)
0–500	0.089	0.2523
500–2800	0.089	0.15708
2800–3600	0.127	0.1469

and

Table 2. Calculation parameters.

Parameter	Unit	Value	Parameter	Unit	Value
Liquid density	kg/m3	1100	Surface temperature	K	293.15
Back pressure	Pa	101325	Temperature gradient	K/(100 m)	2.2
Liquid viscosity	Pa·s	0.02	Surface tension	N/m	0.0761
Gas viscosity	Pa·s	0.0000015	Sound velocity in liquid	m/s	1200

). For more effective control of the downhole pressure during drilling, the application of the existing technology to improve the effectiveness in downhole pressure control is discussed, and an analysis is conducted of the impact of different parameters on flow inside the wellbore.

The wellhead backpressure has a direct effect on the downhole pressure. As shown in Figure 17, lower wellhead backpressure leads to a faster decline in the downhole pressure when overflow arises. When the wellhead backpressure is 0.1 MPa, the downhole pressure declines by 7.5 MPa, which means an increment of 1.9 MPa as compared to the time point when the wellhead backpressure is 2.0 MPa. Therefore, it can be seen that applying pressure to some extent at the wellhead is effective for bringing the reduction in the downhole pressure under control.

Even a slight decline in the density of drilling fluid will cause a sharp drop in the downhole pressure. As shown in Figure 18, the downhole pressure difference is about 3.4 MPa when the difference in drilling fluid density is 0.1 kg/m³. A commonly-seen practice for controlling the downhole pressure is to change the density of the drilling fluid. However, as it takes time for the drilling fluid to reach the downhole from the wellhead, the downhole pressure is incapable of being adjusted promptly, which means that the possibility of various dangerous incidents occurring remains possible. Therefore, a combined approach taken while drilling with backpressure and displacement would be more effective and reasonable to prevent the occurrence of incidents.

Over the course of drilling, displacement has a substantial impact on the downhole pressure. As indicated by Figure 19, a lower displacement leads to a lower downhole pressure, which will decline at a faster pace over time. Therefore, as far as the drilling facilities allow, an increase in the displacement will cause a significant rise in the downhole pressure.

In actual conditions, the pressure difference between the reservoir and the downhole will constantly change, and the invasion amount of gas will also vary synchronously. As shown in Figure 20, a higher initial pressure difference causes a more significant variation in the downhole pressure. Transient gas–liquid two-phase flow is always present in the wellbore, and there is no stable stage. The downhole pressure drop can be divided into two stages: the accelerated stage while gas is in the annulus; and the steady descent stage after the gas reaches the wellhead.

6. Conclusions

In this paper, an MPD wellbore transient flow dynamic model was developed based on the compressible gas–fluid two-phase flow with the drift flux model. The most critical findings are summarized below:

(1)

In laboratory experiments, the flow state is mainly divided into three stages: the liquid phase flow stage (0–25 s), the transient gas–liquid two-phase flow stage (25–32 s), and the stable flow stage (32–60 s). The simulated data are in good agreement with the experimental results, and the error range is within ±10%.

(2)

The pressure shows a gradual decline, which is followed by a sharp fall, before an incremental increase again. The drop in pressure is due to the original fluid being dispelled with the ingress of gas. As the gas moves upwards, the pressure it is subjected to decreases, which causes the gas to expand further. As a result, more fluid is dispelled, and the drop in pressure occurs faster. Up to the point when the gas reaches the wellhead, the downhole pressure is at its minimum.

(3)

The adjustment of wellhead back pressure is mainly realized by throttle valve. The higher the wellhead back pressure is, the smaller the downhole pressure will be. When the gas–liquid two-phase flow in the wellbore reaches an equilibrium state, the downhole pressure will decrease less with the increase of drilling fluid displacement, and the time of gas reaching the wellhead will be earlier.

(4)

The downhole pressure can be controlled by changing the density of drilling fluid. However, the adjustment of drilling fluid density has a serious lag. Considering the variation of gas invasion caused by reservoir pressure difference, there will be no stable gas–liquid two-phase flow equilibrium.

Further work needs to be done to couple the proposed method with deformation near the well at a small scale.