Heat Transfer Mechanisms in Petroleum and Geothermal Wellbores

Abstract

The flow of fluids between wells and reservoirs involves a substantial amount of thermal energy exchange with the formation. Understanding the mechanisms involved in the heat transfer of these processes is crucial to the design of the wells for mechanical integrity. While long term production scenarios may achieve a notional steady state, short term injection scenarios involve an accurate consideration of the thermal transients. With global initiatives towards a transition to clean energy, the design of geothermal wells is becoming an area of great importance these days. Accordingly, correct simulation of the heat transfer in the circulating scenario involved in closed loop wells enables accurate assessments of thermal power generated. This chapter aims to educate the user in how to tackle these problems and explains the physics and mathematics involved in detail.

1. Introduction

1.1 Background

The discovery of reservoirs with hotter in-situ temperatures (above 200° F) over the past several decades has introduced engineering challenges that depend critically on an accurate assessment of wellbore temperatures. In particular, subsea wells are being drilled to deeper horizons these days and are exposed to hotter temperatures than in the past.

These wells have multiple tubulars and fluid-filled annuli as depicted in Figure 1. In addition, many of these wells are prolific producers (of hydrocarbons or geothermally heated water), resulting in high arrival temperatures at the surface. In some instances, the fluid arrival temperatures at the wellhead could, in fact, be hotter than the already high bottomhole temperature, because of the negative Joule-Thomson effect. A problem of equal, if not greater importance, is the effect of the lateral (or radial) heat transfer from the flowing stream to the wellbore layers, resulting in temperature buildup in fluid filled annuli and thermal stresses in the unsupported sections of the bounding tubulars.¹ One of the most serious implications of radial heat transfer is Annular Pressure Buildup (APB). The prediction and mitigation of APB constitutes a vast body of investigation in its own right. Thermal stresses in tubulars influence the structural design of the wellhead, and the control of Wellhead Movement (WHM). The displacement constraints on the tubulars at the wellhead and the tops of cement can cause buckling and the generation of bending stresses during well operation. In a worst case discharge (WCD) scenario, elevated temperatures may potentially dislodge tubulars from the wellhead, and require additional lock down rings to prevent the tubulars from catapulting. All of these phenomena require accurate and reliable estimates of wellbore temperatures. In instances that involve operations with short durations, accurate prediction of the thermal transient response is critical (for example. Drillstem tests, Well Testing to evaluate reservoir performance, Designing APB mitigation mechanisms, wellhead pressure control in platform wells). Injection and circulation scenarios also create temperature changes that generate unsustainable tensile forces in improperly designed wellbore tubulars and tubular connections.

Figure 1.

Schematic of a complex wellbore with multiple annuli and bounding tubulars.

1.2 Heat transfer mechanisms in wellbores

The fundamental mechanisms of heat transfer in a wellbore are indicated in Figure 2. In most of the onshore and offshore locations, the geothermal temperature increases with depth below the surface, at an average of rate of 21–32°C/km. This temperature gradient is the primary driver for all heat exchange phenomena in a wellbore. This is true of wellbores used to extract oil and gas, and of wellbores used to generate geothermal energy.

Figure 2.

Illustration of the various heat transfer phenomena in a wellbore.

For the purposes of thermal and structural analysis, a well can be enclosed in an imaginary volume that encloses the production tubing (i.e. the innermost cylinder and primary flow conduit), and the series of casings and cement sheaths in the intervening annular spaces. The well boundary is located at the interface between the outermost cement sheath and the earth (known hereafter as the formation).

In a wellbore, energy is exchanged between the flow stream(s), the wellbore (i.e. the casing strings and annular contents) and the formation. The thermal analysis of the producing wellbore proceeds in three interlinked steps. The first step is the solution of the balance (mass, momentum, and energy) equations in the tubing. The second step is the assessment of radial heat loss from the tubing to the wellbore. For the purposes of thermal analysis, the wellbore is defined as the region between the outer surface of the tubing and the outer surface of the outermost cement sheath. The third step is the determination of the heat transfer in the formation, i.e. from the wellbore – formation boundary to the earth.

There is forced convection heat transfer between the flowing fluid stream and the conduit boundary. Usually, the uncemented annular sections between tubulars contain incompressible fluids that experience natural convection. Conduction in the radial direction occurs through the walls of the casing, and the cemented sections of the intervening annuli. This is a case of diffusion across in a composite medium. At the well boundary, heat lost by the contents of the wellbore diffuses by conduction into a semi-infinite domain. Sometimes the semi-infinite domain is approximated by a finite domain with a very large farfield radius. In some wells, there is a need to minimise heat loss from the wellbore. In such applications, Vacuum Insulated tubing (VIT) is used.² The heat transfer in this case between the inner and outer pipes is practically by thermal radiation.

1.3 Types of well thermal operations

In terms of thermal interactions, a wellbore is essentially a heat exchanger. Conventional heat exchangers typically involve heat transfer between two counterflowing or parallel streams. In a wellbore however, a single stream flowing up (production) or down (injection) the wellbore, exchanges heat with the formation layers surrounding the wellbore, as indicated in the left two panels of Figure 3. In this figure, the black dotted line represents the geothermal temperature, which prevails in the wellbore until an operation (or operations) induce a thermal disturbance. During production, the hot fluid exits from the reservoir at the bottom of the well and flows upward. During it upward transit, there is loss of fluid enthalpy because of lateral/radial heat transfer. This is responsible for the heating of the tubulars, the annular contents in the well (solid red curve, panel (a)). During injection, cold fluid gets heated during its downward transit (blue curve, panel (b)). The right two panels indicate circulation scenarios which are analogous to classic counterflow heat exchangers. In both cases, qualitative descriptions of the associated temperature profiles are indicated (solid red and blue curves).

Figure 3.

Producing (a), injection (b), forward circulating (c), and reverse circulating (d) scenarios.

In all the three scenarios in Figure 3, the key objective of a thermal analysis is the prediction of the flowing temperature profiles, given appropriate boundary conditions. In the case of transient heat transfer, initial conditions must also be specified. In production and injection scenarios, the (boundary condition) temperatures are either known or stipulated at the reservoir and wellhead locations. In forward circulation, the temperature is specified at the wellhead location of the inner conduit, whereas in reverse circulation the temperature is specified at the wellhead location of the outer (annular) conduit. In either case, the temperature is specified at the inlet to the wellbore of the downward flowing stream. At the bottom of the wellbore, it is typical (but not always) to stipulate the equality of the temperatures of the two flowing streams, as shown in Figure 3 (panels (c) and (d)).

The analysis involves the solution of the transport equations in conjunction with heat transfer in the formation. This requires careful consideration of all relevant fluid and thermal transport phenomena. The subsequent sections will present a systematic analytical approach to the solution of the aforementioned problems.

1.4 Review of relevant literature

The earliest studies of heat transfer in wellbores by Lesem et al. [1] and Moss and White [2] date back to the late 1950s. For a detailed review of the literature on the topic, the reader is referred to Chandrasekhar [3] wherein a comprehensive transient thermal model of a complex wellbore is described in detail. There are several key studies that constitute essential reading and are listed below:

The 1962 study by Ramey [4] was the first systematic study of both flowing and wellbore temperatures. His approach assumed pseudo steady state conditions in the flowing conduit and wellbore, with the transients relegated solely to the formation. This approach is in fact the basis for a very large number of model implementations (the WELLFLO code for example) to this day. The wellbore itself was modelled as a line source in a semi-infinite formation for which a simple expression was used to characterise the transient heat flux. While the approach breaks down for shorter producing intervals, it is valid for time periods corresponding to Fourier numbers in excess of unity.

Willhite [5] extended the approach of Ramey [4] to account for amongst other phenomena, natural convection, and thermal radiation in fluid-filled annuli. An iterative approach is required to calculate the overall heat transfer coefficient linking the temperature of flowing stream to the far field undisturbed geothermal temperature.

It is very likely that Raymond [6] was the first study of the transient circulation problem using a combine Laplace Transform/Finite Difference approach. The key observation of Raymond’s analysis is that the transients are limited to the first few hours of circulation and that the steady state solution was valid for longer periods. The first detailed study of multiple well operating scenarios is that of Wooley [7] in which production, injection, and circulation were considered in the context of a transient analysis using a finite difference approach to couple the well and formation responses.

More recent studies have looked at analytical solutions where possible for coupled wellbore/formation problems. Wu and Pruess [8] considered transient heat transfer between a flowing fluid stream and the formation, but used an overall lumped heat transfer coefficient to model the heat transfer across the wellbore itself. They formulated a more refined formation temperature model using Laplace transforms to model a cylindrical source. The 2018 study of Chandrasekhar et al. [9] is recommended for the reader interested in the application of a circulating model to a complex realistic wellbore considering both hydraulics and thermal phenomena, in addition to several other aspects of actual real-life wellbores.

There are several textbooks in the literature that present a detailed analysis of the fundamentals of wellbore heat transfer. Hasan and Kabir [10] cover several aspects of both heat transfer and fluid flow in wellbores, starting with the governing equations, and several models for multiphase flows in wellbores. In a 2009 SPE monograph, Mitchell and Sathuvalli [11] discuss various phenomena and analytical techniques relevant to temperature prediction in prolific oil and gas producers.

There are a few experimental studies that have investigated aspects of wellbore heat transfer. Jones [12] describes a real time measurement that was quite novel at the time approach to establish circulating temperatures in wellbores during drilling and cementing operations. The performance of Vacuum-Insulated Tubing was studied by Aeschliman et al. [13] in the context of a steam injection well. Their results compared six different commercially available means of achieving thermal insulation by the suppression of convection in the tubing annulus.

2. Governing transport equations

2.1 Mass conservation

Consider a control volume (CV) of length $\mathrm{\Delta }z$ and a fixed radius $R$ as shown in Figure 4. Mass, momentum, and energy enter and leave the CV at locations $z$ and $z+\mathrm{\Delta }z$. In wellbores, usually there is no mass accumulation at a given location in the control volume, so that the constant mass flow rate is given by

Figure 4.

Control volume for mass, momentum, and energy balances.

$$\dot{m}=\rho \left(z\right)V\left(z\right)=\rho \left(z+\mathrm{\Delta }z\right)V\left(z+\mathrm{\Delta }z\right)A$$E1

where $A=\pi R^2$ is the conduit flow area. A mass balance over the control volume in the limit that the size $\mathrm{\Delta }z$ shrinks to zero yields

$$\frac{\partial }{\partial z}\left(\mathit{\rho V}\right)=-\frac{\partial \rho }{\partial t}=0$$E2

so that the instantaneous temporal derivative of the density is zero. From Eqs. (1) and (2), we have

$$\frac{\partial V}{\partial t}=\frac{\partial }{\partial t}\left(\frac{\dot{m}}{\mathit{\rho A}}\right)=\frac{\dot{m}}{A}\frac{\partial }{\partial t}\left(\frac{1}{\rho }\right)=-\frac{V}{\rho }\frac{\partial \rho }{\partial t}=0$$E3

whereupon the temporal derivative of the velocity also vanishes, so that along with no mass accumulation, there is no accumulation of momentum in the control volume either.

2.2 The energy equation

The specific energy $e\left(z,t\right)$ identified in the figure is the sum of the kinetic, potential, and internal energies such that

$$e=u+\frac{1}{2}{V}^2-\mathit{\delta gy}$$E4

where $y$ is the vertical depth relative to some fixed datum, and the negative sign associated with it implies a loss of potential energy as the vertical depth increases. The term $\delta =\pm 1$ defines the orientation of the streamwise coordinate relative to the gravity vector, such that $\delta =1$ and $\delta =-1$ describe injection and production scenarios, respectively. The issue of how to deal with these terms in a circulation scenario will be described later.

Energy enters and leaves the CV in Figure 4 as indicated, with some accumulation at the rate $\mathrm{\Delta }e$ over a period $\mathrm{\Delta }t$. Energy is supplied through the conduit boundary at the rate $\dot{q}\left(t\right)$.

A simple energy balance yields the following expression

$$\dot{m}\left[e\left(z+\mathrm{\Delta }z\right)-e\left(z\right)\right]+\mathit{\rho A}\mathrm{\Delta }z\frac{\mathrm{\Delta }e}{\mathrm{\Delta }t}=2\mathit{\pi R}\underset{s=z}{\overset{s=z+\mathrm{\Delta }z}{\int }}\dot{q}\left(s\right)\mathrm{d}s-\dot{m}\mathrm{\Delta }\left(\frac{P}{\rho }\right)=\left(2\mathit{\pi R}\mathrm{\Delta }z\right)\dot{q}\left(z+\lambda \mathrm{\Delta }z,t\right)$$E5

where the Mean Value Theorem has been used to replace the integral such that $0<\lambda <1$. Dividing by $\mathrm{\Delta }z$ and taking the limit as both $\mathrm{\Delta }z\to 0$ and $\mathrm{\Delta }t\to 0$ results in

$$\begin{array}{c}lim\\ \begin{array}{c}\mathrm{\Delta }z\to 0\\ \mathrm{\Delta }t\to 0\end{array}\end{array}\left\{\dot{m}\left[\frac{e\left(z+\mathrm{\Delta }z\right)-e\left(z\right)}{\mathrm{\Delta }z}\right]+\mathit{\rho A}\frac{\mathrm{\Delta }e}{\mathrm{\Delta }t}\right\}=\begin{array}{c}lim\\ \begin{array}{c}\mathrm{\Delta }z\to 0\\ \mathrm{\Delta }t\to 0\end{array}\end{array}\left\{2\mathit{\pi R}\dot{q}\left(z+\lambda \mathrm{\Delta }z\right)\right\}$$E6

Evaluating the limit and noting that $\dot{m}=\mathit{\rho AV}$ simplifies Eq. (6) to

$$\rho \frac{\partial e}{\partial t}+\mathit{\rho V}\frac{\partial e}{\partial z}=\frac{2}{R}\dot{q}\left(z\right)$$E7

From one of the fundamental thermodynamic relationships relating the enthalpy to internal energy, we have

$$h=u+\frac{P}{\rho }$$E8

Substitution of the above along with Eq. (4) into Eq. (7) yields

$$\rho \left[\frac{\partial u}{\partial t}+V{\cancel{\frac{\partial V}{\partial t}}}^0+\cancel{\frac{\partial }{\partial t}}{\left(\mathit{\delta gy}\right)}^0\right]+\mathit{\rho V}\left(\frac{\partial h}{\partial z}+V\frac{\partial V}{\partial z}-\mathit{\delta g}\right)=\frac{2}{R}\dot{q}\left(z,t\right)$$E9

The derivative of the vertical depth $y$ with the streamwise coordinate $z$ (known as the Measured Depth in wellbore parlance) is the cosine of the local wellbore inclination $\theta$. From Eq. (8), the temporal derivative of the internal energy can be expressed as

$$\frac{\partial u}{\partial t}=\frac{\partial h}{\partial t}-\frac{\partial }{\partial t}\frac{P}{\rho }=\frac{\partial h}{\partial t}-\left[\frac{P}{\rho }{\cancel{\frac{\partial P}{\partial t}}}^0-\frac{P}{{\rho }^2}{\cancel{\frac{\partial \rho }{\partial t}}}^0\frac{P}{\rho }\right]=\frac{\partial h}{\partial t}$$E10

where the temporal derivative of the density vanishes in accordance with the constant mass flow criterion which also stipulates from Eq. (3) that the time derivative of the velocity is zero. In addition, it follows from the momentum conservation equation (which will be presented in the next section), that in the context of a constant mass flow rate the time derivative of the pressure also vanishes. Accordingly, Eq. (9) reduces to

$$\rho \frac{\partial h}{\partial t}+\mathit{\rho V}\left(\frac{\partial h}{\partial z}+V\frac{\partial V}{\partial z}-\mathit{\delta g}cos\theta \right)=\frac{2}{R}\dot{q}\left(z,t\right)$$E11

Note that Eq. (11) contains spatial derivatives of both specific enthalpy and velocity. Closure therefore requires the consideration of the momentum equation which is presented next. It is reiterated here that Eq. (11) as derived is only valid under the assumption of a constant mass flow rate throughout the wellbore.

2.3 The momentum equation

A force balance over the same control volume as in Figure 4 yields the rate of change of momentum. The forces acting on the fluid in the control volume are the static and dynamic pressure forces and the pressure and shear stress as indicated

$$\begin{array}{l}A\left[\rho \left(z+\mathrm{\Delta }z\right)V^2\left(z+\mathrm{\Delta }z\right)-\rho \left(z\right)V^2\left(z\right)\right]+A\left[P\left(z+\mathrm{\Delta }z\right)-P\left(z\right)\right]+\mathit{\rho A}\mathrm{\Delta }z{\cancel{\frac{\mathrm{\Delta }V}{\mathrm{\Delta }t}}}^0\\ =\mathit{\delta \rho }\overline{g}\left(z\right)A\mathrm{\Delta }z+2\mathit{\pi R}\mathrm{\Delta }z\underset{s=z}{\overset{s=z+\mathrm{\Delta }z}{\int }}\tau \left(s\right)\mathrm{d}s=\mathit{\delta \rho }\overline{g}\left(z\right)A\mathrm{\Delta }z+\left(2\mathit{\pi R}\mathrm{\Delta }z\right)\tau \left(z+\lambda \mathrm{\Delta }z\right)\end{array}$$E12

whereupon following the same logic as was used to derive the energy equation and noting from Eq. (3) that the time derivative of the velocity vanishes, we have

$$\mathit{\rho V}\frac{\partial V}{\partial z}+\frac{\partial P}{\partial z}=\mathit{\delta \rho g}cos\theta +\frac{2}{R}\tau \left(z\right)$$E13

which is functionally equivalent to Newton’s Second Law of Motion relating the rate of change of Momentum to the sum of the forces acting on a body of fluid.

Note that the stipulation of no mass accumulation also implies no momentum accumulation, so that the only accumulation in the wellbore is that of energy. Note also, that from Eq. (13), the time derivative of pressure is zero which enables the energy equation to be cast with enthalpy as the sole flux variable on both side of the equation.

2.4 Coupled transport equation system

The kinetic energy term in Eq. (11) is represented by the spatial derivative of the velocity. This term can be expressed in terms of pressure and enthalpy derivatives by invoking the chain rule as follows

$$\begin{array}{ll}\frac{\partial V}{\partial z}& =\frac{\partial }{\partial z}\left(\frac{\dot{m}}{\mathit{\rho A}}\right)=\frac{\dot{m}}{A}\frac{\partial }{\partial z}\left(\frac{1}{\rho }\right)=-\frac{V}{\rho }\frac{\partial \rho }{\partial z}\\ & =-V\left[\left(\frac{1}{\rho }{\left.\frac{\partial \rho }{\partial P}\right|}_h\right)\frac{\partial P}{\partial z}+\left(\frac{1}{\rho }{\left.\frac{\partial \rho }{\partial h}\right|}_P\right)\frac{\partial h}{\partial z}\right]=V\left[{\alpha }_h\frac{\partial h}{\partial z}-\beta \frac{\partial P}{\partial z}\right]\end{array}$$E14

where $\beta =\frac{1}{\rho }{\left.\frac{\partial \rho }{\partial P}\right|}_h$ is the adiabatic compressibility and ${\alpha }_h=-\frac{1}{\rho }{\left.\frac{\partial \rho }{\partial h}\right|}_P$ can be regarded as a two-phase isobaric volume expansivity.³ Substitution of Eq. (14) into Eqs. (11) and (13) yields the system of equations that can be expressed in the compact form

$$\begin{array}{l}\rho \frac{\partial }{\partial t}\left[\begin{array}{c}h\\ 0\end{array}\right]+\left[\begin{array}{cc}\mathit{\rho V}\left(1+V^2{\alpha }_h\right)& -\rho V^3\beta \\ \rho V^2{\alpha }_h& 1-\rho V^2\beta \end{array}\right]\frac{\partial }{\partial z}\left[\begin{array}{c}h\\ P\end{array}\right]\\ =\left[\begin{array}{c}2R^{-1}\dot{q}\left(z,t\right)+\mathit{\text{δρVg}}cos\theta \\ \mathit{\delta \rho g}cos\theta +2R^{-1}\tau \left(z\right)\end{array}\right]\end{array}$$E15

2.5 Extraction of wellbore temperatures

Subject to an initial condition for the enthalpy field in the wellbore and appropriate pressure and enthalpy and boundary conditions, Eq. (15) can be solved in conjunction with the constitutive models for the heat flux ($\dot{q}\left(z,t\right)$) and frictional resistance $\tau \left(z\right)$ terms. Once the enthalpy and pressure distributions are known, the temperature distribution is determined from the appropriate thermophysical property database⁴ or a correlation that has the functional form

$$T=T\left(h,P\right)$$E16

Since the heat flux term itself depends on temperature, the solution involves an iterative sequence at each depth. Furthermore, in a transient multiphase analysis, the coupling of the transport equations with the diffusion in the formation adjacent to the wellbore can present occasional challenges with respect to the latter, as is the case with modelling transient phenomena in steam injector wells.

2.6 Single phase flow (sensible heat)

In Two-Phase flow, the temperature remains constant under an isobaric change in enthalpy. In single phase flow however, an enthalpy change is related to changes in pressure and temperature according to

$$\mathrm{d}h=c_{\mathrm{p}}\mathrm{d}T-c_{\mathrm{p}}{c}_{\mathrm{JT}}\mathrm{d}P$$E17

where $c_{\mathrm{p}}$ and $c_{\mathrm{JT}}$ are the specific heat at constant pressure and the fluid Joule-Thomson Coefficient, respectively. The latter is related to the fluid volume expansivity according to

$$c_{\mathrm{JT}}=\frac{1}{\rho c_{\mathrm{p}}}\left(\mathit{\alpha T}-1\right)$$E18

For liquids with very low expansivity, the Joule-Thomson coefficient is invariably negative. The term “sensible heat” is used to refer to Eq. (17) since a change in enthalpy can be perceived as a change in temperature, which is not possible when the state point is inside the vapour dome of the fluid.

The sensible heat formulation can also be extended to multiphase flow in wellbores in conjunction with the so-called Black-Oil Model, where weighted properties are used for the fluid thermophysical properties in each of the phases. The key advantage of Eq. (17) is that the primary flux variables are now pressure and temperature. Accordingly, the spatial velocity derivative is now expressed as

$$\frac{\partial V}{\partial z}=-V\left[\left(\frac{1}{\rho }{\left.\frac{\partial \rho }{\partial P}\right|}_T\right)\frac{\partial P}{\partial z}+\left(\frac{1}{\rho }{\left.\frac{\partial \rho }{\partial T}\right|}_P\right)\frac{\partial T}{\partial z}\right]=V\left[\alpha \frac{\partial T}{\partial z}-\beta \frac{\partial P}{\partial z}\right]$$E19

where $\alpha =-\frac{1}{\rho }{\left.\frac{\partial \rho }{\partial T}\right|}_P$ is the (single-phase) isobaric volume expansivity, and $\beta =\frac{1}{\rho }{\left.\frac{\partial \rho }{\partial P}\right|}_T$ is the isothermal (not adiabatic) compressibility.

Substitution of Eq. (19) in Eq. (15) results in the system (see Chandrasekhar [3]):

$$\rho c_{\mathrm{p}}\frac{\partial }{\partial t}\left[\begin{array}{c}T\\ 0\end{array}\right]+\left[\begin{array}{cc}\mathit{\rho V}\left(c_{\mathrm{p}}+V^2\alpha \right)& -\mathit{\rho V}\left(c_{\mathrm{p}}{c}_{\mathrm{JT}}+V^2\beta \right)\\ \rho V^2\alpha & 1-\rho V^2\beta \end{array}\right]\frac{\partial }{\partial z}\left[\begin{array}{c}T\\ P\end{array}\right]=\left[\begin{array}{c}Q+\mathit{VH}\\ H+F\end{array}\right]$$E20

where

$$\begin{array}{l}{Q}_{\mathrm{tbg}}=2R^{-1}\dot{q}\left(z,t\right)=2\pi \frac{\overline{U}}{{\dot{m}}_{\mathrm{tbg}}}\left(T_{\mathrm{ann}}-T_{\mathrm{tbg}}\right)\\ H_{\mathrm{tbg}}=-gcos{\mathit{\theta \rho }}_{\mathrm{tbg}}{V}_{\mathrm{tbg}}\\ F_{\mathrm{tbg}}=2R^{-1}\tau \left(z\right)=\frac{1}{2D}\mathit{f\rho }{V}^2\end{array}$$E21

are the Thermal, Hydrostatic, and Frictional forcing functions respectively. If we ignore the transient term for the time being, then the 2×2 system in Eq. (20) can be inverted to yield the following expressions for the streamwise gradients of temperature and pressure as

$$\frac{\mathrm{d}T}{\mathrm{d}z}=\frac{\left(1-\omega \right)Q+\left(\mathit{\alpha TV}\right)H+\left(\mathit{\alpha T}+\omega -1\right)\mathit{VF}}{\mathit{\rho V}\left[\left(1-\omega \right)\left(c_{\mathrm{p}}+\mathit{\eta \alpha }\right)+\eta \left(\mathit{\alpha T}+\omega -1\right)\alpha \right]}$$E22

and

$$\frac{\mathrm{d}P}{\mathrm{d}z}=\frac{c_{\mathrm{p}}H+\left(c_{\mathrm{p}}+\mathit{\eta \alpha }\right)F-\left({g_{\mathrm{c}}}^{-1}\mathit{V\alpha }\right)Q}{\left[\left(1-\omega \right)\left(c_{\mathrm{p}}+\mathit{\eta \alpha }\right)+\eta \left(\mathit{\alpha T}+\omega -1\right)\alpha \right]}$$E23

where

$$\eta =V^2$$E24

and

$$\omega =\rho V^2\beta =\mathit{\eta \rho \beta }$$E25

For an incompressible liquid, $\alpha =\beta =\omega =0$, and accordingly the expression for the temperature gradient reduces to

$$\frac{\mathrm{d}T}{\mathrm{d}z}=\frac{Q-\mathit{VF}}{\rho {\mathit{Vc}}_{\mathrm{p}}}$$E26

and the pressure gradient reduces as it should, to

$$\frac{\mathrm{dP}}{\mathrm{d}z}=H+F$$E27

Before delving into the important physical aspects of Eq. (26), it will be useful to establish the constitutive models for the heat flux and fluid shear stress, which will be presented next.

3. Constitutive models

3.1 Heat flux

The formation adjacent to a wellbore is notionally a semi-infinite cylinder. Accordingly, a true steady state is never reached. However at large times from when a well is put into operation, a notional or pseudo-steady state condition is reached as shown by Ramey [4]. Under these conditions, the heat flux between the flowing fluid stream in the wellbore and the formation can be represented in terms of the local temperature difference between the fluid and the undisturbed formation temperature prevailing at a distance far from the wellbore, at any given depth. Mathematically this can be represented by the very simple form

$$\dot{q}\left(z\right)=-U\left[T\left(z\right)-T_{\mathrm{geo}}\left(z\right)\right]$$E28

where $U$ which will be described in more detail to follow, is an Overall Heat Transfer coefficient that is independent of time, and is associated with the conduit radius $R$, as reflected by the $2\mathit{\pi R}$ term in Eq. (5). Note that the factor of $2\pi$ itself has already been incorporated in Eq. (5) and therefore does not appear in Eq. (28).

3.2 Overall heat transfer coefficient

Consider the section of the wellbore shown in Figure 5 with multiple intervening layers between the fluid and the formation. At some notional steady state typically attained at long elapsed times after a well is put into operation, the fluxes across all of these layers are equal. In addition, this flux is also equal to the flux at the wellbore-formation interface at some frozen time instant $t$. This assumption corresponds to what is termed a Pseudo-Steady-State approach. With respect to the nomenclature of the figure, the heat flux per unit length (not unit area) is

Figure 5.

Wellbore layers between the transport fluid and the formation.

$$\begin{array}{ll}\frac{{\dot{q}}_L\left(z\right)}{2\pi }& =-\mathit{hR}\left(T_0-{\overline{T}}_0\right)=-\frac{k_1}{ln\raisebox{1ex}{${\overline{r}}_1$}\!\left/ \!\raisebox{-1ex}{${\overline{r}}_0$}\right.}\left({\overline{T}}_0-{\overline{T}}_1\right)=-\frac{k_2}{ln\raisebox{1ex}{${\overline{r}}_2$}\!\left/ \!\raisebox{-1ex}{${\overline{r}}_1$}\right.}\left({\overline{T}}_1-{\overline{T}}_2\right)=\cdots \\ & =-\frac{k_j}{ln\raisebox{1ex}{${\overline{r}}_j$}\!\left/ \!\raisebox{-1ex}{${\overline{r}}_{j-1}$}\right.}\left({\overline{T}}_{j-1}-{\overline{T}}_j\right)=\cdots =-\frac{k_N}{ln\raisebox{1ex}{${\overline{r}}_N$}\!\left/ \!\raisebox{-1ex}{${\overline{r}}_{N-1}$}\right.}\left({\overline{T}}_{N-1}-{\overline{T}}_N\right)\\ & =-R_{\mathrm{wb}}{k}_{\mathrm{Geo}}{\left.\frac{\partial T}{\partial r}\left(t\right)\right|}_{r=R_{\mathrm{wb}}}=-k_{\mathrm{Geo}}\left({\overline{T}}_N-T_{\mathrm{Geo}}\left(z\right)\right){\left.\frac{\partial \theta }{\partial \eta }\left(\tau \right)\right|}_{\eta =1}\end{array}$$E29

where the dimensionless time $\tau =\frac{{\alpha }_{\mathrm{Geo}}}{R_{\mathrm{wb}}^2}t$ is a Fourier Number. In Eq. (29) the barred entities refer to interface locations (layer boundaries) and the subscript 0 in the first term on the RHS of Eq. (29) refers to the fluid. Note that this term describes forced convection between the fluid and the conduit. The dimensionless flux in the last term of Eq. (29) is independent of the wellbore outer radius, interface temperature, and formation properties and is obtained from the solution of the diffusion problem in a cylindrical semi-infinite domain. Ramey [4] presented an expression for the dimensionless flux in terms of the Fourier Number⁵ based on an approximation of the line source solution as

$${\left.\frac{\partial \theta }{\partial \eta }\left(\tau \right)\right|}_{\eta =1}=F\left(\tau \right)=-{\left(ln\frac{1}{2\sqrt{\tau }}+0.29\right)}^{-1}$$E30

The constant heat flux per unit length across the wellbore represented by Eq. (29)and out into the formation at some snapshot in time can also be represented in terms of a fluid to formation temperature difference with the use of an Overall Heat Transfer Coefficient as

$$\frac{{\dot{q}}_L\left(z\right)}{2\pi }=\mathit{UR}\left(T-T_{\mathrm{Geo}}\right)$$E31

Eliminating the flux $\dot{q}\left(z\right)$ and the temperatures between Eqs. (29) and (31) results in the following expression for the overall heat transfer coefficient

$$U=\frac{1}{R}{\left[\frac{1}{\mathit{hR}}+\sum _{k=1}^{k=N}\frac{ln\raisebox{1ex}{$r_j$}\!\left/ \!\raisebox{-1ex}{$r_{j-1}$}\right.}{k_j}+\frac{1}{k_{\mathrm{Geo}}F\left(\tau \right)}\right]}^{-1}$$E32

Eq. (32) considers the following phenomena.

• forced convection in the conduit,

• thermal resistances offered by all of the intervening wellbore layers,

• thermal resistance at the wellbore-formation interface.

It is often convenient in wellbore heat transfer analysis to work with an overall conductance rather than a coefficient. In the context of Eq. (32), this is defined as $\overline{U}=\mathit{UR}$, so that Eq. (28) can be rewritten as

$$\dot{q}\left(z\right)R=-\overline{U}\left[T_0\left(z\right)-T_{\mathrm{geo}}\left(z\right)\right]$$E33

where the overall conductance is the reciprocal of the term in brackets in Eq. (32).

3.3 Natural convection in fluid-filled annuli

The thermal conductivity in each wellbore layer depends on the medium of the layer. In the case of tubulars and cemented sections, the thermal conductivity may be regarded as constant. Typical values for steel and cement are 45 W/m-K and 1 W/m-K, respectively. When the layer consists of a fluid however, it is subject to natural convection that must be considered in the analysis. Therefore, the conductivity of a fluid layer as used in Eq. (32) should be replaced by an equivalent thermal conductivity $k_{\mathrm{eq}}$ that accounts for natural convection. In terms of a heat transfer coefficient the flux due to natural convection between the inner and outer walls of the layer is given by

$${\dot{q}}_i\left(z\right)=h_i{\overline{r}}_{i-1}\left({\overline{T}}_{i-1}-{\overline{T}}_i\right)=\frac{k_{\mathrm{eq}}}{ln\raisebox{1ex}{${\overline{r}}_i$}\!\left/ \!\raisebox{-1ex}{${\overline{r}}_{i-1}$}\right.}\left({\overline{T}}_{i-1}-{\overline{T}}_i\right)$$E34

The natural convection correlation used in this context is an extension of the one proposed by Dropkin and Sommerscales [14] as suggested by Willhite [5] such that the equivalent conductivity $k_{\mathrm{eq}}$ of the layer can be obtained by using a multiplier on the thermal conductivity of the static medium that corresponds to the Nusselt Number from the Dropkin-Somerscales correlation according to

$$\frac{k_{\mathrm{eq}}}{k_i}=\frac{h_i{\overline{r}}_{i-1}}{k_i}=\mathrm{Nu}=0.049{\left(\text{GrPr}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}{Pr}^{0.074}$$E35

where the Grashof and Prandtl numbers are defined as

$$\mathrm{Gr}=\frac{{\beta }_{\mathrm{f}}g\left({\overline{T}}_{i-1}-{\overline{T}}_i\right){\left({\overline{r}}_{i-1}-{\overline{r}}_i\right)}^3}{{\nu }^2}$$E36

and

$$Pr=\frac{\nu }{\alpha }$$E37

where in the interest of keeping with the traditional nomenclature used in the literature, the term ${\beta }_{\mathrm{f}}$ in Eq. (36) is the coefficient of volumetric thermal expansion (essentially the isobaric volume expansivity), and should not be mistaken for the isothermal compressibility.⁶ Owing to the exponent of one-third in Eq. (35), an iterative procedure is required to evaluate the overall heat transfer coefficient, with Eqs. (32)–(37) all embedded in the iterative loop.

3.4 Shear stress

For flow in a conduit, the frictional resistance is expressed as a shear stress per unit distance in the streamwise gradient that is related to the flow velocity according to the Darcy–Weisbach model

$$\tau \left(z\right)=-\frac{1}{8}\mathit{f\rho }{V}^2$$E38

where the friction factor $f=f\left(Re,\raisebox{1ex}{$\epsilon $}\!\left/ \!\raisebox{-1ex}{$D$}\right.\right)$ can be obtained in terms of the flow Reynolds number and the pipe roughness ($\epsilon$) to diameter ratio, according to the iterative Colebrook-White model or any one of several noniterative approximations published in the literature. Note that the negative sign in Eq. (38) implies that the shear stress acts in the direction opposing the flow.

4. Steady state temperature profiles

Consider the scenario depicted in Figure 6, in which fluid enters a vertical wellbore from a reservoir at a fixed temperature $T_{\mathrm{BH}}$. This temperature is generally referred to as the Static Bottomhole Temperature. The formation temperature is assumed to decrease linearly with depth down to $T_{\text{Surf}}$ at the wellbore exit, such that

Figure 6.

Production through a wellbore from a reservoir.

$$T_{\mathrm{geo}}\left(z\right)=T_{\mathrm{BH}}-\frac{z}{L}\left(T_{\mathrm{BH}}-T_{\text{Surf}}\right)$$E39

Dimensionless streamwise coordinate, and fluid and geothermal temperatures can be defined according to

$$\xi =\frac{z}{L},\theta =\frac{T-T_{\text{Surf}}}{T_{\mathrm{BH}}-T_{\text{Surf}}}\begin{array}{cc},& \end{array}{\theta }_{\mathrm{geo}}=\frac{T_{\mathrm{geo}}-T_{\text{Surf}}}{T_{\mathrm{BH}}-T_{\text{Surf}}}=1-\xi$$E40

Substitution of Eqs. (28) and (38) for the heat flux and shear stress into Eq. (26) and rearranging, results in the compact form

$$\frac{\mathrm{d}\theta }{\mathrm{d}\xi }=N_{\mathrm{TU}}\left({\theta }_{\mathrm{geo}}-\theta \right)+\mathrm{\Lambda }=N_{\mathrm{TU}}\left(1-\xi -\theta \right)+\mathrm{\Lambda }$$E41

where readers familiar with classical heat exchanger analysis will identify the coefficient of the temperature differential as the Number of Transfer Units defined in this context as

$$N_{\mathrm{TU}}=2\frac{\mathit{URL}}{\rho {\mathit{VR}}^2{c}_{\mathrm{p}}}=2\pi \frac{\overline{U}L}{\dot{m}{c}_{\mathrm{p}}}\equiv \frac{1}{\mathrm{Pe}}$$E42

which as noted above can be expressed as the reciprocal of a⁷ Peclet Number. The dimensionless Frictional Heating parameter in Eq. (41) is defined as

$$\mathrm{\Lambda }=\frac{f}{4}\left(\frac{V^2}{c_{\mathrm{p}}\left(T_{\mathrm{BH}}-T_{\text{Surf}}\right)}\right)\frac{L}{R}$$E43

Note that the dimensionless entity in parenthesis within the expression for $\mathrm{\Lambda }$ is the Eckert Number. If the temperature at the inlet to the wellbore is the same as the reservoir temperature $T_{\mathrm{BH}}$, the boundary condition corresponding to Eq. (41) is

$$\theta \left(0\right)=1$$E44

Subject to the boundary condition above, the solution of Eq. (41) is

$$\theta \left(\xi \right)=1-\xi +{\mathit{Ce}}^{-N_{\mathrm{TU}}\xi }+\frac{1+\mathrm{\Lambda }}{N_{\mathrm{TU}}}$$E45

For $N_{\mathrm{TU}}=1$, the dimensionless temperature profiles are plotted in Figure 7 for various values of the frictional heating parameter $\mathrm{\Lambda }$. What is noteworthy is that as $\mathrm{\Lambda }$ increases from zero, the temperature at the surface (known as the Arrival Temperature) not only approaches the reservoir temperature, but in fact exceeds it, a very common observation in prolific deepwater oil producers. Neglecting the frictional heating term can therefore result in a severe underprediction of temperatures and threaten wellbore integrity if the attendant tubular thermal stresses and annular pressure buildup are accordingly underpredicted.

Figure 7.

Temperature profiles for N_TU = 1 and various values of $\mathrm{\Lambda }$.

The impact of the Number of Transfer Units is shown in Figure 8, and shows that even in the absence of frictional heating, near-isothermal conditions in the wellbore can be achieved as the Number of Transfer Units number becomes very small. In fact, in the limit $N_{\mathrm{TU}}\to 0$, the RHS of Eq. (41) reduces to zero for $\mathrm{\Lambda }$ = 0 implying $\theta \left(\xi \right)=1$ throughout in accordance with the boundary condition of Eq. (44).

Figure 8.

Effect of N_TU on temperature profiles for Λ = 0 (left) and Λ = 0.5 (right).

A word of caution is in order here. The phenomenon of the arrival temperature exceeding the bottomhole temperature as evidenced in Figure 7 is associated only with liquids that almost always have a Negative Joule-Thomson coefficient $c_{\mathrm{JT}}$. For gases, $c_{\mathrm{JT}}$ is negative only below the inversion pressure. Accordingly, for very high gas rate flows in wellbores, it is not uncommon to see the contrary effect of a substantial drop in the fluid temperature towards the surface. Attempting to simulate this effect however, with negative values of $\mathrm{\Lambda }$ (positive $c_{\mathrm{JT}}$) will yield results which while seemingly plausible may not be accurate since Eq. (41) is only valid for incompressible flows, and the assumptions invoked in its derivation tend to break down when the produced fluid is predominantly gaseous.

5. Transient heat transfer in wellbores

Thermal transients in a wellbore are characterised by the fluid exchanging heat with the surrounding formation at a rate that evolves in time. Therefore, there are two adjacent coupled problems that need to be considered – the transient transport equation in the wellbore conduit of radius $R$ and the transient radial diffusion in the formation. These two problems are coupled at the interface between the outer layer of the wellbore and the formation at the radius $\overline{R}$ as shown in Figure 9. Between the radial locations $R$ and $\overline{R}$ are all of the wellbore layers comprised of tubulars and annuli. For the purpose of this illustrative example, it will be assumed that these layers have negligible capacitance, so that they respond instantaneously to the fluid transients.⁸ Unlike in the steady state case where the overall heat transfer coefficient across the wellbore layers was used to link the fluid temperature to the undisturbed geothermal temperature, in a transient analysis, the linkage is between the fluid temperature and the transient temperature at the wellbore-formation interface according to

Figure 9.

Transient production through a wellbore.

$$\dot{q}\left(z,t\right)=-U\left[T\left(z,t\right)-T_{\mathrm{WBF}}\left(z,t\right)\right]=k_{\mathrm{Geo}}{\left.\frac{\partial T_{\text{Form}}}{\partial r}\right|}_{r=\overline{R}}$$E46

With the assumption of an incompressible fluid, the transient transport equation can be extracted from Eq. (20) as

$$\frac{\partial T}{\partial t}+V\frac{\partial T}{\partial z}=\frac{Q-\mathit{VF}}{\rho c_{\mathrm{p}}}$$E47

subject to the same boundary condition as in the steady state case i.e., $T\left(0,t\right)=T_{\mathrm{BH}}$ and the initial condition

$$T\left(z,0\right)=T_{\mathrm{geo}}\left(z\right)$$E48

The temperature at the interface between the wellbore and the formation is not known a-priori, but constitutes one of the radial boundary conditions for the problem of diffusion in the formation. In lieu of a semi-infinite domain, we will regard the formation as a finite cylindrical domain with an outer radius far enough that geothermal conditions prevail therein at the end of the well operational time period of interest. Accordingly, the diffusion in the formation is governed by the partial differential equation

$${\rho }_{\mathrm{geo}}{c}_{{\mathrm{p}}_{\mathrm{geo}}}\frac{\partial T_{\text{Form}}}{\partial t}=k_{\mathrm{geo}}\frac{{\partial }^2{T}_{\text{Form}}}{\partial r^2}$$E49

subject to the initial condition

$$T_{\text{Form}}\left(r,0,z\right)=T_{\mathrm{Geo}}\left(z\right)$$E50

and the boundary conditions

$$T_{\text{Form}}\left(\overline{R},\tau ,z\right)=T_{\mathrm{WBF}}\left(z\right)$$E51

and

$$\frac{\partial T_{\text{Form}}}{\partial r}\left({\overline{R}}_{\infty },\tau ,z\right)=0$$E52

Implicit in Eq. (49) is the assumption that axial diffusion is negligible, which given the length scales of typical wellbores, is eminently justified. As a consequence, the governing equation holds at all depths along the wellbore where the thermal interaction between the wellbore and the formation is described in terms of a purely radial heat transfer mechanism.

The time and spatial coordinates variables are non-dimensionalised as

$$\xi =\frac{z}{L}\begin{array}{cc},& \end{array}\eta =\frac{r}{\overline{R}}\begin{array}{cc},& \end{array}{\eta }_{\infty }=\frac{R_{\infty }}{\overline{R}}\begin{array}{cc},& \end{array}\tau =\frac{{\alpha }_{\mathrm{geo}}}{{\overline{R}}^2}t$$E53

where the dimensionless time is essentially a Fourier Number, and the temperatures are normalised as before with respect to the bottomhole to surface geothermal temperature difference as

$$\begin{array}{l}\theta =\frac{T-T_{\text{Surf}}}{T_{\mathrm{BH}}-T_{\text{Surf}}}\begin{array}{cccc}& & & \end{array}{\theta }_{\text{form}}=\frac{T_{\text{Form}}-T_{\text{Surf}}}{T_{\mathrm{BH}}-T_{\text{Surf}}}\\ \psi =\frac{T_{\mathrm{WBF}}-T_{\text{Surf}}}{T_{\mathrm{BH}}-T_{\text{Surf}}}\begin{array}{cccc}& & & \end{array}\varphi =\frac{T_{\text{Form}}-T_{\text{Surf}}}{T_{\mathrm{BH}}-T_{\text{Surf}}}\end{array}$$E54

Substitution of the modified heat flux constitutive model i.e., Eq. (46), the shear stress model from Eq. (38) and the set of dimensionless variables defined by Eqs. (53) and (54) into Eq. (47) and Eqs. (48)–(52) result in the following dimensionless system of coupled equations:

$$\frac{\partial \theta }{\partial \tau }+\mathrm{Pe}\frac{\partial \theta }{\partial \xi }=\mathrm{\Gamma }\left(\psi -\theta \right)+\mathrm{\Lambda }$$E55

where

$$\mathrm{Pe}=\frac{V\overline{R}}{{\alpha }_{\mathrm{geo}}}\frac{\overline{R}}{L}$$E56

is a Peclet number. In the context of the steady state problem, the term

$$\mathrm{\Gamma }=2\frac{\overline{R}}{R}\frac{U\overline{R}}{\rho c_{\mathrm{p}}{\alpha }_{\mathrm{geo}}}$$E57

is a diffusion coefficient and the frictional heating parameter

$$\mathrm{\Lambda }=\frac{f}{4}\frac{\overline{R}}{R}\frac{V^3\overline{R}}{c_{\mathrm{p}}{\alpha }_{\mathrm{geo}}\left(T_{\mathrm{BH}}-T_{\text{Surf}}\right)}$$E58

as defined above is somewhat different from that described earlier. Eq. (55) is subject to the initial condition

$$\theta \left(\xi ,0\right)=1-\xi$$E59

and the boundary condition

$$\theta \left(0,\tau \right)=1$$E60

which must be solved in conjunction with the formation diffusion problem non-dimensionalised as

$$\frac{\partial \varphi }{\partial \tau }=\frac{{\partial }^2\varphi }{\partial {\eta }^2}$$E61

subject to the initial condition

$$\varphi \left(\eta ,\xi ,0\right)=1-\xi$$E62

along with the a priori unknown boundary condition at the wellbore formation interface

$$\varphi \left(1,\xi ,\tau \right)=\psi \left(\xi ,\tau \right)$$E63

and the farfield boundary condition

$$\frac{\partial \varphi }{\partial \eta }\left({\eta }_{\infty },\xi ,\tau \right)=0$$E64

The value of the farfield radius ratio ${\eta }_{\infty }$ must be choses so as to be consistent with the physics of the problem. While the Neumann condition (zero flux) in Eq. (64) at the drainage radius is by itself adequate from a mathematical standpoint to provide closure to the system of equations, physical realism also requires that the formation temperature asymptotically approach the undisturbed geothermal temperature at the depth in question, at a radial location prior to the drainage radius. Failure to satisfy this criterion due to an insufficiently large value of ${\eta }_{\infty }$ could result in substantially inaccurate calculations. A good rule of thumb for estimating the required drainage ratio is ${\eta }_{\infty }=5\sqrt{\tau }$ where the Fourier number corresponds to the end of the time period of interest.

5.1 Solution of the transient formation diffusion problem

We start with the solution of Eq. (61) subject to the criteria of Eqs. (62)–(64) that following Ozisik [15], involves the use of Duhamel’s Theorem as is customary for problems involving time dependent boundary conditions. The radial temperature profile is not the actual entity of interest. What is necessary to facilitate the coupling of the formation diffusion problem with the fluid transport equation, is the interface flux wherein flux continuity as expressed by Eq. (46) yields in terms of non-dimensional entities

$$\psi -\theta =\gamma {\left.\frac{\mathrm{d}\varphi }{\mathrm{d}\eta }\left(\xi ,\tau \right)\right|}_{\eta =1}=\gamma \sum _{j=1}^{\infty }{\overline{C}}_j{D}_j\left(\xi ,\tau \right)$$E65

where

$$\gamma =\frac{k_{\mathrm{Geo}}}{U\overline{R}}$$E66

and the Duhamel Convolution Integral is

$$D_j\left(\xi ,\tau \right)=-\underset{0}{\overset{\tau }{\int }}{e}^{-{\lambda }_j^2\left(\tau -\beta \right)}{\psi }^{\prime }\left(\xi ,\beta \right)\mathrm{d}\beta$$E67

where the prime denotes a derivative with respect to $\beta$. The eigenvalues ${\lambda }_j$ and the Fourier-Bessel coefficients ${\overline{C}}_j$ are defined in the appendix and depend on the farfield radius ratio ${\eta }_{\infty }$. In practice, the infinite summation in Eq. (65) is obviously restricted to a finite number of Fourier modes. Note that the purpose of the exercise above was to express the temperature difference $\psi -\theta$ in terms of an interface flux.

5.2 Solution of the transient fluid transport equation

As is the case with a large number of problems in transient heat transfer, the first step in the solution of Eq. (55) is Laplace transformation whereupon we have

$$s\mathrm{\Theta }-\theta \left(\xi ,0\right)+\mathrm{Pe}\frac{\mathrm{d\Theta }}{\mathrm{d}\xi }=\mathrm{\Gamma }\left(\mathrm{\Psi }-\mathrm{\Theta }\right)+\frac{\mathrm{\Lambda }}{s}$$E68

Laplace transformation of the interface flux expression of Eq. (65) in conjunction with Eq. (67) and some algebra yields

$$\mathrm{\Psi }-\mathrm{\Theta }=-\gamma \sum _{j=1}^{\infty }{\overline{C}}_j\left[\frac{s\mathrm{\Psi }-\psi \left(0,\xi \right)}{s+{\lambda }_j^2}\right]$$E69

where the Convolution Theorem has been used on the Duhamel Integral. Noting that the initial condition or the interface flux is the undisturbed geothermal temperature, Eq. (69) can be rearranged (see [3]) into the compact form

$$\mathrm{\Psi }-\mathrm{\Theta }=G\left(s\right)\left(1-\xi \right)-\mathit{sG}\left(s\right)\mathrm{\Theta }$$E70

where

$$G\left(s\right)=\frac{W\left(s\right)}{1+\mathit{sW}\left(s\right)}$$E71

and

$$W\left(s\right)=\gamma \sum _{j=1}^{\infty }\frac{{\overline{C}}_j}{s+{\lambda }_j^2}$$E72

Substitution of Eq. (70) into Eq. (68) along with the initial condition of Eq. (59) results in the ordinary differential equation in the frequency domain

$$\frac{\mathrm{d\Theta }}{\mathrm{d}\xi }=\left(\frac{1+\mathrm{\Gamma }G}{\mathrm{Pe}}\right)\left(1-\xi \right)-s\left(\frac{1+\mathrm{\Gamma }G}{\mathrm{Pe}}\right)\mathrm{\Theta }+\left(\frac{\mathrm{\Lambda }}{\mathrm{Pe}}\right)s^{-1}$$E73

subject to the transformed boundary condition

$$\mathrm{\Theta }\left(0,s\right)=\frac{1}{s}$$E74

The solution of Eq. (73) in the frequency domain is

$$\mathrm{\Theta }\left(\xi ,s\right)=F\left(s\right)\left(1-e^{-A\left(s\right)\xi }\right)+\frac{1}{s}\left(e^{-A\left(s\right)\xi }-\xi \right)$$E75

where

$$F\left(s\right)=\frac{1}{s}+\left(\frac{\mathrm{Pe}+\mathrm{\Lambda }}{1+\mathrm{\Gamma }G\left(s\right)}\right)\frac{1}{s^2}$$E76

and

$$A\left(s\right)=s\left(\frac{1+\mathrm{\Gamma }G\left(s\right)}{\mathrm{Pe}}\right)$$E77

5.3 Inversion from the frequency domain

The solution expressed in Eq. (75) in the frequency domain is not particularly useful from a practical point of view. It must therefore be inverted back to physical space using the Inverse Laplace Transform. One approach is to use the Cauchy Residue Theorem by summing the residues over all of the poles of Eq. (75). One of these poles is at the Origin.⁹ The remaining poles lie along the negative real axis of the complex plane and are the zeros of the denominator of the 2nd term in Eq. (76) such that

$$f\left(s\right)=1+\mathrm{\Gamma }G\left(s\right)=0$$E78

which must be solved numerically. An efficient method of doing so involves an asymptotic bracketing technique, the description of which is outside of the scope of this chapter.

The Residue Theorem while attractive from the standpoint of constituting a formal analytical solution can involve some very tedious if otherwise straightforward book-keeping in addition to the requirement of numerical evaluation of the roots of Eq. (78). A far more efficient approach is to use numerical inversion with the Gaver-Stehfest Function-Sampling Algorithm (Stehfest [16]) whereupon the temperatures in physical space are given by

$$\theta \left(\xi ,\tau \right)\approx -\frac{ln2}{\tau }\sum _{k=1}^{2N_{\mathrm{G}}}{\sigma }_k\mathrm{\Theta }\left(\frac{kln2}{\tau }\right)$$E79

where $N_{\mathrm{GS}}$ is the (even) order of the Gaver Summation, and ${\sigma }_k,k=1\cdots N_{\mathrm{GS}}$ are the Stehfest Accelerators, defined and listed in Table 1 in the Appendix.

The evolution of the temperature profiles for two cases –with and without frictional heating, is shown in Figure 10. In both cases, the Gaver-Stehfest function sampling algorithm was used in conjunction with a farfield radius ratio of 200 and 1000 Fourier modes. The Negative Joule-Thomson effect is clearly seen in the right panel of the figure for the case with frictional heating.

Figure 10.

Evolution of the transient temperature profiles for Λ = 0 (left) and Λ = 0.05 (right).

6. Heat transfer in circulating scenarios

Circulation constitutes an important aspect of wellbore operations. This is most commonly encountered in drilling, as well as swapping fluids and hole cleaning. In what is known as forward circulation, fluid is pumped down a drillstring and returns to the surface through the annulus, as depicted in panel (c) of Figure 3. In Reverse Circulation (panel d Figure 3), the flow directions are reversed so that colder fluid is injected down the annulus and hotter fluid returns to the surface. This is not that common in conventional wellbores, but is the primary mode of operation in geothermal wells where hot fluid returns through an insulated or partially insulated inner string known in that context as the tubing.

The thermal interactions in both scenarios are depicted in Figure 11 indicating the known boundary conditions and the a priori unknown return fluid temperatures of interest. At the bottom of the well, a matching condition stipulates that the pipe and annulus temperatures (designated by the subscripts “p” and “a”, respectively) are equal. In the simple well configuration considered, there is an inner pipe of inner radius $R$. The annulus has inner and outer radii $R_{\mathrm{i}}$ and $R_{\mathrm{o}}$ as indicated. The outer casing is cemented with the outer radius of the wellbore $R_{\mathrm{WB}}$ in contact with the formation. A linear geothermal gradient is assumed.

Figure 11.

Thermal interaction in forward (left) and reverse (right) circulation scenarios.

Invoking the previous assumptions of both incompressibility and pseudo steady state heat transfer, the governing equations describing both forward and reverse circulation can be described by the pair of equations

$$\delta \frac{\mathrm{d}{T}_{\mathrm{p}}}{\mathrm{d}z}=\frac{Q_{\mathrm{p}}-V_{\mathrm{p}}{F}_{\mathrm{p}}}{\rho {\mathit{cV}}_{\mathrm{p}}}$$E80

for flow in the pipe (denoted by the subscript “p”), and

$$-\delta \frac{\mathrm{d}{T}_{\mathrm{a}}}{\mathrm{d}z}=\frac{Q_{\mathrm{a}}-V_{\mathrm{a}}{F}_{\mathrm{a}}}{\rho {\mathit{cV}}_{\mathrm{a}}}$$E81

for flow in the annulus (denoted by the subscript “a”). Note that the subscript “p” has been dropped from the specific heat since it is now used to denote the flowing stream in the inner pipe. The direction of the circulation is characterised by the parameter $\delta =\pm 1$, with the positive and negative signs denoting Forward and Reverse circulation, respectively. The associated boundary conditions are

$$\left(1+\delta \right)T_{\mathrm{p}}\left(0\right)+\left(1-\delta \right)T_{\mathrm{a}}\left(0\right)=2T_{\text{in}}$$E82

at the surface inlet and the matching condition

$$T_{\mathrm{a}}\left(L\right)=T_{\mathrm{p}}\left(L\right)$$E83

at the bottom of the well assuming no losses as the fluid leaves one conduit and enters another. The velocities are related through mass conservation such that

$$\dot{m}=\rho A_{\mathrm{p}}{V}_{\mathrm{p}}=\rho A_{\mathrm{a}}{V}_{\mathrm{a}}\Rightarrow V_{\mathrm{a}}=\left(\frac{A_{\mathrm{p}}}{A_{\mathrm{a}}}\right)V_{\mathrm{p}}$$E84

The Heat Transfer terms in Eqs. (80) and (81) are

$$Q_{\mathrm{p}}=2R^{-2}{\overline{U}}_{\mathrm{pa}}\left(T_{\mathrm{a}}-T_{\mathrm{p}}\right)$$E85

for the pipe describing the interaction between the two flowing streams and

$$Q_{\mathrm{a}}=2R^{-2}{\overline{U}}_{\mathrm{pa}}\left(T_{\mathrm{p}}-T_{\mathrm{a}}\right)+2R^{-2}{\overline{U}}_{\mathrm{a}\infty }\left(T_{\mathrm{Geo}}-T_{\mathrm{a}}\right)$$E86

for the annulus describing the interaction between the streams and the interaction between the annulus stream and the formation. In accordance with the formulation of the overall heat transfer conductance, we have with respect to the geometry of Figure 11 the following expressions

$$\frac{{\overline{U}}_{\mathrm{pa}}}{k_{\text{fluid}}}={\left[\frac{k_{\text{fluid}}}{\mathit{hR}}+\frac{k_{\text{fluid}}}{k_{\text{steel}}}ln\frac{R_{\mathrm{i}}}{R}+\frac{k_{\text{fluid}}}{h_{\mathrm{i}}{R}_{\mathrm{i}}}\right]}^{-1}$$E87

and

$$\frac{{\overline{U}}_{\mathrm{a}\infty }}{k_{\text{fluid}}}={\left[\frac{k_{\text{fluid}}}{h_{\mathrm{o}}{R}_{\mathrm{o}}}+\frac{k_{\text{fluid}}}{k_{\text{steel}}}ln\frac{R_{\mathrm{o}}}{R_{\mathrm{C}}}+\frac{k_{\text{fluid}}}{k_{\text{cement}}}ln\frac{R_{\mathrm{C}}}{R_{\mathrm{wb}}}\frac{k_{\text{fluid}}}{k_{\mathrm{Geo}}}\frac{1}{F\left(\tau \right)}\right]}^{-1}$$E88

where $\tau$ is the Fourier number corresponding to the instantaneous snapshot in time at which the temperature profiles correspond to the pseudo steady state. In both expressions above, it should be noted that the division by the constant fluid thermal conductivity (assumed) enables the evaluation of the forced convection Nusselt numbers in terms of known correlations such as the Dittus-Boelter or Sieder-Tate models. The frictional heating terms for the pipe and annulus streams are given by

$$F_{\mathrm{p}}=2\frac{{\tau }_{\mathrm{p}}\left(z\right)}{R}=\frac{f_{\mathrm{p}}}{4}\frac{\rho V_{\mathrm{p}}^2}{R}$$E89

and

$$F_{\mathrm{a}}=2\frac{{\tau }_{\mathrm{a}}\left(z\right)}{R_{\mathrm{o}}-R_{\mathrm{i}}}=\frac{f_{\mathrm{a}}}{4}\frac{\rho V_{\mathrm{a}}^2}{\left(R_{\mathrm{o}}-R_{\mathrm{i}}\right)}$$E90

Normalising the wellbore streamwise coordinate by the length as in the previous exercises, and the pipe and annulus temperatures by the surface to well depth temperature difference as before yields the following coupled system of equations

$$\frac{\mathrm{d}}{\mathrm{d}\xi }\left[\begin{array}{c}{\theta }_{\mathrm{p}}\\ {\theta }_{\mathrm{a}}\end{array}\right]=\left[\begin{array}{cc}-\delta N_{\mathrm{pa}}& \delta N_{\mathrm{pa}}\\ -\delta N_{\mathrm{pa}}& \delta \left(N_{\mathrm{pa}}+N_{\mathrm{a}\infty }\right)\end{array}\right]\left[\begin{array}{c}{\theta }_{\mathrm{p}}\\ {\theta }_{\mathrm{a}}\end{array}\right]+\left[\begin{array}{c}\delta {\mathrm{\Lambda }}_{\mathrm{p}}\\ -\delta {\mathrm{\Lambda }}_{\mathrm{a}}\end{array}\right]+\left[\begin{array}{c}0\\ -\delta N_{\mathrm{a}\infty }\end{array}\right]\xi$$E91

which makes use of the fact that the normalised linear geothermal temperature is ${\theta }_{\mathrm{Geo}}\left(\xi \right)=\xi$. Eq. (91) is subject to the pair of boundary conditions

$$\begin{array}{l}\left(1+\delta \right){\theta }_{\mathrm{p}}\left(0\right)+\left(1-\delta \right){\theta }_{\mathrm{a}}\left(0\right)=2{\theta }_{\text{in}}\\ {\theta }_{\mathrm{a}}\left(1\right)={\theta }_{\mathrm{p}}\left(1\right)\end{array}$$E92

The governing dimensionless parameters are the Number of Transfer Unit parameters

$$N_{\mathrm{pa}}=2\pi \frac{{\overline{U}}_{\mathrm{pa}}L}{\dot{m}c}\begin{array}{ccc}& & \end{array}{N}_{\mathrm{a}\infty }=2\pi \frac{{\overline{U}}_{\mathrm{a}\infty }L}{\dot{m}c}$$E93

and the dimensionless frictional heating parameters

$${\mathrm{\Lambda }}_{\mathrm{p}}=\frac{f_{\mathrm{p}}}{4}\frac{V_{\mathrm{p}}^2}{c\mathrm{\Delta }T}\frac{L}{R}\begin{array}{ccc}& & \end{array}{\mathrm{\Lambda }}_{\mathrm{a}}=\frac{f_{\mathrm{a}}}{4}\frac{V_{\mathrm{a}}^2}{c\mathrm{\Delta }T}\frac{L}{R_{\mathrm{o}}-R_{\mathrm{i}}}$$E94

The analytical solution of Eq. (91) yields the pair of equations for the pipe and annulus temperature profiles as

$${\theta }_{\mathrm{p}}\left(\xi \right)=C_{\lambda }{e}^{\mathit{\lambda \xi }}+C_{\mu }{e}^{\mathit{\mu \xi }}+h_{\mathrm{p}}+\xi$$E95

and

$${\theta }_{\mathrm{a}}\left(\xi \right)=C_{\lambda }{r}_{\lambda }{e}^{\mathit{\lambda \xi }}+C_{\mu }{r}_{\mu }{e}^{\mathit{\mu \xi }}+h_{\mathrm{a}}+\xi$$E96

where $\lambda$ and $\mu$ are the eigenvalues of the matrix in Eq. (91) given by

$$\begin{array}{l}\lambda =\frac{1}{2}\left[\delta N_{\mathrm{a}\infty }+\sqrt{{N_{\mathrm{a}\infty }}^2+4N_{\mathrm{pa}}{N}_{\mathrm{a}\infty }}\right]\\ \mu =\frac{1}{2}\left[\delta N_{\mathrm{a}\infty }-\sqrt{{N_{\mathrm{a}\infty }}^2+4N_{\mathrm{pa}}{N}_{\mathrm{a}\infty }}\right]\end{array}$$E97

The $r$ and $h$ constants in Eqs. (95) and (96) are

$$\begin{array}{l}{r}_{\lambda }=1+\lambda {\left(\delta N_{\mathrm{pa}}\right)}^{-1}\\ r_{\mu }=1+\mu {\left(\delta N_{\mathrm{pa}}\right)}^{-1}\end{array}$$E98

and

$$\begin{array}{l}{h}_{\mathrm{p}}={\left(N_{\mathrm{pa}}{N}_{\mathrm{a}\infty }\right)}^{-1}\left[\left(N_{\mathrm{pa}}+N_{\mathrm{a}\infty }\right){\mathrm{\Lambda }}_{\mathrm{p}}+N_{\mathrm{pa}}{\mathrm{\Lambda }}_{\mathrm{a}}-\delta N_{\mathrm{a}\infty }\right]\\ h_{\mathrm{a}}={N_{\mathrm{a}\infty }}^{-1}\left({\mathrm{\Lambda }}_{\mathrm{p}}+{\mathrm{\Lambda }}_{\mathrm{a}}\right)\end{array}$$E99

The constants of integration are determined from the boundary conditions as

$$\begin{array}{l}{C}_{\lambda }=\frac{\left(e^{\mu }-r_{\mu }{e}^{\mu }\right)\left(2{\theta }_{\text{in}}-\left(1+\delta \right)p_{\mathrm{p}}-\left(1-\delta \right)p_{\mathrm{a}}\right)-\left(p_{\mathrm{a}}-p_{\mathrm{p}}\right)\left[\left(1+\delta \right)+\left(1-\delta \right)r_{\mu }\right]}{\left[\left(1+\delta \right)+\left(1-\delta \right)r_{\lambda }\right]\left(e^{\mu }-r_{\mu }{e}^{\mu }\right)-\left[\left(1+\delta \right)+\left(1-\delta \right)r_{\mu }\right]\left(r_{\lambda }{e}^{\lambda }-e^{\lambda }\right)}\\ \\ C_{\mu }=\frac{\left(r_{\lambda }{e}^{\lambda }-e^{\lambda }\right)\left(2{\theta }_{\text{in}}-\left(1+\delta \right)p_{\mathrm{p}}-\left(1-\delta \right)p_{\mathrm{a}}\right)+\left(p_{\mathrm{a}}-p_{\mathrm{p}}\right)\left[\left(1+\delta \right)+\left(1-\delta \right)r_{\lambda }\right]}{\left[\left(1+\delta \right)+\left(1-\delta \right)r_{\lambda }\right]\left(e^{\mu }-r_{\mu }{e}^{\mu }\right)-\left[\left(1+\delta \right)+\left(1-\delta \right)r_{\mu }\right]\left(r_{\lambda }{e}^{\lambda }-e^{\lambda }\right)}\end{array}$$E100

6.1 Forced convection in the annulus

In Eqs. (87) and (88), the terms $h_{\mathrm{i}}$ and $h_{\mathrm{o}}$ represent the heat transfer coefficients at the inner and outer surfaces of the annulus, respectively. In turbulent flows in annuli with radius ratios approaching unity, the following approximation can be used

$$h_{\mathrm{i}}=h_{\mathrm{o}}=\overline{h}=\frac{k_{\text{fluid}}}{{\overline{D}}_{\mathrm{hyd}}}{\mathrm{Nu}}_{\mathrm{T}}=\frac{k_{\text{fluid}}}{{\overline{D}}_{\mathrm{hyd}}}C{Re}^m{Pr}^n$$E101

where ${\overline{D}}_{\mathrm{hyd}}$ is the hydraulic diameter of the annulus, and $C$, $m$, and $n$ are the constants of the forced convection correlation used.¹⁰ Most annular flows in wellbore circulating scenarios however, tend to be laminar, and associated with annulus radius ratios often well less than unity. In addition, the Non-Newtonian nature of the flow must be considered. Merely replacing the turbulent Nusselt Number ${\mathrm{Nu}}_{\mathrm{T}}$ in Eq. (101) with its laminar analogue ${\mathrm{Nu}}_{\mathrm{L}}$ is not consistent with the physics of the problem. It is recommended that for fluids that obey the Power-Law model, the following correlations from Chandrasekhar [17] be used instead

$$\begin{array}{l}\frac{h_{\mathrm{i}}\overline{D}}{k_{\text{fluid}}}\left(n,\kappa \right)=-\frac{2}{1-{\theta }_{\mathrm{b}}\left(n,\kappa \right)}\left(\frac{1-\kappa }{\kappa ln\kappa }\right)\\ \frac{h_{\mathrm{o}}\overline{D}}{k_{\text{fluid}}}\left(n,\kappa \right)=-\frac{2}{{\theta }_{\mathrm{b}}\left(n,\kappa \right)}\left(\frac{1-\kappa }{ln\kappa }\right)\end{array}$$E102

where the dimensionless bulk temperature is a function of the power law index and radius ratio, and is given by

$$\begin{array}{l}{\theta }_{\mathrm{b}}\left(n,\kappa \right)=\sum _{j=1}^4\left(a_j+{\mathit{nb}}_j\right){\kappa }^{j-1}\\ a_1=0.213\begin{array}{cc},& a_2=0.576\begin{array}{cc},& a_3=-0.439\begin{array}{cc},& a_4=0.152\end{array}\end{array}\end{array}\\ b_1=0.0043\begin{array}{cc},& b_2=-0.0183\begin{array}{cc},& b_3=0.0236\begin{array}{cc},& b_4=-0.0102\end{array}\end{array}\end{array}\end{array}$$E103

6.2 Examples of forward and reverse circulation

For a given set of operational parameters, the intermediate calculations and evaluation of the dimensionless parameters is shown in Figure 12 which represents a case of forward circulation of an oil-based fluid in a fairly typical drilling scenario.

Figure 12.

Estimation of dimension groups from problem data (forward circulation).

The circulating temperature profiles in the Drillpipe and its annulus are shown in Figure 13 for the parameters in Figure 12, but three different flowrates. It is seen that the temperature at the well TD at first decreases with flowrate as would be expected, but at higher flowrates, actually increases due to frictional heating. This is more clearly evident in Figure 14 where the TD and arrival temperatures are plotted over a range of flowrates. The inflexion point in the well TD temperature is where the negative Joule-Thomson effect surpasses the advection effect in the drillpipe. The point at which the (dashed) arrival temperature curve intersects the (solid) well TD curve is where frictional heating is significant even in the annulus.

Figure 13.

Circulating temperature profiles for three different flowrates – Forward circulation of an oil-based drilling fluid.

Figure 14.

Negative joule-Thomson (frictional heating) effect of circulating flow rate on well temperatures.

The temperature profiles shown in Figure 15 for three different mass flow rates of water in a reverse circulating scenario correspond to a geothermal well. In this case, the inner conduit known as the tubing is assumed to be insulated as is common in geothermal wells. The well TD and arrival temperatures quickly tend to become independent of mass flow rate. The thermal power produced by a geothermal well is given by

Figure 15.

Circulating temperature profiles for three different flowrates – Reverse circulation of water in a geothermal well with an insulated tubing.

$${\dot{P}}_{\mathrm{MWT}}=\dot{m}{c}_{\mathrm{P}}\left(T_{\mathrm{arr}}-T_{\text{inlet}}\right)$$E104

and is plotted as a function of mass flow rate as shown in Figure 16. The arrival temperature is seen to become independent of mass flow rate at about 25 kg/sec whereupon the thermal power increases linearly. It is important to note that only a fraction of the thermal power is actually converted into electric power (often at a rate of about 15–20%) that can be transmitted to a grid. Furthermore, a portion of even this converted power has to be used to overcome the parasitic power due to frictional losses in the geothermal wellbore. The thermal power serves however as a useful metric in a parametric sensitivity analysis of a geothermal well.

Figure 16.

Effect of flow rate on thermal power generated and arrival temperature in a geothermal well with an insulated tubing.

The circulating thermal model developed here considers only a simple monobore well for illustrative purposes. For an extension of the methodology to a complex wellbore with multiple wellbore segments coupled via an arbitrary number of interface temperature matching conditions, the reader is referred to [9]. That study also considers curvature and tortuosity effects in deviated wells, variable lithology, multiple geothermal gradients, and the effects of fluid rheology.

7. Evaluation of wellbore and interface temperatures

7.1 Wellbore temperatures in a transient analysis

The methodologies described in the previous sections dealt with the estimation of flowing temperatures in operating and circulating scenarios. An issue of equal - if not often greater - importance is the estimation of temperatures in the intervening wellbore layers (fluid and solid) between the flow conduit and the formation. The description of a fully transient analysis wherein the transient temperatures in all layers are evaluated in tandem with the transient flowing temperature is beyond the scope of this chapter, but the interested reader is referred to [3] where such an analysis is described in near-exhaustive granularity. What will be demonstrated in what follows here is how to estimate the interface thermal conductivities and temperatures which are needed to evaluate the heat fluxes as well as the natural convection multipliers required for the estimation of the nodal thermal conductivities in fluid layers.

Consider the depiction in Figure 17 showing 3 adjacent layers designated $i-1$, $i$ and $i+1$ from left to right. The barred and unbarred symbols refer to interfacial and nodal entities respectively. At the interface $i-1$ at the left of the layer $i$, the flux expressed in terms of the straddling nodal difference and the interface conductivity, can also be expressed in terms of the differences between the nodes and the interface, and the nodal conductivities. This relationship can be expressed as

Figure 17.

Nodal and Interface temperatures and conductivities.

$${\dot{\overline{q}}}_{i-1}\left(z\right)=-\frac{{\overline{k}}_{i-1}}{ln\raisebox{1ex}{$r_i$}\!\left/ \!\raisebox{-1ex}{$r_{i-1}$}\right.}\left(T_{i-1}-T_i\right)=-\frac{k_{i-1}}{ln\raisebox{1ex}{${\overline{r}}_{i-1}$}\!\left/ \!\raisebox{-1ex}{$r_{i-1}$}\right.}\left(T_{i-1}-{\overline{T}}_{i-1}\right)=-\frac{k_i}{ln\raisebox{1ex}{$r_i$}\!\left/ \!\raisebox{-1ex}{${\overline{r}}_{i-1}$}\right.}\left({\overline{T}}_{i-1}-T_i\right)$$E105

from which the interface thermal conductivity at the right and left interfaces with indices $i$ and $i-1$ can be expressed in terms of the nodal values as the weighted harmonic means

$$\begin{array}{l}{\overline{k}}_i=ln\raisebox{1ex}{$r_{i+1}$}\!\left/ \!\raisebox{-1ex}{$r_i$}\right.{\left(\frac{{\lambda }_i}{k_i}+\frac{{\mu }_i}{k_{i+1}}\right)}^{-1}\\ {\overline{k}}_{i-1}=ln\raisebox{1ex}{$r_i$}\!\left/ \!\raisebox{-1ex}{$r_{i-1}$}\right.{\left(\frac{{\lambda }_{i-1}}{k_{i-1}}+\frac{{\mu }_{i-1}}{k_i}\right)}^{-1}\end{array}$$E106

where the geometric coefficients are

$${\lambda }_i=ln\raisebox{1ex}{${\overline{r}}_i$}\!\left/ \!\raisebox{-1ex}{$r_i$}\right.\begin{array}{cccc}& & & \end{array}{\mu }_i=ln\raisebox{1ex}{$r_{i+1}$}\!\left/ \!\raisebox{-1ex}{${\overline{r}}_i$}\right.$$E107

The interface temperatures can be expressed in terms of the geometric coefficients in Eq. (107) and the nodal conductivities as

$$\begin{array}{l}{\overline{T}}_i=\frac{{\lambda }_i^{-1}{k}_i{T}_i+{\mu }_i^{-1}{k}_{i+1}{T}_{i+1}}{{\lambda }_i^{-1}{k}_i+{\mu }_i^{-1}{k}_{i+1}}\\ {\overline{T}}_{i-1}=\frac{{\lambda }_{i-1}^{-1}{k}_{i-1}{T}_{i-1}+{\mu }_{i-1}^{-1}{k}_i{T}_i}{{\lambda }_{i-1}^{-1}{k}_{i-1}+{\mu }_{i-1}^{-1}{k}_i}\end{array}$$E108

When there is at least one fluid layer subject to natural convection, the evaluation of the interface values must be embedded in an iterative sequence within each time step. Note that in a transient analysis, the fluxes on either side of a nodal layer need not be equal, so that in general, ${\overline{q}}_{i-1}\ne {\overline{q}}_i$ with respect to Figure 17.

7.2 Wellbore temperatures in a Pseudo steady state analysis

If the transients are relegated solely to the formation and included as a flux captured at a snapshot in time, then the overall heat transfer coefficient can be calculated from Eq. (32) without any need for explicitly formulating an energy balance for each individual wellbore layer. If at least one layer is a fluid layer, then the interface temperatures are required to model the natural convection which then renders the procedure iterative. At each step of the iteration, the interface temperatures can be calculated as follows - given a value of the fluid temperature evaluated a certain iteration step, the temperature at the conduit wall is estimated from the forced convection component of the overall thermal resistance as

$${\overline{T}}_0=T_0-\frac{\overline{U}}{\mathit{hR}}\left[T\left(z\right)-T_{\mathrm{geo}}\left(z\right)\right]$$E109

Subsequently the temperatures at each of the outer layers is evaluated as

$${\overline{T}}_j={\overline{T}}_{j-1}-ln\left(\frac{{\overline{r}}_j}{{\overline{r}}_{j-1}}\right)\frac{\overline{U}}{k_j}\left[T\left(z\right)-T_{\mathrm{geo}}\left(z\right)\right]\begin{array}{cc},& \end{array}j=1\cdots N$$E110

Once the iteration has converged, the temperature at the layer mid radius is evaluated as

$$T_j=\frac{1}{2}\left({\overline{T}}_{j-1}+{\overline{T}}_j\right)\begin{array}{cc},& \end{array}j=1\cdots N$$E111

Note that in this case, since the fluxes are evaluated using interface temperature differences, there is no need to evaluate interface thermal conductivities as in the transient analysis case.

8. Conclusions

The key aspects of wellbore heat transfer cover the entire gamut of thermal energy transport mechanisms from advection/convection in wellbores, to conduction across wellbore tubular and cement layers, to natural convection in trapped annuli, and diffusion in semi-infinite domains from a wellbore to the surrounding formation layer. A term by term derivation of the transport equation using the enthalpy formulation is crucial to understanding the relative importance of the various energy terms.

The mathematical models developed are applicable to a very wide range of wellbore operations, from production and injection to circulation. While thermal transients can generally be ignored for long term production scenarios, significant errors can result from ignoring them in shorter injection and circulating scenarios. When the flowrates exceed a certain threshold, the seemingly counter-intuitive temperature profiles can be explained in terms of the Negative-Joule-Thomson effect.

The most efficient approach to solving transient wellbore heat transfer problems is by Laplace transformation of the governing equations. An efficient method of inversion back to the physical domain is with the use of the powerful Gaver-Stehfest function sampling algorithm.

Appendix

Fourier-Bessel coefficients

The solution to the one dimensional radial transient diffusion problem defined by Eq. (61) with a constant boundary condition ($\varphi \left(1,\xi ,\tau \right)$ = 1 in Eq. (63)) at the inner boundary and an insulated outer boundary at the radial location ${\eta }_{\infty }$ is

$$\varphi \left(\xi ,\tau \right)=1-\sum _{j=1}^{\infty }{C}_j{F}_j\left(\eta \right)e^{-{\lambda }_j^2\tau }$$E112

in which the radial eigenfunction is

$$F_k\left(\eta \right)=J_1\left({\mu }_k{\eta }_{\infty }\right)Y_0\left({\mu }_k\eta \right)-Y_1\left({\mu }_k{\eta }_{\infty }\right)J_0\left({\mu }_k\eta \right)$$E113

where $J_k\left(\right)$ and $Y_k\left(\right)$ for k = 0, 1 are the Bessel Functions of the First and Second Kind of Order k and the eigenvalues are the zeroes of

$$J_1\left({\mathit{\mu \eta }}_{\infty }\right)Y_0\left(\mu \right)-Y_1\left({\mathit{\mu \eta }}_{\infty }\right)J_0\left(\mu \right)=0$$E114

The Fourier-Bessel coefficients are

$$C_k={\left[\underset{1}{\overset{D}{\int }}{\left[F_k\left(\eta \right)\right]}^2\eta \mathrm{d}\eta \right]}^{-1}\underset{1}{\overset{D}{\int }}{F}_k\left(\eta \right)\eta \mathrm{d}\eta$$E115

The flux at the inner boundary is the entity of interest, so that the coefficient ${\overline{C}}_j$ in Eq. (65) is given by

$${\overline{C}}_j=C_j{\left.\frac{\partial F_j}{\partial \eta }\right|}_{\eta =1}=C_j\left[Y_1\left({\mu }_k{\eta }_{\infty }\right)J_1\left({\mu }_k\right)-J_1\left({\mu }_k{\eta }_{\infty }\right)Y_1\left({\mu }_k\right)\right]$$E116

Stehfest accelerators

The coefficients in the Gaver Summation known as Stehfest accelerators are defined as

$${\lambda }_i={\left(-1\right)}^{\mathrm{Int}\left(\frac{N}{2}+i\right)}\sum _{k=\mathrm{Int}\left(\frac{i+1}{2}\right)}^{k=Min\left(i,\frac{N_{\mathrm{GS}}}{2}\right)}\frac{k^{\frac{N_{\mathrm{GS}}}{2}}\left(2k\right)!}{\left(\frac{n}{2}-k\right)!k!\left(k-1\right)!\left(i-k\right)\left(2k-i\right)!},k=1\cdots N_{\mathrm{GS}}$$E117

and are tabulated below for the first few even orders of the method.

k	Stehfest accelerators
	2	4	6	8	10
1	2	−2	1	−3.333E−01	8.333E−02
2	−2	26	−49	4.833E+01	−3.208E+01
3		−48	366	−9.060E+02	1.279E+03
4		24	−858	5.465E+03	−1.562E+04
5			810	−1.438E+04	8.424E+04
6			−270	1.873E+04	−2.370E+05
7				−1.195E+04	3.759E+05
8				2.987E+03	−3.401E+05
9					1.641E+05
10					−3.281E+04

Table 1.

Stehfest accelerators.

Nomenclature

$A$	Flow Area (m^2)
$c_{\mathrm{JT}}$	Joule-Thomson Coefficient (K-m^2/N)
$c_{\mathrm{p}}$	Specific Heat at Constant Pressure (J/kg-K)
$\overline{D}$	Hydraulic Diameter (m)
$e$	Total Specific Energy (J/kg-K)
$f$	Friction Factor
$g$	Acceleration due to Gravity (m^2/s)
$h$	Specific Enthalpy (J/kg) / heat transfer coefficient (W/m^2-K)
$k$	Thermal Conductivity (W/m-K)
$L$	Conduit Length (m)
$\dot{m}$	Mass Flow Rate (kg/s)
$n$	Power Law Index (dimensionless)
$\mathrm{Nu}$	Nusselt Number (dimensionless)
$\mathrm{Nu}$	Nusselt Number (dimensionless)
$P$	Pressure (N/m^2)
$Pr$	Prandtl Number
$Pr$	Prandtl Number (dimensionless)
$Re$	Reynolds Number (dimensionless)
$\dot{q}$	Heat Fux per Unit Area (W/m^2)
$r$	local radius (m)
$\overline{r}$	Interface Radius (m)
$R$	Conduit Radius (m)
$\overline{R}$	Wellbore Outer Radius (m)
$Re$	Reynolds Number (dimensionless)
$Re$	Reynolds Number (dimensionless)
$s$	Laplace variable in Frequency Domain
$t$	time (s)
$T$	Temperature (K)
$u$	Specific Internal Energy (J/kg)
$U$	Overall Heat Transfer Coefficient (W/m^2-K)
$\overline{U}$	Overall Heat Transfer Conductance (W/m-K)
$V$	Flow Velocity (m/s)
$y$	True Vertical Depth (m)
$z$	Streamwise Coordinate (m)
$\alpha$	Isobaric Volume Expansivity(1/K)
$\beta$	Isothermal Compressibility (m^2/N)
$\kappa$	Annulus Radius Ratio ()
$\delta$	Directional Index
$\varphi$	Dimensionless Formation Temperature
$\eta$	Dimensionless radial Coordinate
$\nu$	Kinematic Viscosity (m^2/s)
$\theta$	Dimensionless Fluid Temperature/ Wellbore Local Inclination Inclination
$\rho$	Fluid Density (kg/m^3)
$\tau$	Dimensionless Time (Fourier Number)
${\tau }_{\mathrm{w}}$	Wall Spear Stress (Pa)
$\xi$	Dimensionless Streamwise coordinate
a	Annulus
Form	Formation
Geo	Geothermal
i	Inner Surface
o	Outer Surface
p	Pipe
wb	Wellbore