# Prediction of Pressure, Temperature, Holdup ... - ACS Publications

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 present...
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 diﬀerential equations (PDEs) concerning pressure, temperature, velocity, and holdup for a gas−oil−water three-phase steady ﬂow 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 diﬀerent 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 eﬀectiveness 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 ﬂow and gas−oil−water three-phase ﬂow 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 ﬂow often occurs in onshore and oﬀshore hydrocarbon production and transportation. Studies on the oil−water two-phase ﬂows and oil−gas−water three-phase ﬂows could solve many important technical problems in the petroleum industry. Studies on the oil−gas−water three-phase ﬂow, especially on the temperature eﬀects and pressure gradients, can solve many important technical problems for the petroleum industry and could also be very important for improvements in multiphase ﬂow 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 signiﬁcantly diﬀerent when outputs are varied (ﬂow velocity) but neither has a simple linear relationship, because the ﬂuid density is not constant. Therefore, for HTHP wells, the prediction of three-phase ﬂows is important to the industry. In this paper, the three-phase ﬂow is treated as a bubbly gas ﬂow. 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 ﬂow patterns had been observed while several researchers described three or four.1 In the past 30 years, there has been signiﬁcant progress in the analysis of two-phase ﬂow patterns. Many new and comprehensive ﬂow patterns have been published.2−9 Because two-phase gas−liquid ﬂows are highly complex, it is apparent that the addition of a third phase will increase this complexity.10 Experimental observations have shown that the ﬂow structures of a three-phase-pipe ﬂow are much more complicated than that of the two-phase-pipe ﬂow. For example, 10 ﬂow patterns were observed by the experiment study.11,12 Seven (7) ﬂow patterns were identiﬁed 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

ﬂow.13 For vertical air−water−oil ﬂow, 8 ﬂow patterns were identiﬁed.14 A two-phase ﬂow is the most common ﬂow for ﬂuid in nature, which widely exists in developed oil wells in the middle and later stages. The understanding, description, and prediction of ﬂow characteristics has become a research focus in the ﬁeld of two-phase ﬂows. Numerical simulation methods have been found to be an eﬀective way to study two-phase ﬂows. A twophase ﬂow 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 ﬂow 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 ﬂow 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 ﬂow-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 ﬁxed 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 ﬂows. A two-phase ﬂow pressure drop computation model covering all patterns in a vertical tube was established and the ﬂow pattern identiﬁcation 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 ﬂow and slug ﬂow focusing on a gas−liquid two-phase ﬂow, where a gas volume factor was introduced into the density and friction losses using a gas− liquid two-phase separating eﬀect.19 The relationship of a liquid hold-up and drag coeﬃcient based on a homogeneous ﬂow pressure gradient equation based on homogeneous ﬂow was experimentally studied.20 Beggs−Brill’s correlations were built and could be used for almost any pipeline inclination. The ﬂow pattern changing conditions from the mechanism was explained and a physical model was proposed.21 Their model had an important impact on two-phase ﬂow 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, ﬂow analysis of steam and water mixtures in vertical ﬂow passages in a general form. They extended the model to applications such as ﬁre 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 ﬂows 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 ﬂow was experimentally studied and found that a droplet size distribution will directly aﬀect the accuracy of

model predictions.34 Hasan et al. made good jobs to the coupled ﬂuid ﬂow model where wellbore/reservoir simulators for modeling single-phase gas, oil and the two-phase gas−oil ﬂow problems were presented. The transient ﬂuid and heat ﬂows models are solved numerically using ﬁnitediﬀerence 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 ﬂow pattern deﬁnitions and some forces are diﬀerent in each pattern. More ﬂow patterns mean more discontinuities and greater complexity in the models. Since there are many papers identifying the patterns of ﬂow, 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 ﬂow 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 ﬂow is to combine oil and water into a single liquid phase. Gas−liquid−liquid three-phase ﬂows can be regarded as a special type of gas−liquid two-phase ﬂow if the two liquids are fully mixed. This is probably true of vertical and steeply inclined ﬂows. 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 ﬂow as a three-layer stratiﬁed ﬂow with gas on the top, oil in the middle, and water at the bottom. This can be done for immiscible liquids ﬂowing in horizontal or slightly inclined pipes with low gas, oil, and water ﬂow 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 stratiﬁed ﬂow in a horizontal pipeline.45 The concept of extended velocity was applied to compute the wall shear stresses. However, the temperature proﬁle was not considered. A three-phase (heavy-oil−water− gas) bubbly ﬂow in upward pipes was simulated using a onedimensional transient two-ﬂuid 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 proﬁles rather than the temperature proﬁles and failed to consider heat transfer.

MODEL BUILDING Problem Description. Within the tubing, ﬂow takes place under turbulent ﬂow conditions. Consider the ﬂow system depicted in Figure 1: a straight cylindrical ﬂow tube with an inclination angle θ, a constant cross-sectional ﬂow area A, a hydraulic diameter d, and a total length Z. Gas ﬂows through the tubing from the bottom to the top with a mass ﬂow rate W. The distance coordinate in the ﬂow 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 ﬂow, 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

Industrial & Engineering Chemistry Research

Article

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

=0

(3)

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

Figure 1. Schematic of the physical ﬁgure.

− 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 ﬂows in the opposite direction: ρkαkg(cos θ)A dz • The friction which ﬂows in the opposite direction: τkbSkb dz • The interfacial shear stress, which ﬂows in the opposite direction: τkjSkj dz The momentum ﬂowing out from dz is given as ρkαkυ2k + A(d(ρkαkυ2k)/dz). The momentum ﬂowing 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 diﬀerent from phase k. From the law of momentum conservation, the momentum equation for each phase is as follows: Water ﬂow momentum equation:

Basic Assumption. Considering the diﬀerential equation model of P, T, v and a, the following assumptions were set. (1) The gas−liquid−liquid ﬂow in the tubing is in one dimension of the ﬂowing 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 ﬂow is with a ﬂat 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 ﬂow 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 ﬂow momentum equation: d(ρoαoυo2) dz

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

Model Formulation. In the paper, a three-phase system in a inﬁnitesimal section is considered as shown in Figure 2. Referring to Figure 2, the ﬂow of three ﬂuids is considered. It is assumed that water is heavier than oil and ﬂows at the bottom. The oil ﬂows 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 diﬀerential depth. Thus, under steady-state conditions, for each phase, it follows the law of ﬂuid 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 ﬂowing 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 ﬂows 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 uniﬁed energy equation is given as

Water ﬂows continuity equation: =0

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

Figure 2. Bubbly gas three phase ﬂow.

d(ρωαωνωA)

= −αo

(2) 6539

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

Industrial & Engineering Chemistry Research

Article

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

The energy ﬂowing out from dz includes the following types:

The term ηk is the Joule−Thomson coeﬃcient, which is deﬁned 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 ﬂuid ﬂowing into the distance element equals the energy sum of losses and ﬂuid ﬂowing 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 ﬂuid and the earth have been discussed in detail.46,47 As shown in Figure 3, the radial heat transfer from the ﬂuid 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 satisﬁes the following relation, denotes the speciﬁc 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

Industrial & Engineering Chemistry Research

Article

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

For one thing, the number of diﬀerential equations for a three-ﬂuid model is twice that for a two-ﬂuid model. For another, the nonlinearity and coupling of the equations are much stronger, because of the interactions between the three ﬂuids. The equations are coupled not only within those for the same phase, but also among those for diﬀerent phases. Because of the coupling of equations, the work to solve the equations directly is very diﬃcult. In this paper, a simpliﬁed 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 ﬂows. 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 diﬀerential 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-ﬂuid model system is obviously more complex to compute than a one-ﬂuid or two-ﬂuid system. 6541

(27)

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

Industrial & Engineering Chemistry Research

Article

Combining eqs 15, 16, 21−23, and 27, the coupled system of diﬀerential 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 modiﬁed 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 coeﬃcient 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

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

Industrial & Engineering Chemistry Research

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 ﬂow boundary conditions are deﬁned 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 ﬂow 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 diﬀerential 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 diﬀerential 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 deﬁned 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, ﬂuid 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

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

Industrial & Engineering Chemistry Research

Article

Table 1. Pipe Parameters diameter (mm)

thickness (mm)

weight

expansion

coeﬃcient

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 speciﬁc 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 eﬀective method that must be analyzed more deeply.10,11,43 To study the suitability of the model, some parameters were adjusted to reﬂect totally diﬀerent 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 ﬂuid temperature and tubing pressure heavily aﬀect the force of the tubular string. Therefore, because, in such an experiment, it is diﬃcult 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 eﬀect of grid ﬁneness 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

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

Industrial & Engineering Chemistry Research

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 proﬁles 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 ﬁxed production. The reason is that the ﬂowing gas phase tends to overcome the gravity and friction. Because the pressure drop is inﬂuenced 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

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

Industrial & Engineering Chemistry Research

Article

Figure 7. Gas production sensitivity analysis.

6546

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

Industrial & Engineering Chemistry Research

Article

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

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

Industrial & Engineering Chemistry Research

Article

Figure 9. Comparison analysis.

fraction is inﬂuenced by the liquid velocity and the decrease in liquid velocity. The eﬀect of an increase in velocity at high gas production is almost balanced out by the increase in the velocity drop eﬀect. 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 eﬀect 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 eﬀect on temperature and the diﬀerence 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 aﬀected the velocity and volume fraction proﬁles shown in Figures 8c−g. Figure 8h shows the gas density proﬁles. At the same depth, the heat transfer had a weak eﬀect on gas density (from 149.1 kg/m3 to 147.7kg/m3). This is because the gas density was inﬂuenced not only by temperature but also by other parameters, according to the density equations. Therefore ,the main inﬂuence on gas density is gas production. The reason why the heat transfer appeared to only aﬀect 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 proﬁles, an energy equation is not necessary; these results reﬂect 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 inﬂuenced the result were taken into consideration. The ﬁrst was gas output with three diﬀerent measures used: 300 000 m3/day, 500 000 m3/day, and 700 000 m3/day. In order to study the diﬀerence 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 Diﬀerent 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 proﬁles with variable gas production. At the same depth, pressure decreases with an increase in production. This is because ﬂowing 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 proﬁles 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

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

Industrial & Engineering Chemistry Research

Article

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 ﬁeld data for a three-phase pipeline ﬂow, a homogeneous model was developed for the pressure prediction, temperature, and velocity distribution based on a two-phase ﬂow 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 diﬀerential 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 diﬀerent parameters such as gas production or heat transfer, are presented. It is physically obvious that the gas production is the main inﬂuential factor. (3) Further investigations, for example with diﬀerent ﬂow-regime and transient, are desirable. (The volume fraction change was so small found from the model and it can be simpliﬁed to ﬁt for large equations in transient ﬂow 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 modiﬁed 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

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

Industrial & Engineering Chemistry Research

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

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

Industrial & Engineering Chemistry Research

Article

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

The partial diﬀerentials 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

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

Industrial & Engineering Chemistry Research

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

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

Industrial & Engineering Chemistry Research

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

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

Industrial & Engineering Chemistry Research

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 diﬀerentials 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

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

Industrial & Engineering Chemistry Research

⎧ ⎪ ∂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

− ρ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

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

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 diﬀerential equations; thus, the Lipschitz condition of F(z;y) was considered ﬁrst. 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

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

Industrial & Engineering Chemistry Research

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 deﬁned 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 deﬁnition 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 deﬁnition of ∥yi − y0∥, we get

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

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

Industrial & Engineering Chemistry Research

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 deﬁnition 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 diﬀerential equations system F is continuous and satisﬁes the Lipschitz conditions. The solution is conﬁrmed 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 ﬁrst-order ordinary diﬀerential 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 eﬃciency, and high precision. However, it has strict conditions for stability. Therefore, stability will be discussed in this section. 6558

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

Industrial & Engineering Chemistry Research

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 diﬀerential 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 inﬂuenced. 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 ﬁxed 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 deﬁnite 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

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

Industrial & Engineering Chemistry Research

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 aﬀected by the Reynolds numbers of ﬂuid 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 coeﬃcient 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 ﬁrst computation point are already known. Then, ρg0 =

MP0 RZ0T0

6560

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

Industrial & Engineering Chemistry Research

Article

ρ = density of ﬂowing 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 ﬂow (gas−liquid) w = water ﬂow phase k = any phase (mixture, gas or water ﬂow) g = gas phase ﬂow o = oil phase ﬂow

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 ﬁnancial interest.

NOMENCLATURE A = cross-sectional area ﬂow conduit (m2) Cp = heat capacity of ﬂuids (J/kg) Cf = speciﬁc 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 ﬂowing liquid (K) Tw = temperature of the wellbore (K) ρf = mass density of surrounding formation (kg/m3) P = pressure of ﬂowing liquid (Pa) Ppc = critical pressure (Pa) Ppr = comparative pressure (dimensionless) ρ = density of ﬂowing liquid (kg/m3) q = heat ﬂux (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 ﬂow and wall (N) Uto = overall-heat-transfer coeﬃcient (W/(m K)) Z = total length of ﬂow conduit (m) Zg = gas compressibility factor z = distance coordinate in the ﬂow 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 eﬀect of oil viscosity, mixture, velocity and water fraction on horizontal oil-water ﬂow. M.Sc. Thesis, University of Tulsa, Tulsa, OK, 1979. (2) Trallero, J. Oil−water ﬂow 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 ﬂow patterns in horizontal pipes. SPE Prod. Facil. 1996 12 (3), 165−172. (4) Nadler, M.; Mewes, D. Flow induced emulsiﬁcation in the ﬂow 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 ﬂows. Int. J. Multiphase Flow 1999, 24, 1183−1203. (6) Shi, H.; Cai, J.; Jepson, W. The eﬀect of surfactants on ﬂow characteristics in oil/water ﬂows 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 ﬂow. Ph.D. Thesis, University of Petroleum, 2001. (8) Shi, H. A study of oil−water ﬂows in large diameter horizontal pipelines. Ph.D. thesis, University of Petroleum, 2001. (9) Shi, H.; Holmes, J.; Diaz, L.; Durlofsky, L.; Aziz, K. Drift−ﬂux parameters for three-phase steady-state ﬂow 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 ﬂow 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 ﬂow 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 ﬂows. AlChE J. 1992, 38, 1049− 1058. (13) Pan, L.; Jayanti, S.; Hewitt, G. Flow patterns, phase inversion and pressure gradient in air−oil−water ﬂow 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 ﬂow. Chem. Eng. Res. Des. 1998, 76, 571− 584. (15) Poettman, F.; Carpenter, P. The multiphase ﬂow of gas, oil, and water through vertical ﬂow strings with application to the design of gas-lift installations. Drill Prod. Pract. 1952, 257−317.

Greek Symbols

θ = inclination angle ﬂow conduit (°) v = velocity of ﬂowing liquid (m/s) 6561

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

Industrial & Engineering Chemistry Research

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 ﬂow 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 ﬂow in inclined pipes. J. Pet. Technol. 1973, 25, 607−617. (21) Taitel, Y.; Bornea, D.; Dukler, A. Modelling ﬂow pattern transitions for steady upward gas−liquid ﬂow in vertical tubes. AIChE J. 1980, 26, 345−354. (22) Markatos, N.; Singhal, A. Numerical analysis of one-dimensional, two-phase ﬂow 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 ﬂow 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 ﬂow and combustion of granular propellants. Int. J. Multiphase Flow 1986, 12, 913−933. (28) Hoﬀman, N.; Galea, E.; Markatos, N. Mathematical modelling of ﬁre 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 proﬁles in a ﬂowing well. SPE Prod. Eng. 1991, 6, 441−448. (31) Alves, I.; Alhanati, F.; Shoham, O. A Uniﬁed 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 ﬂows 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. Simpliﬁed two-phase ﬂow 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 ﬂow 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. Uniﬁed modeling of gas/oil/water pipe ﬂow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 ﬂow in nearly horizontal pipes using the multi-ﬂuid model. Int. J. Multiphase Flow 2003, 29, 1719−1747. (45) Ghorai, S.; Suri, V.; Nigam, K. Numerical modeling of three-phase stratiﬁed ﬂow 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 coeﬃcients 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 Cliﬀs, NJ, 1979. (49) Jiang, Y.; Li, Y.; Feng, S. Numerical simulation of transient twophase ﬂow. Oil Gas Storage 2005, 24, 22−24 (in Chin.). (50) Dekker, K.; Verwer, J. Stability of Runge−Kutta Methods for Stiﬀ Nonlinear Diﬀerential 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 ﬂows in bubble columns: Numerical predictions and measurements. Ind. Eng. Chem. Res. 1998, 37, 2284− 2292. (53) De Henau, V.; Raithby, G. A transient two-ﬂuid model for the simulation of slug ﬂow in pipelines−II. Validation. Int. J. Multiphase Flow 1995, 21, 351−363. (54) Grolman, E.; Fortuin, J. Gas−liquid ﬂow in slightly inclined pipes. Chem. Eng. Sci. 1997, 52, 4461−4471. (55) Birkhoﬀ, G.; Rota, G.-C. Ordinary Diﬀerential Equations; Wiley: New York, 1978. (56) Taitel, Y.; Barnea, D.; Brill, J. Stratiﬁed three phase ﬂow in pipes. Int. J. Multiphase Flow 1995, 21, 53−60.

6562

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