A mechanistic model of heat transfer for gas–liquid flow in vertical wellbore annuli

The most prominent aspect of multiphase flow is the variation in the physical distribution of the phases in the flow conduit known as the flow pattern. Several different flow patterns can exist under different flow conditions which have significant effects on liquid holdup, pressure gradient and heat transfer. Gas–liquid two-phase flow in an annulus can be found in a variety of practical situations. In high rate oil and gas production, it may be beneficial to flow fluids vertically through the annulus configuration between well tubing and casing. The flow patterns in annuli are different from pipe flow. There are both casing and tubing liquid films in slug flow and annular flow in the annulus. Multiphase heat transfer depends on the hydrodynamic behavior of the flow. There are very limited research results that can be found in the open literature for multiphase heat transfer in wellbore annuli. A mechanistic model of multiphase heat transfer is developed for different flow patterns of upward gas–liquid flow in vertical annuli. The required local flow parameters are predicted by use of the hydraulic model of steady-state multiphase flow in wellbore annuli recently developed by Yin et al. The modified heat-transfer model for single gas or liquid flow is verified by comparison with Manabe’s experimental results. For different flow patterns, it is compared with modified unified Zhang et al. model based on representative diameters.


Introduction
As oil and gas development moves from land or shallow water to deep and ultradeep waters, multiphase flow occurs during production and transportation (Chen 2011). The flow normally occurs in horizontal, inclined, or vertical pipes and wells. Gas-liquid two-phase flow in an annulus can be found in a variety of practical situations. In high rate oil and gas production, it may be beneficial to flow fluids vertically through the annulus configuration between well tubing and casing. For surface facilities, some large, underutilized flow lines can be converted to dual-service by putting a second pipe through the large line (i.e., flowing produced water in the inner line and gas in the annulus). During gas production, liquids may accumulate at the bottom of the gas wells during their later life. In order to remove or ''unload'' the undesirable liquids, a siphon tube is often installed inside the tubing string, which would form a gas-liquid two-phase flow in the annulus. Flow-assurance problems, such as hydrate blocking (Jamaluddin et al. 1991;Li et al. 2013;Wang et al. 2014Wang et al. , 2016 and wax deposition (Zhang et al. 2013;Bryan 2016;Theyab and Diaz 2016), are strongly associated with both the hydraulic and thermal behavior. For example, they are related to the fluid velocity, liquid fraction, slug characteristics, pressure gradient and convective-heat-transfer coefficients of different phase and flow patterns in multiphase flow. Therefore, multiphase hydrodynamics and heat transfer in an annulus need to be modeled properly to guide the design and operation of flow systems.
Compared to theoretical studies of multiphase hydrodynamics (Liu et al. 2007; Wang and Sun 2009), and multiphase heat transfer in pipe flow (Zheng et al. 2016;Gao et al. 2016;Karimi and Boostani 2016;Rushd and Sanders 2017), there are very limited research results in the open literature for multiphase heat transfer in wellbore annuli. Davis et al. (1979) obtained a model for calculating the local Nusselt numbers of stratified gas/liquid flow in turbulent liquid/turbulent gas conditions. The model was tested with heat-transfer experiments for air/water flow in a 63.5-mm inside diameter (ID) tube. Guo et al. (2017) established a mathematical model of heat transfer in a gasdrilling system, considering the flowing gas, formation fluid influx, Joule-Thomson cooling and entrained cuttings in the annular space. However, the multiphase flow effect on the heat transfer was not considered. Shoham, et al. (1982) undertook experiments on heattransfer for slug flow in a horizontal pipe. He found a substantial difference in heat-transfer coefficient existed between the top and bottom of the slug. Developing heattransfer correlations of different flow patterns was the aim of most previous studies (Shah 1981). Twenty heat-transfer correlations were compared in Kim's study (Kim et al. 1997). He collected the experimental data from the open literature and recommended the correlations for different flow patterns. However, the errors with the experimental results by Matzain (1999) were large. Later, a comprehensive mechanistic model about heat transfer in gas-liquid pipe flow was obtained (Manabe 2001). It was compared with the experimental data, and the performance was better. However, there were some inconsistencies in annular and slug flow. It needed to be modified. Zhang et al. (2006) developed a unified model of multiphase heat transfer for different flow patterns of gas-liquid pipe flow at all inclinations -90°to ? 90°from the horizontal. The required local flow parameters were predicted by use of the unified hydrodynamic model for gas/liquid pipe flow developed by Zhang et al. (2003a, b). However, it is not fit for the gas-liquid flow in an annulus, because the flow patterns in annuli are different from pipe flow patterns, as seen in Fig. 1 (Caetano 1986). A new heat-transfer model for gas-liquid flow in vertical annuli needs to be established.
A hydraulic model was developed to predict flow patterns, liquid holdup and pressure gradients for steady-state gas-liquid flow in wellbore annuli (Yin et al. 2014). The major advantage of this model compared with previous mechanistic models is that it is developed based on the dynamics of slug flow, and the liquid-film zone is used as the control volume. The effects of the tubing liquid film, casing liquid film and the droplets in the gas core area on the mass and momentum transfers are considered. Multiphase heat transfer depends on the hydrodynamic behavior of the flow. The objective of this study is to develop a heat-transfer model for gas-liquid flow that is consistent with the hydrodynamic model in vertical wellbore annuli.

Single-phase flow
When fluids flow through an annulus, the surrounding temperature is cold, heat is lost from the fluids to the formation, resulting in a decline in temperature, as seen in  In the steady state, there is no heat loss from the annulus to the tube. For the fluids in the annulus, the conservation of mass is: where q l is the density of liquids, kg/cm 3 ; v l is the velocity of liquids in the annulus, m/s; dl is the length of the element. The conservation of momentum is: where p is the unit pressure, MPa; g is the gravity acceleration, m/s 2 ; h is the inclination angle,°; s is the friction; d is the annulus diameter, m; A is the cross area of the annulus, m 2 . The energy equation is: where e is the internal energy of the unit, J/kg; q 1 is the heat loss from the annulus to the formation, W.
Using the mass balance, we can reduce Eqs.
(2) and (3) further: Or where h is the enthalpy, J/kg; w l is the mass flow rate of liquids in the annulus, kg/s. Using the definition of the dimensionless temperature T D proposed by Hasan and Kabir (1991), we may write an expression for heat transfer from the wellbore/formation interface to the formation as where T e is the formation temperature,°C; T wb is the wellbore temperature,°C; k e is the conductivity of formation, W/(m°C); the dimensionless temperature, defined by Hasan and Kabir (1991), T D , can be easily estimated from the following models.
where t D is the dimensionless flowing time, where t is the flowing time, s; r w is the wellbore radius, m. The heat transfer from the annulus to the well wall, q 1 may be expressed in terms of an overall heat coefficient (Hansan and Kabir 1994), where r co is the outside casing radius, m; T a is the fluid temperature in the annulus,°C; U a is the overall heattransfer coefficient of the annulus, W/(m 2°C ).
where r ci is the inside casing radius, m; h a is the convective-heat-transfer coefficient of fluids in the annulus, W/ (m 2°C ); k cas is the conductivity of the casing wall, W/(m°C ); k cem is the conductivity of cement, W/(m°C). Combining Eqs. (7) and (8) gives, where c p is the fluid heat capacity at the constant pressure in the annulus, J/(kg°C).
where A 0 is the local parameter defined in Eq. (10). Combining Eqs. (6), (8) and (10) yields, The enthalpy gradient can be written in terms of the temperature and pressure gradients: where g l is the fluid Joule-Thomson coefficient in the annulus,°C/MPa. Combining Eqs. (11) and (12) gives Defining a dimensionless parameter U a , as We can write Eq. (13) as Equation (14) is the energy conservation equation of fluids in the annulus.
We can write Eq. (14) as Up to this point, only mathematical manipulations have been done to the enthalpy equation, and the analysis has been carried out rigorously without simplification. Now if Eq. (15) is used for an onshore production well, then assuming the surrounding formation temperature is a linear function of depth, it can be expressed as where T ei is the temperature of formation at wellbore intake,°C; g e is the formation thermal gradient,°C/100 m; L is the depth, m. If Eq. (15) is used for the riser of an offshore production well, then the surrounding sea temperature is not a linear function of depth. It will be calculated according to the actual environment (Wang and Sun 2009).
If, for a certain segment of the wellbore, U a , c p , g l , g e , h, v l dv l /dl and dp/dl can be approximately constants, combining Eqs. (15) and (16) and integrating, yields an explicit equation for the temperature:

Bubbly flow and dispersed-bubble flow
In bubbly flow and dispersed-bubble flow, the gas holdup is small and the gas superficial velocity is low, the gas phase is distributed as small discrete bubbles in a continuous liquid phase. So bubbly flow and dispersed-bubble flow can be treated as pseudo-single-phase flow. The fluid physical properties are adjusted based on the liquid holdup. Zhang et al. correction (2006) for bubbly flow will be modified based on ''hydraulic diameter''.
where d R is the hydraulic diameter, m; d ci is the inside diameter of casing, m; d tuo is the outside diameter of tubing, m. Then the convective-heat-transfer coefficient of Eq. (9) for bubbly or dispersed-bubbly flow is obtained from where k m is the conductivity of the mixture, W/(m°C), where H l is the liquid hold up; k l is the conductivity of liquid, W/(m°C); k g is the conductivity of gas, W/(m°C). N m Nu is the mixture Nusselt number (Zhang et al. 2006).
where f m is the friction factor of the mixture; l l is the viscosity of liquid, mPa s; l m is the viscosity of the mixture, mPa s; l g is the viscosity of gas, mPa s, Pr is the Reynolds number and Prandtl number of the mixture.
where q m is the density of the mixture, kg/m 3 ; v m is the velocity of the mixture, m/s; c pm is the specific heat of the mixture, J/(kg°C), where c pl is the specific heat of liquid, J/(kg°C); c pg is the specific heat of gas, J/(kg°C). So the associated parameters will be modified, þ q m g sin h À q m g m c pm dp dl 0 dp dl where the subscript b represents in bubble flow or dispersed-bubble flow; the subscript m represents the mixture properties of gas and liquid, A 0 b is the local parameter defined in Eq. (20).
Equation (15) can be written as Then Eq. (20) can be used in the heat transfer in bubble flow and dispersed-bubble flow based on the modified convective-heat-transfer coefficient h ab .

Annular flow
As shown in Fig. 3, there are tubing and casing films in annular flow in annuli, which is different from the annular flow in pipes. Assume there are no temperature and heattransfer exchanges in the vertical direction. The temperature and heat changes are only caused by the heat transfer in the radial direction. Heat distribution in annuli includes three parts: in casing film, tubing film and gas core. The heat-transfer models in annuli are obtained by a method similar to that above.
For the gas core, the energy conservation equation is as follows, where T gc is the temperature of the gas core,°C; T cf is the temperature of the casing film,°C; T tuf is the temperature of the tubing film,°C; c pgc is the specific heat of the gas core, J/(kg°C); A 0 gc , B 0 gc , U gc are the local parameters defined as follows: U gc ¼ q gc v gc dv gc dl þ q gc g sin h À q gc g gc c pgc dp dl 0 dp dl where w gc is the mass flow rate of the gas core, kg/s; c pgc is the specific heat of the gas core, J/(kg°C); v gc is the velocity of the gas core, m/s; q gc is the density of the gas core, kg/m 3 ; g gc is the Joule-Thomson coefficient of the gas core,°C/MPa; U gc is the overall heat-transfer coefficient of the gas core, W/(m 2°C ), and defined as where r tuo is the outside radius of the tubing, m; d tu is the thickness of the tubing film, m; d c is the thickness of the casing film, m; h gc is the convective-heat-transfer coefficient of the gas core, W/(m 2°C ). U tuf is the overall heat-transfer coefficient of the tubing film, W/(m 2°C ), and defined as where h tuf is the convective-heat-transfer coefficient of the tubing film, W/(m 2°C ); k tu is the conductivity of the tubing wall, W/(m°C).
For the tubing film, the energy conservation equation is as follows, where w tuf is the mass flow rate of the liquid tubing film, kg/s; q tuf is the density of the liquid tubing film, kg/m 3 ; c ptuf is the specific heat of the liquid tubing film, J/(kg°C); U tuf is the local constant and defined as follows: U tuf ¼ q tuf v tuf dv tuf dl þ q tuf g sin h À q tuf g tuf c ptuf dp dl 0 dp dl where v tuf is the velocity of the liquid tubing film, m/s; g tuf is the Joule-Thomson coefficient of the liquid tubing film,°C /MPa. For the casing film, the energy conservation equation is as follows, where w cf is the mass flow rate of liquid casing film, kg/s; q cf is the density of liquid casing film, kg/m 3 ; c pcf is the specific heat of liquid casing film, J/(kg°C); C 0 cf is the local parameter defined in Eq. (23), H gc (l ) H tuf (l ) Tubing liquid film Casing liquid film Fig. 3 Control volume and temperatures in annular flow (q tufc is the heat transfer from the tubing film to the gas core, W; q ccf is the heat transfer from the gas core to the casing film, W; q cff is the heat transfer from the casing film to the formation, W) Pet. Sci. (2018) 15:135-145 139 U cf is the overall heat-transfer coefficient of casing film, W/(m 2°C ), and defined as where h cf is the convective-heat-transfer coefficient of casing film, W/(m 2°C ).
U cf is the local parameter defined in Eq. (23), where v cf is the velocity of the liquid casing film, m/s; g cf is the Joule-Thomson coefficient of the liquid casing film,°C /MPa. Equations (21)-(23) are the heat-transfer models for gas-liquid annular flow in annuli. Then the fluid temperature in the wellbore is calculated by a weighted average method based on holdup.
where H tuf is the holdup of the liquid tubing film; H cf is the holdup of the liquid casing film. The convective-heat-transfer coefficients for the casing film, the tubing film and gas core are obtained by the following equations where k lf , k gc are the thermal conductivities of the liquid film and gas core, W/(m°C); d tufR , d cfR , d gcR are the ''hydraulic diameter'' of the tubing film, the casing film and gas core, m.
The Nusselt numbers for the liquid film and gas core are calculated by using the correlations for single-phase convective heat transfer. The Petukhov correlation (1970) is used for turbulent liquid-film flow: where i is tubing film or casing film, i = tu or c. f if is the friction factor at the wall in contact with the liquid film.
The Dittus and Boelter correlation (1985) is used for turbulent gas core flow, For fully developed laminar flows of the liquid film and gas core, the Nusselt number approaches a constant value. According to Zhang et al. (2006) The associated hydraulic parameters are calculated by the Yin et al. model (Yin et al. 2014), such as the thickness of the liquid film, d c and d tu , the hold up, liquid and gas physical properties .

Film region
The flow characters of the film region are similar to the annular flow. The difference is that the gas core with liquid droplets changes into the Taylor bubble. So, the overall heat-transfer coefficients of Eqs. (21) and (23) should change from U gc to U T . Then the heat transfer of the film region can be calculated.
where U T is the overall heat-transfer coefficient of the Taylor bubble, W/(m 2°C ); h T is the convective-heattransfer coefficient of the Taylor bubble, W/(m 2°C ).
where k T is the thermal conductivities of the Taylor bubble, W/(m°C); d TR is the ''representative diameter'' of the Taylor bubble.

Slug region
There are small discrete bubbles in a continuous liquid phase. The flow characters are similar to bubbly flow. So, the heat transfer can be calculated by the bubbly flow model.

Slug unit
The fluid temperature in the wellbore is calculated by the weighted average method based on holdup.
where T s is the temperature of the fluid in slug flow,°C; T ls is the temperature of the fluid in the slug unit,°C; T T is the temperature of the Taylor bubble,°C; l F is the length of the liquid film, m; l s is the length of the slug unit, m; l U is the length of the whole slug, m.
3 Solution procedure Figure 4 shows the overall solution flowchart for the present model. First, the parameters of fluid properties are provided. Then the flow pattern is determined based on the input variables and then all the flow conditions, such as flow pattern, liquid holdups, local fluid velocities of the liquid film and gas core, and slug characteristics, are predicted by use of the hydraulic model of steady-state multiphase flow in wellbore annuli recently developed by Yin et al. (2014). Based on flow patterns, further calculations are performed in different subroutines. If the flow pattern is single-phase flow, single-phase heat-transfer calculation will be performed; if the flow pattern is bubble or dispersed-bubble flow, the corresponding hydraulic model and heat-transfer model will be called for calculation. For annular flow, corresponding hydraulic and heat-transfer calculations will be made. For slug flow, corresponding hydraulic and heat-transfer calculations will be made, as seen in Fig. 4. Figure 5 is the flowchart for the present annular flow heat-transfer model.

Comparisons with the unified Zhang model
The heat transfer of single-phase gas or liquid flow can be calculated by the present model if the thermal conductivity of the tubing wall is infinite. The composition of simulated fluids is shown in Tables 1 and 2. The simulated annulus outer diameter is 76.2 mm, and the inner diameter is 42.2 mm. The basic parameters are shown in Table 3. Figures 6 and 7 show comparisons between the simulations of new models and experimental measurements (Manabe 2001) of the convective-heat-transfer coefficients for single-phase gas and liquid flows, respectively. The good agreement between the simulations and experimental measurements for single-phase gas and liquid flows indicates that the new models are reliable and the selected correlations are appropriate.
There are few experimental research results in the open literature for multiphase heat transfer in annuli. The unified Zhang et al. model (2006) is verified by comparison with Manabe's experimental results for different flow patterns in a crude-oil/natural gas system, and good agreement has been observed in the comparison. So, the unified Zhang et al. model is modified to calculate the heat transfer of gas-liquid flow in annuli based on the ''hydraulic diameter'' of the annuli, Eq. (17), and the results are compared with the present mechanistic heat-transfer model for gasliquid flow in annuli. Figures 8, 9, and 10 are comparisons of convectiveheat-transfer coefficient for bubble flow, annular flow and slug flow predicted by the present model and the modified unified Zhang et al. model (2006). For the bubble flow, the data points are located inside the 10% error band. The agreement is good. It shows that the influence of annulus geometry is small for the low gas volume fraction. For the annular and slug flow in annuli, most of the data points are located inside the 30% error band and all are overestimated. It may because there is a tubing liquid film and a

Conclusions and discussion
A heat-transfer model for gas-liquid flow in vertical annuli is developed in conjunction with the mechanistic hydrodynamic model of Yin et al. (2014), which can predict flow pattern transitions, liquid holdup, gas void fraction, pressure gradient, and slug characteristics in gas-liquid twophase flow in vertical annuli. The heat-transfer modeling is based on energy-balance equations and analyses of the temperature differences and variations in the tubing liquid film, casing liquid film, gas core, Taylor bubble and slug body.
The heat-transfer model for single gas or liquid flow is verified by comparison with Manabe's experimental results (2001). Good agreement has been observed in the comparison. For different flow patterns, it is compared with unified Zhang et al. model modified based on ''hydraulic diameter''. For bubble or dispersed-bubble flow, the error is lower than 10%. With the gas void fraction and gas flow velocity increasing, the error will be larger but lower than 30%. In other words, the difference between the new model and modified unified Zhang et al. model will be small if the gas void fraction and velocity is small. The difference will be large when it changes. The modified method based on ''hydraulic diameter'' is no longer applicable for slug and annular flow in vertical annuli when the gas void fraction increases. It may be 1.3 times larger than the new model. Experimental investigations of heat transfer in vertical are required to improve the model performance.