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

Authors: Yongjian Pan; Huajun Wang; Qiang Liu
Journal: Geofluids
 (2020)

---

## Abstract

Hydraulic-thermal coupling is a key problem in rock mass engineering, especially in the disposal of nuclear wastes in deep rock mass. To accurately describe the coupling action of rock mass when under the interaction of hydraulic stress and thermal, the hydraulic-thermal coupling solving methods were proposed in this paper; in addition, the corresponding hydraulic-thermal coupling program FRHT3D was complied. Consequently, the numerical simulation was performed, 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. The proposed hydraulic-thermal coupling solving methods and programs can be applied to rock mass engineering practice.


---

## Body

## 1. 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-mechanical-chemical 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 free-face 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–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 coefficienth 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.


## 2. 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


### 2.2. 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.


#### 2.2.1. 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 Figure1.Figure 1
Heat equilibrium of fluid in fissure for a calculation unit.As regard to thermal convection of fluid-fluid, the heat flux variation of unit thickness on thex−y plane can be denoted as
(1)dQconv=dQ1=Q1−Q1′=−∫DρwCvwvx∂Tw∂x+vy∂Tw∂ydxdy.The heat flux variation is(2)dQcond=dQ2=Q2−Q2′=−∫Dλw∂∂x∂Tw∂x+∂∂y∂Tw∂ydxdy.The heat flux variation of rock-fluid thermal convection can be expressed as(3)dQ3a=hTrb−Twdx,(4)dQ3b=hTrb′−Twdx,(5)dQ3=dQ3a+dQ3b,where ρw is the fluid density; Cvw is the fluid heat capacity; vx,vy is the flow velocity in the x, y direction; Tw is the temperature of fluid; λw is the fluid heat conductivity; h is the fluid thermal coefficient; and Trb,Trb′ is the boundary temperature near fissure.Heat conduction coefficienth 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:
(6)h=λwNuL2a=3.549λw2aρwviCpw2a2λwL1/3=3.549λw2aGZL1/3,where vi 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; NuL is the Nusselt coefficient; and GZL is the Graetz coefficient.Figure1 illustrates the heat entry and exit diagram of fluid in a fissure. The heat equilibrium of a fluid unit can be expressed as
(7)ρwcvwdTwdt=dQ1+dQ2+dQ3,(8)cvwρw∂Tw∂t=λw∇2T−cpw⋅ρw⋅∇vi⋅Tw+hTrb−Tw.It obeys Darcy’s law in fissure fluid, which can be expressed as(9)vi=qi=kf∇p,i=x,y,z,(10)kf=βga212v,where vi is the fluid speed, qi is the flow rate, kf is the permeable coefficient, β is the fracture connectivity, and v is the flow viscosity of fluid.Adopting theC0 interpolation function of the Galerkin method, functional equations can be dispersed; then, the finite element discrete equation of the thermal field was attained:
(11)KwΔt+CwTwt+Δt=KwTwt+PwΔt,where Cw is the thermal capacity matrix, Kw is the heat conduction matrix, Pw is the temperature load matrix, Tw is the junction temperature matrix, and Δt is the time step, where the matrix element in Equation (11) can be expressed as follows:
(12)Kwije=∫Ωeλw∂Ni∂x∂Nj∂x+∂Ni∂y∂Nj∂y+∂Ni∂z∂Nj∂zdΩ+∫ΩrCpwρwkfNj∂Ni∂x∂H∂x+∂Ni∂y∂H∂y+∂Ni∂z∂H∂z,(13)Cwije=∫Ωecvw⋅ρw⋅Ni⋅Nj⋅dΩ,(14)Pwie=∫ΩeNi⋅hTrb−TwtdΩ,where Ni,Nj 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:(15)ρrcpr∂Tr∂t=λr∇2Tr+W,where ρr is the rock density, cpr is the specific heat at constant pressure of rock, Tr is the temperature of rock, λr is the thermal conductivity of rock, and W is the heat source sink.C0 interpolation function using the Galerkin method was adopted, the functional equation was dispersed, and the thermal field distribution was attained:
(16)KΔt+CTrt+Δt=KTrt+PΔt,where C is the heat capacity matrix, K is the heat conduction matrix of rock, P is the temperature load matrix of rock, Tr is the junction load temperature matrix of rock, and Δt is the time step. The element matrix in Equation (6) was depicted as follows:
(17)Kije=∫Ωeλr∂Ni∂x∂Nj∂x+∂Ni∂y∂Nj∂y+∂Ni∂z∂Nj∂zdΩ,(18)Cije=∫ΩecprρrNiNjdΩ,(19)Pie=∫ΩeW⋅Ni⋅dΩ+∫Γ2eq⋅Ni⋅dΓ+∫Γ3eh⋅Tra⋅Ni⋅dΓ.


#### 2.2.3. Seepage Field Analysis in Fissure

Based on the law of conservation of mass, the conservation of mass of any unit controlled by fissure can be determined by the following equation:(20)∂∂x′kf∂h∂x′+∂∂y′kf∂h∂y′=s∂h∂t,where h is the water head, s is the specific water storage coefficient, and x′ and y′ are the local coordinates.C0 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:
(21)TΔt+SPt+Δt=SPt+HΔt,where T,S,H is the transmitting matrix, storage matrix, and column vector, respectively. The corresponding unit matrix can be depicted as follows:
(22)Tije=a∫Ωekf∂Ni∂x′∂Nj∂x′+∂Ni∂y′∂Nj∂y′dΩ,(23)Sije=s∫ΩeNiNjdΩ,(24)Hie=∫Γ2eg⋅Ni⋅dΓ.


#### 2.2.4. Fluid-Thermal Coupling Coefficient Transmitting Analysis

With respect to the fluid-thermal coupling analysis of noncontinuum, the main coupling parameters, such as fluid densityρw, viscosity v, seepage coefficient kf, and heat conductivity h could be determined:
(25)ρw=−0.0035T2−0.0817T+1001,(26)v=0.1853T−0.8962,(27)kf=βgb212×0.1853T−0.8962,(28)h=3.549λw2aρwvicpw2a2λwL1/3.


### 2.3. Program for Fluid-Thermal Coupling of Noncontinuum Rock

The fluid-thermal mathematical model of 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 Figure2.Figure 2
Procedure of fluid-thermal coupling.


## 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


## 2.2. 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.


### 2.2.1. 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 Figure1.Figure 1
Heat equilibrium of fluid in fissure for a calculation unit.As regard to thermal convection of fluid-fluid, the heat flux variation of unit thickness on thex−y plane can be denoted as
(1)dQconv=dQ1=Q1−Q1′=−∫DρwCvwvx∂Tw∂x+vy∂Tw∂ydxdy.The heat flux variation is(2)dQcond=dQ2=Q2−Q2′=−∫Dλw∂∂x∂Tw∂x+∂∂y∂Tw∂ydxdy.The heat flux variation of rock-fluid thermal convection can be expressed as(3)dQ3a=hTrb−Twdx,(4)dQ3b=hTrb′−Twdx,(5)dQ3=dQ3a+dQ3b,where ρw is the fluid density; Cvw is the fluid heat capacity; vx,vy is the flow velocity in the x, y direction; Tw is the temperature of fluid; λw is the fluid heat conductivity; h is the fluid thermal coefficient; and Trb,Trb′ is the boundary temperature near fissure.Heat conduction coefficienth 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:
(6)h=λwNuL2a=3.549λw2aρwviCpw2a2λwL1/3=3.549λw2aGZL1/3,where vi 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; NuL is the Nusselt coefficient; and GZL is the Graetz coefficient.Figure1 illustrates the heat entry and exit diagram of fluid in a fissure. The heat equilibrium of a fluid unit can be expressed as
(7)ρwcvwdTwdt=dQ1+dQ2+dQ3,(8)cvwρw∂Tw∂t=λw∇2T−cpw⋅ρw⋅∇vi⋅Tw+hTrb−Tw.It obeys Darcy’s law in fissure fluid, which can be expressed as(9)vi=qi=kf∇p,i=x,y,z,(10)kf=βga212v,where vi is the fluid speed, qi is the flow rate, kf is the permeable coefficient, β is the fracture connectivity, and v is the flow viscosity of fluid.Adopting theC0 interpolation function of the Galerkin method, functional equations can be dispersed; then, the finite element discrete equation of the thermal field was attained:
(11)KwΔt+CwTwt+Δt=KwTwt+PwΔt,where Cw is the thermal capacity matrix, Kw is the heat conduction matrix, Pw is the temperature load matrix, Tw is the junction temperature matrix, and Δt is the time step, where the matrix element in Equation (11) can be expressed as follows:
(12)Kwije=∫Ωeλw∂Ni∂x∂Nj∂x+∂Ni∂y∂Nj∂y+∂Ni∂z∂Nj∂zdΩ+∫ΩrCpwρwkfNj∂Ni∂x∂H∂x+∂Ni∂y∂H∂y+∂Ni∂z∂H∂z,(13)Cwije=∫Ωecvw⋅ρw⋅Ni⋅Nj⋅dΩ,(14)Pwie=∫ΩeNi⋅hTrb−TwtdΩ,where Ni,Nj 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:(15)ρrcpr∂Tr∂t=λr∇2Tr+W,where ρr is the rock density, cpr is the specific heat at constant pressure of rock, Tr is the temperature of rock, λr is the thermal conductivity of rock, and W is the heat source sink.C0 interpolation function using the Galerkin method was adopted, the functional equation was dispersed, and the thermal field distribution was attained:
(16)KΔt+CTrt+Δt=KTrt+PΔt,where C is the heat capacity matrix, K is the heat conduction matrix of rock, P is the temperature load matrix of rock, Tr is the junction load temperature matrix of rock, and Δt is the time step. The element matrix in Equation (6) was depicted as follows:
(17)Kije=∫Ωeλr∂Ni∂x∂Nj∂x+∂Ni∂y∂Nj∂y+∂Ni∂z∂Nj∂zdΩ,(18)Cije=∫ΩecprρrNiNjdΩ,(19)Pie=∫ΩeW⋅Ni⋅dΩ+∫Γ2eq⋅Ni⋅dΓ+∫Γ3eh⋅Tra⋅Ni⋅dΓ.


### 2.2.3. Seepage Field Analysis in Fissure

Based on the law of conservation of mass, the conservation of mass of any unit controlled by fissure can be determined by the following equation:(20)∂∂x′kf∂h∂x′+∂∂y′kf∂h∂y′=s∂h∂t,where h is the water head, s is the specific water storage coefficient, and x′ and y′ are the local coordinates.C0 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:
(21)TΔt+SPt+Δt=SPt+HΔt,where T,S,H is the transmitting matrix, storage matrix, and column vector, respectively. The corresponding unit matrix can be depicted as follows:
(22)Tije=a∫Ωekf∂Ni∂x′∂Nj∂x′+∂Ni∂y′∂Nj∂y′dΩ,(23)Sije=s∫ΩeNiNjdΩ,(24)Hie=∫Γ2eg⋅Ni⋅dΓ.


### 2.2.4. Fluid-Thermal Coupling Coefficient Transmitting Analysis

With respect to the fluid-thermal coupling analysis of noncontinuum, the main coupling parameters, such as fluid densityρw, viscosity v, seepage coefficient kf, and heat conductivity h could be determined:
(25)ρw=−0.0035T2−0.0817T+1001,(26)v=0.1853T−0.8962,(27)kf=βgb212×0.1853T−0.8962,(28)h=3.549λw2aρwvicpw2a2λwL1/3.


## 2.2.1. 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 Figure1.Figure 1
Heat equilibrium of fluid in fissure for a calculation unit.As regard to thermal convection of fluid-fluid, the heat flux variation of unit thickness on thex−y plane can be denoted as
(1)dQconv=dQ1=Q1−Q1′=−∫DρwCvwvx∂Tw∂x+vy∂Tw∂ydxdy.The heat flux variation is(2)dQcond=dQ2=Q2−Q2′=−∫Dλw∂∂x∂Tw∂x+∂∂y∂Tw∂ydxdy.The heat flux variation of rock-fluid thermal convection can be expressed as(3)dQ3a=hTrb−Twdx,(4)dQ3b=hTrb′−Twdx,(5)dQ3=dQ3a+dQ3b,where ρw is the fluid density; Cvw is the fluid heat capacity; vx,vy is the flow velocity in the x, y direction; Tw is the temperature of fluid; λw is the fluid heat conductivity; h is the fluid thermal coefficient; and Trb,Trb′ is the boundary temperature near fissure.Heat conduction coefficienth 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:
(6)h=λwNuL2a=3.549λw2aρwviCpw2a2λwL1/3=3.549λw2aGZL1/3,where vi 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; NuL is the Nusselt coefficient; and GZL is the Graetz coefficient.Figure1 illustrates the heat entry and exit diagram of fluid in a fissure. The heat equilibrium of a fluid unit can be expressed as
(7)ρwcvwdTwdt=dQ1+dQ2+dQ3,(8)cvwρw∂Tw∂t=λw∇2T−cpw⋅ρw⋅∇vi⋅Tw+hTrb−Tw.It obeys Darcy’s law in fissure fluid, which can be expressed as(9)vi=qi=kf∇p,i=x,y,z,(10)kf=βga212v,where vi is the fluid speed, qi is the flow rate, kf is the permeable coefficient, β is the fracture connectivity, and v is the flow viscosity of fluid.Adopting theC0 interpolation function of the Galerkin method, functional equations can be dispersed; then, the finite element discrete equation of the thermal field was attained:
(11)KwΔt+CwTwt+Δt=KwTwt+PwΔt,where Cw is the thermal capacity matrix, Kw is the heat conduction matrix, Pw is the temperature load matrix, Tw is the junction temperature matrix, and Δt is the time step, where the matrix element in Equation (11) can be expressed as follows:
(12)Kwije=∫Ωeλw∂Ni∂x∂Nj∂x+∂Ni∂y∂Nj∂y+∂Ni∂z∂Nj∂zdΩ+∫ΩrCpwρwkfNj∂Ni∂x∂H∂x+∂Ni∂y∂H∂y+∂Ni∂z∂H∂z,(13)Cwije=∫Ωecvw⋅ρw⋅Ni⋅Nj⋅dΩ,(14)Pwie=∫ΩeNi⋅hTrb−TwtdΩ,where Ni,Nj 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:(15)ρrcpr∂Tr∂t=λr∇2Tr+W,where ρr is the rock density, cpr is the specific heat at constant pressure of rock, Tr is the temperature of rock, λr is the thermal conductivity of rock, and W is the heat source sink.C0 interpolation function using the Galerkin method was adopted, the functional equation was dispersed, and the thermal field distribution was attained:
(16)KΔt+CTrt+Δt=KTrt+PΔt,where C is the heat capacity matrix, K is the heat conduction matrix of rock, P is the temperature load matrix of rock, Tr is the junction load temperature matrix of rock, and Δt is the time step. The element matrix in Equation (6) was depicted as follows:
(17)Kije=∫Ωeλr∂Ni∂x∂Nj∂x+∂Ni∂y∂Nj∂y+∂Ni∂z∂Nj∂zdΩ,(18)Cije=∫ΩecprρrNiNjdΩ,(19)Pie=∫ΩeW⋅Ni⋅dΩ+∫Γ2eq⋅Ni⋅dΓ+∫Γ3eh⋅Tra⋅Ni⋅dΓ.


## 2.2.3. Seepage Field Analysis in Fissure

Based on the law of conservation of mass, the conservation of mass of any unit controlled by fissure can be determined by the following equation:(20)∂∂x′kf∂h∂x′+∂∂y′kf∂h∂y′=s∂h∂t,where h is the water head, s is the specific water storage coefficient, and x′ and y′ are the local coordinates.C0 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:
(21)TΔt+SPt+Δt=SPt+HΔt,where T,S,H is the transmitting matrix, storage matrix, and column vector, respectively. The corresponding unit matrix can be depicted as follows:
(22)Tije=a∫Ωekf∂Ni∂x′∂Nj∂x′+∂Ni∂y′∂Nj∂y′dΩ,(23)Sije=s∫ΩeNiNjdΩ,(24)Hie=∫Γ2eg⋅Ni⋅dΓ.


## 2.2.4. Fluid-Thermal Coupling Coefficient Transmitting Analysis

With respect to the fluid-thermal coupling analysis of noncontinuum, the main coupling parameters, such as fluid densityρw, viscosity v, seepage coefficient kf, and heat conductivity h could be determined:
(25)ρw=−0.0035T2−0.0817T+1001,(26)v=0.1853T−0.8962,(27)kf=βgb212×0.1853T−0.8962,(28)h=3.549λw2aρwvicpw2a2λwL1/3.


## 2.3. Program for Fluid-Thermal Coupling of Noncontinuum Rock

The fluid-thermal mathematical model of 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 Figure2.Figure 2
Procedure of fluid-thermal coupling.


## 3. Fluid-Thermal Coupling Effect of Noncontinuum Rock Mass

### 3.1. Numerical Simulation Model

To study the hydrothermal coupling effect of noncontinuum rocks, the calculation example is established in Figure3. The size of the numerical simulation model is x×y×z=47.5m×6m×15m, a total of 1680 grid points and 1254 coupling zones are included in the numerical simulation model, the fissure is placed z=3m, it is parallel with the bottom of the numerical model, and the hydraulic opening of fissure is a=50μm.Figure 3
Numerical simulation model.The initial and boundary condition of seepage is that the fissure is kept with the state of saturation, its initial water head ish0=20m, the water head of the left numerical model is h1=300m, it is kept at this state during the numerical simulation, and the right water head is h2=20m.The initial temperature condition isT0=293K, the left side of the numerical model is injected with cold water with temperature of T1=293K, and the boundary condition of the rock interface and water interface is as follows:
(29)−λw∂Tw∂nΓx,y=hTrb−Tw.The basic mechanical parameters of the hydraulic-thermal coupling model are listed in Table1. During the numerical simulation, the water head and temperature 1-8# were recorded at different time points.Table 1
Basic parameters of the fluid-thermal coupling numerical model.
MediumDensity (kg/m3)Thermal conductivity (W/m/K)Heat absorption capacity coefficient (J/kg/K)Initial permeable coefficient (m/s)Hydraulic opening (μm)Specific water storage coefficient (m-1)Rock2.7×1032.581.0×103Fluid10000.64.0×1032.23×10−4500.007To 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 Table2.Table 2
Calculation parameters of the numerical model without considering hydraulic-thermal coupling.
Water temperature (K)Density (kg/m3)Viscosity (m2/s)Permeable coefficient (m/s)Heat transform coefficient (W/m2/K)2931.0×1030.885×10−60.95×10−42.25×104


### 3.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. Figure4 gives the seepage field distribution when t1=600s and t2=20000s; 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 t2=20000s.Figure 4
Isometric line of the seepage field in the fissure interface.
(a)
t1=600s(b)
t2=20000s


### 3.3. Coupling Effect of Seepage Field to Thermal Field

Figure5 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.Figure 5
Water head evolution at the nodes 2#, 4#, and 7# in the fissure interface.The temperature of injected water (T1=293K) 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 t1=600s and t2=20000s; the thermal field is lower when the coupling effect is considered. Meanwhile, the temperature decreased remarkably near the rock, and the temperature increases towards the rock internal.Figure 6
Comparison of the thermal field inx‐z plane when the coupling effect is considered or not at different times.
(a)
t1=600s(b)
t2=20000s


#### 3.3.1. Coupling Effect of Calculation Parameters

The coupling parameters of hydraulic-thermal coupling of noncontinuum rock mass and the main coupling parameters include fluid densityρw, viscosity v, permeable coefficients kfi, and heat conductivity h.As regards heat conductivityh within a short time at the initial stage of seepage, due to the increase of hydraulic 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×104W/m2/K, and the heat conductivity h of 2#, 4#, and 7# under the action of hydraulic-thermal coupling could be shown in Figure 7.Figure 7
Heat conductivityh evolution of 2#, 4#, and 7# under the condition of hydraulic-thermal coupling.


## 3.1. Numerical Simulation Model

To study the hydrothermal coupling effect of noncontinuum rocks, the calculation example is established in Figure3. The size of the numerical simulation model is x×y×z=47.5m×6m×15m, a total of 1680 grid points and 1254 coupling zones are included in the numerical simulation model, the fissure is placed z=3m, it is parallel with the bottom of the numerical model, and the hydraulic opening of fissure is a=50μm.Figure 3
Numerical simulation model.The initial and boundary condition of seepage is that the fissure is kept with the state of saturation, its initial water head ish0=20m, the water head of the left numerical model is h1=300m, it is kept at this state during the numerical simulation, and the right water head is h2=20m.The initial temperature condition isT0=293K, the left side of the numerical model is injected with cold water with temperature of T1=293K, and the boundary condition of the rock interface and water interface is as follows:
(29)−λw∂Tw∂nΓx,y=hTrb−Tw.The basic mechanical parameters of the hydraulic-thermal coupling model are listed in Table1. During the numerical simulation, the water head and temperature 1-8# were recorded at different time points.Table 1
Basic parameters of the fluid-thermal coupling numerical model.
MediumDensity (kg/m3)Thermal conductivity (W/m/K)Heat absorption capacity coefficient (J/kg/K)Initial permeable coefficient (m/s)Hydraulic opening (μm)Specific water storage coefficient (m-1)Rock2.7×1032.581.0×103Fluid10000.64.0×1032.23×10−4500.007To 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 Table2.Table 2
Calculation parameters of the numerical model without considering hydraulic-thermal coupling.
Water temperature (K)Density (kg/m3)Viscosity (m2/s)Permeable coefficient (m/s)Heat transform coefficient (W/m2/K)2931.0×1030.885×10−60.95×10−42.25×104


## 3.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. Figure4 gives the seepage field distribution when t1=600s and t2=20000s; 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 t2=20000s.Figure 4
Isometric line of the seepage field in the fissure interface.
(a)
t1=600s(b)
t2=20000s


## 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. Figure4 gives the seepage field distribution when t1=600s and t2=20000s; 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 t2=20000s.Figure 4
Isometric line of the seepage field in the fissure interface.
(a)
t1=600s(b)
t2=20000s


## 3.3. Coupling Effect of Seepage Field to Thermal Field

Figure5 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.Figure 5
Water head evolution at the nodes 2#, 4#, and 7# in the fissure interface.The temperature of injected water (T1=293K) 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 t1=600s and t2=20000s; the thermal field is lower when the coupling effect is considered. Meanwhile, the temperature decreased remarkably near the rock, and the temperature increases towards the rock internal.Figure 6
Comparison of the thermal field inx‐z plane when the coupling effect is considered or not at different times.
(a)
t1=600s(b)
t2=20000s


### 3.3.1. Coupling Effect of Calculation Parameters

The coupling parameters of hydraulic-thermal coupling of noncontinuum rock mass and the main coupling parameters include fluid densityρw, viscosity v, permeable coefficients kfi, and heat conductivity h.As regards heat conductivityh within a short time at the initial stage of seepage, due to the increase of hydraulic 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×104W/m2/K, and the heat conductivity h of 2#, 4#, and 7# under the action of hydraulic-thermal coupling could be shown in Figure 7.Figure 7
Heat conductivityh evolution of 2#, 4#, and 7# under the condition of hydraulic-thermal coupling.


## 3.3.1. Coupling Effect of Calculation Parameters

The coupling parameters of hydraulic-thermal coupling of noncontinuum rock mass and the main coupling parameters include fluid densityρw, viscosity v, permeable coefficients kfi, and heat conductivity h.As regards heat conductivityh within a short time at the initial stage of seepage, due to the increase of hydraulic 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×104W/m2/K, and the heat conductivity h of 2#, 4#, and 7# under the action of hydraulic-thermal coupling could be shown in Figure 7.Figure 7
Heat conductivityh evolution of 2#, 4#, and 7# under the condition of hydraulic-thermal coupling.


## 4. 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



---

*Source: 8830862-2020-11-23.xml*