Abstract
Unraveling the seepage and heat transfer coupled process between the working fluid and the fractured rocks in geothermal reservoir is of great significance to the exploitation and utilization of geothermal resources. In this study, based on the numerical modeling of fracture flow in geothermal reservoirs, the seepage and convective heat transfer behavior of fractured granite during geothermal extraction were investigated, and the effects of different fracture patterns, fluid injection temperature, and injection velocity on the temperature evolution of rock mass were comparatively analyzed. The results clearly revealed the influence of five different fracture patterns, i.e., single fracture, echelon fracture, parallel fracture, Y-shaped fracture, and crossfracture, on heat transfer capacity. Generally, the echelon fracture has the highest fluid outlet temperature, while the crossfracture has the largest total heat and shows the strongest heat transfer capacity. Besides the fracture shape, the fluid injection temperature and injection velocity also play significant roles in the heat transfer performance. The increase of fluid injection temperature would improve the total outlet heat of the crossfracture system and also benefit the system life. When the fluid injection temperature raises from 20°C to 35°C, the system life and the total outlet heat would be increased by 58.76% and 22.42%, respectively. However, higher fluid injection velocity would damage the system life of the geothermal reservoir, which obviously results in a decrease in total outlet heat. With the fluid injection velocity increasing from 0.004 m/s to 0.006 m/s, the system life could drop by 72.47%, and the total outlet heat was reduced by 55.72%. This work contributes to the preliminary understanding of the coupled seepage and heat transfer behavior in rock mass with various fracture patterns, and it could provide some practical implications for the rational exploitation of geothermal resources.
1. Introduction
The exploration and utilization of deep geo-energy and geo-resources, such as geothermal, have become a strategic developing tendency to meet the global energy demand [1]. However, as the reservoir depth advances, the commercial development of deep geo-resources is currently faced with a shortage of scientific understanding and applicable theory [2–4]. Generally, geothermal energy is a kind of clean and environmentally friendly renewable resources contained in the interior of the earth. Compared with other renewable resources, e.g., solar energy, wind energy, and ocean energy, geothermal has many advantages including wide distribution, huge reserves, and independence of climate and seasons. Most of geothermal energy exists in the shallow part of the crust in the form of steam, hot water, hot dry rock, and magma, and it has a promising development prospect since the current discovery and utilization rate of geothermal energy is quite low.
As a majority of deep geothermal resources, hot dry rock generally requires artificial reservoir stimulation technologies, e.g., enhanced geothermal system (EGS), to create reservoir fractures to enhance the heat exchange capacity. The exploitation of hot dry rock resources is characterized by large buried depth of reservoir, complex in situ environment, and high-temperature and high-stress coupling conditions [5]. Although the potential reserves of hot dry rock of China are huge, the development process of hot dry rock is relatively slow in China [6]. Technically speaking, the development of hot dry rock still faces many issues, among which the enhancement of seepage-heat transfer performance of different fracture patterns in reservoir rocks is the key topic that affects the efficiency of geothermal energy exploitation. With the increase of reservoir depth, the connection of natural fractures in hot dry rock becomes more complicated, and thus, it is challenging to capture the in situ evolution of fracture network [7–9]. At the same time, the seepage-heat transfer process of deep fractured rocks is extremely complex after the reservoir stimulation, so there always exists a large difference between the laboratory simulation tests and the practical engineering activities [10]. In this background, numerical modeling has proved a feasible method for practical prediction and interpretation of geothermal exploitation if in situ reservoir conditions could be elaborately considered.
Based on numerical simulations, Zhao et al. [11] proposed a THM- (thermo-hydro-mechanical-) coupled mathematical model to investigate the heat extraction process of hot dry rock, and a numerical analysis of geothermal extraction of high-temperature rock mass at the depth of 6000~7000 m was presented, which indicated that the temperature of rock mass, the matrix stress, and the stress around fractures decrease with the depth of geothermal extraction. Xin et al. [12] established a three-dimensional model of a single ideal rough fracture to reveal the distribution of fluid pressure, temperature, stress, and deformation, and a sensitivity analysis was further conducted to evaluate the effects of crack roughness and aperture on the outlet temperature. Sun et al. [13] simulated the EGS engineering example by using the two-dimensional randomly generated fracture model coupled with THM in fractured rock mass and studied the main parameters affecting the outlet temperature of EGS. Gao et al. [14] established a three-dimensional fracture model by scanning a real fractured granite and analyzed the impacts of fluid velocity, fluid temperature, rock temperature, and fracture pore size on heat transfer coefficient systematically, and a new representation of heat transfer coefficient was finally developed.
In summary, the previous studies have demonstrated that the main parameters such as fracture geometry, fluid velocity, fluid temperature, and rock temperature exert significant impacts on the fluid flow and heat transfer process of fractured rocks. To analyze the affecting degree of various parameters, many investigators have also carried out corresponding numerical simulation modeling and laboratory testing [15–17]. The related findings have an important propulsive effect on the development of geothermal exploitation engineering. However, most of the existing researches are focused on the effect of fracture roughness, fluid velocity, and temperature on single fracture flow, while the investigation on heat transfer efficiency in fractured rocks with various fracture patterns is limited.
In nature, the natural or artificial fractures of rock mass in geothermal reservoirs are complicated without a unique distribution [18, 19] (Figure 1). Figure 1 shows the complex fracture network of a granite outcrop, which consists of several main fracture patterns, e.g., single fracture, parallel fracture, crossfracture, echelon fracture, and Y-shaped fracture. The geometric characteristics of joints or fractures largely determine the mechanical properties, hydraulic characteristics, and quality of rock mass. Therefore, by considering five different fracture shapes, a thermo-hydro-coupled numerical model has been established in this work to comparatively reveal the effect of fracture patterns on seepage and heat transfer behaviors in rocks with various fracture patterns. The heat transfer performance of fluid through the crossfractures and the evolution of geothermal reservoir temperature have been further explored under the condition of different fluid velocities and fluid injection temperatures.
2. Modeling Methods
2.1. Model Assumptions
In this study, the common software COMSOL Multiphysics has been employed to unravel the heat transfer mechanism of different types of fractures during geothermal energy extraction. For the sake of simplicity, some model assumptions are made as follows: (1)Given that the reservoir rock mass is a discontinuous medium. The permeability of rock matrix is ignored, and the fluid only can flow through the fracture channel(2)The laminar flow of an incompressible Newtonian fluid is considered, and the fluid phase transition will not occur in the process of migration(3)The instantaneous local heat equilibrium is assumed to be reached in the process of convective conduction(4)The fracture width is much smaller than the fracture length. The fracture width is set as constant everywhere, and the fracture roughness is not considered
2.2. Governing Equations
The fluid flow within the fractures is assumed to follow the Navier-Stokes equation, and the leak off does not happen as the rock matrix is nearly impermeable. For the two-dimensional laminar flow, the governing equation can be expressed as [20–22]: where is the fluid velocity, is the fluid density, is the pressure, and is the dynamic viscosity of water.
For steady-state incompressible fluid, the mass conservation equation can be presented as:
Since this work is focused on the investigation of fluid flow behavior, the mechanical deformation of rock mass is not considered. The local thermal equilibrium theory, assuming that the fracture surface temperature is consistent with the fluid temperature at the local point, is adopted. Meanwhile, the heat transfer in rock matrix is mainly controlled by heat conduction, and thus, the energy conservation equation in steady state for rock matrix is as follows: where is the rock thermal conductivity and is the matrix temperature.
By contrast, the heat transfer in fluid is mainly controlled by thermal convection and conduction, so the corresponding energy conservation equation for fractures is expressed as [23]: where is the fluid density, is the specific heat capacity of the fluid under constant pressure, is the fluid temperature, is fluid heat capacity, and is the normal vector of the fracture surface.
2.3. Computational Domain
To analyze the effect of fracture shapes on fluid flow and heat transfer performance, five fractured granite models with different fracture patterns are established as shown in Figure 2. The numerical model has a length () of 50 mm, and a width () of 50 mm. The fracture width () is set as 0.2 mm. The fluid inflow temperature is denoted as , and the initial temperature of reservoir rock is denoted as . The cold fluid is injected into the fracture at a constant rate from the left side to the right side.
(a) Single fracture
(b) Echelon fracture
(c) Parallel fracture
(d) Y-shaped fracture
(e) Crossfracture
The fluid boundary conditions of the model are velocity inlet at the left side and pressure outlet at the right side, while impermeable boundary conditions are set at the top and bottom. The inlet velocity is fixed at 0.005 mm/s, while the outlet pressure remains zero at the right side. For the temperature boundary conditions, the initial temperature () of the whole system is prescribed as 100°C, and the temperature at the top and bottom of the model is fixed during the heat transfer process. The left and right sides of the model is assumed as thermally insulated, which means that the inlet and outlet boundaries are adiabatic.
Once the model geometry is created, the computational domain can be meshed for finite element computation. The grid resolution is an important factor affecting the precise of numerical simulation results. To speed up computation and ensure simulation accuracy, the physical field control grid method is adopted, and proper sizes of triangle elements are set up. The detailed grid of a single fracture model is illustrated in Figure 3.
3. Effect of Different Shapes of Fractures on Heat Extraction Performance
3.1. Model Configuration
Before computation, the finite element models are configurated with physical parameters. The fracture width is assumed as 0.2 mm everywhere, and the initial temperature of computation domain is prescribed as 100°C. To model geothermal energy extraction, cold water of 20°C is injected into the fracture at a velocity of 0.005 mm/s, and the seepage occurs from left to right in the fracture. Other key parameters of the model are summarized in Table 1. The source of the data can be found in Ref. [12].
3.2. Modeling of a Single Fracture Rock
The seepage and heat transfer coupled behaviors of a single fracture model were simulated first. Figure 4 shows the simulation results of temperature evolution in a single fracture model during the process of fluid injection. In the early stage () of fluid injection, a significant convective heat transfer occurs between the cold fluid and the high-temperature rock due to large temperature difference, and the heat transfer in rock mass maintains dynamic equilibrium before the fluid reaches the outlet. Thus, only the temperature field is locally disturbed near the inlet, where the temperature drops rapidly with significant heat loss. With the continuous injection of cold fluid, the rock temperature around the fracture surfaces decreases continuously, allowing the disturbance area to extend to the both sides of the fracture and the exit. In the late stage of injection, the temperature difference between the rock matrix and the fluid within the fracture decreases gradually, resulting in weaker heat transfer. The convection heat transfer in rock mass finally tends to be stable at .
(a)
(b)
(c)
To quantify the temperature evolution of injected fluid through the fracture channel, three gauge points were placed in the fracture channel along the direction of fracture length. The locations of gauges points are 10 mm, 20 mm, and 30 mm to the injector, respectively. Figure 5 shows the progressive variation of temperatures at the three observation points with the injection time. It can be seen that the temperature at goes down rapidly at the early stage of injection, while other points remain at a high level. With the increase of the distance from the observation point to the inlet, the observed temperature decreases more slowly and eventually reaches stable at a higher value.
To validate the modeling methodology, the simulation results of the seepage and heat transfer coupled behaviors of a single fracture model were compared with the results reported in previous literatures. Figure 6 shows that the simulation results of this work have a good agreement with the data taken from Ref. [24], which demonstrates the validity of the present modeling method.
3.3. Simulation Results and Analysis of Different Fractured Models
To reveal the effect of different fracture patterns on heat extraction, four types of nonsingle fracture models were established for comparison analysis. Figure 7 shows the temperature variation of different fractured rocks during the injection process. It can be seen that the overall tendency of temperature variation of complex fractured rocks is similar to that of the single fracture rock. The temperature field in the vicinity of the injector is disturbed at the initial stage. Afterwards, the disturbed area of rock mass will expand with the continuous injection, and the temperature decreases gradually along both sides of the fracture surface and the direction of seepage. This is due to the heat exchange between the fluid and the rock matrix as the cold fluid flows through the fracture, which causes the temperature of rock matrix near the fracture to decrease continuously.
(a) Echelon fracture
(b) Parallel fracture
(c) Y-shaped fracture
(d) Crossfracture
To evaluate the heat transfer capacity of various fracture models, the fluid outlet temperature and normal total heat flux in different cases were analyzed. The normal total heat flux of the outlet is computed as follows: where is the fluid outlet temperature, is the fluid inlet temperature, is the fluid density, is the fluid specific heat capacity, is the fluid velocity, and is the fracture width.
Figures 8 and 9 show the variations of the outlet temperature and normal total heat flux in fractured granite with five different patterns. Apparently, the fracture pattern has a considerable effect on the fluid outlet temperature. The echelon fracture model has the maximum fluid outlet temperature, followed by the crossfracture model, single fracture model, and Y-shaped fracture model, while the parallel fracture model has the lowest outlet temperature. This is because that the echelon fracture has the longest channel, which is also more curved than the other four. The fluid can absorb more heat through echelon fracture channel under the same boundary conditions, and thus, the average outlet fluid temperature of the echelon fracture model is the highest. Similarly, Figure 9 shows that the normal total heat fluxes at the outlet of fractured rocks with five different patterns are largely different. It can be clearly seen that the normal total heat flux of the outlet decreases rapidly with time in the initial stage and then tends to a stable state at .
Based on Figure 8, the total outlet normal heat during the whole heat transfer process can be further calculated (see Table 2). There is a difference in total outlet normal heat of fractured rocks with different patterns, among which the crossfracture can absorb most heat. Although the average outlet temperature of the crossfracture is not the highest among the five different patterns, the total outlet heat is the maximum because the fluid flow through the crossfracture channels is the largest under the same injection velocity.
4. Effects of Fluid Injection Temperature and Velocity on Heat Extraction Performance
The system life of EGS is important to assess the project benefit, and one of the key indexes of system life is the fluid outlet temperature of the production well. Generally, the operation time when the fluid outlet temperature is reduced to 70%~90% of the initial reservoir temperature is used as the service life of EGS [25, 26]. In this section, the service life is estimated from the time when the fluid outlet temperature drops to 85% of the initial rock temperature, i.e., the operation time when the outlet temperature drops to 85°C in this work.
Since the crossfracture model has the largest total heat flux at the fluid outlet, it is employed as an example to analyze the effects of fluid injection temperature and fluid injection velocity on the total outlet heat during its lifetime. The variation of average outlet temperature in the range of 15~35°C was considered, and different injection velocities in the range of 0.0040~0.0060 m/s were comparatively simulated. Figure 10 shows the variation of average outlet temperature at different fluid injection temperatures. It can be seen that high fluid injection temperature is conductive to the fluid outlet temperature. Besides, the decrease of fluid outlet temperature becomes slower with the increase of fluid injection temperature, and the corresponding system life is longer as it needs more time to reach a steady state. The longest system life is achieved when the fluid injection temperature is 35°C, which is increased by 58.76% compared with that at an injection temperature of 20°C. Figure 11 shows that the total outlet normal heat of various fractured models is positively correlated with the fluid injection temperature. The total normal heat at the outlet reaches the largest when the fluid injection temperature is 35°C, which is increased by 22.42% compared with that at an injection temperature of 20°C.
Figure 12 displays the variation of average fluid outlet temperature at different injection velocities. Apparently, the fluid injection velocity at the inlet has an adverse effect on the average outlet temperature. High fluid injection velocity will damage the average outlet temperature and system life. With the increase of the fluid injection velocity, the decline of outlet temperature is faster, resulting in a shorter lifetime. The system life is the shortest when the fluid injection velocity is 0.0060 m/s, which is reduced by 72.47% compared with that at an injection velocity of 0.004 m/s. Figure 13 shows the relationship between the total outlet normal heat and the fluid injection velocity during the system lifetime. It can be seen that the total outlet heat decreases with the increase of fluid injection velocity. When the fluid injection rate is 0.004 m/s, the total normal outlet heat is reduced by 55.72% compared with that at a fluid injection rate of 0.006 m/s.
5. Conclusions
To reveal the interaction between the rock matrix and fracture during the process of geothermal energy extraction, this work has numerically investigated the seepage and heat transfer coupled behaviors of fractured rocks with different fracture patterns, and the following conclusions can be drawn: (1)During the fluid injection, the temperature of fracture surface near the outlet decreases first, and the disturbance area extends gradually to inside of the rock matrix and the outlet. When the temperature evolution reaches a steady state in the later stage, the convective heat transfer would be weakened progressively(2)The fracture patterns play a significant role on the heat transfer capacity of fractured rocks. The average fluid outlet temperature would reach the maximum in the echelon fracture model, while the parallel fracture model observes the lowest fluid outlet temperature. Under the same conditions, the crossfracture has the largest total outlet heat and relatively high fluid outlet temperature, demonstrating a strong heat transfer capacity(3)The increase of fluid injection temperature helps to not only improve the system life but also increase the total outlet heat during the lifetime. When the fluid injection temperature rises up to 35°C, the system life can be enhanced by 58.76%, and the total outlet heat is improved by 22.42% compared with the results in the 20°C case(4)High fluid injection velocity appears to reduce the system life and the total outlet heat of fractured rocks in geothermal reservoirs during the system life period. Compared with the case of fluid injection velocity at 0.0040 m/s, the system life is reduced by 72.47%, while the total outlet heat is decreased by 55.72% when the fluid injection velocity increases to 0.006 m/s
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the Department of Science & Technology of Guangdong Province (2019ZT08G315), National Natural Science Foundation of China (U2013603), Shenzhen Science and Technology Innovation Commission (JCYJ20190808153416970), the Natural Resources Science and Technology Project of Hunan Province (no. 2020-12), and the State Key Laboratory of Hydraulics and Mountain River Engineering (SKHL1921).