# A Novel Dynamic Model for Predicting Pressure Wave Velocity in Four-Phase Fluid Flowing along the Drilling Annulus

Authors: Xiangwei Kong; Yuanhua Lin; Yijie Qiu; Xing Qi
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

A dynamic pressure wave velocity model is presented based on momentum equation, mass-balance equation, equation of state, and small perturbation theory. Simultaneously, the drift model was used to analyze the flow characteristics of oil, gas, water, and drilling fluid multiphase flow. In addition, the dynamic model considers the gas dissolution, virtual mass force, drag force, and relative motion of the interphase as well. Finite difference and Newton-Raphson iterative are introduced to the numerical simulation of the dynamic model. The calculation results indicate that the wave velocity is more sensitive to the increase of gas influx rate than the increase of oil/water influx rate. Wave velocity decreases significantly with the increase of gas influx. Influenced by the pressure drop of four-phase fluid flowing along the annulus, wave velocity tends to increase with respect to well depth, contrary to the gradual reduction of gas void fraction at different depths with the increase of backpressure (BP). Analysis also found that the growth of angular frequency will lead to an increase of wave velocity at low range. Comparison with the calculation results without considering virtual mass force demonstrates that the calculated wave velocity is relatively bigger by using the presented model.


---

## Body

## 1. Introduction

In petroleum industry, managed pressure drilling (MPD) is considered to be one of the most important techniques, which allows accurate control of bottom hole pressure (BHP) by controlling the flow rate, drilling mud density, and back pressure (BP) at the wellhead [1]. As extensively used in the drilling of huge risk, uneconomical, or even abnormal formation, MPD has been of interest in the literature for at least the past decade, especially the related issues about pressure control inducing the dynamic well and kick detection [2]. Due to those efforts, research of the dynamic pressure wave velocity is of great significance to the detection of gas influx and effective control of the pressure at the bottom of well [3]. During the drilling operation of the so-called “microflux control,” an MPD technique developed by Santos et al. [4], the return flow is monitored and adjusted to control fluid loss or gain. In the light of control principal, simulation studies were performed to determine the most appropriate initial response to kicks arising due to MPD specific complications caused by BHP fluctuations [5]. During the MPD operations, all unsteady operating, such as changing of pumping rate, adjustment of choke, and controls of BP at wellhead, will generate a pressure wave and threaten the drilling equipment [6]. For the same reason, while tripping out of a drill string in the wellbore, bottom hole is submitted to a suddenly decrease in pressure, leading to fluid expansion and movement out of the annulus. The rapidly expanding fluids and dynamic pressure fluctuations can also lead to rock instability in a reservoir [7]. Particularly, the effects are more important in systems in which multiphase flows occur. New kick-detection tools are now available that are based on acoustic principles, which are of great benefit to potentially earlier and more sensitive detection of a gas influx than pit-gain or paddle flow measurements [8]. The study of propagation of pressure wave is also relevant to control of downhole tool, such as intelligent well downhole control valves applied in different field for many purposes. With further development of oilfield, downhole tool technique for special casing wells is receiving much more attention [9]. At the meantime, an important problem in deep drilling is the propagation of measurement-while-drilling (MWD) mud pulse, transmitting real-time various data from sensors located down hole near the drill bit. Hence, the propagation behavior of pressure wave is considered to provide reference for the MPD operations.However, we also noticed that the conventional theories present are difficult to be employed in a systematic and accurate prediction when influxes generates. As the influxes fluid including gas, oil, and water will lead to variations of physic characteristics parameter of fluid in the annulus, the distribution of pressure wave velocity in gas, oil, water, and drilling mud four-phase flow along the annulus will be dynamic changing with respect to time and well depth. These effects may influence the safe operation of devices.Pressure waves are disturbances that transmit energy and momentum along the wellbore through drilling fluid without significant displacement of matter. As migration of dispersed gas, water, and oil towards wellhead is quite complicated, the fundamental characteristics of the four-phase flow are still unknown and modeling results of pressure wave velocity are questionable. The proper treatment of propagation behavior of pressure wave in the four-phase flow in the annulus requires knowledge of description of influxes generation, development, and pressure wave propagation model in a two-phase mixture.In gas/liquid two-phase flow, two-phase medium interaction greatly changes the structure characteristic of flowing fluid, which results in greater compressibility of two-component flow than single-phase gas or liquid and further causes pressure waves propagating speed to be greatly reduced. 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 bottom hole, the density of the drilling mud has little variation while the compressibility increases obviously, which makes the wave velocity be decreased. This phenomenon was first proposed by Mallock [10] and attracts much attention for its important role in the development of science and technology applications. Extensive investigations involved in the issue of pressure wave propagation have been taken and some significant achievements associated with the theories have also been made in the researches. In early stage researching, Wood [11] extended the researches of Mallock and presented a succinct formula by assuming that the compressibility of two-component fluid is a function related to single-phase compressibility and elastic modulusE. Carstensen and Flody [12] measured the velocity of pressure wave under a lake. A dispersion relation for pressure waves propagating through a bubbly fluid was derived by using a linear scattering theory developed by Foldy. Under the hypothesis of homogeneous and adiabatic laminar flow, Thuraisingham [13] took the two-phase media as a homogeneous fluid and derived the solution model concerning the problem of wave velocity model in two-phase flow at low gas void fraction according to the analysis of state equation of mixture. According to the early stage investigations, Hsieh and Plesset [14] and Murray [15] researched the influence of thermal conductivity and viscosity coefficient on pressure wave propagation. Wallis [16] firstly studied the propagation mechanism of pressure wave and derived the propagation velocity in bubbly flow and separated flow using the homogeneous model. The proposed model is based on the respective compressibility of the vapor and the liquid. Also, the two-phase mixture is treated as a compressible fluid with suitably averaged properties. In that expression, mass and heat transfers are neglected, so it can be applied to any gas/liquid mixture, in so far as no major effect due to vaporization or condensation must be taken into account. It may not be valid for such complex mixtures that never reach equilibrium. Assuming that no evaporation or condensation occurs when pressure wave is transporting, Moody [17] developed a simple acoustic wave model for bubbly flow and annular flow and established a relationship between sonic velocity and two-phase critical flow. Similar model was developed by D’Arcy [18]. In the researches of Mcwilliam and Duggins [19], surface tension and compressibility of liquid phase were also considered. Henry et al. [20] calculated the velocity as a function of void fraction using a correlation to account for the change in bubble shape with void fraction. Martin and Padmanabhan [21] extended the simple model proposed by Henry by considering wave reflection and wave transmission at gas-liquid interfaces. Researches of Mori et al. [22, 23] suggest that impact pipe elasticity on pressure wave propagation velocity is limited to the range of gas void of less than 1%. At high gas void range, the pressure wave velocity is in between the two velocities of single-phase. Nguyen et al. [24] proposed another type of model for bubbly flow (diluted gas phase in the liquid). The simple relations for prediction of the propagation of pressure disturbances in liquid-gaseous two-phase systems are presented. The model makes use of the well-known physical behavior that the wave velocity in a single-phase fluid is influenced by the elasticity of the confining walls. The interface of the one phase is considered to act as the elastic wall of the other phase and vice versa. Mecredy and Hamilton [25] used the two-fluid model to predict the pressure wave propagation in vapor-liquid flow in detail. However, the analysis contained the important assumption that the evaporation or condensation was governed by kinetic theory. Michaelides and Zissis [26] developed a computational method which yields the sound velocity in terms of the thermodynamic coordinates of the substance without the use of diagrams. Corresponding velocities of sound for the four substances considered exhibit a certain similarity which is examined statistically. The relationship between the sound velocity and the critical mass flux is also investigated. Thuraisingham [13] studied the wave velocity in bubbly water at megahertz frequencies (1~10 MHZ). Numerical analytical results indicate that volume concentrations and the radius of the bubble relative to the incident wavelength of sound are the important parameters which determine the deviation of sound speed form that of bubble-free water.Currently, the two-fluid continuum model is the most common and reliable method to describe the gas/liquid two-phase flow phenomenon [27]. In the model, the governing equation and phase interface relationship is established based on the assumption that each phase satisfies the continuum conditions. To obtain the practical flow equations, reasonable assumptions and constitutive equations should be introduced. In consequence, the predicted wave velocities were found to depend strongly on the introduced assumptions and equations. In recent years, the two-fluid model was applied in determining the pressure wave propagation characteristics [28]. Ruggles et al. [29, 30] firstly performed the experimental investigation on the dispersion of pressure wave propagation in air-water bubbly flow and studied the propagation of pressure disturbance based on the two-fluid model small perturbation analysis method. Through the comparison with experimental data, they found that the virtual mass force coefficient is a function of gas void. Chung et al. [31] calculated the sonic velocity versus angular frequency form the concept of bubble compressibility in a two-component bubbly flow regime. He also extended such a model to predict the sonic velocity of a vapor-liquid system. Lee et al. [32] constructed the two fluids model to determine the pressure wave propagation speed for two-phase bubbly flow, slug flow, and stratified flow by using pressure disturbance instead of virtual mass and other phase interfacial terms. The results fit well with the experimental in steam water and air water of Henry and theoretical analysis of Nguyen. Zhao and Li [33] derived the general formula of sonic velocity in gas-liquid two-phase flow linear analysis using the linear analysis of the closed fundamental equations of compressible gas-liquid two-phase flow. It is proposed that the appropriate formula for calculating sonic velocity in gas-liquid two-phase flows under usual conditions may be Wood adiabatic sonic velocity formula. By linearizing the conservation equations of two-fluid model, Liu [34] derived a wave number equation of pressure wave for adiabatic gas-liquid two-phase flow. The effects of drag force and virtual mass force on propagation and dispersion of pressure wave were investigated. Xu and Chen [35] used the transient two-fluid model to develop a general relation for acoustic waves with steam-water two-phase mixture in one-dimensional flowing system. Both the mechanical and thermal nonequilibrium are considered. Brennen [36] taken mass and heat exchanges into account and proposed more complete expressions of the speed of sound in two-phase mixture. However, calibration of the mass and heat exchanges requires some further experimental investigations. Yeom and Chang [37] numerically investigated the wave propagation in the two-phase flows. An assessment was made on the effect of interfacial friction terms. Zhang et al. [27] investigated the propagation of the pressure wave in the water-gas two-phase bubbly flow with a one-dimensional two-fluid model and employing small perturbation analysis. The governing equations are simplified and closured according to high-speed aerated flow characteristics in hydraulic engineering. The effects of aerated concentration, liquid pressure, perturbation frequency, and interfacial forces on the acoustic wave velocity and its attenuation in the aerated flow are also explored. With the application of thermal phase change model in computational fluid dynamics code CFX, Li et al. [38] proposed a pressure wave propagation model and investigated the pressure wave propagation characteristics in two-phase fuel systems of liquid-propellant rocket. The propagation of pressure wave during the condensation of R404A and R134A refrigerants in pipe minichannels was given by Kuczyński [39]. Heat exchange between the phases in the condensation process was calculated by using the one-dimensional form of Fourier’s equation.In drilling industry, some scholars have been devoted to this aspect. In the late 1970s, the former Soviet All-Union Drilling Technology Research Institute [40] began to study characteristics of pressure wave velocity in gas-liquid two-phase flow to detect early gas influx and achieved some important results. To study relationship between pressure wave velocity and gas void fraction, Li et al. [41] launched a gas-drilling mud two-phase flow simulation experiment in vertical annulus. It is proved that the two phase flow of gas-liquid patterns and the velocity of gas migration can be determined if the well depth, mud properties, and void fraction in bottom are given The method which is faster in detection time than the method of conventional kick detection was proposed. Starting with the analysis of transient flow, combining with the theory of transmission line, Wang [42] obtained the calculating model of frequency domain for the pulse velocity in drilling fluid built and the impendence and transmission operator of drilling fluid. Alternative initial responses to kicks for various well scenarios during MPD operations were also explored by Davoudi et al. [43]. Li et al. [44] proposed a mathematical model for predicting the attenuation and propagation velocity of measurement while drilling (MWD) pressure pulses in aerated drilling using the two-phase flow model and considering the momentum and energy exchange at the phase interface, gravity of each phase, viscous pipe shear, and other closing conditions. According to the theory of unsteady flow, Xiushan [45] developed the formulas of transmission velocity for mud pulse signal. The formulas which cover all kinds of boundary conditions, including thin wall pipes and thick ones and interaction influence of gas content and solid content on transmission velocity are suitable for positive and negative mud pulse and accord with drilling practice. In a previous work, we proposed a united wave velocity model to predict the pressure wave velocity in gas-drilling mud two-phase steady flow. The effect of well depth, back pressure, gas influx rate, virtual mass force, and angular frequency are all considered. However, under the effects of buoyancy and complicated turbulence interaction, the existing theoretical solutions are not involved in the dynamic model for predicting pressure wave velocity in four-phase fluid flowing along the drilling annulus when influx fluid migrate towards the top of wellbore.In this paper, the drift model was used to analyze the flow characteristics of oil, gas, water, and drilling fluid multiphase flow. As the important characteristics of influx development, the relative motion of the interphase, such as slippage of gas phase and oil phase, is considered. Moreover, to predict pressure wave velocity in gas-oil-water-drilling mud four-phase flow in the annulus during MPD operations, a dynamic mathematical model is presented. By computing, the influence factors of pressure wave velocity, such as back pressure, gas void fraction, oil void fraction, influx time, influx rate, disturbance angular frequency, and virtual mass force, are analyzed.


## 2. Mathematical Model

Before introducing the new dynamic model and to make this point clear, this paper reviews the hydraulic system in MPD operations. The drilling system is a closed circulation with BP at wellhead. The key equipments include the rotating control device, dynamic well control system, conventional pressure control system, industrial personal computer, Coriolis meter, choke, and pressure sensor. First the drilling mud begins to circulate from mud tank, down the drill pipe, and the drill string and returns from the annulus travel back through mud pit where drilling solids are taken away and then to surface mud tank. An important function of the drilling fluid is to provide pressure support to the wellbore wall. The rock formation drilled through has some form of porosity filled with formation fluids. These fluids can be water, or in the case of a reservoir, hydrocarbons. The pressure in these fluids is referred to as the pore pressure. If the pore spaces are connected, these formations will also have permeability. Fluids can flow through them in response to a pressure gradient. The pressure in the annulus is controlled by varying BP to operate the fluid pressure in the wellbore. The aim in MPD is thus to maintain the pressure in the annulus between the two limits of pore pressure and fracture pressure [1].


### 2.1. The United Dynamic Model

When the bottom hole pressure is below the formation pressure, formation fluid will invade into the wellbore and the four-phase flow emerges in the annulus constituted by the drill string and wellbore. As seen in Figure1, take any cross section of the wellbore as an infinitesimal control volume. In the infinitesimal control volume, the four-phase drilling fluid is consisted by drilling mud (considered as a pseudohomogeneous liquid), influx oil (considered as oil phase), influx natural gas (considered as gas phase), and influx water (considered as water phase). Appropriate assumptions and governing equations are critical to simulate realistic four-phase well-control operations. The four-phase model was established based on the following assumptions:(1)
it is unsteady-state four-phase flow;(2)
the flow along the flow path is one-dimensional;(3)
the drilling mud is water-based;(4)
drilling mud is incompressible.Figure 1
Infinitesimal control volume in effective wellbore.In the analysis of the multiphase flow characteristic, oil, water, and drilling mud which are all considered as liquid phase have great differences in physical and chemical properties. As water-based medium, water and drilling mud have no substantive difference. In the flow processes, they blend quickly and have no clear phase boundary. Thus, the water-based fluid phase is discussed as water phase, and the oil is considered as another phase.As the presence of oil and gas interphase mass transfer in the multiphase flow system within the wellbore, the appropriate mass conservation equation can be listed according to oil, gas, and water three components. Given thatxik is the mass fraction of k components in i phase, the mass conservation equations for four-phase mixture are
(1)A∂∂tχokρoϕo+χgkρgϕg+χwiρwϕw+χmiρmϕm+∂∂sχokAρoϕovo+χgkAρgϕgvg+χwiAρwϕwvw+χmiAρmϕmvm=0.The momentum balance equation for the four-phase mixture is(2)∂∂tA∑kρkϕkvk+∂∂sA∑kρkϕkvk2+Ag∑kρkϕk+∂∂sAp+A∂p∂sfr=0.Equation (2) is a general momentum balance equation including hydrostatic pressure gradient, frictional pressure loss gradient, and acceleration loss gradient.Hence(3)∑iχoi=1,∑iχgi=1,∑iχwi=1,∑iχmi=1,
where
(4)χwo=0,χwg=0,χww=1,χwm=0,χmo=0,χmg=0,χmw=0,χmm=1,χgo=0,χgg=1,χgw=0,χgm=0,χoo=momo+mg,χog=mgmo+mg,χow=0,χom=0.As propagation velocity is greatly affected by the gas void fraction and angular frequency of the pressure disturbance, the superficial velocity of flowing medium has almost no effect on the propagation velocity [46]; oil, water, and drilling can be considered as liquid phase for their similarity in mechanics. According to the two-fluid model, the flow can be supposed to be gas-liquid two-phase flow from a macroscopic view.To establish the wave velocity dispersion equation, the mass conservation equations for liquid and gas two phases can be written individually as follows:(5)∂∂tϕgρg+∂∂sϕgρgvg=0,∂∂tϕLρL+∂∂sϕLρLvL=0.
Hence, ρL=ϕoρo+ϕwρw+ϕmρm.The gas momentum conservation equation is(6)∂∂sϕgρgvg+∂∂sϕgρgvg2=-∂∂sϕgρg+∂∂sϕgτgfr+τgRe+Mgi-4τgD.The liquid momentum conservation equation is(7)∂∂tϕLρLvL+∂∂sϕLρLvL2=-∂∂sϕLρL+∂∂sϕLτLfr+τLRe+MLi-4τLD.The transfer of momentumMgi and MLi can be written by the following equations:
(8)Mgi=-MLind-MLid+τLifr+τLiRe∂ϕL∂s+∂ϕgσs∂s+∂ϕgpgi∂s-ϕg∂pLi∂s,MLi=MLind+MLid+pLi∂ϕL∂s-τLifr+τLiRe∂ϕL∂s.Virtual mass force is obtained by the equation in the following form:(9)MLind=cVMϕgρLαVM-0.1ϕgρLvs∂vs∂s-cm1ρLvs2∂ϕg∂s,
where vs=vg-vL and cm1=0.1.The momentum transfer term is described as [47]
(10)MLid=38CDrρLRqvs2.The pressure difference between the liquid interface and liquid can be obtained by (11)(11)pLi-pL=-cpρLvs2,
where cp=0.25.The gas interface pressurepgi is defined as follows:
(12)pgi-pg≈0.The pressure of the liquid is(13)pL=p-0.25ρLϕgvs2.The shear stress and the interphase shear stress can be described as(14)τgfr≈τLifr≈τLfr≈τg≈τgRe≈0.The Reynolds stress and interfacial average Reynolds stress are(15)τLRe=-crρLvs2ϕgϕL,τLiRe=-crρLvs2.
Hence, cr=0.2.The wall shear stress of liquid phase is expressed as [48]
(16)τL=0.5fLρLvL2.The pressure wave velocity of gas phasecg and that of liquid phase cL can be expressed in the following form:
(17)dpLdρL=cL2,dpgdρg=cg2.Based on (17), the hydrodynamic equations of two-fluid model (5)–(7) can be written in the following matrix form:
(18)A∂ξ∂t+B∂ξ∂s=Cξ.
Here,A is the matrix of parameters considered in relation to time;B is the matrix of parameters considered in relation to the spatial coordinate;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
(19)ξ=ξ0+δξexp⁡iwt-kt,
where k is the wave number.According to the solvable condition of the homogenous linear equations that the determinant of the equations is zero, the equation of pressure wave can be expressed in the following form:(20)ρg+cpϕgρLvs2cg2wϕgcg21-cpϕLvs2cL2w-ϕgρgk+2cpϕgϕLρLvscL2w22cpϕgϕLρLvscL2w-ρLw1-ϕgcL2w0-k1-ϕgρLρLvr2k-ϕgcp+cr-ci+cm2-ϕgk1-ϕLcpvs2cL2+civs2cL2ϕgρg+cvmρLw-i34cDrρLϕgvs+4Dfgwρgvg-cvmϕgρLw+i34cDrρLϕgvsρLvs2kϕLcp-2cr-cm2-kϕL+crϕgvs2cL2-cvmϕgρLw+i34cDrρLϕgvsρLϕL+ϕgcvmw-i34cDrρLϕgvs+4DfLρLvL=0,where ci=0.3, cm2=0.1,cr=0.2, and  cp=0.25.The real value of wave number is determined as the pressure wave velocityc, and pressure wave velocity in the four-phase flow is
(21)c=w/R+k-w/R-k2.


### 2.2. Physical Models

To define the velocity of pressure wave propagation in the four-phase flow, the related physical models are required, such as equations of state, temperature distribution model, gas dissolution, and oil phase volume factor.


#### 2.2.1. Equations of State for Gas

The equation of state (EOS) for gas can be expressed as(22)ρg=pZg·R·T.Forp<35 MPa, the compression factor is obtained as follows:
(23)Zg=1+0.3051-1.0467Tr-0.5783Tr3ρr+0.5353-02.6123Tr-0.6816Tr3ρr2,
where Tr=T/Tc, pr=p/pc, ρr=0.27pr/ZGTr.Forp≥35 MPa, the compression factor under the condition of high pressure is [49](24)Zg=0.06125PrTr-1exp⁡-1.21-Tr-12Y,
whereY is given by the follow equations:
(25)-0.06125PrTr-1exp⁡-1.21-Tr-12+Y+Y2+Y3+Y41-Y3=14.76Tr-1-9.76Tr-2+4.58Tr-3Y2-90.7Tr-1-242.2Tr-2+42.4Tr-3Y2.18+2.82Tr-1.


#### 2.2.2. Equations of State for Liquid

WithT<130°C, the density of drilling mud is expressed as follows:
(26)ρm=ρ01+4×10-10pL-4×10-5T-3×10-6T2.
With T≥130°C, the density of drilling mud is expressed as follows:
(27)ρm=ρ0T-130T21+4×10-10pL-4×10-5T-3×10-6T2+0.4T-130T2.


#### 2.2.3. Temperature Distribution Model

The temperature of the drilling mud at any depth of the wellbore is [50]
(28)T=Tei+F1-e(zbh-z)/A-gsinθgcJcpm+N+gTsinθ+e(z_bh-z)/ATfbh-Tebh.


#### 2.2.4. Gas Dissolution

Assuming gas goes into and comes out of solution instantaneously, gas solubility can be obtained by (29)(29)Rs=0.021γgsp+0.1757101.7688/γos-0.001638T1.205.


#### 2.2.5. Oil Phase Volume Factor

The volume factor is calculated as follows:(30)Bo=0.976+0.000125.612γgsγos0.5Rs+2.25T+401.2.


### 2.3. Flow Pattern Prediction Models

The flow regime, the flow pattern, and structure of the flow are some of the important parameters to describe two-phase gas-quid flows, identify two phase gas-quid flow regimes, and calculate the dynamic pressure wave propagation velocity. Thus, the transition among the three main flow regimes (bubbly, slug, and annular) is desirable to be known [51].In the hydraulic calculation of annulus, the effective diameter of annulus [52] should be given as
(31)D=πDo2-Di2/4πDo+Di/4=Do-Di.The effective roughness of annulus can be calculated by(32)ke=k0DoDo+Di+kiDiDo+Di.At low gas flow rate, the liquid is continuous phase, and the gas bubbles are dispersed in the liquid phase. Studies of Taitel [53] give the minimum diameter necessary to form bubbly flow as
(33)Dmin⁡=19ρL-ρgσsgρL20.5.The critical condition for forming bubbly flow is(34)vMcr1.12=5.88D0.48gρL-ρgσs0.5σsρLρMμL0.08,(35)D>Dmin⁡,ϕg≤0.25,vM≤vMcr,ϕg≤0.52,vM>vMcr.For slug flow, the critical balance superficial flow rate of gas carrying droplets needs to meet the condition [54] that
(36)vsgcr=3.1gσρL-ρgρg20.25,(37)ϕg>0.25,vM≤vMcr,ϕg>0.52,vM>vMcr,vsg≤vsgcr.For annular flow, the pattern transition criterions [55] is
(38)vsg>vsgcr.


#### 2.3.1. Bubbly Flow

Gas void fraction of four-phase flow is described as(39)ϕg=vsgSgvso+vsg+vsw+vsm+vgr.The value of the distribution factorSg can be determined by
(40)Sg=1.20+0.371DiDo.Harmathy [56] established the calculation formula of gas slip velocity in bubbly flow based on the study of the migration velocity of the bubble in a stationary liquid as
(41)vgr=1.53gσsρL-ρgρL20.25.The average density of four-phase mixture flow is(42)ρM=ϕLρL+ϕgρg.The oil void fraction for four-phase flow is(43)ϕo=1-ϕgvsoSovso+vsw+vsm+1-ϕgvor.The value of the distribution factor isSo=1.05+0.371(Di/Do).On the basis of the total liquid fluid, establish the oil phase velocity relationship as(44)vo=SovL+vor.According to cross-section flow rate, phase distribution, and slip mechanism of liquid phase, we can the draw the following relationship:(45)vor=1.53gσwo-ρoρwb22,
where
(46)ρwb=ϕwρw+ϕmρm.Water void fraction is(47)ϕw=1-ϕg-ϕovswvsw+vsm.Drilling mud void fraction is(48)ϕm=1-ϕg-ϕo-ϕw.Due to the similar physical properties of water and drilling fluids,vw, vm, and vwb can be expressed by
(49)vw=vm=vwb.The coefficient of virtual mass forceCvm for bubbly flow can be expressed as follows:
(50)Cvm=0.51+2ϕg1-ϕg.The coefficient of resistance coefficientCD for bubbly flow can be expressed by
(51)CD=4Rb3gρL-ρgσs1+17.67ϕL9/718.67ϕL1.52.The friction pressure gradient for bubbly flow can be obtained from the following equation:(52)τf=fρLvL22D.


#### 2.3.2. Slug Flow

Thevoid fraction of the four phases,Φg, Φo, Φw, and Φm, can be determined by (39), (43), (47), and (48) the same as bubbly flow.The value of the distribution factorSg for slug flow can be described as
(53)Sg=1.182+0.9DiDo.For slug flow, the slip velocity can be calculated as follows.From experimental studies, Hasan and Kabir [57] established the calculation formula of drift velocity for slug flow on the basis of research on Taylor bubble migration rule of Davies and Taylor as
(54)vgr=0.35+0.1DiDogDoρL-ρgρL0.5.The coefficient of virtual mass forceCvm for slug flow can be expressed as follows:
(55)Cvm=3.3+1.73Lq-3Rq3Lq-Rq.The coefficient of resistance coefficientCD for slug flow can be expressed as
(56)CD=110ϕL3Rb.


#### 2.3.3. Annular Flow

As for annular flow, due to the miscible flow state of gas at center, the simplification can beVgr=0.The void fraction of gas can be determined by(57)ϕg=1+X0.8-0.378,
where X is defined as
(58)X=dp/dsLfrdp/dsgfr.Oil void fraction is(59)ϕo=1-ϕgvsovso+vsw+vsm.Water void fraction is(60)ϕw=1-ϕgvswvso+vsw+vsm,ϕm=1-ϕg-ϕo-ϕw.The same as slug flow,Cvm and CD can be determined by (55) and (56).


## 2.1. The United Dynamic Model

When the bottom hole pressure is below the formation pressure, formation fluid will invade into the wellbore and the four-phase flow emerges in the annulus constituted by the drill string and wellbore. As seen in Figure1, take any cross section of the wellbore as an infinitesimal control volume. In the infinitesimal control volume, the four-phase drilling fluid is consisted by drilling mud (considered as a pseudohomogeneous liquid), influx oil (considered as oil phase), influx natural gas (considered as gas phase), and influx water (considered as water phase). Appropriate assumptions and governing equations are critical to simulate realistic four-phase well-control operations. The four-phase model was established based on the following assumptions:(1)
it is unsteady-state four-phase flow;(2)
the flow along the flow path is one-dimensional;(3)
the drilling mud is water-based;(4)
drilling mud is incompressible.Figure 1
Infinitesimal control volume in effective wellbore.In the analysis of the multiphase flow characteristic, oil, water, and drilling mud which are all considered as liquid phase have great differences in physical and chemical properties. As water-based medium, water and drilling mud have no substantive difference. In the flow processes, they blend quickly and have no clear phase boundary. Thus, the water-based fluid phase is discussed as water phase, and the oil is considered as another phase.As the presence of oil and gas interphase mass transfer in the multiphase flow system within the wellbore, the appropriate mass conservation equation can be listed according to oil, gas, and water three components. Given thatxik is the mass fraction of k components in i phase, the mass conservation equations for four-phase mixture are
(1)A∂∂tχokρoϕo+χgkρgϕg+χwiρwϕw+χmiρmϕm+∂∂sχokAρoϕovo+χgkAρgϕgvg+χwiAρwϕwvw+χmiAρmϕmvm=0.The momentum balance equation for the four-phase mixture is(2)∂∂tA∑kρkϕkvk+∂∂sA∑kρkϕkvk2+Ag∑kρkϕk+∂∂sAp+A∂p∂sfr=0.Equation (2) is a general momentum balance equation including hydrostatic pressure gradient, frictional pressure loss gradient, and acceleration loss gradient.Hence(3)∑iχoi=1,∑iχgi=1,∑iχwi=1,∑iχmi=1,
where
(4)χwo=0,χwg=0,χww=1,χwm=0,χmo=0,χmg=0,χmw=0,χmm=1,χgo=0,χgg=1,χgw=0,χgm=0,χoo=momo+mg,χog=mgmo+mg,χow=0,χom=0.As propagation velocity is greatly affected by the gas void fraction and angular frequency of the pressure disturbance, the superficial velocity of flowing medium has almost no effect on the propagation velocity [46]; oil, water, and drilling can be considered as liquid phase for their similarity in mechanics. According to the two-fluid model, the flow can be supposed to be gas-liquid two-phase flow from a macroscopic view.To establish the wave velocity dispersion equation, the mass conservation equations for liquid and gas two phases can be written individually as follows:(5)∂∂tϕgρg+∂∂sϕgρgvg=0,∂∂tϕLρL+∂∂sϕLρLvL=0.
Hence, ρL=ϕoρo+ϕwρw+ϕmρm.The gas momentum conservation equation is(6)∂∂sϕgρgvg+∂∂sϕgρgvg2=-∂∂sϕgρg+∂∂sϕgτgfr+τgRe+Mgi-4τgD.The liquid momentum conservation equation is(7)∂∂tϕLρLvL+∂∂sϕLρLvL2=-∂∂sϕLρL+∂∂sϕLτLfr+τLRe+MLi-4τLD.The transfer of momentumMgi and MLi can be written by the following equations:
(8)Mgi=-MLind-MLid+τLifr+τLiRe∂ϕL∂s+∂ϕgσs∂s+∂ϕgpgi∂s-ϕg∂pLi∂s,MLi=MLind+MLid+pLi∂ϕL∂s-τLifr+τLiRe∂ϕL∂s.Virtual mass force is obtained by the equation in the following form:(9)MLind=cVMϕgρLαVM-0.1ϕgρLvs∂vs∂s-cm1ρLvs2∂ϕg∂s,
where vs=vg-vL and cm1=0.1.The momentum transfer term is described as [47]
(10)MLid=38CDrρLRqvs2.The pressure difference between the liquid interface and liquid can be obtained by (11)(11)pLi-pL=-cpρLvs2,
where cp=0.25.The gas interface pressurepgi is defined as follows:
(12)pgi-pg≈0.The pressure of the liquid is(13)pL=p-0.25ρLϕgvs2.The shear stress and the interphase shear stress can be described as(14)τgfr≈τLifr≈τLfr≈τg≈τgRe≈0.The Reynolds stress and interfacial average Reynolds stress are(15)τLRe=-crρLvs2ϕgϕL,τLiRe=-crρLvs2.
Hence, cr=0.2.The wall shear stress of liquid phase is expressed as [48]
(16)τL=0.5fLρLvL2.The pressure wave velocity of gas phasecg and that of liquid phase cL can be expressed in the following form:
(17)dpLdρL=cL2,dpgdρg=cg2.Based on (17), the hydrodynamic equations of two-fluid model (5)–(7) can be written in the following matrix form:
(18)A∂ξ∂t+B∂ξ∂s=Cξ.
Here,A is the matrix of parameters considered in relation to time;B is the matrix of parameters considered in relation to the spatial coordinate;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
(19)ξ=ξ0+δξexp⁡iwt-kt,
where k is the wave number.According to the solvable condition of the homogenous linear equations that the determinant of the equations is zero, the equation of pressure wave can be expressed in the following form:(20)ρg+cpϕgρLvs2cg2wϕgcg21-cpϕLvs2cL2w-ϕgρgk+2cpϕgϕLρLvscL2w22cpϕgϕLρLvscL2w-ρLw1-ϕgcL2w0-k1-ϕgρLρLvr2k-ϕgcp+cr-ci+cm2-ϕgk1-ϕLcpvs2cL2+civs2cL2ϕgρg+cvmρLw-i34cDrρLϕgvs+4Dfgwρgvg-cvmϕgρLw+i34cDrρLϕgvsρLvs2kϕLcp-2cr-cm2-kϕL+crϕgvs2cL2-cvmϕgρLw+i34cDrρLϕgvsρLϕL+ϕgcvmw-i34cDrρLϕgvs+4DfLρLvL=0,where ci=0.3, cm2=0.1,cr=0.2, and  cp=0.25.The real value of wave number is determined as the pressure wave velocityc, and pressure wave velocity in the four-phase flow is
(21)c=w/R+k-w/R-k2.


## 2.2. Physical Models

To define the velocity of pressure wave propagation in the four-phase flow, the related physical models are required, such as equations of state, temperature distribution model, gas dissolution, and oil phase volume factor.


### 2.2.1. Equations of State for Gas

The equation of state (EOS) for gas can be expressed as(22)ρg=pZg·R·T.Forp<35 MPa, the compression factor is obtained as follows:
(23)Zg=1+0.3051-1.0467Tr-0.5783Tr3ρr+0.5353-02.6123Tr-0.6816Tr3ρr2,
where Tr=T/Tc, pr=p/pc, ρr=0.27pr/ZGTr.Forp≥35 MPa, the compression factor under the condition of high pressure is [49](24)Zg=0.06125PrTr-1exp⁡-1.21-Tr-12Y,
whereY is given by the follow equations:
(25)-0.06125PrTr-1exp⁡-1.21-Tr-12+Y+Y2+Y3+Y41-Y3=14.76Tr-1-9.76Tr-2+4.58Tr-3Y2-90.7Tr-1-242.2Tr-2+42.4Tr-3Y2.18+2.82Tr-1.


### 2.2.2. Equations of State for Liquid

WithT<130°C, the density of drilling mud is expressed as follows:
(26)ρm=ρ01+4×10-10pL-4×10-5T-3×10-6T2.
With T≥130°C, the density of drilling mud is expressed as follows:
(27)ρm=ρ0T-130T21+4×10-10pL-4×10-5T-3×10-6T2+0.4T-130T2.


### 2.2.3. Temperature Distribution Model

The temperature of the drilling mud at any depth of the wellbore is [50]
(28)T=Tei+F1-e(zbh-z)/A-gsinθgcJcpm+N+gTsinθ+e(z_bh-z)/ATfbh-Tebh.


### 2.2.4. Gas Dissolution

Assuming gas goes into and comes out of solution instantaneously, gas solubility can be obtained by (29)(29)Rs=0.021γgsp+0.1757101.7688/γos-0.001638T1.205.


### 2.2.5. Oil Phase Volume Factor

The volume factor is calculated as follows:(30)Bo=0.976+0.000125.612γgsγos0.5Rs+2.25T+401.2.


## 2.2.1. Equations of State for Gas

The equation of state (EOS) for gas can be expressed as(22)ρg=pZg·R·T.Forp<35 MPa, the compression factor is obtained as follows:
(23)Zg=1+0.3051-1.0467Tr-0.5783Tr3ρr+0.5353-02.6123Tr-0.6816Tr3ρr2,
where Tr=T/Tc, pr=p/pc, ρr=0.27pr/ZGTr.Forp≥35 MPa, the compression factor under the condition of high pressure is [49](24)Zg=0.06125PrTr-1exp⁡-1.21-Tr-12Y,
whereY is given by the follow equations:
(25)-0.06125PrTr-1exp⁡-1.21-Tr-12+Y+Y2+Y3+Y41-Y3=14.76Tr-1-9.76Tr-2+4.58Tr-3Y2-90.7Tr-1-242.2Tr-2+42.4Tr-3Y2.18+2.82Tr-1.


## 2.2.2. Equations of State for Liquid

WithT<130°C, the density of drilling mud is expressed as follows:
(26)ρm=ρ01+4×10-10pL-4×10-5T-3×10-6T2.
With T≥130°C, the density of drilling mud is expressed as follows:
(27)ρm=ρ0T-130T21+4×10-10pL-4×10-5T-3×10-6T2+0.4T-130T2.


## 2.2.3. Temperature Distribution Model

The temperature of the drilling mud at any depth of the wellbore is [50]
(28)T=Tei+F1-e(zbh-z)/A-gsinθgcJcpm+N+gTsinθ+e(z_bh-z)/ATfbh-Tebh.


## 2.2.4. Gas Dissolution

Assuming gas goes into and comes out of solution instantaneously, gas solubility can be obtained by (29)(29)Rs=0.021γgsp+0.1757101.7688/γos-0.001638T1.205.


## 2.2.5. Oil Phase Volume Factor

The volume factor is calculated as follows:(30)Bo=0.976+0.000125.612γgsγos0.5Rs+2.25T+401.2.


## 2.3. Flow Pattern Prediction Models

The flow regime, the flow pattern, and structure of the flow are some of the important parameters to describe two-phase gas-quid flows, identify two phase gas-quid flow regimes, and calculate the dynamic pressure wave propagation velocity. Thus, the transition among the three main flow regimes (bubbly, slug, and annular) is desirable to be known [51].In the hydraulic calculation of annulus, the effective diameter of annulus [52] should be given as
(31)D=πDo2-Di2/4πDo+Di/4=Do-Di.The effective roughness of annulus can be calculated by(32)ke=k0DoDo+Di+kiDiDo+Di.At low gas flow rate, the liquid is continuous phase, and the gas bubbles are dispersed in the liquid phase. Studies of Taitel [53] give the minimum diameter necessary to form bubbly flow as
(33)Dmin⁡=19ρL-ρgσsgρL20.5.The critical condition for forming bubbly flow is(34)vMcr1.12=5.88D0.48gρL-ρgσs0.5σsρLρMμL0.08,(35)D>Dmin⁡,ϕg≤0.25,vM≤vMcr,ϕg≤0.52,vM>vMcr.For slug flow, the critical balance superficial flow rate of gas carrying droplets needs to meet the condition [54] that
(36)vsgcr=3.1gσρL-ρgρg20.25,(37)ϕg>0.25,vM≤vMcr,ϕg>0.52,vM>vMcr,vsg≤vsgcr.For annular flow, the pattern transition criterions [55] is
(38)vsg>vsgcr.


### 2.3.1. Bubbly Flow

Gas void fraction of four-phase flow is described as(39)ϕg=vsgSgvso+vsg+vsw+vsm+vgr.The value of the distribution factorSg can be determined by
(40)Sg=1.20+0.371DiDo.Harmathy [56] established the calculation formula of gas slip velocity in bubbly flow based on the study of the migration velocity of the bubble in a stationary liquid as
(41)vgr=1.53gσsρL-ρgρL20.25.The average density of four-phase mixture flow is(42)ρM=ϕLρL+ϕgρg.The oil void fraction for four-phase flow is(43)ϕo=1-ϕgvsoSovso+vsw+vsm+1-ϕgvor.The value of the distribution factor isSo=1.05+0.371(Di/Do).On the basis of the total liquid fluid, establish the oil phase velocity relationship as(44)vo=SovL+vor.According to cross-section flow rate, phase distribution, and slip mechanism of liquid phase, we can the draw the following relationship:(45)vor=1.53gσwo-ρoρwb22,
where
(46)ρwb=ϕwρw+ϕmρm.Water void fraction is(47)ϕw=1-ϕg-ϕovswvsw+vsm.Drilling mud void fraction is(48)ϕm=1-ϕg-ϕo-ϕw.Due to the similar physical properties of water and drilling fluids,vw, vm, and vwb can be expressed by
(49)vw=vm=vwb.The coefficient of virtual mass forceCvm for bubbly flow can be expressed as follows:
(50)Cvm=0.51+2ϕg1-ϕg.The coefficient of resistance coefficientCD for bubbly flow can be expressed by
(51)CD=4Rb3gρL-ρgσs1+17.67ϕL9/718.67ϕL1.52.The friction pressure gradient for bubbly flow can be obtained from the following equation:(52)τf=fρLvL22D.


### 2.3.2. Slug Flow

Thevoid fraction of the four phases,Φg, Φo, Φw, and Φm, can be determined by (39), (43), (47), and (48) the same as bubbly flow.The value of the distribution factorSg for slug flow can be described as
(53)Sg=1.182+0.9DiDo.For slug flow, the slip velocity can be calculated as follows.From experimental studies, Hasan and Kabir [57] established the calculation formula of drift velocity for slug flow on the basis of research on Taylor bubble migration rule of Davies and Taylor as
(54)vgr=0.35+0.1DiDogDoρL-ρgρL0.5.The coefficient of virtual mass forceCvm for slug flow can be expressed as follows:
(55)Cvm=3.3+1.73Lq-3Rq3Lq-Rq.The coefficient of resistance coefficientCD for slug flow can be expressed as
(56)CD=110ϕL3Rb.


### 2.3.3. Annular Flow

As for annular flow, due to the miscible flow state of gas at center, the simplification can beVgr=0.The void fraction of gas can be determined by(57)ϕg=1+X0.8-0.378,
where X is defined as
(58)X=dp/dsLfrdp/dsgfr.Oil void fraction is(59)ϕo=1-ϕgvsovso+vsw+vsm.Water void fraction is(60)ϕw=1-ϕgvswvso+vsw+vsm,ϕm=1-ϕg-ϕo-ϕw.The same as slug flow,Cvm and CD can be determined by (55) and (56).


## 2.3.1. Bubbly Flow

Gas void fraction of four-phase flow is described as(39)ϕg=vsgSgvso+vsg+vsw+vsm+vgr.The value of the distribution factorSg can be determined by
(40)Sg=1.20+0.371DiDo.Harmathy [56] established the calculation formula of gas slip velocity in bubbly flow based on the study of the migration velocity of the bubble in a stationary liquid as
(41)vgr=1.53gσsρL-ρgρL20.25.The average density of four-phase mixture flow is(42)ρM=ϕLρL+ϕgρg.The oil void fraction for four-phase flow is(43)ϕo=1-ϕgvsoSovso+vsw+vsm+1-ϕgvor.The value of the distribution factor isSo=1.05+0.371(Di/Do).On the basis of the total liquid fluid, establish the oil phase velocity relationship as(44)vo=SovL+vor.According to cross-section flow rate, phase distribution, and slip mechanism of liquid phase, we can the draw the following relationship:(45)vor=1.53gσwo-ρoρwb22,
where
(46)ρwb=ϕwρw+ϕmρm.Water void fraction is(47)ϕw=1-ϕg-ϕovswvsw+vsm.Drilling mud void fraction is(48)ϕm=1-ϕg-ϕo-ϕw.Due to the similar physical properties of water and drilling fluids,vw, vm, and vwb can be expressed by
(49)vw=vm=vwb.The coefficient of virtual mass forceCvm for bubbly flow can be expressed as follows:
(50)Cvm=0.51+2ϕg1-ϕg.The coefficient of resistance coefficientCD for bubbly flow can be expressed by
(51)CD=4Rb3gρL-ρgσs1+17.67ϕL9/718.67ϕL1.52.The friction pressure gradient for bubbly flow can be obtained from the following equation:(52)τf=fρLvL22D.


## 2.3.2. Slug Flow

Thevoid fraction of the four phases,Φg, Φo, Φw, and Φm, can be determined by (39), (43), (47), and (48) the same as bubbly flow.The value of the distribution factorSg for slug flow can be described as
(53)Sg=1.182+0.9DiDo.For slug flow, the slip velocity can be calculated as follows.From experimental studies, Hasan and Kabir [57] established the calculation formula of drift velocity for slug flow on the basis of research on Taylor bubble migration rule of Davies and Taylor as
(54)vgr=0.35+0.1DiDogDoρL-ρgρL0.5.The coefficient of virtual mass forceCvm for slug flow can be expressed as follows:
(55)Cvm=3.3+1.73Lq-3Rq3Lq-Rq.The coefficient of resistance coefficientCD for slug flow can be expressed as
(56)CD=110ϕL3Rb.


## 2.3.3. Annular Flow

As for annular flow, due to the miscible flow state of gas at center, the simplification can beVgr=0.The void fraction of gas can be determined by(57)ϕg=1+X0.8-0.378,
where X is defined as
(58)X=dp/dsLfrdp/dsgfr.Oil void fraction is(59)ϕo=1-ϕgvsovso+vsw+vsm.Water void fraction is(60)ϕw=1-ϕgvswvso+vsw+vsm,ϕm=1-ϕg-ϕo-ϕw.The same as slug flow,Cvm and CD can be determined by (55) and (56).


## 3. Solution of the Dynamic Model

Now since obtaining the analytic solution of the aforementioned theoretical model directly is impossible, discretization of the model to a numerical model is required [58]. In this paper, the mathematical methods based on finite difference method provide a numerical solution approach for the dynamic model. As for the solution of the pressure wave velocity model, spatial domain includes the entire wellbore and the formation node; time domain is the time period, influx fluid flowing from the bottom hole to the wellhead along the wellbore. Discretizing the domain of determinacy, the entire spatial and time domain can be divided into discrete networked systems.According to finite difference scheme, the four equations (5)–(7) are solved by using the finite difference method with computational cells shown in Figure 2; the difference equation systems of which described the basic principles of four-phase fluid motion in wellbore is presented as follows.Figure 2
Computational cells for semi-implicit difference solution.For the drilling mud phase,(61)Avsmi+1n+1-Avsmin+1Δs=Aϕmin+Aϕmi+1n-Aϕmin+1-Aϕmi+1n+12Δt.For the water phase,(62)Avswi+1n+1-Avswin+1Δs=Aϕwin+Aϕwi+1n-Aϕwin+1-Aϕwi+1n+12Δt.For the oil phase,(63)Avso/Boi+1n+1-Avso/Boin+1Δs=AϕsoBoin+AϕsoBoi+1n-AϕsoBoin+1-AϕsoBoi+1n+1×2Δt-1.For the gas phase,(64)Aρgvsg+ρgsRsvso/Boi+1n+1-Aρgvsg+ρgsRsvso/Boin+1Δs=Aρgϕg+ρgsRsϕo/Boin2Δt+Aρgϕg+ρgsRsϕo/Boi+1n2Δt-Aρgϕg+ρgsRsϕo/Boin+12Δt-Aρgϕg+ρgsRsϕo/Boi+1n+12Δt.For momentum balance equation,(65)Api+1n+1-Apin+1=ξ1+ξ2+ξ3+ξ4,
where(66)ξ1=Δs2ΔtAρmvsm+ρwvsw+ρovso+ρgvsgin+Aρmvsm+ρwvsw+ρovso+ρgvsgi+1n-Aρmvsm+ρwvsw+ρovso+ρgvsgin+1-Aρmvsm+ρwvsw+ρovso+ρgvsgi+1n+1,
and (67)ξ2=Aρmvsm2ϕm+ρwvsw2ϕw+ρgvsg2ϕg+ρovso2ϕoin+1-Aρmvsm2ϕm+ρwvsw2ϕw+ρgvsg2ϕg+ρovso2ϕoi+1n+1,ξ3=-gΔs2AρMin+1+AρMi+1n+1,ξ4=-Δs2A∂p∂sfrin+1+A∂p∂sfri+1n+1.Besides finite-difference method, the characterized method and the Newton-Raphson iterative method are adopted to solve the united dynamic model in four phases along annulus. On the basis of discretization, the solution of models was realized by applying personally complied code on VB. NET and the solution procedure of that are shown in Figure3.Figure 3
Solution procedures for pressure wave velocity in MPD operations.The four-phase flow system is described completely by eight variables including pressure, temperature, gas void fraction and liquid void fraction, gas and liquid densities, and gas and liquid velocities. There are still five unknowns, such as gas and liquid velocities, gas fraction, pressure, and gas density based on the above assumptions. Therefore, five equations are required to compute the unknown variables with boundary conditions. At different annulus depth, we can obtain pressure, temperature, gas velocity, oil velocity, water velocity, drilling mud velocity and gas void fraction, oil void fraction, and water void fraction during MPD operations by application of the finite-difference. At initial time, the wellhead BP, wellbore structure, well depths, wellhead temperature, gas-oil-water-drilling mud properties, and so forth are known. On nodei, flow parameters, such as pressure, temperature, and void fraction of gas-oil-water-drilling mud multiphase, can be obtained by adopting finite difference. Then, the determinant equation (20) is calculated based on the calculated parameters. Omitting the two unreasonable roots, the pressure wave at different depth of annulus in MPD operations can be solved by (21). Upon substitution of the actual magnitudes, v turned out to be the velocity of light. The process is repeated until the pressure wave velocity in the whole annulus has been obtained.The published experimental data presented by Li et al. [41] is used to verify the new developed united dynamic model. In Figure 4, the points represent Li’s experimental data and the lines represent the calculation results. From the comparisons, it can be concluded that the developed united dynamic model agrees well with the experimental data. The averaged error between model and data in Figure 4 is 2.853%.Figure 4
Comparison between calculation results of new dynamic model and the experimental data of Li.


## 4. Analysis and Discussion

One of the potentially most useful capabilities of transient well control model is the ability to simulate various choke control procedures. Typically, the control of choke to maintain constant BHP depends on the experience and training of the hand on the choke. The gas and liquid flow rate measured by Coriolis meter and the BP measured by pressure senor is taken as the initial data for annulus pressure calculation. The well used for calculation is a gas well in Sichuan Region, China. The wellbore structure of this well is shown in Figure5, and information about gas-oil-water-drilling mud properties, well design parameters, and operational conditions of calculation well are listed in Table 1.Table 1
Parameters of the calculation well.
Type
Property
Value
Mud
Density (kg/m3)
1460
Viscosity (Pa⋅s)
20
×
10
-
3
Gas
Relative density
0.65
Viscosity (Pa⋅s)
1.14
×
10
-
5
Oil
Density (kg/m3)
800
Viscosity (Pa⋅s)
40
×
10
-
3
Water
Density (kg/m3)
1000
Viscosity (Pa⋅s)
1.0
×
10
-
3
String
Elastic modulus of string (Pa)
2.07
×
10
11
Poisson ratio of string
0.3
Roughness (m)
1.54
×
10
-
7
Surface condition
Surface temperature (K)
298
Atmosphere pressure (MPa)
0.101Figure 5
Structural parameters of the calculation well.In the following example, the drilling mud mixed with gas, oil, and water is taken as a four-phase flow medium when the influx of gas, oil, and water occurs at the bottom of the well. In the calculation, the BP is 0.2 MPa, and flow rate isQg=3.6 m3/h, Qw=3.6 m3/h, Qo=3.6m3/h, and Qm=72.0m3/h, respectively. The propagation velocity of pressure wave in the four-phase flow is calculated and discussed by using the established model and well parameters.


### 4.1. Effect of Well Depth on the Pressure Wave Velocity

With bottom hole pressure being constant, the effect of well depth on pressure wave velocity is analyzed. When gas-oil-water influx fluid invades into the wellbore at the bottom hole, gas, oil, water, and drilling mud four-phase flow appears in the annulus and migrates from the bottom of the well (H=4000 m) to the wellhead along the flow direction. It is assumed that the initial length of the gas pillar is 500 m under the bottom hole pressure and temperature condition in the analysis of the migration of influx gas. The migration and change of gas column including the solubility and length change lead to the variation of pressure wave velocity in the annulus. The formation and development of influx are a dynamic process, and the flow of fluid in the annulus is very complicated. The flow condition and flow parameter at any time are various.Figure6 illustrates the distribution of pressure wave velocity at different well depth after 5.16 minutes and 18.71 minutes’ development of influx. It is shown that the top of aerated mud column arrives at the position of H=3000 m and H=2000 m, respectively. Due to the circulation of drilling mud and slip of gas phase, the length of aerated mud column that represents the rising height of gas will increase with the gradual decrease of pressure along the annulus. Besides, in a segment of aerated mud column, the gas void fraction of the top is higher than that of the bottom for the slop of gas phase. So, in the declining stage, pressure wave velocity in aerated mud column presents an increasing tendency along the flow direction. As the gas migrates upwards, the casing pressure should be changed dynamically to keep the bottom hole pressure constant, which will lead to an increase of annular pressure and back pressure. As the influx develops for 18.71 minutes, the back pressure will be prompted to 2.39 MPa, increased by 0.29 MPa compared with that of 5.16 minutes.Figure 6
Wave velocity variations at different depth.


### 4.2. Effect of BP on Pressure Wave Velocity

In Figures7 and 8, the three-dimensional figures illustrate the distributions of gas void fraction and variations of pressure wave velocity along the flow direction in the annulus at different time when BP is 0.1 MPa, 3.0 MPa, 6.0 MPa, and 9.0 MPa, respectively. When influx occurs, increase of applied BP is of the most commonly used approach for bottom hole pressure control. Control of back pressure has a great influence on the gas void fraction. In the two graphs, the two full lines give the shape of the void fraction and pressure wave velocity at different back pressure. In each graph, the two clusters of dotted lines at the two-dimensional plane are the projection of the full lines, which, respectively, represent the influence of well depth and time on the distributions of gas void fraction or variations of pressure wave velocity at different back pressure. Increasing the casing pressure will make the compressibility decrease obviously and energy dissipation with the propagation of pressure wave decrease as well. Thus, it can be seen that, with the reduction of back pressure, void fraction shows an increasing trend and that the pressure wave velocity has a decrease at different depth of well. Meanwhile, at the high temperature and high pressure bottom hole, pressure wave velocity change in a comparatively smooth level in a consequence of compressibility of the four-phase fluid changing only in a small degree. As a result, as aerated drilling mud migrating from the bottom hole to the wellhead, gas void fraction significantly increases first slightly and then sharply along the flow direction a drop of annular pressure. Conversely, the pressure wave velocity shows a gradual fist and then remarkable decreasing tendency after the influx gas reaches the position close to the wellhead. Also, with the development of influx, the influence of BP on the gas void fraction and pressure wave velocity is larger over time.Figure 7
Effect of BP on the gas void fraction.Figure 8
Effect of BP on the pressure wave velocity.


### 4.3. Effect of Gas Influx Rate on Pressure Wave Velocity

Figures9 and 10 graphically interpret the effect of gas influx rate on the distribution of gas void fraction and variations of pressure wave velocity along the flow direction in the annulus. The same as Figures 7 and 8, the two clusters of dotted lines at the two-dimensional plane are the projection of the full lines. The pressure wave velocity at different well depth is the transient value with respect to time. It is noted that the changes of pressure wave velocity are obvious with the migration of gas, which is extremely significant at wellhead even at low gas influx rate. With the decreasing of pressure at the position near the wellhead, rapid expansion of gas volume appears, and as a result the gas void fraction increases sharply and the pressure wave velocity decrease obviously at the same time. However, the gas void fraction and pressure wave velocity have no obvious change at the bottom of well due to the high pressure at the bottom of well and insignificant compressibility of gas. When the mainly reflect fundamental factors, initial gas influx rate, is increased, pressure wave velocity shows an obvious descending trend with the increase of initial gas influx rate, the value of pressure wave velocity at the position close to the wellhead is much greater than that at bottom hole. At the gas void fraction range of about 5%~80%, as gas migrating along the annulus, the pressure wave velocity is decreased to a minimum at the wellhead with the gas void fraction increasing to a maximum.Figure 9
Gas void fraction distribution at different gas influx rate.Figure 10
Wave velocity variations at different gas influx rate.


### 4.4. Effect of Oil Influx Rate on Pressure Wave Velocity

Figure11 shows the effect of oil influx on the gas void fraction in gas-oil-water-drilling mud four-phase flow. The increase of oil influx rate has little influence on the gas void fraction distribution, which is because gas void fraction is mainly impacted by pressure, temperature, and gas influx rate. When the oil influx rate increased from Qg=0.36m3/h to Qg=3.60m3/h, the change of mixture density is inconspicuous compared with drilling mud. Therefore, the gas void fraction keeps almost unchanged. Similarly, as shown in Figure 12, with the oil influx rate increase, the change of pressure wave velocity can be neglected due to the little influence of oil influx rate on the gas void fraction. Also, when the water influx invades the wellbore, the influx water will dissolve in drilling mud and reduce the density of drilling mud as the water phase and water base drilling mud have similar physical properties. But the influence of changing water influx rate on pressure wave velocity is too little to be neglected as the oil influx rate does.Figure 11
Gas void fraction distribution at different oil influx rates.Figure 12
Wave velocity variations at different oil influx rates.


### 4.5. Effect of Disturbance Angular Frequency on Pressure Wave Velocity

Figure13 provides an illustration of the influence of disturbance angular frequency on the pressure wave velocity when gas, oil, and water influx generates. It is indicated that an increase in angular frequency results in a rise of pressure wave velocity. In addition, the influence on dynamic pressure wave velocity is enhanced with the continuation of time and the decrease of well depth as the gas migrates from the bottom hole to the wellhead. Accordingly, the disturbance angular frequency has a significantly distinct influence on the pressure wave velocity at the position near the wellhead due to the phenomenon that sharply increases of gas void fraction. For an increase in angular frequency from 1500 Hz to 10000 Hz, the influence on wave velocity is not as great as the increase form 50 Hz to 1500 Hz. This is consistent with the fact that the effect of angular frequency on the pressure wave appears at a low angular frequency.Figure 13
Pressure wave velocity variations at different angular frequency.


## 4.1. Effect of Well Depth on the Pressure Wave Velocity

With bottom hole pressure being constant, the effect of well depth on pressure wave velocity is analyzed. When gas-oil-water influx fluid invades into the wellbore at the bottom hole, gas, oil, water, and drilling mud four-phase flow appears in the annulus and migrates from the bottom of the well (H=4000 m) to the wellhead along the flow direction. It is assumed that the initial length of the gas pillar is 500 m under the bottom hole pressure and temperature condition in the analysis of the migration of influx gas. The migration and change of gas column including the solubility and length change lead to the variation of pressure wave velocity in the annulus. The formation and development of influx are a dynamic process, and the flow of fluid in the annulus is very complicated. The flow condition and flow parameter at any time are various.Figure6 illustrates the distribution of pressure wave velocity at different well depth after 5.16 minutes and 18.71 minutes’ development of influx. It is shown that the top of aerated mud column arrives at the position of H=3000 m and H=2000 m, respectively. Due to the circulation of drilling mud and slip of gas phase, the length of aerated mud column that represents the rising height of gas will increase with the gradual decrease of pressure along the annulus. Besides, in a segment of aerated mud column, the gas void fraction of the top is higher than that of the bottom for the slop of gas phase. So, in the declining stage, pressure wave velocity in aerated mud column presents an increasing tendency along the flow direction. As the gas migrates upwards, the casing pressure should be changed dynamically to keep the bottom hole pressure constant, which will lead to an increase of annular pressure and back pressure. As the influx develops for 18.71 minutes, the back pressure will be prompted to 2.39 MPa, increased by 0.29 MPa compared with that of 5.16 minutes.Figure 6
Wave velocity variations at different depth.


## 4.2. Effect of BP on Pressure Wave Velocity

In Figures7 and 8, the three-dimensional figures illustrate the distributions of gas void fraction and variations of pressure wave velocity along the flow direction in the annulus at different time when BP is 0.1 MPa, 3.0 MPa, 6.0 MPa, and 9.0 MPa, respectively. When influx occurs, increase of applied BP is of the most commonly used approach for bottom hole pressure control. Control of back pressure has a great influence on the gas void fraction. In the two graphs, the two full lines give the shape of the void fraction and pressure wave velocity at different back pressure. In each graph, the two clusters of dotted lines at the two-dimensional plane are the projection of the full lines, which, respectively, represent the influence of well depth and time on the distributions of gas void fraction or variations of pressure wave velocity at different back pressure. Increasing the casing pressure will make the compressibility decrease obviously and energy dissipation with the propagation of pressure wave decrease as well. Thus, it can be seen that, with the reduction of back pressure, void fraction shows an increasing trend and that the pressure wave velocity has a decrease at different depth of well. Meanwhile, at the high temperature and high pressure bottom hole, pressure wave velocity change in a comparatively smooth level in a consequence of compressibility of the four-phase fluid changing only in a small degree. As a result, as aerated drilling mud migrating from the bottom hole to the wellhead, gas void fraction significantly increases first slightly and then sharply along the flow direction a drop of annular pressure. Conversely, the pressure wave velocity shows a gradual fist and then remarkable decreasing tendency after the influx gas reaches the position close to the wellhead. Also, with the development of influx, the influence of BP on the gas void fraction and pressure wave velocity is larger over time.Figure 7
Effect of BP on the gas void fraction.Figure 8
Effect of BP on the pressure wave velocity.


## 4.3. Effect of Gas Influx Rate on Pressure Wave Velocity

Figures9 and 10 graphically interpret the effect of gas influx rate on the distribution of gas void fraction and variations of pressure wave velocity along the flow direction in the annulus. The same as Figures 7 and 8, the two clusters of dotted lines at the two-dimensional plane are the projection of the full lines. The pressure wave velocity at different well depth is the transient value with respect to time. It is noted that the changes of pressure wave velocity are obvious with the migration of gas, which is extremely significant at wellhead even at low gas influx rate. With the decreasing of pressure at the position near the wellhead, rapid expansion of gas volume appears, and as a result the gas void fraction increases sharply and the pressure wave velocity decrease obviously at the same time. However, the gas void fraction and pressure wave velocity have no obvious change at the bottom of well due to the high pressure at the bottom of well and insignificant compressibility of gas. When the mainly reflect fundamental factors, initial gas influx rate, is increased, pressure wave velocity shows an obvious descending trend with the increase of initial gas influx rate, the value of pressure wave velocity at the position close to the wellhead is much greater than that at bottom hole. At the gas void fraction range of about 5%~80%, as gas migrating along the annulus, the pressure wave velocity is decreased to a minimum at the wellhead with the gas void fraction increasing to a maximum.Figure 9
Gas void fraction distribution at different gas influx rate.Figure 10
Wave velocity variations at different gas influx rate.


## 4.4. Effect of Oil Influx Rate on Pressure Wave Velocity

Figure11 shows the effect of oil influx on the gas void fraction in gas-oil-water-drilling mud four-phase flow. The increase of oil influx rate has little influence on the gas void fraction distribution, which is because gas void fraction is mainly impacted by pressure, temperature, and gas influx rate. When the oil influx rate increased from Qg=0.36m3/h to Qg=3.60m3/h, the change of mixture density is inconspicuous compared with drilling mud. Therefore, the gas void fraction keeps almost unchanged. Similarly, as shown in Figure 12, with the oil influx rate increase, the change of pressure wave velocity can be neglected due to the little influence of oil influx rate on the gas void fraction. Also, when the water influx invades the wellbore, the influx water will dissolve in drilling mud and reduce the density of drilling mud as the water phase and water base drilling mud have similar physical properties. But the influence of changing water influx rate on pressure wave velocity is too little to be neglected as the oil influx rate does.Figure 11
Gas void fraction distribution at different oil influx rates.Figure 12
Wave velocity variations at different oil influx rates.


## 4.5. Effect of Disturbance Angular Frequency on Pressure Wave Velocity

Figure13 provides an illustration of the influence of disturbance angular frequency on the pressure wave velocity when gas, oil, and water influx generates. It is indicated that an increase in angular frequency results in a rise of pressure wave velocity. In addition, the influence on dynamic pressure wave velocity is enhanced with the continuation of time and the decrease of well depth as the gas migrates from the bottom hole to the wellhead. Accordingly, the disturbance angular frequency has a significantly distinct influence on the pressure wave velocity at the position near the wellhead due to the phenomenon that sharply increases of gas void fraction. For an increase in angular frequency from 1500 Hz to 10000 Hz, the influence on wave velocity is not as great as the increase form 50 Hz to 1500 Hz. This is consistent with the fact that the effect of angular frequency on the pressure wave appears at a low angular frequency.Figure 13
Pressure wave velocity variations at different angular frequency.


## 5. Conclusions

Solution of dynamic pressure wave velocity in gas-oil-water-drilling mud four-phase flow is a lot more complicated than that in two-phase flow. The dynamic model which takes consideration of gas solubility, oil compression coefficients, perturbation frequency, virtual mass force, and drag force is established. Owing to the solution of the migration of drilling mud which contains influx fluid (gas, water, and oil), the dynamic pressure wave varied with time and well depth is solved, and the main conclusions can be summarized as follows.(
1
) Agreeing well with the previous experimental data, the developed untied dynamic model can be used to calculate the wave velocity in the annulus. The model can accurately give the dynamic at different depth and time during the migration of influx, which will be beneficial to provide a reference for well control and MWD.(
2
) In the dynamic migration process, the pressure wave velocities at different well depth are changing constantly due to the circulation of drilling mud and slip of gas phase. The length of aerated mud column which represents the rising height of gas increases with the gradually decreasing pressure along the annulus. In the segment of aerated mud column, pressure wave velocity presents an increasing tendency along the flow direction.(
3
) The flow condition and flow parameter at any time are various. Both the static pressure and void fraction distribution along the annulus changed with the position in the annulus and time and have direct impact on the pressure wave velocity. The pressure wave velocity happens to decrease after an elapsed period of time in which influx fluid migrates along the annulus. The decreasing tendency will last until the aerated mud column pass through that point of annulus completely. Also, the time required is gradually shortened during process of influx migration and growth.(
4
) At the bottom hole, the void of gas, water, and oil is very low relative to drilling mud. With the migration of influx fluid, the variation of gas void fraction is much greater than other phases; in particular, a sharp increase of gas void fraction presents at the position near the wellhead for the rapid expansion of gas volume, while the void fractions of water, oil, and drilling mud are decreased for the nearly negligible variations of physical properties and drastic expansion of gas at the position near the wellhead.(
5
) When the BP at the wellhead increases, gas void fraction at any position of annulus decreases; instead, the pressure wave velocity increases. Moreover, the influence of BP on the gas void fraction is larger at position close to the wellhead than that at the bottom of well. With the development of influx, the influence of BP on the gas void fraction and pressure wave velocity is larger over time.(
6
) The pressure wave velocity is sensitive to the gas influx rate, in particular at the position close to the wellhead. Even at low range, the rise of gas influx rate also will result in a significantly downward trend of pressure wave velocity. On the contrary, the influence of changing water influx rate and oil influx rate on pressure wave velocity is too little to neglect.(
7
) Ignoring virtual mass force, the pressure wave velocity has an obvious dispersion characteristic at different position of well. The effect of angular frequency on the pressure wave appears at a low angular frequency, and the pressure wave velocity increases together with the growth of the angular frequency.



---

*Source: 134102-2015-01-15.xml*