During the hydraulic fracturing process for shale gas wells, the fracturing fluid is injected into the wellbore with large pump rate. The bottom-hole temperature will decrease sharply, which increasing the risk of casing failure. Based on the actual data of shale gas wells, the influence of rheological parameters of fracturing fluid on heat transfer coefficient is analyzed. A transient thermo-pressure coupling model of casing-cement sheath-formation combination is established. Sensitivity analysis is conducted for different pump rates, pressures and injection temperatures. The results indicate that:

1. The rheological parameters of fracturing fluid can affect the heat transfer coefficient between the fluid and borehole wall, which in turn affects the temperature distribution at the bottom of the wellbore.
2. The increasing of pump rate will drastically reduce the bottom-hole temperature, then thermal stress generating on the casing.
3. The higher the original reservoir temperature, the greater the influence of the temperature.
4. The casing stress decreases then increases with the increasing of fracturing pressure. Therefore, the rheological parameters of fracturing fluid should be chosen reasonably to reduce the heat transfer coefficient. What is more, it is crucial to choose the appropriate fracturing pump rate, fracturing pressure and injection temperature, as far as possible to reduce the bottom-hole temperature difference, which avoiding excessive casing stress to enhance the casing safety for fracturing shale gas wells.

Keywords: rheological parameter thermo-pressure coupling, transient temperature field, casing stress, hydraulic fracturing

The casing deformation problems are presented over 36 wells (among 112 horizontal wells by 2016) during fracturing processes in Weiyuan-Changning shale gas play in China.¹ Subsequent tools could be blocked in the wellbore.^2–4 The segments with serious deformation had to be abandoned before completing fracturing operations.⁵ The casing deformation issues make a great influence on production of shale gas wells.

Shale gas fracturing is carried out with large pump rate and high pressure. When the fracturing fluid is injected into the wellbore, it will have a great influence on the downhole temperature. Zhao et al.⁶ calculated the distribution of wellbore temperature field during fluid injection or fluid production by implicit difference model. Wang et al.,⁷ proposed a new method of heat transfer in unsteady state of wellbore and took a detailed calculation of the wellbore temperature field during fracturing process. Tang et al.⁸ established a mathematical model of heat conduction in porous media, giving the temperature analysis solutions of well in the case of injection and production operations through the mathematical method. Wang et al.⁹ considered the friction heat of fluid during fracturing process, and calculated the wellbore temperature field. He held the point that the friction heat generated on the fracturing fluid temperature could not be ignored. Satman et al.¹⁰ conducted an in-depth discussion of the wellbore heat transfer theory and clarified the essence of static and transient heat transfer. You et al.¹¹ established a fully implicit algorithm for the wellbore to gain the unsteady temperature field of surrounding rock. Cai et al.¹² established a coupled model to study the temperature field distribution during the production of fractured horizontal wells. To sum up, scholars have done a lot of research on the calculation of wellbore temperature field, but most of them only focused on the calculation of temperature field under the condition of tubing injection. Rare researches were conducted on transient temperature field and thermal stress during fracturing process for shale gas wells. Based on rheology and thermo-elasticity theory, a transient temperature-pressure coupling model of casing-cement sheath-formation combination is established in this paper. Considering the influence of fracturing fluid rheological parameters, emphasis is put on the influences of fracturing pump rate, injection temperature and fracturing pressure on casing stress. Corresponding countermeasures are proposed according to the law of the impact.

• Heat transfer coefficient calculation

The fracturing fluid used in the fracturing process is generally non-Newtonian fluid¹³, while its rheological property equation is:

$\tau \text{=}K{\gamma }^n$ ……….. (1)

During fracturing process, the pump rate is usually very large, being in turbulent flow state. The heat transfer coefficient between the fracturing fluid and the inner wall of the casing can be calculated using the Marshall model¹⁴ in Eq. (2-7):

$h=\frac{N_u{k}_m}{D}$ ………………… (2)

$N_u(S_t\text{)=}0.0107R_{eg}^{0.67}{P}_r^{0.33}$  ………3)

${\mu }_{w.app}=K{\left(\frac{3n+1}{4n}\right)}^n{\left(\frac{8v}{D}\right)}^{n-1}$ ………. (4)

$R_{eg}=\frac{{\rho }_a{D}_{eff}v}{{\mu }_{w.app}}$………… (5)

$P_r=\frac{{\mu }_{w.app}{C}_m}{k_m}$………… (6)

$v=\frac{q}{15\pi D^2}$ ………. (7)

 Finite element model establishment

After cementing slurry solidification, the casing, cement sheath and surrounding rock wall will be consolidated as a whole. Assuming that the contact part of the three cementing good, non-slip generated, three parts are uniform homogeneous thermo-elastic,^15,16 based on the thermo-elastic theory, a temperature - pressure coupling combined model is established, shown in Figure 1.

Since the axial dimension of the assembly is much larger than the radial dimension, the assembly can be simplified as a plane strain model.¹⁷ It means that the strain exists only in the radial direction but no strain in the axial direction.

Thermal conduction equations

The transient heat transfer equation is obtained as time-dependent temperature distribution by Eq. (8):

$k_i\left(\frac{{\partial }^{2}T}{\partial x^2}\text{+}\frac{{\partial }^{2}T}{\partial y^2}+\frac{{\partial }^{2}T}{\partial z^2}\right)\text{+}{Q}_i=C_{pi}{\rho }_i\frac{\partial T}{\partial t}$  ………… (8)

Initial condition:$T^{\text{i}}=T^{\text{i}}(x,y,z)|{}_{{}_{t=0}}$

Boundary conditions

Surface S1 boundary condition:$T^{\text{i}}=T^{\text{i}}(S,t)|{}_{S_1}$

Surface S2 boundary condition:$k_i\frac{\partial T^{\text{i}}}{\partial x}{n}_x+k_i\frac{\partial T^{\text{i}}}{\partial y}{n}_y+k_i\frac{\partial T^{\text{i}}}{\partial z}{n}_z=q_i(S,t)|{}_{S_2}$

Surface S3 boundary condition:$k_i\frac{\partial T^{\text{i}}}{\partial x}{n}_x+k_i\frac{\partial T^{\text{i}}}{\partial y}{n}_y+k_i\frac{\partial T^{\text{i}}}{\partial z}{n}_z=h\text{ (}{T}_a-T^{\text{i}}\text{)}|{}_{S_3}$

Based on the Galerkin's weak form of equivalent integrals, the heat conduction equations can be written as Eq. (A-3)

${\Pi }^{\text{i}}(T)={\displaystyle {\int }_{\Omega }k\nabla T^{\text{i}}\cdot \nabla (\delta T^{\text{i}})d\Omega }-{\displaystyle {\int }_{\Omega }\delta T^{\text{i}}(Q_i-{\rho }_i{C}_{pi}{\dot{T}}^i)d\Omega }-{\displaystyle {\int }_{S_2}{q}_i(\delta T^{\text{i}})d\text{S}-}{\displaystyle {\int }_{S_3}h\left(T_a-\frac{1}{2}{T}^{\text{i}}\right)(\delta T^{\text{i}})d\text{S}}$ (9)

The regions are divided into finite elements. Specific element nodes temperatures are used to represent the whole elements temperatures through the appropriate shape function. Finite element equations about nodes temperature are obtained using the boundary conditions. The nodes temperatures can be obtained by solving these finite equations, shown in Eq. (10-13).

$C\dot{T}+KT=P$…….. (10)

$C_{ij}={\displaystyle {\sum }_e{C}_{ij}^e}$…………… (11)

$K_{ij}={\displaystyle {\sum }_e{K}_{ij}^e}+{\displaystyle {\sum }_e{H}_{ij}^e}$………. (12)

$P_i={\displaystyle {\sum }_e{P}_{Qi}^e}+{\displaystyle {\sum }_e{P}_{qi}^e}+{\displaystyle {\sum }_e{P}_{Hi}^e}$…… (13)

 Thermal-pressure coupling equations

The basic equations of mechanics are derived from the theory of elasticity.

Balance equations:

$L^T\sigma +f=0$……… (14)

Geometric equations:

$\epsilon =Lu=LNa{}^e=B{a}^e$……….. (15)

Constitutive equations:

$\sigma =D(\epsilon -{\epsilon }_0)=DBa{}^e=S{a}^e$………. (16)

Boundary conditions in Sσ:

$\sigma n=T$ ……. (17)

Degree of freedom conversion equations

${\Pi }_p={\int }_{\Omega }\frac{1}{2}{\epsilon }^{T}D\epsilon d\Omega -{\int }_{\Omega }{\epsilon }_0^{T}D{\epsilon }_{0}d\Omega -{\int }_{\Omega }{u}^{T}fd\Omega -{\int }_{s_{\sigma }}{u}^{T}Tds$  ……….. (18)

There would be no free space to deform because the casing, cement sheath, and formation consolidated as a whole. Thermal strain would be generated when the temperature changed. The strains in transverse isotropic elastic body are:

${\epsilon }_0\text{=}{\left[\begin{array}{cccccc}\alpha \Delta T& \alpha \Delta T& \alpha \Delta T& 0& 0& 0\end{array}\right]}^T$……… (19)

The functional total of the minimum potential energy is:

${\Pi }_p={\displaystyle {\int }_{\Omega }\frac{1}{2}{\epsilon }^{\text{T}}D\epsilon \text{d}\Omega }-{\displaystyle {\int }_{\Omega }{\epsilon }_0^{\text{T}}D{\epsilon }_0\text{d}\Omega }-{\displaystyle {\int }_{\Omega }{u}^{\text{T}}f\text{d}\Omega }-{\displaystyle {\int }_{S_{\sigma }}{u}^{\text{T}}T\text{d}S}$……… (20)

Solved domains should be discrete firstly. The system potential energy can be represented by the node displacements through the elements total potential energies. According to the variation principle, one order variation of functional is set zero. The available value of the function is the displacement vector sought. Then finite element solution equations can be obtained as the following:

$Ka=P$……… (21)

$K={\displaystyle \sum _e{G}^T\left({\displaystyle {\int }_{{\Omega }_e}{B}^{T}DB\text{d}\Omega }\right)G}$………… (22)

$P=\underset{e}{{\displaystyle \sum G^T}}\left({\int }_{{\Omega }_e}{N}^{T}fd\Omega +{\int }_{S_{\sigma }^e}{N}^{T}Tds\right)+\underset{e}{{\displaystyle \sum {\int }_{{\Omega }_e}}}{B}^{T}D{\epsilon }_{0}d\Omega$……….. (23)

Since the thermal parameters do not change with time within the considered temperature range, it can be superimposed with the mechanical analysis step after the entire transient temperature field analysis is completed. In this way, the stress can be calculated under the condition of thermal-pressure coupling.

Parameter setting

According to the theory of elastic mechanics, stress concentration phenomenon occurs near the circular hole in an infinite plate,²⁰ but when the boundary size exceeds 5-6 times than that of well-hole diameter, the influence of stress concentration is very small. The material and thermal parameters are shown in Table 1.

The heat convection coefficient

The rheological parameters and pump rate of fracturing fluid will affect the heat transfer coefficient²¹. Using equations (2-7), the variation law of heat transfer coefficient under different rheological parameters and pump rates can be calculated, shown in Figure 2.

As shown in Figure 2 when the rheological parameters n or K keeps constant, the heat transfer coefficient increases with the increasing of the pump rate. When the pump rate stays the same, the heat transfer coefficient firstly decreases dramatically then slowly with the increasing of rheological parameters. Therefore, the small heat transfer coefficient can be gained by appropriate rheological parameters and low pump rate.

The transient temperature of fracturing fluid

Setting the rheological parameters of fracturing fluid K=0.08, n=0.8, pump rate Q= 0/3/6/20m³/min, then calculating the corresponding heat transfer coefficient. The results are shown in Table 2.

From the heat transfer coefficient calculated above, the transient bottom-hole temperature change can be calculated by establishing the heat conduction model (10) to calculate the reservoir temperature at 100 and 150, respectively. The results are shown in Figure 3.

It can be seen from Figure 3 that when the pump rate of fracturing fluid is zero, the bottom hole temperature changes more slowly, and when the pump rate increases, the bottom hole temperature will decrease drastically. With a certain pump rate, the higher the temperature at the bottom of the well, the higher the temperature decreases. For example, the maximum temperature difference at the bottomhole temperature of 100°C is close to 80°C and the maximum temperature difference at the bottomhole temperature of 150°C is close to 120°C. Therefore, using large pump rate parameter will have a larger temperature difference.

The influence of pump rate on casing stress

Based on the basic equations of temperature-pressure coupling model, the transient temperature field changes calculated above are substituted into the stress model of the composite model. The temperature-pressure coupling calculation is carried out to calculate the stress variation of the inner wall of the casing body under different pump rates. The results are shown in Figure 4.

It can be seen from the Figure 4 that the fracturing pump rate will exert a certain influence on the casing stress. With the increase of rate, the bottom hole temperature difference increases, then the casing stress increases. At the same time, the higher the formation temperature, the greater the temperature difference of the bottom hole and the larger the casing stress will be. This indicates that during fracturing process, large pump rate will lead to dramatic changes in temperature, thereby increasing the casing stress.

The influence of fluid injection temperature on casing stress

The fracturing fluid temperature is crucial to guarantee the casing safety. The fluid temperature in the model are set as 0/20/40/60/80/100℃. The regularity of the influences of different fluid temperature on casing stress are investigated and the results are shown in Figure 5.

The transient internal casing temperature drops dramatically in the first 5 minutes, then keeps almost the same with the lapse of time. The lower the fluid temperature is, the greater the bottom casing temperature reduction. At the same time, the casing stress will increase with the decreasing of fluid temperature. In view of this, the fracturing fluid temperature is important for the safety of casing. During fracturing, the high temperature fracturing fluid can be adopted to reduce the risk of casing failure.

The influence of fracturing pressure on casing stressThe fracturing fluid with high pressure is adopted to fracture the shale formation. The larger the tectonic stress is, the higher the pressures are. Meanwhile, the high fracturing pressure can pose a great potential challenge to casing deformation. The influences of different pressures on casing stress were simulated. The results are shown in Figure 6.

It can be seen from Figure 6 that the casing stress decreases firstly, then increases with the increasing of fracturing pressure. The maximum stress changes from 90° to 0°. In order to fracture the shale formation, the fracturing pressure mush be high enough. But too large pressure will increase casing stress. That can tell us that some measures needs to be adopted to decrease the fracturing stress of shale formation. Only in this way, can the casing be a safe place.

Figure 1 Transient thermal-pressure coupling model of casing-cement sheath-formation.

Figure 2 Effect of defection angle and phase angle on casing stress.

Figure 3 The influence of the pump rate on the downhole transient temperature.

Figure 4 The influence of the pump rate on casing stress.

Figure 5 The influence of temperature of the fracturing fluid.

Figure 6 The influence of the fracturing pressure on casing stress.

During the hydraulic fracturing process for shale gas wells, temperature drops dramatically in the bottom of the well due to the large pump rate. To evaluate the casing stress during volume fracturing, the heat transfer coefficient was calculated. A finite element model considering the transient thermal-pressure coupling was established. Sensitivity analyses of the influences of pump rate, injection temperature, and fracturing pressure on casing stress were investigated. Conclusions can be drawn as the following:

1. A small convective heat transfer coefficient could be gained by appropriate rheological parameters and low pump rate. Downhole temperature decreased quickly with the increasing of the pump rate.
2. Casing stress increased with the increasing of pump rate. The higher the initial reservoir temperature, the greater the casing temperature decreasing.
3. The bottom-hole temperature decreased with the increasing of fluid injection temperature, then decreasing the casing stress.
4. The casing stress firstly decreased, and then increased with the increasing of fracturing pressure.

During fracturing operation, the warm injection fracturing fluid, appropriate pump rate and fracturing pressure could be adopted to reduce the casing stress.

K is the consistency factor, Pa·s-n;
n is the fluidity index, n<1 presented the pseudo-plastic fluid, n>1 presented the expansive fluid, for the fracturing fluid n<1, meaning the pseudo-plastic fluid;
γ is shear Rate, s-1.
h is the heat transfer coefficient.
W/(m2·°C)
Nu is the Nusselt number
Pr is the Prandtl number
Reg id 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/m3.
v is the fluid velocity.
Q is the fracturing pump rate, m3/min,
km is the coefficient of heat conductivity.
W/(m).
Cm is the fluid specific heat capacity, J/(kg).
ki is the material coefficient of heat conductivity, W/(m·°C).
Cpi is the material specific heat capacity, J/(kg·°C).
ρi is the material density, kg/m3.
Qi is the heat source density, W/m3.
qi is the boundary heat source density.
W/m2
Ti is the temperature, °C
Ta is the fracturing fluid temperature, °C
t is the time ,s, nx, ny, nz are the normal vector of S3 in x, y, z directions.
i=1, 2, 3 are represented the casing, cement sheath, and formation respectively.
α is the coefficient of thermal expansion, °C-1,
ΔT is the temperature difference.
L is the differential operator.
σ is the stress matrix.
ε is the strain matrix.
ε0 is the thermal strain matrix.
f is the body force array.
u is the displacement array.
N is the shape function matrix.
ae is the element displacement array.
D is the elastic matrix.
B is the strain matrix.
S is the stress matrix.
T is the boundary surface force array.
n is the surface normal vector.
a is the structure displacement array
G is the transformation matrix of degree of freedom between element nodes and structure nodes.
K is the Structure overall stiffness matrix.
P is the structural equivalent node load array.
C is the heat capacity matrix, Cij is any element of C.

This work was supported by the National Natural Science Foundation of China, Project No. 51674272, National Science and Technology Department of China, Project No. 2017ZX05009-003, the China University of Petroleum (Beijing), Project No. 2462015QZDX05, and the State Key Laboratory of Shale Oil and Gas Enrichment Mechanisms and Effective Development, Project “Research on casing loading spectrum under multi-factor coupling effect in shale gas wells”

The author declares there is no conflict of interest.

1. Chen Z, Shi L, Xiang D. Mechanism of casing deformation in the Changning–Weiyuan national shale gas demonstration area and countermeasures. Natural Gas Industry B. 2017;4(1):1–6.
2. Liu K, Wang Y, Gao D, et al. Effects of hydraulic fracturing on horizontal wellbore for shale gas. Acta Petrolei Sinica. 2016;37(3):406–414.
3. Jiang K, Li Q, Chen Y, et al. Influence of cementing quality on casing failure in horizontal shale gas wells. Natural Gas Industry. 2015;(12):77–82.
4. Tian Z, Shi L, Qiao L. Problems in the wellbore integrity of a shale gas horizontal well and corresponding countermeasures. Natural Gas Industry. 2015;(09):522–529.
5. Dong W, Shen R, Liang Q, et al. Calculation and analysis of casing thermal stress during stimulated reservoir volume fracturing. Fault-Block Oil & Gas Field. 2016;(05):673–675.
6. Zhao J, Ren S. Numerical calculation of temperature distribution of liquid in wellbore. Oil Drilling and Production Technology. 1986;(03):49–57.
7. Wang H, Li P. Numerical Calculation Method of Wellbore Temperature in Hydraulic Fracturing Process. Acta Petrolei Sinica. 1987;(02):91–99.
8. Lu D, Zeng Y, Guo Y. Analytical solution to wellbore and formation temperature in multilayered stratum. Journal of Hydrodynamics (A). 2002;(03):382–390.
9. Wang L, Li Q. Research on Theory and Calculation Method of Temperature Distribution in Vertical Wellbore Fracturing. Science Technology and Engineering. 2013;(35):10466–10469.
10. Satman A, Tureyen OI. Geothermal wellbore heat transmission: Stabilization times for “static” and “transient” wellbore temperature profiles. Geothermics. 2016;64:482–489.
11. You J, Rahnema H, Mcmillan M D. Numerical modeling of unsteady-state wellbore heat transmission. Journal of Natural Gas Science and Engineering. 2016;34:1062–1076.
12. Cai J, Duan Y. Study on temperature distribution along wellbore of fracturing horizontal wells in oil reservoir. Petroleum. 2015;1(4):358–365.
13. Li D, Jia J, Teng F, et al. Calculation of bottom hole pressure during fracturing. Fault-Block Oil and Gas Field. 2013;(03):384–387.
14. Fan H. Practical Drilling Fluid Mechanics. Beijing: Petroleum Industry Press; 2014.
15. Yin Y, Chen Z, Li P. Analytical Stress Distribution of Casing - Cement Sheath- Formation. Chinese Journal of Theoretical and Applied Mechanics. 2006;(06):835–842.
16. Li J, Chen M, Zhang H. Effects of cement sheath elastic modulus on casing external collapse load. Journal of the University of Petroleum. 2005;(06):41–44.
17. Li Z. Finite element analyses of casing shear failure in Member II of Nenjiang Formation in Daqing Oilfield. J Fault-Block Oil and Gas Field. 2015;(06):812–815.
18. Wang M. Finite Element Method [M]. China: Tsinghua University Press; 2015.
19. Zhang L. Yan X, Yang X. Evaluation of wellbore integrity for HTHP gas wells under solid-temperature coupling using a new analytical model. Journal of Natural Gas Science and Engineering. 2015;25:347–358.
20. Li X, Liu J, Dong S. Design and Application of Local Double Layer Bushing in High Pressure Anomaly. Fault-Block Oil and Gas Field. 1997;5:42–46.
21. Zhang F, Shen Y, Ren T. Application of Similar Reynolds Number Method in Frictional Resistance Test of Slip Water Fracturing String. Science Technology and Engineering. 2016;(25):86–90.