Next Article in Journal
Evaluating the Use of Displacement Ventilation for Providing Space Heating in Unoccupied Periods Using Laboratory Experiments, Field Tests and Numerical Simulations
Next Article in Special Issue
Stage Division of Landslide Deformation and Prediction of Critical Sliding Based on Inverse Logistic Function
Previous Article in Journal
Visualization of Chemical Heavy Oil EOR Displacement Mechanisms in a 2D System
Previous Article in Special Issue
Experimental Study on Infrared Temperature Characteristics and Failure Modes of Marble with Prefabricated Holes under Uniaxial Compression
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Novel Procedure for Coupled Simulation of Thermal and Fluid Flow Models for Rough-Walled Rock Fractures

1
Faculty of Engineering, China University of Geosciences, Wuhan 430074, China
2
School of Civil Engineering, Wuhan University, Wuhan 430072, China
3
School of Earth Sciences and Engineering, Hohai University, Nanjing 210098, China
*
Author to whom correspondence should be addressed.
Energies 2021, 14(4), 951; https://doi.org/10.3390/en14040951
Submission received: 23 December 2020 / Revised: 6 February 2021 / Accepted: 8 February 2021 / Published: 11 February 2021

Abstract

:
An enhanced geothermal system (EGS) proposed on the basis of hot dry rock mining technology has become a focus of geothermal research. A novel procedure for coupled simulation of thermal and fluid flow models (NPCTF) is derived to model heat flow and thermal energy absorption characteristics in rough-walled rock fractures. The perturbation method is used to calculate the pressure and flow rate in connected wedge-shaped cells at pore-scale, and an approximate analytical solution of temperature distribution in wedge-shaped cells is obtained, which assumes an identical temperature between the fluid and fracture wall. The proposed method is verified in Barton and Choubey (1985) fracture profiles. The maximum deviation of temperature distribution between the proposed method and heat flow simulation is 13.2% and flow transmissivity is 1.2%, which indicates the results from the proposed method are in close agreement with those obtained from simulations. By applying the proposed NPCTF to real rock fractures obtained by a 3D stereotopometric scanning system, its performance was tested against heat flow simulations from a COMSOL code. The mean discrepancy between them is 1.51% for all cases of fracture profiles, meaning that the new model can be applicable for fractures with different fracture roughness. Performance analysis shows small fracture aperture increases the deviation of NPCTF, but this decreases for a large aperture fracture. The accuracy of the NPCTF is not sensitive to the size of the mesh.

Graphical Abstract

1. Introduction

Geothermal energy, which is clean, renewable, and widespread, has been developed and utilized in many countries, and considerable attention has been given to exploiting thermal energy from hot dry rock (HDR) 3–10 km underground, using the enhanced geothermal system (EGS) to exchange and loop heat in/out of effectively connected systems of underground fractures [1,2,3,4,5,6,7,8,9,10,11,12,13,14]. For example, the Landau plant in Germany [15], the Soultz plant in France [16], and the recently commissioned Geodynamics Habanero pilot plant [17] demonstrate the feasibility of EGS running. Figure 1 shows the traditional heat extraction process from the geothermal system, where the vast rock fracture system is the main path for heat fluid flow. However, prediction of EGS service life and thermal energy extraction encounters considerable challenges due to limited knowledge of the heat exchange behavior between flowing fluid and high temperature rock [18,19,20,21].
Natural rock fractures are characterized by irregular shape, rough surface, and even contact asperities, which play important roles in the hydraulic and heat transfer properties of fractures. Because of the rough fracture surface, there may be a 1–2 order of magnitude error of permeability calculated by the widely used cubic law, which assumes the fracture is a smooth parallel plate [22,23,24]. Zou et al. [25] and Wang et al. [26] used a 3D Lattice Boltzmann method to predict the fluid flow behavior in fractures and found that the secondary roughness enhanced the local complexity of velocity distribution and led to nonlinear flow behavior. Similar studies were also reported by Yin et al. [27] and Xiong et al. [24]. In addition, the rough fracture surface also causes a decrease in heat transport intensity. Recently, many efforts [28,29,30] paid attention to the study of the effects of fracture roughness on heat transport behaviors in fractures, based on laboratory experiments and numerical simulations. For example, He et al. [28] found that heat transport characteristics mainly depended on fracture surface roughness, followed by aperture and flow rate, combining the experimental and numerical modeling approaches. Huang et al. [29] conducted seepage and convective heat transport experiments, and demonstrated that large roughness in the direction perpendicular to flow would increase the capacity of heat transport. Ma et al. [30] adopted numerical and experimental approaches to improve understanding of the heat transfer characteristics of water flowing through rough fractures, and results showed the total heat transfer coefficient increased with the increase in the value of the joint roughness coefficient (JRC). However, there are no analytical models describing the heat transport properties between fluid and high temperature rock, considering fracture surface geometry. In previous research studies, the parallel plate is traditionally used to simplify the real rough fracture surface, and many models can be observed in the literature, such as the Gehlin and Hellstrom model [2], Gringarten et al. model [31], Cheng et al. model [32], Martinez et al. model [33], Zhao model [34], Yost and Einstein model [35], and Yan and Jiao model [7] which have inevitable relative error comparing with the real heat transport process in rough-walled fractures.
In view of analytical heat transport models for rough fractures being rare, this paper aims at developing a proposed novel procedure for coupled simulation of thermal and fluid flow (NPCTF) models to understand the heat transport process in rough fractures, based on representing the measured fracture geometry with a series of connecting wedges formed using adjacent apertures in longitude direction along the fracture plane, which is a common approximation approach for fracture geometry [36]. The performance of the proposed NPCTF is assessed by comparing predictions with results from numerically solving the Navier–Stokes and energy equations, using COMSOL code in rock fracture profiles obtained by a 3D stereotopometric scanning system, which have different relative surface roughness and aperture.

2. Fluid Flow Model

Flow behavior in natural rock fractures has been studied through considering fracture walls as parallel plates, saw-tooth shaped walls, and sinusoidally varying walls [30]. Among these approximates, fractures are found to be well described by a series of connected wedges. Wang et al. [37,38] used a wedge-shaped fracture to obtain the perturbation solution of pressure ∆p and flowrate Q under the pressure boundary condition. The perturbation parameter ε is selected as the relative aperture variation along the wedge length l in his study. This section provides a brief description of the approach rather than a detailed derivation. A more detailed perturbation derivation can be viewed in Wang et al. [37,38]. By perturbation analysis, the pressure difference can be made dimensionless and expressed as expanded series with a small parameter ε:
Δ P = Δ P 0 + ϵ Δ P 1 + O ϵ 2
Δ P 0 = 1 ω 0 ω 1 B 3 d X
Δ P 1 = 1 ω 0 ω 9 R B 70 B 3 d X
where X is defined as X = εx/bm; ω is the dimensionless absolute aperture variation defined as ω = |a|/bm, where a is the difference between the upper and lower wedge edge; B is the dimensionless aperture defined as B = c/cm, c is the fracture local aperture, and cm is the mean aperture; B′ is the first derivatives of B with respect to X; R is the Reynolds number defined as R = ρQ/μ; Q is flow rate per unit width of fracture and μ is fluid viscosity; ∆P is the dimensionless pressure difference given as ∆P = ∆p/∆pm and ∆pm is the pressure difference of flow through a fracture with a uniform aperture defined from the cubic law as:
Δ p m = 12 μ l Q b m 3
Noted that only the first two order terms in Equation (1) are used which is enough accuracy for nonlinear flow simulation. Because the quadratic polynomial form, called the Forchheimer law, can well describe nonlinear flow behavior in a single fracture [39,40], the highest order term is the second order term.

3. Heat Transport Model

The steady heat fluid transport differential equation in the fracture, including advection, conduction, and convection terms for the fracture walls, can be expressed as [34]:
v T f x , y x K w ρ w c w 2 T f x , y 2 x 2 K r ρ w c w b T r x , y y | b o u n d a r y = 0
where ρw is the fluid density; cw is the specific heat of fluid; Kw is the fluid thermal conductivity; Kr is the rock thermal conductivity; v is the steady flow velocity. b is half aperture of the fracture; Tf is the bulk temperature of the fluid; Tr is the temperature of the reservoir rock matrix. In this equation, the temperature at the rock fracture surface is identified as that of the bulk fluid. The steady heat conduction in the rock is governed by [32,34]:
T 2 r x , y y 2 = 0
which assumes that the heat conduction is two dimensional, perpendicular to the fracture plane.
The boundary conditions associated with the heat flow in the fracture are:
T f 0 = T i n
T f = T 0
T r x , R = T 0
where R is the height of the rock. The identical temperature between the fluid and fracture surfaces is expressed as:
T r x , f t x = T f x
T r x , f l x = T f x
where ft(x) and fl(x) are the upper fracture wall and lower fracture wall, respectively. Zhao et al. [34] used this assumption to obtain a solution of temperature in the plane fracture. This solution has good agreement with the experiment results. However, there are not any temperature solutions for the fracture wedge, which is the basic element in rough fractures. Therefore, a solution for the fracture wedge would be developed based on the above heat fluid transport equations.
(1)
Symmetric fracture wedge
For the symmetric fracture wedge (Figure 2a), the wedge wall (upper fracture ft(x), lower fracture fl(x)) can be described as:
f t x = A x + b 0
f l x = A x b 0
The aperture b can be determined as ft(x) − fl(x) according to Equations (12) and (13). If thermal diffusion in the fluid is not considered, the boundary conditions in Equations (10) and (11) are substituted into Equations (5) and (6); a temperature distribution equation of fluid can be expressed as:
v T f x , y x 2 K r ρ w c w A x + b 0 T f x T 0 A x + b 0 R = 0
The solution to Equation (14) considering the temperature boundary in Equations (7)–(9) is:
T f x = T 0 + T i n T o exp ( k r l n ( | b 0 + R b 0 + R + A x | ) A b 0 c w ρ w v )
(2)
Asymmetric fracture wedge
For the asymmetric fracture wedge (Figure 2b), the wedge wall can be described as:
f t x = A 1 x + b 1
f l x = A 2 x b 2
Similarly, the solutions to Equations (5) and (6) for the asymmetric wedge with boundaries in Equations (16) and (17) are:
T f 1 x = T 0 + T i n T 0 exp 2 k r ln ( | b 1 + R b 1 + R + A 1 x | ) A 1 b 1 c w p w v A 1 b 2 c w p w v
T f 2 x = T 0 + T i n T 0 exp 2 k r ln ( | b 2 R b 2 R + A 2 x | ) A 2 b 1 c w p w v A 2 b 2 c w p w v
T f x = T f 1 x + T f 2 x 2
The temperature mathematic solutions for the fracture wedge were first proposed in the literature based on a few assumptions which have well-defined parameters.

4. Procedure for Coupled Simulation of Thermal and Fluid Flow Models (NPCTF)

4.1. Model Description

As presented above, the fluid flow model (Equation (1)) and heat transport model (Equations (18)–(20)) were obtained in the fracture wedge. Each fracture wedge was connected to another and constituted a rough-walled fracture (Figure 3a). According to the flowrate equilibrium principle and fluid flow model (Equation (1)), the flowrate and pressure distribution in the entire fracture can be determined (seen in Figure 3b). If the flowrate distribution is known, the temperature solution in each connected fracture wedge would be obtained using the heat transport model (as Figure 3c). When the pressure and temperature condition changes, the density of water (ρw) is no longer a constant value; this can be expressed as a function of pressure and temperature [6]:
1 ρ w = 3.086 0.899017 647.25 T f 0.147166 0.39 658.15 T f 1.6 P 225.5 + δ w
where δw is a function of P and Tf, and the value of δw does not exceed 6% of 1/ρw. The temperature of water also affects the kinematic viscosity of water (μ). An empirical formula for kinematic viscosity is [6]:
μ = 0.01775 1 + 0.033 T f + 0.000221 T f 2
Therefore, the function of density and viscosity of fluid can relate the fluid flow model and heat transport model, and a procedure for coupled simulation of thermal and fluid flow (NPCTF) models in a two-dimension fracture profile was proposed. A detailed flowchart for the proposed simulation method is seen in Figure 3d. When fracture profile information (x,y) is obtained, the fluid flow model is used to evaluate flow rate Q and pressure P, and the heat transport model then predicts the temperature distribution Tf(x). An iterative process is used to solve the fluid flow model and heat transport model by modifying the fluid’s physical properties (fluid density ρw and kinematic viscosity μ) using new flowrate, pressure, and temperature before moving to the next step. The iterative process stops once the temperature and flowrate stabilize to a specified tolerance. The tolerances of flowrate and temperature are defined as ||Qj−1Qj|| ≤ 0.00001 and ||Tfj−1-Tfj|| ≤ 0.00001, respectively. j is each iteration step in the calculation process. This calculation method would increase calculated efficiency, compared to the previous finite element method or finite difference method.

4.2. Validation of NPCTF

Barton and Bandis (1985) [41] analyzed the roughness of 136 natural fracture surfaces. According to the roughness of the fracture surface, it is divided into 10 levels, of which values are in the range of 0–20. The roughness index is called the joint roughness coefficient (JRC). The JRC parameter is widely used in rock mechanics and engineering. In order to valid the NPCTF, 5 JRC profiles were selected to calculate temperature and pressure distribution, and those results were compared to numerical results with COMSOL. The JRC values of the selected fracture profiles are 1, 5, 10, 15, and 20, respectively. Table 1 lists the thermal and flow parameters for predicting the temperature in JRC profiles.
Figure 4 shows the comparison between the predicted results of the NPCTF and the calculated results of COMSOL. As illustrated in Figure 4, the calculated results of the NPCTF are consistent with the COMSOL results. With the increase in JRC value, the error becomes larger and a more obvious difference is observed at the nonlinear part of the temperature profile. However, the maximum error value is 12.3%. The reason is that the effect of thermal diffusion is not considered. With the increase in JRC, the thermal diffusion effect is more obvious, and enhances the differences in the nonlinear part of the temperature profile between the proposed NPCTF and COMSOL results. Figure 5 shows the temperature counter in the fracture profile with different JRC. In addition, the results of flow field calculated by proposed NPCTF and COMSOL are compared. The flow field is expressed in terms of transmissivity Ta, which is defined as follows:
T a = Q L P
Figure 6 shows the comparison of transmissivity Ta of fractures under different JRC values. The transmissivities predicted by the proposed NPCTF are consistent with the COMSOL results. As the JRC value increases, the error increases, and the maximum error is 1.2%.
Therefore, the proposed NPCTF can predict temperature and flow velocity distribution in rough-walled fractures.

5. Coupled Hydrothermal Simulation in Three-Dimensional Rock Fracture

5.1. Modelling Setup

To obtain the rock fractures, an intact granite block was split using the Brazilian indirect tensile test into two halves with dimensions of 150 × 150 × 75 mm. This method was also used by Xiong et al. [24] to study fluid flow behavior in rock fractures.
To determine the surface information of the fracture, an advanced 3D CaMega stereotopometric scanning system (BoWeiHengXin Inc., Beijing, China) was used to perform a non-contact 3D scan of the rough fracture surfaces. Then, the coordinates of 22,801 points of the upper half of the fracture were obtained, as well as the lower half fracture. Tatone et al. [42] introduced a novel method to measure fracture geometry based on fracture surface data. Based on this method, a 3D fracture model was obtained, which is shown in Figure 7a.
Five evenly spaced surface profiles parallel to the flow direction with lengths of 150 mm of the upper and lower half of the fracture were taken into consideration, namely slices X = 15, 45, 75, 105, and 135 mm, as shown in Figure 7a. The profile curves of the five slices were obtained by importing the coordinates (x, z) of every point of each slice to Auto CAD. Correspondingly, different values of JRC are calculated according to Tse and Cruden [43], who proposed a relationship of JRC and the root mean square of the first derivative of profile Z2:
JRC = 32.2 + 32.47 log Z 2
Table 2 lists the JRC and mean aperture of the fractures.
The heat flow simulations in this work were conducted by solving a 3D Navier–Stokes equation and an energy equation with an advanced COMSOL Multiphysics code. The mesh scheme in the COMSOL mesh module was optimized to reduce the mesh-related effect on flow simulations. The fracture computational mesh was further refined, with denser mesh in the fracture part to capture the details of the flow. The final fracture mesh is shown in Figure 7b. The mean mesh sizes of the rock matrix and fracture are 35.0 and 5.0 μm, respectively. In general, over 3.2 million tetrahedron mesh elements were constructed for the fracture model. The temperature boundary conductions were the upper and lower surface of the rock matrix. No-flux and non-slip boundary conditions were set for the fracture walls. The inflow boundary was assigned a constant mass flow rate and a constant temperature, and the outflow boundary was set as a constant pressure condition.

5.2. Temperature and Pressure Distribution Along Fracture

To demonstrate the general thermal transport within the fractures, the temperature counter and distribution are plotted in Figure 8 along five fracture profiles tested in this study using COMSOL simulation results. The calculated parameters are listed in Table 1. The corresponding injection pressure was 0.01 Pa. As illustrated in Figure 8, the water has a strong trend to attain equilibrium temperature (363.15 K) close to the fracture inlet since the coupling between the water and the rock is strong and, therefore, the heat transport is more intense. Far away from the fracture inlet, the water temperature slowly attains equilibrium temperature. If the fracture is long enough, the temperature at the fracture outlet position would be same as the equilibrium temperature. Additionally, each temperature along the rough-walled fracture is not a smooth curve; the rougher fracture surface leads to more fluctuation for the temperature curve, indicating the fracture roughness affects the temperature distribution, as shown in Figure 8. It can also be seen that as fracture aperture increases, the temperature at which water is absorbed decreases. The difference between the temperature related to the aperture parameter is readily apparent when the change in temperature that occurred between fractures with bm = 1.2 mm and bm = 1.8 mm (the smallest and largest fracture) aperture is examined. A decrease in temperature of nearly 25% was calculated between the flow in these two fractures. This relatively large change in aperture relates well to the large simulated change in temperature. Hence, both fracture roughness and aperture have significant influences on temperature distribution.
The pressure distribution of the fracture is also plotted in Figure 9. As illustrated in Figure 9, from the inlet to the outlet, the fluid pressure decreases from 1 to 0 Pa. Both fracture roughness and aperture have small influences on pressure distribution.

6. Model Performance

6.1. Comparison with Previous Models from the Literature

To demonstrate the robustness of the proposed NPCTF (the calculation process shown in Figure 3) in describing this thermal transport in realistic 3D rock fractures, the deviation D is defined as the relative error in estimating the temperature distribution:
D = T i m o d e l T i r e a l T i m o d e l × 100 %
where Ti-model is the fracture’s local temperature at i node obtained from the NPCTF and Ti-real is the fracture local temperature at i node obtained from numerical simulations. The temperature deviations D calculated by Equation (25) for the five fracture profiles are plotted in Figure 10 at different fracture positions. In general, the two temperatures predicted by NPCTF and COMSOL simulations are in close agreement for all fracture profiles, with the deviation D ranging from 0.44% to 8.10% with a mean of 1.51%. A highest mean deviation <D> of 2.57% is observed for the fracture profile with JRC = 16.34. The fracture profile with JRC = 11.96 has the lowest <D> of 0.75% in this case and the range of D is from 0.48% to 3.98%. The root mean squares (RMS) of deviation D were calculated. It is interesting that the fracture profile with the higher JRC has the higher RMS. The highest RMS of 0.0173 is observed for the fracture profile with JRC = 16.34. The fracture profile with JRC = 11.96 has the lowest RMS of 0.00858. The deviation D increases with increases in JRC, also suggesting that the prediction accuracy of the proposed NPCTF is related to fracture roughness.
Further assessments of NPCTF were made by comparing the values of D with those obtained from previously proposed heat transport models including the Gehlin and Hellstrom [2] and Zhao models [34]. A steady-state analytical solution describing the temperature distribution in the vertical fracture was developed by Gehlin and Hellstrom [2], as follows:
T f x = T o + T i n T o exp ( α c w v x )
where α is the heat transfer coefficient; the value of 900 W m−2 K−1 was used in the research of Shaik et al. [44]. To is rock temperature and Tin is water temperature in the inlet of the fracture. When the diverging angle or converging angle is 0, the proposed heat transport model can be degenerated as in the Zhao model [34]:
T f x = T o + T i n T o exp ( k r v ρ w c w b R x )
For comparison, the heat transport model took the place of the Gehlin and Hellstrom and Zhao heat transport models in the NPCTF, when applying the Gehlin and Hellstrom [2] and Zhao [34] models to calculate the temperature distribution in all five fracture cases. The mean temperature deviations <D> (see in Equation (25)) for all five fracture cases by comparing the NPCTF, the Gehlin and Hellstrom model [2], and the Zhao model [34] with COMSOL simulation are plotted in Figure 11. Overall, the highest deviation is observed for the Gehlin and Hellstrom model [2], with a mean of 23.27% and a maximum of 27.81%. The Zhao model [34] deviates by less than 10% of <D> when JRC is less than 13; however, the deviation rises to 12.45%, 12.75%, and 14.16%, respectively, as JRC increases to 16.34. The poorer performance of the two models is due to not taking fracture roughness and aperture variation into consideration, compared to the proposed NPCTF. Therefore, it can be observed that the D from the NPCTF is less than other models. Hence, the proposed NPCTF outperformed other models in different JRC with the lowest <D> of 0.44%.

6.2. Sensitivity Analysis

The accuracy of the proposed NPCTF is mainly affected by the fracture geometry (i.e., roughness and aperture variation) and mesh size (one dimensional element). The fracture roughness has been discussed above. Hence, the two factors (aperture and mesh size) need to be tested in the NPCTF.
The fracture profile with maximum roughness (JRC = 16.34 and bm = 1.5 mm) was selected to test, and aperture bm was fixed at 0.6, 0.9, 1.2, 1.5, 1.8, and 2.1 mm by moving the distance of the lower and upper fracture surface. Figure 12 plots the effect of aperture on mean deviation <D> (see in Equation (25)) from NPCTF). As shown in Figure 13, the mean deviation <D> decreases as the value of bm increases. At JRC = 16.34, as the value of bm changes from 0.6 to 2.1 mm, the mean deviation <D> decreases from 1.80 to 0.62. The deviation of cases for aperture less than 1.5 mm is larger than that of cases for an aperture value of 1.8 to 2.1 mm. It is evident that the effect of bm on <D> is larger for a small aperture fracture and smaller for a larger aperture fracture.
In order to analyze the effect of mesh size, the mesh size is set as 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, and 0.9 mm and 1.0, 1.1, and 1.2 mm for the fracture profile with JRC = 16.34 and bm = 1.5 mm, respectively. Figure 13 shows the effect of mesh size on mean deviation <D> (see in Equation (25)) from the NPCTF. As shown in Figure 13, the mean deviation <D> increases as the mesh size increases from 0.6 to 1.2 mm. The maximum deviation for a mesh size value of 1.2 mm is 1.57%. The results tend to converge when the simulation mesh size is less than 0.6 mm. The differences in deviation from using a mesh size of 0.3 and 0.6 mm are approximately 1.38% and 1.41%, respectively. Therefore, this indicates that the NPCTF is not sensitive to the size of mesh defined in the range of 0.2 to 0.8 mm.

6.3. Limitations of NPCTF

The main limitation of this proposed NPCTF from the heat transport model is neglecting the heat diffusion term. This assumption can be violated as the diverging or converging angle increases. This heat diffusion effect can only cause bigger deviation in the distance from the fracture inlet to equilibrium temperature position. Furthermore, the deviation from the proposed heat transport model is reasonable when presented a diverging wedge with an angle that is 63.47°. The feature of the diverging and converging angle is mostly considered and reported for rock fractures [36,37,45] and therefore, the validity of the assumption is warranted in most cases. Additionally, the validity of the used flow perturbation solution in the fracture lies in the assumptions that the variation of aperture along the fracture length direction is small. However, Wang et al. [37] discussed that the assumption can meet the geometry of most natural fractures.

7. Conclusions

We presented an NPCTF for rough-walled fractures, including a fluid flow model and a heat transport model. The fluid flow model describes the relationship of flowrate and pressure in the wedge-shaped cell, based on the perturbation method. The heat transport properties between the flowing fluid and rock fracture wall are captured by the proposed heat transport model, which assumes an identical temperature between the fluid and fracture wall. Both of them connected to each other in a series of connected wedge-shaped fractures. The proposed NPCTF was validated by a commercial COMSOL Multiphysics code in JRC profiles. The temperature relative deviation between them is less than 12.3% and the maximum transmissivity deviation is 1.2%. To examine the performance of the proposed NPCTF, heat flow simulations in real rock fractures with different relative surface roughness, obtained by a 3D stereotopometric scanning system, were performed. It has been shown that the results from the proposed NPCTF agree well with the COMSOL simulation results. In general, the absolute deviation of the proposed model from the simulation results, in terms of the estimated temperature, is in the range of 0.44% to 8.10% with a mean deviation of 1.51% for all cases of the fracture profiles examined. In addition, the NPCTF is sensitive to fracture aperture variation, but not to the size of the mesh. The effect of aperture on deviation of the NPCTF is larger for a small aperture fracture compared with a large aperture fracture.
The proposed NPCTF has good estimation of the heat transport effect of rough-walled fractures. Therefore, in the next work, we would take the proposed NPCTF to use in a two-dimension rock fracture network. On this basis, the heat energy exploitation process in EGS can be simulated.

Author Contributions

Conceptualization, F.X. and Q.J.; methodology, F.X.; software, C.Z.; validation, C.Z.; formal analysis, F.X.; investigation, F.X.; resources, Q.J.; data curation, C.Z.; writing—original draft preparation, F.X.; writing—review and editing, Q.J.; visualization, C.Z.; supervision, Q.J.; project administration, Q.J.; funding acquisition, Q.J. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Science Foundation of China, grant number 1679173.

Institutional Review Board Statement

No applicable.

Informed Consent Statement

No applicable.

Data Availability Statement

All data included in this study are available upon request by contact with the corresponding author.

Acknowledgments

This research was financially supported by the National Science Foundation of China (No. 51679173). This support is gratefully acknowledged.

Conflicts of Interest

The authors declared that they have no conflicts of interest to this work.

Abbreviations

εPerturbation parameter
XDimensionless roughness
ωDimensionless absolute aperture variation
aDifference between the upper and lower wedge edge
BDimensionless aperture
cFracture local aperture
cmMean aperture
BFirst derivatives of B with respect to X
RReynolds number
QFlow rate per unit width of fracture
μFluid viscosity
PDimensionless pressure difference
pmPressure difference of flow through a fracture with a uniform aperture defined from the cubic law
ρwfluid density
cwSpecific heat of fluid;
KwFluid thermal conductivity
KrRock thermal conductivity
vSteady flow velocity
bHalf aperture of the fracture
TfBulk temperature of fluid
TrTemperature of reservoir rock matrix
Torock temperature
Tinwater temperature in inlet of fracture
RHeight of rock
αInclined angle
δwFunction of P and Tf
Tatransmissivity
JRCJoint roughness coefficient
Z2root mean square of first derivative of profile
DRelative deviation between two variables
<D>Mean deviation

References

  1. Kolditz, O.; Clauser, C. Numerical simulation of flow and heat transfer in fractured crystalline rocks: Application to the hot dry rock site in Rosemanowes (UK). Geothermics 1997, 27, 1–23. [Google Scholar] [CrossRef]
  2. Gehlin, S.; Hellstrom, G. Influence on thermal response test by groundwater flow in vertical fractures in hard rock. Renew. Energy 2003, 28, 2221–2238. [Google Scholar] [CrossRef]
  3. Zhao, Y.; Feng, Z.; Feng, Z.; Yang, D.; Liang, W. THM (Thermo-hydro-mechanical) coupled mathematical model of fractured media and numerical simulation of a 3D enhanced geothermal system at 573 K and buried depth 6000-7000 M. Energy 2015, 82, 193–205. [Google Scholar] [CrossRef]
  4. Jiao, Y.Y.; Zhang, X.L.; Zhang, H.Q.; Li, H.B.; Yang, S.Q.; Li, J.C. A coupled thermo-mechanical discontinuum model for simulating rock cracking induced by temperature stresses. Comput. Geotech. 2015, 67, 142–149. [Google Scholar] [CrossRef]
  5. Sun, Z.; Zhang, X.; Yao, J.; Wang, H.; Lv, S.; Sun, Z.; Huang, Y.; Cai, M.; Huang, X. Numerical simulation of the heat extraction in EGS with thermal-hydraulic-mechanical coupling method based on discrete fractures model. Energy 2017, 120, 20–33. [Google Scholar] [CrossRef]
  6. Yao, J.; Zhang, X.; Sun, Z.; Huang, Z.; Liu, J.; Li, Y.; Xin, Y.; Yan, X.; Liu, W. Numerical simulation of the heat extraction in 3D-EGS with thermal-hydraulic-mechanical coupling method based on discrete fractures model. Geothermics 2018, 74, 19–34. [Google Scholar] [CrossRef]
  7. Yan, C.; Jiao, Y.Y. FDEM-TH3D: A three-dimensional coupled hydrothermal model for fractured rock. Int. J. Numer. Anal. Methods Geomech. 2019, 43, 415–440. [Google Scholar] [CrossRef] [Green Version]
  8. Ehyaei, M.; Ahmadi, A.; Rosen, M.A.; Davarpanah, A. Thermodynamic Optimization of a Geothermal Power Plant with a Genetic Algorithm in Two Stages. Processes 2020, 8, 1277. [Google Scholar] [CrossRef]
  9. Davarpanah, A.; Shirmohammadi, R.; Mirshekari, B.; Aslani, A. Analysis of hydraulic fracturing techniques: Hybrid fuzzy approaches. Arab. J. Geosci. 2019, 12, 402. [Google Scholar] [CrossRef]
  10. Sun, S.; Zhou, M.; Lu, W.; Davarpanah, A. Application of Symmetry Law in Numerical Modeling of Hydraulic Fracturing by Finite Element Method. Symmetry 2020, 12, 1122. [Google Scholar] [CrossRef]
  11. Zhu, M.; Yu, L.; Zhang, X.; Davarpanah, A. Application of Implicit Pressure-Explicit Saturation Method to Predict Filtrated Mud Saturation Impact on the Hydrocarbon Reservoirs Formation Damage. Mathematics 2020, 8, 1057. [Google Scholar] [CrossRef]
  12. Hu, X.; Xie, J.; Cai, W.; Wang, R.; Davarpanah, A. Thermodynamic effects of cycling carbon dioxide injectivity in shale reservoirs. J. Pet. Sci. Eng. 2020, 195, 107717. [Google Scholar] [CrossRef]
  13. Xiong, F.; Jiang, Q.; Xu, C. Fast Equivalent Micro-Scale Pipe Network Representation of Rock Fractures Obtained by Computed Tomography for Fluid Flow Simulations. Rock Mech. Rock Eng. 2020, 1–17. [Google Scholar] [CrossRef]
  14. Xiong, F.; Wei, W.; Xu, C.; Jiang, Q. Experimental and numerical investigation on nonlinear flow behaviour through three dimensional fracture intersections and fracture networks. Comput. Geotech. 2020, 121, 103446. [Google Scholar] [CrossRef]
  15. Heimlich, C.; Gourmelen, N.; Masson, F.; Schmittbuhl, J.; Kim, S.W.; Azzola, J. Uplift around the geothermal power plant of Landau (Germany) as observed by InSAR monitoring. Geotherm. Energy 2015, 3, 2. [Google Scholar] [CrossRef]
  16. Tenzer, H.; Park, C.-H.; Kolditz, O.; McDermott, C.I. Application of the geomechanical facies approach and comparison of exploration and evaluation methods used at Soultz-sous-Forêts (France) and Spa Urach (Germany) geothermal sites. Environ. Earth Sci. 2010, 61, 853–880. [Google Scholar] [CrossRef]
  17. Xu, C.; Dowd, P.; Tian, Z. A simplified coupled hydro-thermal model for enhanced geothermal systems. Appl. Energy 2015, 140, 135–145. [Google Scholar] [CrossRef]
  18. Kolditz, O. Modelling flow and heat transfer in fractured rocks: Conceptual model of a 3-D deterministic fracture network. Geothermics 1995, 24, 451–470. [Google Scholar] [CrossRef]
  19. Natarajan, N.; Xu, C.; Dowd, P.A.; Hand, M. Numerical modelling of thermal transport and quartz precipitation/dissolution in a coupled fracture–skin–matrix system. Int. J. Heat Mass Transf. 2014, 78, 302–310. [Google Scholar] [CrossRef] [Green Version]
  20. Bai, B.; He, Y.; Li, X.; Hu, S.; Huang, X.; Li, J.; Zhu, J. Local heat transfer characteristics of water flowing through a single fracture within a cylindrical granite specimen. Environ. Earth Sci. 2016, 75, 1460. [Google Scholar] [CrossRef]
  21. Li, Z.; Feng, X.; Zhang, Y.; Zhang, C.; Xu, T.; Wang, Y. Experimental research on the convection heat transfer characteristics of distilled water in manmade smooth and rough rock fractures. Energy 2017, 133, 206–218. [Google Scholar] [CrossRef]
  22. Tsang, Y. The effect of tortuosity on fluid flow through a single fracture. Water Resour. Res. 1984, 20, 1209–1215. [Google Scholar] [CrossRef]
  23. Brown, S. Fluid flow through rock Joints-the effect of surface roughness. J. Geophys. Res. 1987, 92, 1337–1347. [Google Scholar] [CrossRef]
  24. Xiong, F.; Jiang, Q.; Ye, Z.; Zhang, X. Nonlinear flow behavior through rough-walled rock fractures: The effect of contact area. Comput. Geotech. 2018, 102, 179–195. [Google Scholar] [CrossRef]
  25. Zou, L.; Jing, L.; Cvetkovic, V. Roughness decomposition and nonlinear fluid flow in a single rock fracture. Int. J. Rock Mech. Min. Sci. 2015, 75, 102–118. [Google Scholar] [CrossRef]
  26. Wang, M.; Chen, Y.; Ma, G.; Zhou, J.; Zhou, C. Influence of surface roughness on nonlinear flow behaviors in 3D self-affine rough fractures: Lattice Boltzmann simulations. Adv. Water Resour. 2016, 96, 373–388. [Google Scholar] [CrossRef]
  27. Yin, Q.; Ma, G.; Jing, H.; Wang, H.; Su, H.; Wang, Y.; Liu, R.C. Hydraulic properties of 3D rough-walled fractures during shearing: An experimental study. J. Hydrol. 2017, 55, 169–184. [Google Scholar] [CrossRef]
  28. He, Y.; Bai, B.; Hu, S.; Li, X. Effects of surface roughness on the heat transfer characteristics of water flow through a single granite fracture. Comput. Geotech. 2016, 80, 312–321. [Google Scholar] [CrossRef]
  29. Huang, Y.; Zhang, Y.; Yu, Z.; Ma, Y.; Zhang, C. Experimental investigation of seepage and heat transfer in rough fractures for enhanced geothermal systems. Renew. Energy 2019, 135, 846–855. [Google Scholar] [CrossRef]
  30. Ma, Y.; Zhang, Y.; Yu, Z.; Huang, Y.; Zhang, C. Heat transfer by water flowing through rough fractures and distribution of local heat transfer coefficient along the flow direction. Int. J. Heat Mass 2018, 119, 139–147. [Google Scholar] [CrossRef]
  31. Gringarten, A.C.; Witherspoon, P.A.; Ohnishi, Y. Theory of heat extraction from fractured hot dry rock. J. Geophys. Res. 1975, 80, 1120–1124. [Google Scholar] [CrossRef]
  32. Cheng, A.; Ghassemi, A.; Detournay, E. Integral equation solution of heat extraction from a fracture in hot dry rock. Int. J. Numer. Anal. Methods Geomech. 2001, 25, 1327–1338. [Google Scholar] [CrossRef]
  33. Martínez, A.R.; Roubinet, D.; Tartakovsky, D.M. Analytical models of heat conduction in fractured rocks. J. Geophys. Res. 2014, 119, 83–98. [Google Scholar] [CrossRef]
  34. Zhao, Z. On the heat transfer coefficient between rock fracture walls and flowing fluid. Comput. Geotech. 2014, 59, 105–111. [Google Scholar] [CrossRef]
  35. Yost, K.; Valentin, A.; Einstein, H. Estimating cost and time of wellbore drilling for Engineered Geothermal Systems (EGS)–Considering uncertainties. Geothermics 2015, 53, 85–99. [Google Scholar] [CrossRef]
  36. Brush, D.J.; Thomson, N.R. Fluid flow in synthetic rough-walled fractures: Navier-Stokes, Stokes, and local cubic law simulations. Water Resour. Res. 2003, 39, 1085. [Google Scholar] [CrossRef]
  37. Wang, Z.; Xu, C.; Dowd, P. Perturbation Solutions for Flow in a Slowly Varying Fracture and the Estimation of Its Transmissivity. Transp. Porous Media 2019, 128, 97–121. [Google Scholar] [CrossRef]
  38. Wang, Z.; Xu, C.; Dowd, P.; Xiong, F.; Wang, H. A non-linear version of the Reynolds equation for flow in rock fractures with complex void geometries. Water Resour. Res. 2020, 56, e2019WR026149. [Google Scholar] [CrossRef]
  39. Zimmerman, R.; Al-Yaarubi, A.; Pain, C.; Grattoni, C. Non-linear regimes of fluid flow in rock fractures. Int. J. Rock Mech. Min. Sci. 2004, 41, 163–169. [Google Scholar] [CrossRef]
  40. Chen, Y.F.; Zhou, J.Q.; Hu, S.H.; Hu, R.; Zhou, C.B. Evaluation of Forchheimer equation coefficients for non-Darcy flow in deformable rough-walled fractures. J. Hydrol. 2015, 529, 993–1006. [Google Scholar] [CrossRef]
  41. Barton, N.; Bandis, S.; Bakhtar, K. Strength, deformation and conductivity coupling of rock joints. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 1985, 22, 121–140. [Google Scholar] [CrossRef]
  42. Tatone, B.; Grasselli, G. Quantitative measurements of fracture aperture and directional roughness from rock cores. Rock Mech. Rock Eng. 2012, 45, 619–629. [Google Scholar] [CrossRef]
  43. Tse, R.; Cruden, D.M. Estimating joint roughness coefficients. Int. J. Min. Sci. Geomech. Abstr. 1979, 16, 303–307. [Google Scholar] [CrossRef]
  44. Shaik, A.R.; Rahman, S.S.; Tran, N.H.; Tran, T. Numerical simulation of fluid-rock coupling heat transfer in naturally fractured geothermal system. Appl. Therm. Eng. 2011, 31, 1600–1606. [Google Scholar] [CrossRef]
  45. Sisavath, S.; Al-Yaarubi, A.; Pain, C.C.; Zimmerman, R.W. A Simple Model for Deviations from the Cubic Law for a Fracture Undergoing Dilation or Closure. Pure Appl. Geophys. 2003, 160, 1009–1022. [Google Scholar] [CrossRef]
Figure 1. Schematic illustration of heat extraction from a geothermal reservoir.
Figure 1. Schematic illustration of heat extraction from a geothermal reservoir.
Energies 14 00951 g001
Figure 2. A conceptual diagram of heat flow through a rock fracture. (a) Symmetric fracture wedge; (b) Asymmetric fracture wedge.
Figure 2. A conceptual diagram of heat flow through a rock fracture. (a) Symmetric fracture wedge; (b) Asymmetric fracture wedge.
Energies 14 00951 g002
Figure 3. Computational flow chart of the proposed novel procedure for coupled simulation of thermal and fluid flow models (NPCTF). (a) Roughness-walled fracture; (b) Fluid flow model; (c) Heat transport model; (d) Flowchart of NPCT.
Figure 3. Computational flow chart of the proposed novel procedure for coupled simulation of thermal and fluid flow models (NPCTF). (a) Roughness-walled fracture; (b) Fluid flow model; (c) Heat transport model; (d) Flowchart of NPCT.
Energies 14 00951 g003
Figure 4. Temperature distribution along 5 fracture profiles calculated by the proposed NPCTF and COMSOL with different joint roughness coefficient (JRC) values. (a) JRC = 1; (b) JRC = 5; (c) JRC = 10; (d) JRC = 15; (e) JRC = 20.
Figure 4. Temperature distribution along 5 fracture profiles calculated by the proposed NPCTF and COMSOL with different joint roughness coefficient (JRC) values. (a) JRC = 1; (b) JRC = 5; (c) JRC = 10; (d) JRC = 15; (e) JRC = 20.
Energies 14 00951 g004aEnergies 14 00951 g004b
Figure 5. Temperature contour of fracture profile with different JRC (1–20).
Figure 5. Temperature contour of fracture profile with different JRC (1–20).
Energies 14 00951 g005
Figure 6. The transmissivity comparison calculated by the proposed NPCTF and COMSOL simulation in 5 fracture profiles.
Figure 6. The transmissivity comparison calculated by the proposed NPCTF and COMSOL simulation in 5 fracture profiles.
Energies 14 00951 g006
Figure 7. The 3D fracture model setup scanned by stereotopometric scanning system. (a) 3D fracture geometry model. (b) 3D fracture mesh model.
Figure 7. The 3D fracture model setup scanned by stereotopometric scanning system. (a) 3D fracture geometry model. (b) 3D fracture mesh model.
Energies 14 00951 g007
Figure 8. Temperature distribution of 3D fracture for (a) spatial temperature variation and (b) temperature change along fracture profiles.
Figure 8. Temperature distribution of 3D fracture for (a) spatial temperature variation and (b) temperature change along fracture profiles.
Energies 14 00951 g008
Figure 9. Pressure distribution of 3D fracture.
Figure 9. Pressure distribution of 3D fracture.
Energies 14 00951 g009
Figure 10. Deviation from heat flow simulations Ds for the proposed NPCTF, where the mean deviations of the model with different JRC are 2.57%, 1.85%, 1.29%, 1.10%, and 0.75%, respectively, as represented by dashed lines.
Figure 10. Deviation from heat flow simulations Ds for the proposed NPCTF, where the mean deviations of the model with different JRC are 2.57%, 1.85%, 1.29%, 1.10%, and 0.75%, respectively, as represented by dashed lines.
Energies 14 00951 g010
Figure 11. Mean deviation from heat flow simulations <D> for the proposed NPCTF, Gehlin and Hellstrom (2003) [2], and Zhao (2014) [35] models.
Figure 11. Mean deviation from heat flow simulations <D> for the proposed NPCTF, Gehlin and Hellstrom (2003) [2], and Zhao (2014) [35] models.
Energies 14 00951 g011
Figure 12. The effect of fracture aperture on mean deviation <D> from NPCTF for JRC = 16.34 fracture.
Figure 12. The effect of fracture aperture on mean deviation <D> from NPCTF for JRC = 16.34 fracture.
Energies 14 00951 g012
Figure 13. The effect of mesh size on mean deviation <D> from NPCTF or JRC = 16.34 fracture.
Figure 13. The effect of mesh size on mean deviation <D> from NPCTF or JRC = 16.34 fracture.
Energies 14 00951 g013
Table 1. The used parameters for heat transport modeling.
Table 1. The used parameters for heat transport modeling.
ParametersKr/W m−1 K−1ρw/kg m−3μ/pa scw/J kg−1 K−1Tin/KT0/K
Value3.510000.0014200293.15363.15
Table 2. Fracture geometry parameters.
Table 2. Fracture geometry parameters.
Fracture PositionJRCMean Aperture/mm
X = 1514.861.2
X = 4516.341.5
X = 7513.261.5
X = 10512.261.8
X = 13511.961.6
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Xiong, F.; Zhu, C.; Jiang, Q. A Novel Procedure for Coupled Simulation of Thermal and Fluid Flow Models for Rough-Walled Rock Fractures. Energies 2021, 14, 951. https://doi.org/10.3390/en14040951

AMA Style

Xiong F, Zhu C, Jiang Q. A Novel Procedure for Coupled Simulation of Thermal and Fluid Flow Models for Rough-Walled Rock Fractures. Energies. 2021; 14(4):951. https://doi.org/10.3390/en14040951

Chicago/Turabian Style

Xiong, Feng, Chu Zhu, and Qinghui Jiang. 2021. "A Novel Procedure for Coupled Simulation of Thermal and Fluid Flow Models for Rough-Walled Rock Fractures" Energies 14, no. 4: 951. https://doi.org/10.3390/en14040951

APA Style

Xiong, F., Zhu, C., & Jiang, Q. (2021). A Novel Procedure for Coupled Simulation of Thermal and Fluid Flow Models for Rough-Walled Rock Fractures. Energies, 14(4), 951. https://doi.org/10.3390/en14040951

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

Article Metrics

Back to TopTop