Next Article in Journal
Social Networking Service as a Marketing Technology Tool and Sustainable Business in the Lodging Industry: Investigating the Difference across Older and Younger Age Groups among Tourists
Next Article in Special Issue
Assessment of Emission Reduction and Meteorological Change in PM2.5 and Transport Flux in Typical Cities Cluster during 2013–2017
Previous Article in Journal
AI Technology and Online Purchase Intention: Structural Equation Model Based on Perceived Value
Previous Article in Special Issue
Day-Ahead Wind Power Forecasting Based on Wind Load Data Using Hybrid Optimization Algorithm
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Investigation on Wellbore Temperature Prediction during the CO2 Fracturing in Horizontal Wells

1
State Key Laboratory of Petroleum Resources and Engineering, China University of Petroleum, Beijing 102249, China
2
National Computer Network Emergency Response Technical Team/Coordination Center of China, Beijing 100029, China
3
State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Chengdu 610500, China
*
Author to whom correspondence should be addressed.
Sustainability 2021, 13(10), 5672; https://doi.org/10.3390/su13105672
Submission received: 31 March 2021 / Revised: 2 May 2021 / Accepted: 5 May 2021 / Published: 18 May 2021
(This article belongs to the Special Issue Environmental and Economic Analysis of Low-Carbon Energy Technologies)

Abstract

:
A novel model is established to predict the temperature field in the horizontal wellbore during CO2 fracturing. The pressure work and viscous dissipation are considered, and the transient energy, mass and momentum equations as well as the CO2 physical properties are solved fully coupled. The model passes the convergence test and is verified through a comparison using the COMSOL software. Then, a sensitivity analysis is performed to study the effects of the treating parameters. Results illustrate that the relationship between the injection rate and the stable bottom-hole temperature (hereinafter referred to as BHT) is non-monotonic, which is different from the hydraulic fracturing. The existence of the horizontal section will increase the BHT at 2 m3/min condition but reduce the BHT at 10 m3/min condition. The problem of high wellbore friction can be alleviated through tube size enhancement, and the ultimate injection rate allowed increased from 2.7 m3/min to 29.6 m3/min when the tube diameter increased from 50.3 mm to 100.3 mm. Additionally, the open-hole completion method of the horizontal section can increase the BHT to 2.7 °C but reduce the near formation temperature to 24.5 °C compared with the casing completion method.

1. Introduction

Shale gas is one of the most important unconventional oil and gas resources [1,2]. Slick-water fracturing and horizontal well drilling are the basic technologies for the development of the shale gas reservoirs [3]. However, the usage of the water-based fracturing fluid can generate formation damage and consume a large amount of water resources. CO2 fracturing can simultaneously avoid the above disadvantages and achieve the effects of carbon sequestration (CCS) [4,5] and gas recovery enhancement (EGR) [6,7,8,9], which is a potential alternative technology to the slick-water fracturing. Unlike the water-based fracturing fluid, the physical properties of CO2 vary significantly with temperature and pressure [10,11], which can affect the sand-carrying capacity and the performance of the stimulation. Additionally, it is impractical to modify the CO2 temperature during the fracturing process based on the temperature monitor data because of the expense and time consumed. Therefore, accurately predicting the temperature and pressure of the horizontal well CO2 fracturing process is necessary in terms of both the environmental and economic aspects.
The study of the wellbore temperature field can be traced back to 1937 [12]. After a long-term development, the wellbore temperature prediction methods can be divided into three categories. The first category is the quasi-steady numerical method represented by the Ramey’s model [13], which assumes that the fluid and wellbore zones are steady-state heat transfer, and the formation zone is transient heat transfer. The quasi-steady numerical method is applicable to production and injection (flooding process) wells, and many researchers have extended the model to various wellbore conditions [14,15,16,17,18,19,20,21,22,23,24,25,26,27]. However, it is not suitable for the fracturing process [28]. The second category is the transient numerical method represented by the Eickmeier’s model [29], which assumes that all the calculation zones are transient heat transfer. Considering the transient heat transfer condition, this kind of method is widely used in the temperature prediction of the hydraulic fracturing and acidification treatments [30,31,32,33], in which the injection rates are typically high. The third category is the analytic method. The assumptions of the steady [34] and transient [35,36,37,38,39,40,41] heat transfers are both used in the analytic method. However, although the analytic method is faster and more stable than the numerical method, it requires further simplified conditions to obtain the solution, which restricts its applications. The three kinds of wellbore temperature prediction methods discussed above are all applied to the CO2 injection process.
As for the quasi-steady numerical method, Lu and Connell [42] used the quasi-steady heat transfer equation and steady-state mass and momentum equations to simulate the non-isothermal flow of a vertical well during the injection of a multicomponent mixture containing CO2. The physical properties of the mixture were calculated by the Pen–Robinson Equations of State (hereinafter referred to as EoS). Wang [43] combined the Span–Wagner model [44], Fenghour model [45], and Vesovic model [46] with the quasi-steady temperature equation and the steady mass and momentum equations to analyze the wellbore temperature and pressure distribution in a vertical well during supercritical CO2 drilling. Dou [47] applied a similar method to simulate the wellbore temperature field of the CO2 flooding process and compared it with the measured data to verify the accuracy of the model. Li [48] expanded the application of Wang and Dou’s work to CO2 sequestration and CO2 fracturing. Sun [49] applied the quasi-steady numerical method to the CO2 injection heat recovery process of a geothermal U-pipe system and simulated the temperature and pressure distribution of CO2 in a closed cycle system.
Although the quasi-steady numerical method performed well in low injection rate (typically < 0.1 m3/min) conditions such as CO2 flooding, CO2 sequestration, and CO2 heat recovery, it is not suitable for a short time and fast injection processes such as fracturing owing to the adaptation of the steady heat transfer assumption in the fluid energy equation. As stated above, the transient numerical method can effectively solve this problem. Pan [50] developed a vertical wellbore simulator following the assumption of the transient non-isothermal flow of CO2 brine mixture and using the drift flux model (DFM) to describe momentum conservation. Then, Pan [51,52] combined their wellbore model with the multiphase flow to analyze the non-isothermal flow of the mixture in the formation. In their wellbore heat transfer model, the fluid heat transfer equation is transient. By contrast, the heat transfer of the wellbore zone is characterized by the comprehensive heat transfer coefficient, which is derived on the basis of the steady-state heat transfer assumption. Guo [53,54] and Gong [55,56] considered the influence of the frictional heat and the Joule–Thomson effect based on the Eickmeier’s model to establish a wellbore temperature field model for CO2 fracturing in a vertical well, in which the CO2 physical property model and the transient continuity and momentum equation are coupled.
Compared with the numerical method, research scarcely focuses on the analytic method. Singhe [57] considered the Joule–Thomson effect on the basis of Hagoort’s model [40] to build an analytical wellbore temperature model for CO2 injection wells and verified it through the data of the Ketzin area in Germany; however, the analytical solution is based on the quasi-steady heat transfer assumption, which is also not suitable for the fracturing process. Additionally, some researchers [58,59] used 2D turbulence models such as the k-ε model to investigate the flow and heat transfer characteristics of the CO2 injection process in vertical wells. However, considering the computation and the convergence requirements, utilizing it is difficult in the simulation of the horizontal wells. In general, although the above literature brought important achievements for the wellbore temperature prediction of CO2 injection process, it has focused on the vertical wells. Moreover, an efficient wellbore temperature prediction model remains lacking for horizontal wells, which are extensively used in unconventional oil and gas reservoirs.
In this article, a transient fully coupled model of wellbore temperature for CO2 fracturing in horizontal wells is established, which combines the 1D non-isothermal flow in the axial direction with the 1D transient heat transfer in the radial direction. The effects of the pressure work and viscous dissipation are considered in the energy equation, and the transient mass and momentum equations as well as the CO2 physical property models are solved simultaneously. The 1D + 1D model is verified through a convergence test and a comparison using the COMSOL software. Finally, the parameter analysis is performed to provide guidance for the treatment of the CO2 fracturing in horizontal wells.

2. Theory

2.1. Physical Model

Figure 1 illustrates the physical model. The wellbore is divided into vertical, curved, and horizontal sections. In the vertical and curved sections, the wellbore structure includes the tubing, annulus, casing, and the cement sheath. In the horizontal section, as the open-hole completion method is considered here, the wellbore structure only includes tubing and annulus.

2.2. Basic Assumption

The horizontal well CO2 fracturing has three characteristics. First, the aspect ratio of the wellbore of oil and gas wells is usually very large (1 × 103 to 1 × 104) given the enormous gap between the length of the axial direction and the width of the radial direction. Second, the injection rate is typically above 3 m3/min, which results in a high Reynolds number (Re > 1 × 106) and a turbulent flow during fracturing. Last but not least, the L-shaped structure of the horizontal wellbore cannot be simplified through the method of axisymmetric modeling as the vertical well in 2D or 3D space. All of the above necessitate the meshes in the 2D, or the 3D Computational Fluid Dynamics (hereinafter referred to as CFD) model becomes a huge number to obtain stability and convergence. An example is the k-ε model. The dimensionless wall distance function y+ of the first boundary layer grids should be 30–200 [58,60], and the dimensionless wall distance function can be calculated as y+ = vΔy/μ, where y is the thickness of the first boundary layer (please notice that parameters that are not further explained in the main text hereinafter are listed at the Nomenclature section). For CO2 fracturing, this process requires the thickness of the first boundary layer to be less than 0.001 m (as shown in Figure 2).
To solve the above problems, considering the following two points: (1) the radial dimension is pretty small compared with the axial dimension, (2) the effect of the axial heat conduction on the wellbore temperature can be ignored [53,54], the following assumptions are used:
  • The CO2 zone is considered a one-dimensional axial non-isothermal flow;
  • The wellbore and formation zones are considered a one-dimensional radial transient heat transfer;
  • The heterogeneity and anisotropy of the formation are ignored;
  • At the initial moment, the wellbore is full of fluid, and the fluid is stationary.

2.3. Mathematical Model

The calculation zone includes the CO2, wellbore, and formation zones. The CO2 zone is considered a non-isothermal flow, and the continuity equation, momentum equation, as well as the energy equation must be solved simultaneously. The wellbore zone is considered a fluid natural convective heat transfer, and the formation zone is considered a solid heat transfer process, for which the energy equation must be solved.
Continuity equation of the CO2 zone:
v ρ f z + ρ f v z + ρ f t = 0
Momentum equation of the CO2 zone:
p z = ρ f g sin θ 1 2 f D ρ f v 2 d ti ρ f v t
Energy equation of the CO2 zone:
A ρ f C p f T f t + A ρ f C p f v T f z = Q wall + Q fric + Q pres .
where fD is the friction coefficient. The fitting formula (R2 > 0.999) [61] of the CO2 fracturing wellbore friction plate [62] is used here to obtain the expression of fD.
P f = 0.156 Q inj 2 + 1.639 Q inj 0.228 d ti = 50.3 mm P f = 0.037 Q inj 2 + 0.777 Q inj 0.169 d ti = 62 mm P f = 0.017 Q inj 2 + 0.268 Q inj 0.041 d ti = 76 mm P f = 0.003 Q inj 2 + 0.089 Q inj 0.021 d ti = 100.3 mm ,
f D = 2 P f d ti / 100 ρ v 2 .
Qwall, Qfric and Qpres are the heat flows of the unit axial length generated by the heat transfer, viscous dissipation, and pressure work, respectively, W/(m2·K). They are defined as follows:
Q wall = π d ti h t T f T t ,
Q fric = f D v ρ f A v 2 / 2 d ti ,
Q pres = A α p T f p / t + v p / z .
where ht is the heat transfer coefficient of the forced heat convection. The Gnielinski correlation, which is valid for 0.5 ≤ Pr ≤ 2000 and 3000 ≤ Re ≤ 5 × 106, is applied here to calculate ht [63]:
h t = λ f f D / 8 Re 1000 Pr 2 r ti 1 + 12.7 f D / 8 1 / 2 Pr 2 / 3 1 .
Pr is the Prandtl number, defined as Pr = Cpμ/λ. Unless otherwise specified, the heat flow mentioned in the following refers to the heat flow of the unit axial length.

2.3.1. Energy Equation of the Wellbore Zone

The wellbore zone includes the tubing, annulus, casing, and the cement sheath. Except for the annulus, the remaining regions are all solid, and only the heat conduction requires consideration. For the annulus, the internal fluid is in a state of natural convection, which is characterized here by the natural convection heat transfer coefficient. The overall governing equation of the wellbore zone is as follows:
1 r r λ r T r = ρ C p T t .  
In the formula, ρ, Cp and λ take the values of different materials according to different spatial positions, and the temperature at the interface of different zones is continuous. The specific form of the governing equation of each part is different, and the wellbore structure in the horizontal section is also different from the vertical and curved sections. The governing equations for each part of the wellbore are given below.

Vertical Section (Curved Section)

Figure 3 presents the radial grids of the vertical section (curved section), in which the white point with blue edge indicates the location of the various radius of the tubing, casing, cement sheath, and formation, while the red point with black edge indicates the location of the calculation temperature.
The energy equation of the tubing:
2 r to h an T an T t r to 2 r ti 2 2 r ti h t T t T f r to 2 r ti 2 = ρ t C p t T t t .  
The energy equation of the annulus:
2 r ci h an T c T an r ci 2 r to 2 2 r to h an T an T t r ci 2 r to 2 = ρ an C p an T an t .  
The energy equation of the casing:
4 r co λ c , h T h T c r co 2 r ci 2 r h r ci 2 r ci h an T c T an r co 2 r ci 2 = ρ c C p c T c t .
The energy equation of the cement sheath:
4 r h λ h , r T r 1 T h r h 2 r co 2 r r 1 r co 4 r co λ c , h T h T c r h 2 r co 2 r h r ci = ρ h C p h T h t .  
where han is the natural convection heat transfer coefficient in annulus, W/(m2·°C), defined as
h an = 0.049 Gr Pr 1 / 3 Pr 0.074 λ an r to ln r ci / r to .
Gr is the Grashof number of the annulus fluid, defined as
Gr = r ci r to 3 g ρ an 2 β an T f T c μ an 2 .  
where λc,h and λh,r are the composite thermal conductivity, which can have the following definition based on Fourier law:
λ c , h = λ c λ h ln r co + r h r co + r ci λ c ln r co + r h 2 r co + λ h ln 2 r co r co + r ci ,  
λ h , r = λ h λ r ln r h + r r 1 r h + r co λ h ln r h + r r 1 2 r h + λ r ln 2 r h r h + r co .  

Horizontal Section

Figure 4 illustrates the radial grids of the horizontal section. The implications of the two kinds of points are the same as in Figure 3.
The wellbore in the horizontal section only includes the tubing and the annulus. The energy governing equation of the tubing is the same as the vertical (curved) section. The annulus energy governing equation is as follows:
2 r h h an T r 1 T an r h 2 r to 2 2 r ti h an T an T t r h 2 r to 2 = ρ an C p an T an t .

2.3.2. Energy Equation of the Formation Zone

No flow occurs in the formation. Accordingly, the energy governing equation in the formation can be simplified to the heat conduction equation:
4 r r i λ r T r i + 1 T r i r r i 2 r r i 1 2 r r i + 1 r r i 1 4 r r i 1 λ r T r i T r i 1 r r i 2 r r i 1 2 r r i r r i 2 = ρ r C p r T r i t .
Owing to the difference between the radial grids of the vertical (curved) and horizontal sections, the equations at the boundary of the formation and the wellbore are different.
The first layer grids of the vertical (curved) section:
4 r r 1 λ r T r 2 , j T r 1 , j r r 1 2 r h 2 r r 2 r h 4 r h λ h , r T r 1 , j T h , j r r 1 2 r h 2 r r 1 r co = ρ r C p r T r 1 t .  
The first layer grids of the horizontal section:
4 r r 1 λ r T r 2 , j T r 1 , j r r 1 2 r h 2 r r 2 r h 2 r h h an , j T r 1 , j T an , j r r 1 2 r h 2 = ρ r C p r T r 1 t .  

2.4. Initial and Boundary Conditions

2.4.1. Initial Conditions

The initial wellbore temperature is calculated from the original geothermal gradient. As for the pressure, considering that the wellbore is assumed to be full of CO2, the initial wellbore hydrostatic pressure can be obtained through an iterative method combined with the CO2 density model. The governing equations for the initial conditions are as follows:
T t = 0 = T surf + g G z p t = 0 = 0 z ρ T , p g d z

2.4.2. Boundary Conditions

For the flow field in the CO2 zone, a given injection rate at the wellhead is used as the velocity boundary, and a pressure boundary is required at the bottom of the well. The current study uses the fracture closure pressure. For the temperature field of the CO2 zone, the injection temperature at the wellhead is given as the inlet boundary, and the bottom-hole of the well is insulated. However, the thermal convection is considered at the bottom-hole. The governing equations for the boundary conditions are as follows:
v z = 0 = Q inj / π r ti 2 p z = z bh = p w T z = 0 = T inj T z z = z bh = 0 .  

2.5. Solving Method

In the 1D + 1D transient fully coupled model, the continuity, momentum, and energy equations are coupled together through the interaction among pressure, temperature, and CO2 physical properties to form a transient nonlinear system, which requires an iteration solving method. The Span–Wagner EoS [44], Fenghour model [45] and Vesovic model [46] are used here as the physical properties models, which are proven to be highly accurate [10,11]. Two iteration loops are set up during the solution process: the pressure-speed and temperature iterations. Figure 5 shows the model solution flowchart.

3. Validation

3.1. Test of the Convergence

This paper used the iterative method to solve the velocity, pressure, and temperature. To test the convergence of the numerical algorithm, the changes in the calculation results must be analyzed during the iteration. Given that the velocity-pressure iteration is nested in the temperature iteration, the convergence test was performed on the temperature calculation results. Table 1 details the parameters of the benchmark case.
The convergence condition of the temperature calculation is defined as
T f k T f k 1 2 0.01
where ‖X2 represents the 2-norm of vector X. Figure 6a presents the T f k T f k 1 2 with various iteration steps and the number of iterations used to reach the convergence condition at different time steps.
Figure 6b indicates that as the number of iteration steps increased, the error continued to decrease; as the time step increased, the number of iteration steps used to achieve convergence condition also continued to decrease. Thus, the numerical solution method used in this paper met the requirements of convergence.

3.2. Comparison of the Non-Isothermal Flow Simulation in a Horizontal Tube with COMSOL

Currently, measured data of CO2 fracturing are lacking in horizontal wells. Accordingly, a comparison method with commercial software was used to verify the calculation results of the 1D + 1D transient fully coupled model in the horizontal section. COMSOL is a multi-physics coupled finite element software. The non-isothermal tube flow module can simulate thermal convection and heat conduction in the tube. However, in COMSOL, the wellbore wall layer (tubing, annulus, casing, cement sheath, etc.) notably adopts the steady-state heat transfer assumption and does not consider the natural convection in the annulus [64].
The calculation parameters were based on Table 1. The vertical section was ignored, and the horizontal section was extended to 1000 m. Simultaneously, the initial values were set as the parameters of the depth of 1600 m, implying that the initial temperature was 68 °C, and the initial pressure was 32 MPa, based on Table 1. As the horizontal wellbore was completed with the open-hole completion method, the casing and cement sheath were removed. The natural convection in the annulus was also ignored in the 1D + 1D fully coupled model for comparative analysis.
Figure 7a shows the variation of the BHT with time. Apparently, after ignoring the natural convection in the annulus, the results of the 1D + 1D transient fully coupled model were close to the results of COMSOL. However, at the early stage, the change rate of BHT calculated by COMSOL was faster than that of the 1D + 1D transient fully coupled model because of the steady-state wall heat transfer assumption used in COMSOL. The latter ignored the temperature stabilization process of each wall surface of the wellbore and underestimated the CO2 temperature. However, due to the large injection rate and the neglection of the natural convection, the effect of the wall heat transfer on the CO2 temperature was small; the temperature difference between the two was only 0.2 °C after 120 min. Considering that the treatment of wall heat transfer was the main difference between COMSOL and the 1D + 1D transient fully coupled model, another comparison was performed in which the wall heat transfer term was removed both from COMSOL and the 1D + 1D transient fully coupled model. As shown in Figure 7b, the calculation results of the two were highly consistent after neglecting the effect of wall heat transfer, and the deviation did not exceed 0.03 °C. This finding validated the model accuracy in the horizontal section.

4. Results and Discussions

This section presents an analysis of the effects of the treating parameters using the 1D + 1D fully coupled model. Table 1 details the parameters of the benchmark case.

4.1. Analysis of the Injection Rate

The cases were considered when the injection rates were 2 m3/min, 6 m3/min, and 10 m3/min, based on the benchmark case.

4.1.1. Analysis of the Time Domain

Figure 8 shows the variation of the BHT with time under various injection rates. When the injection time was less than 50 min, an increase in the injection rate speeded up the decrease in the temperature. This finding was due to the change in the convective heat transfer efficiency with the CO2 velocity. Figure 9 shows the average CO2 velocity under different injection rates. As the size of the tubing was unchanged, the velocity was basically linear with the injection rate. However, faster convective heat transfer will also hasten the stability of the CO2 temperature. Using 0.03 °C/min as the criterion of the temperature stability, Figure 10 shows that the temperatures of the cases of 2 m3/min, 4 m3/min, and 6 m3/min did not reach the stable state after injection for 120 min. The stable temperature of CO2 was related to the heat transfer mechanism determined by the combination of the pressure work, viscous dissipation, and wall heat transfer. The analysis of space domain will illustrate this process.
Figure 11 shows the variation of the wellhead pressure (hereinafter referred to as WHP) with time under various injection rates. Given that no flow occurred at the initial time, the initial WHP was a constant with the value of 17.5 MPa. When fracturing started, CO2 flowed into the wellbore and generated friction. To keep the BHP constant, the WHP increased. The enhancement was determined by the friction, which increased significantly as the injection rate increased. Figure 12 shows that the wellbore friction reached 67.7 MPa at 10 m3/min. In addition, Figure 11 shows that the pressure increased to a maximum within a short period of time, and then the pressure decreased slightly. This phenomenon was caused by the variation of the CO2 density, which was affected by the CO2 temperature. As shown in Figure 13, after 120 min, the CO2 density at the bottom-hole increased by 185 kg/m3 compared with that after 1 min.

4.1.2. Analysis of the Space Domain

Figure 14 shows the temperature profiles of CO2 with various injection rates after 120 min. With the increase in the measured depth, the axial direction and the wellbore structure changed. This was divided into three sections: vertical, curving, and horizontal sections. In the vertical section, the temperature of CO2 decreased first and then increased with the increase in the injection rate. This phenomenon was caused by the heating effect of the viscous dissipation. As shown in Figure 15a, the heat flow generated by the viscous dissipation was only 120 W/m at 2 m3/min but reached 4796 W/m at 10 m3/min. In the curved section, the well inclination angle slowly changed from 0° to 90°, and the gradient of the hydrostatic pressure gradually disappeared. Thus, the total pressure drop in the curved section was significantly increased, which resulted in an increase in the cooling effect of pressure work. As shown in Figure 15b, from the bottom of the vertical section to the bottom of the curving section, the variations of the heat flow of the pressure work were as follows:
2 m 3 / min : 196 W / m Curved section 116 W / m 6 m 3 / min : 317 W / m Curved section 1202 W / m 9 m 3 / min : 2228 W / m Curved section 3773 W / m
Therefore, in the curved section, the cooling effect of pressure work was obviously enhanced, leading to a decrease in the temperature gradients of all three cases. In the horizontal section, the absolute value of the Qpres under the injection rate of 10 m3/min continued to increase, whereas the Qpres under the injection rate of 2 m3/min did not change much. Simultaneously, as shown in Figure 15c, compared with the vertical and curved sections, the Qwall in the horizontal section increased significantly, which was due to the open-hole completion method. Section 4.4 discusses the effects of horizontal completions. Generally, with the combined effects of the pressure work and wall heat transfer, the temperature gradient in the horizontal section was higher than that of the vertical section under the 2 m3/min case but lower than that of the vertical section under the 10 m3/min case. The 6 m3/min case fell in between the other two cases. Therefore, compared with the vertical section, the existence of the horizontal section increased the BHT under a low injection rate but reduced the BHT under a high injection rate.
Figure 16 shows the pressure profiles of CO2 with various injection rates after 120 min. CO2 pressure was determined by the hydrostatic pressure and wellbore friction. In the vertical section, when the injection rate was 2 m3/min, the hydrostatic column pressure was higher than the wellbore friction, resulting in a positive pressure gradient. Additionally, the cases of 6 m3/min and 10 m3/min were contrary to the 2 m3/min case. In the horizontal section, as the hydrostatic pressure no longer increased, the CO2 pressure decreased along the axial direction under any injection rates. The curved section was the transition between the vertical and horizontal sections, during which the hydrostatic column pressure gradient slowly decreased to zero.

4.2. Analysis of the Injection Temperature

The cases where the injection temperatures were −20 °C, 0 °C and 20 °C were considered on the bases of the benchmark case.

4.2.1. Analysis of the Time Domain

Figure 17 shows the variation of the BHT with time under various injection temperatures. Evidently, the increase in the injection temperature effectively enhanced the BHT. When the injection temperature increased to 20 °C, the BHT maintained above the supercritical temperature during the entire fracturing process. Figure 18 shows the variation of the WHP with time under various injection temperatures. As the lower temperatures increased the CO2 density, which enhanced the hydrostatic pressure, the WHP of the −20 °C case was 1.8 MPa lower than that of the 20 °C case.

4.2.2. Analysis of the Space Domain

Figure 19 shows the temperature profiles of CO2 with various injection temperatures after 120 min. When the injection temperature increased from −20 °C to 20 °C, the temperature of CO2 increased significantly. However, the temperature gradient of CO2 decreased continuously, especially in the horizontal section. Figure 20 indicates that when the temperature increased from −20 °C to 20 °C, the temperature difference between the wellhead and the bottom-hole decreased by 5.2 °C. This phenomenon was caused by the effects of the injection temperature on Qpres and Qwall. As shown in Figure 21a,b, respectively, the heating effect of Qwall reduced, whereas the cooling effect of Qpres increased with the enhancement of the injection temperature. Therefore, with the increase in the injecting temperature, the CO2 temperature gradient decreased.
Figure 22 illustrates the pressure profiles of CO2 with various injection temperatures after 120 min. As stated above, the hydrostatic pressure increased with the reduction of the temperature. Therefore, the lower the temperature, the higher the pressure in the vertical section. However, given that the hydrostatic pressure gradient disappeared, the pressure profiles of the three cases in the horizontal section were slightly different.

4.3. Analysis of the Tube Inner Diameters

The cases where the tube inner diameters were 62, 76, and 100.3 mm were considered on the bases of the benchmark case. Table 2 shows the specific size of the tube.

4.3.1. Analysis of the Time Domain

Figure 23 shows the variation of the BHT with time under various tube inner diameters. The smaller the inner diameter of the tube, the faster the cooling rate. However, the CO2 temperature of the 62-mm case was 11.8 °C higher than that of the 100.3-mm case after 120 min. This phenomenon was caused by the variation of the convective heat transfer efficiency and Qfric. As shown in Figure 24, the flow velocity enhanced obviously with the decrease in the tube inner diameter, which resulted in a higher convective heat transfer effectiveness. However, the extra friction from the higher flow velocity Qfric increased by 3489 W/m, which made the lowest temperature higher in the 62-mm case. Figure 25 shows the variation of the WHP with time under various tube inner diameters. Its effects on WHP were similar to that of the injection rate. When the inner diameter decreased from 100.3 mm to 62 mm, the wellhead pressure increased by 79.6 MPa.

4.3.2. Analysis of the Space Domain

Figure 26 shows the temperature profiles of CO2 with various tube inner diameters after 120 min. In the vertical section, when the inner diameter increased, the effect of viscous dissipation decreased and the CO2 temperature decreased. However, in the horizontal section, the temperature gradient of the 62-mm case decreased significantly. Although the viscous dissipation increased with the increase in the wellbore friction, the cooling effect of the pressure work also became increasingly obvious. The combined effect of the two factors determined the temperature gradient of CO2. Figure 27a provides the sum of Qfric and Qpres under various tube inner diameters after 120 min, indicating that the smaller the inner diameter, the greater the heat flow gradient, especially in the horizontal section. The heat flow reduction rate was very fast in the 62-mm case: at a depth of 2300 m, the cooling effect of the pressure work exceeded the warming effect of the viscous dissipation, which created a negative sum of the two heat flows. Although the total heat flow remained positive throughout the wellbore due to the effect of the wall heat transfer (as shown in Figure 27b), the temperature gradient was anticipated to be negative along the wellbore if the horizontal section continued to extend in the 62-mm case.
Figure 28 indicates the pressure profiles of CO2 with various tube inner diameters after 120 min. In the vertical section, the hydrostatic pressure of the 100.3-mm case was greater than the wellbore friction, and the pressure gradient was positive; the wellbore friction of the 76- and 62-mm cases was greater than the hydrostatic pressure, and the pressure gradient was negative. In the horizontal section, as the hydrostatic pressure gradient was 0, the pressure gradients of the three cases were controlled by the wellbore friction, which made them all negative. The simulation results showed that the wellbore friction could be reduced significantly by increasing the inner diameter. Therefore, by increasing the inner diameter, the maximum of the CO2 injection rate could be increased, which is beneficial for improving the sand-carrying capacity of the CO2 fracturing fluid. At present, the newly developed 3000 type fracturing pump can stably provide 140 MPa pumping pressure [65]. Under the parameters of the benchmark case, using the wellhead pressure of 140 MPa as the ultimate pressure, the ultimate injection rate corresponding to different inner diameters was obtained as shown in Figure 29. Apparently, when the inner diameter was enlarged from 50.3 mm to 100.3 mm, the ultimate injection rate was increased from 2.7 m3/min to 29.6 m3/min.

4.4. Analysis of the Completion Method of the Horizontal Section

The cases where the completion methods of the horizontal section were casing completion and open-hole completion were considered on the bases of the benchmark case. The main difference between the two methods was that casing and cement sheath existed between the annulus and the formation in the casing completion, which were otherwise absent in the open-hole completion.

4.4.1. Analysis of the Time Domain

Figure 30 shows the variation of the BHT with time under various completion methods. The BHT of the open-hole completion in the horizontal section was 2.7 °C higher than that of the casing completion after 120 min. The variation of the annulus natural convection explained this phenomenon. In the open-hole completion case, the radius of the annulus was much larger than that of the casing case, which enhanced the efficiency of the natural convection. As shown in Figure 31a, compared with the casing completion case, the average annulus natural convection heat transfer coefficient of open-hole completion was 30 W/(m2·°C) higher. The enhancement of the natural convection in the annulus improved the wall heat transfer efficiency. Figure 31b indicates that the average Qwall with open-hole completion method was 313 W/m higher than that of the casing completion method. Figure 32 shows the variation of the WHP with time under various completion methods. Apparently, the completion method had little effect on the WHP. This finding is further discussed in the analysis of the space domain.

4.4.2. Analysis of the Space Domain

Figure 33 shows the temperature profiles of CO2 with various completion methods after 120 min. In the vertical and curved sections, the temperature profiles coincided, and the difference between the two cases only appeared after entering the horizontal section. As discussed above, the difference in the Qwall was the main cause of the temperature variation. Figure 34 provides the Qwall profiles with different horizontal section completion methods, which indicated that Qwall was the same in the vertical and curved sections, whereas a jump (ΔQwall = 817 W/m) occurred at the boundary between the curved and horizontal sections, which increased the CO2 temperature gradient in the horizontal section in the open-hole completion method. Therefore, open-hole completion for the horizontal section was beneficial to increase the BHT.
Figure 35 shows the variation of the CO2 pressure profiles with time under various completion methods after 120 min. Firstly, as the completion method did not change the flow characteristics, the wellbore friction was the same at different wellbore sections. Secondly, the variation of CO2 density occurred at the horizontal section, in which the hydrostatic pressure gradient was 0. Therefore, the completion method had little effect on CO2 pressure in both the time and space domains.
In addition, the usage of the open-hole completion method also significantly increased the cooling effect of CO2 on the formation. As shown in Figure 36a,b, after injecting 120 min, the formation temperature near the wellbore of the open-hole completion case was lower by 24.5 °C than in the casing completion case.

5. Conclusions

A 1D + 1D numerical model was established to investigate the horizontal wellbore temperature field in CO2 fracturing, which can provide an inexpensive and efficient way to investigate and guide the CO2 fracturing design. The convergence test and the comparison with the COMSOL software were performed to verify the model. The following conclusions were drawn on the bases of the assumptions and conditions of this work:
  • The higher the injection rate, the faster the reduction of temperature and the shorter the time for the BHT to reach stability. The relationship between the stable BHT and the injection rates is determined by the heat transfer mechanisms, which are non-monotonic. Under the parameters of the benchmark case, after 120 min, the corresponding injection rates were 2 m3/min, 10 m3/min, and 6 m3/min, according to the sequence of BHT from large to small.
  • In the horizontal section, the cooling effect of the pressure work and the heating effect of the wall heat transfer are enhanced. Under the combined effect of the two, the temperature gradient in the horizontal section was higher than that in the vertical section when the injection rate was 2 m3/min. The opposite was the case where the injection rate was 10 m3/min. Therefore, compared with the vertical section, the existence of the horizontal section will increase the BHT when the injection rate is low and reduce the BHT when the injection rate is high.
  • Wellbore friction increases exponentially with the increase in the injection rate. Enhancing the tube diameter can effectively reduce the wellbore friction, thus extra attention to casing fracturing is necessary. Based on the benchmark case and taking WHP at 140 MPa as the ultimate pressure, the ultimate injection rate increased from 2.7 m3/min to 29.6 m3/min when the tube diameter increased from 50.3 mm to 100.3 mm.
  • The open-hole completion method in the horizontal section enlarges the annulus space more than the casing completion method, which improves the natural convection heat transfer efficiency of the annulus. Based on the benchmark case, the BHT of the open-hole completion in the horizontal section was 2.7 °C higher than that of the casing completion after 120 min, whereas its effects on pressure were negligible. In addition, the open-hole completion in the horizontal section will significantly enhance the cooling effect of CO2 on the formation near the wellbore. After 120 min, the formation temperature of the bottom-hole with the open-hole method was 24.5 °C lower than that in the casing completion method.

Author Contributions

Conceptualization, S.Z.; Data curation, C.Z.; Funding acquisition, Y.H.; Investigation, X.L.; Methodology, X.L.; Software, X.L. and C.Z.; Supervision, Y.H. and Z.Z.; Writing—original draft, X.L.; Writing—review and editing, S.Z., C.Z. and Z.M. All authors have read and agreed to the published version of the manuscript.

Funding

The authors would like to acknowledge the National Natural Science Foundation of China (No. 51504266), National Science and Technology Major Project (No. 2016ZX05014-005-012), National Science and Technology Major Project (No. 2017ZX05049-006), “An Experimental Study on Prepad CO2 Energized Fracturing Stimulation Project” and “A Research on Prepad CO2 Energized Fracturing Technology of Horizontal Wells in Hahu Glutenite Reservoir Project” for their financial support.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

Across sectional area of tubing, m2
BHTbottom-hole temperature, °C
BHPbottom-hole pressure, Pa
Cpisobaric heat capacity, J/(kg·°C)
Cpcisobaric heat capacity of casing, J/(kg·°C)
Cpfisobaric heat capacity of CO2, J/(kg·°C)
Cptisobaric heat capacity of tubing, J/(kg·°C)
Cpanisobaric heat capacity of annulus fluid, J/(kg·°C)
Cphisobaric heat capacity of cement, J/(kg·°C)
Cprisobaric heat capacity of rock, J/(kg·°C)
dtitubing internal diameter, m
fDcoefficient of wellbore friction
ggravity acceleration, m/s2
gGgeothermal gradient, °C/m
GrGrasse number
hanheat transfer coefficient of natural convection, W/(m2·°C)
htheat transfer coefficient of forced convection, W/(m2·°C)
pCO2 pressure, Pa
pwbottom-hole pressure, Pa
Pfwellbore friction of CO2 fracturing, MPa/100 m
rradial distance, m
rtitubing internal radius, m
rtotubing external radius, m
rcicasing internal radius, m
rcocasing external radius, m
rhcement sheath radius, m
rr,iformation radius of unit i, m
ttime, s
Ttemperature, °C
vCO2 velocity, m/s
Tinjinjection temperature, °C
Tiniinitial temperature, °C
TfCO2 temperature, °C
Tttubing temperature, °C
Tanannulus temperature, °C
Tccasing temperature, °C
Thcement sheath temperature, °C
Tr,iformation temperature of unit i, °C
WHPwellhead pressure, Pa
zmeasured depth, m
zbhmeasured depth of the bottom-hole, m
αpcoefficient of CO2 thermal expansion, 1/°C, take −1/ρf(∂ρ/∂Tf) here
βancoefficient of annulus thermal expansion, 1/°C, take 2.5 × 10−4 here
θdeviation angle
ρdensity, kg/m3
ρfdensity of CO2, kg/m3
ρtdensity of tubing, kg/m3
ρandensity of annulus fluid, kg/m3
ρcdensity of casing, kg/m3
ρhdensity of cement, kg/m3
ρrdensity of rock, kg/m3
μfviscosity of CO2, Pa·s
μanviscosity of annulus, Pa·s
λthermal conductivity, W/(m·°C)
λfthermal conductivity of CO2, W/(m·°C)
εanthermal conductivity of annulus fluid, W/(m·°C)
λcthermal conductivity of casing, W/(m·°C)
λhthermal conductivity of cement, W/(m·°C)
λrthermal conductivity of rock, W/(m·°C)

References

  1. Jiang, J.; Yunis, R.M. Numerical study of complex fracture geometries for unconventional gas reservoirs using a discrete fracture-matrix model. J. Nat. Gas Sci. Eng. 2015, 26, 1174–1186. [Google Scholar] [CrossRef]
  2. Zhang, L.; Li, D.; Li, L.; Lu, D. Development of a new compositional model with multi-component sorption isotherm and slip flow in tight gas reservoirs. J. Nat. Gas Sci. Eng. 2014, 21, 1061–1072. [Google Scholar] [CrossRef]
  3. Yuan, J.; Luo, D.; Feng, L. A review of the technical and economic evaluation techniques for shale gas development. Appl. Energy 2015, 148, 49–65. [Google Scholar] [CrossRef]
  4. IEA. Prospects for CO2 Capture and Storage; IEA/OECD: Paris, France, 2004. [Google Scholar]
  5. Aminu, M.D.; Nabavi, S.A.; Rochelle, C.A.; Manovic, V. A review of developments in carbon dioxide storage. Appl. Energy 2017, 208, 1389–1419. [Google Scholar] [CrossRef] [Green Version]
  6. Liu, F.; Ellett, K.; Xiao, Y.; Rupp, J.A. Assessing the feasibility of CO2 storage in the New Albany Shale (Devo-nian–Mississippian) with potential enhanced gas recovery using reservoir simulation. Int. J. Greenh. Gas Control 2013, 17, 111–126. [Google Scholar] [CrossRef]
  7. Nuttall, B.C. Reassessment of CO2 sequestration capacity and enhanced gas recovery potential of Middle and Upper Devonian Black Shales in the Appalachian Basin. In MRCSP Phase II Topical Report, October 2005–October 2010; Kentucky Geological Survey: Lexington, KY, USA, 2010. [Google Scholar]
  8. Liu, J.; Xie, L.Z.; He, B.; Zhao, P.; Ding, H.-Y. Performance of free gases during the recovery enhancement of shale gas by CO2 injection: A case study on the depleted Wufeng–Longmaxi shale in northeastern Sichuan Basin, China. Pet. Sci. 2021, 18, 530–545. [Google Scholar] [CrossRef]
  9. Wang, T.Y.; Tian, S.C.; Liu, Q.L.; Li, G.S.; Sheng, M.; Ren, W.X.; Zhang, P.P. Pore structure characterization and its effect on methane adsorption in shale kerogen. Pet. Sci. 2021, 18, 565–578. [Google Scholar] [CrossRef]
  10. Lyu, X.; Zhang, S.; Ma, X.; Wang, F.; Mou, J. Numerical study of non-isothermal flow and wellbore heat transfer characteristics in CO2 fracturing. Energy 2018, 156, 555–568. [Google Scholar] [CrossRef]
  11. Lyu, X.; Zhang, S.; Ma, X.; Wang, F.; Mou, J. Numerical investigation of wellbore temperature and pressure fields in CO2 fracturing. Appl. Therm. Eng. 2018, 132, 760–768. [Google Scholar] [CrossRef]
  12. Schlumberger, M.; Perebinossoff, A.A.; Doll, H.G. Temperature measurements in oil wells. J. Pet. Technol. 1937, 23, 159. [Google Scholar]
  13. Ramey, H.J. Wellbore heat transmission. J. Pet. Technol. 1962, 14, 427–435. [Google Scholar] [CrossRef]
  14. Satter, A. Heat Losses During Flow of Steam Down a Wellbore. J. Pet. Technol. 1965, 17, 845–851. [Google Scholar] [CrossRef]
  15. Willhite, G.P. Over-all Heat Transfer Coefficients in Steam and Hot Water Injection Wells. J. Pet. Technol. 1966, 19, 607–615. [Google Scholar] [CrossRef]
  16. Romero-Juarez, A. A Simplified Method for Calculating Temperature Changes in Deep Wells. J. Pet. Technol. 1979, 31, 763–768. [Google Scholar] [CrossRef]
  17. Shiu, K.C.; Beggs, H.D. Predicting Temperatures in Flowing Oil Wells. J. Energy Resour. Technol. 1980, 102, 2–11. [Google Scholar] [CrossRef]
  18. Fontanilla, J.P.; Aziz, K. Prediction of Bottom-hole Conditions for Wet Steam Injection Wells. J. Can. Pet. Technol. 1982, 21, 82–88. [Google Scholar] [CrossRef]
  19. Alves, I.N.; Alhanati, F.J.S.; Shoham, O. A Unified Model for Predicting Flowing Temperature Distribution in Wellbores and Pipelines. SPE Prod. Eng. 1992, 7, 363–367. [Google Scholar] [CrossRef]
  20. Izgec, B.; Kabir, C.S.; Zhu, D.; Hasan, A.R. Transient Fluid and Heat Flow Modeling in Coupled Well-bore/Reservoir Systems. SPE Res. Eval. Eng. 2007, 10, 294–301. [Google Scholar] [CrossRef] [Green Version]
  21. McSpadden, A.R.; Coker, O.D. Multiwell Thermal Interaction: Predicting Wellbore and Formation Tempera-tures for Closely Spaced Wells. SPE Drill. Completion 2010, 25, 448–457. [Google Scholar] [CrossRef]
  22. Gu, H.; Cheng, L.; Huang, S.; Du, B.; Hu, C. Prediction of thermophysical properties of saturated steam and wellbore heat losses in concentric dual-tubing steam injection wells. Energy 2014, 75, 419–429. [Google Scholar] [CrossRef]
  23. Gu, H.; Cheng, L.; Huang, S.; Bo, B.; Zhou, Y.; Xu, Z. Thermophysical properties estimation and performance analysis of superheated-steam injection in horizontal wells considering phase change. Energy Convers. Manag. 2015, 99, 119–131. [Google Scholar] [CrossRef]
  24. Gu, H.; Cheng, L.; Huang, S.; Li, B.; Shen, F.; Fang, W.; Hu, C. Steam injection for heavy oil recovery: Modeling of wellbore heat efficiency and analysis of steam injection performance. Energy Convers. Manag. 2015, 97, 166–177. [Google Scholar] [CrossRef]
  25. Cheng, W.-L.; Han, B.-B.; Nian, Y.-L.; Wang, C.-L. Study on wellbore heat loss during hot water with multiple fluids injection in offshore well. Appl. Therm. Eng. 2016, 95, 247–263. [Google Scholar] [CrossRef]
  26. Nian, Y.-L.; Cheng, W.-L. Insights into heat transport for thermal oil recovery. J. Pet. Sci. Eng. 2017, 151, 507–521. [Google Scholar] [CrossRef]
  27. Cai, J.; Duan, Y. Study on temperature distribution along wellbore of fracturing horizontal wells in oil reservoir. Petroleum 2015, 1, 358–365. [Google Scholar] [CrossRef] [Green Version]
  28. Raymond, L.R. Temperature Distribution in a Circulating Drilling Fluid. J. Pet. Technol. 1969, 21, 333–341. [Google Scholar] [CrossRef]
  29. Eickmeier, J.R.; Ersoy, D.; Ramey, H.J. Wellbore Temperatures and Heat Losses During Production or Injection Operations. J. Can. Pet. Technol. 1970, 9, 115–121. [Google Scholar] [CrossRef]
  30. You, J.; Rahnema, H.; McMillan, M.D. Numerical modeling of unsteady-state wellbore heat transmission. J. Nat. Gas Sci. Eng. 2016, 34, 1062–1076. [Google Scholar] [CrossRef]
  31. Dong, W.; Shen, R.; Liang, Q. Model Calculations and Factors Affecting Wellbore Temperatures During SRV Fracturing. Arab. J. Sci. Eng. 2018, 43, 6475–6480. [Google Scholar] [CrossRef]
  32. Abdelhafiz, M.M.; Hegele, L.A.; Oppelt, J.F. Numerical transient and steady state analytical modeling of the wellbore temperature during drilling fluid circulation. J. Pet. Sci. Eng. 2020, 186, 106775. [Google Scholar] [CrossRef]
  33. Zhang, Z.; Xiong, Y.; Mao, L.; Lu, J.; Wang, M.; Peng, G. Transient temperature prediction models of wellbore and formation in well-kick condition during circulation stage. J. Pet. Sci. Eng. 2019, 175, 266–279. [Google Scholar] [CrossRef]
  34. Wu, Y.-S.; Pruess, K. An Analytical Solution for Wellbore Heat Transmission in Layered Formations (includes associated papers 23410 and 23411). SPE Reserv. Eng. 1990, 5, 531–538. [Google Scholar] [CrossRef] [Green Version]
  35. Hasan, A.R.; Kabir, C.S. Aspects of Wellbore Heat Transfer During Two-Phase Flow (includes associated papers 30226 and 30970). SPE Prod. Facil. 1994, 9, 211–216. [Google Scholar] [CrossRef]
  36. Hasan, A.R.; Kabir, C.S.; Lin, D. Analytic Wellbore Temperature Model for Transient Gas-Well Testing. SPE Reserv. Eval. Eng. 2005, 8, 240–247. [Google Scholar] [CrossRef]
  37. Hasan, A.R.; Kabir, C.S.; Wang, X. A Robust Steady-State Model for Flowing-Fluid Temperature in Complex Wells. SPE Prod. Oper. 2009, 24, 269–276. [Google Scholar] [CrossRef]
  38. Hasan, A.; Kabir, C. Modeling two-phase fluid and heat flows in geothermal wells. J. Pet. Sci. Eng. 2010, 71, 77–86. [Google Scholar] [CrossRef]
  39. Hasan, A.; Kabir, C. Wellbore heat-transfer modeling and applications. J. Pet. Sci. Eng. 2012, 86-87, 127–136. [Google Scholar] [CrossRef]
  40. Hagoort, J. Ramey’s Wellbore Heat Transmission Revisited. SPE J. 2004, 9, 465–474. [Google Scholar] [CrossRef]
  41. Lin, C.; Zhongren, Y. Construction and application of a wellbore temperature and pressure model for unsteady flow in gas wells. Nat. Gas Ind. 2017, 37, 70–76. [Google Scholar]
  42. Lu, M.; Connell, L. Non-isothermal flow of carbon dioxide in injection wells during geological storage. Int. J. Greenh. Gas Control 2008, 2, 248–258. [Google Scholar] [CrossRef]
  43. Haizhu, W.; Zhonghou, S.; Gensheng, L. Wellbore temperature and pressure coupling calculation of drilling with supercritical carbon dioxide. Pet. Explor. Dev. 2011, 38, 97–102. [Google Scholar] [CrossRef]
  44. Span, R.; Wagner, W. A New Equation of State for Carbon Dioxide Covering the Fluid Region from the Triple-point Temperature to 1100 K at Pressures up to 800 MPa. J. Phys. Chem. Ref. Data. 1996, 25, 151–159. [Google Scholar] [CrossRef] [Green Version]
  45. Fenghour, A.; Wakeham, W.A.; Vesovic, V. The Viscosity of Carbon Dioxide. J. Phys. Chem. Ref. Data 1998, 27, 31–44. [Google Scholar] [CrossRef]
  46. Vesovic, V.; Wakeham, W.A.; Olchowy, G.A.; Sengers, J.V.; Watson, J.T.R.; Millat, J. The Transport Properties of Carbon Dioxide. J. Phys. Chem. Ref. Data 1990, 19, 763–808. [Google Scholar] [CrossRef]
  47. Liangbin, D.; Gensheng, L.; Zhonghou, S.; Chunfang, W.; Gang, B. Wellbore Pressure and Temperature Prediction Model and Its Affecting Factors for CO2 Injection Wells. Pet. Explor. Dev. 2013, 41, 76–81. [Google Scholar]
  48. Li, X.; Li, G.; Wang, H.; Tian, S.; Song, X.; Lu, P.; Wang, M. A unified model for wellbore flow and heat transfer in pure CO2 injection for geological sequestration, EOR and fracturing operations. Int. J. Greenh. Gas Control 2017, 57, 102–115. [Google Scholar] [CrossRef]
  49. Sun, F.; Yao, Y.; Li, G.; Li, X. Geothermal energy development by circulating CO2 in a U-shaped closed loop geothermal system. Energy Convers. Manag. 2018, 174, 971–982. [Google Scholar] [CrossRef]
  50. Pan, L.; Oldenburg, C.M.; Wu, Y.-S.; Pruess, K. Wellbore flow model for carbon dioxide and brine. Energy Procedia 2009, 1, 71–78. [Google Scholar] [CrossRef] [Green Version]
  51. Pan, L.; Oldenburg, C.M. T2Well—An integrated wellbore–reservoir simulator. Comput. Geosci. 2014, 65, 46–55. [Google Scholar] [CrossRef]
  52. Pan, L.; Oldenburg, C.M.; Pruess, K.; Wu, Y.-S. Transient CO2 leakage and injection in wellbore-reservoir systems for geologic carbon sequestration. Greenh. Gases Sci. Technol. 2011, 1, 335–350. [Google Scholar] [CrossRef]
  53. Jianchun, G.; Ji, Z. A coupling model for wellbore transient temperature and pressure of fracturing with su-percitical carbon dioxide. Acta Pet. Sin. 2015, 36, 203–209. [Google Scholar]
  54. Jianchun, G.; Ji, Z.; Ran, Z.; Changlin, Z. A dual transient coupling model for wellbore of carbon dioxide injection well. Acta Pet. Sin. 2015, 36, 977–982. [Google Scholar]
  55. Gong, Q.; Xu, Z.; Wang, M.; Qin, J. Numerical investigation on wellbore temperature and pressure during carbon dioxide fracturing. Appl. Therm. Eng. 2019, 157, 113675. [Google Scholar] [CrossRef]
  56. Gong, Q.; Qin, J.; Lan, J.; Zhao, C.; Xu, Z. Numerical investigation of key parameter effects on temperature and pressure in wellbore during carbon dioxide fracturing. Heat Transf. Res. 2020, 51, 115–128. [Google Scholar] [CrossRef]
  57. Singhe, A.T.; Ursin, J.R.; Henninges, J.; Pusch, G.; Ganzer, L. Modeling of Temperature Effects in CO2 Injection Wells. Energy Procedia 2013, 37, 3927–3935. [Google Scholar] [CrossRef] [Green Version]
  58. Ruan, B.; Xu, R.; Wei, L.; Ouyang, X.; Luo, F.; Jiang, P. Flow and thermal modeling of CO2 in injection well during geological sequestration. Int. J. Greenh. Gas Control 2013, 19, 271–280. [Google Scholar] [CrossRef]
  59. Jiang, P.; Li, X.; Xu, R.; Wang, Y.; Chen, M.; Wang, H.; Ruan, B. Thermal modeling of CO2 in the injection well and reservoir at the Ordos CCS demonstration project, China. Int. J. Greenh. Gas Control 2014, 23, 135–146. [Google Scholar] [CrossRef]
  60. Cussler, E.L. Diffusion—Mass Transfer in Fluid Systems, 3rd ed.; Cambrige Press University: Cambridge, UK, 2009. [Google Scholar]
  61. Weishi, Z.; Zhenyun, S.; Weidong, S. To calculate and analyze friction along tubing string during liquid CO2 fracturing. Dring Prod. Technol. 2017, 40, 53–56. [Google Scholar]
  62. Campbell, S.M.; Fairchild, N.R., Jr.; Arnold, D.L. Liquid CO2 and Sand Stimulations in the Lewis Shale, San Juan Basin, New Mexico: A Case Study. SPE Rocky Mountain Regional/Low-Permeability Reservoirs Symposium and Exhibition. In Proceedings of the SPE Rocky Mountain Regional/Low-Permeability Reservoirs Symposium and Exhibition, Denver, CO, USA, 12–15 March 2020. [Google Scholar]
  63. Incopera, F.P.; DeWitt, D.P. Fundamentals of Heat and Mass Transfer, 4th ed.; John Wiley and Sons Inc: Hoboken, NJ, USA, 2011. [Google Scholar]
  64. COMSOL. Multiphysics User’s Guide: Pipe Flow Module User’s Guide; COMSOL: Burlington, MA, USA, 2018; pp. 75–77. [Google Scholar]
  65. Peng, J.; Zhou, Q.; Dai, Q.; Wang, L.; Liu, Y.; He, Y.; Huang, X. Development Status and Analysis of Domestic Large-scale Fracturing Equipment. China Pet. Mach. 2016, 44, 82–86. [Google Scholar]
Figure 1. Physical model.
Figure 1. Physical model.
Sustainability 13 05672 g001
Figure 2. Grids of the 2D model.
Figure 2. Grids of the 2D model.
Sustainability 13 05672 g002
Figure 3. Radial grids of the vertical (curved) section.
Figure 3. Radial grids of the vertical (curved) section.
Sustainability 13 05672 g003
Figure 4. Radial grids of the horizontal section.
Figure 4. Radial grids of the horizontal section.
Sustainability 13 05672 g004
Figure 5. Physical model.
Figure 5. Physical model.
Sustainability 13 05672 g005
Figure 6. Test of the convergence; left: (a) T f k T f k 1 2 with various iteration steps; right: (b) Number of iterations with various time steps.
Figure 6. Test of the convergence; left: (a) T f k T f k 1 2 with various iteration steps; right: (b) Number of iterations with various time steps.
Sustainability 13 05672 g006
Figure 7. Comparison with commercial software; left: (a) Comparison of BHT; right: (b) Comparison of temperature profiles.
Figure 7. Comparison with commercial software; left: (a) Comparison of BHT; right: (b) Comparison of temperature profiles.
Sustainability 13 05672 g007
Figure 8. Variation of BHT with time under various Qinj.
Figure 8. Variation of BHT with time under various Qinj.
Sustainability 13 05672 g008
Figure 9. Mean velocities.
Figure 9. Mean velocities.
Sustainability 13 05672 g009
Figure 10. Bottom-hole temperature change rates.
Figure 10. Bottom-hole temperature change rates.
Sustainability 13 05672 g010
Figure 11. Variation of the WHP with time under various Qinj.
Figure 11. Variation of the WHP with time under various Qinj.
Sustainability 13 05672 g011
Figure 12. Wellbore friction.
Figure 12. Wellbore friction.
Sustainability 13 05672 g012
Figure 13. CO2 density profiles.
Figure 13. CO2 density profiles.
Sustainability 13 05672 g013
Figure 14. CO2 temperature profiles with various Qinj.
Figure 14. CO2 temperature profiles with various Qinj.
Sustainability 13 05672 g014
Figure 15. Heat flow profiles with various Qinj; left: (a) Profiles of Qfric; middle: (b) Profiles of Qpres; right: (c) Profiles of Qwall.
Figure 15. Heat flow profiles with various Qinj; left: (a) Profiles of Qfric; middle: (b) Profiles of Qpres; right: (c) Profiles of Qwall.
Sustainability 13 05672 g015aSustainability 13 05672 g015b
Figure 16. CO2 pressure profiles with various Qinj.
Figure 16. CO2 pressure profiles with various Qinj.
Sustainability 13 05672 g016
Figure 17. Variation of BHT with time under various Tinj.
Figure 17. Variation of BHT with time under various Tinj.
Sustainability 13 05672 g017
Figure 18. Variation of WHP with time under various Tinj.
Figure 18. Variation of WHP with time under various Tinj.
Sustainability 13 05672 g018
Figure 19. CO2 temperature profiles with various Tinj.
Figure 19. CO2 temperature profiles with various Tinj.
Sustainability 13 05672 g019
Figure 20. Difference between BHT and wellhead temperature (hereinafter referred to as WHT) with various Tinj.
Figure 20. Difference between BHT and wellhead temperature (hereinafter referred to as WHT) with various Tinj.
Sustainability 13 05672 g020
Figure 21. Heat flow profiles with various Tinj; left: (a) Profiles of Qpres; right: (b) Profiles of Qwall.
Figure 21. Heat flow profiles with various Tinj; left: (a) Profiles of Qpres; right: (b) Profiles of Qwall.
Sustainability 13 05672 g021
Figure 22. CO2 pressure profiles with various Tinj.
Figure 22. CO2 pressure profiles with various Tinj.
Sustainability 13 05672 g022
Figure 23. Variation of BHT with time under various dti.
Figure 23. Variation of BHT with time under various dti.
Sustainability 13 05672 g023
Figure 24. Mean v and mean Qfric under various dti.
Figure 24. Mean v and mean Qfric under various dti.
Sustainability 13 05672 g024
Figure 25. Variation of WHP with time under various dti.
Figure 25. Variation of WHP with time under various dti.
Sustainability 13 05672 g025
Figure 26. CO2 temperature profiles with various dti.
Figure 26. CO2 temperature profiles with various dti.
Sustainability 13 05672 g026
Figure 27. Heat flow profiles with various dti; left: (a) Profiles of Qpres + Qfric; right: (b) Profiles of Qpres + Qfric + Qwall.
Figure 27. Heat flow profiles with various dti; left: (a) Profiles of Qpres + Qfric; right: (b) Profiles of Qpres + Qfric + Qwall.
Sustainability 13 05672 g027
Figure 28. CO2 pressure profiles with various dti.
Figure 28. CO2 pressure profiles with various dti.
Sustainability 13 05672 g028
Figure 29. Ultimate injection rates with various dti.
Figure 29. Ultimate injection rates with various dti.
Sustainability 13 05672 g029
Figure 30. Variation of BHT with time under various completion methods of the H section.
Figure 30. Variation of BHT with time under various completion methods of the H section.
Sustainability 13 05672 g030
Figure 31. Effects of the completion method of the horizontal section on mean han and mean Qwall; left: (a) Mean han; right: (b) Mean Qwall.
Figure 31. Effects of the completion method of the horizontal section on mean han and mean Qwall; left: (a) Mean han; right: (b) Mean Qwall.
Sustainability 13 05672 g031aSustainability 13 05672 g031b
Figure 32. Variation of WHP with time under various completion methods of the H section.
Figure 32. Variation of WHP with time under various completion methods of the H section.
Sustainability 13 05672 g032
Figure 33. CO2 temperature profiles with various completion methods.
Figure 33. CO2 temperature profiles with various completion methods.
Sustainability 13 05672 g033
Figure 34. Profiles with various completion methods.
Figure 34. Profiles with various completion methods.
Sustainability 13 05672 g034
Figure 35. CO2 pressure profiles with various completion methods.
Figure 35. CO2 pressure profiles with various completion methods.
Sustainability 13 05672 g035
Figure 36. Radial temperature distribution with various well completion methods after 120 min; left: (a) Casing completion; right: (b) Open-hole completion.
Figure 36. Radial temperature distribution with various well completion methods after 120 min; left: (a) Casing completion; right: (b) Open-hole completion.
Sustainability 13 05672 g036aSustainability 13 05672 g036b
Table 1. Parameters of the basic computation case.
Table 1. Parameters of the basic computation case.
ParametersValuesParametersValues
Depth1600 mSurface temperature20 °C
Injection temperature0 °CClosure pressure (hereinafter referred to as BHP)32 MPa
Injection rate6 m3/minGeothermal gradient0.03 °C/m
Tubing inner diameter76 mmInjecting time120 min
Tubing outer diameter89 mmTubing (Casing) heat capacity460 J/(kg·°C)
Casing inner diameter157.8 mmAnnulus capacity 4180 J/(kg·°C)
Casing outer diameter177.8 mmCement capacity 880 J/(kg·°C)
Cement sheath diameter237.8 mmFormation capacity833 J/(kg·°C)
Tubing (Casing) density7800 kg/m3Tubing (Casing) heat conductivity 53 W/(m·°C)
Annulus density1000 kg/m3Annulus heat conductivity0.557 W/(m·°C)
Cement density2000 kg/m3Cement heat conductivity0.627 W/(m·°C)
Formation density2505.53 kg/m3Formation heat conductivity2.5 W/(m·°C)
Curved section length200 mHorizontal section600 m
Completion method of the vertical (curved) section Casing completionCompletion method of the horizontal section Open-hole completion
Table 2. Size of the tube.
Table 2. Size of the tube.
Inner DiameterOuter DiameterOuter Diameter
Mmmmin
6273 2 7 8
7688.9 3 1 2
100.3114.3 4 1 2
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Lyu, X.; Zhang, S.; He, Y.; Zhuo, Z.; Zhang, C.; Meng, Z. Numerical Investigation on Wellbore Temperature Prediction during the CO2 Fracturing in Horizontal Wells. Sustainability 2021, 13, 5672. https://doi.org/10.3390/su13105672

AMA Style

Lyu X, Zhang S, He Y, Zhuo Z, Zhang C, Meng Z. Numerical Investigation on Wellbore Temperature Prediction during the CO2 Fracturing in Horizontal Wells. Sustainability. 2021; 13(10):5672. https://doi.org/10.3390/su13105672

Chicago/Turabian Style

Lyu, Xinrun, Shicheng Zhang, Yueying He, Zihan Zhuo, Chong Zhang, and Zhan Meng. 2021. "Numerical Investigation on Wellbore Temperature Prediction during the CO2 Fracturing in Horizontal Wells" Sustainability 13, no. 10: 5672. https://doi.org/10.3390/su13105672

APA Style

Lyu, X., Zhang, S., He, Y., Zhuo, Z., Zhang, C., & Meng, Z. (2021). Numerical Investigation on Wellbore Temperature Prediction during the CO2 Fracturing in Horizontal Wells. Sustainability, 13(10), 5672. https://doi.org/10.3390/su13105672

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop