Abstract
The study of fluid-heat coupling in deep fractured surrounding rock is the basis of design, safety, and extraction of geothermal energy of deep underground spaces. The heat transfer and fractured media seepage theories were employed to establish a three-dimensional unsteady model for fluid-heat coupling heat transfer in fractured surrounding rock. Using COMSOL multiphysics simulation software, the temperature field of the fractured surrounding rock was determined. Furthermore, the influences of ventilation time, Darcy’s velocity, fracture aperture, and thermal conductivity coefficient of the surrounding rock on the fractured surrounding rock temperature field distribution were investigated. The results of the numerical simulation show that the ventilation time and fracture have a major impact on the temperature field distribution of the fractured surrounding rock. As ventilation time is 200 days, an average water temperature in centerline of the fracture decreases 9.4 K as Darcy’s velocity increased from 3e-4m/s to 2e-3m/s. As ventilation time is 200 days, an average water temperature in centerline of the fracture decreases 5.3 K as fracture aperture increased from 3 mm to 9 mm. A set of experimental devices for fluid-heat coupling heat transfer in surrounding rock with a single fracture was designed and built to validate the numerical simulation results. Numerical simulation results are, in general, in agreement with the experimental results.
1. Introduction
Shallow mineral resources and spaces are gradually being reduced and exhausted. Deep mining and space use have become a trend in resource and space development [1]. The temperature of original rock and deep underground space increases with increasing depth of the underground space [2]. The high temperature restricts the exploitation efficiency of the deep resources, effective use of deep underground space, and comfort of personnel. High temperature can also induce rock mass collapse and gas explosion accidents in deep engineering, which seriously threaten the safety of deep underground spaces [3].
High temperature in deep underground spaces occurs mainly because of heat transfer from deep surrounding rock [3]. Geothermal energy from deep rock masses is a renewable energy. Extraction and use of geothermal energy from deep underground spaces have great potential for development and, compared with ground-source heat pump systems, have reduced drilling costs. Ghoreishi-Madiseh et al. [4] extracted geothermal energy from backfilled mines; heat from the underground surrounding rock was transformed into a sustainable geothermal heat source to provide low-cost and clean geothermal energy for the mining area, thereby improving the sustainability of the mining industry. The study of heat transfer in deep surrounding rock is very important for the exploitation of deep resources, the use of deep underground spaces, and the exploitation of mining geothermal resources [5].
Wang et al. [6] simplified the heat conduction of surrounding rock to one-dimensional unsteady heat conduction. The analytical solution for dimensionless temperature was obtained by using the method of separation of variables. Zhang et al. [7] simplified the temperature field of the surrounding rock to a one-dimensional unsteady heat conduction problem. These studies assumed no fluid flow in the surrounding rock and took into account only the conductive heat transfer of the surrounding rock.
In fact, most deep surrounding rocks are fractured. In fractured surrounding rocks, heat transfer and seepage flow are strongly coupled processes. According to seepage models for rock mass, the fluid-heat coupling heat transfer in surrounding rock can be organized in two categories: an equivalent continuum model and a discrete fracture network model [8].
The equivalent continuum model is suitable for high fracture development; a representative elementary volume (REV) of fractured rock mass exists and is not too large for the study area. In this model, the fractured rock mass is regarded as a porous medium for research, without considering the physical structure of a single fracture. The fluid-heat coupling heat transfer in the surrounding rock is established by using porous medium seepage and heat transfer [9–11].
Discrete fracture models are suitable for a rock mass with fracture network and sparse fractures where there are no representative elementary volumes. The discrete fracture model considers that the rock mass is impervious and that the entire movement of water is accomplished through the fracture network. A single fracture is the basic unit of the fracture network model. Considerable amount of research has been undertaken on the fluid-heat coupling heat transfer in the surrounding rock with a single fracture [12–14].
An analytical model for heat transfer in fractured surrounded rock was developed as a discrete fracture network model. The results demonstrated that longitudinal thermal diffusivity is a critical parameter that determines temperature distribution in the fracture [15]. A discrete fracture network geothermal reservoir model based on a parallel plate model was developed, and the model was validated for synthetic fracture systems using a discrete fracture network model with four vertical fractures and three intersections. In this model, each fracture is modeled explicitly as a parallel plate [16]. A three-dimensional coupled hydrothermal model for fractured rock based on the finite discrete element method was proposed. Analytical solutions were provided to verify the model. The effects of fracture aperture, fluid viscosity, and pressure difference between the fluid and rock were studied [17]. The fluid-heat coupling heat transfer based on a dual porosity model considering heat transfer and mass transfer was solved numerically, and the temperature field distribution of the fractured rock mass was determined [18]. A coupled thermal-hydraulic-mechanical model with discrete fractures was developed to simulate fracture fluid flow, heat transfer, and shearing dilation behaviors in a hot volcanic reservoir system. Fluid flows in the fracture were calculated based on the cubic law. Heat transfer modes within the investigated fracture were thermal conduction, thermal advection, and thermal dispersion [19].
Although scholars have performed much work on the fluid-thermal coupling of fractured rock mass with respect to nuclear waste burial, hot dry rock, oil and gas resource exploitation, and underground sequestration of carbon dioxide, they have paid less attention to deep fractured surrounding rock. Research on the fluid-heat coupling heat transfer in surrounding rock is therefore urgent.
In this study, to study the fluid-thermal coupling in surrounding rock with a single fracture, the discrete fracture network model was applied. Regarding deep fractured surrounding rock as a fractured and matrix rock block and considering the interaction between heat transfer and seepage flow, a transient heat transfer mathematical model for a single fracture of surrounding rock in deep underground space was proposed to calculate the temperature distribution at any time. COMSOL Multiphysics is direct coupling multiphysics analysis finite element simulation software. The software can be used on its own or expanded with functionality from any combination of add-on modules for simulating electromagnetics, structural mechanics, acoustics, fluid flow, heat transfer, and chemical engineering. COMSOL Multiphysics was used to numerically solve the fluid-heat coupling heat transfer in deep surrounding rock and obtain the heat transfer law of deep fractured surrounding rock. The influence of ventilation time, seepage velocity, fracture apertures, and thermal conductivity of the surrounding rock on the temperature field distribution of the surrounding rock was studied. A corresponding experimental simulation platform for fluid-thermal coupling of deep fractured surrounding rock with a single fracture was built for verification.
2. Heat Transfer Model for Deep Fractured Surrounding Rock
2.1. Basic Assumptions
The following assumptions were made for the heat transfer of deep fractured surrounding rock:
(1) Fractured surrounding rock with large fractures is a deformable rock mass composed of a matrix block and fractured rock mass with negligible water-holding capacity and permeability. The rock is homogeneous and isotropic. There is only fracture in the rock, and the aperture of the fracture is much smaller than the length of the fracture.
(2) The permeability of the fractured surrounding rock mass can be ignored and the mass is regarded as a pure solid. Groundwater can flow only in the fracture, the flow direction is constant, phase change in the fluid can be ignored, flow velocity is not affected by density and viscosity, and the seepage obeys Darcy’s law linearly [20].
(3) The heat in the rock matrix is transmitted by conduction and convection, and the influence of thermal radiation can be ignored.
(4) The temperature of the fluid between the rock matrix and adjacent points of the fracture is not equal to that of the rock matrix. Heat exchange between the fluid and surrounding rock matrix is through convection.
2.2. Governing Equation of Heat Transfer in Deep Fractured Surrounding Rock
Based on the above assumptions, the fluid-heat coupling model for surrounding rock included three parts: the temperature field governing equation of the surrounding rock matrix, the seepage field governing equation of the fracture water, and the temperature field governing equation of the fracture water.
2.2.1. Governing Equation of Temperature Field in Surrounding Rock Matrix
According to the heat conservation of the surrounding rock matrix, the heat absorbed by a microelement of surrounding rock matrix in a unit of time should be equal to the heat conducted by the microelement of the surrounding rock matrix itself. The governing equation of the temperature field of the surrounding rock matrix is given bywhere is the density of surrounding rock matrix, in kg/m3; is the specific heat of the surrounding rock matrix, in kJ/kg∙K; is the temperature of the surrounding rock matrix, in K; is the heat transfer time, in s; is the thermal conductivity of the surrounding rock matrix, in W/(m·K); is the convective heat transfer coefficient of the fracture water to the surrounding rock wall, in W/(m2·K); is the temperature of the fracture water, in K; and is the aperture degree of the fracture, in m.
2.2.2. Control Equation of Fracture Water Seepage Field
The fracture seepage satisfies the following continuity equation:where is the density of the fracture water, in kg/m3; and , , and are the components of the seepage velocity vectors in the , , and directions, respectively, in m/s.
The flow in the fracture follows Darcy’s laws, which are The seepage coefficient can be deduced from the cubic law [21], which is where is the head loss, in m, is the gravitational acceleration, in m/s2, and is the kinematic viscosity of water, in m2/s.
2.2.3. Governing Equation of Temperature Field of Fracture Water
According to the heat conservation of fracture water, the sum of the net energy carried by microelement fracture water through the interface in a unit of time and the net heat transferred by heat conduction at the interface is equal to the change rate of total energy of the microelement fracture water with time. This can be expressed aswhere is the thermal conductivity of the fracture water, in W/(m∙K); and is the specific heat of the fracture water at constant pressure, in kJ/kg∙K.
The initial and boundary conditions of the fracture water seepage field were set as follows: the initial velocity (pressure) and velocity (pressure) boundary of the fracture water and the peripheral contact surfaces of the fracture and surrounding rock matrix are set with a nonslip condition. The initial and boundary conditions for the temperature field of the fracture water were set as follows: the initial temperature of the fracture water and the inlet temperature. The coupled boundary condition is that the heat fluxes at the junction of fracture water and surrounding rock matrix are equal.
Coupling Eqs. (1)–(7) and supplementing with initial and boundary conditions, the fluid-heat coupling equation for deep fractured surrounding rock is obtained.
3. Experimental Verification
To verify the results of the numerical model, an experimental system of fluid-heat coupling heat transfer in a surrounding rock with a single fracture was designed and built. The experimental system consisted of a rock model, tunnel model, fracture model, air flow simulation system, thermal boundary simulation system, seepage condition simulation system, and measurement and monitoring system.
In the air flow simulation system, an air conditioning unit for controlled temperature and humidity was used to reproduce the air flow in deep underground space with different temperature, humidity, and wind speeds. The thermal boundary simulation system was a multiheating belt and temperature controller on the wall of the surrounding rock and was used to simulate geothermal heat. The seepage condition simulation system included a constant-temperature water tank, constant flow pump, solenoid valve, and flowmeter. The monitoring system controlled the air flow temperature, humidity, velocity, internal temperature of the surrounding rock, and temperature and velocity at the fracture entrance and exit and collected the data.
Experimental schematic and physical diagrams are shown in Figures 1 and 2. The dimensions of the experimental surrounding rock were 2.4 m long, 1 m wide, and 1 m high. The space section was set as a circle and located in the center of the surrounding rock. The space radius was 0.1 m, length and width of the fracture were 0.5 m each, and aperture was 8 mm. Thermal properties of the experimental surrounding rock matrix are listed in Table 1. As shown in Table 1, the thermal properties of the surrounding rock in the experimental model are smaller than the properties of the artificial rock. This is because the equivalent diameter in the experimental model is small compared to the equivalent diameter in the prototype, so, with being equal, smaller thermal conductivity is used in the experimental model. In addition, it is ensured that the bottom of the surrounding rock model is a constant-temperature thermal boundary for a sufficient period of heat transfer time. The measurement parameters used for the experimental setup and the type and accuracy of the sensor are shown in Table 2. The sensor position is shown in Figure 1.
(a) Schematic diagram
(b) Cross-sectional view of a central portion of the fracture with layout of measuring points
The numerical simulation, which was based on the experimental model, provided results for fluid-heat coupling heat transfer in the surrounding rock with a single fracture. For this analysis, the original rock temperature was 313 K, the initial and entrance temperatures of the air flow in the space were 293 K, and the entrance temperature of the fracture water was 293 K. The wind velocity was 2 m/s. The layout of the measuring points in the surrounding rock model is shown in Figure 1(b). The simulated temperatures of the nodes at different ventilation times (0–60 min) were calculated using COMSOL and compared with the experimental values, as shown in Figure 3.
As shown in Figure 3, as time increases, temperature of all measuring points decreases. Within the ventilation time range, measuring point 4 has the lowest temperature. This is because measuring point 4 is closest to the fracture where the temperature is lowest. Figure 3 indicates that the simulated values and experimental values are in good agreement. The longer the ventilation time is, the better the consistency is between the simulated and experimental temperatures.
4. Numerical Simulation
4.1. Calculation Model
In the calculation model, the deep fractured surrounding rock is composed of a single fracture. To accurately study the influence of the fracture on the temperature field of the fractured surrounding rock, rock with a single fracture in deep underground space is taken as the research object. The computational domain is a single-fracture surrounding rock with dimensions of 40 m × 20 m × 200 m. The shape of the deep underground space is a semicircular arch. The height of the space is 3.8 m, the width is 4.8 m, the height of the vertical wall is 1.4 m, the arch height is 2.4 m, and the bottom of the space is 4 m from the bottom of the surrounding rock. The fracture length is 20 m, the width is 50 m, the aperture of the fracture with horizontal arrangement is d, the angle between the fracture and the horizontal plane is 0°, the bottom of the fracture is 18 m from the surrounding rock bottom, and the seepage is two-dimensional unsteady flow. The seepage flows in at y = 0 m and out at y = 20 m. The physical model is shown in Figure 4.
4.2. Calculation Parameters
4.2.1. Physical Parameters
Thermal properties of the simulated surrounding rock matrix with a single fracture are listed in Table 3.
4.2.2. Initial and Boundary Conditions
The initial and boundary conditions of the surrounding rock matrix temperature field were set as follows: the initial temperature of the surrounding rock is 313 K; the bottom of the surrounding rock model is a constant-temperature thermal boundary with a temperature of 313 K; the other edges are thermally insulated; the initial temperature of the air flow is 293 K; the velocity of the air flow is 2 m/s; and the surface heat transfer coefficient () is given by the following equation [22]:where Nu is Nusselt number; Re is Reynolds number; Pr is Prandtl number; is the surface heat transfer coefficient, in W/(m2∙K); D is the diameter of space, in m.
The seepage in the calculation model neglects the flow in the z direction, and the flow is two-dimensional unsteady flow. This simplifies the continuity equation of the seepage field in fracture water to a two-dimensional equation. The initial and boundary conditions of the fracture water seepage field were set as follows: the initial velocity of the fracture water is 0 m/s; and the peripheral contact surfaces of the fracture and surrounding rock matrix are set with a nonslip condition. The entrance of the fracture seepage is set as the velocity boundary, and the velocity is given.
The initial and boundary conditions for the temperature field of the fracture water were set as follows: the initial temperature of the fracture water and the inlet temperature are 293 K.
The coupled boundary condition is that the heat fluxes at the junction of fracture water and surrounding rock matrix are equal.
4.3. Mesh Generation
The calculation model was meshed with extremely fine grids using COMSOL. A total of 311329 domain elements, 29200 boundary elements, and 1662 edge elements were obtained.
4.4. Results and Analysis
To study the influence of main factors on fluid-heat coupling heat transfer in surrounding rock with a single fracture, the effects of fracture, ventilation time, seepage velocity, fracture aperture, and thermal conductivity of the surrounding rock matrix on the temperature field distribution in the surrounding rock were studied by COMSOL.
4.4.1. Effect of Fracture on Surrounding Rock Temperature Field
Figure 5 shows the temperature distribution of the surrounding rock in different sections with a ventilation time of 300 days, seepage velocity of 6e-4 m/s, fracture aperture of 5 mm, and thermal conductivity of surrounding rock matrix of 2.035 W/(m·K).
Figure 5 indicates that the existence of the fracture changes the distribution of the temperature field in the surrounding rock. The closer the section is to the fracture, the higher the temperature of fracture water is and the lower the temperature of surrounding rock matrix is, owing to the heat exchange between fracture water and surrounding rock. The isothermal surface presents an asymmetric distribution related to the fracture. The fracture water brings heat energy from the surrounding rock to the downstream area of seepage, which causes the isotherms of the surrounding rock shown in Figure 5 to move to the right. The heat transfer of the fracture water does not affect the surrounding rock far from the fracture. The heat transfer mode for the surrounding rock is only heat conduction. The isothermal surface of the surrounding rock cross-section is similar in shape to that of the space cross-section, showing a regular and symmetrical distribution.
4.4.2. Effect of Ventilation Time
Figure 6 is a temperature field map of a deep single-fracture surrounding rock at x = 100 m with a seepage velocity of 6e-4 m/s, fracture aperture of 5 mm, and thermal conductivity of surrounding rock matrix of 2.035 W/(m·K); ventilation times of 50, 150, 200, and 250 days are shown, and the abscissa is the displacement in the y direction.
(a) 50 days
(b) 150 days
(c) 200 days
(d) 250 days
At a ventilation time of 50 days, separate temperature fields are formed between the fracture and the surrounding rock that do not affect each other. At 150 days of ventilation time, the two temperature fields have fused and interfere with each other. A yellow isotherm at a temperature of 306 K can be seen clearly in the figure. The area surrounded by the yellow isotherm is defined as the cooling zone. With increasing time, the cooling zone formed by seepage and ventilation of fracture water gradually expands and the cooling zone gradually offsets in the direction of the fracture water flow; when the ventilation time is 200 days, the distance between the right boundary of the cooling zone and fracture intersection and the left boundary of the model is approximately 22 m. When the ventilation time is 250 days, the distance between the right boundary of the cooling zone and fracture intersection and the left boundary of the model is approximately 27 m. The cooling zone gradually expands and offsets to the right with increasing ventilation time. Therefore, the longer the ventilation time is, the greater the disturbance of fracture water is to the temperature field of the surrounding rock.
4.4.3. Effect of Seepage Velocity
Figures 7 and 8 are the temperature distributions at cross-sectional line A (x = 100 m, z = 18.0025 m) and at cross-sectional line B (x = 100 m, z = 16 m), respectively. Figure 4 shows illustration cross-sectional lines A and B. In terms of conditions, the ventilation time is 200 days, fracture aperture is 5 mm, and thermal conductivity of the surrounding rock matrix is 2.035 W/(m·K); seepage velocities are 3e-4, 6e-4, 9e-4, and 2e-3 m/s, and the abscissa is the displacement in the y direction. In the figures, the figures water is at x = 100 m, z = 18.0025 m and surrounding rock matrix is at x = 100 m, z = 16 m. Figure 7 is the temperature distribution map of the fracture water, and Figure 8 is the temperature map of the surrounding rock matrix.
As shown in Figure 7, when the inlet temperature of fracture water is the same, the fracture water temperature decreases with increasing seepage velocity of fracture water. The main reason is that, with the increase in seepage velocity, the amount of fracture water and the amount of cooling provided by the fracture water increase, and the temperature of fracture water therefore increases slowly. Figure 8 indicates that the temperature distribution for the surrounding rock matrix is asymmetric owing to the influence of the fracture. With the increase in seepage velocity of the fracture water, the temperature of surrounding rock matrix on the line decreases. The main reason is that, with the increase in seepage velocity, the cooling capacity of the fracture water increases, the convective heat transfer effect between fracture water and surrounding rock matrix becomes more intense, and the temperature of surrounding rock matrix decreases faster.
4.4.4. Effect of Fracture Aperture
Figures 9 and 10 are the temperature distributions at cross-sectional lines A and B, respectively. In terms of conditions, the ventilation time is 200 days, the seepage velocity is 6e-4 m/s, and the thermal conductivity of the surrounding rock matrix is 2.035 W/(m·K); fracture apertures are 3, 5, 7, and 9 mm, and the abscissa is the displacement in the y direction.
As shown in Figure 9, the temperature of fracture water decreases with increasing fracture aperture. The aperture increase causes an increase in fracture water flow; more low-temperature fluid enters the fracture to exchange heat with the surrounding rock, which reduces the temperature of the seepage water in the fracture.
As shown in Figure 10, the temperature of surrounding rock decreases with increasing fracture aperture. The increase in fracture width causes increasingly more water to enter the fracture, and increasingly more heat is exchanged with the surrounding rock matrix, which steadily reduces the temperature of the surrounding rock matrix.
4.4.5. Effect of Thermal Conductivity of Surrounding Rock Matrix
Figures 11 and 12 are the temperature distributions at cross-sectional lines A and B, respectively. In terms of conditions, the ventilation time is 100 days, the seepage velocity is 6e-4 m/s, the fracture aperture is 5 mm, and the thermal conductivity of the surrounding rock matrix is 2.035 W/(m·K); thermal conductivities are 1.409, 2.035, 2.608, and 3.288 W/(m·K), and the abscissa is the displacement in the y direction.
Figures 11 and 12 indicate that the thermal conductivity of the surrounding rock matrix has little effect on the temperature of the fracture water or the surrounding rock matrix. With the increase in thermal conductivity, the temperature of fracture water increases and the temperature of the surrounding rock decreases. The reason is that, as the thermal conductivity increases, the heat transferred from the surrounding rock to the fracture water increases; the fracture water receives more heat and the surrounding rock loses much heat.
5. Conclusions
(1) Considering the surrounding deep fractured rock as composed of fracture and rock matrix blocks and considering the mutual coupling of heat transfer and seepage flow, a three-dimensional unsteady mathematical model for fluid-heat coupling heat transfer in the surrounding rock was established.
(2) Corresponding experimental simulation platforms for fluid-thermal coupling of deep fractured surrounding rock with a single fracture were built, and the numerical simulation results were verified by the experiments. The numerical simulation results and the experimental results are consistent.
(3) A geometric model for heat transfer of a deep surrounding rock fracture with dimensions of 40 m × 20 m × 200 m was established. The initial and boundary conditions of the model were set up. The coupling calculation was performed using COMSOL. The effects of ventilation time, seepage velocity, fracture aperture, and thermal conductivity of the surrounding rock on the temperature field distribution of the fracture and surrounding rock were studied.
(4) The numerical results show that the existence of the fracture reduces the temperature of the adjacent surrounding rock matrix, and the isotherm of the surrounding rock temperature field moves in the direction of the fracture water flow. At a ventilation time of 50 days, the fracture and the surrounding rock do not affect each other. At 150 days of ventilation time, the two temperature fields have fused and interfere with each other. Results show that the longer the ventilation time, the larger the cooling zone formed by seepage of fracture water and ventilation of the space.
When the seepage velocity of fracture water increases from 3e-4 to 2e-3 m/s, the water temperature at the fracture centerline decreases by 9.4 K; the average decrease is 3.13%. When the fracture aperture increases from 3 to 9 mm, the water temperature at the fracture centerline decreases by 5.3 K; the average decrease is 1.75%. When the thermal conductivity of surrounding rock matrix increases from 1.409 to 3.288 W/(m·K), the water temperature at the fracture centerline increases by 1.9 K; the average decrease is 0.62%.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study was supported by the National Natural Science Foundations of China (No. 51404191, 51774229).