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

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


Introduction
In the oil and gas industry, the continuously increasing demand for underground resources has led to a shift in the focus on exploration and exploitation to deep or ultra-deep reservoirs. Due to the complicated geological conditions of these reservoirs, drilling technology aimed at creating a wellbore to reach the gas or oil reservoir encounters a multitude of challenges, particularly within narrow pore-fracture pressure windows [1,2]. Over the course of drilling, incidents occur frequently and cause various disadvantages, such as longer drilling times, higher costs, and vast risks. Conventional drilling is incapable of addressing these issues, and advanced technology-managed pressure drilling (MPD) is proposed as an alternative, which is useful for resolving complicated situations such as leakage, gas-kick, collapse, and trapping within a narrow density window through reducing the non-productive time and avoiding incidents [3,4]. Regarding offshore drilling, MPD shows remarkable advantages in coping with complex situations and has been widely used in drilling engineering, considering the continued exploration of challenging reservoirs [5,6]. Moreover, for the drilling and production of natural gas hydrate, MPD also creates incomparable advantages, including precise monitoring and control of pressure, avoidable formation damage and unwanted fractures, and rare gas kick and lost circulation [7,8].
Managing annular pressure in wellbores is considered a challenging task, as flow materials are likely to contain mud, oil, gas, and cuttings [9]. Typically, a gas-liquid two-phase flow is included in Figure 1. Relation between the models [11]. PDE: partial differential equation; DFM: drift flow model; and ODE: ordinary differential equation.
Concerning the DFM, accuracy and efficiency are crucial for the mathematical algorithm. The DFM is a typical hyperbolic equation that can be solved using the finite volume method. The Godunov and Roe schemes have been suggested to solve the DFM, but they depend on algebraic Riemann solvers or analytic Jacobian matrices of conservation equations [25,26]. A new pressurebased method was developed to address the convergence problem of the drift flux model by Wang [27]. Among the different schemes, it has been confirmed that the advection upwind splitting method (AUSM) schemes are effective solutions to the DFM in terms of their simplicity, accuracy, and robustness [28][29][30]. Not only do these schemes display accuracy in flux difference splitting (FDS) in the boundary layer, they also exhibit robustness in flux vector splitting (FVS) when capturing strong discontinuity. AUSM schemes are relatively late in performing the calculation of gas-liquid twophase flow, and nearly all of them were studied after the year 2000. Niu [31] applied the AUSM based on flux difference (AUSMD) scheme to solve the flow model of seven equations. Nevertheless, it was discovered that Niu's method was susceptible to numerical oscillation. Subsequently, Paillère et al. LOL Figure 1. Relation between the models [11]. PDE: partial differential equation; DFM: drift flow model; and ODE: ordinary differential equation.
Concerning the DFM, accuracy and efficiency are crucial for the mathematical algorithm. The DFM is a typical hyperbolic equation that can be solved using the finite volume method. The Godunov and Roe schemes have been suggested to solve the DFM, but they depend on algebraic Riemann solvers or analytic Jacobian matrices of conservation equations [25,26]. A new pressure-based method was developed to address the convergence problem of the drift flux model by Wang [27]. Among the different schemes, it has been confirmed that the advection upwind splitting method (AUSM) schemes are effective solutions to the DFM in terms of their simplicity, accuracy, and robustness [28][29][30]. Not only do these schemes display accuracy in flux difference splitting (FDS) in the boundary layer, they also exhibit robustness in flux vector splitting (FVS) when capturing strong discontinuity. AUSM schemes are relatively late in performing the calculation of gas-liquid two-phase flow, and nearly all of them were studied after the year 2000. Niu [31] applied the AUSM based on flux difference (AUSMD) scheme to solve the flow model of seven equations. Nevertheless, it was discovered that Niu's method was susceptible to numerical oscillation. Subsequently, Paillère et al. [32] achieved success in applying AUSM+ to solve the TFM of six equations. Evje et al. [28,33] also used the AUSM+ and AUSM based on flux vector (AUSMV) methods on a flow model of four equations and a drift flux model of three equations. Moreover, this method is applicable to engineering computing in drilling [34,35].
According to a general drift flux model of compressed gas-liquid two-phase flow, a model to describe the transient flow dynamics in the wellbore for the MPD is established in this paper. Pressure variations of the throttle valve at the wellhead and flow change resulting from the pressure difference at the bottom hole are considered as well. The governing equation of the DFM is numerically solved by the AUSM scheme with a second-order spatial accuracy and Van Leer slope limiter. The data obtained from published papers are used to validate the established model. In addition, a field application case is researched to support data analysis to make operation decisions.

Governing Equations
The typical hydraulic cycle used in the process of MPD is shown in Figure 2. The mud from the rig pump flows into the drill pipe through the surface equipment and riser, and flows through the bit nozzle at the bottom of the drill pipe. Then it enters the annular space and flows upward. After reaching the surface, it is processed by the solids control equipment and returns to the mud pit, then enters the rig pump again. The hydraulic model is the basis of the design and calculation of pressure control. MPD requires a real-time estimate based on the measured data, and it inputs the calculated results into the control system to achieve automatic control. It is costly to use pressure while drilling (PWD) tools to measure the downhole pressure, so the hydraulic model is of great significance. Additionally, even in the case of PWD, the hydraulic model is an essential auxiliary analytical tool.  [32] achieved success in applying AUSM+ to solve the TFM of six equations. Evje et al. [28,33] also used the AUSM+ and AUSM based on flux vector (AUSMV) methods on a flow model of four equations and a drift flux model of three equations. Moreover, this method is applicable to engineering computing in drilling [34,35]. According to a general drift flux model of compressed gas-liquid two-phase flow, a model to describe the transient flow dynamics in the wellbore for the MPD is established in this paper. Pressure variations of the throttle valve at the wellhead and flow change resulting from the pressure difference at the bottom hole are considered as well. The governing equation of the DFM is numerically solved by the AUSM scheme with a second-order spatial accuracy and Van Leer slope limiter. The data obtained from published papers are used to validate the established model. In addition, a field application case is researched to support data analysis to make operation decisions.

Governing Equations
The typical hydraulic cycle used in the process of MPD is shown in Figure 2. The mud from the rig pump flows into the drill pipe through the surface equipment and riser, and flows through the bit nozzle at the bottom of the drill pipe. Then it enters the annular space and flows upward. After reaching the surface, it is processed by the solids control equipment and returns to the mud pit, then enters the rig pump again. The hydraulic model is the basis of the design and calculation of pressure control. MPD requires a real-time estimate based on the measured data, and it inputs the calculated results into the control system to achieve automatic control. It is costly to use pressure while drilling (PWD) tools to measure the downhole pressure, so the hydraulic model is of great significance. Additionally, even in the case of PWD, the hydraulic model is an essential auxiliary analytical tool. As discussed in the previous section, the DFM is primarily a simplification of high-fidelity models using a complex conservation equation. In DFM, a mixture momentum equation is adopted rather than using an isolated equation for each phase. To achieve a high level of precision, some As discussed in the previous section, the DFM is primarily a simplification of high-fidelity models using a complex conservation equation. In DFM, a mixture momentum equation is adopted rather Energies 2019, 12, 3930 4 of 21 than using an isolated equation for each phase. To achieve a high level of precision, some closure relations are applied into the equations. The basic DFM includes two conservation equations for the gas and liquid and a mixture momentum equation [28]: Equation (1) can be written in vector form: where: where ρ, p, and υ are the density, pressure, and velocity. α is the volume fraction. The subscripts l and g denote the liquid and gas phases, respectively. F G represents the gravity term: F f is the friction term: where: f m can be calculated by [26]: where D and d denote the outer diameter and inner diameter of the annulus, respectively. Re is the Reynolds number; and µ m is the viscosity.

Closure Equations
For a constant mass flow, the following equations must be satisfied [20]: a is the speed of sound in the fluid. Neglect the last term: Similarly, for the momentum: Therefore, the density of the liquid is described as: Similarly, The following relationship is satisfied: In addition, the slip velocity equation proposed by Zuber and Findlay [16] is represented by where C 0 is the distribution parameter and υ t is the gas drift velocity.

Shi Slip Relation
In this research, a Shi slip relation was adopted, which can be found in [36]: where: 1.53 with a 1 = 0.2 and a 2 = 0.4. K u is the critical Kutateladze number: where C K u and C w denote the fitted parameters, which are set to 142 and 0.008, respectively. N B is the Bond number: where σ gl is the gas-liquid interfacial tension. V c is the characteristic velocity: The distribution parameter C 0 can be determined by: C max represents the value of the gas distribution parameter, which is usually set to 1.2. γ is calculated by: where: where F υ is a multiplier of the flooding velocity fraction that defaults to 1.0.
The above relation is based on the derivation of the vertical flow. In an inclined pipe, the terminal rise velocity can be indicated as follows: where θ is the angle of deviation and m is the drift velocity multiplier. m(θ) represents the influence of the wellborn inclination: For the water-gas flow, m 0 = 1.28, n 1 = 0.5, n 2 = 1.7 (θ < 88 • ). When 88 • ≤ θ ≤ 90 • , the relation has a large deviation for horizontal wells.

Primitive Variables
Equation (1) can be expressed as: where u 1 = α l ρ l , u 2 = α g ρ g and u 3 = α l ρ l v l + α g ρ g v g . Equation (21) can be expressed in terms of u 1 and u 2 : Substitute Equations (19) and (20) into Equation (35): It can be obtained by simple transformation: where: The pressure is then given by: The volume fractions are determined as: One can get the following equation by Equation (34): According to Equation (22): A system of linear equations from Equations (43) and (44) can be expressed in matrix form: Hence, we can get expressions for υ l and υ g :

Boundary Conditions
At the inlet, the convective flux is determined by the flow of the gas and liquid. The pressure flux can be obtained by extrapolation: It is necessary to analyze the two operating conditions of a shut-in and open well at the outlet boundary. During shut-in, the flow rates of gas and liquids at the wellhead are zero, and then the convective flux is zero. The pressure flux can be calculated as: While the well is open, the convective flux can be dealt with using the original variables at the outlet (velocity, density, volume fraction). The pressure flux is equal to the ground backpressure.
The two-phase flow in the choke valve can be represented as [20]: where K c is the valve discharge coefficient, z c is the choke valve opening index, and γ is the gas expansion factor. P choke and P atm represent the pressure of the choke and the atmosphere, respectively. The inflow rate resulting from the pressure difference between the reservoir and the bottom hole can be expressed as: P r and P wf are the reservoir pressure and bottom-hole pressure, respectively. q is the inflow rate. A and B are regression coefficients, which are mainly determined by the reservoir.

Advection Upwind Splitting Method Scheme
The flux at the cell interface is composed of the convective flux and pressure flux, which can be represented as [26,31]: F c l and F c g represent the liquid and gas phase convective flux, respectively, and F p represents the pressure flux.
For the convective terms: where Ψ is the speed function, l and g are liquid and gas phases, respectively, and L and R represent left and right cells, respectively. where: χ is the weight coefficient and can be expressed as: where v is the mixture velocity and c is the sound velocity.
ε is a small parameter to ensure smooth transition. The pressure flux is updated by:

Second-Order Accuracy
In this paper, we use a classical second-order monotonic upwind scheme for conservation laws (MUSCLs) technique to improve the accuracy: S is the Van Leer slope limiter:

Solution Method
When Equation (2) is discretized explicitly: Whereas Equation (2) is discretized implicitly: The Newton-Raphson method is used for solving the implicit equation, so that the process of the calculation is stable.

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

Solution Method
When Equation (2) is discretized explicitly: Whereas Equation (2) is discretized implicitly: The Newton-Raphson method is used for solving the implicit equation, so that the process of the calculation is stable.

Workflow
The primary calculation steps are shown in Figure 3

Laboratory Test
To complete the verification of the established model, we designed a simple experimental device, as shown in Figure 4. The device mainly includes a wellbore simulator, liquid circulation system, and air supply unit, as shown in Figures 5-8. The wellbore is a 6 m long transparent glass tube with a diameter of 180 mm, inside of which is an 80 mm blue plastic tube. The fluid is injected into the blue string, and the air is injected at the bottom. The air enters the annulus, and a gas-liquid two-phase flow is formed. Multiple pressure sensors are set in the wellbore to record the changes in pressure.
device, as shown in Figure 4. The device mainly includes a wellbore simulator, liquid circulation system, and air supply unit, as shown in Figures 5-8. The wellbore is a 6 m long transparent glass tube with a diameter of 180 mm, inside of which is an 80 mm blue plastic tube. The fluid is injected into the blue string, and the air is injected at the bottom. The air enters the annulus, and a gas-liquid two-phase flow is formed. Multiple pressure sensors are set in the wellbore to record the changes in pressure.   device, as shown in Figure 4. The device mainly includes a wellbore simulator, liquid circulation system, and air supply unit, as shown in Figures 5-8. The wellbore is a 6 m long transparent glass tube with a diameter of 180 mm, inside of which is an 80 mm blue plastic tube. The fluid is injected into the blue string, and the air is injected at the bottom. The air enters the annulus, and a gas-liquid two-phase flow is formed. Multiple pressure sensors are set in the wellbore to record the changes in pressure.   To complete the verification of the established model, we designed a simple experimental device, as shown in Figure 4. The device mainly includes a wellbore simulator, liquid circulation system, and air supply unit, as shown in Figures 5-8. The wellbore is a 6 m long transparent glass tube with a diameter of 180 mm, inside of which is an 80 mm blue plastic tube. The fluid is injected into the blue string, and the air is injected at the bottom. The air enters the annulus, and a gas-liquid two-phase flow is formed. Multiple pressure sensors are set in the wellbore to record the changes in pressure.   During the experiment (Figure 9), the annulus was filled with liquid at the beginning with a flow rate of 0.0028 m 3 /s and then injected with gas at 25 s, of which the flow rate was 0.01 kg/s. (32-60 s). Generally, the pressure in the wellbore increases first and then decreases, and finally, tends to be stable. During the experiment (Figure 9), the annulus was filled with liquid at the beginning with a flow rate of 0.0028 m 3 /s and then injected with gas at 25 s, of which the flow rate was 0.01 kg/s. Figure 8 shows the annulus flow dynamics. The flow state is mainly divided into three stages: the liquid phase flow stage (0-25 s), the transient gas-liquid two-phase flow stage (25-32 s), and the stable flow stage (32-60 s). Generally, the pressure in the wellbore increases first and then decreases, and finally, tends to be stable.   During the experiment (Figure 9), the annulus was filled with liquid at the beginning with a flow rate of 0.0028 m 3 /s and then injected with gas at 25 s, of which the flow rate was 0.01 kg/s. Figure 8 shows the annulus flow dynamics. The flow state is mainly divided into three stages: the liquid phase flow stage (0-25 s), the transient gas-liquid two-phase flow stage (25-32 s), and the stable flow stage (32-60 s). Generally, the pressure in the wellbore increases first and then decreases, and finally, tends to be stable.   During the experiment (Figure 9), the annulus was filled with liquid at the beginning with a flow rate of 0.0028 m 3 /s and then injected with gas at 25 s, of which the flow rate was 0.01 kg/s. Figure 8 shows the annulus flow dynamics. The flow state is mainly divided into three stages: the liquid phase flow stage (0-25 s), the transient gas-liquid two-phase flow stage (25-32 s), and the stable flow stage (32-60 s). Generally, the pressure in the wellbore increases first and then decreases, and finally, tends to be stable.

Full-Scale Experiment
Lage and Fjelde [37] conducted their injection experiments in a real well. The vertical test well was 1275 m deep with a parasitic pipe used for gas injection at a depth of 760 m and a pressure sensor fitted at 605 m. At the start, the wellbore was filled with water at a flow rate of 10.11 m 3 /s prior to the air being injected into the wellbore. The change in gas flow is presented in Figure 11. Over the course of gas injection, there was a significant change in the gas volume. The total time taken to conduct this experiment amounted to 1500 s, and the injection depth was relatively deep, which was used to exhibit the characteristics of the gas invasion process in practice.  Figure 12 presents a diagram that shows the changes in pressure over time at a latitude of 605 m. In this paper, the original measurement data as well as the calculation results obtained by Meng et al. [35] and Lage and Fjelde [37] are compared. As revealed by this figure, there is a decent consistency between the simulation results and the on-site testing data. The curve obtained in this study is smoother than that obtained by Lage and Fjelde. Meanwhile, the measurement data were contrasted with the calculation results, which indicates that the resulting error is within ±5%.
In order to further understand the flow dynamics in the wellbore, we analyzed the variation of other parameters according to the simulated data. Figure 13 reflects the changes in pressure with time and depth. As can be seen from the figure, in both the measured and simulated data, the pressure

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

Full-Scale Experiment
Lage and Fjelde [37] conducted their injection experiments in a real well. The vertical test well was 1275 m deep with a parasitic pipe used for gas injection at a depth of 760 m and a pressure sensor fitted at 605 m. At the start, the wellbore was filled with water at a flow rate of 10.11 m 3 /s prior to the air being injected into the wellbore. The change in gas flow is presented in Figure 11. Over the course of gas injection, there was a significant change in the gas volume. The total time taken to conduct this experiment amounted to 1500 s, and the injection depth was relatively deep, which was used to exhibit the characteristics of the gas invasion process in practice.  Figure 12 presents a diagram that shows the changes in pressure over time at a latitude of 605 m. In this paper, the original measurement data as well as the calculation results obtained by Meng et al. [35] and Lage and Fjelde [37] are compared. As revealed by this figure, there is a decent consistency between the simulation results and the on-site testing data. The curve obtained in this study is smoother than that obtained by Lage and Fjelde. Meanwhile, the measurement data were contrasted with the calculation results, which indicates that the resulting error is within ±5%.
In order to further understand the flow dynamics in the wellbore, we analyzed the variation of other parameters according to the simulated data. Figure 13 reflects the changes in pressure with time and depth. As can be seen from the figure, in both the measured and simulated data, the pressure  Figure 12 presents a diagram that shows the changes in pressure over time at a latitude of 605 m. In this paper, the original measurement data as well as the calculation results obtained by Meng et al. [35] and Lage and Fjelde [37] are compared. As revealed by this figure, there is a decent consistency between the simulation results and the on-site testing data. The curve obtained in this study is smoother than that obtained by Lage and Fjelde. Meanwhile, the measurement data were contrasted with the calculation results, which indicates that the resulting error is within ±5%.
In order to further understand the flow dynamics in the wellbore, we analyzed the variation of other parameters according to the simulated data. Figure 13 reflects the changes in pressure with time and depth. As can be seen from the figure, in both the measured and simulated data, the pressure shows a gradual decline, which is followed by a sharp fall, before an incremental increase again. The drop in pressure is due to the original fluid within the ambient environment being dispelled with the ingress of gas. It is known that gas has a lower density than fluid, causing the downhole pressure to decline on a continued basis. As the gas moves upwards, the pressure it is subjected to decreases, which causes the gas to expand further. As a result, more fluid is dispelled, and the drop in pressure occurs faster. Up to the point when the gas reaches the wellhead, the downhole pressure is at its minimum. With the continued ingress of gas, there will be an increase in downhole pressure to some extent until the gas-fluid two-phase flow stabilizes. Figure 14 shows the changes in the gas holdup. As revealed by this figure, there is a constant change in gas holdup at varying depths over time, which is consistent with the changes in the gas mass flow. This is most obvious at the wellhead. The gas holdup reaches its maximum at the wellhead when the maximum gas mass flow reaches the wellhead. shows a gradual decline, which is followed by a sharp fall, before an incremental increase again. The drop in pressure is due to the original fluid within the ambient environment being dispelled with the ingress of gas. It is known that gas has a lower density than fluid, causing the downhole pressure to decline on a continued basis. As the gas moves upwards, the pressure it is subjected to decreases, which causes the gas to expand further. As a result, more fluid is dispelled, and the drop in pressure occurs faster. Up to the point when the gas reaches the wellhead, the downhole pressure is at its minimum. With the continued ingress of gas, there will be an increase in downhole pressure to some extent until the gas-fluid two-phase flow stabilizes. Figure 14 shows the changes in the gas holdup. As revealed by this figure, there is a constant change in gas holdup at varying depths over time, which is consistent with the changes in the gas mass flow. This is most obvious at the wellhead. The gas holdup reaches its maximum at the wellhead when the maximum gas mass flow reaches the wellhead. shows a gradual decline, which is followed by a sharp fall, before an incremental increase again. The drop in pressure is due to the original fluid within the ambient environment being dispelled with the ingress of gas. It is known that gas has a lower density than fluid, causing the downhole pressure to decline on a continued basis. As the gas moves upwards, the pressure it is subjected to decreases, which causes the gas to expand further. As a result, more fluid is dispelled, and the drop in pressure occurs faster. Up to the point when the gas reaches the wellhead, the downhole pressure is at its minimum. With the continued ingress of gas, there will be an increase in downhole pressure to some extent until the gas-fluid two-phase flow stabilizes. Figure 14 shows the changes in the gas holdup. As revealed by this figure, there is a constant change in gas holdup at varying depths over time, which is consistent with the changes in the gas mass flow. This is most obvious at the wellhead. The gas holdup reaches its maximum at the wellhead when the maximum gas mass flow reaches the wellhead.  Figure 15 presents an analysis of the gas mass flow. As revealed in the figure, there is a significant change in gas mass flow, with the maximum at 0.265 kg/s. The variation in gas mass flow accounts for the changes to all other characterizing parameters for the wellbore. In practice, due to the difference between the geological pressure and that at the downhole, there will be a constant change in the gas mass flow, as a result of which the characteristics exhibited by the transient multiphase flow will change inside the wellbore. Figure 16 reflects the changes in fluid mass flow. The maximum fluid mass at the wellhead is 63 kg/s, and the minimum is 4 kg/s, which is because as the amount of filled gas is on the rise and the fluid inside the wellbore is dispelled, the fluid mass flow is significant when there is less gas ingress. When the gas mass increases, the fluid mass flow will drop noticeably.  Figure 15 presents an analysis of the gas mass flow. As revealed in the figure, there is a significant change in gas mass flow, with the maximum at 0.265 kg/s. The variation in gas mass flow accounts for the changes to all other characterizing parameters for the wellbore. In practice, due to the difference between the geological pressure and that at the downhole, there will be a constant change in the gas mass flow, as a result of which the characteristics exhibited by the transient multiphase flow will change inside the wellbore. Figure 16 reflects the changes in fluid mass flow. The maximum fluid mass at the wellhead is 63 kg/s, and the minimum is 4 kg/s, which is because as the amount of filled gas is on the rise and the fluid inside the wellbore is dispelled, the fluid mass flow is significant when there is less gas ingress. When the gas mass increases, the fluid mass flow will drop noticeably.  Figure 15 presents an analysis of the gas mass flow. As revealed in the figure, there is a significant change in gas mass flow, with the maximum at 0.265 kg/s. The variation in gas mass flow accounts for the changes to all other characterizing parameters for the wellbore. In practice, due to the difference between the geological pressure and that at the downhole, there will be a constant change in the gas mass flow, as a result of which the characteristics exhibited by the transient multiphase flow will change inside the wellbore. Figure 16 reflects the changes in fluid mass flow. The maximum fluid mass at the wellhead is 63 kg/s, and the minimum is 4 kg/s, which is because as the amount of filled gas is on the rise and the fluid inside the wellbore is dispelled, the fluid mass flow is significant when there is less gas ingress. When the gas mass increases, the fluid mass flow will drop noticeably.

Sensitivity Analysis
In the previous sections, the model proposed in this paper was validated and shown to be suited to industrial application in the oil and gas industry. One example is a well located in the north of Sichuan Province, where the MPD technique is applied. To complete the calculation, it was supplemented with the following parameters (Tables 1 and 2). For more effective control of the downhole pressure during drilling, the application of the existing technology to improve the effectiveness in downhole pressure control is discussed, and an analysis is conducted of the impact of different parameters on flow inside the wellbore.

Sensitivity Analysis
In the previous sections, the model proposed in this paper was validated and shown to be suited to industrial application in the oil and gas industry. One example is a well located in the north of Sichuan Province, where the MPD technique is applied. To complete the calculation, it was supplemented with the following parameters (Tables 1 and 2). For more effective control of the downhole pressure during drilling, the application of the existing technology to improve the effectiveness in downhole pressure control is discussed, and an analysis is conducted of the impact of different parameters on flow inside the wellbore. The wellhead backpressure has a direct effect on the downhole pressure. As shown in Figure 17, lower wellhead backpressure leads to a faster decline in the downhole pressure when overflow arises. When the wellhead backpressure is 0.1 MPa, the downhole pressure declines by 7.5 MPa, which means an increment of 1.9 MPa as compared to the time point when the wellhead backpressure is 2.0 MPa. Therefore, it can be seen that applying pressure to some extent at the wellhead is effective for bringing the reduction in the downhole pressure under control.
The wellhead backpressure has a direct effect on the downhole pressure. As shown in Figure 17, lower wellhead backpressure leads to a faster decline in the downhole pressure when overflow arises. When the wellhead backpressure is 0.1 MPa, the downhole pressure declines by 7.5 MPa, which means an increment of 1.9 MPa as compared to the time point when the wellhead backpressure is 2.0 MPa. Therefore, it can be seen that applying pressure to some extent at the wellhead is effective for bringing the reduction in the downhole pressure under control. Even a slight decline in the density of drilling fluid will cause a sharp drop in the downhole pressure. As shown in Figure 18, the downhole pressure difference is about 3.4 MPa when the difference in drilling fluid density is 0.1 kg/m 3 . A commonly-seen practice for controlling the downhole pressure is to change the density of the drilling fluid. However, as it takes time for the drilling fluid to reach the downhole from the wellhead, the downhole pressure is incapable of being adjusted promptly, which means that the possibility of various dangerous incidents occurring remains possible. Therefore, a combined approach taken while drilling with backpressure and displacement would be more effective and reasonable to prevent the occurrence of incidents.
Over the course of drilling, displacement has a substantial impact on the downhole pressure. As indicated by Figure 19, a lower displacement leads to a lower downhole pressure, which will decline at a faster pace over time. Therefore, as far as the drilling facilities allow, an increase in the displacement will cause a significant rise in the downhole pressure.
In actual conditions, the pressure difference between the reservoir and the downhole will constantly change, and the invasion amount of gas will also vary synchronously. As shown in Figure  20, a higher initial pressure difference causes a more significant variation in the downhole pressure. Even a slight decline in the density of drilling fluid will cause a sharp drop in the downhole pressure. As shown in Figure 18, the downhole pressure difference is about 3.4 MPa when the difference in drilling fluid density is 0.1 kg/m 3 . A commonly-seen practice for controlling the downhole pressure is to change the density of the drilling fluid. However, as it takes time for the drilling fluid to reach the downhole from the wellhead, the downhole pressure is incapable of being adjusted promptly, which means that the possibility of various dangerous incidents occurring remains possible. Therefore, a combined approach taken while drilling with backpressure and displacement would be more effective and reasonable to prevent the occurrence of incidents.
Over the course of drilling, displacement has a substantial impact on the downhole pressure. As indicated by Figure 19, a lower displacement leads to a lower downhole pressure, which will decline at a faster pace over time. Therefore, as far as the drilling facilities allow, an increase in the displacement will cause a significant rise in the downhole pressure.
In actual conditions, the pressure difference between the reservoir and the downhole will constantly change, and the invasion amount of gas will also vary synchronously. As shown in Figure 20, a higher initial pressure difference causes a more significant variation in the downhole pressure. Transient gas-liquid two-phase flow is always present in the wellbore, and there is no stable stage. The downhole pressure drop can be divided into two stages: the accelerated stage while gas is in the annulus; and the steady descent stage after the gas reaches the wellhead. Transient gas-liquid two-phase flow is always present in the wellbore, and there is no stable stage. The downhole pressure drop can be divided into two stages: the accelerated stage while gas is in the annulus; and the steady descent stage after the gas reaches the wellhead.

Conclusions
In this paper, an MPD wellbore transient flow dynamic model was developed based on the compressible gas-fluid two-phase flow with the drift flux model. The most critical findings are summarized below: (1) In laboratory experiments, the flow state is mainly divided into three stages: the liquid phase flow stage (0-25 s), the transient gas-liquid two-phase flow stage (25-32 s), and the stable flow stage (32-60 s). The simulated data are in good agreement with the experimental results, and the error range is within ±10%.
(2) The pressure shows a gradual decline, which is followed by a sharp fall, before an incremental increase again. The drop in pressure is due to the original fluid being dispelled with the ingress of gas. As the gas moves upwards, the pressure it is subjected to decreases, which causes the gas to expand further. As a result, more fluid is dispelled, and the drop in pressure occurs faster. Up to the point when the gas reaches the wellhead, the downhole pressure is at its minimum. (3) The adjustment of wellhead back pressure is mainly realized by throttle valve. The higher the wellhead back pressure is, the smaller the downhole pressure will be. When the gas-liquid two-phase flow in the wellbore reaches an equilibrium state, the downhole pressure will decrease less with the increase of drilling fluid displacement, and the time of gas reaching the wellhead will be earlier. (4) The downhole pressure can be controlled by changing the density of drilling fluid. However, the adjustment of drilling fluid density has a serious lag. Considering the variation of gas invasion caused by reservoir pressure difference, there will be no stable gas-liquid two-phase flow equilibrium.
Further work needs to be done to couple the proposed method with deformation near the well at a small scale.