Volume fracturing of horizontal well has become the key technology in shale gas development, and core techniques, represented by "slick water at high pumping rate+soluble bridge plug" and "intensive staged clustering+temporary plugging and diverting", have been widely used^[1,2,3]. However, casing deformation (CD) is likely to occur in volume fracturing of horizontal wells^[4]. By the end of 2019, 187 wells had been fractured in Changning- Weiyuan area (C-W area), of them, 75 had CD, accounting for 40.1%. Consequently, fracturing was unable to implement in 174 stages^[5]. CD will affect the wellbore integ-rity, hinder running down of tools, and even cause some invalid stages, greatly affecting the stimulation efficiency^[6,7].

Some studies suggest that the main factors inducing CD in shale gas horizontal wells include natural fracture slippage, casing yield collapse and thermal stress damage^[8]. In terms of research on CD caused by natural fracture (NF) slippage shear, Liao et al. analyzed CD characteristics of wells in C-W area based on field stimulation data, and concluded that NF slippage shear was the main cause of CD^[9,10,11]. Gao et al. analyzed the influences of fracturing zone, formation slippage, cementing quality, natural fracture and other factors on CD with simulation software^{[12,13,14,15,16]}. Fu et al. analyzed the influence of NF parameters on CD through large-scale physical modeling experiments^[17,18]. In general, the previous researches on NF slippage shear induced CD mostly focused on the analysis of the phenomenon and laws of CD, but few studies on NF slipage, stress on casing, deformation mechanism, CD rapid prediction, and response to CD risk during fracturing have been conducted. In this work, in light of NF slip shear induced CD, a formation-fracture-casing stress unit has been established to analyze the mechanical mechanism of CD, a CD calculation model based on the complex function has been worked out to analyze influence factors of CD. The research results can provide an analytical tool for rapid prediction of the deformation position and deformation severity of horizontal well in fracturing engineering design.

High-angle natural fractures are common in shale reservoirs of southern Sichuan Basin^{[11, 19-21]}. During the fracturing process of shale reservoir, due to the existence of flow channels such as hydraulic fractures or voids in cement, fracturing fluid may flow into the NF intersecting with the wellbore and make the fluid in fractures rise in pressure^[5]. When fluid pressure in fractures increases to a critical value, the resultant force of fracture surface friction and wellbore shear force (the shear force of the wellbore against fracture slippage) will be less than the shear force of formation slippage, and fracture will slip and cause casing deformation. Based on the above physical process, as shown in Fig. 1, the stress model of formation-fracture-casing system was established to study the mechanical mechanism of CD caused by the fracture slippage-shear. In the figure, the angle between the direction of the maximum horizontal principal stress (Oe direction) and the fracture surface in the clockwise direction (referred to as the fracture approach angle) is θ, and the angle between the Oe direction and the horizontal wellbore in the clockwise direction (referred to as the wellbore approach angle) is α. The model assumes that: (1) The horizontal wellbore (in the same horizontal plane) passes through a vertical NF, and fracturing fluid flows into the NF from the wellbore through a flow channel. (2) The horizontal wellbore is composed of casing and cementing sheath; the whole cementing sheath and casing bear the shear stress together, but damaged cementing sheath cannot bear the shear stress. (3) The NF is uncemented, the fluid fills natural fracture when it enters; (4) Fracturing fluid leak-off and fracture propagation are not considered, fluid pressure is equal everywhere in the NF.

Since the NF in the model is assumed to be vertical, the formation vertical stress has no shear stress component on the fracture surface, and the effect of horizontal principal stress on the fracture surface is the main interest of the study. A horizontal square section of casing and formation (Fig. 1b) in the formation-fracture-casing system stress unit model (Fig. 1a) is selected for stress analysis, and 4 sides of the stress unit are indicated by the dotted lines a', b', c', and d' with the length equal to fracture length. In the figure, the x axis is perpendicular to the fracture surface, and the y axis is parallel to the fracture surface. According to the plane stress analysis, the normal stress in the x direction and shear stress in the y direction of the stress unit before fracturing operation can be expressed as follows:

1.1.1. Partially opened fracture model

Fracturing fluid enters and supports NFs during hydraulic fracturing. When the fracture is partially opened, as shown in Fig. 2a, the fracture surface in the x direction is mainly subjected to the contact force of rocks and the fracture fluid pressure, and the main forces in the y direction are the fracture surface friction force and the wellbore shear force. Assuming that the stress unit is stable, the force balance established in the x direction can be expressed as:

The force balance in the y direction can be expressed as:

Since the near-wellbore stress has been released before cementing, and the cement solidification stress is ignored, when the fracture surface friction force is greater than the shear force on formation in the y direction (b' surface), from Eq. (3), the actual friction force (f) on the fracture surface is equal to the shear force on formation, and the wellbore shear force is zero.

When the fracture surface friction force is less than the shear force on formation in the y direction, the formation will slip along the fracture surface if there is no wellbore in the formation. And if existing, the wellbore will bear part of the shear force to resist the formation slippage, in this case, f is the maximum static friction force (f_max). The shear stress of the wellbore against fracture slippage under current conditions (referred to as the wellbore shear stress) can be calculated from Eq. (3):

In the above equation, f_max can be calculated by the fracture surface friction coefficient μ, and the values of which can be 0.6-1.0^[5]:

Assuming that the stress unit is stable, fracture geometric parameters are fixed, the relationship between the fluid pressure inside fracture and the forces in the tangential direction of fracture is shown in Fig. 3. In the figure, with the increase of p_f, f_max decreases linearly, f is equal to the shear force on formation in the y direction ${{\tau }_{xy}}{{A}_{\text{f}}}$, and the wellbore shear force is zero. When f_max decreases to ${{\tau }_{xy}}{{A}_{\text{f}}}$, the wellbore will support part of the shear force (${{\tau }_{\text{c}}}{{A}_{\text{c}}}$) to keep the stress unit static, and f is equal to f_max. When p_f increases to${{\sigma }_{x}}$, the fracture is fully opened by the fracturing fluid, the rocks on both sides of the fracture are no longer in contact, and f is zero, at this point, the wellbore shear force is equal to${{\tau }_{xy}}{{A}_{\text{f}}}$.

1.1.2. Model of fully opened fracture

When the NF is fully opened by fracturing fluid, the rocks on both sides of fracture are no longer in contact

(Fig. 2b), the fracture surface is only subjected to p_f in the x direction, and f is zero in the y direction, and the shear force in this direction is borne by the wellbore. Therefore, the wellbore shear stress under current conditions can be calculated from Eq. (3):

A_c is the area enclosed by the cementing sheath on the section of wellbore and fracture:

A_c is the area enclosed by the casing (${{{A}'}_{\text{c}}}$) when the cementing sheath is damaged:

Tables 1 and 2 show the basic welllbore parameters and geological engineering parameters of shale reservoir stimulation in a block of Weiyuan area. In this area, the direction of ground stress is E-W, the maximum horizontal principal stress (σ_H) is about 75 MPa, the minimum horizontal principal stress (σ_h) is about 65 MPa, and the direction of horizontal wellbore is approximately parallel to the direction of σ_h. Comprehensive interpretation results of the ant-body and logging data show that NFs in this block are mostly in N-E strike and about 150 m long on average. The effects of factors, such as approach angle, fracture area, fracture fluid pressure, and fracture friction coefficient on wellbore shear stress were analyzed based on the above data.

1.2.1. The approach angle

Fig. 4 shows the relationship curves of fracture approach angle (θ), wellbore approach angle (α) and wellbore shear stress (τ_c) under the conditions of A_f=10 m², p_f=75 MPa, cases of intact cementing sheath and damaged cementing sheath. When the fracture trend is parallel or perpendicular to the direction of the maximum horizontal principal stress (σ_H), τ_c is equal to zero, which is independent of α. In other cases, the value of τ_c is affected by both θ and α. Wellbores in Weiyuan area are along the minimum horizontal stress direction in general (α is approximately equal to 90°). The corresponding θ of high shear stress area (τ_c≥1 GPa) in this case is 20°-55° (or its supplementary angle) when α=90°, which is consistent with the NF approach angle at CD position in Table 2. Thus, the value of α can be appropriately adjusted according to the chart to avoid the impact of high shear stress on wellbore at a given θ.

Besides, the integrity of the cementing sheath has a significant impact on the wellbore shear stress. In this case, the maximum shear stress on wellbore (α=135°, θ=45°) with damaged cementing sheath is 5.97 times that with intact cementing sheath. Fig. 4a shows that the value of τ_c is much higher than the shear strength (S_s) of the cementing sheath. Thus, it can be inferred that the cementing sheath integrity will be difficult to maintain when the formation slips along the fracture surface, so the focus of the next analysis is the CD in the case with damaged cementing sheath. Moreover, high steel-grade and thick casing is used in this case, the τ_c value in the case of damaged cement sheath is compared with the casing shear strength (S_c). It is found that the τ_c value is much higher than the casing shear strength (Fig. 4b), which shows that enhancing casing strength or cementing quality has limited effect on reducing the risk of CD^[22].

1.2.2. Fracture area

Fig. 5 shows the relationship curves between fracture area (A_f) and τ_c under α in different values of a fully opened fracture with θ=30°. In the figure, the larger the A_f, the larger the τ_c is; and the value of τ_c first increases and then decreases with the increase of α. From Eq. (6), at given τ_xy and A_c, τ_c is mainly affected by A_f, namely, the larger the NF area, the higher the risk of CD.

1.2.3. Fluid pressure in fractures

Fracture fluid pressure will affect natural fracture slippage when the fracture is partially opened. Fig. 6 shows the relationship curves between θ and τ_c under different p_f conditions when α=90°, μ=0.6, and A_f=50 m². When the p_f value is higher than σ_x, τ_c reaches the maximum, and when p_f is lower than σ_x-μ^-1τ_xy, τ_c becomes zero at any θ conditions. Taking the case of θ=30° as an example, σ_x is calculated to be 67.5 MPa and σ_x-μ^-1τ_xy is 60.3 MPa. In the figure, τ_c increases with the increase of p_f, the τ_c curves at p_f of 70, 75, and 80 MPa overlap at the position of θ=30°, and τ_c reaches the maximum value when the p_f value is higher than 67.5 MPa. In contrast, when the p_f value is lower than 60.3 MPa, τ_c is zero.

1.2.4. Fracture friction coefficient

Fig. 7 shows the relationship between θ and τ_c under different fracture friction coefficient (μ) conditions of a partially opened fracture when α=90°, p_f=65 MPa, and A_f=50 m². It can be seen at given ground stress and fracture geometric parameters, f_max increases with the increase of μ, resulting in decrease of τ_c. However, because of the limited range of μ, the influence of the variation of μ on wellbore shear stress is not as significant as that of p_f.

Based on the above analysis, it is believed that NF is the main factor inducing CD in horizontal wells. Therefore, strengthening the monitoring and identifica-tion of NF in the block will help prevent CD in engineering: (1) Placing the wellbore parallel or perpendicular to the trend of natural fractures and avoid large NFs can reduce the wellbore shear force and lower the risk of CD based on a deep understanding of the geological conditions of the block; (2) Reducing net pressure in the fracture through engineering measures is an effective method to reduce the risk of CD^[9].

“Multi-cluster perforation+temporary plugging and diverting” technology applied in the horizontal well fracturing in Weiyuan area has achieved good results in reducing the risk of CD. The proportion of CD wells treated by this technology is 21%, which is far less than that of wells treated by the coventional technology of 48%. In blocks with known geological engineering parameters, the method proposed in this paper can be used to establish a chart and quickly estimate the risk of CD. Fluid pressure in the fracture can be reasonably controlled through engineering measures, such as optimizing the wellbore azimuth, controlling pumping rate, and applying temporary plugging, to reduce the risk of CD.

The rapid determination of CD is of great significance to the engineering design of CD well. Since the fracture slip shear stress is much greater than wellbore shear strength, it can be assumed that the CD is approximately equal to the fracture slippage in the tangential direction. Meanwhile, assuming that: (1) NF is the “I+II” composite type fracture^[23] (the opened fracture with tensile stress is defined as type I, and the slip fracture with shear stress is defined as type II), as shown in Fig. 8, and the fluid pressure in the fracture is p_f; (2) the fracture is fully opened under p_f, fluid pressure is equal everywhere inside the fracture, and fluid leak-off and fracture propagation are ignored. Based on the above assumptions, a CD calculation model is established based on the complex function method.

The Westergaard stress function of “I+II” type fracture can be expressed as^[23]:

Based on the Muskhelishvili fracture displacement calculation method^[23], the displacement field of the “I+II” type fracture can be calculated by Eq. (10):

The fracture relative displacement can be derived from Eq. (9) and Eq. (10):

And the fracture relative displacement under the additional fluid pressure in fracture can be derived as:

Thus, the fracture relative displacement of “I+II” type fracture can be derived from Eq. (11) and (12) as:

During the perforation operation of Well HX-1 in Wei 202 block, the bridge plug was obstructed when pumped to 3090 m downhole. The result of ant-body interpretation showed that a NF (150 m long, 50° N by E in azimuth) across the wellbore at the sticking point, the wellbore was about 30 m away from the fracture midpoint. Well temperature logging data showed that the temperature at the sticking point dropped greatly, and caliper logging data showed that grade-3 CD occurred at the sticking point (Tables 3 and 4). As shown in Fig. 9, the calculated CD (Δu_f, relative displacement in tangential direction of “I+II” type fracture) of Well HX-1 is 33.5 mm based on the above data. And the caliper logging data showed that the maximum CD of this well is 31.8 mm. The relative error of calculated CD is 5.3%, indicating that the calculation model is applicable.

2.3.1. Rock mechanical parameters

Fig. 10 shows the relationship curves of k'/a and Δu_f under different elasticity moduli at θ=45°, υ=0.25, and a=50 m. The figure shows the calculated CD varies in parabolic shape, and the maximum CD occurs at the fracture midpoint; under the same stress, the higher the elasticity modulus, the smaller the rock strain, so the larger the elasticity modulus, the smaller the CD will be.

With θ and a fixed, the influence of rock Poisson’s ratio on CD at E=30 GPa is shown in Fig. 11. The figure shows the CD is negatively correlated with Poisson’s ratio, but Poisson’s ratio has a limited effect on CD.

2.3.2. Fracture geometric parameters

Fig. 12 shows the influence of fracture half-length on CD at θ and υ with fixed values, and E of 30 GPa. The figure shows fracture half-length has a great impact on CD. This is because the larger the a, the larger the fracture area, the greater the shear force on formation, and the greater the formation displacement will be, resulting in larger CD.

With υ and a fixed, the influence of fracture approach angle on CD at E of 30 GPa is shown in Fig. 13. Since the shear stress is a sine function of θ, only the variation of the CD at θ∈[0°, 90°] is discussed. The calculation results show that the CD first increases and then decreases with the increase of θ, and the maximum value appears at θ of 45°.

The fracture and wellbore approach angles have significant impacts on the wellbore shear stress. Under the condition of common wellbore approach angles in the Weiyuan shale gas field, under fracture approach angles of 20° to 55° (or supplementary angles), the wellbore is subjected to large shear stress and high risk of CD.

When the fracture is partially opened by fracturing fluid, the wellbore shear stress is positively correlated with the fluid pressure in fracture, and negatively correlated with the fracture friction coefficient. When the fracture is fully opened, the wellbore shear stress is positively correlated with the natural fracture area.

The lower the rock elasticity modulus and the longer the fracture length, the more serious the CD will be. Poisson's ratio has a weak influence on the CD. The CD first increases and then decreases with the increase of the fracture approach angle, and reaches the maximum at the fracture approach angle of 45°.

For a given fracture approach angle, appropriately adjust the wellbore approach angle can avoid high shear stress on the wellbore. Reasonable control of the fluid pressure in fracture can reduce the CD risk. The shear stress on casing is usually much greater than casing shear strength, increasing casing strength or cementing quality has limited effect on reducing the risk of CD.

Caliper logging data has verified the CD prediction model, and the model can be used to establish risk analysis charts and calculate CD, which provides reference for quick CD risk prediction in fracturing design.

a—NF half length, m;

a', b', c', d'—the 4 sides of the stress unit;

A_c,${{{A}'}_{\text{c}}}$—the area enclosed by the cementing sheath and casing on the section of wellbore and fracture, m²;

A_f—NF area, m²;

d, d₁, d₂—the distance from z point to the midpoint and two endpoints of the long axis of the fracture in “I+II” type fracture plain stress model, m;

E—rock elasticity modulus, MPa;

E'—rock elasticity modulus under the plain strain condition, MPa;

f—actual friction force on fracture surface, 10⁶ N;

f_max—the maximum static friction force on fracture surface, 10⁶ N;

h—NF half height, m;

Im—imaginary part of complex function;

k—any point on the X axis;

k°—X coordinate value of k point, m;

p_f—fluid pressure in fracture, MPa;

r₁, r₂—inner diameter and outer diameter of casing, m;

R—wellbore diameter, m;

Re—real part of complex function;

S_c—casing shear strength, MPa;

S_s—cementing sheath shear strength, MPa;

u, v—tangential and normal displacement of “I+II” type fracture, m;

x, y—the direction perpendicular and parallel to the fracture surface;

X—1D coordinate axis along the major axis of “I+II” type fracture;

z—the complex variable at any point in the “I+II” type fracture plain strain model;

$\bar{z}$—conjugate complex number of complex variable z;

Z(z),${{Z}_{\text{I}}}\left( z \right)$—Westergaard stress function of “I+II” type fracture and I type fracture;

$\widetilde{Z}\left( z \right)$,$\widetilde{{{Z}_{\text{I}}}}\left( z \right)$—integral of$Z\left( z \right)$, ${{Z}_{\text{I}}}\left( z \right)$ with respect to z;

α— wellbore approach angle, (°);

β, β₁, β₂—the angles between the line from z point to the midpoint and two endpoints of the long axis of the fracture with fracture long axis, (°);

Δu, Δv—relative displacement in tangential and normal directions of “I+II” type fracture, m;

Δu_p, Δv_p—relative displacement in tangential and normal directions of fracture with the influence of fluid pressure, m;

Δu_f, Δv_f—relative displacement in tangential and normal directions of “I+II” type fracture with influence of fluid pressure, m;

θ—NF approach angle, (°);

μ—fracture friction coefficient, dimensionless;

σ_a, σ_b—normal stress components of ground stress on the long and short axis of “I+II” type fracture, MPa;

σ_H, σ_h—the maximum and minimum horizontal principal stress, MPa;

σ_n—contact normal stress of matrix on fracture surface, MPa;

σ_v—vertical stress, MPa;

σ_x, σ_y—the normal stress component of ground stress in the x and y direction of the stress unit, MPa;

τ—the shear stress component of ground stress of the stress unit, MPa;

τ_c—wellbore shear stress, MPa;

τ_yx, τ_xy—the shear stress component of ground stress in the x and y direction of the stress unit, MPa;

υ—rock Poisson’s ratio, dimensionless;

υ°—rock Poisson’s ratio under plain strain condition, dimensionless;

Ψ—intermediate variables, (°).