Coupled Hydraulic-Thermal Modelling and Related Numerical Analysis on Rock Fractures

Zhejiang Engineering Survey and Design Institute Group Co. Ltd., Ningbo, 315010 Zhejiang Province, China Work Safety Key Lab on Prevention and Control of Gas and Roof Disasters for Southern Coal Mines, Hunan University of Science and Technology, Xiangtan, 411201 Hunan Province, China School of Resource Environment and Safety Engineering, Hunan University of Science and Technology, Xiangtan, 411201 Hunan Province, China


Introduction
Hydraulic-thermal coupling is one of the most important subsystems of the thermal-hydraulic-mechanical (THM) system; it is a new research area with wide application. The thermal-hydraulic coupling process and mechanism of rock mass are conducive to petroleum excavation, natural gas excavation, underground water excavation, geothermal resource development, and disposal of nuclear wastes in deep rock mass. Besides, there exist problems in large water conservation and hydropower project in the river and gorge; the related study is critically important for the stability of water conservation and hydropower project. At present, the study on the hydraulic-thermal coupling mechanism of rock mass with fissures is mainly based on the continuum and equivalent continuum seepage theory. With regard to the study of THM coupling, how to simulate the process of hydraulic-thermal coupling of rock mass with fissures with noncontinuum theory has become a hot topic.
McDermott et al. [1] presented such an approach for the simulation of fluid flow through a fracture validated against experimental data and cross-comparison with results of other modelling teams within the DECOVALEX 2015 project by replacing the mechanical behavior and chemical transport process within physical models, and some results of thermal-hydraulic-mechanical-chemical conditions for fracture rock were attained. Cacace and Jacquey [2] developed theory and numerical implementation describing groundwater flow and the transport of heat and solute mass in fully saturated rocks with elastoplastic mechanical feedbacks; in the model, fractures were considered being of lower dimension than the hosting deformable porous rock. Ogata et al. [3] developed a multiphysics numerical model to predict the fluid flow and mass transport behavior of rock fractures under coupled thermal-hydraulic-mechanicalchemical conditions. Particularly, the model was employed for the purpose of describing the evolution of permeability and reactive transport behavior within rock fractures by taking into account the geochemical processes of the freeface dissolution and the pressure dissolution. Liu et al. [4] proposed a coupled thermal-hydraulic-chemical model to study the influence of rock heterogeneity and the coupling effect of temperature, groundwater, and hydrochemistry on rock damage, and an actual deep coal mine model was established to predict water inrush. Furthermore, Shao et al. [5,6], Olsson and Barton [7], Xie et al. [8], Zhang and Nemcik [9], Tsang and Witherspoon [10], and many other scholars [11][12][13][14][15][16][17][18][19][20][21][22][23][24] studied the rheological fracture behavior of rock cracks subjected and the coupling effect of rock, which promoted the development of mechanical characteristics of rocks under complicated conditions.
Though some numerical simulation software, such as FLAC and 3DEC, can simulate the THM coupling, this coupling is one-way coupling, which only considered the influence of the thermal field to the seepage field. Meanwhile, the heat transform coefficient, heat conductivity coefficient, and width of cracks are constant; however, many researches proved that the heat transform coefficient h is related to the width of fissures.
In this paper, the fluid mechanics and boundary layer theory of heat conductivity theory were applied to the study of water-rock heat transform, and the heat transform between single fissure interface and water fluid was analyzed, its hydraulic-thermal model was established, and the effect of hydraulic-thermal coupling was systematically studied.

Hydrothermal Theory and
Its Implementation 2.1. Basic Assumptions. The seepage field and thermal field of discrete fissure net interact. On the one hand, the variation of fissure seepage of rock mass could influence the transmitting and exchange of water fluid and fissure interface, which further changes the distribution of the thermal field; on the other hand, the variation of the thermal field could change the fluid viscosity and density, influencing the distribution of the crack seepage field. For convenience of establishing the noncontinuum hydraulic-thermal coupling model of rock mass, some basic assumptions are as follows.
(1) Rock mass consists of impermeable rock and rock fissures; the fissure is the channel of rock mass seepage (2) The fluid in fissure seepage cannot be compressed, and it obeys Darcy's law and cubic law (3) The heat in rock is transmitted by heat conduction; the heat of water can be transmitted by heat conduction and thermal convection, which occurs on the interface of rock and water

Hydraulic-Thermal Coupling Equation.
Hydraulic-thermal analysis systems of noncontinuum rock mass include thermal field analysis of rock-fissure water, seepage field of fissures, and hydraulic-thermal coupling parameter transmitting analysis.

Thermal Field Analysis of Water in Fissures.
Heat transmission of thermal field analysis of rock seepage can be summarized: fluid-fluid heat transmission, thermal convection of fluid-fluid, and thermal convection of rock-fluid, which is displayed in Figure 1.
As regard to thermal convection of fluid-fluid, the heat flux variation of unit thickness on the x − y plane can be denoted as The heat flux variation is The heat flux variation of rock-fluid thermal convection can be expressed as where ρ w is the fluid density; C vw is the fluid heat capacity; v x , v y is the flow velocity in the x, y direction; T w is the temperature of fluid; λ w is the fluid heat conductivity; h is the fluid thermal coefficient; and T rb , T rb ′ is the boundary temperature near fissure.

Geofluids
Heat conduction coefficient h is related to characteristics of fissure interface, fluid speed, viscosity, and boundary condition. Midoux reckoned that the heat conduction coefficient h can be expressed as follows: where v i is the fluid speed along the fissure direction; a is the width of crack; L is the characteristic length related to the dimension of the problem; N u L is the Nusselt coefficient; and G Z L is the Graetz coefficient. Figure 1 illustrates the heat entry and exit diagram of fluid in a fissure. The heat equilibrium of a fluid unit can be expressed as It obeys Darcy's law in fissure fluid, which can be expressed as where v i is the fluid speed, q i is the flow rate, k f is the permeable coefficient, β is the fracture connectivity, and v is the flow viscosity of fluid.
Adopting the C 0 interpolation function of the Galerkin method, functional equations can be dispersed; then, the finite element discrete equation of the thermal field was attained: where C w is the thermal capacity matrix, K w is the heat conduction matrix, P w is the temperature load matrix, T w is the junction temperature matrix, and Δt is the time step, where the matrix element in Equation (11) can be expressed as follows: where N i , N j are the interpolation functions and H is the water head.
2.2.2. Thermal Field Analysis. The heat in rock is transmitted by the way of heat conduction; the heat theory can be expressed as follows:

Geofluids
where ρ r is the rock density, c pr is the specific heat at constant pressure of rock, T r is the temperature of rock, λ r is the thermal conductivity of rock, and W is the heat source sink.
C 0 interpolation function using the Galerkin method was adopted, the functional equation was dispersed, and the thermal field distribution was attained: where C is the heat capacity matrix, K is the heat conduction matrix of rock, P is the temperature load matrix of rock, T r is the junction load temperature matrix of rock, and Δt is the time step. The element matrix in Equation (6) was depicted as follows: where h is the water head, s is the specific water storage coefficient, and x ′ and y ′ are the local coordinates. C 0 interpolation function using the Galerkin method was adopted, the functional equation is dispersed, and the finite element discrete equation of the seep field is attained: where T, S, H is the transmitting matrix, storage matrix, and column vector, respectively. The corresponding unit matrix can be depicted as follows:   4 Geofluids noncontinuum rock is a nonlinear three-dimensional parabola equation, these equations consist of many coupling parameters and variables, they are hard to be solved even for a one-dimensional problem, and we complied the numerical simulation program FRHT3D of fluid-coupling of noncontinuum rock mass. A planar four-node unit was used for the planar fissure; eight nodes in space are used for the rock. The solving systems contain two subsystems: fissure seepage and rock system and thermal seepage field system. The main solving procedure is illustrated in Figure 2.

Numerical Simulation Model.
To study the hydrothermal coupling effect of noncontinuum rocks, the calculation example is established in Figure 3. The size of the numerical simulation model is x × y × z = 47:5 m × 6 m × 15 m, a total of 1680 grid points and 1254 coupling zones are included in the numerical simulation model, the fissure is placed z = 3 m, it is parallel with the bottom of the numerical model, and the hydraulic opening of fissure is a = 50 μm. The initial and boundary condition of seepage is that the fissure is kept with the state of saturation, its initial water head is h 0 = 20 m, the water head of the left numerical model is h 1 = 300 m, it is kept at this state during the numerical simulation, and the right water head is h 2 = 20 m.
The initial temperature condition is T 0 = 293 K, the left side of the numerical model is injected with cold water with temperature of T 1 = 293 K, and the boundary condition of the rock interface and water interface is as follows: The basic mechanical parameters of the hydraulicthermal coupling model are listed in Table 1. During the numerical simulation, the water head and temperature 1-8# were recorded at different time points.

Geofluids
To study the hydraulic-thermal coupling effect of noncontinuum rock mass, the coupling and uncoupling numerical simulation would be analyzed; the uncoupling parameters are listed in Table 2.

Hydraulic-Thermal Coupling Effect
3.2.1. Coupling Effect of Thermal Field to Seepage Field. The coupling effect of the thermal field to the seepage field is mainly due to the fact that fluid density, viscosity, and permeable coefficient are the function of fluid temperature; the variation of fluid temperature leads to the variation of seepage control equations, which influence the distribution of the fluid field. Figure 4 gives the seepage field distribution when t 1 = 600 s and t 2 = 20000 s; the solid line represents the coupling calculation results, while the dotted line represents the uncoupling numerical simulation results. From the analysis of Figure 4, it could be concluded that there exists a big difference for numerical simulation results when the coupling effects are considered or not. By comparison of Figures 4(a) and 4(b), the water head increased faster when the coupling effects are considered; the permeable water head solution is much larger than that of the uncoupling solution at the unstable seepage stage. With the development of seepage, the difference of coupling solution and uncoupling solution becomes smaller; in particular, the calculation results are close when t 2 = 20000 s. Figure 5 illustrates the seepage water head evolution difference of nodes 2#, 4#, and 7# when the coupling effect was considered or not; it could be concluded that the coupling results and uncoupling results show the same increased trend; at the initial seepage stage, the increased speed of coupling solution is faster than the later, the time for seepage tending to a constant value. When the coupling effect is not considered, the water temperature is kept at 293 K; when the coupling effect is considered, the temperature of fissure water changed with the development of seepage; during the whole coupling numerical simulation, the temperature of water flow is higher than 293 K when the coupling effect is considered. The higher the temperature is, the larger the viscosity is, which results in the increase of the permeable coefficient. Therefore, the flow speed is faster when the coupling effect is considered at the unstable seepage stage, and the seepage coupling solution is larger than that of uncoupling.

Coupling Effect of Seepage Field to Thermal Field.
The temperature of injected water (T 1 = 293 K) is lower than that of rock, the heat of rock is transmitted to water, and the flow water reduces the heat, which further changes the thermal field distribution of rock. The influence of the seepage field to the thermal field is mainly caused by two aspects: (1) thermal forced convection of the seepage flow and (2) thermal convection of seepage flow and rock surface. With the increase of water flow speed, two kinds of thermal convection accelerate both. Due to the increase in water flow speed, the heat transfer coefficient h increases; the influence of the seepage field coupling effect could be more apparent. At the unstable seepage stage, water flow is larger when the hydraulic-thermal coupling effect is considered; the cooling effect of water to rock is more significant when the coupling effect is considered, which can be observed in Figure 6; there exists an obvious difference in the thermal field when coupling effect is considered or not for t 1 = 600 s and t 2 = 20000 s; the thermal field is lower when the coupling effect is considered. Meanwhile, the temperature decreased     7 Geofluids gradient of every node, the flow speed increases; then, the heat conductivity increases subsequently; with the decrease of permeable coefficient, the flow speed slows down and the heat conductivity decreases. At the unstable stage, 2# is quite close to the location of water injection; the hydraulic gradient of 2# is larger than that of 4# and 7#. Due to the fact that flow speed of 2# is higher than 4# and 7#, the heat conductivity h of 2# is larger than that of 4# and 7#. With the development of seepage, the seepage tends to be stable, the rate of flow for the nodes tends to be a constant value, the heat conductivity all tends to be 2:35 × 10 4 W/m 2 /K, and the heat conductivity h of 2#, 4#, and 7# under the action of hydraulic-thermal coupling could be shown in Figure 7.

Conclusions
The hydraulic-thermal coupling phenomenon is widely existing in rock mass engineering,to better predict the rock mass engineering stability, the hydraulic-thermal coupling solving methods were proposed, based on the methods, the corresponding program was complied; finally, the numerical simulation of hydraulic-thermal coupling was performed. The results of this paper can be summarized as below.
(1) The influence between the seepage field and the thermal field was analyzed; meanwhile, the corresponding theory was applied to implement the influence between the seepage field and the thermal field. Meanwhile, the fortran program FRHT3 was complied to simulate the hydraulic-thermal coupling (2) Through the numerical simulation, it can be concluded that the flow speed is faster when the coupling effect is considered at the unstable seepage stage, and the seepage coupling solution is larger than that of uncoupling. Furthermore, when the coupling effects are considered, the permeable water head solution is much larger than that of the uncoupling solution at the unstable seepage stage

Data Availability
All data and models generated or used during the study appear in the submitted article.

Conflicts of Interest
The authors declare that they have no conflicts of interest.