Prediction of Pressure, Temperature, Holdup, Velocity, and Density

Article pubs.acs.org/IECR

Prediction of Pressure, Temperature, Holdup, Velocity, and Density Distribution for Steady-State Bubbly Gas Three-Phase Flow in High-Temperature−High-Pressure (HTHP) Wells Jiuping Xu,*,† Yunqiang Liu,† Shize Wang,‡ and Bin Qi‡ †

Uncertainty Decision-Making Laboratory, Sichuan University, Chengdu 610064, People’s Republic of China Research School of Engineering Technology, The Southwest Petroleum and Gas Corp, China Petroleum and Chemical Corp, Deyang 618000, People’s Republic of China

‡

ABSTRACT: In this paper, a coupled system model of partial differential equations (PDEs) concerning pressure, temperature, velocity, and holdup for a gas−oil−water three-phase steady flow in high-temperature−high-pressure (HTHP) wells is presented. A solution framework shared by the basic models is built, making it convenient to perform to solve. An algorithm solution model using the fourth-order Runge−Kutta method is given. Basic data from “X Well” (HTHP well, 7110 m deep, in Sichuan, PRC) is used as a case study for calculations and a sensitivity analysis is performed on the model. Pressure, temperature, velocity, and hold-up curve graphs, along with the depth of the well, are plotted at different depths. A trend and sensitivity analysis are conducted to test the parameter characteristics. In addition, a comparison between the present model and other models is conducted to illustrate the effectiveness of the model and the algorithm. The results provide both technical reliability for well test design in HTHP gas wells and provides a dynamic analysis of production.

■

INTRODUCTION Oil−water two-phase flow and gas−oil−water three-phase flow are commonly encountered in the petroleum, chemical process, energy, and nuclear industries. Oil is one of the most important energy sources in the world, and it is closely related to other industrial modernization processes. Oil extraction is often not a single phase, but rather involves both crude oil and associated gas. At the same time, there is a large amount of connate water in the reservoirs, especially in the late stages of oil exploration. Therefore, an oil−water−-gas three-phase flow often occurs in onshore and offshore hydrocarbon production and transportation. Studies on the oil−water two-phase flows and oil−gas−water three-phase flows could solve many important technical problems in the petroleum industry. Studies on the oil−gas−water three-phase flow, especially on the temperature effects and pressure gradients, can solve many important technical problems for the petroleum industry and could also be very important for improvements in multiphase flow theory and the associated practical application. Therefore, it is vital to obtain accurate data for temperature, velocity, and pressure in the tubular string. In a high-temperature−high-pressure (HTHP) deep well, the bottom hole temperature is ∼160 °C with a pressure of ∼35 MPa. The temperature distribution and pressure on the tubing are significantly different when outputs are varied (flow velocity) but neither has a simple linear relationship, because the fluid density is not constant. Therefore, for HTHP wells, the prediction of three-phase flows is important to the industry. In this paper, the three-phase flow is treated as a bubbly gas flow. At the same time, a coupled, steady, mathematical model for HTHP wells is presented. The model can be widely used for oil and gas wells, including horizontal, vertical, and deviated wells. The mass, momentum, and energy equations are solved simultaneously to predict temperature, pressure, density, gas © 2012 American Chemical Society

velocity, oil velocity, water velocity, and holdup. Because of the complexity of the model, a new frame is designed to solve the equation. An algorithm solution model, combined with a fourth-order Runge−Kutta (RK4) method, is recommended to solve the coupled system model. Basic data from X well (a HTHP well) in China, with a depth of 7100 m, is used for case history calculations, and a sensitivity analysis is conducted on the model. The well-posedness analysis of the system and the algorithm stability help to make the system compact and well-organized.

■

LITERATURE REVIEW Generally, two-phase research methods are extended to threephase problems, so research covering both two-phase and three-phase problems are reviewed together. Flow patterns are key when deciding a phase strength situation. About 14 flow patterns had been observed while several researchers described three or four.1 In the past 30 years, there has been significant progress in the analysis of two-phase flow patterns. Many new and comprehensive flow patterns have been published.2−9 Because two-phase gas−liquid flows are highly complex, it is apparent that the addition of a third phase will increase this complexity.10 Experimental observations have shown that the flow structures of a three-phase-pipe flow are much more complicated than that of the two-phase-pipe flow. For example, 10 flow patterns were observed by the experiment study.11,12 Seven (7) flow patterns were identified for horizontal gas−oil−water Received: Revised: Accepted: Published: 6537

September 2, 2011 February 11, 2012 March 8, 2012 March 8, 2012 dx.doi.org/10.1021/ie2019917 | Ind. Eng. Chem. Res. 2012, 51, 6537−6562

Industrial & Engineering Chemistry Research

Article

flow.13 For vertical air−water−oil flow, 8 flow patterns were identified.14 A two-phase flow is the most common flow for fluid in nature, which widely exists in developed oil wells in the middle and later stages. The understanding, description, and prediction of flow characteristics has become a research focus in the field of two-phase flows. Numerical simulation methods have been found to be an effective way to study two-phase flows. A twophase flow in vertical tube research in a gas-lift operation design research was launched.15 The mixture was treated as homogeneous and single phase. The density was constant and the viscosity was ignored. This method was found to be suitable for a high flow and low gas−liquid ratio, but made large errors in other types. A pressure drop computation method was proposed.16 The oil−liquid two-phase flow was described using dimensional analysis and a mobility pattern distribution map was presented with the experiment, which proved to have very accurate engineering but could not be used for deep or great pressure drop wells. Empirical slip models suitable for a long tube flow-based gas−liquid two-phase slip was built.17 It was established that most energy loss was caused by friction, which could be correlated with a Reynolds number. The friction factor needed to be fixed by calculating the liquid holdup, and, although the equation involved empirical correlations, it had high application accuracy and has been widely used in horizontal multiphase flows. A two-phase flow pressure drop computation model covering all patterns in a vertical tube was established and the flow pattern identification method was developed.18 The technique proved to be one of most reliable application methods in the oil industry. A pressure drop model was developed for bubbly flow and slug flow focusing on a gas−liquid two-phase flow, where a gas volume factor was introduced into the density and friction losses using a gas− liquid two-phase separating effect.19 The relationship of a liquid hold-up and drag coefficient based on a homogeneous flow pressure gradient equation based on homogeneous flow was experimentally studied.20 Beggs−Brill’s correlations were built and could be used for almost any pipeline inclination. The flow pattern changing conditions from the mechanism was explained and a physical model was proposed.21 Their model had an important impact on two-phase flow research, in that the pressure drop calculation method was from empirical research to mechanism study. Early correlations were empirical and were based on experimental studies, and the results were generally satisfactory for the conditions under which each model was developed. Marktos et al. had made outstanding contributions to the two-phase work. They presented a one-dimensional, isothermal, flow analysis of steam and water mixtures in vertical flow passages in a general form. They extended the model to applications such as fire sprinkler systems, batch sedimentation, granular propellants, pressurized-water-reactor design and airlift pumps simulation, etc.22−28 Because of the complex nature of the gas−liquid two-phase flows in wells and pipes, many attempts have been made to develop predicable techniques through empirical or semiempirical methods (such methods have dominated practical design procedures). However, most research calculates the distribution prediction of temperature and pressure separately,29−33 and their interdependence is ignored. In fact, it is well-known that there is interdependence in pressure, temperature, density, and velocity in injection wells and production wells. The characteristics of the oil−water twophase pipe flow was experimentally studied and found that a droplet size distribution will directly affect the accuracy of

model predictions.34 Hasan et al. made good jobs to the coupled fluid flow model where wellbore/reservoir simulators for modeling single-phase gas, oil and the two-phase gas−oil flow problems were presented. The transient fluid and heat flows models are solved numerically using finitedifference methods to obtain some parameters such as pressure and velocity. The numerical algorithms apply a double-iterative procedure on both temperature and pressure to solve the three conservation equations simultaneously.35−40 The hydrodynamic three-phase modeling is based on flow pattern definitions and some forces are different in each pattern. More flow patterns mean more discontinuities and greater complexity in the models. Since there are many papers identifying the patterns of flow, a model using one pattern can be used to analyze the forces and can then be applied to the other patterns. Among the research on limited coupled models, there were two main ideas. One treatment for a three-phase flow is to consider the system of oil, water, and gas as one liquid phase with mixture properties.41,42 However, the model can predict neither the parameter mentioned nor the volumetric fraction for each phase. The other treatment for a three-phase flow is to combine oil and water into a single liquid phase. Gas−liquid−liquid three-phase flows can be regarded as a special type of gas−liquid two-phase flow if the two liquids are fully mixed. This is probably true of vertical and steeply inclined flows. The physical properties of the liquid mixture can be calculated based on the fractions and the individual physical properties of the two liquids. In this treatment, the slip between the oil and water is ignored and a homogeneous mixture is assumed for the liquid phase.9,10,43,44 The other extreme is to a treat three-phase flow as a three-layer stratified flow with gas on the top, oil in the middle, and water at the bottom. This can be done for immiscible liquids flowing in horizontal or slightly inclined pipes with low gas, oil, and water flow rates. The solution is further complicated for two phases and becomes much more complex for three phases. A model was developed to predict the value of the hold-up and pressure gradient for a three-phase stratified flow in a horizontal pipeline.45 The concept of extended velocity was applied to compute the wall shear stresses. However, the temperature profile was not considered. A three-phase (heavy-oil−water− gas) bubbly flow in upward pipes was simulated using a onedimensional transient two-fluid model in which the continuity and momentum equations for the two liquids (heavy oil and water) were combined to obtain a new equation for the liquid mixture quantities.10 However, the parameters those models predicted are the pressure profiles rather than the temperature profiles and failed to consider heat transfer.

■

MODEL BUILDING Problem Description. Within the tubing, flow takes place under turbulent flow conditions. Consider the flow system depicted in Figure 1: a straight cylindrical flow tube with an inclination angle θ, a constant cross-sectional flow area A, a hydraulic diameter d, and a total length Z. Gas flows through the tubing from the bottom to the top with a mass flow rate W. The distance coordinate in the flow direction along the tubing is denoted by z. Mass, momentum, and energy balances, along with the pressure, temperature, and velocity, related to the steady gas−liquid−liquid three-phase flow, as well as the stated equation for gas, are used to generate the constitutive equation. 6538

dx.doi.org/10.1021/ie2019917 | Ind. Eng. Chem. Res. 2012, 51, 6537−6562


Article

Gas flows continuity equation: d(ρg αg νg A) dz

=0

(3)

Momentum Equation. A unified momentum equation can be presented as follows: sum vector of the forces = the momentum flowing out from dz

Figure 1. Schematic of the physical figure.

− the momentum flowing into dz

The force on dz include the following: • The pressure on cross section is given as [αkPk − (αkPk + αk(dPk/dz)dz)]A = −αk(dPk/dz)A dz • The gravity which flows in the opposite direction: ρkαkg(cos θ)A dz • The friction which flows in the opposite direction: τkbSkb dz • The interfacial shear stress, which flows in the opposite direction: τkjSkj dz The momentum flowing out from dz is given as ρkαkυ2k + A(d(ρkαkυ2k)/dz). The momentum flowing into dz is given as ρkαkυ2kA. The subscript k represents gas (g), oil (o), and water (w) and implies a phase that is different from phase k. From the law of momentum conservation, the momentum equation for each phase is as follows: Water flow momentum equation:

Basic Assumption. Considering the differential equation model of P, T, v and a, the following assumptions were set. (1) The gas−liquid−liquid flow in the tubing is in one dimension of the flowing direction, with negligible heat or mass transfer between phases. (2) When the oil−gas−water three-phase system reaches thermodynamic equilibrium, the temperatures of all points are equal in the cross-section of the arbitrary position. (3) The oil−water two phases flow is with a flat interface. (4) Oil and water can be supposed to be incompressible phases, because of their small compressibility. (5) The phase pressures are assumed to be equal at a given axial location. (6) The gas flow is approximately considered to be that of an ideal gas.

d(ρωαωυ2ω) dz

= −αω

τ S τ S dP − ρωαωg cos θ − wb wb − wo wo dz A A (4)

Oil flow momentum equation: d(ρoαoυo2) dz

Gas flow momentum equation: d(ρg αg υ2g )

Model Formulation. In the paper, a three-phase system in a infinitesimal section is considered as shown in Figure 2. Referring to Figure 2, the flow of three fluids is considered. It is assumed that water is heavier than oil and flows at the bottom. The oil flows at the top. Material Balance. The top of the well is considered to be the origin of the coordinate axis and vertical up is the positive direction. Let dz denote the differential depth. Thus, under steady-state conditions, for each phase, it follows the law of fluid dynamics.

dz

dz

dz

=0

dP − ρg αg g cos θ dz

(6)

translate energy the second net = energy flowing out from dz − energy flowing into dz

The types of energy considered in this paper include inner energy, pressure energy, kinetic energy, and potential energy. The energy flowing into dz includes the following types: • The internal energy: Wk(z) • The kinetic energy: mkυ2k(z)/2 • The potential energy: mkgz cos θ • The pressure energy: Pk(z)υk(z)

(1)

Oil flows continuity equation: d(ρoαoνoA)

= −αg

Here, τwbSwb and τobSob denote the friction shear stress between the water and the wall and between the oil and the wall, respectively. τ woSwo and τowSow (τwoSwo = τowSow) represent the shear force between the oil phase and the water phase. Energy Equation. The unified energy equation is given as

Water flows continuity equation: =0

τ S τ S dP − ρoαog cos θ − ob ob − ow ow dz A A (5)

Figure 2. Bubbly gas three phase flow.

d(ρωαωνωA)

= −αo

(2) 6539



Article

Figure 3. Schematic depiction of the radial transfer of heat.

The energy flowing out from dz includes the following types:

The term ηk is the Joule−Thomson coefficient, which is defined as10

• The internal energy: Wk(z + Δz) • The kinetic energy: mkυ2k(z + Δz)/2 • The potential energy: mkg(z + Δz) cos θ • The pressure energy: Pk(z + Δz)υk(z + Δz) where the inner energy and pressure energy are collectively referred to as enthalpy (Hk(z) = Wk(z) + Pk(z)υk(z)) and dQ represents the radial heating of the tubing. According to the energy balance rule, the energy of fluid flowing into the distance element equals the energy sum of losses and fluid flowing out from the distance element:

⎧η = 0 k=g ⎪ k ⎨ 1 k = o, w ⎪ ηk = − C pk ρk ⎩

In eq 7, dqk denotes the radial heat transfer between the gas and the surrounding earth tube. The radial transfer of heat between the fluid and the earth have been discussed in detail.46,47 As shown in Figure 3, the radial heat transfer from the fluid to the cement/earth interface is described by

mk υk2(z) − mk gz cos θ 2 m υ2(z + Δz) = Hk(z + Δz) + k k − mk g (z + Δz) cos θ 2

Hk(z) +

+ dQ

2πrt 0Ut 0 (T − Tk) dz w

dqfe =

(7)

2πKe (Tk − Te) dz wf (tD)

dqes =

dz

+ υk

d υk dh − g cos θ + =0 dz dz

dq =

2πrt 0Ut 0Ke (T − Te) dz w(rt 0Ut 0f (tD) + Ke)

(14)

Let a = (2πrt0Ut0Ke)/[w(rt0Ut0 f(tD) + Ke)], substituting eqs 10, 11, and 14 into eq 9. Then the energy equation for the water phase is as follows:

(9)

The term hk, which satisfies the following relation, denotes the specific enthalpy. dhk dT dP = Cpk k − ηk Cpk k dz dz dz

(13)

Thus, the radial heat transfer between the gas and the surrounding earth tube is

(8)

Division by mk yields dqk

(12)

and the radial heat transfer from the cement/earth interface to the surrounding earth is described by

Equation 7 can be written as dHk dυ dQ + mk υk k − mk g cos θ + =0 dz dz dz

(11)

ρwCpw (10)

dTw dυ dP + ρwυw w + w − ρwg cos θ + aρw(T − Te) dz dz dz

=0 6540

(15) dx.doi.org/10.1021/ie2019917 | Ind. Eng. Chem. Res. 2012, 51, 6537−6562


Article

Similarly, the oil flow energy equation can be written as ρoCpo

For one thing, the number of differential equations for a three-fluid model is twice that for a two-fluid model. For another, the nonlinearity and coupling of the equations are much stronger, because of the interactions between the three fluids. The equations are coupled not only within those for the same phase, but also among those for different phases. Because of the coupling of equations, the work to solve the equations directly is very difficult. In this paper, a simplified model is applied. First, considering that the gas density is only a function of pressure (ρg = f(Pg)),48 then

dTo dυ dP + ρoυo o + o − ρog cos θ + aρo(T − Te) dz dz dz (16)

=0

From Figure 2, the gases are supposed to be fully dispersed in oil−water two-phase flows. Therefore, the gas energy equation is not considered here. Finally, a gas state equation can be added: ρg =

MPg RZgT

dρg (17)

dx

Combining eqs 1−6 and 15−17, the coupled system model of differential equations is as follows: ⎧ d(ρ αωνωA) ⎪ ω =0 dz ⎪ ⎪ ⎪ d(ρoαoνoA) = 0 ⎪ dz ⎪ ⎪ d(ρg αg νg A) =0 ⎪ dz ⎪ ⎪ d(ρ α υ2 ) ⎪ ω ω ω = −α dPw − ρ α g cos θ ω ω ω ⎪ dz dz τwbSwb τ S ⎪ − − wo wo ⎪ A A ⎪ ⎪ d(ρ αoυo2) τ S dP ⎪ o = −αo o − ρoαog cos θ − ob ob dz A ⎪ dz τ S ⎪ − ow ow ⎪ A ⎪ 2 ⎨ d(ρg αg υg ) dPg ⎪ = −αg − ρg αg g cos θ dz dz ⎪ ⎪ ⎪ρ C dTw + ρ υ d υw + dPw − ρ g cos θ w w w ⎪ w pw dz dz dz ⎪ + aρw(Tw − Te) = 0 ⎪ ⎪ dT dυ dP ⎪ρoCpo o + ρoυo o + o − ρog cos θ dz dz dz ⎪ + aρo(To − Te) = 0 ⎪ ⎪ MPg ⎪ ⎪ρg = RZgT ⎪ ⎪ ⎪ Pw = Po = Pg = P , ⎪ αw + αo + αg = 1, ⎪ ⎪ Tw = To = Tg = T , ⎪ y(z 0) = φ(z 0) ⎩

dρg dP 1 dP = 2 dP dx c dx

=

(19)

where c is the sound velocity, which is calculated using the following equation:49 M ZRT

c=

(20)

So eqs 1−3 can be written as follows, respectively: d υw dα + υw w = 0 dz dz

(21)

d υo dα + υo o = 0 dz dz

(22)

αw αo

υg αg P

d υg dα dα dP − ρ2g υg o − ρ2g υg w + ρ2g αg =0 dz dz dz dz

(23)

Substituting eqs 2, 3, and 23 into the momentum equation, respectively, the following equations are obtained: 2ρwαwυw

d υw dP + αw dz dz

= −ρwg αw cos θ −

2ρoαoυo

τwbSwb τ S − wo wo A A

(24)

d υo τ S τ S dP + αo = −ρog αo cos θ − ob ob − ow ow dz dz A A (25)

2ρ2g αg υg =

d υg

+ αg (υ2g P + ρg )

dz

−ρ2g g αg

d αg dP + ρ2g υ2g dz dz

cos θ

(26)

Adding eqs 24 and 25 to eliminate the shear force between two phases, the following equation is acquired: (18)

■

d υw dυ dP + 2ρoαoυo o + (αw + αo) dz dz dz τwbSwb τ S = −g cos θ(ρwαw + ρoαo) − − ob ob A A

2ρwαwυw

SOLVING MODEL Model Analysis. A three-fluid model system is obviously more complex to compute than a one-fluid or two-fluid system. 6541

(27)



Article

Combining eqs 15, 16, 21−23, and 27, the coupled system of differential equations on pressure, velocity and temperature in mathematical form is obtained and the well-posedness of the model is analyzed in Appendix A. ⎧ d υg dα dα dP ⎪ υg αg P − ρ2g υg o − ρ2g υg w + ρ2g αg =0 dz dz dz dz ⎪ ⎪ dυ ⎪ αw w + υw d αw = 0 ⎪ dz dz ⎪ ⎪ α d υo + υ d αo = 0 o ⎪ o dz dz ⎪ ⎪(α + α ) dP + 2ρ α υ d υw + 2ρ α υ d υo o w w w o o o ⎪ w dz dz dz ⎪ τwbSwb + τobSob ⎪ = − g cos θ(ρwαw + ρoαo) − A ⎪ ⎪ d υg dα dP ⎪ 2ρ2g αg υg + αg (υ2g P + ρg ) − ρ2g υ2g w dz dz dz ⎪ ⎪ 2 2 d αo 2 g cos − ρ υ = −ρ α θ g g g g ⎪ dz ⎪ dυ ⎪ dP dT + ρwυw w + ρwCpw ⎪ dz dz ⎨ dz ⎪ = ρwg cos θ − aρw(T − Te) ⎪ d υo ⎪ dP dT ⎪ dz + ρoυo dz + ρoCpo dz ⎪ ⎪ = ρog cos θ − aρo(T − Te) ⎪ ⎪ρ = MP ⎪ g RZgT ⎪ ⎪ P(z ) = P , 0 ⎪ 0 υg (z 0) = υg 0 , ⎪ ⎪ υo(z 0) = υo0 , ⎪ υw(z 0) = υw0 , ⎪ ⎪ ⎪ αw(z 0) = αw0 , ⎪ αo(z 0) = αo0 , ⎪ Tw(z 0) = T0 , ⎪ ⎪ αw + αo + αg = 1 ⎩

(28)

The modified system of equations that can be written in a compact form as

A

dU =B dz

Figure 4. Flowchart showing the overall procedure of the model.

where A represents the coefficient matrices, B is a vector containing all algebraic terms, and U is the solution vector.

(29)

⎡ ρ2g (1 − αw − αo) −ρ2g υg −ρ2g υg vg (1 − αw − αo)P 0 0 0 ⎤ ⎥ ⎢ ⎢ αw υw 0 0 0 0 0 ⎥ ⎥ ⎢ ⎢ αo υo 0 0 0 0 0 ⎥ ⎥ ⎢ αw + αo 2ρwυwαw 2ρoυoαo 0 0 0 0 ⎥ A=⎢ ⎥ ⎢ ⎢(1 − αw − αo)(υ2g P + ρg ) 0 0 2ρ2g (1 − αw − αo)υg −ρ2g υ2g −ρ2g υ2g 0 ⎥ ⎥ ⎢ ⎢ ρwυw ρwCpw ⎥ 1 0 0 0 0 ⎥ ⎢ ⎢⎣ ρoυo ρoCpo ⎥⎦ 1 0 0 0 0

6542



Article

⎡ ⎤ 0 ⎢ ⎥ 0 ⎢ ⎥ ⎢ ⎥ 0 ⎢ ⎥ ⎢−g cos θ(ρ α + ρ α ) − τwbSwb + τobSob ⎥ w w o o ⎥ A B=⎢ ⎢ ⎥ 2 −ρg g (1 − αw − αo) cos θ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ρwg cos θ − aρw(T − Te) ⎢ ⎥ ⎢⎣ ⎥⎦ ρog cos θ − aρo(T − Te)

⎡P ⎤ ⎢υ ⎥ ⎢ w⎥ ⎢ υo ⎥ ⎢ ⎥ U = ⎢ υg ⎥ ⎢α ⎥ ⎢ w⎥ ⎢ αo ⎥ ⎢ ⎥ ⎣T ⎦

aforementioned discussion, the fourth-order Runge−Kutta (RK4) method was used to simulate the model. The coupled models for three-phase flow boundary conditions are defined at the bottom of the well. Some parameters needed are shown in Appendix C. Solution Process. The proposed algorithm procedure for solving the problem is designed as follows. The overall algorithma̧rś program flow diagram is presented in Figure 4. Step 1: Set the step length of depth. In addition, the relatively tolerant error is denoted by ε. The smaller Δλ and ε, the more accurate the results. However, it can lead to an rapid increase in the calculation time. Here, set Δh = 250(m), Δλ = 1, and ε = 5%. Step 2: Compute each point’s inclination. Step 3: Give the initial conditions and let h = 0. Step 4: Compute the parameters under the initial conditions or the last depth variables and let λ = 0. Step 5: Let the differential equations to be functions f i, where (i = 1, 2, 3, 4, 5, 6, 7). A system of coupled functions then is derived as follows:

Several numerical methods can be used to solve ordinary differential equations (ODEs), such as the Runge−Kutta methods, linear multistep methods, and predictor-correcting methods. The Runge−Kutta technique is far more widely used than any of the other techniques for the defined initial conditions.50 The stability of the numerical algorithm is given by Appendix B. Equation 29 represents large linear equations, and there are many ways to solve them. Considering the large and sometimes ill-posed linear systems, in this paper, a Singular Value Decomposition (SVD) Method, which can make the solutions stable, is adopted.51 Numerical Solution. To simplify the calculations, the wells were divided into several short segments of the same length h. The length of a segment varies, depending on variations of well thickness, hole diameter, fluid density inside and outside the pipe, and the geometry of the wells. The velocity, pressure, and temperature calculations then are performed for each successive “segment” of the pipe up to the surface. According to the

⎧ υ (1 − α − α )Pf − ρ2 υ f − ρ2 υ f + ρ2 (1 − α − α )f = 0 w o w o 4 g g6 g g5 g 1 ⎪ g ⎪α f + υ f = 0 w5 ⎪ w2 ⎪α f + υ f = 0 o6 ⎪ o3 ⎪ ⎪ τ S τ S ⎨ 2ρwαwυwf2 + 2ρoαoυof3 + (αw + αo)f1 + P(f5 + f6 ) = −ρwg αw cos θ − ρog αo cos θ − wb wb − ob ob A A ⎪ ⎪ 2 2 2 2 2 2 2 ρ − α − α υ + − α − α υ + α − ρ υ − ρ υ = −ρ − α − α θ 2 (1 ) f (1 )( P ) f f f g (1 ) cos w o g 4 w o g g 1 w o g g6 g g5 g ⎪ g ⎪ ⎪ρwCpwf7 + ρwυwf2 + f1 − ρwg cos θ + a(T − Te) = 0 ⎪ ⎪ρ Cpof + ρ υof + f − ρ g cos θ + a(T − Te) = 0 ⎩ o 7 o 3 o 1

Step 6: Assume that P, υw, υo, υg, αw, αo, T to be yi(i = 1, 2, 3, 4, 5, 6, 7, respectively). By solving the equations using SVD methods, some basic parameters are obtained as follows: ⎧ ai = f (y , y , y , y , y , y , y ) i 1 2 3 4 5 6 7 ⎪ ⎪b = f (y + λa /2, y + λa /2, y + λa /2, y + λa /2, y + λa /2, y + λa /2, y + λa /2) 1 2 3 4 5 6 7 ⎪ i i 1 2 3 4 5 6 7 ⎨ ⎪ ci = fi (y1 + λb1/2, y2 + λb2 /2, y3 + λb3/2, y4 + λb4 /2, y5 + λb5/2, y6 + λb6 /2, y7 + λb7 /2) ⎪ ⎪ di = f (y + λc1, y + λc 2 , y + hc3 , y + λc4 , y + λc5 , y + λc6 , y + λc 7) ⎩ i 1 2 3 4 5 6 7

yi j + 1 = yi j +

Step 7: Calculate the gas velocity, water velocity, oil velocity, oil cut, water cut, pressure, and temperature at point (j + 1):

λ(ai + 2bi + 2ci + di) 6

i = 1, 2, 3, 4, 5, 6, 7; j = 1, 2 , ..., n 6543



Article

Table 1. Pipe Parameters diameter (mm)

thickness (mm)

weight

expansion

coefficient

Young

Modulus

88.9 88.9 88.9 73 73

9.53 7.34 6.45 7.82 5.51

18.9 15.18 13.69 12.8 9.52

0.0000115 0.0000115 0.0000115 0.0000115 0.0000115

215 215 215 215 215

0.3 0.3 0.3 0.3 0.3

1400 750 4200 600 150

Step 8: Compute the gas density:

calculations are performed for each successive segment of the pipe up to the surface.

■

MPkj

ρ gj =

RZgTkj

PARAMETERS In this simulation, a pipe in X well, China, is used as a case study. All the needed parameters are as follows: depth of the well, 7100 m; critical pressure, 4.968 MPa; gas specific weight, 0.6434; ground thermal conductivity parameter, 2.06; ground temperature, 16 °C; ground temperature gradient, 0.0218 °C/m; the roughness of the inner surface of the well, 0.000015; the parameters of pipes, inclined well, inclination and vertical depth are given in Tables 1−3. Some data, needed by the model, are taken from Sagar et al.,30 as shown in Table 4.

Step 9: λ = λ + Δλ. Repeat steps 6−8 until λ > Δh/Δλ. Step 10: h = h + Δh. Repeat steps 4−9 until h < hmax is calculated.

■

NUMERICAL SIMULATION AND RESULTS DISCUSSION Some models were validated using the experimental data, and it was found to be an effective method that must be analyzed more deeply.10,11,43 To study the suitability of the model, some parameters were adjusted to reflect totally different operating conditions. However, as a result of the test conditions, it was not possible to totally simulate the real situation, so the error was large. Up until now, a 4570-m-deep well has been considered a deep well, and a 7100-m-deep well has been considered a superdeep well, so a completion test for a deep well has become an emerging new problem. In the applied basic theory for deep well testing research, tubular string mechanical analysis has been shown to be extremely complex as fluid temperature and tubing pressure heavily affect the force of the tubular string. Therefore, because, in such an experiment, it is difficult to model the real tubular sting, here the real data are derived by comparing the consistency between the model and the real situation. As was described previously, the algorithm begins with a calculation for the pipe at the bottom of the pipe. The

Table 4. Model Parameters

depth (mm)

internal

4325.69 6301.7 7100

168.56 168.3 121.42

193.7 193.7 146.1

value

CPo Ke kang CPw kcem w

0.485 Btu/(lbm °F) 33.6 Btu/(lbm °F) 0.504 Btu/(lbm °F) 1.0 Btu/(lbm °F) 96.5 Btu/(lbm °F) 0.08 lbm/s

■

STEP ANALYSIS The effect of grid fineness had to be investigated. First, the adequacy of the number of the grids in the fourth-order Runge−Kutta (RK4) method to obtain independent results is tested. The maximum number of steps should be no less than 10 (see Appendix B). As shown in Figure 5, grids of 125, 250, 500, and 2500 steps were chosen where the pressure along the depth becomes practically grid-independent for the 250- and 500-step grid solutions. An interesting phenomenon was found that the more-precise results were not gained with a higher number of grids, mainly because errors tended to increase from accumulation.

Table 2. Casing Parameters measured (mm)

parameter

Table 3. Parameters of Azimuth, Inclination, and Vertical Depth number

measured

inclination

vertical depth

number

measured

inclination

vertical depth

1 2 3 4 5 6 7 8 9 10 11 12

0 303 600 899 1206 1505 1800 2105 2401 2669 3021 3299

0 1.97 1.93 0.75 1.25 1.04 0.49 2.49 1.27 2.44 0.14 1.18

0 302.87 599.73 898.59 1205.45 1504.32 1799.18 2104.04 2399.91 2667.79 3019.63 3297.5

13 14 15 16 17 18 19 20 21 22 23 24 25

3605 3901 4183 4492 4816.07 5099.07 5394.07 5706.07 5983.07 6302.07 6597.07 6911.12 7100

2.05 0.16 2.92 2.73 1.98 2.74 0.13 0.63 2.09 2.69 2.45 0.15 1.15

3603.36 3899.22 4181.09 4489.95 4813.87 5096.74 5391.61 5703.47 5980.34 6299.19 6594.06 6907.96 7085.88

6544



Article

Figure 5. Step analysis.

Figure 6. Trend analysis.

■

TREND ANALYSIS

with increased depth, as shown by Figure 5b. The temperature was 162.02 °C at the wellhead. Figure 6a shows the volumetric fraction profiles for each phase. The gas volumetric fraction minimum value reached was almost 79.52%, which is larger than both the oil (maximum value = 8.1%) and the water (maximum value = 12.38%), which seems to indicate that this well is a well that yields natural gas, which is obviously in a dominant continuous phase. The decrease in the gas volume fraction is more pronounced in the wellhead than at the well bottom, because the gas phase accelerates rapidly in the wellhead. Figure 6b shows the velocity simulation results. In accordance with mass conservation theory, the three-phase velocity changes with depth in a trend opposite to those of

A trend analysis was performed to determined whether the developed model is physically correct or not. To test the developed model, the distribution of oil velocity, gas velocity, water velocity, water cut, oil cut, pressure, temperature, and density in the wells were determined. Figure 5a shows that the pressure multiplies as the depth increases with the wellhead pressure at 24.6 MPa at a fixed production. The reason is that the flowing gas phase tends to overcome the gravity and friction. Because the pressure drop is influenced by liquid velocity and holdup, the pressure drop grows larger and larger, so the larger the liquid velocity, the larger the pressure drop. Meanwhile, the phase also transmits heat to the stratum, which causes the temperature to decrease 6545



Article

Figure 7. Gas production sensitivity analysis.

6546



Article

Figure 8. Sensitivity analysis for samples with and without heat transfer. 6547



Article

Figure 9. Comparison analysis.

fraction is influenced by the liquid velocity and the decrease in liquid velocity. The effect of an increase in velocity at high gas production is almost balanced out by the increase in the velocity drop effect. At the same depth, the greater the gas production, the larger the gas velocity, as shown by Figure 7e. Maximum gas velocity reached 8.613 m/s at the wellhead. However, the damage in nodular cast iron happens when the velocity is more than a certain amount. Generally, from an engineering point, the velocity should be controlled to less than 30 m/s. The same trend can be seen for the water and oil phases from Figure 7f and 7g, respectively (the wellhead water and oil velocity were 6.82 m/s and 6.65 m/s). The gas density decreased with an increase in production, as found in Figure 7h. When the gas production was from 300 000 m3/day to 700 000 m3/day, the wellhead pressure reached a range from 157.02 kg/m3 to 122.1 kg/m3. This is because the pressure was decreasing as the temperature was increasing. The Values with and without Heat Transfer. The other parameter chosen was heat transfer, because heat transfer has been neglected in some models;10,41 therefore, the purpose was to investigate how important this factor was. Figure 8a shows the pressure distributions along the depth when considering both measurements with heat transfer taken into account and those without heat transfer. From this, it can be seen that the heat transfer had almost no effect on the pressure (the pressure at the wellhead decreased from 29.19 MPa with heat transfer to 29.01 MPa without heat transfer). However, heat transfer had a strong effect on temperature and the difference became greater as the depth decreased, as is shown in Figure 8b. The wellhead temperature decreased nearly 3 °C. Yet, from the present model’s calculation accuracy, heat transfer barely affected the velocity and volume fraction profiles shown in Figures 8c−g. Figure 8h shows the gas density profiles. At the same depth, the heat transfer had a weak effect on gas density (from 149.1 kg/m3 to 147.7kg/m3). This is because the gas density was influenced not only by temperature but also by other parameters, according to the density equations. Therefore ,the main influence on gas density is gas production. The reason why the heat transfer appeared to only affect temperature is that the energy equation was not coupled with material or momentum equations. If this equation is only used for only pressure and velocity profiles, an energy equation is not necessary; these results reflect those of many researchers.44,52−54

pressure and liquid holdup. The primary reason for the decline in the oil or water velocity with depth is that an increase in liquid (gas or water) volumetric fractions requires a decrease in velocity to satisfy mass balance. The gas velocity is the largest at the same initial conditions. The gas velocity at the wellhead was 6.388 m/s, while that of the oil and water were 5.001 m/s and 5.205 m/s, respectively. This gas velocity trend was also observed by Cazarez et al.,10 who concluded that this was because the gas density was the smallest. However, here, the oil velocity was not as small, which could be because the volumetric fractions made the velocity change later. As shown in Figure 6c, gas density decreases along with depth and the minimum density reached 153.0781 kg/m3 at the wellhead. At a depth of 3750 m, the parameters had a sudden change (water and oil were not obvious because there were only slight variations in volumetric fractions). The reason was that the interior diameter of tubing changed from 73 mm to 88.9 mm. Sensitivity Analysis. During this analysis, two parameters that influenced the result were taken into consideration. The first was gas output with three different measures used: 300 000 m3/day, 500 000 m3/day, and 700 000 m3/day. In order to study the difference between a model with a radial transfer of heat to the tubing and a model without it, the results allowed for the two conditions. Through the algorithm and by simulation, a series of results was obtained. The Values at Different Gas Outputs. For this purpose, synthetic sets were prepared where, in each set, only the gas output parameter was changed while other parameters were kept constant. Figure 7a shows the pressure profiles with variable gas production. At the same depth, pressure decreases with an increase in production. This is because flowing velocity increases with an increase in output resulting in a drop in increased frictional pressure. When the gas production ranged from 300 000 m3/day to 700 000 m3/day, the wellhead pressure ranged from 31.01 MPa to 24.6 MPa. Figure 7b shows the temperature profiles with varied gas production. At the same depth, temperature increases with an increase in production. However, the degree of increase decreases. This is because an increase in production results in an increase in velocity; thus, the heat loss decreases. When the gas production was from 300 000 m3/day to 700 000 m3/day, the wellhead temperature ranged from 160.045 °C to 162.75 °C. The oil and water volumetric fraction, respectively, saw almost no change with the varied gas productions, as can be seen in Figures 7c and 7d. The liquid (oil or water) volumetric 6548



Article

Table 5. Comparative Results wellhead

pressure (MPa)

temperature (°C)

gas velocity (m/s)

gas density (kg/m3)

the present model Octavio’s homogeneous model relative error

29.19 26.4 9.5%

159.08 156.5 1.6%

6.388 6.302 1.3%

149.1 147.2 1.2%

Comparison Analysis. In the comparison step, because of the lack of experimental or actual field data for a three-phase pipeline flow, a homogeneous model was developed for the pressure prediction, temperature, and velocity distribution based on a two-phase flow in oil wells. Therefore, the performances of pressure and temperature were adopted for contrasting with the present model. Figures 9a and 9b show that the pressure and temperature measures close to the well bottom were much more accurate than those near the wellhead. This was probably because the gas volume fraction at the well bottom was very high and therefore could be considered as a single phase. As the gas volume fraction decreases along the depth, errors may accumulate. Therefore, the maximum discrepancy over the depth is at wellhead and the max error in pressure and temperature are 9.5% and 1.6%, respectively. Other computational parameters are shown in Table 5. Computer-Time Requirements. The computation time on a Lenovo G470 microcomputer with an i5 CPU was 23 s for a typical 10-node calculation.

■

CONCLUSION A coupled system model of pressure, gas velocity, liquid velocity, and temperature differential equations in high-temperature−high-pressure (HTHP) wells was presented according to the mass, momentum, and energy balances. The features of the models are summarized below. (1) The wellposed-ness of a three-phase system was investigated, and the fourth-order Runge−Kutta (RK4) methods stable step was also computed. Therefore, a solution process was used to predict the model. (2) The predictions, based on different parameters such as gas production or heat transfer, are presented. It is physically obvious that the gas production is the main influential factor. (3) Further investigations, for example with different flow-regime and transient, are desirable. (The volume fraction change was so small found from the model and it can be simplified to fit for large equations in transient flow model) (4) Although the paper is inclined toward pipelines, the relevant model is already capable of handling horizontal or vertical wells.

■

APPENDIX

A. The Well-Posedness Analysis of the System

The task here is to check the nature of the proposed model, and to establish whether the system is stable. Like Cazarez et al.,10 the equations for the two liquids (oil and water) were combined to obtain the mixture. The modified model is given as follows. ⎧ ⎛ f ρ υ2 ⎞ −αg (ρ2g + ρg υ2g P − 2P υ2g )⎜ 2ld l + ρlg cos θ(1 − αg )⎟ + αg (1 − αg )ρ3g g cos θ ⎪ dα ⎝ ⎠ ⎪ g = ⎪ dz 2ρ2g υ2g − ρ3g υ2g − 2αg ρ2g υ2g + αg ρ3g υ2g + 2αg ρ2g ρlυl2 − 4P αg υ2g ρlυl2 + 2P αg ρg ρlυ2g υl2 ⎪ ⎪ d αg υl dz ⎪ d υl = ⎪ 1 − αg ⎪ dz ⎪ f ρ υ2 ⎪ dυ −ρlg cos θ(1 − αg ) − 2ld l − 2ρl(1 − αg )υl dzl ⎨ dP = ⎪ 1 − αg ⎪ dz ⎪ α d g dP ⎪ dυ αg υg P dz + ρ2g υg dz ⎪ g =− 2 ⎪ dz ρ g αg ⎪ ⎪ dυ dP ρlg cos θ − α(T − Te) − dz − ρlυl dzl ⎪ dT = ⎪ dz ρlCPl ⎩

where αl = αo + αw ρl =

ρw αw + ρo αo αl

νl = νo = νw CPl =

ρw CPw + ρo CPo αlρl

Tl = To = Tw

6549



Article

Let dαg dz

(

= f1 z ; αg , υl , P , υg , T

dυl = f2 z ; αg , υl , P , υg , T dz dPg dz dυg dz

)

(

)

(

)

(

)

= f3 z ; αg , υl , P , υg , T = f4 z ; αg , υl , P , υg , T

dT = f5 z ; αg , υl , P , υg , T dz

(

) T

Denote F = ( f1,f 2,f 3,f4,f5)T, y = (αg,υl,P,υg,T)T, y(0) = y(z0) = (αg(z0),υl(z0),P(z0),υg(z0),T(z0)) ; then, the systems of ODEs can be written y′ = F(z ; y),

y(0) = y(z 0)

The norm of vector f was given as follows: F || = max(|f1 | , |f2 | , |f3 | , |f4 | , |f5 |)

For f1, f 2, f 3, f4, f5, ⎛ f ρ υ2 ⎞ −αg (ρ2g + ρg υ2g P − 2P υ2g )⎜ 2ld l + ρlg cos θ(1 − αg )⎟ + αg (1 − αg )ρ3g g cos θ ⎝ ⎠

|f1 | =

|2ρ2g υ2g − ρ3g υ2g − 2αg ρ2g υ2g + αg ρ3g υ2g + 2αg ρ2g ρlυl2 − 4P αg υ2g ρlυl2 + 2P αg ρg ρlυ2g υl2|

≤

|K1||K 2| + |K3| |K 4 |

where |K1| = |− αg(ρ2g + ρgυ2g P − 2Pυ2g ) | ≤ | ρ2g |+| ρg ||υ2g ||P| + 2|P||υ2g |, |K2| = |( fρgυ2g /(2d)) + ρlg cos θ (1 − αg) | ≤ 2 ((|f ||ρl||υ2l |)/|d|)| + |ρl||g|, |K3| = |αg(1 − αg)ρ2g g cos θ | ≤ | |ρ2g ||g|, |K4| = |2ρ2g υ2g − ρ3g υ2g − 2αgρ2g υ2g + αgρ3g υ2g + 2αgρ2g ρlυ2l − 4Pαgρ2g ρlυ2l + 2Pαgρgρlυ2g υ2l |. All parameters are bounded quantities; therefore, |K1|, |K2|, |K3|, |K4| are bounded. Let N1 = sup{((|K1||K2| + |K3|)/(|K4|))}; then, |f1| ≤ N1. Similarly,

|f2 | = |

|f3 | ≤

|f4 | ≤

|f5 | ≤

υl υl f1 | ≤ | |N1 1 − αg 1 − αg

|ρl||g | + |

|f | |ρl| |υl2| 2 |d |

| + 2|ρl||υl|N2

|1 − α g | |υg ||P|N3 + |ρ2g ||υg |N1 |ρ2g |

⎧ ⎫ ⎪ ⎪ υl ⎬, then |f2 | ≤ N2 Letting N2 = sup⎨ N 1 ⎪ ⎪ ⎩ 1 − αg ⎭ ⎫ ⎧ |f | |ρ | |υl2| ⎪ |ρl||g | + | 2 |ld| | + 2|ρl||υl|N2 ⎪ ⎬, then |f3 | ≤ N3 Letting N3 = sup⎨ |1 − α g | ⎪ ⎪ ⎭ ⎩ 2 ⎧ ⎫ ⎪ |υg ||P|N3 + |ρ g ||υg |N1 ⎪ ⎬ , then |f4 | ≤ N4 Letting N4 = sup⎨ ⎪ ⎪ |ρ2g | ⎩ ⎭

⎧ |g | |f3 | |f3 | |υl||f2 | ⎫ |υl||f2 | |g | |a||T − Te| |a||T − Te| ⎬ , then |f5 | ≤ N5 + Letting N5 = sup⎨ + + + + + |ρl||CPl| |CPl| ⎭ |ρl||CPl| |CPl| |ρl||CPl| |CPl| |ρl||CPl| ⎩ |CPl|

6550

⎪

⎪

⎪

⎪



Article

Therefore, || F || ≤ max{N1 , N2 , N3 , N4 , N5}

The partial differentials of f1, f 2, f 3, f4, and f5 about αg,υl,P,υg,T are bounded, as described by the next series of equations. ∂f1

=

∂αg

∂P

=

∂f1

=

∂υg ∂f1 ∂T

d

=

∂υl

−

K4 f ρlυl

∂f1 ∂f1

− K 2(ρ2g + ρg υ2g P − 2υ2g ) + (1 − 2αg )ρ3g g cos θ

K1

K4

−

K42

(4αg ρ2g ρlυl − 8αg υ2g ρlυl + 4P αg ρg ρlυ2g υl)(K1K 2 + K3) K42

K1 + αg ρg g cos θ(1 − αg ) K4 −αg ρg υ2g K1 K4

(− 2ρ2g υ2g + 2ρ2g ρlυl2 − 4υ2g ρlυl2 + 2Pρg ρlυ2g υl2)(K1K 2 + K3)

−

−

(υ2g − ρg υ2g + αg ρlυl2)(K1K 2 + K3) K42

(4ρ2g υg − 2ρ3g υg − 4αg ρ2g υg − 8αg υg ρlυl2 + 4P αg ρg ρlυg υl)(K1K 2 + K3) K42

=0

Then, ∂f1 ∂αg

=

≤

− K 2(ρ2g + ρg υ2g P − 2υ2g ) + (1 − 2αg )ρ3g g cos θ K4 (|ρ2g | + |ρg ||υ2g ||P| + 2|υ2g |)|K 2| + |ρ3g ||g | |K 4 |

+

−

(− 2ρ2g υ2g + 2ρ2g ρlυl2 − 4υ2g ρlυl2 + 2Pρg ρlυ2g υl2)(K1K 2 + K3) K42

(2|ρ2g ||υ2g | + 2|ρ2g ||ρl||υl2|)+ 4|υ2g ||ρl||υl2| + 2P|ρg ||ρl||υ2g ||υl2|) + (|K1||K 2| + |K3|) |K42|

Letting 2 2 3 ⎧ 2 (2|ρ2g ||υ2g | + 2|ρ2g ||ρl||υl2| + 4|υ2g ||ρl||υl2| + 2|P||ρg ||ρl||υ2g ||υl2|)(|K1||K 2| + |K3|) ⎫ ⎪ (|ρ g | + |ρg ||υg ||P| + 2|υg |)|K 2| + |ρg ||g | ⎪ ⎬ M11 = sup⎨ + 2 ⎪ ⎪ |K 4 | |K 4 | ⎭ ⎩

then, ∂f1 ∂αg

∂f1 ∂υl

f ρlυl

=

≤

d

K1

K4

−

(4αg ρ2g ρlυl − 8αg υ2g ρlυl + 4P αg ρg ρlυ2g υl)(K1K 2 + K3)

|K1| + |ρg ||g | |K 4 |

< M11

K42 +

(4|ρ2g ||ρl||υl| + 8|υ2g ||ρl||υl| + 4|P||ρg ||ρl||υ2g ||υl|)(|K1||K 2| + |K3|) |K42|

then, ∂f1 ∂υl

< M12

6551



Article

Letting ⎧ (4|ρ2g ||ρl||υl| + 8|υ2g ||ρl||υl| + 4|P||ρg ||ρl||υ2g ||υl|)(|K1||K 2| + |K3|) ⎫ ⎪ |K1| + |ρg ||g | ⎪ ⎬ M12 = sup⎨ + 2 ⎪ ⎪ K | | K | | 4 4 ⎩ ⎭

then, ((∂f1 )/(∂υl )) ≤ M12

∂f1 ∂P

=

K1 + αg ρg g cos θ(1 − αg ) K4

−

(υ2g − ρg υ2g + αg ρlυl2)(K1K 2 + K3) K42

≤

|K1| + |ρg ||g | |K 4 |

+

(|υ2g | + |ρg ||υ2g | + |ρl||υl2|)(|K1||K 2| + |K3|) |K42|

Letting ⎧ (|υ2g | + |ρg ||υ2g | + |ρl||υl2|)(|K1||K 2| + |K3|) ⎫ ⎪ |K1| + |ρg ||g | ⎪ ⎬ M13 = sup⎨ + 2 ⎪ ⎪ K | | K | | 4 4 ⎩ ⎭

then, ∂f1 ∂P

∂f1 ∂υg

< M13

⎛ −α ρ υ2 K (4ρ2g υg − 2ρ3g υg − 4αg ρ2g υg − 8αg υg ρlυl2 + 4P αg ρg ρlυg υl)(K1K 2 + K3) g g g 1 =⎜ − ⎜ K4 K42 ⎝ ≤

|ρg ||υ2g ||K1| |K 4 |

+

⎞ ⎟ ⎟ ⎠

(4|ρ2g ||υg | + 2|ρ3g ||υg | + 4|ρ2g ||υg | + 8|υg ||ρl||υl2| + 4|P||ρg ||ρl||υg ||υl|)(|K1||K 2| + |K3|) |K42|

Letting ⎧ ⎫ 2 4|ρ2g ||υg | + 2|ρ3g ||υg | + 4|ρ2g ||υg | + 8|υg ||ρl||υl2| + 4|P||ρg ||ρl||υg ||υl|)(|K1||K 2| + |K3|) ⎪ ⎪ |ρg ||υg ||K1| ⎬ M14 = sup⎨ + |K42| ⎪ ⎪ |K 4 | ⎭ ⎩

(

then ∂f1 ∂υg

≤ M14

and ∂f1 ∂T

= 0 ≤ M15

6552



Article

According to a similar method, the following system was obtained: ⎧ ∂f ⎪ 2 ⎪ ∂αg ⎪ ⎪ ⎪ ∂f2 ⎪ ∂υ ⎪ l ⎪ ∂f ⎨ 2 ⎪ ∂P ⎪ ⎪ ∂f ⎪ 2 ⎪ ∂υg ⎪ ⎪ ∂f ⎪ 2 ⎩ ∂T

∂f1 υlf1 υl + 1 − αg ∂αg (1 − αg )2

=

∂f

f1 + υl ∂υ1

l

=

1 − αg

=

υl ∂f1 1 − αg ∂P

=

υl ∂f1 1 − αg ∂υg

=0

Thus, ∂f2 ∂αg

υl

≤

1 − αg

υl N1

M11 +

2

(1 − αg )

Letting ⎧ ⎪ υ l M 21 = sup⎨ M11 + ⎪ 1 − αg ⎩

υl N1 2

(1 − αg )

⎫ ⎪ ⎬ ⎪ ⎭

then ∂f2 ∂αg

≤ M 21 +

∂f2 ∂υl

≤

υl N1 2

(1 − αg )

N1 + υl M12 1 − αg

Letting ⎧ ⎫ ⎪ N + υl M12 ⎪ ⎬ M 22 = sup⎨ 1 ⎪ 1 − αg ⎪ ⎭ ⎩

then ∂f2 ∂υl

≤ M 22

6553



Article

∂f2

≤

∂P

υl M13 1 − αg

Letting ⎧ ⎪ υM M 23 = sup⎨ l 13 ⎪ 1 − αg ⎩

⎫ ⎪ ⎬ ⎪ ⎭

then ∂f2

≤ M 23

∂P

∂f2

υl M14

≤

∂υg

1 − αg

Letting ⎧ ⎪ υM M 24 = sup⎨ l 14 ⎪ 1 − αg ⎩

⎫ ⎪ ⎬ ⎪ ⎭

then ∂f2

≤ M 24

∂υg

and ∂f2 ∂T

= 0 ≤ M 25

The partial differentials of f 3, f4, and f5, about αg, υl, P, υg, and T, may be written as follows, respectively. ⎧ ⎪ ∂f3 ⎪ ⎪ ∂αg ⎪ ⎪ ⎪ ∂f3 ⎪ ∂υ l ⎪ ⎪ ⎨ ∂f ⎪ 3 ⎪ ∂P ⎪ ⎪ ∂f3 ⎪ ⎪ ∂υg ⎪ ⎪ ∂f3 ⎪ ⎩ ∂T

∂f

=

ρlg cos θ + 2ρlυlf2 − 2ρl(1 − αg )υl ∂α2

g

1 − αg −f ρlυl d

=

+

ρlg cos θ(1 − αg ) +

f ρlυl2

+ 2ρl(1 − αg )υlf2

2d

2

(1 − αg )

∂f2

− 2ρl(1 − αg )f2 − 2ρl(1 − αg )υl ∂υ

l

1 − αg

=

− 2ρl(1 − αg )υl ∂f2

=

− 2ρl(1 − αg )υl ∂f2

1 − αg

1 − αg

∂P

∂υg

=0

6554



⎧ ⎪ ∂f4 ⎪ ⎪ ∂αg ⎪ ⎪ ⎪ ∂f4 ⎪ ∂υ ⎪ l ⎪ ⎪ ⎨ ∂f ⎪ 4 ⎪ ∂P ⎪ ⎪ ⎪ ∂f4 ⎪ ∂υ ⎪ g ⎪ ∂f ⎪ 4 ⎪ ⎩ ∂T

Article

∂f

∂f

υg Pf3 + αg υg P ∂α3 + ρ2g υg ∂α1 g

=

g

ρ2g αg ∂f

−

ρ2g (αg υg Pf3 + ρ2g υg f1 ) ρ4g α 2g

∂f

αg υg P ∂υ3 + ρ2g υg ∂υ1 l

=

l

ρ2g αg ∂f

∂f

αg Pf3 + αg υg P ∂υ3 + ρ2g f1 + ρ2g υg ∂υ1 g

=

g

ρ2g αg ∂f

∂f

αg υg f3 + αg υg P ∂P3 + ρ2g υg ∂P1

=

ρ2g αg

=0

and

⎧ ⎪ ∂f5 ⎪ ⎪ ∂αg ⎪ ⎪ ∂f ⎪ 5 ⎪ ∂υl ⎪ ⎪ ⎪ ∂f ⎨ 5 ⎪ ∂P ⎪ ⎪ ⎪ ∂f5 ⎪ ⎪ ∂υg ⎪ ⎪ ∂f ⎪ 5 ⎪ ∂T ⎩

∂f

=

∂f

− ∂α3 − ρlυl ∂α2 g

g

ρlCPl ∂f

=

∂f

− ∂υ3 − ρlf2 − ρlυl ∂υ2 l

l

ρlCPl ∂f

=

∂f

− ∂P3 − ρlυl ∂P2 ∂f

=

=

−

ρlCPl

⎡ρ ∂CPl ⎤[ρ g cos θ − a(T − T ) − f − ρ υ f ] e ⎣ l ∂P ⎦ l l l2 3 ρl2CPl2

∂f

− ∂υ3 − ρlυl ∂υ2 g

g

ρlCPl a−

∂f3 ∂T

∂C

∂f

− ρlυl ∂T2

ρlCPl

−

(

ρl ∂TPl ρlg cos θ − α(T − Te) −

dP dz

dυ

− ρlυl dzl

)

ρl2CPl2

Repeating the above method, the following conditions can be obtained: ∂f3 ∂αg

∂f4 ∂αg

≤ M41 ,

≤ M31,

∂f4 ∂υl

∂f3 ∂υl

≤ M32 ,

≤ M42 ,

∂f3 ∂P

∂f4 ∂P

≤ M33 ,

≤ M43 ,

6555

∂f3 ∂υg

≤ M34 ,

∂f4 ∂υg

∂f3 ∂T

≤ M44 ,

≤ M35

∂f4 ∂T

≤ M45


Industrial & Engineering Chemistry Research ∂f5 ∂αg

≤ M51 ,

∂f5 ∂υl

Article

≤ M52 ,

∂f5 ∂P

∂f5

≤ M53 ,

∂υg

≤ M54 ,

∂f5 ∂T

≤ M55

The Lipschitz condition is very important in discussing the solution of the system of differential equations; thus, the Lipschitz condition of F(z;y) was considered first. The primal problem was written as follows: d αg

= f1 (z ; αg , υl , P , υg , T ), dz d υl = f2 (z ; αg , υl , P , υg , T ), dz dP = f3 (z ; αg , υl , P , υg , T ), dz d υg = f4 (z ; αg , υl , P , υg , T ), dz dT = f5 (z ; αg , υl , P , υg , T ) dz

This can be written as αg ′ = f1 (z ; αg , υl , P , υg , T ), υl′ = f2 (z ; αg , υl , P , υg , T ), P′ = f3 (z ; αg , υl , P , υg , T ), υg ′ = f4 (z ; αg , υl , P , υg , T ), T ′ = f5 (z ; αg , υl , P , υg , T )

The primal condition is αg (z 0) = αg 0 , υl(z 0) = υl 0 , P(z 0) = P0 , υg (z 0) = υg 0 , T (z 0) = T0

Using the Euler method, for i = 0, 1, 2, ..., αg(i + 1) = αgi + (zi + 1 − zi)f1 (zi ; αgi , υli , Pi , υgi , Ti ), υl(i + 1) = υli + (zi + 1 − zi)f2 (zi ; αgi , υli , Pi , υgi , Ti ), Pi + 1 = Pi + (zi + 1 − zi)f3 (zi ; αgi , υli , Pi , υgi , Ti ),

υg(i + 1) = υgi + (zi + 1 − zi)f4 (zi ; αgi , υli , Pi , υgi , Ti ),

Tg(i + 1) = Ti + (zi + 1 − zi)f5 (zi ; αgi , υli , Pi , υgi , Ti )

Here, αgi, υli, Pi, υgi, and Ti are intended to approximate αg(zi), υl(zi), P(zi), υg(zi), and T(zi), where z0 < z1 < z2< ... is a subdivision of the interval of integration. Let yi = (αgi, υli, Pi, υgi, Ti)T; then, yi + 1 = yi + (zi + 1 − zi)F(zi ; yi )

i = 0, 1, 2, ..., n − 1

If set hi = zi+1 − zi, then, for the subdivision above, the following can be written: h = (h0 , h1 , ..., hn − 1)

6556



Article

If we connect y0 and y1, y1 and y2, ..., etc. using straight lines, we obtain the Euler polygon: yh (z) = yi + (z − zi)f (zi ; yi )

for zi ≤ z ≤ zi + 1

Theorem 1. For ∥F(z;y)∥ ≤ N = max{N1,N2,N3,N4,N5}, then, for αgi, υli, Pi, υgi, Ti defined as described above, the estimate

|| yi − y0 || ≤ N |zi − z 0|

where yi = (αgi, υli, Pi, υgi, Ti)T. For ∂fk ∂αg

≤ Mk1 ,

∂fk ∂υl

≤ Mk 2 ,

∂fk ∂P

≤ Mk 3 ,

∂fk ∂υg

≤ Mk 4 ,

∂fk ∂T

≤ Mk 5

then || F(z ; y) − F(z ; y ̂)|| ≤ L|y − y ̂|

where k = 1, 2, 3, 4, 5 and L = maxk(∑5i=1Mki). Proof. (1) From αg(i+1) = αgi + (zi+1 − zi) f1(zi;αgi,υli,Pi,υgi,Ti) and definition of ∥F(z;αg, υl, P, υg, T)∥,

|αg(i + 1) − αgi| = |zi + 1 − zi||f1 (zi ; αgi , υli , Pi , υgi , Ti )| ≤ N (zi + 1 − zi)

Therefore, |αgi − αg(i − 1)| ≤ N (zi − zi − 1), ..., |αg 2 − αg1| ≤ N (z 2 − z1), |αg1 − αg 0| ≤ N (z1 − z 0)

Thus, |αgi − αg(i − 1)| + ... + |αg 2 − αg1| + |αg1 − αg 0| ≤ N (zi − z 0)

Since |αgi − αg(i − 1) + ... + αg 2 − αg1 + αg1 − αg 0| ≤ |αgi − αg(i − 1)| + ... + |αg 2 − αg1| + |αg1 − αg 0|

so, |αgi − αg 0| ≤ N (zi − z 0)

Similarly, |υli − υl0| ≤ N(zi − z0), |Pi − P0| ≤ N(zi − z0), |υgi − υg0| ≤ N(zi − z0), |Ti − T0| ≤ N(zi − z0). From the definition of ∥yi − y0∥, we get

|| yi − y0 || ≤ N (zi − z 0) 6557



Article

(2) For f1(z;y), f 2(z;y), f 3(z;y), f4(z;y), f5(z;y), y = (αg, υl, P, υg, T)T,

f1 (z ; y ̂) − f1 (z ; y) =

∂f1 ∂αg

(α̂ g − αg ) +

∂f1 ∂υl

(υ̂ l − υl) +

∂f1 ∂P

(P ̂ − P ) +

∂f1 ∂υg

(υ̂g − υg ) +

∂f1 ∂T

(T ̂ − T )

Thus,

f1 (z ; y ̂) − f1 (z ; y) ≤

∂f1 ∂αg

α̂ g − αg +

∂f1 ∂υl

∂f1

υ̂ l − υl +

∂P

P̂ − P +

∂f1 ∂υg

υ̂g − υg +

∂f1 ∂T

T̂ − T

Let Δy = max{|α̂ g − αg|,|υ̂l − υl|,|P̂ − P|,|υ̂l − υl||T̂ − T|}, then ⎛ ∂f |f1 (z ; y ̂) − f1 (z ; y)| ≤ ⎜⎜ 1 ⎝ ∂αg

+

⎛ ∂f |f2 (z ; y ̂) − f2 (z ; y)| ≤ ⎜⎜ 2 ⎝ ∂αg

+

⎛ ∂f |f3 (z ; y ̂) − f3 (z ; y)| ≤ ⎜⎜ 3 ⎝ ∂αg

+

⎛ ∂f |f4 (z ; y ̂) − f4 (z ; y)| ≤ ⎜⎜ 4 ⎝ ∂αg

+

⎛ ∂f |f5 (z ; y ̂) − f5 (z ; y)| ≤ ⎜⎜ 5 ⎝ ∂αg

+

∂f1 ∂υl

+

∂f1 ∂P

+

∂f1 ∂υg

+

∂f1 ⎞ ⎟Δy ∂T ⎟⎠

+

∂f2 ⎞ ⎟Δy ∂T ⎟⎠

+

∂f3 ⎞ ⎟Δy ∂T ⎟⎠

+

∂f4 ⎞ ⎟Δy ∂T ⎟⎠

+

∂f5 ⎞ ⎟Δy ∂T ⎟⎠

Similarly, we can get ∂f2 ∂υl ∂f3 ∂υl ∂f4 ∂υl ∂f5 ∂υl

+

∂f2

+

∂f3

∂P

∂P

+

∂f4

+

∂f5

∂P

∂P

+

+

+

+

∂f2 ∂υg ∂f3 ∂υg ∂f4 ∂υg ∂f5 ∂υg

From definition of norm, we have || F(z ; y )̂ − F(z ; y)|| = max{|f1 (z ; y ̂) − f1 (z ; y)| , |f2 (z ; y )̂ − f2 (z ; y)| , |f3 (z ; y ̂) − f3 (z ; y)| , |f4 (z ; y ̂) − f4 (z ; y)| , |f5 (z ; y ̂) − f5 (z ; y)|}

Let L = maxk(∑5i=1Mki) and ∥ŷ − y∥ = Δy, then we get || F(z ; y )̂ − F(z ; y)|| ≤ L || y ̂ − y ||

Therefore, the differential equations system F is continuous and satisfies the Lipschitz conditions. The solution is confirmed and gain the only from the original problem.55 B. The Stability Analysis of Fourth-Order Runge−Kutta Method

According to the solution process above, at every computation point, the model was a linear first-order ordinary differential equation (ODE). Therefore, it cannot be solved analytically; however, in science and engineering, a numeric approximation of the solution is often suitable enough. There are many numerical methods, such as the Euler Method, the Trapezoidal Method, the θ Method, the Adams Method, and the Runge−Kutta Method. The Runge−Kutta method is an important method for the approximation of solutions from ODEs. The fourth-order explicit Runge−Kutta (RK4) method has the following characteristics of generality; less calculation, high velocity, high efficiency, and high precision. However, it has strict conditions for stability. Therefore, stability will be discussed in this section. 6558



Article

For simplicity, equations for the liquid mixture were also used to analyze the stability. The system can be translated to the following form: D

dU = FU + G dz

where ⎡ vg αg P ρ2g υg ρ2g υg 0 0 ⎤ ⎢ ⎥ ⎢ −υl 0 1 − αg 0 ⎥ 0 ⎢ ⎥ ⎢ 0 0 2αlρlυl 0 ⎥ D = ⎢ 1 − αg ⎥ ⎢ ⎥ 2 2 2 2 2 0 0 ⎥ ⎢ αg (υg P + ρg ) ρg υg ρg υg ⎢ ⎥ υlρl ρlCPl ⎦ 1 0 0 ⎣

⎡0 0 ⎢ 0 ⎢0 ⎢ 0 g cos θρ l F=⎢ ⎢ 0 −g ρ2 cos θ g ⎢ ⎢0 0 ⎣

⎤ ⎡ 0 ⎥ ⎢ 0 ⎥ ⎢ G = ⎢ −g cos θρl ⎥ ⎥ ⎢ 0 ⎥ ⎢ ⎥ ⎢ g cos θ + a ρ T ⎣ l e⎦

0 ⎤ ⎥ 0 ⎥ 0 f ρl /(2d) 0 ⎥ ⎥ 0 0 0 ⎥ ⎥ −aρl ⎥⎦ 0 0 0 0

0 0

⎡P⎤ ⎢α ⎥ ⎢ g⎥ U = ⎢ υg ⎥ ⎢ ⎥ ⎢ υl ⎥ ⎢ ⎥ ⎣T ⎦

Considering the characteristic equation, the mathematical characters for a set of ODEs can be found from the solving of the following eigenvalue system: det[λD − F] = 0

λ1 = λ 2 = λ 3 = 0,

λ5 =

λ4 = −

a CPl

2ρ2g g cos θd(1 − 2αg ) + 2α 2g ρg g cos θd(ρg − ρl) + αg ρg f υl(ρg − P υ2g ) + 2αg ρlg cos θd(ρg − P υ2g + P αg υ2g ) 2dρg [ρg υ2g (αg − 1) − 2α 2g ρlυl2] + 4Pρld αg (αg − 1)υ2g υl2 + 4d αg ρg (ρlυl2 − ρg υ2g )

In engineering practice, the values of λ4 and λ5 can be found (λ4 < 0 and λ5 > 0). Therefore, the differential system solution tended toward divergence. However, the value of the positive eigenvalue is very small (in the range of 0.009−0.1). For a solution cycle, the approximate solution is weakly influenced. These models belong to a well-conditioned model. The absolute value of the negative characteristic parameter range is very small but has an important impact on the RK4 methods. For ODEs, if the numerical results are stable with a fixed grid size h, the product of λ and h must fall into an absolutely stable interval. The RK4 interval is [−2.78,0].55 Therefore, the system is stable at grid h < 28. C. Calculation of Some Parameters and the Initial Conditions

In order to solve model eq 29, some definite conditions (i.e., initial conditions) must be added. The initial conditions comprise the pressure distribution, temperature, and velocity along the wall at the top of the well. In this section, the calculation method for some parameters is given. (1) Each point’s inclination

θj = θj − 1 + (θk − θk − 1)

6559

Δsj Δsk



Article

where j represents the calculation segment point of calculation, Δsk represents the measurement depth between inclination angle θk and θk−1, and Δsj represents the step length of calculation. (2) Transient heat transfer f unction (Hasan et al.38) ⎧1.128 t (1 − 0.3 t ) tD ≤ 1.5 D D ⎪ ⎪ f (tD) = ⎨ ⎛ 0.6 ⎞ ⎪(0.4063 + 0.5 ln tD)⎜1 + ⎟ tD > 1.5 ⎪ tD ⎠ ⎝ ⎩

(3) Frictional force between the liquid phase and the pipe wall. In the present model, the gas−wall friction and the oil−wall friction are neglected, and the friction force is represented by the mixture−wall interaction, which is given by56 f ρ υk2 τkbSkb = k k A 2d

Because the friction factor may be affected by the Reynolds numbers of fluid in pipeline,the Blasius equation should be adopted. Many researchers have used this equation to compute the frictional force between the oil (or water) and the pipe wall. ⎧ Re Re ≤ 2000 ⎪ ⎪ 64 fk = ⎨ ⎪ 0.3164 ⎪ 0.25 Re > 1.5 ⎩ Re

Rek =

ρkυk d μk

(4) Overall heat-transfer coef f icient. Under the same assumptions,30 the overall heat-transfer coefficient can be obtained using the following equation:

1 Uto−1 = + hc + hr

( ) + r ln( )

rti ln

rcem rco

kcem

ti

rci rto

kang

(5) Gas condensing parameter: ⎧ ⎛ ⎞ ⎞ ⎛ ρ2 ⎪1 + ⎜0.31506 − 1.0467 − 0.5783 ⎟ρ + ⎜0.053 − 0.6123 ⎟ρ2 + 0.6815 pr if P < 35 MPa ⎜ ⎟ pr ⎜ 3 3 ⎪ Tpr Tpr ⎟⎠ pr T pr T pr ⎝ ⎝ ⎠ ⎪ ⎪ Zg = ⎨ ⎪(90.7x − 242.2x 2 + 42.4x 3)y1.18 + 2.82x − (14.76x − 9.76x 2 + 4.58x 3)y otherwise ⎪ (1 + y + y 2 + y 3 ) ⎪ + ⎪ (1 − y)3 ⎩

F(y) = − 0.06125Pprx exp[− 1.2(1 − x)2 ] + (90.7x − 242.2x 2 + 42.4x 3)y 2.18 + 2.82x +

1 + y + y2 + y3 − (14.76x − 9.76x 2 + 4.58x 3)y 2 = 0 (1 − y)3

x=

1 Tpr

(6) Calculation of the conditions of the gas velocity, water velocity, oil velocity, oil cut, water cut, pressure and temperature at the initial point. Assume that the pressure and temperature at the first computation point are already known. Then, ρg0 =

MP0 RZ0T0

6560



Article

ρ = density of flowing liquid (kg/m3) α = volume fraction

D. SI Metric Conversion Factors

Factors used to convert experimental values to SI units are given in Table D1.

Subscripts and Superscripts

Table D1. SI Metric Conversion Factorsa non-SI unit Btu °F ft lbm lbf psi ft2 eV hp a

factor × × × × × × × ×

1.055056 (°F − 32)/1.8 3.048 × 10−1 4.535924 × 10−1 4.44822 6.894757 9.290304 × 10−2 1.60219 × 10−19 7.46043 × 10−1

SI unit = = = = = = = = =

kJ °C m kg N kPa m2 J kW

■

i = grid discretization index in the z-direction m = mixture phase flow (gas−liquid) w = water flow phase k = any phase (mixture, gas or water flow) g = gas phase flow o = oil phase flow

AUTHOR INFORMATION

Corresponding Author

*Tel.: +86 28 85418522. Fax: +86 28 85400222. E-mail: [email protected]. Notes

Conversion factors are exact.

The authors declare no competing financial interest.

■

NOMENCLATURE A = cross-sectional area flow conduit (m2) Cp = heat capacity of fluids (J/kg) Cf = specific heat of formation (J/(kg J)) d = hydraulic diameter (m) f tD = dimensionless transient heat conduction time function for formation rD = dimensionless radius (dimensionless) h = enthalpy (J2/s2) g = acceleration constant of gravity (m/s2) tD = dimensionless time Re = Reynolds number Ke = formation conductivity (J/(m K)) rti = inner radius of conduit (m) T = temperature of flowing liquid (K) Tw = temperature of the wellbore (K) ρf = mass density of surrounding formation (kg/m3) P = pressure of flowing liquid (Pa) Ppc = critical pressure (Pa) Ppr = comparative pressure (dimensionless) ρ = density of flowing liquid (kg/m3) q = heat flux (J/(m s)) M = molar mass of the gas rto = outer radius of conduit (m) R = gas constant Si = wetted perimeter at depth i (m) Tr = temperature of the formation/earth (K) Te = initial temperature of formation (K) Tk = second interface temperature (K) Tpc = critical temperature (K) Tref = second interface temperature (K) Tpr = comparative temperature (dimensionless) t = production time of gas (s) τi = shearing friction force between flow and wall (N) Uto = overall-heat-transfer coefficient (W/(m K)) Z = total length of flow conduit (m) Zg = gas compressibility factor z = distance coordinate in the flow direction along the conduit (m)

■

ACKNOWLEDGMENTS

■

REFERENCES

This research was supported by the Key Program of NSFC (Grant No. 70831005), and the Key Project of China Petroleum and Chemical Corporation (Grant No. GJ-73-0706).

(1) Oglesby, K. An experimental study on the effect of oil viscosity, mixture, velocity and water fraction on horizontal oil-water flow. M.Sc. Thesis, University of Tulsa, Tulsa, OK, 1979. (2) Trallero, J. Oil−water flow patterns in horizontal pipes. Ph.D. Thesis, University of Tulsa, Tulsa, OK, 1995. (3) Trallero, J.; Sarica, C.; Brill, J. A study of oil−water flow patterns in horizontal pipes. SPE Prod. Facil. 1996 12 (3), 165−172. (4) Nadler, M.; Mewes, D. Flow induced emulsification in the flow of two immiscible liquids in horizontal pipes. Int. J. Multiphase Flow 1997, 23, 55−68. (5) Angeli, P.; Hewitt, G. Pressure gradient in horizontal liquid− liquid flows. Int. J. Multiphase Flow 1999, 24, 1183−1203. (6) Shi, H.; Cai, J.; Jepson, W. The effect of surfactants on flow characteristics in oil/water flows in large diameter horizontal pipelines. Presented at the 1999 Multiphase Technology Conference, Cannes, France, June 1999. (7) Mu, H. Experimental research on oil-water horizontal pipe flow. Ph.D. Thesis, University of Petroleum, 2001. (8) Shi, H. A study of oil−water flows in large diameter horizontal pipelines. Ph.D. thesis, University of Petroleum, 2001. (9) Shi, H.; Holmes, J.; Diaz, L.; Durlofsky, L.; Aziz, K. Drift−flux parameters for three-phase steady-state flow in wellbores. SPE J. 2005, 10, 130−137. (10) Cazarez, O.; Montoya, D.; Vital, A.; Bannwart, A. Modeling of three-phase heavy oil−water−gas bubbly flow in upward vertical pipes. Int. J. Multiphase Flow 2010, 36, 439−448. (11) Açikgöz, M.; Franca, F.; Lahey, R., Jr. An experimental study of three-phase flow regimes. Int. J. Multiphase Flow 1992, 18, 327−336. (12) Lahey, R., Jr.; Açikgöz, M.; Franca, F. Global volumetric phase fractions in horizontal three-phase flows. AlChE J. 1992, 38, 1049− 1058. (13) Pan, L.; Jayanti, S.; Hewitt, G. Flow patterns, phase inversion and pressure gradient in air−oil−water flow in a horizontal pipe. Presented at the 2nd International Conference on Multiphase, Kyoto, Japan, 1995. (14) Woods, G.; Spedding, P.; Watterson, J.; Raghunathan, R. Threephase oil/water/air vertical flow. Chem. Eng. Res. Des. 1998, 76, 571− 584. (15) Poettman, F.; Carpenter, P. The multiphase flow of gas, oil, and water through vertical flow strings with application to the design of gas-lift installations. Drill Prod. Pract. 1952, 257−317.

Greek Symbols

θ = inclination angle flow conduit (°) v = velocity of flowing liquid (m/s) 6561



Article

(16) Duns, H.; Ros, N. Vertical Flow of Gas and Liquid Mixtures in Well. In Proceedings of the Sixth World Petroleum Congress, June 19−26, 1963; p 451. (17) Hagedorn, A.; Brown, K. Experimental study of pressure gradients occurring during continuous two-phase flow in smalldiameter vertical conduits. J. Pet. Technol. 1965, 17, 475−484. (18) Orkiszewski, J. Predicting two-phase pressure drops in vertical pipe. J. Pet. Technol. 1967, 19, 829−838. (19) Aziz, K.; Govier, G.; Govier, G. Pressure drop in wells producing oil and gas. J. Can. Pet. Technol. 1972, 11. (20) Beggs, D.; Brill, J. A study of two-phase flow in inclined pipes. J. Pet. Technol. 1973, 25, 607−617. (21) Taitel, Y.; Bornea, D.; Dukler, A. Modelling flow pattern transitions for steady upward gas−liquid flow in vertical tubes. AIChE J. 1980, 26, 345−354. (22) Markatos, N.; Singhal, A. Numerical analysis of one-dimensional, two-phase flow in a vertical cylindrical passage. Adv. Eng. Software 1982, 4, 99−106. (23) Nenes, A.; Assimacopoulos, D.; Markatos, N.; Mitsoulis, E. Simulation of airlift pumps for deep water wells. Can. J. Chem. Eng. 1996, 74, 448−456. (24) Latsa, M.; Assimacopoulos, D.; Stamou, A.; Markatos, N. Twophase modeling of batch sedimentation. Appl. Math. Modell. 1999, 23, 881−897. (25) Markatos, N.; Kirkcaldy, D. Analysis and computation of threedimensional, transient flow and combustion through granulated propellants. Int. J. Heat Mass Transfer 1983, 26, 1037−1053. (26) Markatos, N.; Pericleous, K. Laminar and turbulent natural convection in an enclosed cavity. Int. J. Heat Mass Transfer 1984, 27, 755−772. (27) Markatos, N. Modelling of two-phase transient flow and combustion of granular propellants. Int. J. Multiphase Flow 1986, 12, 913−933. (28) Hoffman, N.; Galea, E.; Markatos, N. Mathematical modelling of fire sprinkler systems. Appl. Math. Modell. 1989, 13, 298−306. (29) Xiao, Z. The calculation of oil temperature in a well. SEP 1987, 17125, 1−14. (30) Sagar, R.; Doty, D.; Schmidt, Z. Predicting temperature profiles in a flowing well. SPE Prod. Eng. 1991, 6, 441−448. (31) Alves, I.; Alhanati, F.; Shoham, O. A Unified Model for Predicting Flowing Temperature Distribution. in Wellbores and Pipelines. SPE Prod. Eng. 1992, 7, 363−367. (32) Kirkpatrick, C. Advances in gas-lift technology. Drill. Prod. Prac., API 1959, 24−60. (33) Hagoort, J. Ramey’s wellbore heat transmission revisited. SPE J. 2004, 9, 465−474. (34) Vielma, M. A.; Atmaca, S.; Sarica, C.; Zhang, H.-Q. Characterization of oil/water flows in horizontal pipes. SPE Proj. Facil. Construct. 2007, 3 (4), 1−21. (35) Kabir, C.; Hasan, A.; Jordan, D.; Wang, X. A wellbore/reservoir simulator for testing gas wells in high-temperature reservoirs. SPE Form. Eval. 1996, 11, 128−134. (36) Hasan, A.; Kabir, C.; Wang, X. Development and application of a wellbore/reservoir simulator for testing oil wells. SPE Form. Eval. 1997, 12, 182−188. (37) Hasan, A.; Kabir, C.; Wang, X. Wellbore Two-Phase Flow and Heat Transfer During Transient Testing. SPE J. 1998, 3, 174−180. (38) Hasan, A.; Kabir, C.; Sarica, C. Fluid Flow and Heat Transfer in Wellbores; Society of Petroleum Engineers, 2002. (39) Hasan, A.; Kabir, C.; Lin, D. Analytic Wellbore-Temperature Model for Transient Gas−Well Testing. SPE Reservoir Eval. Eng. 2005, 8 (1), 240−247. (40) Hasan, A.; Kabir, C.; Sayarpour, M. Simplified two-phase flow modeling in wellbores. J. Pet. Sci. Technol. 2010, 72, 42−49. (41) Cazarez-Candia, O.; Vasquez-Cruz, M. Prediction of pressure, temperature, and velocity distribution of two-phase flow in oil wells. J. Pet. Sci. Technol. 2005, 46, 195−208.

(42) Xu, J.; Wu, Z.; Wang, S.; Qi, B.; Chen, K.; Li, X.; Zhao, X. Prediction of temperature, pressure, density, velocity distribution in HT-H-P gas wells. Can. J. Chem. Eng. 2011, DOI: 10.1002/cjce.20689. (43) Zhang, H.; Sarica, C. Unified modeling of gas/oil/water pipe flowBasic approaches and preliminary validation. Presented at the SPE Annual Technical Conference and Exhibition, Tulsa, OK, 2005. (44) Bonizzi, M.; Issa, R. On the simulation of three-phase slug flow in nearly horizontal pipes using the multi-fluid model. Int. J. Multiphase Flow 2003, 29, 1719−1747. (45) Ghorai, S.; Suri, V.; Nigam, K. Numerical modeling of three-phase stratified flow in pipes. Chem. Eng. Sci. 2005, 60, 6637− 6648. (46) Ramey, J. R.; Mobil, H.; Wellbore, O. Heat Transmission. J. Pet. Technol. 1962, 14, 427−435. (47) Willhite, G. Over-all heat transfer coefficients in steam and hot water injection wells. J. Pet. Technol. 1967, 19, 607−615. (48) Bertin, J. J.; Smith, M. L. Aerodynamics for Engineers; Prentice− Hall: Englewood Cliffs, NJ, 1979. (49) Jiang, Y.; Li, Y.; Feng, S. Numerical simulation of transient twophase flow. Oil Gas Storage 2005, 24, 22−24 (in Chin.). (50) Dekker, K.; Verwer, J. Stability of Runge−Kutta Methods for Stiff Nonlinear Differential Equations; North−Holland: Amsterdam, 1984; Vol. 984. (51) Chang, L. Singular Value Decomposition Method and Its Application for Solution of Linear Pathological Equations; Geological Publishing House: Beijing, 1990 (in Chin.). (52) Mitra-Majumdar, D.; Farouk, B.; Shah, Y.; Macken, N.; Oh, Y. Two- and three-phase flows in bubble columns: Numerical predictions and measurements. Ind. Eng. Chem. Res. 1998, 37, 2284− 2292. (53) De Henau, V.; Raithby, G. A transient two-fluid model for the simulation of slug flow in pipelines−II. Validation. Int. J. Multiphase Flow 1995, 21, 351−363. (54) Grolman, E.; Fortuin, J. Gas−liquid flow in slightly inclined pipes. Chem. Eng. Sci. 1997, 52, 4461−4471. (55) Birkhoff, G.; Rota, G.-C. Ordinary Differential Equations; Wiley: New York, 1978. (56) Taitel, Y.; Barnea, D.; Brill, J. Stratified three phase flow in pipes. Int. J. Multiphase Flow 1995, 21, 53−60.

6562


Prediction of Pressure, Temperature, Holdup, Velocity, and Density

Recommend Documents