A Transient Analytical Model for Predicting Wellbore/Reservoir Temperature and Stresses during Drilling with Fluid Circulation

Abstract

Accurate characterization of heat transfer in a wellbore during drilling, which includes fluid circulation, is important for wellbore stability analysis. In this work, a pseudo-3D model is developed to simultaneously calculate the heat exchange between the flowing fluid and the surrounding media (drill pipe and rock formation) and the in-plane thermoelastic stresses. The cold drilling fluid descends through the drill pipe at constant injection rates and returns to the ground surface via the annulus. The fluid circulation will decrease the wellbore bottom temperature and reduce the near-wellbore high compressive stress, potentially leading to tensile fracturing of the well. The governing equations for the coupled heat transfer stress problem are formulated to ensure that the most important parameters are taken into account. The wellbore is subject to a non-hydrostatic in situ far-field stress field. In modeling heat exchange between fluid and surrounding media, the heat transfer coefficients are dependent on fluid properties and flow behavior. Analytical solutions in the Laplace space are obtained for the temperatures of the fluid in both the drill pipe and annulus and for the temperature and stress changes in the formation. The numerical results in the time domain are obtained by using an efficient inversion approach. In particular, the near-well stresses are compared for the cases with fixed and time-dependent cooling wellbore conditions. This comparison indicates that the using a fixed temperature wellbore conditions may over-estimate or under-estimate the bottom-hole stress change, potentially leading to wellbore stability problems.

1. Introduction

Drilling and completion of wells provides a direct way to extract the resources such as oil, gas and geothermal energy from beneath the earth’s surface. Assessment of the wellbore stability issues is an important but very challenging topic related to safe and efficient drilling operations. Instability of wellbores historically costs the oil and gas industry billions of dollars worldwide annually [1]. During drilling, fluid is circulated through the drill pipe/annulus system and this fluid circulation process serves several fundamental functions: (1) removal of the large amount of drilled cuttings from the wellbore [2]; (2) prevention of breakdown (tensile failure) and breakout (shear failure) of the wellbore by adjusting the mud weight and circulation rate, thus controlling the pressure [3,4]; and (3) cooling and lubricating of the drill bit [5]. In addition, the heat exchange between drilling fluids and surrounding formation results in thermoelastic stress redistribution near the wellbore.

Wellbore stability is affected by a variety of factors, such as mechanical impact as a result of drillstring vibration during drilling [6], formation material properties [7,8] and change in loading boundary conditions. Especially, with the exploration and development of unconventional energy resources such as shale gas/oil and geothermal energy, wells are being drilled deeper and at orientations with respect to the stress field that result in higher stress concentrations around the well. In such cases, the temperature and pressure conditions also become more extreme. For example, at a depth of 5 km in geothermal reservoirs, the rock temperature can be more than 270 °C [9]. Under high temperature and high pressure (HTHP) conditions, with cold fluid injected into the wellbore, the heat transfer behavior in both the fluid and rock formation must be accounted for in estimating the complex stress redistribution. Therefore, understanding the transient thermal behavior in the fluid and accounting for the interaction between the wellbore and reservoir rock must be included in the wellbore stability analysis.

A fully coupled well/reservoir (W/R) drilling system analysis generally contains three parts, i.e., fluid circulation (heat and mass) in the wellbore, energy exchange between the wellbore and the reservoir, and flow and deformation in rock formations. A comprehensive model for predicting the wellbore integrity should address these three parts at the same time. Most previous studies have focused on one or two of these parts. As for the first two processes, many models and approaches have been developed to describe heat diffusion, fluid flow and exchange behaviors in the system.

During the past several decades, many models and approaches have been developed to describe the heat exchange process and stress-associated wellbore stability in the oil, gas and geothermal industry. For example, Bullard [10] and Moss and White [11] used a line/source model to study the wellbore/reservoir (W/R) temperature. Edwardson [12] and Tragesser [13] obtained the formation temperature during mud circulation with an “exact method”, which solved the differential equation of heat conduction. Approximate solutions were presented by Ramey [14] for the case where steady-state heat transfer in the wellbore is coupled with transient heat flow in the formation. Later, Raymond [15] derived the complete equations for the heat transfer in the W/R system, and since then a great deal of further work has been carried out. For example, Holmes and Swift [16], Keller [17], Sump [18], Arnold [19], Kabir et al. [20], Fomin et al. [21] and Wu et al. [22] studied similar problems using different approaches or extending Raymond’s model to a more general case.

In addition, many models have also been developed to describe the rock deformation and wellbore stability. According to the effects to be considered, which are mechanical (M), hydro- (H), and thermal (T), these models are mainly divided into three types, i.e., classical elasticity (CM), poro-elasticity (HM) and thermo-poro-elasticity (THM). The well-known Kirsch solutions [23] for 2D cases, and those obtained by Hiramatsu and Oka [24] for 3D cases, belong to the first type. The second type of the wellbore models applied the HM theory developed by Biot [25]. For example, Carter and Booker [26] analyzed the consolidation of a linear elastic soil when a deep circular tunnel is cut. Detournay and Cheng [27] obtained the poroelastic (HM) response of a wellbore in a non-hydrostatic stress field. Rajapakse [28], Ekbote et al. [29] and Chen and Yu [30] studied similar two-dimensional cases in an isotropic HM media. In addition, Abousleiman and Cui [31] extended the loading decomposition scheme (LDS) and obtained analytical solutions for the cases where inclined wellbores are in transversely isotropic HM medium.

The third type takes into account the thermal effect on stress and deformation for analysis of wellbore stability and hydraulic fracturing, and is especially applicable to deep hot reservoirs. For example, Mctigue [32], Kurashige [33] and Charlez [34] derived similar theoretical formulations that can be used for studying coupled thermo-poro-elasticity wellbore problems. By using the above theories, Zhou et al. [35], Ghassemi and Diek [36], Wang and Dusseault [37], Choi et al. [38] and Wu et al. [39,40] studied the THM behaviors around a wellbore or sphere. Li et al. [41] and Abousleiman and Ekbote [42] and Gao et al. [43] applied the LDS to an inclined wellbore in an isotropic and transversely isotropic THM medium.

The motivation of the present work arises for two reasons. First, most of the previous fluid circulation models, such as [10,11,12,13,14,15,16,17,18,19,20,21,22], consider only the heat transfer between the wellbore and the reservoir, without considering any mechanical responses of the system. In addition, most existing analytical solutions for the wellbore circulation system are based on the assumption that the heat transfer in the wellbore or rock formation is steady state [8,9,10,11,12,13,14]. Use of complete transient solutions requires resorting to fully numerical approaches [15].

Second, in nearly all of the existing wellbore stability models, such as [23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43], the boundary conditions at wellbore wall are assumed to be fixed (especially with a fixed wellbore temperature) or given in an arbitrary way. Obviously, this does not reflect the real wellbore conditions where the temperature changes with fluid circulation or drilling time. This simplification can lead to under- or over-estimation on the near-well stress field. Circulation of a fluid with a temperature lower than the formation temperature leads to a gradual reduction in both the bottom-hole temperature and the near-wellbore compressive stresses. The temperature at most locations along the wellbore also reduces with time. As different temperatures at the wellbore wall generate different heat transfer efficiencies in the formation, the time-dependent response of the wellbore wall temperature should be considered.

The purpose of the present work is both to obtain semi-analytical solutions for the temperature evolution in the injection string, annulus and rock formation during fluid circulation and to provide an accurate evaluation for near-well stress changes under the variable wellbore temperature, which is a more reasonable auxiliary condition for wellbore stability analysis during fluid circulation.

2. Problem Formulation

2.1. Problem Description

The geometry configuration for the present model of a W/R system is shown in Figure 1. The origin of the cylindrical coordinate system is set to coincide with the centre of the well at the surface. The system is initially in an equilibrium state with the formation temperature distribution T₀^* = A₀ z + B₀ [21], where A₀ is the geothermal gradient and B₀ is the surface soil temperature. At an arbitrary depth of z, the cross-sectional wellbore is subjected to maximum horizontal principal stress σ_H(z), minimum horizontal principal stress σ_h(z), vertical principal stress σ_v(z) and wellbore pressure p_w(z). The temperatures of the fluid inside the annulus and tubing are denoted as T_a^* and T_d^*, respectively, independent of the radius. The rock temperature is denoted as T_r^*. The temperature on the borehole wall is denoted as T_w^*. At time t > 0, a fluid with a given temperature T_in^* is injected into the tubing at a mass rate Q_d and returns to the surface along the annulus when it reaches the bottom of the tubing string. The return mass flow rate is denoted as Q_a.

The wellbore radius is denoted by r_w and the wellbore depth is H. The inner and outer radii of the tubing are r₀ and r_d, respectively. Therefore, the tubing wall thickness is denoted as δ₀ = r_d − r₀ and the effective width of the annulus is w_a = r_w − r_d.

For simplicity, the following assumptions are made: (1) the fluid is incompressible and the rock formation is impermeable, homogeneous, isotropic, linear and elastic; (2) the material properties are constants, independent of temperature change; and (3) the heat conduction along the vertical direction in the formation is neglected.

2.2. Governing Equations

2.2.1. Wellbore Heat Transfer

The fluid temperature in the tubing is determined by the heat exchange between the tubing and the annulus and the heat advection down the tubing. The fluid temperature in the annulus is determined by the heat exchange between the annular fluid and the surrounding rock, the heat exchange between the tubing and the annulus and the heat advection up along the annulus. According to Raymond [15], the advective heat transfer equations for fluids inside the tubing string and in the annulus are expressed as follows.

$$\begin{array}{l}{\rho }_l{c}_l{A}_d{v}_d\frac{\partial T_d^{*}}{\partial z}+2\pi r_d{h}_{ad}(T_d^{*}-T_a^{*})=-{\rho }_l{A}_d{c}_l\frac{\partial T_d^{*}}{\partial t},\\ {\rho }_l{c}_l{A}_a{v}_a\frac{\partial T_a^{*}}{\partial z}+2\pi r_d{h}_{ad}(T_d^{*}-T_a^{*})+2\pi r_w{h}_w(T_w^{*}-T_a^{*})={\rho }_l{A}_a{c}_l\frac{\partial T_a^{*}}{\partial t},\end{array}$$

(1)

where c_l is the fluid specific heat capacity; A_m and T_m^* (m = d refers to tubing and m = a refers to annulus) denote the cross sectional area and the fluid temperature, respectively; h_ad is the overall heat transfer coefficient (HTC) characterizing heat exchange between the fluids in the tubing and in the annulus through the tubing wall; and h_w is the film HTC between the fluid in the annulus and the surrounding rock formation. The velocities v_m (m = d or a) in the tubing and the annuls are expressed as

$$v_a=\frac{Q_a}{A_a{\rho }_l},v_d=\frac{Q_d}{A_d{\rho }_l},A_a=\pi \left(r_w^2-r_d^2\right),A_d=\pi r_0^2.$$

2.2.2. Heat Conduction in Rock

As for the heat transfer in the rock formation, two simplifications can be made. First, as the temperature gradient in the radial direction is much greater than that in the vertical direction, the derivative of the formation temperature with respect to z can be ignored [15]. Then, the governing equation to describe the heat flow in the formation is written as

$$k_r{\nabla }^2{T}_r^{*}={\rho }_r{c}_r\frac{\partial T_r^{*}}{\partial t}+\gamma T^R\frac{\partial {\epsilon }_v}{\partial t}$$

(2)

where ε_v is the volumetric strain. Second, according to the analysis by Coussy [44] and Wu et al. [39], the deformation induced temperature change is very small so that it can be ignored. Therefore, the second term on the right side of Equation (2) is omitted in the present work.

2.2.3. Rock Deformation

The following governing equations in a polar coordinate system (r, θ) are used to describe the mechanical responses of the rock formation, whose temperature distribution is disturbed by the circulating fluid.

• Equations of equilibrium

  $$\frac{\partial {\sigma }_{rr}}{\partial r}+\frac{1}{r}\frac{\partial {\sigma }_{r\theta }}{\partial \theta }+\frac{{\sigma }_{rr}-{\sigma }_{\theta \theta }}{r}=0,\frac{\partial {\sigma }_{r\theta }}{\partial r}+\frac{1}{r}\frac{\partial {\sigma }_{\theta \theta }}{\partial \theta }+2\frac{{\sigma }_{r\theta }}{r}=0$$

  (3)
• Strain-displacement relations

  $${\epsilon }_{rr}=\frac{\partial u_r}{\partial r},{\epsilon }_{\theta \theta }=\frac{u_r}{r}+\frac{1}{r}\frac{\partial u_{\theta }}{\partial \theta },{\epsilon }_{r\theta }=\frac{1}{2}\left[\frac{1}{r}\frac{\partial u_r}{\partial \theta }+\frac{\partial u_{\theta }}{\partial r}-\frac{u_{\theta }}{r}\right]$$

  (4)
• Linear TM constitutive equation for isotropic medium

  $${\sigma }_{ij}=2G\left[\frac{v}{1-2v}{\delta }_{ij}{\epsilon }_v+{\epsilon }_{ij}\right]-{\delta }_{ij}\gamma T_r$$

  (5)

In the above equations, σ_ij is the stress change tensor (positive is tension) and is equal to σ_ij^* − σ_ij^R, the difference between the total stress σ_ij^* and the initial stress σ_ij^R; ε_v is the volumetric strain and ε_ij the strain tensor; T_r = T_r^* − T₀^* denotes the rock temperature change; G, v and K = 2G(1 + v)/3(1 − 2v) are shear modulus, Poisson’s ratio and bulk modulus for the rock formation, respectively; α is the coefficient of volumetric thermal expansion of the rock formation; γ = αK; and δ_ij is Kronecker’s delta.

2.3. Boundary Conditions

The fluid temperature at the inlet of the tubing is specified during drilling. In order to satisfy the continuity condition, the temperature of fluids in the tubing and in the annulus at the bottomhole should be identical. Thus, we have

$$\begin{array}{ll}{T}_d^{*}(z,t)=T_{\mathrm{in}}^{*},& \mathrm{at}z=0,\\ T_d^{*}(z,t)=T_a^{*}(z,t),& \mathrm{at}z=H.\end{array}$$

(6)

In addition, at the wellbore wall, the heat flux from the formation is equal to the heat flux into the wellbore, i.e., the following continuity condition must be met

$$2\pi r_w{h}_w(T_r^{*}-T_a^{*})=2\pi r_w{k}_r\frac{\partial T_r^{*}}{\partial r},\mathrm{at}r=r_w.$$

(7)

The pressure at the wellbore wall at a depth z is assumed to be time-independent and is approximated by

$$p_w(z)={\rho }_{l}gz,$$

(8)

where g is the gravitational acceleration.

The initial condition of the system is given as

$$T_d(z,t)=T_a(z,t)=T_r(z,r,t)=0,\mathrm{on}t=0.$$

(9)

where T_m = T_m^* − T^R (m = d refers to tubing, m = a refers to annulus and m = r refers to formation).

The remote boundary conditions are expressed in terms of changes of temperature and stresses

$$\begin{array}{ll}{T}_r=0& \mathrm{on}r\to \infty ,\\ {\sigma }_{rr}=0,& \mathrm{on}r\to \infty ,\\ {\sigma }_{r\theta }=0,& \mathrm{on}r\to \infty ,\end{array}$$

(10)

while at the wellbore wall, the stress boundary conditions are

$$\begin{array}{ll}{\sigma }_{rr}=-p_w-(p_0+s_0\cos 2\theta )& \mathrm{on}r=r_w,\\ {\sigma }_{r\theta }=S_0\sin 2\theta & \mathrm{on}r=r_w,\end{array}$$

(11)

where

$$p_0=\frac{{\sigma }_{\mathrm{H}\max }+{\sigma }_{\mathrm{H}\min }}{2},s_0=\frac{{\sigma }_{\mathrm{H}\max }-{\sigma }_{\mathrm{H}\min }}{2}.$$

2.4. Heat Transfer Coefficients (HTCs)

In the present paper, the value of the HTC for forced advection along a pipe is calculated by

$$\mathrm{HTC}=\frac{N_u{k}_l}{D},$$

(12)

where k_l is the fluid thermal conductivity, D is the hydraulic diameter (a characteristic length) and N_u is the dimensionless Nusselt number. If the Nusselt number is known, the HTC value can be obtained.

When Re ≤ 2320, the flow is laminar and N_u is equal to be 3.66. When Re ≥ 10,000, the flow is turbulent and the Nusselt number is calculated with the Mikheev equation [45] which is reported to be the reliable formula for estimating the HTC for the flow of Newtonian fluids in a tube. In particular, the Nusselt numbers for turbulent flow at the inner tubing and wellbore walls, respectively, are given by Isachenko and Osipova [46] as follows

$$\begin{array}{l}{N}_d=\frac{h_d(2r_0)}{k_l}=0.021{({\mathrm{Re}}_d)}^{0.8}{({\mathrm{Pr}}_l)}^{0.43}{({\mathrm{Pr}}_l/{\mathrm{Pr}}_d)}^{0.25},\\ N_a=\frac{h_a(2w_a)}{k_l}=0.021{({\mathrm{Re}}_a)}^{0.8}{({\mathrm{Pr}}_l)}^{0.43}{({\mathrm{Pr}}_l/{\mathrm{Pr}}_r)}^{0.25},\end{array}$$

(13)

where μ is the fluid dynamic viscosity, and k_j, ρ_j and c_j are the thermal conductivity, mass density and specific heat capacity (here the symbol j = l, d or r denotes the liquid, tubing or rock), respectively. The Reynolds number for fluid flow in the tubing and annulus and the Prandtl number are defined below

$${\mathrm{Re}}_d=\frac{{\rho }_l{v}_{d}2r_0}{\mu },{\mathrm{Re}}_a=\frac{{\rho }_l{v}_{a}2w_a}{\mu },{\mathrm{Pr}}_d=\frac{\mu c_d}{k_d},{\mathrm{Pr}}_l=\frac{\mu c_l}{k_l},{\mathrm{Pr}}_r=\frac{\mu c_r}{k_r}.$$

When 2320 < Re < 10,000 (for transition flow), a linear interpolation method is used, as show in Figure 2.

The overall HTC across the tubing wall is determined by the empirical formula derived by Willhite [47], which gives

$$\frac{1}{h_{ad}}=\frac{1}{h_d}+\frac{1}{k_d/{\delta }_0}+\frac{1}{h_a},$$

(14)

3. Dimensional Analysis

By using the following transformation

$$\begin{array}{l}{\mathsf{\Theta }}_a=\frac{T_a^{*}-T^R}{A_{0}H},{\mathsf{\Theta }}_d=\frac{T_d^{*}-T^R}{A_{0}H},{\mathsf{\Theta }}_r=\frac{T_r^{*}-T^R}{A_{0}H},{\mathsf{\Theta }}_{\mathrm{in}}=\frac{T_{\mathrm{in}}^{*}-B_0}{A_{0}H},\\ {\mathsf{\Sigma }}_{ij}=\frac{{\sigma }_{ij}}{G},{\mathsf{\Omega }}_i=\frac{u_i}{r_w},R=\frac{r}{r_w},\tau =\frac{d_{r}t}{r_w^2},Z=\frac{z}{H},B_i=\frac{h_a{r}_w}{k_r},\\ d_r=\frac{k_r}{{\rho }_r{c}_r},\beta =\sqrt{\frac{2(1-v)}{1-2v}},c=\frac{\gamma A_{0}H}{G},a_1=\frac{\gamma }{{\rho }_r{c}_r},a_2=\frac{B_0}{A_{0}H}\end{array}$$

(15)

the governing Equation (1) for heat transfer in the wellbore and heat exchange with reservoir is rewritten as

$$\begin{array}{l}{\chi }_d\frac{\partial {\mathsf{\Theta }}_d}{\partial \tau }+\frac{\partial {\mathsf{\Theta }}_d}{\partial z}=b[{\mathsf{\Theta }}_a-{\mathsf{\Theta }}_d]-1,\\ {\chi }_a\frac{\partial {\mathsf{\Theta }}_a}{\partial \tau }-\frac{\partial {\mathsf{\Theta }}_a}{\partial z}=-d[{\mathsf{\Theta }}_a-{\mathsf{\Theta }}_d]-e[{\mathsf{\Theta }}_a-{\mathsf{\Theta }}_w]+1,\end{array}$$

(16)

where the coefficients are defined as

$$\begin{array}{l}{\chi }_d=\frac{H{d}_r}{r_w^2{v}_d},{\chi }_a=\frac{H{d}_r}{r_w^2{v}_a},b=\frac{2N_{ad}{\zeta }_H}{\mathrm{P}{e}_d{\zeta }_0},d=\frac{2N_{ad}{\zeta }_H}{\mathrm{P}{e}_a(1+{\zeta }_d)},e=\frac{2N_w{\zeta }_H}{\mathrm{P}{e}_a(1-{\zeta }_d^2)},\\ N_{ad}=\frac{h_{ad}2r_d}{k_l},{\zeta }_H=\frac{H}{r_w},{\zeta }_0=\frac{r_0}{r_w},{\zeta }_d=\frac{r_d}{r_w},\mathsf{\Phi }=\frac{Q_a}{Q_d},{\epsilon }_d=\frac{k_d}{k_l},{\epsilon }_r=\frac{k_r}{k_l}.\end{array}$$

(17)

The stress–strain correlation (5) is given as

$${\mathsf{\Sigma }}_{ij}=\left({\beta }^2-2\right){\delta }_{ij}{\mathsf{\Xi }}_v+2{\mathsf{\Xi }}_{ij}-{\delta }_{ij}c{\mathsf{\Theta }}_r$$

(18)

and the strain-displacement relationship is rewritten as

$${\mathsf{\Xi }}_{rr}=\frac{\partial {\mathsf{\Omega }}_R}{\partial R},{\mathsf{\Xi }}_{\theta \theta }=\frac{{\mathsf{\Omega }}_R}{R}+\frac{1}{R}\frac{\partial {\mathsf{\Omega }}_{\theta }}{\partial \theta },{\mathsf{\Xi }}_r{}_{\theta }=\frac{1}{2}\left[\frac{1}{R}\frac{\partial {\mathsf{\Omega }}_R}{\partial \theta }+\frac{\partial {\mathsf{\Omega }}_{\theta }}{\partial R}-\frac{{\mathsf{\Omega }}_{\theta }}{R}\right].$$

(19)

In addition, heat conduction in the rock formation becomes

$${\nabla }^2{\mathsf{\Theta }}_r=\frac{\partial {\mathsf{\Theta }}_r}{\partial \tau }$$

The boundary and initial conditions for the wellbore become

$$\begin{array}{ll}\frac{\partial {\mathsf{\Theta }}_r}{\partial R}={\mathrm{B}}_{\mathrm{i}}[{\mathsf{\Theta }}_r-{\mathsf{\Theta }}_a],& \mathrm{on}R=1,\\ {\mathsf{\Theta }}_d={\mathsf{\Theta }}_{\mathrm{in}},& \mathrm{on}Z=0,\\ {\mathsf{\Theta }}_d={\mathsf{\Theta }}_a,& \mathrm{on}Z=1,\\ {\mathsf{\Theta }}_r={\mathsf{\Theta }}_d={\mathsf{\Theta }}_a=0,& \mathrm{on}\tau =0.\end{array}$$

(20)

and the boundary conditions for the reservoir are rewritten as

$$\begin{array}{ll}{\mathsf{\Theta }}_r=0& \mathrm{on}R\to \infty ,\\ {\mathsf{\Sigma }}_{RR}=0,& \mathrm{on}R\to \infty ,\\ {\mathsf{\Sigma }}_{R\theta }=0,& \mathrm{on}R\to \infty ,\end{array}$$

(21)

While, at the wellbore wall, the stress boundary conditions are

$$\begin{array}{ll}{\mathsf{\Sigma }}_{RR}=-{\mathsf{\Pi }}_w-(P_0+S_0\cos 2\theta )& \mathrm{on}r=r_w,\\ {\mathsf{\Sigma }}_{R\theta }=S_0\sin 2\theta & \mathrm{on}r=r_w,\end{array}$$

(22)

where

$${\mathsf{\Pi }}_w=\frac{p_w}{G},P_0=\frac{p_0}{G},S_0=\frac{s_0}{G}.$$

4. Solution Method

From Equations (1)–(11) we find that the heat transfer and deformation process represent a linear system. By using the same approach as that in [27,39], the non-hydrostatic remote stresses are decomposed into isotropic and deviatoric components, i.e., p₀ and s₀. According to the superposition principle, the present problem can be solved by studying two loading modes, i.e., isotropic (Mode 1) and anti-symmetric (Mode 2). In Mode 1, the fluid circulates in the wellbore that is subjected to the isotropic far-field stress p₀, wellbore pressure and initial temperature conditions, while in Mode 2 the wellbore is subjected to the anti-symmetric far-field stress s₀ and homogenous boundary conditions. The final solutions of near-well stress fields are the combination of the above results. It should be mentioned that the wellbore temperature varies with circulation time. In addition, the temperature evolution in the pipe and annulus can be obtained analytically in Laplace space.

The Laplace transformation is defined as

$$\widehat{F}(s)={\displaystyle {\int }_{-\infty }^{+\infty }{e}^{-s\tau }f(\tau )d\tau }$$

(23)

where s is a complex number and the symbol ^ denotes the Laplace transform of a quantity (${\widehat{\mathsf{\Theta }}}_a$, ${\widehat{\mathsf{\Theta }}}_d$, ${\widehat{\mathsf{\Theta }}}_r$, ${\widehat{\mathsf{\Sigma }}}_{RR}$, ${\widehat{\mathsf{\Sigma }}}_{\theta \theta }$ and ${\widehat{\mathsf{\Sigma }}}_{R\theta }$ for example).

The stress boundary condition for the Mode 1 and Mode 2 problems are

$${\mathsf{\Sigma }}_{RR}=-{\mathsf{\Pi }}_w-P_0\mathrm{on}R=1,{\mathsf{\Sigma }}_{RR}=0\mathrm{on}R\to \infty ,$$

and

$$\begin{array}{llll}{\mathsf{\Sigma }}_{RR}=-S_0\cos 2\theta & \mathrm{on}R=1,& {\mathsf{\Sigma }}_{RR}=0& \mathrm{on}R\to \infty ,\\ {\mathsf{\Sigma }}_{R\theta }=S_0\sin 2\theta & \mathrm{on}R=1,& {\mathsf{\Sigma }}_{R\theta }=0& \mathrm{on}R\to \infty ,\end{array}$$

(24)

respectively.

In Mode 1, it should be noted that ${\mathsf{\Omega }}_{\theta }$ = 0 and ${\mathsf{\Sigma }}_{R\theta }$ = 0. Substitution of Equation (16) into the dimensionless forms of Equations (2) and (3) leads to the volumetric strain

$${\delta }_v=\frac{c}{{\beta }^2}{\mathsf{\Theta }}_r+C_0$$

(25)

where the constant C₀ is equal to zero as δ_v vanishes at infinity. After putting δ_v into Equation (19), we obtain

$${\nabla }^2{\mathsf{\Theta }}_r=\frac{\partial {\mathsf{\Theta }}_r}{\partial \tau }$$

By applying the Laplace transformation to the above equation, the following ordinary differential equation (ODE) is obtained for the rock temperature change

$$s{\widehat{\mathsf{\Theta }}}_r=\frac{1}{R}\frac{d{\widehat{\mathsf{\Theta }}}_r}{dR}+\frac{d^2{\widehat{\mathsf{\Theta }}}_r}{d{R}^2}$$

(26)

whose solution in Laplace space is obtained as

$${\widehat{\mathsf{\Theta }}}_r=F_1(Z,s)K_0(\sqrt{s}R)+F_2(Z,s)I_0(\sqrt{s}R).$$

where I_n and K_n denote the modified Bessel functions of the first and second kind of order n, respectively; and F₁(s) and F₂(s) are unknowns to be determined by the boundary conditions. Considering zero rock temperature change at infinity, the function F₂(Z, s) = 0 and then

$${\widehat{\mathsf{\Theta }}}_r=F_1(Z,s)K_0(\sqrt{s}R),{\widehat{\mathsf{\Theta }}}_w=F_1(Z,s)K_0(\sqrt{s})$$

(27)

By using the first boundary condition in Equation (20), the annulus temperature change is found

$${\widehat{\mathsf{\Theta }}}_a=\mathsf{\Lambda }{F}_1(Z,s),\mathrm{where}\mathsf{\Lambda }=K_1(\sqrt{s})\sqrt{s}/{\mathrm{B}}_{\mathrm{i}}+K_0(\sqrt{s})$$

(28)

After applying the Laplace transformation to Equation (16), we obtain the following ODEs

$$\begin{array}{l}f{\widehat{\mathsf{\Theta }}}_d+\frac{d{\widehat{\mathsf{\Theta }}}_d}{dZ}=b\mathsf{\Lambda }{F}_1-\frac{1}{s},\\ h{F}_1-\mathsf{\Lambda }\frac{d{F}_1}{dZ}=d{\widehat{\mathsf{\Theta }}}_d+\frac{1}{s},\end{array}$$

(29)

whose solutions are found to be

$$\begin{array}{l}{F}_1(Z,s)=Q_1(s)e^{-{\lambda }_{1}Z}+Q_2(s)e^{-{\lambda }_{2}Z}+C_1,\\ {\widehat{\mathsf{\Theta }}}_d(Z,s)=C_2{Q}_1(s)e^{-{\lambda }_{1}Z}+C_3{Q}_2(s)e^{-{\lambda }_{2}Z}+C_4,\end{array}$$

(30)

and then we have

$$\begin{array}{l}{\widehat{\mathsf{\Theta }}}_a(Z,s)=\left[Q_1(s)e^{-{\lambda }_{1}Z}+Q_2(s)e^{-{\lambda }_{2}Z}+C_1\right]\mathsf{\Lambda },\\ {\widehat{\mathsf{\Theta }}}_r(R,Z,s)=\left[Q_1(s)e^{-{\lambda }_{1}Z}+Q_2(s)e^{-{\lambda }_{2}Z}+C_1\right]K_0(\sqrt{s}R),\end{array}$$

(31)

where the coefficients are defined as

$$\begin{array}{l}{C}_1=\frac{f-d}{s(fh-d\mathsf{\Lambda }b)},C_2=\frac{(\mathsf{\Lambda }{\lambda }_1+h)}{d},C_3=\frac{(\mathsf{\Lambda }{\lambda }_2+h)}{d},C_4=\frac{\mathsf{\Lambda }b-h}{s(fh-d\mathsf{\Lambda }b)},\\ {\lambda }_1=\frac{-h+f\mathsf{\Lambda }+\sqrt{{(f\mathsf{\Lambda }+h)}^2-4db{\mathsf{\Lambda }}^2}}{2\mathsf{\Lambda }},{\lambda }_2=\frac{-h+f\mathsf{\Lambda }-\sqrt{{(f\mathsf{\Lambda }+h)}^2-4db{\mathsf{\Lambda }}^2}}{2\mathsf{\Lambda }},\\ f=({\chi }_{d}s+b),h=({\chi }_{a}s+d+e)\mathsf{\Lambda }-e{K}_0(\sqrt{s}).\end{array}$$

Based on the boundary conditions, the coefficients Q₁(s) and Q₂(s) are calculated to be

$$\begin{array}{l}{Q}_1(s)=\frac{(s{C}_4-{\mathsf{\Theta }}_{\mathrm{in}})(\mathsf{\Lambda }-C_3)e^{-{\lambda }_2}+s{C}_3(C_4-\mathsf{\Lambda }{C}_1)}{s\left[C_3{e}^{-{\lambda }_1}(\mathsf{\Lambda }-C_2)-C_2{e}^{-{\lambda }_2}(\mathsf{\Lambda }-C_3)\right]},\\ Q_2(s)=-\frac{(-{\mathsf{\Theta }}_{\mathrm{in}}+C_4)(-C_2+\mathsf{\Lambda })e^{-{\lambda }_1}+s{C}_2(C_4-\mathsf{\Lambda }{C}_1)}{s\left[C_3{e}^{-{\lambda }_1}(\mathsf{\Lambda }-C_2)-C_2{e}^{-{\lambda }_2}(\mathsf{\Lambda }-C_3)\right]}.\end{array}$$

From Equation (18), we obtain the induced stresses for Mode 1 if the radial displacement is known

$${\mathsf{\Sigma }}_{RR}=-\frac{2{\mathsf{\Omega }}_R}{R},{\mathsf{\Sigma }}_{\theta \theta }=-\frac{2c}{{\beta }^2}{\mathsf{\Theta }}_r+\frac{2{\mathsf{\Omega }}_R}{R},{\mathsf{\Sigma }}_{zz}=-\frac{2c}{{\beta }^2}{\mathsf{\Theta }}_r$$

(32)

Based on Equations (3), (25) and the second of Equation (31), we obtain the radial displacement

$${\mathsf{\Omega }}_R=-\frac{\mathsf{\Lambda }{K}_1(\sqrt{s}R)}{\sqrt{s}}+\frac{W(s)}{R}$$

(33)

where W(s) is an unknown function and is calculated by the stress boundary condition Equation (22) at the wellbore wall

$$W(s)=\frac{2\mathsf{\Lambda }s{K}_1\left(\sqrt{s}\right)-\left(-{\mathsf{\Pi }}_w-P_0\right)\sqrt{s}}{2s\sqrt{s}}.$$

The solution to the Mode 2 loading provides the temperature, deformation and stress induced by the anti-symmetric far-field stress.

Again, substitution of Equation (16) into the dimensionless forms of Equations (2) and (3) leads to the following equations

$$\begin{array}{l}\frac{\partial {\delta }_v}{\partial R}-\frac{1}{{\beta }^2}\frac{1}{R}\frac{\partial \varpi }{\partial \theta }=\frac{c}{{\beta }^2}\frac{\partial \mathsf{\Theta }}{\partial R},\\ \frac{1}{R}\frac{\partial {\delta }_v}{\partial \theta }+\frac{1}{{\beta }^2}\frac{\partial \varpi }{\partial R}=\frac{c}{{\beta }^2}\frac{1}{R}\frac{\partial \mathsf{\Theta }}{\partial \theta }.\end{array}$$

(34)

Here the rotation displacement $\varpi$ is defined as

$$\varpi =\frac{\partial {\mathsf{\Omega }}_{\theta }}{\partial R}+\frac{{\mathsf{\Omega }}_{\theta }}{R}-\frac{1}{R}\frac{\partial {\mathsf{\Omega }}_R}{\partial \theta }.$$

(35)

The Laplace transforms of the variables ${\delta }_v$, ${\mathsf{\Omega }}_R$, ${\mathsf{\Sigma }}_{RR}$, ${\mathsf{\Sigma }}_{\theta \theta }$, ${\mathsf{\Theta }}_r$, $\varpi$ and ${\mathsf{\Omega }}_{\theta }$ are defined in the same way in Equation (23) in terms of ${\overline{\delta }}_v$, ${\overline{\mathsf{\Omega }}}_R$, ${\overline{\mathsf{\Sigma }}}_{RR}$, ${\overline{\mathsf{\Sigma }}}_{\theta \theta }$, ${\overline{\mathsf{\Theta }}}_r$, $\overline{\varpi }$ and ${\overline{\mathsf{\Omega }}}_{\theta }$. Under Mode 2 loading, these Laplace transformed quantities are expressed as follows

$$\begin{array}{l}({\overline{\delta }}_v,{\overline{\mathsf{\Omega }}}_R,{\overline{\mathsf{\Sigma }}}_{RR},{\overline{\mathsf{\Sigma }}}_{\theta \theta },{\overline{\mathsf{\Theta }}}_r)=({\widehat{\delta }}_v,{\widehat{\mathsf{\Omega }}}_R,{\widehat{\mathsf{\Sigma }}}_{RR},{\widehat{\mathsf{\Sigma }}}_{\theta \theta },{\widehat{\mathsf{\Theta }}}_r)\cos 2\theta ,\\ (\overline{\varpi },{\overline{\mathsf{\Omega }}}_{\theta },{\overline{\mathsf{\Sigma }}}_{R\theta })=(\widehat{\varpi },{\widehat{\mathsf{\Omega }}}_{\theta },{\widehat{\mathsf{\Sigma }}}_{R\theta })\sin 2\theta ,\end{array}$$

(36)

where ${\delta }_v$, ${\mathsf{\Omega }}_R$, ${\mathsf{\Sigma }}_{RR}$, ${\mathsf{\Sigma }}_{\theta \theta }$, ${\mathsf{\Theta }}_r$, $\varpi$ and ${\mathsf{\Omega }}_{\theta }$ are all functions of R and s only.

By Laplace transformation, Equation (34) are changed into the following two ODEs

$$\begin{array}{l}\frac{d{\widehat{\delta }}_v}{dR}-\frac{2}{{\beta }^2}\frac{\widehat{\varpi }}{R}=\frac{c}{{\beta }^2}\frac{d{\widehat{\mathsf{\Theta }}}_r}{dR},\\ \frac{{\widehat{\delta }}_v}{R}-\frac{1}{2{\beta }^2}\frac{d\widehat{\varpi }}{dR}=\frac{c}{{\beta }^2}\frac{1}{R}{\widehat{\mathsf{\Theta }}}_r.\end{array}$$

(37)

The solutions of above ODE are found to be

$${\widehat{\delta }}_v=\frac{c}{{\beta }^2}{\mathsf{\Theta }}_r+\frac{F(s)}{R^2},\widehat{\varpi }=-{\beta }^2\frac{F(s)}{R^2},$$

(38)

where F(s) is an unknown function to be calculated and thus the induced stresses are obtained from Equation (18)

$$\begin{array}{ll}{\widehat{\mathsf{\Sigma }}}_{RR}={\beta }^2\frac{F(s)}{R^2}-2{\widehat{\mathsf{\Xi }}}_{\theta \theta },& {\widehat{\mathsf{\Sigma }}}_{\theta \theta }={\beta }^2\frac{F(s)}{R^2}-2{\widehat{\mathsf{\Xi }}}_{RR},\\ {\widehat{\mathsf{\Sigma }}}_{R\theta }=2{\widehat{\mathsf{\Xi }}}_{R\theta },& {\widehat{\mathsf{\Sigma }}}_{zz}=v({\widehat{\mathsf{\Sigma }}}_{RR}+{\widehat{\mathsf{\Sigma }}}_{\theta \theta }).\end{array}$$

(39)

For Mode 2, the Laplace transforms of the governing Equation (19) is simplified to be

$$\frac{d^2{\widehat{\mathsf{\Theta }}}_r}{\partial R^2}+\frac{1}{R}\frac{d{\widehat{\mathsf{\Theta }}}_r}{\partial R}-\frac{4{\widehat{\mathsf{\Theta }}}_r}{R^2}=s{\widehat{\mathsf{\Theta }}}_r+as\frac{F(s)}{R^2},$$

(40)

whose solution, combined with the far-field boundary conditions, is found to be

$${\widehat{\mathsf{\Theta }}}_r=K_2(\sqrt{s}R)G(s)-\frac{aF(s)}{R^2}.$$

(41)

As ${\widehat{\mathsf{\Theta }}}_r$ = 0 at R = 1, the function G(s) is

$$G(s)=\frac{aF(s)}{K_2(\sqrt{s}R)}.$$

Therefore, the induced temperature and volumetric strain are

$${\widehat{\mathsf{\Theta }}}_r=F(s)a\left[\frac{K_2(\sqrt{s}R)}{K_2(\sqrt{s})}-\frac{1}{R^2}\right],{\widehat{\mathsf{\Xi }}}_v=F(s)m_1\frac{K_2(\sqrt{s}R)}{K_2(\sqrt{s})}+\frac{m_{3}F(s)}{R^2},$$

(42)

where

$$m_1=\frac{ca}{{\beta }^2+ac},m_3=1-m_1.$$

In addition, ${\widehat{\mathsf{\Xi }}}_v$ and $\widehat{\varpi }$ can be expressed through the geometric relationships in Laplace space as follows

$${\widehat{\mathsf{\Xi }}}_v=\frac{\partial {\widehat{\mathsf{\Omega }}}_R}{\partial R}+\frac{{\widehat{\mathsf{\Omega }}}_R}{R}+\frac{2{\widehat{\mathsf{\Omega }}}_{\theta }}{R},\widehat{\varpi }=\frac{\partial {\widehat{\mathsf{\Omega }}}_{\theta }}{\partial R}+\frac{{\widehat{\mathsf{\Omega }}}_{\theta }}{R}+\frac{2{\widehat{\mathsf{\Omega }}}_R}{R}.$$

(43)

By using Equations (38), (42) and (43), the displacements are obtained to be

$$\begin{array}{l}{\widehat{\mathsf{\Omega }}}_R=-\frac{m_{1}F(s)}{\sqrt{sq}{K}_2(\eta )}\left[K_1(\xi )+\frac{2}{\xi }{K}_2(\xi )\right]+\frac{m_{3}F(s)}{2R}+\frac{Q(s)}{R^3},\\ {\widehat{\mathsf{\Omega }}}_{\theta }=-\frac{2m_{1}F(s)}{\sqrt{sq}{K}_2(\eta )}\frac{K_2(\xi )}{\xi }+\frac{m_{4}F(s)}{2R}+\frac{Q(s)}{R^3},\end{array}$$

(44)

where m₄ = −β². Thus, from the constitutive Equation (18), the stresses are obtained

$$\begin{array}{l}{\widehat{\mathsf{\Sigma }}}_{RR}=\frac{2m_{1}F(s)}{K_2(\eta )}\left[\frac{1}{\xi }{K}_1(\xi )+\frac{6}{{\xi }^2}{K}_2(\xi )\right]+2\left({\beta }^2-m_3\right)\frac{F(s)}{R^2}-\frac{6Q(s)}{R^4},\\ {\widehat{\mathsf{\Sigma }}}_{\theta \theta }=-\frac{2m_{1}F(s)}{K_2(\eta )}\left[\frac{1}{\xi }{K}_1(\xi )+\left(1+\frac{6}{{\xi }^2}\right)K_2(\xi )\right]+\frac{6Q(s)}{R^4},\\ {\widehat{\mathsf{\Sigma }}}_{R\theta }=\frac{4m_{1}F(s)}{K_2(\eta )}\left[\frac{1}{\xi }{K}_1(\xi )+\frac{3}{{\xi }^2}{K}_2({\xi }_1)\right]+\left({\beta }^2-m_3\right)\frac{F(s)}{R^2}-\frac{6Q(s)}{R^4},\end{array}$$

(45)

where the function F(s) and Q(s) are calculated via the stress boundary conditions at R = 1

$$F(s)=\frac{4S_0\eta K_2(\eta )}{s{H}_1(s)},Q(s)=-\frac{S_0}{6s}\frac{H_2(s)}{H_1(s)},$$

(46)

where

$$\begin{array}{l}{H}_1(s)=-4K_1(\eta )m_1-2\eta K_2(\eta )\left(m_3-{\beta }^2\right),\\ H_2(s)=12K_3(\eta )m_1-6\eta K_2(\eta )\left(m_3-{\beta }^2\right).\end{array}$$

As the temperature, displacements, and stresses induced by both modes are known, complete solutions are easily obtained by superposition.

5. Results and Discussions

In this section, the physical solutions to the temperature and stresses are obtained by applying the Stehfest method [48]. The parameters used in the following calculations are listed in

Table 1. Parameters for the calculations.

Parameters	Value
Pipe internal radius r0 and thickness δ0 (m)	0.0462, 0.01
Wellbore radius rw and height H (m)	0.0762, 4131
Injection rate Qd and pump out rate Qa (Kg/s)	23.0, 23.0
Injection and surface temperature Tin*, B0 (°C)	36.5, 27
Formation geothermal gradient A0 (°C/m)	0.047
Initial temperature T0* (°C)	A0z + B0
Fluid and rock specific heat cl, cr (J/(kg∙K))	4200, 790
Pipe specific heat cd (J/(kg∙K))	460
Fluid, rock ther. conductivity kl, kr (W/(m∙K))	0.68, 2.2
pipe thermal conductivity kd (W/(m∙K))	50
Fluid and rock mass density ρl, ρr (Kg/m3)	900, 2700
Pipe mass density ρd (Kg/m3)	7800
Fluid viscosity μ (Pa∙s)	0.0004
Thermal expansion coefficient α (1/K)	5.0 × 10−6
Max horizontal, principal stress σH (MPa)	−170
Min horizontal, principal stress σh (MPa)	−130
Vertical, principal stress σv (MPa)	−110
Shear modulus G (Pa)	1.5 × 1010
Poisson ratio v	0.25

unless otherwise specified.

5.1. Wellbore Temperature Responses

The temperature responses in the wellbore/reservoir system can be calculated based on Equations (30) and (31).

Figure 3 shows the outlet and bottom-hole temperature changes for a range of injection rates. In Figure 3a, it can be seen that the outlet temperature change increases quickly (in about 2 h) to a certain value and then change very slowly after that time. In addition, when the injection rate increases from 3 Kg/s to 23 Kg/s, there is only a less than 3 °C difference between the outlet temperature changes. However, by increasing the injection rate, the bottom-hole temperature can be decreased quickly and by a significant amount as shown in Figure 3b. Therefore, changing the injection rate can be an efficient way to control the bottom-hole temperature and therefore the stress conditions near the wellbore, as shown by these simulations.

In Figure 4, the outlet and bottom-hole temperature changes for different injection temperatures are displayed. The injection rate is kept as a constant of 23 Kg/s. The bottom-hole temperature change decreases by 6 °C when the injection fluid temperature varies from 56.5 °C to 16.5 °C, even after a significant injection period, as shown in Figure 4b. This demonstrates that, to manage bottom-hole temperature and bottom-hole stress conditions, it is more effective to control the injection rate, rather than control injection fluid temperature for a deep wellbore.

In contrast, for a constant injection rate, the outlet fluid temperature is determined to a certain degree by the injection fluid temperature. In addition to the heat carried by fluid flowing up the annulus, the region near the outlet undergoes a quick heat exchange with the injection fluid. Therefore, if the injection temperature is lower than the surface soil/rock temperature, the outlet region is in a cooled. This can be seen in Figure 4a, where the outlet fluid temperature change varies from about −6 °C to 31 °C when the injection fluid temperature changes from 16.5 °C to 56.5 °C and the initial surface temperature B₀ is fixed as 27 °C.

The fluid temperature profiles along the wellbore for two injection fluid temperatures (16.5 °C and 56.5 °C) are plotted in Figure 5. It should be mentioned that Z in the y axis denotes the dimensionless depth, i.e., Z = z/H. In both cases, the rock temperature in the lower portion of the wellbore decreases with time while that in the upper portion of the wellbore depends greatly on the injection fluid temperature. This makes the wellbore behave in a relatively complicated thermal manner and leads to the stress increasing or decreasing at different depths along the wellbore through heat exchange with the surrounding rock. Therefore, a transient analysis of coupled heat exchange for the W/R system is required.

5.2. Temperature Change Near the Wellbore in the Rock

Figure 6 shows the contours for the temperature changes around the circular well at different times. Here Z and R denote the dimensionless depth and radial distance from the center of the wellbore, respectively, i.e., Z = z/H and R = r/r_w. It can be seen that the near-well region at the bottom of the well is greatly impacted by the circulation. The low temperature area near the bottom of the well extends into the formation to a depth of one wellbore radius after 120 h. This cooled zone will induce thermo-elastic stresses that may be useful in helping initiate fractures and to enhance the permeability of certain natural fracture sets during stimulation.

It should be noted that, at the initial time of circulation, the upper part of the well is heated due to the hotter fluid flowing upwards, thus leading to enhanced compressive stress in the near-wellbore region. This may cause wellbore failure in terms of breakout. Therefore, some measure needs to be taken to prevent heat from loss into the formation by using casing with low thermal conductivity during the initial stage.

5.3. Comparisons of Near-Well Temperature and Stresses for a System Subject to Fixed and Variable Wellbore Temperature

In this section, a comparison of near-well temperature and stress change is made for two cases, i.e., one with constant wellbore temperature (Case 1, the wellbore wall experiences a step temperature change) and the other with varying temperature that is dependent on the fluid circulation history (Case 2).

Figure 7 shows the temperature and hoop stress variations in the σ_Hmax direction at the wellbore bottom for Case 1 where the bottom-hole temperature is assumed to take a step change and be equal to the bottom-hole temperature after circulating for five days, and Case 2 where the circulation rate is equal to 28 Kg/s. Obviously, for Case 1, the rock temperature starts from the same point as the temperature on the wall that is fixed and increases with the distance away from the wellbore at a specific time. However, for Case 2, the rock temperature at the near-well region decreases with circulation time. There exists a region where the temperature predicted from the circulation case is higher than that from the fixed wellbore temperature case before the time t = 5 days. If we increase the circulation time (to say 25 days), the rock temperature will further decrease, thus leading to a region where the temperature predicted from the circulation case will be lower than that from the constant wellbore temperature case.

Due to different temperature responses predicted from the above two cases, the stress responses in the rock formation near the wellbore are also different. In Figure 7b, it can be seen that for a very small circulation time, the difference of the hoop stresses near the wellbore predicted from the two cases is significant. With time elapsed, the region with significant stress difference expands. At the time t = 5 days, the hoop stresses predicted from both cases are very close as the temperature changes in both cases have a similar profile. After five days, the difference of the hoop stresses for both cases becomes larger as the bottom-hole temperature in the circulation case continues to decrease with time.

Therefore, based on the circulation time, the conventional wellbore stability model with constant well temperature assumptions may over- or under-estimate the breakdown pressure (or mud weight) of the wellbore.

6. Conclusions

This paper presents a transient analysis of temperature and stress changes inside the wellbore and in the formation during fluid circulation. In particular, the solutions are derived for simultaneously solving fluid circulation along the wellbore and elastic deformation in the formation, through a unique temperature profile. By using the Laplace transformation and superposition method, the analytical solutions in the Laplace space are obtained. Numerical inversion is carried out to obtain the solutions in the time domain. The following conclusions are drawn based on the numerical results:

• The fluid circulation rate plays a dominant role in the temperature evolution of the W/R system. The higher circulation rates (say cooling), the larger bottom-hole temperature change and thus the larger induced tensile stresses around the wellbore.
• The effect of the injection fluid temperature on the outlet rock temperature change occurs rapidly. The outlet rock temperature change reaches a value quickly and is almost unchanged after that time.
• The rock temperature of the upper part of the wellbore is determined mainly by the injection condition. It is possible for the upper open-hole section to develop breakouts due to the thermal stresses induced by the heating. Therefore, both wellbore cooling and heating should be taken into account during wellbore stability analysis.
• This work provides a more consistent prediction on the temperature evolution and stress distribution along the wellbore resulting from the variable well temperature profiles associated with fluid circulation, thus making possible a more accurate wellbore stability analysis. The analysis for fixed boundary condition may over- or under-estimate the stress conditions around the wellbore, thus leading to inaccurate prediction of the mud weight density required to maintain a stable well.
• Based on the first two points, varying circulation rates may be a more efficient way to manage bottom-hole temperature and bottom-hole stress conditions, rather than changing injection fluid temperature.