A Three-Dimensional Analytical Solution of Stress Field in Casing-Cement-Stratum System Considering Initial Stress State

Abstract

Accurate stress field calculation of the casing-cement-stratum system is crucial for evaluating wellbore integrity. Previous models treated in-situ stress as boundary pressure loads, leading to unrealistic infinite displacements at infinity. This study presents a three-dimensional (3D) analytical solution for the stress field within the casing-cement-stratum system in inclined wells, considering in-situ stress and hydrostatic stress in cement as the initial stress state and taking into account stress components related to the axial direction. Assuming a plane strain condition and superimposing the in-plane plane strain problem, elastic uni-axial stress problem and anti-plane shear problem, a 3D analytical solution is obtained. Comparisons with previous models indicate that the existing model overestimates the absolute values of stress components and failure potential of casing and cement in both 2D and 3D scenarios. The presence of initial stress in cement greatly increases the absolute value of the compressive stress state but decreases the failure potential in cement, which has not been well studied. Additionally, a low Young’s modulus and high initial stress state of the cement benefits the cement’s integrity since the maximum Mises stress significantly decreases. The new 3D analytical solution can provide a benchmark for 3D numerical simulation and quick assessment for wellbore integrity.

1. Introduction

Zonal isolation often fails throughout the life of oil/gas wells due to mechanical failure of casing and cement [1]. Moreover, the oil/gas leaks from the reservoir to the shallow aquifers and abandonment of the well occurs, and environmental and safety concerns are inevitable [2,3]. Hence, precise stress field calculations of the casing-cement-stratum system play a crucial role in evaluating mechanical wellbore integrity.

Aside from numerical modeling works, analytical solutions have been developed for determining the stress field in the casing-cement-stratum system [4,5,6]. Yin et al. [7] developed a theoretical solution for in-plane stress distribution in the casing-cement-stratum system. Moreover, Yin and Gao [8] derived the in-plane stress field in directional (inclined) cemented wells. Then, in the directional wells [9] and vertical section of the wells [10], the stress field and collapse resistance of multilayer cemented casing were examined. The wellbore integrity of HPHT gas wells [11], carbon sequestration projects [12], well testing and production process [13] have all been analytically analyzed by taking into consideration the thermal effect. Moreover, the impact of casing eccentricity on the failure of the cement sheath was analytically examined [14]. Xia et al. [15] obtained a transient response to the casing-cement-stratum system under dynamic radial tractions. Li et al. [16] obtained an analytical thermo-poroelastic response of the casing-cement-stratum system under internal hydraulic and thermal perturbation. However, the aforementioned models only obtained the in-plane stress field in vertical and inclined wells, neglecting the 3D stress field, which takes into account stress components with the axial direction. Atkinson and Eftaxiopoulos [17] employed complex potentials to determine the stress field in the in-plane and anti-plane issues for an inclined, cased and cemented wellbore. Vitali et al. [18] combined classical Kirsch and Einstein-Schwartz solutions to obtain an analytical solution for tunnels not aligned with geostatic principal stress directions, which provides the stress and displacement distribution induced by far-field axial shear stresses. Jo [19] derived the 3D stress field for inclined wells under the assumption of generalized plane strain condition and taking thermal effect into account.

However, the majority of models [7,8,9,10,11,12,13,14,15,16,17,18,19] considered the in-situ stress that already existed prior to excavation as the boundary pressure load rather than as the initial stress; this results in infinite displacements at infinity in the stratum and significant error for calculating an appropriate stress field [20]. Li et al. [21] considered the in-situ stress as the initial stress state to assess the casing integrity during the hydraulic fracturing of shale gas wells. Wei et al. [22] simulated the factors influencing the cement sheath integrity in hydraulic fractured wells, discovering that injection pressure significantly affects cement debonding. Yan et al. [23] simulated the debonding degree of cementing interfaces during fracturing and found that improving the interfacial bonding strength and reducing the elastic modulus of cement can promote wellbore integrity. Simone et al. [24] developed an analytical solution to calculate the in-plane stress field of the vertical casing-cement-stratum system with uniform in-situ stress, taking into account the drilling, construction, and production phases. To the best of the authors’ knowledge, a 3D stress field has not yet been established for the inclined casing-cement-stratum system that treats in-situ stress as the initial stress state.

A further factor that must be overlooked is the cement’s initial stress state. According to Saint-Marc et al. [25], attaining long-term cement-sheath integrity depends on the initial condition of stress. Petersen and Ulm [26] examined the impact of early-age stress and pressure developments in the setting of cement on cement failure. When Meng et al. [27] proposed a model to estimate the drop in total interface stress in cement during hydration, they discovered that increasing the cement’s initial stress will protect it. For the sake of simplification, the hydrostatic stress of cement is frequently taken to be equal to the initial stress state [28,29,30].

This work establishes a 3D analytical model of the stress field in the casing-cement-stratum system for inclined wells to assess wellbore integrity. The 3D problem is separated into three parts on the plane strain condition: in-plane plane strain problem, elastic uni-axial stress problem and anti-plane shear problem. Moreover, the in-situ stress in the stratum and the hydrostatic stress in the cement are regarded as the initial stress state. The stress field in the casing-cement-stratum system for inclined wells is then obtained by applying continuity of stress and displacement at the interfaces and boundary conditions. Comparisons with the existing model are made to demonstrate the necessity of treating the in-situ stress as the initial stress. In addition, a sensitivity analysis is conducted to estimate the influence of different factors on the integrity of the casing and cement. In the end, a 3D analytical solution can provide a benchmark for numerical simulation and a rapid assessment of wellbore integrity.

2. Formulation

2.1. Problem Description, Decomposition and Basic Equations

The system of steel casing, cement sheath, and stratum constitutes the near wellbore region after wellbore drilling and completion, as depicted in Figure 1. Assuming perfect bonding at the interfaces of casing-cement and cement-stratum, the continuity of the corresponding stress and displacement components is satisfied. The inner radii of casing, cement and stratum are a, b, and c, respectively. It should be noted that in-situ compressive stresses in the stratum have been applied before well construction, so they represent the initial stress state for the casing-cement-stratum system. Then, the three principle stresses in the stratum at the Cartesian coordinate can be decomposed by the vertical stress component ${\sigma }_v$ and unequal horizontal stress components ${\sigma }_{\mathrm{H}}\mathrm{and}{\sigma }_{\mathrm{h}}$, respectively. Moreover, the gravity of the casing and cement (hydrostatic pressure) has also been exerted before wellbore completion. In this study, the initial stress state in the cement is simulated by the hydraulic pressure of the slurry. Meanwhile, a uniform wellbore pressure $p_w$ acts on the inner surface of the casing. Given that the cross-sectional dimension is trivial compared to the length dimension of the system, this 3D problem can be treated as a plane strain problem [31,32].

The transformation from the Cartesian coordinates to the borehole coordinate is shown in Figure 2, where an inclined or horizontal wellbore in the petroleum industry is considered. We assume that the borehole is inclined to the X-Y-Z far-field principal stress axes. The Cartesian coordinate system coincides with the principal axes of the in-situ far-field stress, which are designated as ${\sigma }_x, {\sigma }_y, {\sigma }_z$. The borehole coordinate is defined by a local coordinate system x-y-z, with the borehole axes coinciding with the z-axis. In the borehole system, $\mathsf{\alpha }$ is an azimuth angle to the Z-axis, and $\mathsf{\beta }$ is a zenith angle toward the Z-axis.

Then, in the borehole coordinates, the stress tensor reads:

$$\sigma =\left[\begin{array}{ccc}{\sigma }_{xx}& {\sigma }_{xy}& {\sigma }_{xz}\\ {\sigma }_{xy}& {\sigma }_{yy}& {\sigma }_{zy}\\ {\sigma }_{xz}& {\sigma }_{zy}& {\sigma }_{zz}\end{array}\right]=\mathit{L}\left[\begin{array}{ccc}{\sigma }_H& 0& 0\\ 0& {\sigma }_h& 0\\ 0& 0& {\sigma }_v\end{array}\right]{\mathit{L}}^{\mathit{T}}$$

(1)

where

$$\mathit{L}=\left[\begin{array}{ccc}\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\beta }\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\alpha }& -\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\beta }\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\alpha }& \mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\beta }\\ \mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\alpha }& \mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\alpha }& 0\\ -\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\beta }\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\alpha }& \mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\beta }\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\alpha }& \mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\beta }\end{array}\right]$$

(2)

To obtain the principal horizontal stress components, a rotation of coordinates around the z-axis is performed. The angle of rotation is given as:

$${\theta }_r=\frac{1}{2}{\tan }^{-1}\frac{2{\sigma }_{xy}}{{\sigma }_{xx}-{\sigma }_{yy}}$$

(3)

.

Then, the initial in-situ stress in the polar coordinate reads:

$${\sigma }_{rr0}^{\mathrm{s}}=-P-S\cos (2{\theta }^{\prime }),{\sigma }_{\theta \theta 0}^{\mathrm{s}}=-P+S\cos (2{\theta }^{\prime }),\hspace{1em}{\tau }_{r\theta 0}^{\mathrm{s}}=S\sin (2{\theta }^{\prime }).$$

(4)

where $P=({\sigma }_{xx}+{\sigma }_{yy})/2,S=({\sigma }_{xx}-{\sigma }_{yy})/2,{\theta }^{\prime }=\theta -{\theta }_r.$ Variables with superscripts containing ‘i = ca, ce, s’ are related to corresponding casing, cement and stratum, respectively. Moreover, the subscript 0 denotes the initial stress state.

Owing to the linearity of elasticity theory and the in-situ far-field conditions, the problem can be decomposed into three sub-problems: (I) in-plane plane strain problem; (II) elastic uni-axial stress problem; (III) elastic anti-plane shear problem [31,32,33]. Then, the total solution of the inclined casing-cement-stratum system can be obtained by superposition. In addition, the tensile stress is taken as a positive.

The boundary conditions of three Sub-problems I, II and III can be listed as:

(1)

In-plane plane strain Sub-problem I: this sub-problem can be divided into an axisymmetric and asymmetric mode, which is widely applied.

Mode one, axisymmetric loading:

$$\begin{array}{l}{\left.{\sigma }_{rr}^{\mathsf{I}(1)}\right|}_{r=r_w}=-p_w,{\left.{\tau }_{r\theta }^{\mathsf{I}(1)}\right|}_{r=r_w}=0,\\ {\sigma }_{rr0}^{\mathrm{c}\mathrm{e},\mathsf{I}(1)}={\sigma }_{\theta \theta 0}^{\mathrm{c}\mathrm{e},\mathsf{I}(1)}=-{\sigma }_0^{\mathrm{c}\mathrm{e}}=-{\rho }_{\mathrm{c}\mathrm{e}}gh,\\ {\sigma }_{rr0}^{\mathrm{s},\mathsf{I}(1)}={\sigma }_{\theta \theta 0}^{\mathrm{s},\mathsf{I}(1)}=-P,\end{array}$$

(5)

Mode two, asymmetric loading:

$$\begin{array}{l}{\left.{\sigma }_{rr}^{\mathsf{I}(2)}\right|}_{r=r_w}=0,{\left.{\tau }_{r\theta }^{\mathsf{I}(2)}\right|}_{r=r_w}=0,\\ {\sigma }_{rr0}^{\mathrm{s},\mathsf{I}(2)}=-Scos(2{\theta }^{\prime }),{\sigma }_{\theta \theta 0}^{\mathrm{s},\mathsf{I}(2)}=Scos(2{\theta }^{\prime }),{\tau }_{r\theta 0}^{\mathrm{s},\mathsf{I}(2)}=Ssin(2{\theta }^{\prime }),\end{array}$$

(6)

(2)

Elastic uni-axial stress problem II:

$${\sigma }_{zz0}^{\mathrm{c}\mathrm{e},\mathsf{I}\mathsf{I}}=-{\sigma }_0^{\mathrm{c}\mathrm{e}}=-{\rho }_{\mathrm{c}\mathrm{e}}gh,{\sigma }_{zz0}^{\mathrm{s},\mathsf{I}\mathsf{I}}=-{\sigma }_{zz0}^{\mathrm{s}},$$

(7)

(3)

Elastic anti-plane shear problem III:

$$\begin{array}{l}{\left.{\sigma }_{rz}^{\mathrm{c}\mathrm{a},\mathsf{I}\mathsf{I}\mathsf{I}}\right|}_{r=r_w}=0,{\left.{\sigma }_{\theta z}^{\mathrm{c}\mathrm{a},\mathsf{I}\mathsf{I}\mathsf{I}}\right|}_{r=r_w}=0,\\ {\sigma }_{rz0}^{\mathrm{s},\mathsf{I}\mathsf{I}\mathsf{I}}=-({\sigma }_{xz0}^{\mathrm{s}}\cos \theta +{\sigma }_{yz0}^{\mathrm{s}}\sin \theta ),\\ {\sigma }_{\theta z0}^{\mathrm{s},\mathsf{I}\mathsf{I}\mathsf{I}}=({\sigma }_{xz0}^{\mathrm{s}}\sin \theta -{\sigma }_{yz0}^{\mathrm{s}}\cos \theta ),\end{array}$$

(8)

2.2. Solution Formulation

2.2.1. Solution of Mode One in Sub-Problem I

The classical solution for an elastic hollow cylinder subjected to axisymmetric loads is used for the casing, cement, and stratum. Chemical shrinkage during the cement setting process, which can affect the internal stress state in cement, is not considered in this paper since the stress state in the cement is caused by fluid column pressure and the deepest cement sets first. Then, the general solution for the stress and displacement field in increment form is given as [19]:

$$\begin{array}{l}\mathsf{\Delta }{\sigma }_{\mathrm{r}\mathrm{r}}^{\mathrm{i},\mathsf{I}(1)}=\frac{{\mathrm{A}}^{\mathrm{i}}}{2}+\frac{B^{\mathrm{i}}}{r^2},\\ \mathsf{\Delta }{\sigma }_{\theta \theta }^{\mathrm{i},\mathsf{I}(1)}=\frac{{\mathrm{A}}^{\mathrm{i}}}{2}-\frac{{\mathrm{B}}^{\mathrm{i}}}{r^2},\\ u_{rr}^{\mathrm{i},\mathsf{I}(1)}=\frac{1}{2G_{\mathrm{i}}}\left[(\frac{1}{2}-{\nu }_{\mathrm{i}}){\mathrm{A}}^{\mathrm{i}}\mathrm{r}-\frac{{\mathrm{B}}^{\mathrm{i}}}{r}\right],\end{array}$$

(9)

Then the total stress in Mode one can be written as:

$$\begin{array}{l}{\sigma }_{rr}^{\mathrm{i},\mathsf{I}(1)}={\sigma }_{rr0}^{\mathrm{i}}+\mathsf{\Delta }{\sigma }_{rr}^{\mathrm{i},\mathsf{I}(1)},\\ {\sigma }_{\theta \theta }^{\mathrm{i},\mathsf{I}(1)}={\sigma }_{\theta \theta 0}^{\mathrm{i}}+\mathsf{\Delta }{\sigma }_{\theta \theta }^{\mathrm{i},\mathsf{I}(1)}.\end{array}$$

(10)

Due to the infinity of displacement at infinity in the stratum, ${\mathrm{A}}_{\mathrm{s}}=0$ applies. In contrast, the existing models [19] treat the in-situ stress as boundary pressure loads instead of an initial stress state, thus Equations (9) and (10) can be rewritten as:

$$\begin{array}{l}{\sigma }_{rr}^{\mathrm{s},\mathsf{I}(1)}=\frac{{\mathrm{A}}^{\mathrm{s}}}{2}+\frac{B^{\mathrm{s}}}{r^2},\\ {\sigma }_{\theta \theta }^{\mathrm{i},\mathsf{I}(1)}=\frac{{\mathrm{A}}^{\mathrm{s}}}{2}-\frac{B^{\mathrm{s}}}{r^2},\end{array}$$

(11)

In this case, according to ${\left.{\sigma }_{rr}^{\mathrm{s},\mathsf{I}(1)}\right|}_{\mathrm{r}=\infty }=-P$ in Equation (5), ${\mathrm{A}}^{\mathrm{s}}=-2P$ applies. Consequently, the corresponding ultimate stress field will differ from that of the present model with ${\mathrm{A}}^{\mathrm{s}}=0$.

The normal traction and displacement continuity and boundary conditions at the casing-cement-stratum interfaces in Mode one are:

$$\begin{array}{l}{\left.{\sigma }_{rr}^{\mathrm{c}\mathrm{a},\mathsf{I}(1)}\right|}_{\mathrm{r}={\mathrm{R}}_1}=-P_{\mathrm{w}},\\ {\left.{\sigma }_{rr}^{\mathrm{c}\mathrm{a},\mathsf{I}(1)}\right|}_{\mathrm{r}={\mathrm{R}}_2}={\left.{\sigma }_{rr}^{\mathrm{c}\mathrm{e},\mathsf{I}(1)}\right|}_{\mathrm{r}={\mathrm{R}}_2},{\left.{\mathbf{u}}_{rr}^{\mathrm{c}\mathrm{a},\mathsf{I}(1)}\right|}_{\mathrm{r}={\mathrm{R}}_2}={\left.{\mathbf{u}}_{rr}^{\mathrm{c}\mathrm{e},\mathsf{I}(1)}\right|}_{\mathrm{r}={\mathrm{R}}_2},\\ {\left.{\sigma }_{rr}^{\mathrm{c}\mathrm{e},\mathsf{I}(1)}\right|}_{\mathrm{r}={\mathrm{R}}_3}={\left.{\sigma }_{rr}^{\mathrm{s},\mathsf{I}(1)}\right|}_{\mathrm{r}={\mathrm{R}}_3},{\left.{\mathbf{u}}_{rr}^{\mathrm{c}\mathrm{e},\mathsf{I}(1)}\right|}_{\mathrm{r}={\mathrm{R}}_3}={\left.{\mathbf{u}}_{rr}^{\mathrm{s},\mathsf{I}(1)}\right|}_{\mathrm{r}={\mathrm{R}}_3},\end{array}$$

(12)

Substituting Equations (9) and (10) into Equation (12), the linear equations can be given as:

$${\mathit{M}}^{\mathsf{I}(1)}\cdot {\mathit{X}}^{\mathsf{I}(1)}={\mathit{N}}^{\mathsf{I}(1)}$$

(13)

where ${\mathit{M}}^{\mathsf{I}(1)}$ is a five-by-five matrix whose non-zero entries are given in Appendix A, and:

$$\begin{array}{l}{\mathit{X}}^{\mathsf{I}(1)}={\left[\begin{array}{ccccc}{\mathrm{A}}^{\mathrm{c}\mathrm{a}}& {\mathrm{B}}^{\mathrm{c}\mathrm{a}}& {\mathrm{A}}^{\mathrm{c}\mathrm{e}}& {\mathrm{B}}^{\mathrm{c}\mathrm{e}}& {\mathrm{B}}^{\mathrm{s}}\end{array}\right]}^{\mathrm{T}},\\ {\mathit{N}}^{\mathsf{I}(1)}={\left[\begin{array}{ccccc}-P_w& -{\sigma }_0^{\mathrm{c}\mathrm{e}}& 0& {\sigma }_0^{\mathrm{c}\mathrm{e}}-P_0& 0\end{array}\right]}^{\mathrm{T}}.\end{array}$$

(14)

From Equation (13), the array ${\mathit{X}}^{\mathsf{I}(1)}$ consisting of the unknown coefficients can be solved straightforwardly. Then, the corresponding stress and displacement distributions within the casing-cement-stratum system in Mode one can be yielded. The vertical stress components induced by horizontal stress, assuming plane strain condition, in Mode one are given as:

$$\mathsf{\Delta }{\sigma }_{\mathrm{z}\mathrm{z}}^{\mathrm{i},\mathsf{I}(1)}=v_{\mathrm{i}}(\mathsf{\Delta }{\sigma }_{rr}^{\mathrm{i},\mathsf{I}(1)}+\mathsf{\Delta }{\sigma }_{\theta \theta }^{i,\mathsf{I}(1)})$$

(15)

2.2.2. Solution of Mode two in Sub-Problem I

In Mode two, the in-situ stress in the stratum varies with the polar angle ${\theta }^{\prime }$ in Equation (6). By solving the stress function, the solution of a hollow cylinder subjected to cosinoidal and sinusoidal distributions of normal and shear traction can be given as [20,21]:

$$\left\{\begin{array}{c}{\sigma }_{rr}^{i,\mathsf{I}(2)}={\sigma }_{rr0}^{i,\mathsf{I}(2)}+\mathsf{\Delta }{\sigma }_{rr}^{i,\mathsf{I}(2)}={\sigma }_{rr0}^{i,\mathsf{I}(2)}-2\left(B^i+\frac{3C^i}{r^4}+\frac{2D^i}{r^2}\right)\cos (2{\theta }^{\prime }),\\ {\sigma }_{\theta \theta }^{i,\mathsf{I}(2)}={\sigma }_{\theta \theta 0}^{i,\mathsf{I}(2)}+\mathsf{\Delta }{\sigma }_{\theta \theta }^{i,\mathsf{I}(2)}={\sigma }_{\theta \theta 0}^{i,\mathsf{I}(2)}+2\left(6A^i{r}^2+B^i+\frac{3D^i}{r^4}\right)\cos (2{\theta }^{\prime }),\\ {\sigma }_{r\theta }^{i,\mathsf{I}}={\sigma }_{r\theta 0}^{i,\mathsf{I}(2)}+\mathsf{\Delta }{\sigma }_{r\theta }^{i,\mathsf{I}(2)}={\sigma }_{r\theta 0}^{i,\mathsf{I}(2)}+2\left(3A^i{r}^2+B^i-\frac{C^i}{r^2}-\frac{3D^i}{r^4}\right)\sin (2{\theta }^{\prime }),\end{array}\right.$$

(16)

Moreover, the displacement components are given as [18,19]:

$$\left\{\begin{array}{c}{u}_{rr}^{i,\mathsf{I}(2)}=-\frac{2(1+{\nu }_i)}{E_i}(2v_i{A}^i{r}^3+B^{i}r-(1-v_i)\frac{2C^i}{r}-\frac{D^i}{r^3})\cos (2{\theta }^{\prime }),\\ u_{\theta \theta }^{i,\mathsf{I}(2)}=\frac{2(1+{\nu }_i)}{E_i}((3-2v_i)A^i{r}^3+B^{i}r-(1-2v_i)\frac{C^i}{r^3}+\frac{D^i}{r^3})\mathrm{s}in(2{\theta }^{\prime }),\end{array}\right.$$

(17)

Similarly, the vertical stress components in Mode two under the plane strain condition are given as:

$$\mathsf{\Delta }{\sigma }_{\mathrm{z}\mathrm{z}}^{\mathrm{i},\mathsf{I}(2)}=v_{\mathrm{i}}(\mathsf{\Delta }{\sigma }_{rr}^{\mathrm{i},\mathsf{I}(2)}+\mathsf{\Delta }{\sigma }_{\theta \theta }^{\mathrm{i},\mathsf{I}(2)})$$

(18)

In Mode two, the initial stress components are zero in the casing and cement, while the initial stress components are non-zero for the stratum. Since the increment of stress and the displacement of the stratum in Equations (16) and (17) at infinity is zero, we obtain:

$$A^{\mathrm{s}}=B^{\mathrm{s}}=0$$

(19)

The continuity of stresses and displacements at the interfaces and boundary conditions in Equation (6) give the following linear aligns:

$${\mathit{M}}^{\mathsf{I}(2)}\cdot {\mathit{X}}^{\mathsf{I}(2)}={\mathit{N}}^{\mathsf{I}(2)},$$

(20)

where ${\mathit{M}}^{\mathsf{I}(2)}$ is a 10-by-10 matrix whose non-zero entries are given in Appendix A, and:

$$\begin{array}{l}{\mathit{X}}^{\mathsf{I}(2)}={\left[A^{\mathrm{c}\mathrm{a}},B^{\mathrm{c}\mathrm{a}},C^{\mathrm{c}\mathrm{a}},D^{\mathrm{c}\mathrm{a}},A^{\mathrm{c}\mathrm{e}},B^{\mathrm{c}\mathrm{e}},C^{\mathrm{c}\mathrm{e}},D^{\mathrm{c}\mathrm{e}},C^{\mathrm{s}},D^{\mathrm{s}}\right]}^{\mathrm{T}},\\ {\mathit{N}}^{\mathsf{I}(2)}={\left[0,0,0,0,S/2,S/2,0,0,0,0\right]}^{\mathrm{T}},\end{array}$$

(21)

By substituting the unknown coefficients in ${\mathit{X}}^{\mathsf{I}(2)}$, the incremental and total stress distribution in Mode two of the system can be obtained. In contrast, the existing model [19] treats the in-situ stress as boundary pressure loads instead of an initial stress state; Equation (16) can be rewritten as:

$$\left\{\begin{array}{c}{\sigma }_r^{i,\mathsf{I}(2)}=-2\left(B^i+\frac{3C^i}{r^4}+\frac{2D^i}{r^2}\right)\cos 2{\theta }^{\prime },\\ {\sigma }_{\theta }^{i,\mathsf{I}(2)}=2\left(6A^i{r}^2+B^i+\frac{3D^i}{r^4}\right)\cos 2{\theta }^{\prime },\\ {\tau }_{r\theta }^{i,\mathsf{I}(2)}=2\left(3A^i{r}^2+B^i-\frac{C^i}{r^2}-\frac{3D^i}{r^4}\right)\sin 2{\theta }^{\prime },\end{array}\right.$$

(22)

According to the boundary conditions at infinity for the stratum in Equation (6), $A^s=0, B^s=S/2$ differing from those in corresponding Equation (19) applies. Correspondingly, the matrix ${\mathit{N}}^{\mathsf{I}(2)}$ in Equation (21) can be rewritten as:

$${\mathit{N}}_{\mathrm{o}\mathrm{l}\mathrm{d}}^{\mathsf{I}(2)}={\left[0,0,0,0,S/2,S/2,0,0,S\mathrm{c}/2/{\mathrm{G}}_{\mathrm{s}},S\mathrm{c}/2/{\mathrm{G}}_{\mathrm{s}}\right]}^{\mathrm{T}},$$

(23)

2.2.3. Comparison between the Present and Existing Model

Previous models considered in-situ stress as the boundary pressure load at infinity for the stratum, which led to displacements of the stratum at infinity being infinite, not reflecting the actual situation [20,21]. The present model improves upon this by considering in-situ stress as the initial stress state instead. For instance, the incremental radial displacement in the stratum in Sub-problem I includes the terms in Equations (9) and (17):

$$\left\{\begin{array}{c}{u}_{rr}^{s,\mathsf{I}(1)}=-\frac{(1-2v_{\mathrm{s}})}{4G_s}{A}^{\mathrm{s},\mathsf{I}(1)}r,\\ u_{rr}^{s,\mathsf{I}(2)}=\frac{2(1+{\nu }_{\mathrm{s}})}{E_s}{B}^{s,\mathsf{I}(2)}r\cos (2{\theta }^{\prime }),\end{array}\right.$$

(24)

In Equation (24) $A^{s,\mathsf{I}(1)}=-2P,B^{s,\mathsf{I}(2)}=S/2$ applies in the existing model while $A^{s,\mathsf{I}(1)}=B^{s,\mathsf{I}(2)}=0$ applies in the present model. So the displacements at infinity in the stratum will be infinity rather than zero in the existing model. This difference between the previous and present model ensures that the displacement and stress distribution in the stratum are modeled more accurately, and this applies to Sub-problem III.

2.2.4. Solution of Sub-Problem II: Elastic Uni-Axial Stress Problem

Based on the plane strain assumption, vertical stress can be given as:

$${\sigma }_{\mathrm{z}\mathrm{z}}^i={\sigma }_{\mathrm{z}\mathrm{z}0}^i+{\nu }_{\mathrm{i}}(\mathsf{\Delta }{\sigma }_{\mathrm{z}\mathrm{z}}^{\mathrm{i},\mathsf{I}(1)}+\mathsf{\Delta }{\sigma }_{\mathrm{z}\mathrm{z}}^{\mathrm{i},\mathsf{I}(2)})$$

(25)

2.2.5. Solution of Sub-Problem III: Anti-Plane Shear Problem

In this problem, the Navier equation of equilibrium is [16]:

$${\nabla }^2{u}_{\mathrm{z}}^{\mathrm{i}}=0,$$

(26)

where $u_{\mathrm{z}}^{\mathrm{i}}$ is the vertical displacement component. Given that the wellbore is long, $\partial u_{\mathrm{i}}/\partial z=0$ applies, and we can assume that $u_{\mathrm{z}}$ is a function of the radial coordinate $r$ and polar angle $\theta$. Then, the solution of Equation (26) has the following forms [18]:

$$u_{\mathrm{z}}^{\mathrm{i}}={\displaystyle \sum _{\mathrm{n}=1}^{\infty }({\mathrm{A}}_{\mathrm{n}}^{\mathrm{i}}\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{n}\theta )+{\mathrm{B}}_{\mathrm{n}}^{\mathrm{i}}\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{n}\theta ))}{r}^{\mathrm{n}}+{\displaystyle \sum _{\mathrm{n}=1}^{\infty }({\mathrm{C}}_{\mathrm{n}}^{\mathrm{i}}\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{n}\theta )+{\mathrm{D}}_{\mathrm{n}}^{\mathrm{i}}\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{n}\theta ))}{r}^{-\mathrm{n}}$$

(27)

Using symmetry considerations, it can be suggested that the dependence of the stress upon the polar angle $\theta$ in Equation (8) makes $\mathrm{n}=1$. Therefore, the stress components are derived as follows:

$$\begin{array}{l}{\sigma }_{rz}^{\mathrm{i}}={\sigma }_{rz0}^{\mathrm{i}}+\mathsf{\Delta }{\sigma }_{rz}^{\mathrm{i}}={\sigma }_{rz0}^{\mathrm{i}}+G_{\mathrm{i}}((A^{\mathrm{i}}-\frac{C^{\mathrm{i}}}{r^2})\sin \theta +(B^{\mathrm{i}}-\frac{D^{\mathrm{i}}}{r^2})\cos \theta ),\\ {\sigma }_{\theta z}^{\mathrm{i}}={\sigma }_{\theta z0}^{\mathrm{i}}+\mathsf{\Delta }{\sigma }_{\theta z}^{\mathrm{i}}={\sigma }_{\theta z0}^{\mathrm{i}}+G_{\mathrm{i}}((A^{\mathrm{i}}+\frac{C^{\mathrm{i}}}{r^2})\cos \theta -(B^{\mathrm{i}}+\frac{D^{\mathrm{i}}}{r^2})\mathrm{s}\mathrm{i}n\theta ),\end{array}$$

(28)

According to the boundary and continuity conditions, linear equations can be obtained for the stress distributions induced by the initial stress in Sub-problem III. Additionally, since the displacement of the stratum at infinity must be zero, ${\mathrm{A}}^{\mathrm{s}}={\mathrm{B}}^{\mathrm{s}}=0$ applies in Equations (27) and (28). Therefore, the corresponding linear aligns are given as:

$$\begin{array}{l}\mathit{K}\cdot {\mathit{X}}_1^{\mathsf{I}\mathsf{I}\mathsf{I}}={\mathit{L}}_1,\\ \mathit{K}\cdot {\mathit{X}}_2^{\mathsf{I}\mathsf{I}\mathsf{I}}={\mathit{L}}_2,\end{array}$$

(29)

where

$$\begin{array}{l}{\mathit{X}}_1^{\mathsf{I}\mathsf{I}\mathsf{I}}={\left[A^{\mathrm{c}\mathrm{a}},C^{\mathrm{c}\mathrm{a}},A^{\mathrm{c}\mathrm{e}},C^{\mathrm{c}\mathrm{e}},C^{\mathrm{s}}\right]}^{\mathrm{T}},{\mathit{X}}_2^{\mathsf{I}\mathsf{I}\mathsf{I}}={\left[B^{\mathrm{c}\mathrm{a}},D^{\mathrm{c}\mathrm{a}},B^{\mathrm{c}\mathrm{e}},D^{\mathrm{c}\mathrm{e}},D^{\mathrm{s}}\right]}^{\mathrm{T}},\\ {\mathit{L}}_1={\left[0,0,S_{yz0}^{\mathrm{s}},0,0\right]}^{\mathrm{T}},{\mathit{L}}_2={\left[0,0,S_{xz0}^{\mathrm{s}},0,0\right]}^{\mathrm{T}},\end{array}$$

(30)

and $\mathit{K}$ is listed in Appendix A.

In contrast to the present model, the existing model assumes that the sum displacement and stress at the interface is continuous. The coefficients in Equations (27) and (28) can be given as $A_{old}^s=S_{yz}/G_s,B_{old}^s=S_{xz}/G_s$, according to the initial stress state in Equation (8), following a similar approach as in Sub-problem I. However, this leads to the displacement at infinity of the stratum being infinite due to the term $(A_{old}^s\sin \theta +B_{old}^s\cos \theta )r$, which is not realistic and will be discussed later. Correspondingly, the modified parameters ${\mathit{L}}_1,{\mathit{L}}_2$ in the existing model can be given as:

$${\mathit{L}}_{1\mathrm{o}\mathrm{l}\mathrm{d}}={\left[0,0,S_{yz0}^{\mathrm{s}},0,S_{yz0}^{\mathrm{s}}\mathrm{c}/{\mathrm{G}}_{\mathrm{s}}\right]}^{\mathrm{T}},{\mathit{L}}_{2\mathrm{o}\mathrm{l}\mathrm{d}}={\left[0,0,S_{xz0}^{\mathrm{s}},0,S_{xz0}^{\mathrm{s}}\mathrm{c}/{\mathrm{G}}_{\mathrm{s}}\right]}^{\mathrm{T}}$$

(31)

2.3. The 3D Total Stress Field within the Casing-Cement-Stratum System

The ultimate stress components can be obtained as a superposition solution of three Sub-problems: in-plane plane strain, uni-axial stress and anti-plane shear problems. Additionally, the initial stress of cement is considered, which affects the stress field and failure potential. The ultimate stresses can be given as:

$$\begin{array}{l}{\sigma }_{rr}^i={\sigma }_{rr0}^{\mathrm{i}}+\mathsf{\Delta }{\sigma }_{rr}^{i,\mathsf{I}(1)}+\mathsf{\Delta }{\sigma }_{rr}^{i,\mathsf{I}(2)},\\ {\sigma }_{\theta \theta }^i={\sigma }_{\theta \theta 0}^i+\mathsf{\Delta }{\sigma }_{\theta \theta }^{i,\mathsf{I}(1)}+\mathsf{\Delta }{\sigma }_{\theta \theta }^{i,\mathsf{I}(2)},\\ {\sigma }_{r\theta }^i={\sigma }_{r\theta 0}^i+\mathsf{\Delta }{\sigma }_{r\theta }^{i,\mathsf{I}(2)},\\ {\sigma }_{zz}^i={\sigma }_{zz0}^i+{\nu }_i(\mathsf{\Delta }{\sigma }_{zz}^{\mathrm{i},\mathsf{I}(1)}+\mathsf{\Delta }{\sigma }_{zz}^{\mathrm{i},\mathsf{I}(2)}),\\ {\sigma }_{rz}^i={\sigma }_{rz0}^i+\mathsf{\Delta }{\sigma }_{rz}^i,\\ {\sigma }_{\theta z}^i={\sigma }_{\theta z0}^i+\mathsf{\Delta }{\sigma }_{\theta z}^i.\end{array}$$

(32)

3. Results and Discussions

3.1. Validation of Present Analytical Model

The purpose of this paper is to introduce a straightforward and comprehensive 3D analytical solution for the stress field of the casing-cement-stratum system, taking into account the initial stress state of cement and the stratum. The input parameter values used in the simulation, listed in

Table 1. Input parameters (Data from Liu et al. [20]).

Component	Casing	Cement	Stratum
Elastic coefficients and geometry
Young’s modulus $E/\mathrm{G}\mathrm{P}\mathrm{a}$	210	12	10
Poisson’s ratio $\nu$	0.3	0.25	0.2
Inner radius/mm	59.31	69.85	107.95
Stress state
Initial vertical stress ${\sigma }_V/\mathrm{M}\mathrm{P}\mathrm{a}$	0	0	−100 (assumed)
Initial horizontal stress ${\sigma }_H,{\sigma }_h/\mathrm{M}\mathrm{P}\mathrm{a}$	0	0	−82, −55
Internal pressure $p_w/\mathrm{M}\mathrm{P}\mathrm{a}$	−34

, are obtained from the literature of Liu et al. [20]. The new analytical solution to the present model in Sub-problem I is compared with the Liu model [20], which considers the initial displacement induced by the in-situ stresses in the stratum and ignores the hydrostatic stress in the cement. In contrast, this paper uses the initial stress in the stratum and cement directly to obtain the ultimate stress field, resulting in a more straightforward and understandable approach.

The accuracy of the model results in Sub-problem I, which considers the in-situ stress as the initial stress state (present model) or as boundary pressure loads (existing model); this is validated by excellent agreement with existing results and finite element analysis (FEA) [20], as shown in Figure 3. The existing model considers the initial stress in the stratum as the boundary pressure loads, which causes great error compared to the present model. A detailed analysis has been conducted by Liu et al. [20] and Li et al. [21], which is not discussed here.

Furthermore, the existing model leads to an infinite radial displacement at infinity in the stratum, which is clearly unrealistic, as shown in Figure 4. In contrast, in the present model, as the radius increases, the radial displacement approaches zero. The reason lies in the fact that the initial stress in the stratum has been exerted before wellbore completion and cannot be regarded as the boundary pressure loads, which is always assumed in the existing model.

3.2. Model Comparison between the Present Solution and Existing One in Antiplane Shear Problem III

In the case of an inclined wellbore with a deviation azimuth of 30° and a deviation angle of 60°, the in-situ anti-plane shear stresses in the stratum can be given as ${\sigma }_{xz0}^s=16.6 \mathrm{M}\mathrm{P}\mathrm{a},{\sigma }_{yz0}^s=5.8 \mathrm{M}\mathrm{P}\mathrm{a}$. The induced anti-plane stresses using both the present and existing analytical models are shown in Figure 5. Similar to Sub-problem I, the existing analytical model leads to significant errors compared to the corresponding present model in Sub-problem III.

Equations (29) and (30) imply that the coefficients ${\mathit{X}}_1^{\mathsf{I}\mathsf{I}\mathsf{I}},{\mathit{X}}_2^{\mathsf{I}\mathsf{I}\mathsf{I}}$ are proportional, resulting in the same nondimensionless anti-plane shear stress as the ${\sigma }_{rz}/{\sigma }_{xz0}^s,{\sigma }_{rz}/{\sigma }_{yz0}^s$ distributions when considering the dependence of cosinoidal and sinusoidal anti-plane shear stress in Equation (8). However, the existing model greatly amplifies the induced anti-plane shear stresses ${\sigma }_{rz},{\sigma }_{\theta z}$ compared to the corresponding present model, especially in casing and cement. Therefore, the present analytical model is recommended to accurately calculate the stress field in the casing-cement-stratum system in all Sub-problems. Additionally, the existing model in Sub-problem III causes linearly increasing radial displacement in the stratum in Figure 6, which is also unrealistic at infinity.

Similar to the result of Sub-problem I, a low Young’s modulus of cement will decrease the anti-plane shear stress ${\sigma }_{rz},{\sigma }_{\theta z}$ in the casing and cement in both the present and existing models, as shown in Figure 7. This reduction will benefit wellbore integrity. However, Poisson’s ratio of cement has little influence on the anti-plane shear stress.

3.3. The Effect of Initial Stress State of Cement on the Stress Field

In previous models [7,8,9,10,11,12,13,14,15,16,17,18,19,20,21], the initial stress of the cement was often ignored. However, the existence of initial stress in cement has a significant influence on the stress distribution, as shown in Figure 8. Due to the isotropic initial stress state in the cement, the effect of the initial stress state of cement is studied in Mode one of Sub-problem I. The initial stress state in cement promotes the uniformity of radial and hoop stress along the radial direction and increases the absolute values of compressive radial and hoop stresses. Moreover, for the given parameters in this case, the initial stress state in cement has little influence on the stress field in the casing.

The effect of Young’s modulus and Poisson’s ratio on cement in the stress field in Mode one of Sub-problem I, considering the initial stress state of cement, is studied in Figure 9. Considering the initial stress state of cement, a low Young’s modulus also induces low absolute values of radial and hoop stress. Ignoring the initial stress state of cement will amplify the decreasing degree of absolute values in radial and hoop stresses in the casing and cement under low Young’s modulus of cement condition.

4. Estimation of the Casing and Cement Integrity

The Mohr-Coulomb failure criterion has been widely used to evaluate cement failure potential, while the von Mises yield criterion is used to assess longitudinal zonal isolation in the cement [34]. In this paper, the 3D von Mises yield criterion is adopted to evaluate the failure potential of casing and cement:

$${\sigma }_{\mathrm{Mises}}=\sqrt{\frac{1}{2}\left[{({\sigma }_r-{\sigma }_{\theta })}^2+{({\sigma }_r-{\sigma }_z)}^2+{({\sigma }_z-{\sigma }_{\theta })}^2\right]+3{\tau }_{r\theta }^2+3{\tau }_{rz}^2+3{\tau }_{\theta z}^2}$$

(33)

Correspondingly, the 2D von Mises yield criterion is simplified as:

$${\sigma }_{\mathrm{Mises},2\mathrm{D}}=\sqrt{\frac{1}{2}({({\sigma }_r-{\sigma }_{\theta })}^2+{\sigma }_r^2+{\sigma }_{\theta }^2+6{\tau }_{r\theta }^2)}$$

(34)

The Mises stress in the 2D and 3D solution at the inner surfaces of casing and cement are evaluated by the present and existing models, as shown in Figure 10. In this case, an inclined wellbore with a deviation azimuth of 30°, a deviation angle of 60° and zero initial stress of cement is considered. The symmetry of von Mises stress in the 2D solution, ${\theta }_r$, validates the accuracy of the model results in Sub-problem I. In contrast, the existence of an anti-plane shear stress distribution in Sub-problem III breaks the symmetry of the Mises stress distribution, ${\theta }_r$, in the casing and cement.

Consistent with the law above, the existing model greatly signifies the values of Mises stress in casing and cement both in the 2D and 3D solution. The 3D solutions of the present and existing models differ from the corresponding 2D solutions, which highlights the necessity of establishing the 3D analytical solution to evaluate wellbore integrity. Counter-intuitively, the 2D Mises stress in both the present and existing models is higher than that in the 3D case, particularly for the cement. This suggests that the common 2D solution overestimates the failure potential of cement compared to the 3D model. This could be because ignoring the existence of stress components with axial direction increases the nonuniform degree of stress components.

Figure 11 demonstrates how Young’s modulus and the initial stress state of cement affect the maximum Mises stress in casing and cement in the present 3D model. As Young’s modulus increases, the maximum Mises stress in the casing, and cement approaches a constant value. Simultaneously, the presence of initial stress in cement significantly decreases the maximum Mises stress in the cement with moderate Young’s modulus values of cement. This occurs because the initial hydrostatic compressive stress state in cement moves Mohr’s circle of stress in cement away from the failure envelope compared to the case of zero initial stress. Additionally, under a low Young’s modulus of cement conditions, the presence of initial stress in cement aggravates the failure potential of the casing, and vice versa.

Figure 12 illustrates the influence of wellbore pressure on the maximum Mises stress in casing and cement in the present 3D model. As the wellbore pressure increases, the failure potential of the casing decreases while that of the cement increases. Furthermore, the presence of initial stress in cement significantly decreases the increased failure potential of cement with increasing internal wellbore pressure. The initial stress state has little influence on the maximum Mises stress in the casing since radial and hoop stresses hardly change with increasing initial stress in cement for the given specific parameters, as shown in Figure 8.

5. Conclusions

This paper presents a new 3D analytical model for evaluating the stress field in inclined oil/gas wells, considering in-situ stress as the initial stress instead of boundary pressure loads and the existence of the initial stress of cement. The 3D analytical solution is obtained by dividing the model into an in-plane plane strain problem, elastic uni-axial stress problem and anti-plane shear problem. Moreover, the accuracy of the present model is validated by the existing model results and finite element analysis. The following main conclusions can be made:

(1)

Compared to the existing model, the present model provides a more accurate stress level in the anti-plane shear problem and eliminates unrealistic infinite displacements at infinity.

(2)

The maximum Mises stress in the present 3D model differs from that in the corresponding 2D model, highlighting the importance of considering the uni-axial stress problem and anti-plane shear problem. For the specific conditions studied in this paper, the maximum Mises stress in the cement in the present 3D model is lower than that in the corresponding 2D model, resulting in a more conservative failure potential in the 3D circumstance.

(3)

The existence of initial stress in cement has a significant effect on the cement failure potential. For the specific conditions studied in this paper, as the wellbore pressure increases, the failure potential of the casing decreases while that of cement increases. However, the presence of initial stress in cement significantly decreases the increased failure potential of cement with increasing internal wellbore pressure.