# New Analytical Method to Evaluate Casing Integrity during Hydraulic Fracturing of Shale Gas Wells

Authors: Jun Li; Xueli Guo; Gonghui Liu; Shuoqiong Liu; Yan Xi
Journal: Shock and Vibration
 (2019)

---

## Abstract

An accurate analysis of casing stress distribution and its variation regularities present several challenges during hydraulic fracturing of shale gas wells. In this paper, a new analytical mechanical-thermal coupling method was provided to evaluate casing stress. For this new method, the casing, cement sheath, and formation (CCF) system was divided into three parts such as initial stress field, wellbore disturbance field, and thermal stress field to simulate the processes of drilling, casing, cementing, and fracturing. The analytical results reached a good agreement with a numerical approach and were in-line with the actual boundary condition of shale gas wells. Based on this new model, the parametric sensitivity analyses of casing stress such as mechanical and geometry properties, operation parameters, and geostress were conducted during multifracturing. Conclusions were drawn from the comparison between new and existing models. The results indicated that the existing model underestimated casing stress under the conditions of the geostress heterogeneity index at the range of 0.5–2.25, the fracturing pressure larger than 25 MPa, and a formation with large elastic modulus or small Poisson’s ratio. The casing stress increased dramatically with the increase of in situ stress nonuniformity degree. The stress decreased first and then increased with the increase of fracturing pressure. Thicker casing, higher fluid temperature, and cement sheath with small modulus, large Poisson’s ratio, and thinner wall were effective to decrease the casing stress. This new method was able to accurately predict casing stress, which can become an alternative approach of casing integrity evaluation for shale gas wells.


---

## Body

## 1. Introduction

During the multistage fracturing process, fracturing fluids are pressed into a borehole with a high pump rate and pressure. The complex downhole environments—high pressure and large temperature variation—increase the risk of casing deformation. The volume fracturing technique effectively reconstructs shale reservoirs; however, frequent and serious casing deformation failures occur [1]. There were over 36 wells with casing deformation (including 112 horizontal wells by 2017) during fracturing process in some shale gas plays. Drilling tools were blocked, and serious deformation sections were abandoned before fracturing operation completion [2–4]. Therefore, casing integrity has always been the main issue in designing shale gas wells under harsh conditions, and accurately calculating casing stress becomes primary to guarantee the casing safety. A casing stress analysis presents several challenges regarding the drilling, completion, and fracturing phases of shale gas wells. Many scholars have predicted the casing stress by simplifying the actual situations.Analytical solutions have been developed under different conditions for casing-cement sheath-formation system. In-plane and out-of-plane analyses for the stress field around an internally pressurized, cased, cemented, and remotely loaded circular hole were developed [5]. Taking the in situ stresses and well trajectory into account, the mechanical model of casing in the directional well under in situ stresses is established [6]. Other scholars gave a comprehensively analytical solution of the stress distribution in a casing-cement-formation system under anisotropic in situ stresses. Fang et al. [7] developed a multilayer cemented casing system in the directional well under anisotropic formation to investigate the collapse resistance of casing under nonuniform in situ stress and anisotropic formation. Based on the stress function method, a three-dimensional model of the casing-cement sheath-formation (CCF) system was proposed subjected to linear crustal stress, and then an analytical solution of the model was obtained [8].A finite element method (FEM) had been proposed to achieve a better understanding of the ultimate collapse strength of casing [9]. Using this method, casing stress was better analyzed, subjected to external and nonuniform loading caused by void and pressure [10]. The results of a series of finite element studies of the cemented casing under a variety of stress conditions for both burst and collapse were presented, demonstrating the inadequacy of accepted design equations under many cemented conditions [11]. A purpose-built finite element model was applied to simulate the radial displacement of a casing string constrained within an outer wellbore [12]. Casing stress on the inner wall under the condition of elliptic casing was calculated and analyzed to improve the designing and correction of casing strength [13]. The main controlling factor to the plastic limit load of defective casing was analyzed by the combination of high temperature and high steam injection pressure [14].However, the above analytical or FEM solutions are obtained by setting the casing and cement together instantly at the beginning of the analysis. The strains induced by the initial stress are included in the model, which is not in-line with the actual situation. In fact, initial stress has already existed in the formation before wellbore excavation. The wellbore stress redistributed after removing the rocks that originally occupied the borehole volume and just affected the stress and displacement near wellbore zones. In view of this, more sophisticated solutions have been developed with additional parameters and appropriate assumptions regarding the drilling, completion, and fracturing phases.Pattillo and Kristiansen [15] developed a finite model starting with the virgin reservoir, which considered prewellbore depletion, drilling the wellbore, installation of casing and cement, and subsequent draw down. Behavior of various configurations during subsequent draw down permitted them to be ranked according to the life expectancy of the resulting completion. Gray et al. [16] presented a staged-finite element procedure during the well construction which considered sequentially the stress states and displacements at and near the wellbore. But, temperature, flow, and poroelasticity effects were not included. Mackay and Fontoura [17] carried a numerical analysis to determine the effect of salt creep and contacts amongst the materials before, during and after drilling the well, focusing on the drilling of the wellbore and on the hardening of the cement. Zhang et al. [18] used the staged-finite element modeling approach to simulate the well construction processes and injection cycle, using a “realistic” bottom-hole state of stress to simulate the microannuli generation by the tensile debonding of the cement-formation interface. Simone et al. [19] developed an analytical solution of single and double casing configurations to assess the stresses during the drilling, construction, and production phases. But the nonuniform boundary stress was not considered in this model. Liu et al. [20] presented an analytical method for evaluating the stress field within a casing-cement-formation system of oil/gas wells under anisotropic in situ stresses in the rock formation and uniform pressure within the casing. However, the temperature was not considered in this model.In this work, to evaluate the thermal and mechanical stresses of casing, an analytical model of casing-cement sheath-formation system was established considering wellbore construction. The boundary stresses in the wellbore coordinate system were obtained through three-dimensional rotations from principal in situ stresses. The casing stress was obtained by dividing the model into three parts such as initial stress field, disturbance stress field, and thermal stress field. The continuous homogeneous isotropy and linear elasticity were taken into consideration. Solutions were validated by a finite element method. Sensitivity analyses were conducted to estimate the influence of different factors on casing stress such as property of cement sheath and casing, fracturing pressure, fluid temperature, and initial geostress. Useful countermeasures were put forward to decrease casing stress during the fracturing operation.


## 2. Method and Basic Conditions Analysis

### 2.1. Method Comparison

#### 2.1.1. Existing Method

For existing method, the casing, cement sheath, and formation were set together at the beginning of analysis. Then, temperature boundary and loads were added to the system to calculate wellbore stress distribution as shown in Figure1. This method did not consider the process of wellbore construction and hydraulic fracturing. So, the wellbore stress distribution could not accurately be calculated under the condition of formation strain induced by the initial in situ stress.Figure 1
Loading process of the existing method. Outer boundary temperature,Tf; inner boundary temperature, Tn; normal stress, σx and σy; shear stress τxy.


#### 2.1.2. New Method

To exactly predict the wellbore stress distribution, a casing-cement sheath-formation model using a new analysis method was provided. The analysis process was divided into four steps (Figure2). Before drilling, the initial stresses such as normal stress and shear stress were loaded in the model and initial strains were produced. The constrained boundary conditions were assigned to the entire model which reached initial mechanical equilibrium. During drilling, the rock was excavated, causing stress concentration around the borehole. The drilling mud pressure Pm was applied to the internal face of the wellbore. During casing and cementing, casing and cement sheath elements were added simultaneously to the model and the cement hardening procedure was not considered. Initial stress and strain in casing and cement sheath were ignored. The cement slurry pressure Pc was applied to the internal casing wall. During fracturing, a low temperature Tn and the fracturing pressure Pf were assigned in the internal casing wall. The wellbore stress field was obtained under the condition of thermal-mechanical coupling.Figure 2
Loading process of the new method. Outer boundary temperature,Tf; inner boundary temperature, Tn; normal stress, σx and σy; shear stress, τxy; drilling mud pressure, Pm; cement slurry pressure, Pc; fracturing pressure, Pf.


### 2.2. Stress Transformation

The stress state and coordinate transformation system are shown in Figure3. The coordinate rotation processes from the principal in situ stress coordinate system to the wellbore coordinate system are shown as XHYhZy⟶X′Y′Z′⟶X″Y″Z″⟶X″′Y″′Z″′⟶XYZ: first, rotating anticlockwise φ around the Zv-axis and rotating clockwise β around the Zv-axis; second, rotating anticlockwise α around the Y″ axis (after the second time rotation); and finally, rotating anticlockwise 90° around the Z″′ axis (after the third time rotation).Figure 3
The coordinate rotation processes: (a) stress state of a horizontal well and (b) coordinate transformation system. Principal in situ stress coordinate system,XHYhZy; geodetic coordinate system, X′Y′Z′; coordinate system, X″Y″Z″ after the second time rotation; coordinate system, X″′Y″′Z″′ after the third time rotation; local wellbore coordinate system, XYZ; wellbore deviation angle, α; wellbore azimuth angle, β; the maximum horizontal stress azimuth angle, φ.In the principal stress coordinate system, the principal horizontal stress matrix is σ0, whose components are maximum principal stress, σH; the minimum principal stress, σh; and the overburden pressure, σv, shown in equation (2). The stress matrix in the wellbore coordinate system is σ. According to the right-hand rule, the direction cosine matrices rotating around X-axis, Y-axis, and Z-axis are shown in equation (1) [21]. The wellbore boundary stress is obtained by the three-dimensional rotations from the principal in situ stress coordinate system shown in equation (3):(1)Cxαx=1000cosαxsinαx0−sinαxcosαx,Cyαy=cosαy0−sinαy010sinαy0cosαy,Czαz=cosαzsinαz0−sinαzcosαz0001,(2)σ0=−σH000−σh000−σv,(3)σ=Cz90CyαCzβCzjσ0CzjTCzβTCyαTCz90T,where Cxαx, Cyαy, and Czαz are the coordinate rotation matrices of the x, y, and z directions; σ0 and σ are the stress matrices in the principal stress coordinate and local wellbore coordinate systems; and αx, αy, and αz are the rotation angles in a counterclockwise direction when looking towards the origin coordinate.


### 2.3. Basic Hypotheses

A thermo-pressure coupling model of casing-cement sheath-formation (CCF) system was established (Figure4). The boundary stresses of σx, σy, and τxy were obtained by using equation (3). Compared to the longitude of the well, the radial dimension was very small and a long cylindrical model was loaded by forces that were perpendicular to the axial line and did not vary in length. Both the geometric form of the object and the external loads exerted on the object did not change along the longitudinal (z-axis) direction and that the length of the object might be treated as an infinite one. There were no additional restrictions on the external loads. The wellbore stress was obtained under this kind of stress state called the plane strain problem [22, 23].Figure 4
CCF composite assembly. Formation boundary temperature,Tf; internal casing temperature, Tn; internal casing pressure, Pi; radius, Ri, i = 1, 2, 3, 4 present the radii of the internal casing wall, outer casing wall, internal wellbore, and formation boundary, respectively; the counterclockwise angle from the x-direction to the calculated point, θ.For simplicity, some assumptions have been made [24]:(1)
Geometry: the casing, cement, and borehole were concentric circles, which were assumed to be perfectly bonded to each other at each interface. The perfect bonding mathematically indicated that the continuity of radial stress and displacement was satisfied at each interface.(2)
Thermal effect: the stress induced by wellbore temperature variation was assumed to be steady state and the time effect was ignored.(3)
Material: to simplify the complex property of strong anisotropy and well-developed bedding planes of shale formation [25], the formation was assumed to be a linear elastic material with an infinite radius (R4 ⟶ ∞). The cement sheath was also a complex material. The 3D images revealed the evolution of a large connected pore network with characteristic widths on the micrometer scale as hydration proceeded [26]. It was assumed to be an elastic material neglecting the complex microstructure. It is generally accepted that the casing was an elastic-plastic material, and the casing yielding had nothing to do with hydrostatic pressure. In the model, it was assumed that the deformation of the casing was within the elastic range.


### 2.4. Stress Superposition

The boundary compression normal stresses (−σx, −σy) were decomposed into uniform stress (p0) and deviator stress (s) expressed as equation (4) in the Cartesian coordinate system. The other boundary stress was shear stress (τxy):(4)−σx00−σy=−p000−p0+−s00s,where the uniform stress p0=σx+σy/2 and the deviator stress s=σx−σy/2.According to the basic hypotheses in Section2.2, the deformation history of all phases in the CCF system was independent with each other. The principal of the linear superposition for the stress was applied as shown in the following equation:(5)σstressi=σUniform−stress′i+σDeviator−stress″i+σShear−stress″′i+σThermal−stressTi.The stress distribution of the thermal-pressure coupling model around a wellbore was decomposed into four parts as shown in Figure5. The stress field of the first part was induced by the uniform stress and inner casing pressure. The second one was induced by the deviator stress, and the third one was induced by the shear stress. Thermal stress of the fourth part was induced by the temperature variation.Figure 5
Stress decompositions. Inner casing pressure,Pi; thermal stress, σT.It was convenient to convert the Cartesian coordinate system into the polar coordinate system to calculate wellbore stress. The normal boundary stresses in the polar coordinate system were expressed as follows under the conditions of the infinite outer boundary radius ofR4:(6)σr′30r=∞=−p0−scos2θ,τrθ′30r=∞=ssin2θ,where σr′30 and τrθ′30 are the initial normal and shear stresses in the formation, respectively.Since temperature and stress were coupled, the stress distribution around a cased wellbore induced by temperature variation was hard to solve in the closed form. However, the steady-state condition made the temperature and stress decouple and the problem analytically solvable [27].


## 2.1. Method Comparison

### 2.1.1. Existing Method

For existing method, the casing, cement sheath, and formation were set together at the beginning of analysis. Then, temperature boundary and loads were added to the system to calculate wellbore stress distribution as shown in Figure1. This method did not consider the process of wellbore construction and hydraulic fracturing. So, the wellbore stress distribution could not accurately be calculated under the condition of formation strain induced by the initial in situ stress.Figure 1
Loading process of the existing method. Outer boundary temperature,Tf; inner boundary temperature, Tn; normal stress, σx and σy; shear stress τxy.


### 2.1.2. New Method

To exactly predict the wellbore stress distribution, a casing-cement sheath-formation model using a new analysis method was provided. The analysis process was divided into four steps (Figure2). Before drilling, the initial stresses such as normal stress and shear stress were loaded in the model and initial strains were produced. The constrained boundary conditions were assigned to the entire model which reached initial mechanical equilibrium. During drilling, the rock was excavated, causing stress concentration around the borehole. The drilling mud pressure Pm was applied to the internal face of the wellbore. During casing and cementing, casing and cement sheath elements were added simultaneously to the model and the cement hardening procedure was not considered. Initial stress and strain in casing and cement sheath were ignored. The cement slurry pressure Pc was applied to the internal casing wall. During fracturing, a low temperature Tn and the fracturing pressure Pf were assigned in the internal casing wall. The wellbore stress field was obtained under the condition of thermal-mechanical coupling.Figure 2
Loading process of the new method. Outer boundary temperature,Tf; inner boundary temperature, Tn; normal stress, σx and σy; shear stress, τxy; drilling mud pressure, Pm; cement slurry pressure, Pc; fracturing pressure, Pf.


## 2.1.1. Existing Method

For existing method, the casing, cement sheath, and formation were set together at the beginning of analysis. Then, temperature boundary and loads were added to the system to calculate wellbore stress distribution as shown in Figure1. This method did not consider the process of wellbore construction and hydraulic fracturing. So, the wellbore stress distribution could not accurately be calculated under the condition of formation strain induced by the initial in situ stress.Figure 1
Loading process of the existing method. Outer boundary temperature,Tf; inner boundary temperature, Tn; normal stress, σx and σy; shear stress τxy.


## 2.1.2. New Method

To exactly predict the wellbore stress distribution, a casing-cement sheath-formation model using a new analysis method was provided. The analysis process was divided into four steps (Figure2). Before drilling, the initial stresses such as normal stress and shear stress were loaded in the model and initial strains were produced. The constrained boundary conditions were assigned to the entire model which reached initial mechanical equilibrium. During drilling, the rock was excavated, causing stress concentration around the borehole. The drilling mud pressure Pm was applied to the internal face of the wellbore. During casing and cementing, casing and cement sheath elements were added simultaneously to the model and the cement hardening procedure was not considered. Initial stress and strain in casing and cement sheath were ignored. The cement slurry pressure Pc was applied to the internal casing wall. During fracturing, a low temperature Tn and the fracturing pressure Pf were assigned in the internal casing wall. The wellbore stress field was obtained under the condition of thermal-mechanical coupling.Figure 2
Loading process of the new method. Outer boundary temperature,Tf; inner boundary temperature, Tn; normal stress, σx and σy; shear stress, τxy; drilling mud pressure, Pm; cement slurry pressure, Pc; fracturing pressure, Pf.


## 2.2. Stress Transformation

The stress state and coordinate transformation system are shown in Figure3. The coordinate rotation processes from the principal in situ stress coordinate system to the wellbore coordinate system are shown as XHYhZy⟶X′Y′Z′⟶X″Y″Z″⟶X″′Y″′Z″′⟶XYZ: first, rotating anticlockwise φ around the Zv-axis and rotating clockwise β around the Zv-axis; second, rotating anticlockwise α around the Y″ axis (after the second time rotation); and finally, rotating anticlockwise 90° around the Z″′ axis (after the third time rotation).Figure 3
The coordinate rotation processes: (a) stress state of a horizontal well and (b) coordinate transformation system. Principal in situ stress coordinate system,XHYhZy; geodetic coordinate system, X′Y′Z′; coordinate system, X″Y″Z″ after the second time rotation; coordinate system, X″′Y″′Z″′ after the third time rotation; local wellbore coordinate system, XYZ; wellbore deviation angle, α; wellbore azimuth angle, β; the maximum horizontal stress azimuth angle, φ.In the principal stress coordinate system, the principal horizontal stress matrix is σ0, whose components are maximum principal stress, σH; the minimum principal stress, σh; and the overburden pressure, σv, shown in equation (2). The stress matrix in the wellbore coordinate system is σ. According to the right-hand rule, the direction cosine matrices rotating around X-axis, Y-axis, and Z-axis are shown in equation (1) [21]. The wellbore boundary stress is obtained by the three-dimensional rotations from the principal in situ stress coordinate system shown in equation (3):(1)Cxαx=1000cosαxsinαx0−sinαxcosαx,Cyαy=cosαy0−sinαy010sinαy0cosαy,Czαz=cosαzsinαz0−sinαzcosαz0001,(2)σ0=−σH000−σh000−σv,(3)σ=Cz90CyαCzβCzjσ0CzjTCzβTCyαTCz90T,where Cxαx, Cyαy, and Czαz are the coordinate rotation matrices of the x, y, and z directions; σ0 and σ are the stress matrices in the principal stress coordinate and local wellbore coordinate systems; and αx, αy, and αz are the rotation angles in a counterclockwise direction when looking towards the origin coordinate.


## 2.3. Basic Hypotheses

A thermo-pressure coupling model of casing-cement sheath-formation (CCF) system was established (Figure4). The boundary stresses of σx, σy, and τxy were obtained by using equation (3). Compared to the longitude of the well, the radial dimension was very small and a long cylindrical model was loaded by forces that were perpendicular to the axial line and did not vary in length. Both the geometric form of the object and the external loads exerted on the object did not change along the longitudinal (z-axis) direction and that the length of the object might be treated as an infinite one. There were no additional restrictions on the external loads. The wellbore stress was obtained under this kind of stress state called the plane strain problem [22, 23].Figure 4
CCF composite assembly. Formation boundary temperature,Tf; internal casing temperature, Tn; internal casing pressure, Pi; radius, Ri, i = 1, 2, 3, 4 present the radii of the internal casing wall, outer casing wall, internal wellbore, and formation boundary, respectively; the counterclockwise angle from the x-direction to the calculated point, θ.For simplicity, some assumptions have been made [24]:(1)
Geometry: the casing, cement, and borehole were concentric circles, which were assumed to be perfectly bonded to each other at each interface. The perfect bonding mathematically indicated that the continuity of radial stress and displacement was satisfied at each interface.(2)
Thermal effect: the stress induced by wellbore temperature variation was assumed to be steady state and the time effect was ignored.(3)
Material: to simplify the complex property of strong anisotropy and well-developed bedding planes of shale formation [25], the formation was assumed to be a linear elastic material with an infinite radius (R4 ⟶ ∞). The cement sheath was also a complex material. The 3D images revealed the evolution of a large connected pore network with characteristic widths on the micrometer scale as hydration proceeded [26]. It was assumed to be an elastic material neglecting the complex microstructure. It is generally accepted that the casing was an elastic-plastic material, and the casing yielding had nothing to do with hydrostatic pressure. In the model, it was assumed that the deformation of the casing was within the elastic range.


## 2.4. Stress Superposition

The boundary compression normal stresses (−σx, −σy) were decomposed into uniform stress (p0) and deviator stress (s) expressed as equation (4) in the Cartesian coordinate system. The other boundary stress was shear stress (τxy):(4)−σx00−σy=−p000−p0+−s00s,where the uniform stress p0=σx+σy/2 and the deviator stress s=σx−σy/2.According to the basic hypotheses in Section2.2, the deformation history of all phases in the CCF system was independent with each other. The principal of the linear superposition for the stress was applied as shown in the following equation:(5)σstressi=σUniform−stress′i+σDeviator−stress″i+σShear−stress″′i+σThermal−stressTi.The stress distribution of the thermal-pressure coupling model around a wellbore was decomposed into four parts as shown in Figure5. The stress field of the first part was induced by the uniform stress and inner casing pressure. The second one was induced by the deviator stress, and the third one was induced by the shear stress. Thermal stress of the fourth part was induced by the temperature variation.Figure 5
Stress decompositions. Inner casing pressure,Pi; thermal stress, σT.It was convenient to convert the Cartesian coordinate system into the polar coordinate system to calculate wellbore stress. The normal boundary stresses in the polar coordinate system were expressed as follows under the conditions of the infinite outer boundary radius ofR4:(6)σr′30r=∞=−p0−scos2θ,τrθ′30r=∞=ssin2θ,where σr′30 and τrθ′30 are the initial normal and shear stresses in the formation, respectively.Since temperature and stress were coupled, the stress distribution around a cased wellbore induced by temperature variation was hard to solve in the closed form. However, the steady-state condition made the temperature and stress decouple and the problem analytically solvable [27].


## 3. Stress Distribution around Wellbore

### 3.1. Stress Induced by Uniform Stress

Under the condition of the uniform internal pressure and external stress, the stress and displacement in a thin wall cylinder were obtained by using the following equations shown in Figure6:(7)ur′i=12Gi1−2μiAi′r+Ci′1rqi+12Gi1−2μiBi′r+Ci′1rqi+1−rεr′i0,(8)σr′i=Ai′−Ci′1r2qi+Bi′−Ci′1r2qi+1,σθ′i=Ai′+Ci′1r2qi+Bi′+Ci′1r2qi+1,where σr′i,σθ′i,andur′i are the radial stress, tangential stress, and radial displacement, respectively; σr′i0,σθ′i0,andur′i0 are the initial radial stress, tangential stress, and displacement; Ei is the material elastic modulus; μi is the material Poisson’s ratio; Gi=Ei/1+μi/2 is the material shear modulus; qi, qi+1 were the interfacial pressure, positive in the radial increase direction, i = 1, 2, 3 represented the casing, cement sheath, and formation, respectively; Ai′=Ri2/Ri+12−Ri2; Bi′=Ri+12/Ri+12−Ri2; Ci′=Ai′Ri+12 are the constants, Ri (i = 1, 2, 3, 4) is the radii of internal casing wall, external casing wall, external cement sheath wall, and formation boundary, respectively.Figure 6
Stress induced by the uniform stress.


#### 3.1.1. Formation Stress

Before drilling the borehole, the initial geostress field already existed in the formation. When the rock was removed from the borehole, the wellbore stress field redistributed to produce a disturbance field, which only affected the near-wellbore zones [28]. So, the model was decomposed into two parts such as the original field and the disturbance field. The original field had initial stress and displacement. The disturbance field was induced by drilling a wellbore and mud pressure. In view of this, the actual stress field of F1 in the formation induced by the uniform stress, and internal pressure was decomposed into three parts as shown in Figure 7. They were the original stress field of A1, the excavation disturbance field of B1 induced by drilling of a wellbore and the interface disturbance field of C1 induced by the fluid column pressure.Figure 7
Wellbore stress components under the condition of uniform stress.In the polar coordinate system, the initial conditions ofA1 were σr30=σθ30=−p0, εr′30=p0/2G31−2μ3, and boundary stress conditions of B1 and C1 were σr′3r=R3,B1=p0 and σr′3r=R3,C1=−p3′. Substituting the initial and boundary conditions into equations (7) and (8), the displacement and stress in formation were obtained as shown below. Because R4 approached to infinity, A3′=0 and C3′=R32 were obtained:(9)ur′3r=12G3R32rp3′−p0,σr′3r=−p0−R32r2p3′−p0,σθ′3r=−p0+R32r2p3′−p0,where ur′3 was the radial displacement in formation and ur′3 and σθ′3 were the radial and tangential stresses in formation, respectively.


#### 3.1.2. Casing and Cement Sheath Stress

The pressures at casing-cement sheath interface and cement sheath-formation interface werep2′ and p3′, respectively, (Figure 8). The initial stresses of the casing and the cement sheath were σr′i0=σθ′i0=0 and εr′i0=0. The boundary stress conditions were σr′2r=R3=−p3′ and σr′2r=R2=−p2′, σr′1r=R1=−pi.Figure 8
Interface pressures induced by uniform stress.Substituting these initial and boundary conditions in equations (7) and (8), the displacement and stress in casing and cement sheath were obtained as follows. Subscripts 1 and 2 represent the casing and cement sheath, respectively:(10)ur′1=12G11−2μ1A1′r+C1′1rpi+12G11−2μ1B1′r+C1′1r−p2′,ur′2=12G21−2μ2A2′r+C2′1rp2′+12G21−2μ2B2′r+C2′1r−p3′,(11)σr′1=A1′−C1′1r2pi+B1′−C1′1r2−p2′,σθ′1=A1′+C1′1r2p2′+B1′+C1′1r2−p2′,σr′2=A2′−C2′1r2p2′+B2′−C2′1r2−p3′,σθ′2=A2′+C2′1r2p2′+B2′+C2′1r2−p3′.According to the hypotheses that cement sheath-formation interface and casing-cement sheath interface were perfectly bonded to each other, the interfacial displacement continuity conditions were expressed in the following equation:(12)ur′1r=R2=ur′2r=R2,ur′2r=R3=ur′3r=R3.Substituting equation (10) into equation (12), the binary equations were obtained as(13)AAp2′−BBp3′=CCpi,DDp2′−EEp3′=FFp0,where(14)AA=12G11−2μ1B1′R2+C11R2+12G21−2μ2A2′R2+C2′1R2,BB=12G11−2μ1B1′R2+C11R2,CC=12G11−2μ1A1′r+C1′1R2,DD=12G21−2μ2A2′R3+C2′1R3,EE=12G11−2μ1B1′R+3C11R3+12G3R3,FF=−R32G3.The interfacial pressuresp2′ and p3′ could be calculated by using equation (14). Substituting them into equation (11), the stresses induced by uniform stress were obtained subsequently.


### 3.2. Stress Induced by Deviator Stress

The deviator stress boundary conditions are shown in Figure9. To calculate the stress distribution induced by deviator stress, the stress function was defined as(15)ϕ=Ai″r4+Bi″r2+Ci″+Di″r2cos2θ.Figure 9
Stress induced by deviator stress. Outer stress,scos2θ; interface pressure, pcos2θ.The stress and strain under the condition of nonuniform stress are(16)σr″i=−2Bi″+4Ci″r2+6Di″r4cos2θ,σθ″i=12Ai″r2+2Bi″+6Di″r4cos2θ,τrθ″i=6Ai″r2+2Bi″−2Ci″r2−6Di″r4sin2θ,εr″i=−1+μiEi12υiAi″r2+2Bi″+1−μi4Ci″r2+6Di″r4cos2θ−εr″i0,εθ″i=1+μiEi121−μiAi″r2+2Bi″+μi4Ci″r2+6Di″r4cos2θ−εθ″i0,where σr″iandσθ″i are the radial and tangential stresses; εr″i0andεθ″i0 are the initial radial and tangential strains; and Ai″,Bi″,Ci″,andDi″ were the constants, i = 1, 2, 3 represented the casing, cement sheath, and formation.From the geometric equations,(17)εr″i=∂ur″i∂r,εθ″i=1r∂uθ″i∂θ+ur″ir.The radial displacementur″i and tangential displacement uθ″i were obtained as(18)ur″i=−1+μiEi4μiAi″r3+2Bi″r−1−μi4Ci″r−2Di″r3cos2θ−rεr″i0,uθ″i=1+μi2Ei43−2μiAi″r3+4Bi″r−1−2μi4Ci″r+4Di″r3sin2θ+r∫εr″i0−εθ″i0dθ,εr″i0=1+μiEi1−μiσr″i0−μiσθ″i0,εθ″i0=1+μiEi1−μiσθ″i0−μiσr″i0,where σr″i0andσθ″i0 were the initial radial and tangential stresses.


#### 3.2.1. Formation Stress

Similar to that of uniform stress, the actual stress field,F2 in the strata induced by the nonuniform stress, was decomposed into three parts: the original stress field, A2; the disturbance field, B2 induced by the wellbore excavation; and the interface pressure, C2 induced by the interface pressure (Figure 10).Figure 10
Formation stress components under the nonuniform stress condition.In the polar coordinate system, initial stresses wereσr″30=−scos2θ, σθ″30=scos2θ, and τrθ″30=ssin2θ; initial strains were εr″30=−1+μ3/E3⋅scos2θ and εθ″30=1+μ3/E3⋅scos2θ; and the boundary stresses were σr″3r=∞=−scos2θ, σθ″3=scos2θ, and τrθ″3r=∞=ssin2θ. Substituting the initial and boundary conditions into (14) and (15), it was obtained that A3″=0andB3″=S/2. The displacements and stresses in formation were expressed as shown in the following equations:(19)ur″3=−1G3−1−μ32C3″r−D3″r3cos2θ,uθ″3=1G3−1−2μ3C3″r+D3″r3sin2θ,(20)σr″3=−s+4C3″r2+6D3″r4cos2θ,σθ″3=s+6D3″r4cos2θ,τrθ″3=s−2C3″r2−6D3″r4sin2θ.


#### 3.2.2. Casing and Cement Sheath Stress

For casing and cement sheath in the polar coordinate system, initial stresses wereσr″i0=σθ″i0=0 and initial strains were εr″i0=εθ″i0=0. Substituting the initial and boundary conditions into equations (14) and (15), the displacements and stresses were obtained as follows:(21)ur″i=−1Gi2μiAi″r3+Bi″r−1−μi2Ci″r−Di″r3cos2θ,uθ″i=1Gi3−2μiAi″r3+Bi″r−1−2μiCi″r+Di″r3sin2θ,(22)σr″i=−2Bi″+4Ci″r2+6Di″r4cos2θ,σθ″i=12Ai″r2+2Bi″+6Di″r4cos2θ,τrθ″i=6Ai″r2+2Bi″−2Ci″r2−6Di″r4sin2θ.The interfacial displacement and stress continuity and boundary conditions were expressed in the following equation:(23)σr″1r=R1=0,τrθ″1r=R1=0,σr″1r=R2=σr″2r=R2,τrθ″1r=R2=τrθ″2r=R2,σr″2r=R3=σr″3r=R3,τrθ″2r=R3=τrθ″3r=R3,ur″1r=R2=ur″2r=R2,uθ″1r=R2=uθ″2r=R2,ur″2r=R3=ur″3r=R3,uθ″2r=R3=uθ″3r=R3.Substituting equations (20)–(23) into the following equation, equations were obtained as(24)−2B1″−4C1″R12−6D1″R14cos2θ=0,6R12A1″+2B1″−2R12C1″−6R14D1″sin2θ=0,(25)−2B1″+4C1″R22+6D1″R24+2B2″+4C2″R22+6D2″R24=0,6A1″R22+2B1″−2C1″R22−6D1″R24−6A2″R22+2B2″−2C2″R22−6D2″R24=0,(26)−B2″+2C2″R32+3D2″R34+2C3″R32+3D3″R34=−s,3A2″R32+B2″−C2″R32−3D2″R34+C3″R32+3D3″R34=s,(27)−1G12μ1R23A1″+R2B1″−1−μ12R2C1″−1R23D1″+1G22μ2R23A2″+R2B2″−1−μ22R2C2″−1R23D2″=0,1G13−2μ1R23A1″+R2B1″−1−2μ11R2C1″+1R23D1″−1G23−2μ2R23A2″+R2B2″−1−2μ21R2C2″+1R23D2″=0,(28)−1G22μ2R33A2″+R3B2″−1−μ22R3C2″−1R33D2″+1G3−1−μ32R3C3″−1R33D3″=0,1G23−2μ2R33A2″+R3B2″−1−2μ21R3C2″+1R33D2″−1G3−1−2μ31R3C3″+1R33D3″=0,where the constants of A1″,B1″,C1″,D1″,A2″,B2″,C2″,D2″,C3″,andD3″ were calculated by the total 10 equations in equations (24)–(28). Then, wellbore stress distribution induced by deviator stress was obtained by substituting these 10 constants and A3″andB3″ into equations (19)–(22).


### 3.3. Stress Induced by Shear Stress

The stress induced by shear stress wasur″′i,uθ″′i,σr″′i,σθ″′i,andτrθ″′i, i = 1, 2, 3 represented the casing, cement sheath, and formation, respectively, (Figure 11). The angle of Ω between σx and x-direction was calculated by using equation (29). Then, the principal stresses were obtained as follows [29]:(29)Ω=12arctan−2τxyσx−σy=π4,(30)σπ/4=τxy,σ−π/4=−τxy.Figure 11
Stress distribution induced by shear stress.It could be found that the stress distribution induced by shear stress was similar with that by deviator stress when counterclockwise rotating the angle ofπ/4. The stresses and displacements were obtained by substituting the reference variable θ=θ′−π/4 into the stress induced by deviator stress discussed in Section 3.2.


### 3.4. Stress Induced by Temperature Variation

The thermal field was obtained by using the steady temperature distribution model to calculate the thermal stress. When fracturing fluids were pumped into a wellbore with a high pump rate, they were always in the turbulent state. The heat transfer coefficient between casing and fluid was calculated using the Marshall model [30] shown in the following equation:(31)h=StkmD=0.0107kmDρaDeff4Q/πD2K3n+1/4nn32Q/πD3n−10.67K3n+1/4nn32Q/πD3n−1Cmkm0.33,where h is the heat transfer coefficient (W·m−2·°C−1), St is the Stanton number, Pr is the Prandtl number, Reg is the Reynolds number, μw,app is the fluid apparent viscosity, D is the inner diameter (m), Deff is the equivalent diameter (m), ρa is the fluid density (kg·m−3), n is the liquidity index, K is the consistency coefficient (Pa·sn), v is the fluid velocity, Q is the fracturing pump rate (m3·min−1), km is the coefficient of heat conductivity (W·m−1·°C−1), and Cm is the fluid specific heat capacity (J·kg−1·°C−1).The temperature distribution among casing, cement sheath, and formation is shown in Figure12. In the cylindrical coordinate system of CCF, the differential equation of steady heat conduction of the cylinder is expressed as [31](32)d2Tdr2+1rdTdr=0.Figure 12
The distribution of interface temperatures.Temperature field distribution solutions were obtained according to integral and boundary conditionskidT/dr=hnTi−Tn,Tr=Ri=Ti,Tr=Ri+1=Ti+1 shown in the following equation:(33)Tir=AiTlnr+BiT,(34)AiT=Ti+1−TnlnRi+1/Ri+i2/2−5i/2+3ki/Rih,BiT=TnlnRi+1−Ti+1lnRi+i2/2−5i/2+3ki/RihTi+1lnRi+1/Ri+i2/2−5i/2+3ki/Rih,where Ti is the temperature (°C); Ti is the temperature at the interface (°C); Tn is the fluid temperature (°C); ki is the material thermal conductivity (W·m−1·°C−1); Ri is the radius (m); and AiTandBiT were the constants, i = 1, 2, 3 represented casing, cement sheath, and formation, respectively.The heat flow density continuity conditions were expressed as(35)kidTirdrr=Ri+1=ki+1dTi+1rdrr=Ri+1.The temperatures at interfaces of casing-cement sheath and cement sheath-formation system were defined asT2 and T3 and were calculated by using the following equation:(36)1+β1T2−β1T3=T1,−T2+1+β2T3=β2T4,where(37)β1=k2k1lnR2/R1+k1/R1hlnR3/R2,β2=k3k2lnR3/R2lnR4/R3.Interfacial temperature ofTi was obtained by solving equation (36). The steady-state temperature field around the wellbore could be calculated by substituting Ti into equations (33) and (34). According to thermal elastic mechanics, constitutive equations for a plane strain problem were expressed as(38)εrT=1+μiEi1−μiσrT−μiσθT+1+μiαiT,εθT=1+μiEi1−μiσθT−μiσrT+1+μiαiT,εzT=0,γrθT=21+μiEiτrθT.The actual thermal stress field,F3 in the strata induced by the temperature changes, was decomposed into two parts: the original stress field, A3, and the disturbance field, B3 induced by the temperature variation shown in Figure 13.Figure 13
Thermal stress field. (a) Formation stress components. (b) Interface pressures.pi is the interface pressure, i = 2, 3 represented the casing-cement sheath interface and cement sheath-formation interface.
(a)
(b)The initial stresses wereσrTi0=σθTi0=0, and the initial strains were εrTi0=εθTi0=0. The stresses and displacements induced by thermal variations were expressed as(39)urTi=1+μi1−μiαir∫RirrΔTidr+C1Tir+C2Tir−rεrTi0,(40)σrTi=−αiEi1−μi1r2∫RirrΔTidr+Ei1+μiC1Ti1−2μi−C2Tir2,σθTi=αiEi1−μi1r2∫RirrΔTidr+Ei1+μiC1Ti1−2μi+C2Tir2−αiEi1−μiΔTi,where C1TiandC2Ti are the constants, σrTiandσθTi are the radial and tangential stresses (Pa), urTi is the radial displacement (m), ΔTi is the temperature changes (°C), pi is the interface pressure (Pa), and αi is the material thermal expansion coefficient, i = 1, 2, 3 represented casing, cement sheath, and formation, respectively.The temperatures were known, and the boundary was free at internal casing and external formation. So, radial stress at inner and outer boundaries equals to zero, and radial displacement at the outer boundary equals to zero as well. The boundary and interfacial displacement continuity conditions were expressed as(41)urT1r=R2=urT2r=R2,urT2r=R3=urT3r=R3,σrT1r=R2=σrT2r=R2,σrT2r=R3=σrT3r=R3,σrT1r=R1=0,σrT3r=R4=0.Substituting equations (39) and (40) into the following equation, the equations were obtained as(42)C1T1R2+C2T1R2−C1T2R2−C2T2R2=−1+μ11−μ1α1R2∫R1R2rT1dr,C1T2R3+C2T2R3−C1T3R3−C2T3R3=1+μ21−μ2α2R3∫R2R3rT2dr,(43)E11+μ1C1T11−2μ1−C2T1R22−E21+μ2C1T21−2μ2−C2T2R22=α1E11−μ11R22∫R1R2rΔT1dr,E21+μ2C1T21−2μ2−C2T2R32−E31+μ3C1T31−2μ3−C2T3R32=α2E21−μ21R32∫R2R3rΔT2dr,(44)C11R1+C21R1=0,C13R4+C23R4=−1+μ21−μ2α2R4∫R3R4rTdr.The constants ofC1T1,C2T1,C1T2,C2T2,C1T3,andC2T3 were obtained by equations (42)–(44). The wellbore stress was obtained by substituting these constants into equation (40).The total stresses were obtained using the following equation:(45)σri=σr′i+σr″i+σr″′i+σrTi,σθi=σθ′i+σθ″i+σθ″′i+σθTi,σzi=μiσri+σθi,τrθi=τrθ″i+τrθ″′i,where σri is the radial stress, σθi is the tangential stress, σzi is the axial stress, and τrθi is the shear stress.


### 3.5. Estimation of Wellbore Integrity

It is generally accepted that the yield of isotropic material such as casing has nothing to do with hydrostatic pressure, while hydrostatic pressure is not considered in von Mises yield criterion. So, this criterion was adopted to determine the casing failure:(46)fJ2,k=J2−k=0,J2=16σ11−σ222+σ22−σ332+σ33−σ112+σ122+σ232+σ312,where J2 is the second stress partial tensor; k is the critical value of failure; and σij is the stress components, i, j = 1, 2, 3 represented the three directions of the system, respectively.For uniaxial tension,J2=σ/3, the von Mises stress could be expressed as follows in the polar coordinate:(47)σMises=12σr−σθ2+σθ−σz2+σz−σr2+3τrθ2+3τθz2+3τzr2.


## 3.1. Stress Induced by Uniform Stress

Under the condition of the uniform internal pressure and external stress, the stress and displacement in a thin wall cylinder were obtained by using the following equations shown in Figure6:(7)ur′i=12Gi1−2μiAi′r+Ci′1rqi+12Gi1−2μiBi′r+Ci′1rqi+1−rεr′i0,(8)σr′i=Ai′−Ci′1r2qi+Bi′−Ci′1r2qi+1,σθ′i=Ai′+Ci′1r2qi+Bi′+Ci′1r2qi+1,where σr′i,σθ′i,andur′i are the radial stress, tangential stress, and radial displacement, respectively; σr′i0,σθ′i0,andur′i0 are the initial radial stress, tangential stress, and displacement; Ei is the material elastic modulus; μi is the material Poisson’s ratio; Gi=Ei/1+μi/2 is the material shear modulus; qi, qi+1 were the interfacial pressure, positive in the radial increase direction, i = 1, 2, 3 represented the casing, cement sheath, and formation, respectively; Ai′=Ri2/Ri+12−Ri2; Bi′=Ri+12/Ri+12−Ri2; Ci′=Ai′Ri+12 are the constants, Ri (i = 1, 2, 3, 4) is the radii of internal casing wall, external casing wall, external cement sheath wall, and formation boundary, respectively.Figure 6
Stress induced by the uniform stress.


### 3.1.1. Formation Stress

Before drilling the borehole, the initial geostress field already existed in the formation. When the rock was removed from the borehole, the wellbore stress field redistributed to produce a disturbance field, which only affected the near-wellbore zones [28]. So, the model was decomposed into two parts such as the original field and the disturbance field. The original field had initial stress and displacement. The disturbance field was induced by drilling a wellbore and mud pressure. In view of this, the actual stress field of F1 in the formation induced by the uniform stress, and internal pressure was decomposed into three parts as shown in Figure 7. They were the original stress field of A1, the excavation disturbance field of B1 induced by drilling of a wellbore and the interface disturbance field of C1 induced by the fluid column pressure.Figure 7
Wellbore stress components under the condition of uniform stress.In the polar coordinate system, the initial conditions ofA1 were σr30=σθ30=−p0, εr′30=p0/2G31−2μ3, and boundary stress conditions of B1 and C1 were σr′3r=R3,B1=p0 and σr′3r=R3,C1=−p3′. Substituting the initial and boundary conditions into equations (7) and (8), the displacement and stress in formation were obtained as shown below. Because R4 approached to infinity, A3′=0 and C3′=R32 were obtained:(9)ur′3r=12G3R32rp3′−p0,σr′3r=−p0−R32r2p3′−p0,σθ′3r=−p0+R32r2p3′−p0,where ur′3 was the radial displacement in formation and ur′3 and σθ′3 were the radial and tangential stresses in formation, respectively.


### 3.1.2. Casing and Cement Sheath Stress

The pressures at casing-cement sheath interface and cement sheath-formation interface werep2′ and p3′, respectively, (Figure 8). The initial stresses of the casing and the cement sheath were σr′i0=σθ′i0=0 and εr′i0=0. The boundary stress conditions were σr′2r=R3=−p3′ and σr′2r=R2=−p2′, σr′1r=R1=−pi.Figure 8
Interface pressures induced by uniform stress.Substituting these initial and boundary conditions in equations (7) and (8), the displacement and stress in casing and cement sheath were obtained as follows. Subscripts 1 and 2 represent the casing and cement sheath, respectively:(10)ur′1=12G11−2μ1A1′r+C1′1rpi+12G11−2μ1B1′r+C1′1r−p2′,ur′2=12G21−2μ2A2′r+C2′1rp2′+12G21−2μ2B2′r+C2′1r−p3′,(11)σr′1=A1′−C1′1r2pi+B1′−C1′1r2−p2′,σθ′1=A1′+C1′1r2p2′+B1′+C1′1r2−p2′,σr′2=A2′−C2′1r2p2′+B2′−C2′1r2−p3′,σθ′2=A2′+C2′1r2p2′+B2′+C2′1r2−p3′.According to the hypotheses that cement sheath-formation interface and casing-cement sheath interface were perfectly bonded to each other, the interfacial displacement continuity conditions were expressed in the following equation:(12)ur′1r=R2=ur′2r=R2,ur′2r=R3=ur′3r=R3.Substituting equation (10) into equation (12), the binary equations were obtained as(13)AAp2′−BBp3′=CCpi,DDp2′−EEp3′=FFp0,where(14)AA=12G11−2μ1B1′R2+C11R2+12G21−2μ2A2′R2+C2′1R2,BB=12G11−2μ1B1′R2+C11R2,CC=12G11−2μ1A1′r+C1′1R2,DD=12G21−2μ2A2′R3+C2′1R3,EE=12G11−2μ1B1′R+3C11R3+12G3R3,FF=−R32G3.The interfacial pressuresp2′ and p3′ could be calculated by using equation (14). Substituting them into equation (11), the stresses induced by uniform stress were obtained subsequently.


## 3.1.1. Formation Stress

Before drilling the borehole, the initial geostress field already existed in the formation. When the rock was removed from the borehole, the wellbore stress field redistributed to produce a disturbance field, which only affected the near-wellbore zones [28]. So, the model was decomposed into two parts such as the original field and the disturbance field. The original field had initial stress and displacement. The disturbance field was induced by drilling a wellbore and mud pressure. In view of this, the actual stress field of F1 in the formation induced by the uniform stress, and internal pressure was decomposed into three parts as shown in Figure 7. They were the original stress field of A1, the excavation disturbance field of B1 induced by drilling of a wellbore and the interface disturbance field of C1 induced by the fluid column pressure.Figure 7
Wellbore stress components under the condition of uniform stress.In the polar coordinate system, the initial conditions ofA1 were σr30=σθ30=−p0, εr′30=p0/2G31−2μ3, and boundary stress conditions of B1 and C1 were σr′3r=R3,B1=p0 and σr′3r=R3,C1=−p3′. Substituting the initial and boundary conditions into equations (7) and (8), the displacement and stress in formation were obtained as shown below. Because R4 approached to infinity, A3′=0 and C3′=R32 were obtained:(9)ur′3r=12G3R32rp3′−p0,σr′3r=−p0−R32r2p3′−p0,σθ′3r=−p0+R32r2p3′−p0,where ur′3 was the radial displacement in formation and ur′3 and σθ′3 were the radial and tangential stresses in formation, respectively.


## 3.1.2. Casing and Cement Sheath Stress

The pressures at casing-cement sheath interface and cement sheath-formation interface werep2′ and p3′, respectively, (Figure 8). The initial stresses of the casing and the cement sheath were σr′i0=σθ′i0=0 and εr′i0=0. The boundary stress conditions were σr′2r=R3=−p3′ and σr′2r=R2=−p2′, σr′1r=R1=−pi.Figure 8
Interface pressures induced by uniform stress.Substituting these initial and boundary conditions in equations (7) and (8), the displacement and stress in casing and cement sheath were obtained as follows. Subscripts 1 and 2 represent the casing and cement sheath, respectively:(10)ur′1=12G11−2μ1A1′r+C1′1rpi+12G11−2μ1B1′r+C1′1r−p2′,ur′2=12G21−2μ2A2′r+C2′1rp2′+12G21−2μ2B2′r+C2′1r−p3′,(11)σr′1=A1′−C1′1r2pi+B1′−C1′1r2−p2′,σθ′1=A1′+C1′1r2p2′+B1′+C1′1r2−p2′,σr′2=A2′−C2′1r2p2′+B2′−C2′1r2−p3′,σθ′2=A2′+C2′1r2p2′+B2′+C2′1r2−p3′.According to the hypotheses that cement sheath-formation interface and casing-cement sheath interface were perfectly bonded to each other, the interfacial displacement continuity conditions were expressed in the following equation:(12)ur′1r=R2=ur′2r=R2,ur′2r=R3=ur′3r=R3.Substituting equation (10) into equation (12), the binary equations were obtained as(13)AAp2′−BBp3′=CCpi,DDp2′−EEp3′=FFp0,where(14)AA=12G11−2μ1B1′R2+C11R2+12G21−2μ2A2′R2+C2′1R2,BB=12G11−2μ1B1′R2+C11R2,CC=12G11−2μ1A1′r+C1′1R2,DD=12G21−2μ2A2′R3+C2′1R3,EE=12G11−2μ1B1′R+3C11R3+12G3R3,FF=−R32G3.The interfacial pressuresp2′ and p3′ could be calculated by using equation (14). Substituting them into equation (11), the stresses induced by uniform stress were obtained subsequently.


## 3.2. Stress Induced by Deviator Stress

The deviator stress boundary conditions are shown in Figure9. To calculate the stress distribution induced by deviator stress, the stress function was defined as(15)ϕ=Ai″r4+Bi″r2+Ci″+Di″r2cos2θ.Figure 9
Stress induced by deviator stress. Outer stress,scos2θ; interface pressure, pcos2θ.The stress and strain under the condition of nonuniform stress are(16)σr″i=−2Bi″+4Ci″r2+6Di″r4cos2θ,σθ″i=12Ai″r2+2Bi″+6Di″r4cos2θ,τrθ″i=6Ai″r2+2Bi″−2Ci″r2−6Di″r4sin2θ,εr″i=−1+μiEi12υiAi″r2+2Bi″+1−μi4Ci″r2+6Di″r4cos2θ−εr″i0,εθ″i=1+μiEi121−μiAi″r2+2Bi″+μi4Ci″r2+6Di″r4cos2θ−εθ″i0,where σr″iandσθ″i are the radial and tangential stresses; εr″i0andεθ″i0 are the initial radial and tangential strains; and Ai″,Bi″,Ci″,andDi″ were the constants, i = 1, 2, 3 represented the casing, cement sheath, and formation.From the geometric equations,(17)εr″i=∂ur″i∂r,εθ″i=1r∂uθ″i∂θ+ur″ir.The radial displacementur″i and tangential displacement uθ″i were obtained as(18)ur″i=−1+μiEi4μiAi″r3+2Bi″r−1−μi4Ci″r−2Di″r3cos2θ−rεr″i0,uθ″i=1+μi2Ei43−2μiAi″r3+4Bi″r−1−2μi4Ci″r+4Di″r3sin2θ+r∫εr″i0−εθ″i0dθ,εr″i0=1+μiEi1−μiσr″i0−μiσθ″i0,εθ″i0=1+μiEi1−μiσθ″i0−μiσr″i0,where σr″i0andσθ″i0 were the initial radial and tangential stresses.


### 3.2.1. Formation Stress

Similar to that of uniform stress, the actual stress field,F2 in the strata induced by the nonuniform stress, was decomposed into three parts: the original stress field, A2; the disturbance field, B2 induced by the wellbore excavation; and the interface pressure, C2 induced by the interface pressure (Figure 10).Figure 10
Formation stress components under the nonuniform stress condition.In the polar coordinate system, initial stresses wereσr″30=−scos2θ, σθ″30=scos2θ, and τrθ″30=ssin2θ; initial strains were εr″30=−1+μ3/E3⋅scos2θ and εθ″30=1+μ3/E3⋅scos2θ; and the boundary stresses were σr″3r=∞=−scos2θ, σθ″3=scos2θ, and τrθ″3r=∞=ssin2θ. Substituting the initial and boundary conditions into (14) and (15), it was obtained that A3″=0andB3″=S/2. The displacements and stresses in formation were expressed as shown in the following equations:(19)ur″3=−1G3−1−μ32C3″r−D3″r3cos2θ,uθ″3=1G3−1−2μ3C3″r+D3″r3sin2θ,(20)σr″3=−s+4C3″r2+6D3″r4cos2θ,σθ″3=s+6D3″r4cos2θ,τrθ″3=s−2C3″r2−6D3″r4sin2θ.


### 3.2.2. Casing and Cement Sheath Stress

For casing and cement sheath in the polar coordinate system, initial stresses wereσr″i0=σθ″i0=0 and initial strains were εr″i0=εθ″i0=0. Substituting the initial and boundary conditions into equations (14) and (15), the displacements and stresses were obtained as follows:(21)ur″i=−1Gi2μiAi″r3+Bi″r−1−μi2Ci″r−Di″r3cos2θ,uθ″i=1Gi3−2μiAi″r3+Bi″r−1−2μiCi″r+Di″r3sin2θ,(22)σr″i=−2Bi″+4Ci″r2+6Di″r4cos2θ,σθ″i=12Ai″r2+2Bi″+6Di″r4cos2θ,τrθ″i=6Ai″r2+2Bi″−2Ci″r2−6Di″r4sin2θ.The interfacial displacement and stress continuity and boundary conditions were expressed in the following equation:(23)σr″1r=R1=0,τrθ″1r=R1=0,σr″1r=R2=σr″2r=R2,τrθ″1r=R2=τrθ″2r=R2,σr″2r=R3=σr″3r=R3,τrθ″2r=R3=τrθ″3r=R3,ur″1r=R2=ur″2r=R2,uθ″1r=R2=uθ″2r=R2,ur″2r=R3=ur″3r=R3,uθ″2r=R3=uθ″3r=R3.Substituting equations (20)–(23) into the following equation, equations were obtained as(24)−2B1″−4C1″R12−6D1″R14cos2θ=0,6R12A1″+2B1″−2R12C1″−6R14D1″sin2θ=0,(25)−2B1″+4C1″R22+6D1″R24+2B2″+4C2″R22+6D2″R24=0,6A1″R22+2B1″−2C1″R22−6D1″R24−6A2″R22+2B2″−2C2″R22−6D2″R24=0,(26)−B2″+2C2″R32+3D2″R34+2C3″R32+3D3″R34=−s,3A2″R32+B2″−C2″R32−3D2″R34+C3″R32+3D3″R34=s,(27)−1G12μ1R23A1″+R2B1″−1−μ12R2C1″−1R23D1″+1G22μ2R23A2″+R2B2″−1−μ22R2C2″−1R23D2″=0,1G13−2μ1R23A1″+R2B1″−1−2μ11R2C1″+1R23D1″−1G23−2μ2R23A2″+R2B2″−1−2μ21R2C2″+1R23D2″=0,(28)−1G22μ2R33A2″+R3B2″−1−μ22R3C2″−1R33D2″+1G3−1−μ32R3C3″−1R33D3″=0,1G23−2μ2R33A2″+R3B2″−1−2μ21R3C2″+1R33D2″−1G3−1−2μ31R3C3″+1R33D3″=0,where the constants of A1″,B1″,C1″,D1″,A2″,B2″,C2″,D2″,C3″,andD3″ were calculated by the total 10 equations in equations (24)–(28). Then, wellbore stress distribution induced by deviator stress was obtained by substituting these 10 constants and A3″andB3″ into equations (19)–(22).


## 3.2.1. Formation Stress

Similar to that of uniform stress, the actual stress field,F2 in the strata induced by the nonuniform stress, was decomposed into three parts: the original stress field, A2; the disturbance field, B2 induced by the wellbore excavation; and the interface pressure, C2 induced by the interface pressure (Figure 10).Figure 10
Formation stress components under the nonuniform stress condition.In the polar coordinate system, initial stresses wereσr″30=−scos2θ, σθ″30=scos2θ, and τrθ″30=ssin2θ; initial strains were εr″30=−1+μ3/E3⋅scos2θ and εθ″30=1+μ3/E3⋅scos2θ; and the boundary stresses were σr″3r=∞=−scos2θ, σθ″3=scos2θ, and τrθ″3r=∞=ssin2θ. Substituting the initial and boundary conditions into (14) and (15), it was obtained that A3″=0andB3″=S/2. The displacements and stresses in formation were expressed as shown in the following equations:(19)ur″3=−1G3−1−μ32C3″r−D3″r3cos2θ,uθ″3=1G3−1−2μ3C3″r+D3″r3sin2θ,(20)σr″3=−s+4C3″r2+6D3″r4cos2θ,σθ″3=s+6D3″r4cos2θ,τrθ″3=s−2C3″r2−6D3″r4sin2θ.


## 3.2.2. Casing and Cement Sheath Stress

For casing and cement sheath in the polar coordinate system, initial stresses wereσr″i0=σθ″i0=0 and initial strains were εr″i0=εθ″i0=0. Substituting the initial and boundary conditions into equations (14) and (15), the displacements and stresses were obtained as follows:(21)ur″i=−1Gi2μiAi″r3+Bi″r−1−μi2Ci″r−Di″r3cos2θ,uθ″i=1Gi3−2μiAi″r3+Bi″r−1−2μiCi″r+Di″r3sin2θ,(22)σr″i=−2Bi″+4Ci″r2+6Di″r4cos2θ,σθ″i=12Ai″r2+2Bi″+6Di″r4cos2θ,τrθ″i=6Ai″r2+2Bi″−2Ci″r2−6Di″r4sin2θ.The interfacial displacement and stress continuity and boundary conditions were expressed in the following equation:(23)σr″1r=R1=0,τrθ″1r=R1=0,σr″1r=R2=σr″2r=R2,τrθ″1r=R2=τrθ″2r=R2,σr″2r=R3=σr″3r=R3,τrθ″2r=R3=τrθ″3r=R3,ur″1r=R2=ur″2r=R2,uθ″1r=R2=uθ″2r=R2,ur″2r=R3=ur″3r=R3,uθ″2r=R3=uθ″3r=R3.Substituting equations (20)–(23) into the following equation, equations were obtained as(24)−2B1″−4C1″R12−6D1″R14cos2θ=0,6R12A1″+2B1″−2R12C1″−6R14D1″sin2θ=0,(25)−2B1″+4C1″R22+6D1″R24+2B2″+4C2″R22+6D2″R24=0,6A1″R22+2B1″−2C1″R22−6D1″R24−6A2″R22+2B2″−2C2″R22−6D2″R24=0,(26)−B2″+2C2″R32+3D2″R34+2C3″R32+3D3″R34=−s,3A2″R32+B2″−C2″R32−3D2″R34+C3″R32+3D3″R34=s,(27)−1G12μ1R23A1″+R2B1″−1−μ12R2C1″−1R23D1″+1G22μ2R23A2″+R2B2″−1−μ22R2C2″−1R23D2″=0,1G13−2μ1R23A1″+R2B1″−1−2μ11R2C1″+1R23D1″−1G23−2μ2R23A2″+R2B2″−1−2μ21R2C2″+1R23D2″=0,(28)−1G22μ2R33A2″+R3B2″−1−μ22R3C2″−1R33D2″+1G3−1−μ32R3C3″−1R33D3″=0,1G23−2μ2R33A2″+R3B2″−1−2μ21R3C2″+1R33D2″−1G3−1−2μ31R3C3″+1R33D3″=0,where the constants of A1″,B1″,C1″,D1″,A2″,B2″,C2″,D2″,C3″,andD3″ were calculated by the total 10 equations in equations (24)–(28). Then, wellbore stress distribution induced by deviator stress was obtained by substituting these 10 constants and A3″andB3″ into equations (19)–(22).


## 3.3. Stress Induced by Shear Stress

The stress induced by shear stress wasur″′i,uθ″′i,σr″′i,σθ″′i,andτrθ″′i, i = 1, 2, 3 represented the casing, cement sheath, and formation, respectively, (Figure 11). The angle of Ω between σx and x-direction was calculated by using equation (29). Then, the principal stresses were obtained as follows [29]:(29)Ω=12arctan−2τxyσx−σy=π4,(30)σπ/4=τxy,σ−π/4=−τxy.Figure 11
Stress distribution induced by shear stress.It could be found that the stress distribution induced by shear stress was similar with that by deviator stress when counterclockwise rotating the angle ofπ/4. The stresses and displacements were obtained by substituting the reference variable θ=θ′−π/4 into the stress induced by deviator stress discussed in Section 3.2.


## 3.4. Stress Induced by Temperature Variation

The thermal field was obtained by using the steady temperature distribution model to calculate the thermal stress. When fracturing fluids were pumped into a wellbore with a high pump rate, they were always in the turbulent state. The heat transfer coefficient between casing and fluid was calculated using the Marshall model [30] shown in the following equation:(31)h=StkmD=0.0107kmDρaDeff4Q/πD2K3n+1/4nn32Q/πD3n−10.67K3n+1/4nn32Q/πD3n−1Cmkm0.33,where h is the heat transfer coefficient (W·m−2·°C−1), St is the Stanton number, Pr is the Prandtl number, Reg is the Reynolds number, μw,app is the fluid apparent viscosity, D is the inner diameter (m), Deff is the equivalent diameter (m), ρa is the fluid density (kg·m−3), n is the liquidity index, K is the consistency coefficient (Pa·sn), v is the fluid velocity, Q is the fracturing pump rate (m3·min−1), km is the coefficient of heat conductivity (W·m−1·°C−1), and Cm is the fluid specific heat capacity (J·kg−1·°C−1).The temperature distribution among casing, cement sheath, and formation is shown in Figure12. In the cylindrical coordinate system of CCF, the differential equation of steady heat conduction of the cylinder is expressed as [31](32)d2Tdr2+1rdTdr=0.Figure 12
The distribution of interface temperatures.Temperature field distribution solutions were obtained according to integral and boundary conditionskidT/dr=hnTi−Tn,Tr=Ri=Ti,Tr=Ri+1=Ti+1 shown in the following equation:(33)Tir=AiTlnr+BiT,(34)AiT=Ti+1−TnlnRi+1/Ri+i2/2−5i/2+3ki/Rih,BiT=TnlnRi+1−Ti+1lnRi+i2/2−5i/2+3ki/RihTi+1lnRi+1/Ri+i2/2−5i/2+3ki/Rih,where Ti is the temperature (°C); Ti is the temperature at the interface (°C); Tn is the fluid temperature (°C); ki is the material thermal conductivity (W·m−1·°C−1); Ri is the radius (m); and AiTandBiT were the constants, i = 1, 2, 3 represented casing, cement sheath, and formation, respectively.The heat flow density continuity conditions were expressed as(35)kidTirdrr=Ri+1=ki+1dTi+1rdrr=Ri+1.The temperatures at interfaces of casing-cement sheath and cement sheath-formation system were defined asT2 and T3 and were calculated by using the following equation:(36)1+β1T2−β1T3=T1,−T2+1+β2T3=β2T4,where(37)β1=k2k1lnR2/R1+k1/R1hlnR3/R2,β2=k3k2lnR3/R2lnR4/R3.Interfacial temperature ofTi was obtained by solving equation (36). The steady-state temperature field around the wellbore could be calculated by substituting Ti into equations (33) and (34). According to thermal elastic mechanics, constitutive equations for a plane strain problem were expressed as(38)εrT=1+μiEi1−μiσrT−μiσθT+1+μiαiT,εθT=1+μiEi1−μiσθT−μiσrT+1+μiαiT,εzT=0,γrθT=21+μiEiτrθT.The actual thermal stress field,F3 in the strata induced by the temperature changes, was decomposed into two parts: the original stress field, A3, and the disturbance field, B3 induced by the temperature variation shown in Figure 13.Figure 13
Thermal stress field. (a) Formation stress components. (b) Interface pressures.pi is the interface pressure, i = 2, 3 represented the casing-cement sheath interface and cement sheath-formation interface.
(a)
(b)The initial stresses wereσrTi0=σθTi0=0, and the initial strains were εrTi0=εθTi0=0. The stresses and displacements induced by thermal variations were expressed as(39)urTi=1+μi1−μiαir∫RirrΔTidr+C1Tir+C2Tir−rεrTi0,(40)σrTi=−αiEi1−μi1r2∫RirrΔTidr+Ei1+μiC1Ti1−2μi−C2Tir2,σθTi=αiEi1−μi1r2∫RirrΔTidr+Ei1+μiC1Ti1−2μi+C2Tir2−αiEi1−μiΔTi,where C1TiandC2Ti are the constants, σrTiandσθTi are the radial and tangential stresses (Pa), urTi is the radial displacement (m), ΔTi is the temperature changes (°C), pi is the interface pressure (Pa), and αi is the material thermal expansion coefficient, i = 1, 2, 3 represented casing, cement sheath, and formation, respectively.The temperatures were known, and the boundary was free at internal casing and external formation. So, radial stress at inner and outer boundaries equals to zero, and radial displacement at the outer boundary equals to zero as well. The boundary and interfacial displacement continuity conditions were expressed as(41)urT1r=R2=urT2r=R2,urT2r=R3=urT3r=R3,σrT1r=R2=σrT2r=R2,σrT2r=R3=σrT3r=R3,σrT1r=R1=0,σrT3r=R4=0.Substituting equations (39) and (40) into the following equation, the equations were obtained as(42)C1T1R2+C2T1R2−C1T2R2−C2T2R2=−1+μ11−μ1α1R2∫R1R2rT1dr,C1T2R3+C2T2R3−C1T3R3−C2T3R3=1+μ21−μ2α2R3∫R2R3rT2dr,(43)E11+μ1C1T11−2μ1−C2T1R22−E21+μ2C1T21−2μ2−C2T2R22=α1E11−μ11R22∫R1R2rΔT1dr,E21+μ2C1T21−2μ2−C2T2R32−E31+μ3C1T31−2μ3−C2T3R32=α2E21−μ21R32∫R2R3rΔT2dr,(44)C11R1+C21R1=0,C13R4+C23R4=−1+μ21−μ2α2R4∫R3R4rTdr.The constants ofC1T1,C2T1,C1T2,C2T2,C1T3,andC2T3 were obtained by equations (42)–(44). The wellbore stress was obtained by substituting these constants into equation (40).The total stresses were obtained using the following equation:(45)σri=σr′i+σr″i+σr″′i+σrTi,σθi=σθ′i+σθ″i+σθ″′i+σθTi,σzi=μiσri+σθi,τrθi=τrθ″i+τrθ″′i,where σri is the radial stress, σθi is the tangential stress, σzi is the axial stress, and τrθi is the shear stress.


## 3.5. Estimation of Wellbore Integrity

It is generally accepted that the yield of isotropic material such as casing has nothing to do with hydrostatic pressure, while hydrostatic pressure is not considered in von Mises yield criterion. So, this criterion was adopted to determine the casing failure:(46)fJ2,k=J2−k=0,J2=16σ11−σ222+σ22−σ332+σ33−σ112+σ122+σ232+σ312,where J2 is the second stress partial tensor; k is the critical value of failure; and σij is the stress components, i, j = 1, 2, 3 represented the three directions of the system, respectively.For uniaxial tension,J2=σ/3, the von Mises stress could be expressed as follows in the polar coordinate:(47)σMises=12σr−σθ2+σθ−σz2+σz−σr2+3τrθ2+3τθz2+3τzr2.


## 4. Model Validation

From 2009 to 2017, PetroChina has drilled 141 fracturing wells (including 112 horizontal wells) in the Changning-Weiyuan National Shale Gas Demonstration Area. The geometrical dimensions of the CCF model were a wellbore diameter of 8.5 in, casing diameter of 5.5 in, and casing thickness of 9.17 mm. According to the Saint-Venant principle, a formation boundary dimension should be five to six times larger than that of the wellbore geometry to avoid the influence of boundary effect on wellbore stress. In view of this, the model geometry was 2,000 × 2,000 mm, while the corresponding wellbore diameter was 215.9 mm. The direction of horizontal in situ stress was N120°E. The well deviation angle was 90°, and the wellbore azimuth was N30°E, indicating that the horizontal trajectory was along the minimum in situ stress direction. The internal casing pressure was calculated from the pump pressure plus the downhole hydrostatic fluid pressure. The external boundary stress was obtained from the geostress data of the shale reservoir. The thermal and mechanical properties of different materials are presented in Table1. The casing stress and displacement were calculated and analyzed considering thermal-pressure coupling.Table 1
Thermal and mechanical Parameters of fluid-casing-cement sheath-formation system.
Property
Casing
Cement sheath
Formation
Fluid
Elastic modulus,Ei (GPa)
210
5
35
—
Poisson’s ratio,μi
0.3
0.15
0.25
—
Coefficient of thermal expansion,αi (10−5·°C−1)
1.5
1.0
1.0
—
Thermal conductivity,ki (W·m−1·°C−1)
58.2
1.0
1.0
1.73
Specific heat,Cpi (J·kg−1·°C−1)
460
1830
1043
3935
Density,ρi (kg·m−3)
7850
1800
2500
1080
Note: properties in parenthesis were used in the parametric study.The applied maximum horizontal stressσH was 82 MPa, the minimum horizontal stress σh was 55 MPa, the vertical stress σv was 57 MPa, the inner casing pressure Pi was 75 MPa, the boundary temperature T4 was 100°C, the fluid temperature Ta was 20°C, and the convective heat transfer coefficient was obtained by using equation (20) (1 890 W·m−2·°C−1) at the pump rate of 20 m3/min.The finite element analysis method was adopted to validate the results of the analytical models. A steady-state thermal analysis followed by a static structural analysis was conducted to calculate the stress considering thermal-pressure coupling. The solutions of radial stress, circumferential stress, and Mises stress are compared in Figure14.Figure 14
Comparison of numerical and analytical solutions. (a) Radial stress along the radial directions of 0° and 90°. (b) Radial stress at the internal casing wall. (c) Circumferential stress along the radial directions of 0° and 90°. (d) Circumferential stress at the internal casing wall. (e) Mises stress at inner and outer casing walls.
(a)
(b)
(c)
(d)
(e)The analytical solutions of radial stress, circumferential stress, and Mises stress were in good agreement with the results obtained by a finite element method, which indicates the validity of the analytical method. The maximum deviation between analytical and finite element results was 1.4–13.9%, indicating that the analytical model could provide an accurate calculation of stress distribution for the CCF system.From Figures14(a) and 14(b), the radial stress increased with the increase of radius in casing and cement sheath, but decreased in the formation. The absolute value of radial stress calculated by the new model was smaller than that of the existing model. This was mainly because the new model excluded the strain induced by the initial stress. From Figures 14(c) and 14(d), the circumferential stress decreased with the increase of radius in the casing and cement sheath and increased slowly to a constant value in the formation. The interfacial stress at the internal casing wall was larger than that at the external casing wall. The solutions calculated by the new model were larger than those by the existing model. From Figure 14(e), casing Mises stress obtained by the new model was larger than that of the existing model. It could be explained that circumferential stress was larger than radial stress and had a main influence on Mises stress.The radial displacements along the 0° direction calculated by the new model and existing model under the same conditions were shown in Figure15. There was an obvious difference for two models, especially at the outer boundary. The displacements of new model approached zero when the outer boundary was infinite, which reached an agreement with the actual boundary condition. However, the displacements obtained by the existing model increased linearly in the formation. So, only the new model could reflect the actual situation.Figure 15
Radial displacements of the wellbore assembly along the 0° direction.


## 5. Sensitivity Analysis

The sensitivity analyses were carried out to study the influences of cement sheath properties, geostress, fracturing pressure, fluid temperature, casing thickness, and cement sheath thickness on casing stress. During analyzing, only one parameter was variable and others were constants. Unless otherwise mentioned, the parameters were set as mentioned in Section4.


### 5.1. Influence of Elastic Modulus

Cement sheath properties is crucial for casing safety. To evaluate the effect of elastic modulus on casing stress, the cement sheath elastic modulus ofE2 was set at the range from 2 GPa to 50 GPa, and the formation elastic modulus of E3 was set as 5 and 35 GPa to simulate a soft and hard formation. The Mises stresses at internal casing are shown in Figure 16.Figure 16
Casing Mises stress. (a, b)E3 = 5 GPa and (c, d) E3 = 35 GPa.
(a)
(b)
(c)
(d)From Figures16(a) and 16(c), the maximum Mises stress appeared at the angles of 0° and 180° for the new model and 90° and 270° for the existing model when the formation modulus was small. However, the maximum stress all appeared at the angles of 0° and 180° for the new and existing models when the formation modulus was large. From Figure 16(b), in a soft formation (a modulus of 5 GPa), with the increase of the cement sheath modulus, the maximum casing stress increased first, and then decreased for existing model, while decreasing all the time for the new model. From Figure 16(b), in a hard formation (modulus of 35 GPa), the maximum casing stress always decreased with the increase of the cement sheath modulus for two models. In the soft formation, the stress calculated by the new model was smaller than that by the existing model. However, the stress obtained by the new model was larger than that by the existing model in a hard formation. According to the fact that shale formation had a large elastic modulus, the existing model underestimated casing stress during the fracturing operation.


### 5.2. Influence of Poisson’s Ratio

To evaluate the effect of Poisson’s ratio on casing stress, cement sheath Poisson’s ratio,μ2, with a range from 0.05 to 0.45, was adopted and the formation Poisson’s ratio, μ3, was set as 0.05 and 0.45 to simulate a hard and soft formation. The casing Mises stresses are shown in Figure 17.Figure 17
Casing Mises stresses for different Poisson’s ratios. (a)μ3 = 0.05. (b) μ3 = 0.45.
(a)
(b)From Figures17(a) and 17(b), the maximum Mises stress decreased with the increase of cement sheath Poisson’s ratio for two models. In a hard formation (Poisson’s ratio of 0.05), the maximum stress obtained by the new model was larger than that by the existing model. However, in a soft formation (Poisson’s ratio of 0.45), it was a little smaller than that by the existing model. According to the fact that shale formation had a small Poisson’s ratio, the existing model underestimated casing stress during the fracturing process.


### 5.3. Influence of In Situ Stress Nonuniformity

During the multifracturing operation for shale gas wells, the fracturing fluid was pressed into the formation and the in situ stress field changed abruptly to increase the nonuniformity of the stress around the wellbore. To evaluate the effect of in situ stress nonuniformity on casing stress, the nonuniformity index was defined asδ=σH/σv. Different δ with a range of 0.1–3.0 was adopted. The casing Mises stresses calculated by two models are shown in Figure 18.Figure 18
Casing Mises stress for different stress nonuniformity indexes: (a) casing internal stress and (b) maximum casing stress.
(a)
(b)As seen from Figure18(a), for δ smaller than 1.0, the maximum Mises stresses appeared at 90° and 270° directions. However, for δ larger than 1.0, the maximum Mises stresses appeared at 0° and 180° directions. For δ of 1.0, the casing Mises stress around the wellbore was at a uniform state. From Figure 18(b), the maximum casing stress increased dramatically with the increase of δ from 1.0 or decrease of δ from 1.0. The solution obtained by the new model was larger than that by the existing model for δ between 0.5 and 2.25. When δ was larger than 2.25 or smaller than 0.5, the casing stress obtained by the existing model was larger than that by the new model.


### 5.4. Influence of Fracturing Pressure

A fracturing fluid with high pressure was used to fracture a shale formation. The high pressure depended on the formation regional tectonic stress; the larger the tectonic stress, the higher the pressure. Moreover, a high fracturing pressure posed a great potential challenge to casing failure. Different fracturing pressures with a range of 5–105 MPa were adopted to evaluate the effect of fracturing pressure on casing stress. The maximum casing Mises stresses are shown in Figure19.Figure 19
Maximum casing Mises stress for different pressures.It can be seen from Figure19 that the casing stresses calculated by the two models decreased first and then increased with the increase of fracturing pressure. The minimum stress appeared at 15 MPa for the new model; however, it appeared at about 25 MPa for the existing model. In addition, the casing Mises stress obtained by the existing model was larger than that by the new model for pressure lower than 25 MPa and smaller than that by the new model for pressure higher than 25 MPa. During fracturing operation, pressure must be large enough to fracture the formation, so the existing model underestimated the casing stress.


### 5.5. Influence of Fluid Temperature

During the cycle injection of fracturing fluid, the heat transfer coefficienth was calculated using equation (20) with a pump rate of 20 m3/min. The corresponding casing internal Mises stress was calculated under different fluid temperatures at a range of 10–100°C to evaluate the effect on casing stress. Figure 20 presented the maximum casing stress over temperature and the comparison of the results obtained by the existing model and new model.Figure 20
Casing Mises stress for different temperatures.From Figure20, the maximum Mises stress decreased with the increase of the injection fluid temperature, indicating that a fracturing fluid with high temperature was effective to decrease casing stress. Furthermore, the stress obtained by the existing model was smaller than that by the new model. It revealed that the existing model underestimated the casing Mises stress.


### 5.6. Influence of Thickness

The thickness of cement sheath and casing was curial for casing safety. To evaluate the effect of thickness on the casing stress, the cement thickness was set at a range of 2–50 mm and the casing thickness was set at a range of 5–15 mm. The comparisons of maximum casing Mises stress obtained by the two models are shown in Figure21.Figure 21
Casing Mises stress for different thicknesses of casing and cement sheath. (a) Different cement thickness. (b) Different casing thickness.
(a)
(b)As shown in Figure21, the maximum casing Mises stress increased with the increase of cement sheath thickness and, however, decreased with the increase of casing thickness. So, a thicker casing wall and thinner cement sheath were effective to ensure the casing integrity. Meanwhile, casing stress obtained by the existing model was smaller than that by the new model indicating that the existing model underestimated casing stress.


## 5.1. Influence of Elastic Modulus

Cement sheath properties is crucial for casing safety. To evaluate the effect of elastic modulus on casing stress, the cement sheath elastic modulus ofE2 was set at the range from 2 GPa to 50 GPa, and the formation elastic modulus of E3 was set as 5 and 35 GPa to simulate a soft and hard formation. The Mises stresses at internal casing are shown in Figure 16.Figure 16
Casing Mises stress. (a, b)E3 = 5 GPa and (c, d) E3 = 35 GPa.
(a)
(b)
(c)
(d)From Figures16(a) and 16(c), the maximum Mises stress appeared at the angles of 0° and 180° for the new model and 90° and 270° for the existing model when the formation modulus was small. However, the maximum stress all appeared at the angles of 0° and 180° for the new and existing models when the formation modulus was large. From Figure 16(b), in a soft formation (a modulus of 5 GPa), with the increase of the cement sheath modulus, the maximum casing stress increased first, and then decreased for existing model, while decreasing all the time for the new model. From Figure 16(b), in a hard formation (modulus of 35 GPa), the maximum casing stress always decreased with the increase of the cement sheath modulus for two models. In the soft formation, the stress calculated by the new model was smaller than that by the existing model. However, the stress obtained by the new model was larger than that by the existing model in a hard formation. According to the fact that shale formation had a large elastic modulus, the existing model underestimated casing stress during the fracturing operation.


## 5.2. Influence of Poisson’s Ratio

To evaluate the effect of Poisson’s ratio on casing stress, cement sheath Poisson’s ratio,μ2, with a range from 0.05 to 0.45, was adopted and the formation Poisson’s ratio, μ3, was set as 0.05 and 0.45 to simulate a hard and soft formation. The casing Mises stresses are shown in Figure 17.Figure 17
Casing Mises stresses for different Poisson’s ratios. (a)μ3 = 0.05. (b) μ3 = 0.45.
(a)
(b)From Figures17(a) and 17(b), the maximum Mises stress decreased with the increase of cement sheath Poisson’s ratio for two models. In a hard formation (Poisson’s ratio of 0.05), the maximum stress obtained by the new model was larger than that by the existing model. However, in a soft formation (Poisson’s ratio of 0.45), it was a little smaller than that by the existing model. According to the fact that shale formation had a small Poisson’s ratio, the existing model underestimated casing stress during the fracturing process.


## 5.3. Influence of In Situ Stress Nonuniformity

During the multifracturing operation for shale gas wells, the fracturing fluid was pressed into the formation and the in situ stress field changed abruptly to increase the nonuniformity of the stress around the wellbore. To evaluate the effect of in situ stress nonuniformity on casing stress, the nonuniformity index was defined asδ=σH/σv. Different δ with a range of 0.1–3.0 was adopted. The casing Mises stresses calculated by two models are shown in Figure 18.Figure 18
Casing Mises stress for different stress nonuniformity indexes: (a) casing internal stress and (b) maximum casing stress.
(a)
(b)As seen from Figure18(a), for δ smaller than 1.0, the maximum Mises stresses appeared at 90° and 270° directions. However, for δ larger than 1.0, the maximum Mises stresses appeared at 0° and 180° directions. For δ of 1.0, the casing Mises stress around the wellbore was at a uniform state. From Figure 18(b), the maximum casing stress increased dramatically with the increase of δ from 1.0 or decrease of δ from 1.0. The solution obtained by the new model was larger than that by the existing model for δ between 0.5 and 2.25. When δ was larger than 2.25 or smaller than 0.5, the casing stress obtained by the existing model was larger than that by the new model.


## 5.4. Influence of Fracturing Pressure

A fracturing fluid with high pressure was used to fracture a shale formation. The high pressure depended on the formation regional tectonic stress; the larger the tectonic stress, the higher the pressure. Moreover, a high fracturing pressure posed a great potential challenge to casing failure. Different fracturing pressures with a range of 5–105 MPa were adopted to evaluate the effect of fracturing pressure on casing stress. The maximum casing Mises stresses are shown in Figure19.Figure 19
Maximum casing Mises stress for different pressures.It can be seen from Figure19 that the casing stresses calculated by the two models decreased first and then increased with the increase of fracturing pressure. The minimum stress appeared at 15 MPa for the new model; however, it appeared at about 25 MPa for the existing model. In addition, the casing Mises stress obtained by the existing model was larger than that by the new model for pressure lower than 25 MPa and smaller than that by the new model for pressure higher than 25 MPa. During fracturing operation, pressure must be large enough to fracture the formation, so the existing model underestimated the casing stress.


## 5.5. Influence of Fluid Temperature

During the cycle injection of fracturing fluid, the heat transfer coefficienth was calculated using equation (20) with a pump rate of 20 m3/min. The corresponding casing internal Mises stress was calculated under different fluid temperatures at a range of 10–100°C to evaluate the effect on casing stress. Figure 20 presented the maximum casing stress over temperature and the comparison of the results obtained by the existing model and new model.Figure 20
Casing Mises stress for different temperatures.From Figure20, the maximum Mises stress decreased with the increase of the injection fluid temperature, indicating that a fracturing fluid with high temperature was effective to decrease casing stress. Furthermore, the stress obtained by the existing model was smaller than that by the new model. It revealed that the existing model underestimated the casing Mises stress.


## 5.6. Influence of Thickness

The thickness of cement sheath and casing was curial for casing safety. To evaluate the effect of thickness on the casing stress, the cement thickness was set at a range of 2–50 mm and the casing thickness was set at a range of 5–15 mm. The comparisons of maximum casing Mises stress obtained by the two models are shown in Figure21.Figure 21
Casing Mises stress for different thicknesses of casing and cement sheath. (a) Different cement thickness. (b) Different casing thickness.
(a)
(b)As shown in Figure21, the maximum casing Mises stress increased with the increase of cement sheath thickness and, however, decreased with the increase of casing thickness. So, a thicker casing wall and thinner cement sheath were effective to ensure the casing integrity. Meanwhile, casing stress obtained by the existing model was smaller than that by the new model indicating that the existing model underestimated casing stress.


## 6. Conclusions

A new analytical model considering drilling construction was established to assess the casing stress under different conditions considering thermal-pressure coupling. The solutions were obtained by dividing the model into three parts such as initial stress field, wellbore disturbance field, and thermal stress field. Sensitivity analyses of different factors were conducted to evaluate the influences on casing stress. Some conclusions were drawn from the comparisons between new model and existing model:(1)
The results of radius stress, tangential stress, and casing Mises stress calculated by the analytical method were in good agreement with the solutions by a finite element solution. The minor deviations did not exceed 13.9%. Moreover, the analytical solutions were in-line with the actual boundary conditions of shale gas wells.(2)
The casing stress calculated by the existing model was smaller than that by the new model for hard formation with larger modulus or low Poisson’s ratio, geostress heterogeneity index at a range of 0.5–2.25, and fracturing pressure larger than 25 MPa.(3)
The casing stress increased with the increase of the in situ stress nonuniformity index. With the increase of fracturing pressure, casing stress decreased first and then increased.(4)
Cement sheath with appropriate modulus and larger Poisson’s ratio, thinner cement sheath, thicker casing, and higher fluid temperature were effective to decrease the casing stress.In conclusion, the new analytical model can accurately predict casing stress and become an alternative method of casing integrity evaluation for shale gas wells. It is a useful and efficient method for a preliminary design, being capable of simulating the actual situations in order to assess the casing stresses and integrity.



---

*Source: 4253241-2019-04-24.xml*