# Calculation Analysis of Pressure Wave Velocity in Gas and Drilling Mud Two-Phase Fluid in Annulus during Drilling Operations

Authors: Yuanhua Lin; Xiangwei Kong; Yijie Qiu; Qiji Yuan
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Investigation of propagation characteristics of a pressure wave is of great significance to the solution of the transient pressure problem caused by unsteady operations during management pressure drilling operations. With consideration of the important factors such as virtual mass force, drag force, angular frequency, gas influx rate, pressure, temperature, and well depth, a united wave velocity model has been proposed based on pressure gradient equations in drilling operations, gas-liquid two-fluid model, the gas-drilling mud equations of state, and small perturbation theory. Solved by adopting the Runge-Kutta method, calculation results indicate that the wave velocity and void fraction have different values with respect to well depth. In the annulus, the drop of pressure causes an increase in void fraction along the flow direction. The void fraction increases first slightly and then sharply; correspondingly the wave velocity first gradually decreases and then slightly increases. In general, the wave velocity tends to increase with the increase in back pressure and the decrease of gas influx rate and angular frequency, significantly in low range. Taking the virtual mass force into account, the dispersion characteristic of the pressure wave weakens obviously, especially at the position close to the wellhead.


---

## Body

## 1. Introduction

One of the future trends of the petroleum industry is the exploration and development of high pressure, low permeability reservoirs [1]. Drilling-related issues such as excessive mud cost, wellbore ballooning/breathing, kick-detection limitations, difficulty in avoiding gross overbalance conditions, differentially stuck pipe, and resulting well-control issues together contribute to the application of managed pressure drilling (MPD) technology [2]. MPD technology has the ability to quickly react to the expected drilling problems and formations pressures uncertainties, reduce nonproductive time and mitigating drilling hazards, and offer a considerable amount of tangible benefits while drilling in extremely narrow fracture/pore pressure windows [3]. It allows drilling operations to proceed where conventional drilling is easy to cause formation damage or considered uneconomical, of  high risk, or even impossible [4]. Although drilling operators try to avoid the risk of influxes, occasionally there are influxes for various reasons. Gas influx occurs whenever the pressure of a gas-bearing formation exceeds the pressure at the bottom of a wellbore. Since the subsequent intrusion of gas displaces drilling mud, it decreases the pressure in the wellbore and makes gas enter even faster [5, 6]. If it is not counteracted in time, the unstable effect can escalate into a blowout creating severe financial losses, environmental contamination, and potential loss of human lives. The basic principle of MPD well control is to keep the bottomhole pressure (BHP) as constant as possible at a value that is at least equal to the formation pressure [7]. MPD is a class of techniques that allow precise management of BHP under both static and dynamic conditions through a combination of controlling the flow rate, mud density, and back pressure (or wellhead pressure) on the fluid returns (choke manifold) of the enclosed and pressurized fluid system [8, 9]. As a result of the ability to manipulate the back pressure, MPD offers the capability to improve safety and well control through early detection of influxes and losses in a microflux level, reduces the risk of influx and thus the chance of blowouts, and controls an influx dynamically without conventionally shutting-in [10].During the managed pressure drilling process, all the unsteady operations such as adjustment of choke, wellhead back pressure controls [11], tripping in/out [12], shutting-in [13], and mud pump rate changes [14] will cause generation and propagation of pressure waves, which would threaten the whole drilling system from the wellhead facilities to the bottomhole drilling equipment and the formation [15]. In the analysis of an influx well and formulation of well control scheme, the dynamic effects of these operations, appeared as pressure fluctuations, should be accounted for [16]. As a basic parameter of pressure fluctuation, the pressure wave velocity has close relevance with the determination of transient pressure and safe operating parameters for well control. However, the gas influx is more troublesome for the higher compressibility and lower density of influx gas than the single phase drilling mud [17]. The works of Bacon et al. [18] had demonstrated that compressibility effects of a gas influx can be significant during an applied-back-pressure, dynamic, MPD well control response and can impact the well control process. Hence, this paper considered the propagation behavior of pressure wave in gas-drilling mud two-phase flow in the annulus to provide reference for the MPD operations.Well control includes not only the handling but also early detection of a gas influx. Besides the transient pressure problem mentioned previously, investigation of the propagation characteristics of pressure wave is of great significance to the early detection of gas influx [19]. Many scholars believe that the pressure fluctuations procedure contains a wealth of information about the flow. Thus, characteristics of pressure wave can be easily used to measure some important parameters in the two-phase flow [20]. In the late 1970s, the former Soviet All-Union Drilling Technology Research Institute began to study characteristics of pressure wave propagation velocity in gas-liquid two-phase flow to detect early gas influx and achieved some important results [21]. According to the functional relationship between the pressure wave propagation velocity in gas-liquid two-phase flow and the gas void fraction, Li et al. [22] presented a method of detecting the gas influxes rate and the height of gas migration early after gas influxes into the wellbore. Furthermore, mud pulse telemetry is the most common method of data transmission used by measurement-while-drilling, and the transmission velocity of the pulse is a basic parameter for this kind of data transmission mode [23].The pressure wave discussed in this paper can be transmitted as a pressure perturbation along the direction of flow in wellbore, which propagates with the speed of sound in the mud and gas two-phase drilling fluid. Due to the compressibility of the gas phase, the changes in interface between the gas and drilling mud, and the momentum and energy transfer between two phases, it is complicated to predict the pressure wave velocity in gas and drilling mud two-phase flow. Since the 1940s, many experimental and theoretical studies have been performed. Experimental tests were conducted to inspect the contributions of fluctuation and flow characteristics on pressure wave. Ruggles et al. [24] firstly performed the experimental investigation on the dispersion property of pressure wave propagation in air-water bubbly flow. It was demonstrated that the propagation speed of pressure wave varies over a range of values for the given state, depending on the angular frequency of the pressure wave. Legius et al. [25] tested the propagation of pressure waves in bubbly and slug flow. The experimental result is similar to the calculation result of the Nguyen model and simulation result of the Sophy-2 package. Concluded form Miyazaki and Nakajima experiments [26] in Nitrogen-Mercury two-phase system, the slip between the phases plays a very important role in the mechanism of pressure wave propagation. From experimental investigation, Bai [27] found that the fractal dimension, correlation dimension, and the Kolmogorov entropy have close relationship with flow regime, and the fractal dimension will be greater than 1.5 when the flow is annular with high gas velocity. The characteristics of pressure wave propagation in bubbly and slug flow in a vertical pipe were investigated experimentally by Huang et al. in detail [28]. It confirmed that the propagation velocity is greatly affected by the gas void fraction and angular frequency of the pressure disturbance, and the superficial velocity of flowing medium has almost no effect on the propagation velocity. Also, there are some widely accepted models including elasticity model, homogeneous mixture model, and continuum model for pressure wave propagation in gas-liquid two-phase flow. Tangren et al. [29] took the two-phase media as a homogeneous fluid and presented the solution model concerning the problem of pressure wave velocity in two-phase flow at low gas void fraction. Wallis [30] studied the propagation mechanism of pressure waves and derived the propagation velocity in a bubbly flow and separated flow using the homogeneous model, in which the two-phase mixture is treated as a compressible fluid with suitably averaged properties. Nguyen et al. [31] applied the elastic theory to predict the propagation velocity of pressure waves in several different flow regimes. The comparison between the calculation results and available experimental data suggests its success at low void fraction. Mecredy and Hamilton [32] derived a detailed continuum model for sound wave propagation in gas-liquid flow by using six separate conservation equations to describe the flow of the vapor and liquid phases. This so-called two-fluid representation allows for nonequilibrium mass, heat, and momentum transfer between the phases. Results indicated that in a bubbly flow, high angular frequency waves travelled an order of magnitude faster than low angular frequency waves. With the development of hydrodynamics, the two-fluid model is widely used to determine the propagation velocity of the pressure wave in two-phase flow as it can provide a general dispersion relation valid for arbitrary flow regimes including effects of the interphase mass, heat, and momentum transfer. Ardron and Duffey [33] developed a model for sound-wave propagation in nonequilibrium vapor-liquid flows which predicts sound speeds and wave attenuations dependent only on measurable flow properties on the basis of two-fluid conservation equations. Ruggles et al. [24], Xu and Chen [34], Chung et al. [35], Huang et al. [28], and Bai et al. [27] investigated the propagation velocity behavior of pressure wave via two-fluid model and small perturbation theory, and the predicted results show good agreement with the experimental data. In recent years, some new researches and models that are especially important in this area were developed. Xu and Gong [36] proposed a united model to predict wave velocity for different flow patterns. In this united model, the effect of a virtual mass coefficient was taken into consideration. The propagation of pressure wave during the condensation of R404A and R134a refrigerants in pipe mini-channels that undergo periodic hydrodynamic disturbances was given by Kuczynski [37, 38]. Li et al. [39] simulated the condensation of gas oxygen in subcooled liquid oxygen and the corresponding mixing process in pump pipeline with the application of thermal phase change model in Computational Fluid Dynamics code CFX and investigated the pressure wave propagation characteristics in two-phase flow pipeline for liquid-propellant rocket based on a proposed pressure wave propagation model and the predicted flow parameters. Based on the unified theory of Kanagawa et al. [40], the nonlinear wave equation for pressure wave propagation in liquids containing gas bubbles is rederived. On the basis of numerical simulation of the gaseous oxygen and liquid oxygen condensation process with the thermal phase model in ANSYS CFX, Chen et al. [41] solved the pressure wave propagation velocity in pump pipeline via the dispersion equation derived from ensemble average two-fluid model. Li and He [42] developed an improved slug flow tracking model and analyzed the variation rule of the pressure wave along the pipeline and influence of the variation of initial inlet liquid flow rate and gas flow rate in horizontal air-water slug flow with transient air flow rate. Meanwhile, the compressibility effects of gas had been noticed in the research field of propagation of pressure waves in the drilling industry, and some efforts have been made. Li et al. [22] established the relationship between wave velocity and gas void fraction according to the empirical formula of the homogeneous mixture model presented by Martin and Padmanatbhan [43] and frequency response model presented by Henry [44]. By applying the unsteady flow dynamic theory, Liu et al. [45] derived the pressure wave velocity calculation formula for gas-drilling mud-solid three-phase flow based on the continuity equation. Wang and Zhang [46] studied the pressure pulsation in mud and set up a model for calculating the amplitude of pressure pulsation when pressure wave is transmitted in drilling-fluid channel especially drilling hose with different inside diameters. However some efforts have been made; the pressure wave velocity is usually determined by empirical formula. In the past researches, the influencing factors for pressure wave propagation were simulated and analyzed with the mathematical model; however, the variation of wave velocity and gas void at different depth of wellbore was not considered. In addition, the current researches are limited in their assumption and neglect the flow pattern translation and interphase forces along the annulus. Up to now, no complete mathematical model of pressure wave in an annulus with variations of gas void, flow pattern, temperature, and back pressure during MPD operations has been derived.The object of the present work is to study the velocity for the transmission of pressure disturbance in the two-phase drilling fluid in the form of a pressure wave in annulus during MPD operations. In this paper, in addition to the pressure, temperature, and the void fraction in the annulus, the compressibility of the gas phase, the virtual mass force, and the changes of interface in two phases are also taken into consideration. By introducing the pressure gradient equations in MPD operations, gas-liquid two-fluid model, the gas-drilling mud equations of state (EOS), and small perturbation theory, a united model for predicting pressure wave velocity in gas and drilling mud in an annulus is developed. The model can be used to predict the wave velocity of various annulus positions at different influx rates, applied back pressures, and angular frequencies with a full consideration of drilling mud compressibility and interphase forces.


## 2. The Mathematical Model

### 2.1. The Basic Equation

In this paper, the two-fluid model and the pressure gradient equation along the flow direction in the annulus are combined to study the pressure wave velocity in MPD operations. Drilling fluid contains clay, cuttings, barite, other solids, and so forth. The solid particles are small and uniformly distributed; therefore, drilling fluid is considered to be a pseudohomogeneous liquid, and the influx natural gas is considered to be the gas phase. The following assumptions are made:(i)
the two-phase flow is treated as one dimension;(ii)
no mass transfers between the gas and drilling mud;(iii)
the flow pattern in an annulus is either bubbly or slug flow.As shown in Figure1, the gas and drilling mud two-phase fluid travels along the annulus in the drilling process. The fluid flows along the annulus in “-z” direction, and the annulus is formed by the casing and drill string.Figure 1
The schematic of gas and drilling mud two-phase flow in annulus.The momentum conservation equation for gas and liquid two-phase flow is(1)∑Fz=ρmAdzdvmdt.The mass conservation equations are(2)∂∂t(ϕGρG)+∂∂z(ϕGρGvG)=0,∂∂t(ϕLρL)+∂∂z(ϕLρLvL)=0.The momentum conservation equations for gas are(3)∂∂t(ϕGρGvG)+∂∂z(ϕGρGvG2)+∂∂z(ϕGρG)-∂∂z[ϕG(τGfr+τGRe)+MGi-4τGD]=0.The momentum conservation equations for liquid are(4)∂∂t(ϕLρLvL)+∂∂z(ϕLρLvL2)+∂∂z(ϕLρL)-∂∂z[ϕL(τLfr+τLRe)+MLi-4τLD]=0.The interphase forces include virtual mass force, drag force, and the wall shear stress. The transfer of momentum between the gas and drilling mud phasesMGi and MLi can be written as follows:
(5)MGi=-MLind-MLid+(τLifr+τLiRe)∂ϕL∂z+∂(ϕGσs)∂z+∂(ϕGpGi)∂z-ϕG∂(pLi)∂z,MLi=MLind+MLid+pLi∂(ϕL)∂z-(τLifr+τLiRe)∂ϕL∂z.The momentum transfer term caused by the gas-liquid interfacial relative acceleration motion (i.e., the virtual mass force) can be expressed as(6)MLind=cVMϕGρLαVM-0.1ϕGρLvs∂vs∂z-cm1ρLvs2∂ϕG∂z,
where cm1=0.1 and vs is the slip velocity; vs=vG-vL.The momentum transfer term caused by the drag force provided by Park et al. [47] can be described as
(7)MLid=38CDrρLRGvs2.The pressure difference between the liquid interface and liquid is defined by the following formula presented by Park et al.:(8)pLi-pL=-cpρLvs2,
where cp=0.25.The proportion of gas phase in the interface between the gas and drilling mud is rather small in that the pressure difference between the gas interface and gas is not very high. Omitting the pressure difference, the gas interface pressurePGi can be written as
(9)pGi-pG≈0.According to the equations of Arnold [48], the pressure of the drilling mud is described as follows:
(10)pL=p-0.25ρLϕGvs2.In fact, the shear stress and the interphase shear stress are very small relative to the Reynolds stress. Also, the Reynolds stress in gas phase can be omitted relative to the interphase force. Hence, it can be described as(11)τGfr≈τLifr≈τLfr≈τG≈τGRe≈0.The Reynolds stress and interfacial average Reynolds stress can be obtained by(12)τLRe=-crρLvs2ϕGϕL,τLiRe=-crρLvs2,
where cr=0.2.The wall shear stress of liquid phase is given by Wallis in the following form [30]:
(13)τL=0.5fLρLvL2.


### 2.2. Physical Equations

#### 2.2.1. Equations of State for Gas

The equation of state (EOS) for gas can be expressed as follows:(14)ρG=p(ZG·R·T),
where ZG is the compression factor of gas.The formula presented by Dranchuk and Abou-Kassem has been used to solve the gas deviation factor under the condition of low and medium pressure(p<35MPa) [49]:
(15)ZG=1+(0.3051-1.0467Tr-0.5783Tr3)ρr+(0.5353-02.6123Tr-0.6816Tr3)ρr2,
where Tr=T/Tc, pr=p/pc, and ρr=0.27pr/ZGTr.The formula presented by Yarborough and Hall has been adopted to solve the gas deviation factor under the condition of high pressure(p≥35MPa) [50]:
(16)ZG=0.06125PrTr-1exp(-1.2(1-Tr-1)2)Y,
where Y is given by
(17)-0.06125prTr-1exp(-1.2(1-Tr-1)2)+Y+Y2+Y3+Y4(1-Y)3=(14.76Tr-1-9.76Tr-2+4.58Tr-3)Y2-(90.7Tr-1-242.2Tr-2+42.4Tr-3)Y(2.18+2.82Tr-1).


#### 2.2.2. Equations of State for Liquid

Under different temperatures and pressures, the density of drilling mud can be obtained by the empirical formulas.IfT<130C°, the density of drilling mud can be obtained by the following equation:
(18)ρL=ρ0(1+4×10-10pL-4×10-5T-3×10-6T2).IfT≥130C°, the density of drilling mud is
(19)ρL=ρ0((T-130T)21+4×10-10PL-4×10-5T-3×10-6T2+0.4(T-130T)2).


#### 2.2.3. Correlation of Temperature Distribution

The temperature of the drilling mud at different depths of the annulus can be determined by the relationship presented by Hasan and Kabir [51]:
(20)T=Tei+F[1-e(zbh-z)/A](-gsinθgcJcpm+N+gTsinθ)+e(zbh-z)/A(Tfbh-Tebh).


### 2.3. Flow Pattern Analysis

Based on the analysis of flow characteristics in the enclosed drilling system, it can be safely assumed that the flow pattern in an annulus is either bubbly or slug flow [52]. The pattern transition criteria for bubbly flow and slug flow given by Orkiszewski are used to judge the flow pattern in the gas-drilling mud two-phase flow [53].For bubbly flow,(21)qGqm<LB.For slug flow,(22)qGqm>LB,NGV<LS,
where qm is the volumetric flow rate of two-phase flow, qm=qG+qL.The dimensionless numbersLs and LB are defined as
(23)Ls=50+36NGVqLqG,LB=1.071-0.7277vm2D,
where
(24)NGV=vs(ρLgσs)0.25.Flow parameters such as void fraction, mixture density, and virtual mass force coefficient are discussed for a specific flow pattern.The correlation between void fraction and liquid holdup is expressed as(25)ϕG+ϕL=1.


#### 2.3.1. Bubbly Flow

As for bubbly flow, the density of gas and drilling mud mixture is described as the following formula:(26)ρm=ϕLρL+ϕGρG.The void fraction is determined by the following formula:(27)ϕG=12[1+qmvsA-(1+qmvsA)2-4qGvsA].The coefficient of virtual mass forceCvm for bubbly flow can be expressed as follows [54]:
(28)Cvm=0.51+2ϕG1-ϕG.The coefficient of resistance coefficientCD for bubbly flow can be expressed as
(29)CD=4Rb3g(ρL-ρG)σs[1+17.67(1-ϕG)9/718.67(1-ϕG)1.5]2.The friction pressure gradient for bubbly flow can be obtained from the following equation:(30)τf=fρLvL22D.


#### 2.3.2. Slug Flow

As for slug flow, distribution coefficient of gas in the liquid phase is(31)C0=0.00252lg(103μL)A1.38-0.782+0.232lgvm-0.428lgA.The average density of the mixture for slug flow is determined by(32)ρm=Wm+ρLvsAqm+vsA+C0ρL.The void fraction for slug flow is determined by the following formula:(33)ϕG=qGqG+qL.The coefficient of virtual mass forceCvm for slug flow can be expressed as follows:
(34)Cvm=3.3+1.73Lq-3Rq3Lq-Rq.The coefficient of resistance coefficientCD for slug flow can be expressed as
(35)CD=110ϕL3Rb.The friction pressure gradient of bubbly flow can be obtained from the following equation:(36)τf=fρLvm22D(qL+vsAqm+vsA+C0).


### 2.4. Annulus Characteristic Analysis

The annulus effective diameter proposed by Sanchez [55] is used in the hydraulic calculation of annulus
(37)D=π(Do2-Di2)/4π(Do+Di)/4=Do-Di.The effective roughness of the annulus can be expressed as(38)ke=koDoDo+Di+kiDiDo+Di,
where the Di and Do are the diameters of inner pipe and outer pipe, respectively; the k0 and ki are the roughnesses of outer pipe and the inner pipe, respectively.


## 2.1. The Basic Equation

In this paper, the two-fluid model and the pressure gradient equation along the flow direction in the annulus are combined to study the pressure wave velocity in MPD operations. Drilling fluid contains clay, cuttings, barite, other solids, and so forth. The solid particles are small and uniformly distributed; therefore, drilling fluid is considered to be a pseudohomogeneous liquid, and the influx natural gas is considered to be the gas phase. The following assumptions are made:(i)
the two-phase flow is treated as one dimension;(ii)
no mass transfers between the gas and drilling mud;(iii)
the flow pattern in an annulus is either bubbly or slug flow.As shown in Figure1, the gas and drilling mud two-phase fluid travels along the annulus in the drilling process. The fluid flows along the annulus in “-z” direction, and the annulus is formed by the casing and drill string.Figure 1
The schematic of gas and drilling mud two-phase flow in annulus.The momentum conservation equation for gas and liquid two-phase flow is(1)∑Fz=ρmAdzdvmdt.The mass conservation equations are(2)∂∂t(ϕGρG)+∂∂z(ϕGρGvG)=0,∂∂t(ϕLρL)+∂∂z(ϕLρLvL)=0.The momentum conservation equations for gas are(3)∂∂t(ϕGρGvG)+∂∂z(ϕGρGvG2)+∂∂z(ϕGρG)-∂∂z[ϕG(τGfr+τGRe)+MGi-4τGD]=0.The momentum conservation equations for liquid are(4)∂∂t(ϕLρLvL)+∂∂z(ϕLρLvL2)+∂∂z(ϕLρL)-∂∂z[ϕL(τLfr+τLRe)+MLi-4τLD]=0.The interphase forces include virtual mass force, drag force, and the wall shear stress. The transfer of momentum between the gas and drilling mud phasesMGi and MLi can be written as follows:
(5)MGi=-MLind-MLid+(τLifr+τLiRe)∂ϕL∂z+∂(ϕGσs)∂z+∂(ϕGpGi)∂z-ϕG∂(pLi)∂z,MLi=MLind+MLid+pLi∂(ϕL)∂z-(τLifr+τLiRe)∂ϕL∂z.The momentum transfer term caused by the gas-liquid interfacial relative acceleration motion (i.e., the virtual mass force) can be expressed as(6)MLind=cVMϕGρLαVM-0.1ϕGρLvs∂vs∂z-cm1ρLvs2∂ϕG∂z,
where cm1=0.1 and vs is the slip velocity; vs=vG-vL.The momentum transfer term caused by the drag force provided by Park et al. [47] can be described as
(7)MLid=38CDrρLRGvs2.The pressure difference between the liquid interface and liquid is defined by the following formula presented by Park et al.:(8)pLi-pL=-cpρLvs2,
where cp=0.25.The proportion of gas phase in the interface between the gas and drilling mud is rather small in that the pressure difference between the gas interface and gas is not very high. Omitting the pressure difference, the gas interface pressurePGi can be written as
(9)pGi-pG≈0.According to the equations of Arnold [48], the pressure of the drilling mud is described as follows:
(10)pL=p-0.25ρLϕGvs2.In fact, the shear stress and the interphase shear stress are very small relative to the Reynolds stress. Also, the Reynolds stress in gas phase can be omitted relative to the interphase force. Hence, it can be described as(11)τGfr≈τLifr≈τLfr≈τG≈τGRe≈0.The Reynolds stress and interfacial average Reynolds stress can be obtained by(12)τLRe=-crρLvs2ϕGϕL,τLiRe=-crρLvs2,
where cr=0.2.The wall shear stress of liquid phase is given by Wallis in the following form [30]:
(13)τL=0.5fLρLvL2.


## 2.2. Physical Equations

### 2.2.1. Equations of State for Gas

The equation of state (EOS) for gas can be expressed as follows:(14)ρG=p(ZG·R·T),
where ZG is the compression factor of gas.The formula presented by Dranchuk and Abou-Kassem has been used to solve the gas deviation factor under the condition of low and medium pressure(p<35MPa) [49]:
(15)ZG=1+(0.3051-1.0467Tr-0.5783Tr3)ρr+(0.5353-02.6123Tr-0.6816Tr3)ρr2,
where Tr=T/Tc, pr=p/pc, and ρr=0.27pr/ZGTr.The formula presented by Yarborough and Hall has been adopted to solve the gas deviation factor under the condition of high pressure(p≥35MPa) [50]:
(16)ZG=0.06125PrTr-1exp(-1.2(1-Tr-1)2)Y,
where Y is given by
(17)-0.06125prTr-1exp(-1.2(1-Tr-1)2)+Y+Y2+Y3+Y4(1-Y)3=(14.76Tr-1-9.76Tr-2+4.58Tr-3)Y2-(90.7Tr-1-242.2Tr-2+42.4Tr-3)Y(2.18+2.82Tr-1).


### 2.2.2. Equations of State for Liquid

Under different temperatures and pressures, the density of drilling mud can be obtained by the empirical formulas.IfT<130C°, the density of drilling mud can be obtained by the following equation:
(18)ρL=ρ0(1+4×10-10pL-4×10-5T-3×10-6T2).IfT≥130C°, the density of drilling mud is
(19)ρL=ρ0((T-130T)21+4×10-10PL-4×10-5T-3×10-6T2+0.4(T-130T)2).


### 2.2.3. Correlation of Temperature Distribution

The temperature of the drilling mud at different depths of the annulus can be determined by the relationship presented by Hasan and Kabir [51]:
(20)T=Tei+F[1-e(zbh-z)/A](-gsinθgcJcpm+N+gTsinθ)+e(zbh-z)/A(Tfbh-Tebh).


## 2.2.1. Equations of State for Gas

The equation of state (EOS) for gas can be expressed as follows:(14)ρG=p(ZG·R·T),
where ZG is the compression factor of gas.The formula presented by Dranchuk and Abou-Kassem has been used to solve the gas deviation factor under the condition of low and medium pressure(p<35MPa) [49]:
(15)ZG=1+(0.3051-1.0467Tr-0.5783Tr3)ρr+(0.5353-02.6123Tr-0.6816Tr3)ρr2,
where Tr=T/Tc, pr=p/pc, and ρr=0.27pr/ZGTr.The formula presented by Yarborough and Hall has been adopted to solve the gas deviation factor under the condition of high pressure(p≥35MPa) [50]:
(16)ZG=0.06125PrTr-1exp(-1.2(1-Tr-1)2)Y,
where Y is given by
(17)-0.06125prTr-1exp(-1.2(1-Tr-1)2)+Y+Y2+Y3+Y4(1-Y)3=(14.76Tr-1-9.76Tr-2+4.58Tr-3)Y2-(90.7Tr-1-242.2Tr-2+42.4Tr-3)Y(2.18+2.82Tr-1).


## 2.2.2. Equations of State for Liquid

Under different temperatures and pressures, the density of drilling mud can be obtained by the empirical formulas.IfT<130C°, the density of drilling mud can be obtained by the following equation:
(18)ρL=ρ0(1+4×10-10pL-4×10-5T-3×10-6T2).IfT≥130C°, the density of drilling mud is
(19)ρL=ρ0((T-130T)21+4×10-10PL-4×10-5T-3×10-6T2+0.4(T-130T)2).


## 2.2.3. Correlation of Temperature Distribution

The temperature of the drilling mud at different depths of the annulus can be determined by the relationship presented by Hasan and Kabir [51]:
(20)T=Tei+F[1-e(zbh-z)/A](-gsinθgcJcpm+N+gTsinθ)+e(zbh-z)/A(Tfbh-Tebh).


## 2.3. Flow Pattern Analysis

Based on the analysis of flow characteristics in the enclosed drilling system, it can be safely assumed that the flow pattern in an annulus is either bubbly or slug flow [52]. The pattern transition criteria for bubbly flow and slug flow given by Orkiszewski are used to judge the flow pattern in the gas-drilling mud two-phase flow [53].For bubbly flow,(21)qGqm<LB.For slug flow,(22)qGqm>LB,NGV<LS,
where qm is the volumetric flow rate of two-phase flow, qm=qG+qL.The dimensionless numbersLs and LB are defined as
(23)Ls=50+36NGVqLqG,LB=1.071-0.7277vm2D,
where
(24)NGV=vs(ρLgσs)0.25.Flow parameters such as void fraction, mixture density, and virtual mass force coefficient are discussed for a specific flow pattern.The correlation between void fraction and liquid holdup is expressed as(25)ϕG+ϕL=1.


### 2.3.1. Bubbly Flow

As for bubbly flow, the density of gas and drilling mud mixture is described as the following formula:(26)ρm=ϕLρL+ϕGρG.The void fraction is determined by the following formula:(27)ϕG=12[1+qmvsA-(1+qmvsA)2-4qGvsA].The coefficient of virtual mass forceCvm for bubbly flow can be expressed as follows [54]:
(28)Cvm=0.51+2ϕG1-ϕG.The coefficient of resistance coefficientCD for bubbly flow can be expressed as
(29)CD=4Rb3g(ρL-ρG)σs[1+17.67(1-ϕG)9/718.67(1-ϕG)1.5]2.The friction pressure gradient for bubbly flow can be obtained from the following equation:(30)τf=fρLvL22D.


### 2.3.2. Slug Flow

As for slug flow, distribution coefficient of gas in the liquid phase is(31)C0=0.00252lg(103μL)A1.38-0.782+0.232lgvm-0.428lgA.The average density of the mixture for slug flow is determined by(32)ρm=Wm+ρLvsAqm+vsA+C0ρL.The void fraction for slug flow is determined by the following formula:(33)ϕG=qGqG+qL.The coefficient of virtual mass forceCvm for slug flow can be expressed as follows:
(34)Cvm=3.3+1.73Lq-3Rq3Lq-Rq.The coefficient of resistance coefficientCD for slug flow can be expressed as
(35)CD=110ϕL3Rb.The friction pressure gradient of bubbly flow can be obtained from the following equation:(36)τf=fρLvm22D(qL+vsAqm+vsA+C0).


## 2.3.1. Bubbly Flow

As for bubbly flow, the density of gas and drilling mud mixture is described as the following formula:(26)ρm=ϕLρL+ϕGρG.The void fraction is determined by the following formula:(27)ϕG=12[1+qmvsA-(1+qmvsA)2-4qGvsA].The coefficient of virtual mass forceCvm for bubbly flow can be expressed as follows [54]:
(28)Cvm=0.51+2ϕG1-ϕG.The coefficient of resistance coefficientCD for bubbly flow can be expressed as
(29)CD=4Rb3g(ρL-ρG)σs[1+17.67(1-ϕG)9/718.67(1-ϕG)1.5]2.The friction pressure gradient for bubbly flow can be obtained from the following equation:(30)τf=fρLvL22D.


## 2.3.2. Slug Flow

As for slug flow, distribution coefficient of gas in the liquid phase is(31)C0=0.00252lg(103μL)A1.38-0.782+0.232lgvm-0.428lgA.The average density of the mixture for slug flow is determined by(32)ρm=Wm+ρLvsAqm+vsA+C0ρL.The void fraction for slug flow is determined by the following formula:(33)ϕG=qGqG+qL.The coefficient of virtual mass forceCvm for slug flow can be expressed as follows:
(34)Cvm=3.3+1.73Lq-3Rq3Lq-Rq.The coefficient of resistance coefficientCD for slug flow can be expressed as
(35)CD=110ϕL3Rb.The friction pressure gradient of bubbly flow can be obtained from the following equation:(36)τf=fρLvm22D(qL+vsAqm+vsA+C0).


## 2.4. Annulus Characteristic Analysis

The annulus effective diameter proposed by Sanchez [55] is used in the hydraulic calculation of annulus
(37)D=π(Do2-Di2)/4π(Do+Di)/4=Do-Di.The effective roughness of the annulus can be expressed as(38)ke=koDoDo+Di+kiDiDo+Di,
where the Di and Do are the diameters of inner pipe and outer pipe, respectively; the k0 and ki are the roughnesses of outer pipe and the inner pipe, respectively.


## 3. The United Model Developed

The total pressure drop gradient is the sum of pressure drop gradients due to potential energy change and kinetic energy and frictional loss. From (1), the equation used to calculate the pressure gradient of gas and drilling mud flow within the annulus can be written as
(39)dpdz=ρmgsinθ-τwπDA-ρmvmdvmdz.Assuming the compressibility of the gas is only related to the pressure in the annulus, the kinetic energy or acceleration term in the previous equation can be simplified to(40)ρmvmdvmdz=-ρmvmvsgpdpdz=-WmqgA2pdpdz.Substituting (39) into (40), the total pressure drop gradient along the flow direction within the annulus can be expressed as
(41)dpdz=ρmg+τf1-WmqG/(A2p).It is assumed that the gas obeys to the EOS (14), and the compressibility of drilling mud can be obtained by adopting the simplified EOS ((18) and (19)) which neglects the thermal expansion of liquid. The sonic speed of gas phase cG and that of liquid phase cL can be presented in the following form:
(42)dpLdρL=cL2,dpGdρG=cG2.By introducing (42), the hydrodynamic equations of two-fluid model ((2)–(4)) can be written in the matrix form
(43)A∂ξ∂t+B∂ξ∂z=Cξ,
where A is the matrix of parameters considered in relation to time, B is the matrix of parameters considered in relation to the spatial coordinate, and C is the vector of extractions.By introducing the small disturbance theory, the disturbance of the state variableξ(ϕG,p,vG,vL)T can be written as
(44)ξ=ξ0+δξexp[i(wt-kt)].In (44), k is the wave number, and w is the angular frequency of the disturbances. Substituting (44) into (43) gives homogenous linear equations concerning the expression (δϕG,δp,δvG,δvL)T. According to the solvable condition of the homogenous linear equations that the determinant of the equations is zero, dispersion equation of pressure wave can be expressed in the following form:
(45)|M1M2M3M4-ρLw1-ϕGcL2w0-k(1-ϕG)ρLM5M6M7M8M9M10M8M11|=0,
where M1–M11 can be illustrated by
(46)M1=(ρG+cpϕGρLvs2cG2)w,M2=ϕGcG2[1-cpϕL]vs2cL2w,M3=-[ϕGρGk+2cpϕGϕLρLvscL2w],M4=22cpϕGϕLρLvscL2w,M5=ρLvr2k(-ϕGcp+cr-ci+cm2),M6=-ϕGk[1-ϕLcpvs2cL2+civs2cL2],M7=ϕG(ρG+cvmρL)w-i(34cDrρLϕGvs+4DfGwρGvG),M8=-cvmϕGρLw+i(34cDrρLϕGvs),M9=ρLvs2k(ϕLcp-2cr-cm2),M10=-k(ϕL+crϕGvs2cL2),M11=ρL[ϕL+ϕGcvm]w-i(34cDrρLϕGvs+4DfLρLvL),
where cm2=0.1, cp=0.25, ci=0.3, and cr=0.2.By solving the dispersion equation indicated previously, four roots can be obtained. Because the wavelengths associated with two of the roots are too short to allow the two-fluid medium to be treated as a continuum, the two roots should be omitted. As for the two remaining roots, one of them represents a pressure wave that transmited along the axisz, and the other represents a pressure wave transmits in the opposite direction in accordance with the direction of the flow along the annulus. The real value of wave number can be used to determine the wave velocity c. The wave velocity in the two-phase fluid can be determined by the following model:
(47)c=|(w/R+(k))-(w/R-(k))|2.


## 4. Solution of the United Model

Obtaining the analytical solution of the mathematical models concerned with flow pattern, void fraction, characteristic parameters, and pressure drop gradient is generally impossible for two-phase flow. In this paper, the Runge-Kutta method (R-K4) is used to discretize the theoretical model.We can obtain pressure, temperature, gas velocity, drilling mud velocity, and void fraction at different annulus depths by R-K4. The solution of pressure drop gradient equation (41) can be seen as an initial-value problem of the ordinary differential equation:
(48)dpdz=F(z,p),p(z0)=P0.With the initial value (z0,p0) and the function F(z,p), (50)–(53) can be obtained by
(49)k1=F(z0,p0),(50)k2=F(z0+h2,p0+h2k1),(51)k3=F(z0+h2,p0+h2k2),(52)k4=F(z0+h,p0+hk3),
where h is the step of depth. The pressure on the nod i=i+1 can be obtained by
(53)p1=p0+Δp=p0+h6(k1+2k2+2k3+k4).In the present work, the mathematical model and pressure wave velocity calculation model are solved by a personally compiled code on VB.NET (Version 2010). The solution procedure for the wave velocity in the annulus is shown in Figure1. At initial time, the wellhead back pressure, wellhead temperature, wellbore structure, well depth, gas and drilling mud properties, and so forth, are known. On the node i, the pressure gradient, temperature, and the void fraction can be obtained by adopting R-K4. Then, the determinant (45) is calculated based on the calculated parameters. Omitting the two unreasonable roots, the pressure wave velocity at different depths of the annulus in MPD operations can be solved by (47). The process is repeated until the pressure wave velocity of every position in the wellbore is obtained  as shown in Figure 2.Figure 2
Solution procedure for wave velocity in MPD operations.The developed model takes full consideration of the interfacial interaction and the virtual mass force. Owing to the complex conditions of the annulus in MPD operations, measurement of wave velocity in the actual drilling process is very difficult. In order to verify the united model, the predicted pressure waves are compared with the results of previous simulated experimental investigations presented for gas and drilling mud by Liu et al. [45] in Figure 3(a) and by Li et al. [22] in Figure 3(b). The lines represent the calculation results, and the points represent the experimental data.Experimental verification by comparison with previous experimental data.
(a)
Comparison with the experimental data of Liu
(b)
Comparison with the experimental data of LiThe comparisons reveal that the developed united model fits well with the experimental data. Thus, the united model can be used to accurately predict wave velocity at different wellhead back pressures and gas influx rates (the influx rate of gas in the bottomhole) in MPD operations.


## 5. Analysis and Discussion

The drilling system described is an enclosed system. The schematic diagram of the gas influx process is illustrated in Figure4. The drilling mud is pumped from surface storage down the drill pipe. Returns from the wellbore annulus travel back through surface processing, where drilling solids are removed, to surface storage. The key equipments include the following.(i)
The Rotating Control Device (RCD). The rotating control device provides the seal between atmosphere and wellbore, while allowing pipe movement and diverting returns flow. In conjunction with the flow control choke, the RCD provides the ability to apply back pressure on the annulus during an MPD operation.(ii)
Choke. The MPD choke manifold provides an adjustable choke system which is used to dynamically control the required BHP by means of applying surface BP.(iii)
Coriolis Meter. A Coriolis meter is used to accurately measure the mass flow rate of fluid exiting the annulus. The ability to measure return flow accurately is essential for the applied back pressure.(iv)
Pressure Sensor. A pressure sensor is used to measure surface back pressure on the wellhead.Figure 4
Schematic diagram of the gas influx process.The gas and drilling mud flow rate measured by the Coriolis meter and the back pressure measured by the pressure sensor are the initial data for annulus pressure calculation. The well used for calculation is a gas well in Xinjiang Uygur Autonomous Region, Northwest China. The wellbore structure, well design parameters (depths and diameters), gas and drilling mud properties (density and viscosity), and operational conditions of calculation well are displayed in Table1.Table 1
Parameters of calculation well.
Type
Property
Value
Mud
Dynamic viscosity (Pa · s)
0.056
Density (kg/cm3)
1460
Gas
Relative density
0.65
Viscosity (Pa · s)
1.14 × 10−5
String
Elastic modulus of string (Pa)
2.07 × 1011
Poisson ratio of string
0.3
Roughness (m)
1.54 × 10−7
Surface condition
Surface temperature (K)
298
Atmosphere pressure (MPa)
0.101The drilling mud mixed with gas is taken as a two-phase flow medium. The propagation velocity of pressure wave in the gas-drilling mud flow is calculated and discussed by using the established model and well parameters.


### 5.1. Effect of Back Pressure on Wave Velocity

Increase in applied back pressure is the most common approach for dynamic well control. The wellbore can be seen as an enclosed pressurized system. The pressure at different depths of the annulus varied with the change of back pressure. According to the EOS of gas and drilling mud, the influence of pressure on gas volume is much greater than that of drilling mud for the greater compressibility of gas. So, the gas void fraction changes with the variation of gas volume at the nearly constant flow rate of drilling mud at different annular pressures. Meanwhile, the pressure wave velocity is sensitive to the gas void fraction. As a result, when the back pressure at the wellhead is changed by adjusting the choke, the wave velocity in the two-phase drilling fluid and the distribution of void fraction at different depth of the annulus will diverge. The calculation results affected by the back pressure are presented.Figures5 and 6 show the distributions of void fraction and variations of wave velocity along the flow direction in the annulus when the back pressure at the wellhead is 0.1 MPa, 0.8 MPa, 1.5 MPa, 2.5 MPa, 5.0 MPa, and 6.5 MPa, respectively. It can be seen that the void fraction significantly increases along the flow direction in the annulus. Conversely, the wave velocity shows a remarkable decrease tendency. For instance, BP=1.5 MPa, at the wellhead, the wave velocity is 90.12 m/s, and the void fraction is 0.611, while in the bottomhole, the wave velocity reaches 731.42 m/s, and void fraction is reduced to 0.059. However, a sudden increase in wave velocity is observed in Figure 5 when the back pressure is 0.1 MPa at the position close to wellhead as the void fraction is further increased. Moreover, the influence of back pressure on the void fraction is larger at the position close to the wellhead than that in the bottomhole.Figure 5
Void fraction distribution at different back pressures.Figure 6
Wave velocity variations at different back pressures.This can be explained from the viewpoints of mixture density and compressibility of two-phase fluid and the pressure drop along the flow direction in the wellbore. In the low void fraction range, the gas phase is dispersed in the liquid as bubble, so the wave velocity is influenced greatly by the added gas phase. According to the EOS, if gas invades into the wellbore with a small amount in the bottomhole, the density of the drilling mud has little variation while the compressibility increases obviously, which makes the wave velocity decreased. Then, the two-phase drilling mud flow from the bottomhole to the wellhead along the annulus with a drop of pressure caused by, potential energy change, kinetic energy and frictional loss, which leads to an increase in void fraction. When the void fraction is increased, both the density and the compressibility of the two-phase fluid change slightly, resulting in a flat decrease in wave velocity along the flow direction. With the further increasing in void fraction, the bubbly flow turns into slug flow according to the flow pattern transient criteria of Orkiszewski which shows that the flow pattern is dependent on the void fraction. Generally speaking, the gas-drilling mud slug flow is composed of  liquid and gas slugs. In the preliminary slug flow, the liquid slugs are much longer than the gas slugs, so the wave velocity is determined primarily by the wave velocities in the liquid slugs which are hardly affected by the void fraction. At the position close to the wellhead, the pressure of two-phase flow fluid can be reduced to a low value which approaches the back pressure. The compressibility of gas will be improved at the low pressure. It results in significant increase in void fraction. As the void fraction increases greatly at the position close to the wellhead, gas slugs become much longer than the liquid slugs. It is assumed that the actual wave velocity in slug flow happens to play a leading role when the void fraction increases to some extent. For the decrease in liquid holdup in the gas slug with the increasing in void fraction, the wave velocity in the gas slug is infinitely close to that in pure gas, which is equivalent to an increase in wave velocity in the gas slug. Meanwhile, the gas content in liquid slug also increases and results in a lower wave velocity in the liquid slug. As a result, a slight increase in wave velocity appears.Figures7 and 8 present the influence of back pressure on the void fraction and wave velocity with respect to the parameters of well depth. As the back pressure increases, the void fraction at different depth of the annulus is reduced gradually, while the wave velocity in the two-phase flow tends to increase. Analytical results show that the increased back pressure is equivalent to be applied to the entire enclosed drilling fluid cyclical system. The pressure transmits from the wellhead to the bottomhole; therefore, the annular pressure in the entire wellbore is increased. According to the EOS, the density of gas increases and the compressibility of gas decreases with the increasing of gas pressure. So, the loss of interphase momentum and energy exchange is reduced and the interphase momentum exchange is promoted. It contributes to the increase in wave velocity with the increase in gas pressure. In addition, due to lower compressibility of two-phase flow medium under high pressure, the increase tendency of pressure wave velocity and the decrease tendency of void fraction are slowed down in the high back pressure range.Figure 7
Effect of back pressure on the void fraction.Figure 8
Effect of back pressure on the pressure wave velocity.


### 5.2. Effect of Gas Influx Rate on Wave Velocity

Figures9 and 10 graphically interpret the distributions of void fraction and variations of wave velocity along the flow direction in the annulus. When gas influx occurs in the bottomhole, gas invades into the wellbore and migrates from the bottomhole to the wellhead along the flow direction. At a low gas influx rate, it is extremely obvious that the void fraction and wave velocity first slightly change in a comparatively smooth value then change sharply. It is because of the rapid expansion of gas volume with the decreasing in pressure near the wellhead that the void fraction increases sharply, and the wave velocity decreases obviously at the same time. But under the high bottomhole pressure (up to 50 MPa), the compressibility of the gas is low. This results in a slight change in void fraction and wave velocity at the position far away from the wellhead. Since the compressible component increases with the increase in the gas influx rate, the compressibility of the gas and drilling mud two-phase fluid is improved. So the variations of void fraction and wave velocity become more prominent. Also, the void fraction still shows an increase tendency that is steady first and then sharp. It acts in accordance with the variation of void fraction at a low gas influx rate. At a low gas influx rate, if the void fraction at the position close to the wellhead can not increase to a high extent, the wave velocity always shows a decrease tendency. At a high gas influx rate, such as the Q=8.312 m3/h, the wave velocity tends to increase because the void fraction in the wellhead is increased to a high extent. In conclusion, the wave pressure is sensitive to the void fraction, and the void fraction is dominated by influx rate and pressure in the annulus, especially the influx rate.Figure 9
Void fraction distribution at different gas influx rates.Figure 10
Wave velocity variations at different gas influx rates.With the influx of gas, the mixing of gas and drilling mud occurs in the annulus and the corresponding interfacial transfer of momentum and mass causes an increase in gas phase void fraction and a decrease in pressure wave velocity, as shown in Figures11 and 12. Within the range of low gas influx rate, the wave velocity decreases significantly. It is because of this that the compressibility of the gas increases remarkably, and the medium appears to be of high elasticity, though the density of gas-drilling mud two-phase flow changes slightly. With the increase in the gas influx rate and corresponding increase in the void fraction in the annulus, the compressibility of the two-phase unceasingly increases, which promotes the momentum and energy exchange in the interface. So, the pressure wave continuously decreases. When the void fraction increases to some extent following the increase in the gas influx rate, the decrease in wave velocity in the liquid slug is gradually less than the increase in wave velocity in the gas slug; thus, the decrease of wave velocity is slowed down for the growth of gas slug. Especially in the wellhead, a slight increase in wave velocity is observed.Figure 11
Effect of gas influx rate on the void fraction.Figure 12
Effect of gas influx rate on the pressure wave velocity.


### 5.3. Effect of Virtual Mass Force on Wave Velocity

Figure13 illustrates the effect of virtual mass force on the wave velocity with the parameter of gas influx rate. As shown in Figure 15, under the same gas influx rate (0.594 m3/h, 1.098 m3/h, and 1.768 m3/h, resp.) the wave velocity curve of  Cvm=Re diverges from the curve of  Cvm=0 at the position close to the wellhead, whereas the wave velocity curve of the two types is almost coincided below the position of H=500 m. This divergence between the two types of curves is connected with the fact that significance of influence of interphase virtual mass forces increases together with the increase of the relative motion in the interphase. It is evident that the bottomhole pressure is hundreds of times higher than the wellhead pressure. As a result, the void fraction at the position close to the wellhead increases sharply, meanwhile the interphase momentum and energy exchanges are promoted. The virtual mass force can be described as the transfer of momentum between the gas phase and drilling mud phase caused by the relative motion in the interphase. If the relative motion in the interphase is quite weak, the value of virtual mass force will intend to approach 0, and the influence on the wave velocity can be ignored. However, if the relative motion is rather intense, the effect of virtual mass force on the wave velocity should not be ignored. Furthermore, taking the virtual mass force into account, the dispersion characteristic of the pressure wave weakens obviously. Compared with the pressure wave velocity calculated by ignoring the effect of virtual mass force, the calculated pressure wave velocity with a consideration of the virtual mass force is lower. Therefore, it is necessary to analyze the effect of the virtual mass force on the wave velocity at the position close to the wellhead in MPD operations.Figure 13
Effect of virtual mass force on the pressure wave velocity.


### 5.4. Effect of Disturbance Angular Frequency on Wave Velocity

Figure14 presents the pressure wave velocity in the annulus at different angular frequencies. According to Figures 5 and 9, the void fraction first slightly increases in the bottomhole and then sharply increases along the flow direction at the position close to the wellhead. At a fixed angular frequency, this phenomenon results in an overall decrease in the pressure wave velocity except for the section near the wellhead. At the section close to the wellhead, an opposite change trend of pressure wave velocity is observed for the transition of flow pattern from bubbly flow to slug flow due to the continuous increase in void fraction along the flow direction. In addition, it can be clearly seen from the curves that the wave velocity increases accompanying with the increase in the angular frequency above the position of H=500 m. This property is not very distinct at the position below H=500 m for the low void fraction.Figure 14
Wave velocity variations at different angular frequencies.Figure 15
Effect of angular frequency on wave velocity.Figure15 shows the effect of frequency on the wave velocity in the gas-drilling mud flow (Cvm=0). The curves of the wave velocity of different positions reveal that the propagation velocity of pressure disturbances increases together with the growth of the angular frequency (0<w<500 Hz). It proves that the pressure wave has an obvious dispersion characteristic in the two-phase flow. As the angular frequency increases in the range of less than 500 Hz, there is sufficient time to carry out the exchange of energy and momentum between two phases. It achieves a state of mechanical and thermodynamic equilibrium between the two phases. So the wave velocity increases gradually with the increase in the angular frequency at different depths of the annulus. It is considered that the wave velocity is mainly affected by the interphase mechanical and thermodynamic equilibrium at low anglular frequencies. When the angular frequency reaches the value of w=500 Hz, the pressure wave velocity achieves a constant value and remains on this level regardless of the further growth in angular frequency w. With the increase in angular frequency, there is not enough time for energy and momentum exchange between the gas-drilling mud two phases to reach the mechanical and thermodynamic equilibrium state, and thus the wave dispersion does not exist. At a high angular frequency, the wave velocity is mainly dominated by the mechanical and thermal nonequilibrium in the flow and keeps almost unchanged when the angular frequency is further increased. This is consistent with the influences of the angular frequency in the horizontal pipe [28]. Moreover, the effect of virtual mass force is also shown in this figure (Cvm=Re). It is observed that the wave velocity is significantly reduced when the virtual mass force is taken into consideration in the calculation of the wave velocity at the position close to the wellhead such as H=0 m. It can be explained by the reduction of dispersion characteristic of a pressure wave. Along the wellbore from the wellhead to the bottomhole, the distinction gradually decreases.


## 5.1. Effect of Back Pressure on Wave Velocity

Increase in applied back pressure is the most common approach for dynamic well control. The wellbore can be seen as an enclosed pressurized system. The pressure at different depths of the annulus varied with the change of back pressure. According to the EOS of gas and drilling mud, the influence of pressure on gas volume is much greater than that of drilling mud for the greater compressibility of gas. So, the gas void fraction changes with the variation of gas volume at the nearly constant flow rate of drilling mud at different annular pressures. Meanwhile, the pressure wave velocity is sensitive to the gas void fraction. As a result, when the back pressure at the wellhead is changed by adjusting the choke, the wave velocity in the two-phase drilling fluid and the distribution of void fraction at different depth of the annulus will diverge. The calculation results affected by the back pressure are presented.Figures5 and 6 show the distributions of void fraction and variations of wave velocity along the flow direction in the annulus when the back pressure at the wellhead is 0.1 MPa, 0.8 MPa, 1.5 MPa, 2.5 MPa, 5.0 MPa, and 6.5 MPa, respectively. It can be seen that the void fraction significantly increases along the flow direction in the annulus. Conversely, the wave velocity shows a remarkable decrease tendency. For instance, BP=1.5 MPa, at the wellhead, the wave velocity is 90.12 m/s, and the void fraction is 0.611, while in the bottomhole, the wave velocity reaches 731.42 m/s, and void fraction is reduced to 0.059. However, a sudden increase in wave velocity is observed in Figure 5 when the back pressure is 0.1 MPa at the position close to wellhead as the void fraction is further increased. Moreover, the influence of back pressure on the void fraction is larger at the position close to the wellhead than that in the bottomhole.Figure 5
Void fraction distribution at different back pressures.Figure 6
Wave velocity variations at different back pressures.This can be explained from the viewpoints of mixture density and compressibility of two-phase fluid and the pressure drop along the flow direction in the wellbore. In the low void fraction range, the gas phase is dispersed in the liquid as bubble, so the wave velocity is influenced greatly by the added gas phase. According to the EOS, if gas invades into the wellbore with a small amount in the bottomhole, the density of the drilling mud has little variation while the compressibility increases obviously, which makes the wave velocity decreased. Then, the two-phase drilling mud flow from the bottomhole to the wellhead along the annulus with a drop of pressure caused by, potential energy change, kinetic energy and frictional loss, which leads to an increase in void fraction. When the void fraction is increased, both the density and the compressibility of the two-phase fluid change slightly, resulting in a flat decrease in wave velocity along the flow direction. With the further increasing in void fraction, the bubbly flow turns into slug flow according to the flow pattern transient criteria of Orkiszewski which shows that the flow pattern is dependent on the void fraction. Generally speaking, the gas-drilling mud slug flow is composed of  liquid and gas slugs. In the preliminary slug flow, the liquid slugs are much longer than the gas slugs, so the wave velocity is determined primarily by the wave velocities in the liquid slugs which are hardly affected by the void fraction. At the position close to the wellhead, the pressure of two-phase flow fluid can be reduced to a low value which approaches the back pressure. The compressibility of gas will be improved at the low pressure. It results in significant increase in void fraction. As the void fraction increases greatly at the position close to the wellhead, gas slugs become much longer than the liquid slugs. It is assumed that the actual wave velocity in slug flow happens to play a leading role when the void fraction increases to some extent. For the decrease in liquid holdup in the gas slug with the increasing in void fraction, the wave velocity in the gas slug is infinitely close to that in pure gas, which is equivalent to an increase in wave velocity in the gas slug. Meanwhile, the gas content in liquid slug also increases and results in a lower wave velocity in the liquid slug. As a result, a slight increase in wave velocity appears.Figures7 and 8 present the influence of back pressure on the void fraction and wave velocity with respect to the parameters of well depth. As the back pressure increases, the void fraction at different depth of the annulus is reduced gradually, while the wave velocity in the two-phase flow tends to increase. Analytical results show that the increased back pressure is equivalent to be applied to the entire enclosed drilling fluid cyclical system. The pressure transmits from the wellhead to the bottomhole; therefore, the annular pressure in the entire wellbore is increased. According to the EOS, the density of gas increases and the compressibility of gas decreases with the increasing of gas pressure. So, the loss of interphase momentum and energy exchange is reduced and the interphase momentum exchange is promoted. It contributes to the increase in wave velocity with the increase in gas pressure. In addition, due to lower compressibility of two-phase flow medium under high pressure, the increase tendency of pressure wave velocity and the decrease tendency of void fraction are slowed down in the high back pressure range.Figure 7
Effect of back pressure on the void fraction.Figure 8
Effect of back pressure on the pressure wave velocity.


## 5.2. Effect of Gas Influx Rate on Wave Velocity

Figures9 and 10 graphically interpret the distributions of void fraction and variations of wave velocity along the flow direction in the annulus. When gas influx occurs in the bottomhole, gas invades into the wellbore and migrates from the bottomhole to the wellhead along the flow direction. At a low gas influx rate, it is extremely obvious that the void fraction and wave velocity first slightly change in a comparatively smooth value then change sharply. It is because of the rapid expansion of gas volume with the decreasing in pressure near the wellhead that the void fraction increases sharply, and the wave velocity decreases obviously at the same time. But under the high bottomhole pressure (up to 50 MPa), the compressibility of the gas is low. This results in a slight change in void fraction and wave velocity at the position far away from the wellhead. Since the compressible component increases with the increase in the gas influx rate, the compressibility of the gas and drilling mud two-phase fluid is improved. So the variations of void fraction and wave velocity become more prominent. Also, the void fraction still shows an increase tendency that is steady first and then sharp. It acts in accordance with the variation of void fraction at a low gas influx rate. At a low gas influx rate, if the void fraction at the position close to the wellhead can not increase to a high extent, the wave velocity always shows a decrease tendency. At a high gas influx rate, such as the Q=8.312 m3/h, the wave velocity tends to increase because the void fraction in the wellhead is increased to a high extent. In conclusion, the wave pressure is sensitive to the void fraction, and the void fraction is dominated by influx rate and pressure in the annulus, especially the influx rate.Figure 9
Void fraction distribution at different gas influx rates.Figure 10
Wave velocity variations at different gas influx rates.With the influx of gas, the mixing of gas and drilling mud occurs in the annulus and the corresponding interfacial transfer of momentum and mass causes an increase in gas phase void fraction and a decrease in pressure wave velocity, as shown in Figures11 and 12. Within the range of low gas influx rate, the wave velocity decreases significantly. It is because of this that the compressibility of the gas increases remarkably, and the medium appears to be of high elasticity, though the density of gas-drilling mud two-phase flow changes slightly. With the increase in the gas influx rate and corresponding increase in the void fraction in the annulus, the compressibility of the two-phase unceasingly increases, which promotes the momentum and energy exchange in the interface. So, the pressure wave continuously decreases. When the void fraction increases to some extent following the increase in the gas influx rate, the decrease in wave velocity in the liquid slug is gradually less than the increase in wave velocity in the gas slug; thus, the decrease of wave velocity is slowed down for the growth of gas slug. Especially in the wellhead, a slight increase in wave velocity is observed.Figure 11
Effect of gas influx rate on the void fraction.Figure 12
Effect of gas influx rate on the pressure wave velocity.


## 5.3. Effect of Virtual Mass Force on Wave Velocity

Figure13 illustrates the effect of virtual mass force on the wave velocity with the parameter of gas influx rate. As shown in Figure 15, under the same gas influx rate (0.594 m3/h, 1.098 m3/h, and 1.768 m3/h, resp.) the wave velocity curve of  Cvm=Re diverges from the curve of  Cvm=0 at the position close to the wellhead, whereas the wave velocity curve of the two types is almost coincided below the position of H=500 m. This divergence between the two types of curves is connected with the fact that significance of influence of interphase virtual mass forces increases together with the increase of the relative motion in the interphase. It is evident that the bottomhole pressure is hundreds of times higher than the wellhead pressure. As a result, the void fraction at the position close to the wellhead increases sharply, meanwhile the interphase momentum and energy exchanges are promoted. The virtual mass force can be described as the transfer of momentum between the gas phase and drilling mud phase caused by the relative motion in the interphase. If the relative motion in the interphase is quite weak, the value of virtual mass force will intend to approach 0, and the influence on the wave velocity can be ignored. However, if the relative motion is rather intense, the effect of virtual mass force on the wave velocity should not be ignored. Furthermore, taking the virtual mass force into account, the dispersion characteristic of the pressure wave weakens obviously. Compared with the pressure wave velocity calculated by ignoring the effect of virtual mass force, the calculated pressure wave velocity with a consideration of the virtual mass force is lower. Therefore, it is necessary to analyze the effect of the virtual mass force on the wave velocity at the position close to the wellhead in MPD operations.Figure 13
Effect of virtual mass force on the pressure wave velocity.


## 5.4. Effect of Disturbance Angular Frequency on Wave Velocity

Figure14 presents the pressure wave velocity in the annulus at different angular frequencies. According to Figures 5 and 9, the void fraction first slightly increases in the bottomhole and then sharply increases along the flow direction at the position close to the wellhead. At a fixed angular frequency, this phenomenon results in an overall decrease in the pressure wave velocity except for the section near the wellhead. At the section close to the wellhead, an opposite change trend of pressure wave velocity is observed for the transition of flow pattern from bubbly flow to slug flow due to the continuous increase in void fraction along the flow direction. In addition, it can be clearly seen from the curves that the wave velocity increases accompanying with the increase in the angular frequency above the position of H=500 m. This property is not very distinct at the position below H=500 m for the low void fraction.Figure 14
Wave velocity variations at different angular frequencies.Figure 15
Effect of angular frequency on wave velocity.Figure15 shows the effect of frequency on the wave velocity in the gas-drilling mud flow (Cvm=0). The curves of the wave velocity of different positions reveal that the propagation velocity of pressure disturbances increases together with the growth of the angular frequency (0<w<500 Hz). It proves that the pressure wave has an obvious dispersion characteristic in the two-phase flow. As the angular frequency increases in the range of less than 500 Hz, there is sufficient time to carry out the exchange of energy and momentum between two phases. It achieves a state of mechanical and thermodynamic equilibrium between the two phases. So the wave velocity increases gradually with the increase in the angular frequency at different depths of the annulus. It is considered that the wave velocity is mainly affected by the interphase mechanical and thermodynamic equilibrium at low anglular frequencies. When the angular frequency reaches the value of w=500 Hz, the pressure wave velocity achieves a constant value and remains on this level regardless of the further growth in angular frequency w. With the increase in angular frequency, there is not enough time for energy and momentum exchange between the gas-drilling mud two phases to reach the mechanical and thermodynamic equilibrium state, and thus the wave dispersion does not exist. At a high angular frequency, the wave velocity is mainly dominated by the mechanical and thermal nonequilibrium in the flow and keeps almost unchanged when the angular frequency is further increased. This is consistent with the influences of the angular frequency in the horizontal pipe [28]. Moreover, the effect of virtual mass force is also shown in this figure (Cvm=Re). It is observed that the wave velocity is significantly reduced when the virtual mass force is taken into consideration in the calculation of the wave velocity at the position close to the wellhead such as H=0 m. It can be explained by the reduction of dispersion characteristic of a pressure wave. Along the wellbore from the wellhead to the bottomhole, the distinction gradually decreases.


## 6. Conclusions

With full consideration of the important factors such as virtual mass force, drag force, gas void fraction, pressure, temperature, and angular frequency, a united wave velocity model has been proposed based on pressure drop gradient equations in MPD operations, gas-liquid two-fluid model, the gas-drilling mud equations of state, and small perturbation theory. Solved by the fourth-order explicit Runge-Kutta method, the model is used to predict wave velocity for different back pressures and gas influx rates in MPD operations. The main conclusions can be summarized as follows.(1)
With the introduction of virtual mass force and drag force, the united model agrees well with the previous experimental data. The united model can be used to accurately calculate the wave velocity in the annulus. The application of the model will be beneficial to further study the wave velocity at different gas influx rates and back pressures in MPD operations, reduce nonproductive times, and provide a reference for the drilling operations in extremely narrow pore/fracture windows existing conditions.(2)
The wave velocity and void fraction have different values with respect to well depth. In the annulus, the drop of pressure causes an increase in the void fraction along the flow direction. The void fraction increases first slightly and then sharply. Correspondingly, the wave velocity first gradually decreases in th bubbly flow and preliminary slug flow, and then the wave velocity slightly increases accompanying with the increase in the relative length ratio of gas slug to the liquid slug for the continuous increase in void fraction at the position close to the wellhead. The minimum wave velocity appears in the long gas slug flow.(3)
When the back pressure in the MPD operations increases, the void fraction at different depths of the annulus is reduced gradually, while the wave velocity in the two-phase flow tends to increase. Moreover, the influence of back pressure on the void fraction is greater at a position close to the wellhead than that in the bottomhole because of the great decrease in the pressure from the bottomhole to the wellhead. Also, the effect of back pressure on void fraction and wave velocity is decreased in the high back pressure range.(4)
The wave velocity is sensitive to the void fraction, but the void fraction is dominated by gas influx rate and pressure in the annulus, especially the gas influx rate. Since the compressibility of the gas and drilling mud two-phase fluid is improved with the increase in the influx rate, the void fraction increases greatly, and the wave velocity decreases significantly within the low gas influx rate range. When the void fraction is increased to some extent following the increase in the gas influx rate, the decrease of wave velocity is slowed for the growth of gas slug. Especially at the wellhead, a slight increase in wave velocity is observed at a high gas influx rate for the sharp increase in void fraction.(5)
It is necessary to analyze the effect of the virtual mass force on the wave velocity in MPD operations. Especially, at the position close to the wellhead the effect of virtual mass force is more prominent for the intense phase interactions. Taking the virtual mass force into account, the dispersion characteristic of the pressure wave weakens obviously. Compared with the results calculated by ignoring the effect of virtual mass force, the calculated pressure wave velocity with a consideration of the virtual mass force is lower.(6)
The effect of angular frequency on the propagation velocity of pressure wave appears at a low angular frequency. The propagation velocity of pressure disturbances increases together with the growth of the angular frequency (0<w<500 Hz). When the angular frequency reaches the value of w=500 Hz, the pressure wave velocity achieves a constant value and remains on this level regardless of the further growth in angular frequency w.



---

*Source: 318912-2013-09-04.xml*