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

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.

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 Figure 1, the gas and drilling mud two-phase fluid travels along the annulus in the drilling process. The fluid flows along the annulus in “−𝑧” 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)

The mass conservation equations are 𝜕𝜙𝜕𝑡^𝐺𝜌_𝐺+𝜕𝜙𝜕𝑧^𝐺𝜌_𝐺𝑣_𝐺𝜕=0,𝜙𝜕𝑡^𝐿𝜌_𝐿+𝜕𝜙𝜕𝑧^𝐿𝜌_𝐿𝑣_𝐿=0. (2)

The momentum conservation equations for gas are 𝜕𝜙𝜕𝑡^𝐺𝜌_𝐺𝑣_𝐺+𝜕𝜙𝜕𝑧^𝐺𝜌_𝐺𝑣²_𝐺+𝜕𝜙𝜕𝑧^𝐺𝜌_𝐺−𝜕𝜙𝜕𝑧^𝐺𝜏_𝐺^𝑓𝑟+𝜏_𝐺^Re+𝑀_Gi𝜏−4^𝐺𝐷=0. (3)

The momentum conservation equations for liquid are 𝜕𝜙𝜕𝑡^𝐿𝜌_𝐿𝑣_𝐿+𝜕𝜙𝜕𝑧^𝐿𝜌_𝐿𝑣²_𝐿+𝜕𝜙𝜕𝑧^𝐿𝜌_𝐿−𝜕𝜙𝜕𝑧^𝐿𝜏_𝐿^𝑓𝑟+𝜏_𝐿^Re+𝑀_Li𝜏−4^𝐿𝐷=0. (4)

The interphase forces include virtual mass force, drag force, and the wall shear stress. The transfer of momentum between the gas and drilling mud phases 𝑀_Gi and 𝑀_Li can be written as follows: 𝑀_Gi=−𝑀^nd_Li−𝑀^𝑑_Li+𝜏^𝑓𝑟_Li+𝜏^Re_Li𝜕𝜙_𝐿+𝜕𝜙𝜕𝑧_𝐺𝜎_𝑠+𝜕𝜙𝜕𝑧^𝐺𝑝_Gi𝜕𝑧−𝜙_𝐺𝜕𝑝_Li,𝑀𝜕𝑧_Li=𝑀^nd_Li+𝑀^𝑑_Li+𝑝_Li𝜕𝜙_𝐿−𝜏𝜕𝑧^𝑓𝑟Li+𝜏^Re_Li𝜕𝜙_𝐿.𝜕𝑧 (5)

The momentum transfer term caused by the gas-liquid interfacial relative acceleration motion (i.e., the virtual mass force) can be expressed as 𝑀^nd_Li=𝑐_VM𝜙_𝐺𝜌_𝐿𝛼_VM−0.1𝜙_𝐺𝜌_𝐿𝑣_𝑠𝜕𝑣_𝑠𝜕𝑧−𝑐_𝑚1𝜌_𝐿𝑣²_𝑠𝜕𝜙_𝐺,𝜕𝑧 (6) where 𝑐_𝑚1=0.1 and 𝑣_𝑠 is the slip velocity; 𝑣_𝑠=𝑣_𝐺−𝑣_𝐿.

The momentum transfer term caused by the drag force provided by Park et al. [47] can be described as 𝑀^𝑑_Li=38𝐶_𝐷𝑟𝜌_𝐿𝑅_𝐺𝑣²_𝑠. (7)

The pressure difference between the liquid interface and liquid is defined by the following formula presented by Park et al.: 𝑝_Li−𝑝_𝐿=−𝑐_𝑝𝜌_𝐿𝑣²_𝑠, (8) where 𝑐_𝑝=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 pressure 𝑃_Gi can be written as 𝑝_Gi−𝑝_𝐺≈0. (9)

According to the equations of Arnold [48], the pressure of the drilling mud is described as follows: 𝑝_𝐿=𝑝−0.25𝜌_𝐿𝜙_𝐺𝑣²_𝑠. (10)

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 𝜏_𝐺^𝑓𝑟≈𝜏^𝑓𝑟_Li≈𝜏_𝐿^𝑓𝑟≈𝜏_𝐺≈𝜏_𝐺^Re≈0. (11)

The Reynolds stress and interfacial average Reynolds stress can be obtained by 𝜏_𝐿^Re=−𝑐_𝑟𝜌_𝐿𝑣²_𝑠𝜙_𝐺𝜙_𝐿,𝜏^Re_Li=−𝑐_𝑟𝜌_𝐿𝑣²_𝑠, (12) where 𝑐_𝑟=0.2.

The wall shear stress of liquid phase is given by Wallis in the following form [30]: 𝜏_𝐿=0.5𝑓_𝐿𝜌_𝐿𝑣²_𝐿. (13)

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) where 𝑍_𝐺 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 (𝑝<35MPa) [49]: 𝑍_𝐺=1+0.3051−1.0467𝑇_𝑟−0.5783𝑇³_𝑟𝜌_𝑟+0.5353−02.6123𝑇_𝑟−0.6816𝑇³_𝑟𝜌²_𝑟, (15) where 𝑇_𝑟=𝑇/𝑇_𝑐, 𝑝_𝑟=𝑝/𝑝_𝑐, and 𝜌_𝑟=0.27𝑝_𝑟/𝑍_𝐺𝑇_𝑟.

The formula presented by Yarborough and Hall has been adopted to solve the gas deviation factor under the condition of high pressure (𝑝≥35MPa) [50]: 𝑍_𝐺=0.06125𝑃_𝑟𝑇_𝑟⁻¹exp−1.21−𝑇_𝑟⁻¹2𝑌, (16) where 𝑌 is given by −0.06125𝑝_𝑟𝑇_𝑟⁻¹exp−1.21−𝑇_𝑟⁻¹²+𝑌+𝑌²+𝑌³+𝑌⁴(1−𝑌)³=14.76𝑇_𝑟⁻¹−9.76𝑇_𝑟⁻²+4.58𝑇_𝑟⁻³𝑌²−90.7𝑇_𝑟⁻¹−242.2𝑇_𝑟⁻²+42.4𝑇_𝑟⁻³𝑌^{(2.18+2.82𝑇𝑟−1)}. (17)

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.

If  𝑇<130^∘C, the density of drilling mud can be obtained by the following equation: 𝜌_𝐿=𝜌₀1+4×10⁻¹⁰𝑝_𝐿−4×10⁻⁵𝑇−3×10⁻⁶𝑇². (18)

If 𝑇≥130^∘C, the density of drilling mud is 𝜌_𝐿=𝜌₀1+4×10⁻¹⁰𝑃_𝐿−4×10⁻⁵𝑇−3×10⁻⁶𝑇²+0.4𝑇−130𝑇². (19)

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]: 𝑇=𝑇_𝑒𝑖+𝐹1−𝑒^{(𝑧𝑏ℎ−𝑧)/𝐴}−𝑔sin𝜃𝑔_𝑐𝐽𝑐_𝑝𝑚+𝑁+𝑔_𝑇sin𝜃+𝑒^{(𝑧𝑏ℎ−𝑧)/𝐴}𝑇_𝑓𝑏ℎ−𝑇_𝑒𝑏ℎ. (20)

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)

For slug flow, 𝑞_𝐺𝑞_𝑚>𝐿_𝐵,𝑁_𝐺𝑉<𝐿_𝑆, (22) where 𝑞_𝑚 is the volumetric flow rate of two-phase flow, 𝑞_𝑚=𝑞_𝐺+𝑞_𝐿.

The dimensionless numbers 𝐿_𝑠 and 𝐿_𝐵 are defined as 𝐿_𝑠=50+36𝑁_𝐺𝑉𝑞_𝐿𝑞_𝐺,𝐿_𝐵=1.071−0.7277𝑣²_𝑚𝐷, (23) where 𝑁_𝐺𝑉=𝑣_𝑠𝜌_𝐿𝑔𝜎_𝑠^0.25. (24)

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 𝜙_𝐺+𝜙_𝐿=1. (25)

2.3.1. Bubbly Flow

As for bubbly flow, the density of gas and drilling mud mixture is described as the following formula: 𝜌_𝑚=𝜙_𝐿𝜌_𝐿+𝜙_𝐺𝜌_𝐺. (26)

The void fraction is determined by the following formula: 𝜙_𝐺=12𝑞1+^𝑚𝑣_𝑠𝐴−𝑞1+^𝑚𝑣_𝑠𝐴²−4𝑞_𝐺𝑣_𝑠𝐴. (27)

The coefficient of virtual mass force 𝐶_vm for bubbly flow can be expressed as follows [54]: 𝐶_vm=0.51+2𝜙_𝐺1−𝜙_𝐺. (28)

The coefficient of resistance coefficient 𝐶_𝐷 for bubbly flow can be expressed as 𝐶_𝐷=4𝑅_𝑏3𝑔𝜌_𝐿−𝜌_𝐺𝜎_𝑠1+17.671−𝜙_𝐺^9/718.671−𝜙_𝐺^1.52. (29)

The friction pressure gradient for bubbly flow can be obtained from the following equation: 𝜏_𝑓𝜌=𝑓^𝐿𝑣²_𝐿.2𝐷 (30)

2.3.2. Slug Flow

As for slug flow, distribution coefficient of gas in the liquid phase is 𝐶₀=0.00252lg10³𝜇_𝐿𝐴^1.38−0.782+0.232lg𝑣_𝑚−0.428lg𝐴. (31)

The average density of the mixture for slug flow is determined by 𝜌_𝑚=𝑊_𝑚+𝜌_𝐿𝑣_𝑠𝐴𝑞_𝑚+𝑣_𝑠𝐴+𝐶₀𝜌_𝐿. (32)

The void fraction for slug flow is determined by the following formula: 𝜙_𝐺=𝑞_𝐺𝑞_𝐺+𝑞_𝐿. (33)

The coefficient of virtual mass force 𝐶_vm for slug flow can be expressed as follows: 𝐶_vm=3.3+1.73𝐿_𝑞−3𝑅_𝑞3𝐿_𝑞−𝑅_𝑞. (34)

The coefficient of resistance coefficient 𝐶_𝐷 for slug flow can be expressed as 𝐶_𝐷=110𝜙³_𝐿𝑅_𝑏. (35)

The friction pressure gradient of bubbly flow can be obtained from the following equation: 𝜏_𝑓𝜌=𝑓^𝐿𝑣²_𝑚𝑞2𝐷^𝐿+𝑣_𝑠𝐴𝑞_𝑚+𝑣_𝑠𝐴+𝐶₀. (36)

2.4. Annulus Characteristic Analysis

The annulus effective diameter proposed by Sanchez [55] is used in the hydraulic calculation of annulus 𝜋𝐷𝐷=^2𝑜−𝐷²_𝑖/4𝜋𝐷_𝑜+𝐷_𝑖/4=𝐷_𝑜−𝐷_𝑖. (37)

The effective roughness of the annulus can be expressed as 𝑘_𝑒=𝑘_𝑜𝐷_𝑜𝐷_𝑜+𝐷_𝑖+𝑘_𝑖𝐷_𝑖𝐷_𝑜+𝐷_𝑖, (38) where the 𝐷_𝑖 and 𝐷_𝑜 are the diameters of inner pipe and outer pipe, respectively; the 𝑘₀ and 𝑘_𝑖 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 𝑑𝑝𝑑𝑧=𝜌_𝑚𝜏𝑔sin𝜃−^𝑤𝜋𝐷𝐴−𝜌_𝑚𝑣_𝑚𝑑𝑣_𝑚.𝑑𝑧 (39)

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)

Substituting (39) into (40), the total pressure drop gradient along the flow direction within the annulus can be expressed as 𝑑𝑝=𝜌𝑑𝑧^𝑚𝑔+𝜏_𝑓1−𝑊_𝑚𝑞_𝐺/𝐴²𝑝. (41)

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 𝑐_𝐺 and that of liquid phase 𝑐_𝐿 can be presented in the following form: 𝑑𝑝_𝐿𝑑𝜌_𝐿=𝑐²_𝐿,𝑑𝑝_𝐺𝑑𝜌_𝐺=𝑐²_𝐺. (42)

By introducing (42), the hydrodynamic equations of two-fluid model ((2)–(4)) can be written in the matrix form 𝐴𝜕𝜉𝜕𝑡+𝐵𝜕𝜉𝜕𝑧=𝐶𝜉, (43) where 𝐴 is the matrix of parameters considered in relation to time, 𝐵 is the matrix of parameters considered in relation to the spatial coordinate, and 𝐶 is the vector of extractions.

By introducing the small disturbance theory, the disturbance of the state variable 𝜉(𝜙_𝐺,𝑝,𝑣_𝐺,𝑣_𝐿)^𝑇 can be written as 𝜉=𝜉₀+𝛿𝜉exp[𝑖(𝑤𝑡−𝑘𝑡)]. (44)

In (44), 𝑘 is the wave number, and 𝑤 is the angular frequency of the disturbances. Substituting (44) into (43) gives homogenous linear equations concerning the expression (𝛿𝜙_𝐺,𝛿𝑝,𝛿𝑣_𝐺,𝛿𝑣_𝐿)^𝑇. 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: 𝑀₁𝑀₂𝑀₃𝑀₄−𝜌_𝐿𝑤1−𝜙_𝐺𝑐²_𝐿𝑤0−𝑘1−𝜙_𝐺𝜌_𝐿𝑀₅𝑀₆𝑀₇𝑀₈𝑀₉𝑀₁₀𝑀₈𝑀₁₁=0, (45) where 𝑀₁–𝑀₁₁ can be illustrated by 𝑀₁=𝜌_𝐺+𝑐_𝑝𝜙_𝐺𝜌_𝐿𝑣²_𝑠𝑐²_𝐺𝑀𝑤,₂=𝜙_𝐺𝑐²_𝐺1−𝑐_𝑝𝜙_𝐿𝑣²_𝑠𝑐²_𝐿𝑀𝑤,₃=−𝜙_𝐺𝜌_𝐺𝑘+2𝑐_𝑝𝜙_𝐺𝜙_𝐿𝜌_𝐿𝑣_𝑠𝑐²_𝐿𝑤,𝑀₄=22𝑐_𝑝𝜙_𝐺𝜙_𝐿𝜌_𝐿𝑣_𝑠𝑐²_𝐿𝑀𝑤,₅=𝜌_𝐿𝑣²_𝑟𝑘−𝜙_𝐺𝑐_𝑝+𝑐_𝑟−𝑐_𝑖+𝑐_𝑚2,𝑀₆=−𝜙_𝐺𝑘1−𝜙_𝐿𝑐_𝑝𝑣²_𝑠𝑐²_𝐿+𝑐_𝑖𝑣²_𝑠𝑐²_𝐿,𝑀₇=𝜙_𝐺𝜌_𝐺+𝑐_vm𝜌_𝐿3𝑤−𝑖4𝑐_𝐷𝑟𝜌_𝐿𝜙_𝐺𝑣_𝑠+4𝐷𝑓_𝐺𝑤𝜌_𝐺𝑣_𝐺,𝑀₈=−𝑐_vm𝜙_𝐺𝜌_𝐿3𝑤+𝑖4𝑐_𝐷𝑟𝜌_𝐿𝜙_𝐺𝑣_𝑠,𝑀₉=𝜌_𝐿𝑣²_𝑠𝑘𝜙_𝐿𝑐_𝑝−2𝑐_𝑟−𝑐_𝑚2,𝑀₁₀=−𝑘𝜙_𝐿+𝑐_𝑟𝜙_𝐺𝑣²_𝑠𝑐²_𝐿,𝑀₁₁=𝜌_𝐿𝜙_𝐿+𝜙_𝐺𝑐_vm3𝑤−𝑖4𝑐_𝐷𝑟𝜌_𝐿𝜙_𝐺𝑣_𝑠+4𝐷𝑓_𝐿𝜌_𝐿𝑣_𝐿, (46) where 𝑐_𝑚2=0.1, 𝑐_𝑝=0.25, 𝑐_𝑖=0.3, and 𝑐_𝑟=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 axis 𝑧, 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 𝑐. The wave velocity in the two-phase fluid can be determined by the following model: 𝑐=𝑤/𝑅⁺(𝑘)−(𝑤/𝑅⁻(𝑘))2. (47)

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)

With the initial value (𝑧₀,𝑝₀) and the function 𝐹(𝑧,𝑝), (50)–(53) can be obtained by 𝑘₁=𝐹𝑧₀,𝑝₀,𝑘₂=𝐹𝑧₀+ℎ2,𝑝₀+ℎ2𝑘₁𝑘₃,=𝐹𝑧₀+ℎ2,𝑝₀+ℎ2𝑘₂𝑘₄,=𝐹𝑧₀+ℎ,𝑝₀+ℎ𝑘₃, (49)(50)(51)(52) where ℎ is the step of depth. The pressure on the nod 𝑖=𝑖+1 can be obtained by 𝑝₁=𝑝₀+Δ𝑝=𝑝₀+ℎ6𝑘₁+2𝑘₂+2𝑘₃+𝑘₄. (53)

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 Figure 1. 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.

(a) Comparison with the experimental data of Liu

(b) Comparison with the experimental data of Li

Figure 3 

Experimental verification by comparison with previous experimental data.

The 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 Figure 4. 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 Table 1.

Table 1 

Parameters of calculation well.

The 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.

Figures 5 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.

Figures 7 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

Figures 9 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 𝑄=8.312 m³/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 Figures 11 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

Figure 13 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 m³/h, 1.098 m³/h, and 1.768 m³/h, resp.) the wave velocity curve of  𝐶_vm=Re diverges from the curve of  𝐶_vm=0 at the position close to the wellhead, whereas the wave velocity curve of the two types is almost coincided below the position of 𝐻=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

Figure 14 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 𝐻=500 m. This property is not very distinct at the position below 𝐻=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.

Figure 15 shows the effect of frequency on the wave velocity in the gas-drilling mud flow (𝐶_vm=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<𝑤<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 𝑤=500 Hz, the pressure wave velocity achieves a constant value and remains on this level regardless of the further growth in angular frequency 𝑤. 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 (𝐶_vm=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 𝐻=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<𝑤<500 Hz). When the angular frequency reaches the value of 𝑤=500 Hz, the pressure wave velocity achieves a constant value and remains on this level regardless of the further growth in angular frequency 𝑤.

Nomenclature

𝐴:	Annulus effective cross area (m2)
𝑐𝐺:	Wave velocity in gas (m/s)
𝑐𝐿:	Wave velocity in Drilling mud (m/s)
𝐶𝐷:	The coefficient of drag force
𝐷:	Annulus effective diameter (m)
𝐷𝑖:	Diameter of the inner pipe (m)
𝐷𝑜:	Diameter of the outer pipe (m)
𝑓𝐺:	Shear stresses coefficient of gas interface
𝑓𝐿:	Shear stresses coefficient of drilling mud interface
𝐹𝑧:	Pressure applied to the two-phase flow (N)
𝑔:	Acceleration due to gravity (m2/s)
𝑔𝑐:	Conversion factor
𝑔𝑇:	Geothermal temperature gradient (K/m)
ℎ:	One discrete length of annulus (m)
𝑘𝑒:	Annulus effective roughness (m)
𝑘𝑜:	Roughness of outer pipe (m)
𝑘𝑖:	Roughness of inner pipe (m)
𝐻:	The depth of annulus (m)
𝐿𝑞:	Air bubbles length of slug (m)
𝑀Gi:	Momentum transfer in gas interface (N/m3)
𝑀Li:	Momentum transfer in liquid interface (N/m3)
𝑀ndLi:	The none-drag-force (virtual mass force) (N/m3)
𝑀𝑑Li:	The momentum transfer term caused by the drag force (N/m3)
𝑁:	Parameter combining Thompson energy effects
𝑝Gi:	Gas interface pressure (MPa)
𝑝Li:	Liquid interface pressure (MPa)
𝑝:	Pressure (MPa)
𝑝0:	Initial pressure of annulus (MPa)
𝑝𝑐:	Critical pressure (kPa)
𝑝𝐿:	Pressure of drilling mud (MPa)
𝑝𝑟:	Reduced pressure
𝑞𝑚:	Volumetric flow rate of gas and drilling mud (m3/s)
𝑞𝐺:	Velocity of the gas (m3/s)
𝑞𝐿:	Velocity of the drilling mud (m3/s)
𝑄𝐿:	Drilling mud rate within the annulus (m3/h)
𝑟:	Average diameter of the bubble (m)
𝑅:	Constant of EOS (J/KgK)
𝑅𝑞:	Air bubbles radius of slug (m)
𝑅𝑏:	The bubble diameter (m)
𝑇:	Temperature (K)
𝑇𝑐:	Critical temperature (K)
𝑇𝑟:	Reduced temperature
𝑇𝑒𝑖:	Undisturbed formation temperature at a depth (K)
𝑇𝑒𝑏ℎ:	Undisturbed formation temperature of wellhead (K)
𝑇𝑓𝑏ℎ:	Undisturbed formation temperature at the bottomhole (K)
𝑣𝑚:	Gas and drilling mud flow velocity (m/s)
𝑣𝐺:	Gas flow velocity (m/s)
𝑣𝑠:	Slip velocity (m/s)
𝑣𝑠𝑔:	Superficial gas velocity (m/s)
𝑣𝐿:	Drilling mud flow velocity (m/s)
𝑊𝑚:	Mass flow velocity (kg/m3)
𝑍𝐺:	Gas deviation factor
𝑧0:	Initial calculation height of annulus (m)
𝑧𝑏ℎ:	Total well depth from surface (m)
𝑧:	Variable of distance from surface (m).

Greek Letters

𝜌𝐺:	Gas density (kg/m3)
𝜏𝑓𝑟Li:	Shear stresses of drilling mud noninterfacial (N/m2)
𝜏𝐿𝑓𝑟:	Shear stresses of drilling mud interface (N/m2)
𝜏𝐿Re:	Reynolds stress of drilling mud (N/m2)
𝜏𝐿Re:	Reynolds stress of drilling mud interface (N/m2)
𝜏𝐺Re:	Reynolds stress of gas interface (N/m2)
𝜏𝐺𝑓𝑟:	Shear stresses of gas interface (N/m2)
𝛼VM:	Acceleration due to virtual mass force coefficient (m/s2)
𝜏𝐿:	Shear stresses of drilling mud along wall (N/m2)
𝜏𝐺:	Shear stresses of gas along wall (N/m2)
𝜏𝐺Re:	Reynolds stress of gas (N/m2)
𝜏𝑓:	Frictional pressure gradient (Pa/m)
𝜏𝑤:	Frictional pressure coefficient
𝜌𝑚:	Gas and drilling density (kg/m3)
𝜌0:	Density under standard atmospheric pressure (kg/m3)
𝜌𝐿:	Drilling mud density (kg/m3)
𝜎𝑠:	Surface tension (N/m2)
Φ𝐺:	Gas void fraction
Φ𝐿:	Drilling mud holdup
𝜌𝑟:	Reduced density
𝜇𝐿:	Viscosity of drilling mud (Pa · s).

Subscripts

BHP:	Bottomhole pressure
BP:	Back pressure
EOS:	Equations of state
MPD:	Managed pressure drilling
RCD:	The rotating control device
R-K4:	The fourth-order explicit Runge-Kutta.

Subscripts of Graph

BP:	Back pressure (MPa)
𝑐:	Wave velocity in gas and drilling mud two-phase flow (m/s)
𝐶vm:	The coefficient of virtual mass force
𝐻:	Depth of annulus (m)
𝑄:	Gas influx rate at the bottomhole (m3/h)
Re:	The value of virtual mass force coefficient according to (28) or  (34)
𝑤:	Angle frequency (Hz)
Φ𝐺:	Void fraction (%).

Acknowledgment

Research work was cofinanced by the National Natural Science Foundation of China (nos. 51074135 and 51274170). Without their support, this work would not have been possible.