Abstract
The complex thermal-hydraulic-mechanical (THM) coupling is the key issue to the energy extraction from a geothermal reservoir, where fractures are the main channels for fluid circulation and heat transfer. However, the effects of matrix deformation-induced aperture variation and fracture roughness on heat recovery efficiency are unclear. In this paper, a fully coupling THM model based on a discrete fracture network is proposed to explore these coupling effects. First, the fracture roughness and the fracture aperture variation with effective stress are introduced. Second, the water flow and heat transfer in the matrix and fractures as well as the deformation of the geothermal reservoir are individually formulated for a fractured geothermal reservoir. Third, the model is validated with analytical solution for its thermal-hydraulic (TH) coupling effect and literature data for its hydraulic-mechanical (HM) coupling effect. Finally, the features of heat transfer and fluid flow in the fractured geothermal reservoir are comparatively analyzed through four scenarios. The simulation results indicate that the discrete fracture network severely impacts the pressure distribution and temperature advance. The aperture variation induced by solid deformation can enhance heat transfer efficiency, and the fracture roughness can reduce the heat transfer efficiency.
1. Introduction
A general operation method for energy extraction from geothermal reservoirs is to inject cold fluid in and to pump hot fluid out [1]. The fluid circulation and heat exchange are the key processes in energy extraction. The rock matrix in geothermal reservoirs is usually dense and nearly impervious; thus, natural or artificial fractures or joints are the main channels for fluid circulation and heat exchange [2–4]. This energy extraction process involves the interactions among thermal (T), hydro (H), and mechanical (M) processes [5, 6]. The injection and pump of working fluid will redistribute in situ stress, resulting in the variation of fracture aperture and directly impacting the transmissivity of the fractured geothermal reservoir. This may affect the energy extraction efficiency. Therefore, the THM coupling in a fractured geothermal reservoir is the key issue to the heat extraction performance of this reservoir.
Different models have been proposed for the THM coupling in an enhanced geothermal system (EGS). First, TH coupling models were developed for the forced convection in geothermal exploitations. They investigated the heterogeneity of porous media [7, 8], the system parameter optimization [9, 10], the heat exchange mechanism [11–13], and the effects of reservoir fracture treatments [14]. These TH coupling models can well evaluate the heat extraction in a rigid reservoir because they ignore the impact of reservoir deformation in the heat extraction process. Second, THM coupling models were developed based on continuum porous media [15–18]. Recent THM coupling models for geothermal reservoirs consider the impact of several isolated discrete fractures on heat extraction performance [19–21]. The fluid circulation in a complex fracture network is still unknown. The fluid circulation and heat transfer in a complex fracture network may heavily affect the extraction efficiency from a fractured geothermal reservoir, but such a study has rarely been observed in public literature.
A fractured geothermal reservoir can be modeled by either continuum-based models (CM) or discrete fracture models (DFM) [22, 23]. In the CM, the fractured reservoir is treated as a continuum composed of the fractures and matrix. Typical continuum-based models include an equivalent porous medium model [11, 12, 17] and dual-porosity model [24, 25]. These continuum-based models cannot explicitly express the fractures and thus cannot include the details of fracture characteristics such as location, length, and orientation. In the DFM, the geothermal reservoir is composed of the discrete fracture network (DFN) and rock matrix [6, 26, 27]. These DFM models consider the seepage properties of fractures in a different way from the rock matrix. The complex features of discrete fractures can be individually considered within geothermal reservoirs. However, the complexity of the fracture network largely increases the difficulty in both fracture description and computational cost. A real geothermal reservoir is fractured, and the properties of fractures may play a key role in energy extraction. Thus, the THM coupling analysis is necessary to include the discrete fracture network in the fractured geothermal reservoir.
The heat extraction performance of a geothermal system has been evaluated based on DFN models [13, 14, 28]. These models mainly focused on the heat exchange process but overlooked the effect of solid deformation and the fracture roughness. The impacts of aperture evolution and fracture roughness on heat extraction efficiency are still unknown. Hence, a DFN-based THM coupling model is necessary. This model can include the aperture evolution and the fracture roughness in the heat extraction from a geothermal system.
This paper develops a THM coupling model to consider the effects of aperture variation and fracture roughness for the heat energy extraction from a geothermal reservoir. Firstly, the discrete fracture network is generated through the Monte Carlo method. The fracture characteristics such as location, length, and orientation are included. Secondly, the aperture variation is described by a function of fracture aperture with effective stress. Further, a roughness factor is introduced into Darcy’s law to describe the effect of fracture roughness on fluid flow in a fracture. Thirdly, a full THM coupling model is established for the heat extraction process from a fractured geothermal reservoir and is validated by two degraded submodels. Finally, the distribution of pressure and temperature in the fractured reservoir is numerically simulated. Their coupling effects are comparatively analyzed.
2. Generation of the Fracture Network
2.1. Monte Carlo Method for Fracture Generation
A fracture network is constructed by randomly generating fractures through the algorithm of the Monte Carlo method. The trace length, orientation, and location of fractures are the three key parameters to a two-dimensional fracture network. Compared with trace length, fracture aperture is too tiny to be expressed on the graph of the fracture network although it severely impacts the fluid seepage. These three key parameters are determined through the following procedure.
2.2. Fracture Trace Lengths
Previous studies have shown that the trace length of fractures in geothermal reservoirs obeys a power-law distribution [29, 30]. This study expresses the trace length of each fracture by where and are the minimum and maximum trace lengths, is the fractal dimension of length distribution of fractures, and is a random number uniformly distributed in the range of 0 to 1.
2.3. Orientations of Fractures
The fractures in the geothermal reservoir can be grouped into different sets according to their orientations [29]. A fracture network may have one, two, or three sets of fractures. Each set of fractures has a mean orientation angle and a standard deviation . The orientation angle of fractures follows a Gaussian distribution. This fracture orientation angle after the Box-Muller transform can be generated by where is a fracture orientation angle in a set of fractures and and are random numbers uniformly distributed in the range of 0 to 1.
2.4. Location of Fractures
The midpoint of a fracture is used to indicate the location of a fracture. The coordinates of the fracture center is fixed by generating random numbers based on uniform distribution. The fracture centers are sequentially generated. If one fracture is too close to other fractures (<2 m in the domain of ), the fracture center will be regenerated. The coordinate generation function for the midpoints of fractures is expressed as where and are the length and width of the generation domain, respectively.
2.5. Generated Fracture Networks
Two types of fracture networks are generated in Figure 1. Table 1 lists the parameters used for the generation of these DFNs. Type (a) consists of a set of horizontal fractures and a set of vertical fractures. The fracture density of (a-1) is much larger than that of (a-2), and the connectivity of (a-1) is better than that of (a-2). Type (b) has two sets of oblique fractures, and each set has a fixed orientation. The fracture density of (b-2) is half of that of (b-1). All the four DFNs are imported into the COMSOL Multiphysics for numerical simulation. The simulation results of the four DFNs are discussed in Section 5 of this paper.
3. Coupling Model for Energy Extraction from a Geothermal Reservoir
3.1. Basic Assumptions
This fully coupling model has five separate yet interacting governing equations: a thermo-elastic mechanical equation, two flow equations (for the fractures and rock matrix), two heat transfer equations in the fractures and rock matrix. This model is mainly based on the following assumptions: (1) The rock matrix is continuous, isotropic, and homogeneous. (2) Surrounding rocks are impermeable. The fluid loss is not considered. (3) Both the matrix and fractures are saturated with single-phase liquid. (4) The fluid flow in both the matrix and fractures follows Darcy’s law. (5) The local thermal equilibrium is assumed at the interface between the rock matrix and the working fluid.
3.2. Governing Equations for Fluid Flow
3.2.1. Fluid Flow in the Matrix
In the saturated porous matrix, the mass conservation of fluid is where is the fluid density, is the porosity of the matrix, denotes time, is the flow velocity of the Darcy seepage, and is the source or sink term.
The seepage velocity is expressed by Darcy’s law as
Furthermore, the storage term at the first term of Equation (4) is
The porosity dilation-induced fluid variation is treated as the source term which is shown as where is the intrinsic permeability of the rock matrix, is the dynamic viscosity of water, is the storage coefficient of the porous matrix, is the water pressure, is Biot’s coefficient of the porous matrix, and is the volumetric strain.
Thus, the governing equation of water flow in the matrix is
3.2.2. Fluid Flow in Fractures
The continuity equation of fluid flow in discrete fractures is where is the fracture aperture, is the flow velocity in the fracture, denotes the gradient operator restricted to the tangential plane of a fracture, and is the source or sink term.
A fracture consists of two smooth parallel plates; thus, the volumetric flow rate is and the permeability of a single fracture is then obtained as where is the fracture width and is the fracture aperture. is the pressure gradient along the fracture.
We introduce of the roughness factor to consider the effect of fracture roughness on flow:
The following empirical expression is used in our computation: where is the standard deviation of the fracture surface height and is a constant ( by Lomize [31] and by Louis [32]). Single fracture flow found that varies from 1.04 to 1.65 [33]. Hence, the fracture flow velocity is
In this study, fractures are treated as the internal boundaries of the matrix seepage model. The fluid pressure in fractures is the Dirichlet boundary condition of fluid flow in the matrix, while the pore pressure of the matrix is the motivation of the source/sink term of fracture flow. The source or sink term in the governing equation for fracture flow is where and are the normal direction of up and down surfaces of the fracture, respectively.
Finally, the governing equation for fracture flow is obtained as
3.3. Governing Equations for Heat Transfer in the Matrix and Fractures
As a kind of porous media, the rock matrix has a higher specific surface area. Fluid can exchange heat with the matrix instantly. Hence, the temperature of pore water is assumed to be identical to that of the rock matrix. Then, the heat transfer equation in the matrix can be written as where is the effective heat capacity, is the porosity of the rock matrix, is the matrix density, is the specific heat capacity of the matrix, is the specific heat capacity of fluid, is the matrix temperature, and is the heat flux in the rock matrix.
The heat flux can be expressed through Fourier’s law as where is the effective heat conductivity, is the heat conductivity of the rock matrix, and is the heat conductivity of fluid.
The heat transfer equation in discrete fractures is where is the heat flux in the fracture and is
is the source or sink term in fracture heat transfer. It is the heat exchange between the matrix and fracture. Based on the local thermal equilibrium assumption, can be expressed as
3.4. Mechanical Equilibrium
The constitutive equation of poroelasticity for THM coupling is
The Navier equation for displacement can be expressed as where is the stress tensor (positive for tension), is the shear modulus, is the strain tensor, is Poisson’s ratio, and is the unit tensor. is Biot’s coefficient, which depends on the compressibility; is the pore pressure (negative for suction); is the bulk modulus of porous media, and is Young’s modulus. is the coefficient of thermal volumetric expansion of porous media. is the initial temperature of the rock mass. is the body force.
3.5. Multiphysical Couplings
3.5.1. Variation of Fracture Aperture with Effective Normal Stress
A geothermal reservoir is composed of a rock matrix and fracture network. The rock matrix is usually of low porosity and low permeability. The slight change of porosity has little effect on water flow and heat transfer; thus, the rock matrix is assumed to have constant porosity. The fracture network is the primary channel for fluid flow. Any change of fracture aperture may severely affect water flow. The fracture aperture is related to effective normal stress as [34] where is the fracture closure, where is the initial aperture with zero effective normal stress; is the normal component of the effective stress tensor ; is the maximum possible fracture closure; and is an empirical parameter representing the initial normal stiffness. and can be determined with the tested samples under the specific testing conditions. Baghbanan and Jing [35] found that the above two parameters are difficultly evaluated for a set of fractures. They presented some empirical expressions for use: .
3.5.2. Change of Fluid Property with Temperature
The fluid property includes water density and dynamic viscosity and affects water flow and heat transfer. It is found that pressure change does not significantly modify water density, but temperature severely modifies both water density and dynamic viscosity. A piecewise function is used for water density and dynamic viscosity as [36]
4. Simulation Implementation and Model Validation
4.1. Numerical Simulation Procedure
These governing equations with the initial and boundary conditions formulate a complete THM coupling model. Generally, these partial differential equations are too complex to derive analytical solutions. This study uses a finite element method (FEM) to solve this problem. Particularly, the FEM solutions are numerically obtained through the COMSOL Multiphysics, a powerful solver for partial differential equations. The solution procedure has four steps. Firstly, rock fractures are generated with the Monte Carlo algorithm. Then, the geometric model of fractured rock mass is imported into COMSOL Multiphysics. Thirdly, both governing equations and boundary conditions are assigned to the geometric model, and model parameters and initial conditions are set. Finally, the distributions of temperature and pressure are solved and visualized.
4.2. Model Parameters
A realistic geothermal reservoir usually has a three-dimensional (3D) fracture network. A complex 3D fracture network is extremely difficult in implementation and largely increases computation cost. In this study, a 2D fracture network in Figure 2 is used to model a geothermal reservoir. This geothermal reservoir is at 3500 m deep and has a computational domain of . The simulation time is 50 years from the beginning of cool water injection. Working fluid of water is injected into the geothermal reservoir. The left side boundary of the reservoir is the injection well, and the right side boundary is the production well. The initial and boundary conditions are described below.
For the reservoir seepage, the hydraulic pressure is kept as 80 MPa at the injection well and as 66 MPa at the right boundary. The top and bottom boundaries are impermeable. The initial pressure of this reservoir is 70 MPa. For the heat transfer in the reservoir, the water temperature is 30°C at the injection well. The top and bottom boundaries are thermally insulated. The initial temperature of the reservoir is 180°C. For the mechanical equilibrium model, the top and right boundaries are applied a constant normal compression of 90 MPa, which is equivalent to the overburden stress at 3500 m deep. The other two boundaries have zero normal displacement and free tangential displacement. Table 2 lists the model parameters used in simulations.
4.3. Model Validation
This fully coupling THM model is divided into two submodels: a TH coupling model for the heat transfer in a single fracture and an HM coupling model for the variation of fracture aperture. Figure 3 is the geometry model for the heat transfer in a single fracture. Cool water is injected into the fracture and heated by the impermeable rock matrix. For a problem similar to Figure 3 but with different boundary conditions, Ye and Wang [37] derived an analytical solution through the Laplace transform and Stehfest method. Following the same method, the temperature distribution of heat transfer in a single fracture can be obtained analytically. The deduction process is presented in the appendix, which is used for model validation.
The distribution of fracture temperature and the observation line at the 100th day are plotted in Figure 4(a), and the variation of temperature at the fixed points is plotted in Figure 4(b). The temperature distribution obtained from the current numerical simulation is in good agreement with the analytical solution. Only small deviation is observed at the observation point, but within the acceptable range.
(a) Temperature at the 100th day
(b) Temperature at fixed points
Figure 5 presents a geometry model with two horizontal wells and four fractures. The domain size and boundary conditions are marked. The two wells are located at the left and right boundaries. The four fractures are located inside the domain: fracture (a) has a length of 200 m and is connected with two wellbores, fracture (b) has a length of 100 m and is only connected with the injection wellbore, fracture (c) has a length of 100 m and is only connected with the production wellbore, and fracture (d) has a length of 100 m and is isolated from the two wellbores.
Under the same parameters, both the simulation results in this study and those from Moradi et al. [38] are plotted in Figure 6. Generally, both solutions are in good agreement. Our numerical model can correctly describe the variation of fracture aperture induced by the HM coupling effect.
(a)
(b)
(c)
(d)
5. Simulation Results and Discussion
5.1. Pore Pressure Distribution
The distribution of pore pressure has similar characteristics in the four fracture networks. Thus, only the pressure of the fracture network (a-1) is plotted for discussion. As shown in Figure 7, the pressure distribution is mainly along the fractures. That is because fractures are the main channels for fluid seepage. In the initial stage of water injection (before 3 years), the fractures beside the injection well have higher pressure than the matrix around the fracture. The fractures beside the production well have lower pressure than the matrix around the fracture. When the pressure of fracture water reaches a stable state, the pressure of matrix water is still changing continuously. After 15 years of water injection, the pore pressure of the reservoir tends to be stable. The area with connected fractures has a lower pressure gradient. That is, the water pressure decreases gently (the dash circle area of Figure 7). The area with nonconnected fractures has a higher pressure gradient, and the water pressure decreases sharply (the solid circle area of Figure 7).
5.2. Reservoir Temperature Distribution
The temperature distributions of the four fracture networks at the 30th year are presented in Figure 8. After cool water is injected into the geothermal reservoir, the water as working fluid is immediately heated by the rock matrix. At the same time, the temperature of the rock matrix decreases quickly. The flow velocity is much higher in fractures than in the matrix. Fracture temperature first decreases and then extends to the matrix. As more and more heat is extracted, the cooling area gradually expands from the injection well to the production well. However, the temperature at the production well does not decrease sharply before 30 years. This implies that the geothermal system can maintain a steady heat production for a quite long period.
If the fractures form a network which connects the injection well with the production well, the cooling area will advance faster along these fractures than in other zones (the dash circle area of Figure 8). Higher fracture density increases the probability of fracture intersection, resulting in better connectivity. Thus, the connectivity of (a-1) and (b-1) is better than that of (a-2) and (b-2). Actually, fracture networks (a-2) and (b-2) only have one cluster of fractures that connected both injection and production wells (dashed arrow line in Figure 8). Thus, their thermal front propagation is much slower than that of (a-1) and (b-1). The connected fractures are the main flow channels for fluid circulation. Their fracture flow flux is high, carrying away a majority of heat energy.
5.3. Temperature and Flow Velocity at the Outlet
The temperature and flow velocity at the production well are two important parameters to measure EGS performance. Based on the simulation results of the fracture network (a-1), the water temperature in the matrix and fractures along the production well is plotted in Figure 9(a). The line markers are for the matrix, and the scatter markers are for fractures. Since the permeability of fractures is much higher than that of the matrix, the flow in fractures is much faster than that in the matrix. The flow velocity in the matrix (lines) and fractures (scatter markers) along the production well is plotted in Figure 9(b).
(a) Temperature
(b) Flow velocity
The temperature at the production well is basically kept at about 180°C in the first 15 years of water injection. During this period, the EGS can maintain a high-efficiency heat production. After 30 years of cool water injection, the outlet temperature could remain above 140°C. Then, the cooling area gradually expands to the production well. The temperature drop along the production well is not uniform. It is faster in the zones with dense fractures where the faster flow of working fluid can more quickly transfer the thermal heat of the reservoir.
Under THM coupling, the aperture of fractures is in the range of 60 μm to 90 μm. Thus, the permeability in fractures is several orders of magnitude higher than that in the matrix. As shown in Figure 9(b), the outlet flow velocity in fractures (scatter markers) is about 7 orders of magnitude higher than that in the matrix (line markers). At the initial stage of heat extraction (before 1 year), the flow velocity in the matrix away from fractures is greater than that near fractures. The water in the matrix near fractures firstly flows into the fractures and continues to flow into the production well.
5.4. Comparison of Different Coupling Models
The TH coupling is the core for heat extraction. In the EGS operation, solid deformation and some other factors affect the TH coupling, thus modifying the heat extraction efficiency. In this section, four scenarios in the fully coupling THM model are investigated to comparatively analyze the different coupling effects on the outlet temperature and flow velocity. Table 3 lists the coupling effect in each scenario.
The initial stress equilibrium refers to the initial state of fluid-solid static coupling before cool water injection. Hence, the first step is to calculate the initial stress equilibrium. Then, the four models with different couplings are implemented to investigate their different coupling effects on heat transfer.
Based on the simulation results, the temperature distributions of the four scenarios at 30 years are presented in Figure 10. The temperature difference of the four scenarios mainly lies in the locations of thermal fronts. The thermal front of the THM model is more ahead than that of the other three models. This is mainly due to the increase in heat transfer efficiency caused by the enlargement of fracture aperture. Since fracture roughness results in the decrease in reservoir permeability, the thermal front of the TH+roughness model is in the most lag behind place.
The outlet temperature and flow velocity are plotted in Figure 11. Compared to the TH model, the THM model with smooth fracture has higher flow velocity but lower temperature at the production well. This implies that solid deformation can enhance reservoir permeability and finally increases heat transfer efficiency. If the fracture roughness is considered, the outlet velocity decreases but the outlet temperature increases. Hence, the fracture roughness leads to a decrease in reservoir permeability and finally decreases the heat transfer efficiency.
(a) Temperature 25th year
(b) Flow velocity 25th year
The average outlet temperature is calculated by
The summation term refers to the fractures while the integration term refers to the rock matrix. The average outlet temperature is a comprehensive indicator for the evaluation of EGS performance. The outlet temperatures at 25 years are plotted in Figure 12 for the four scenarios. It is observed that the TH model with fracture roughness has the highest average outlet temperature. The THM model with smooth fracture has the lowest average outlet temperature. The variation of fracture aperture induced by matrix deformation leads to the enhancement of heat transfer efficiency. However, fracture roughness leads to a decrease in heat transfer efficiency. It is observed that the average outlet temperature of the fully coupling model is similar to that of the simple TH model in this case.
6. Conclusions
This study developed a fully coupling thermal-hydro-mechanical model based on a discrete fracture network. The effects of deformation-induced aperture variation and fracture roughness on heat transfer efficiency in geothermal reservoirs were numerically investigated. Both TH and HM coupling effects in this coupling model were validated by either analytical solution or literature data. The different coupling effects on the distributions of pressure and temperature and heat recovery efficiency were comparatively analyzed through four scenarios. Based on these investigations, the following conclusions can be drawn. (1)The pressure distribution and temperature advance in the geothermal reservoir are severely impacted by the discrete fracture network. The lower pressure gradient and faster temperature advance usually occur in the zone with denser interconnected fractures(2)The deformation of the matrix induced by cooler water increases the aperture of the fracture and its permeability. This will enhance the heat transfer efficiency. However, the fracture roughness reduces the fracture flow velocity, thus reducing the heat transfer efficiency(3)The simple TH model can predict the similar outlet temperature to the fully coupling THM model in our calculation case. This may be due to the fact that the combining coupling effects of deformation and roughness are ignorable in that case. On this sense, a TH model is a simple and practical tool for heat extraction prediction
Appendix
A. Analytical Solution for Heat Transfer in a Smooth Fracture
A.1. TH Coupling Model for Heat Transfer in a Single Fracture
The conceptual geometric model is shown in Figure 3. All the properties of rock and water are kept as constants, and the heat conduction of water is omitted. For the convenience of symmetrical derivation, the fracture aperture is equivalent to .
Based on the cubic law of fracture flow, the flow velocity in the fracture is
Since all the parameters in Equation (A.1) are kept constant, the flow velocity is constant in this context. The heat transfer in the fracture can be described as
As the fracture rock geometry model is symmetrical on the midline and is too tiny compared to , Equation (A.1) can be expressed as
Further, the heat transfer equation in rock can be rewritten as where and are the temperatures of rock and water, respectively.
The initial conditions are
The boundary conditions are
The governing equations of Equations (A.3) and (A.4), combining the initial conditions of Equations (A.5) and (A.6) and boundary conditions of Equations (A.7), (A.8), and (A.9), form the simplified heat transfer model of the fractured rock. Its analytical solution can be obtained after variable dimensionless and Laplace transform.
A.2. Analytical Solution in the Laplace Domain
The following dimensionless variables are defined as
Thus, the governing equations can be expressed through these dimensionless variables as where
The dimensionless initial conditions are
The dimensionless boundary conditions are
The time-domain variables are transformed into the Laplace domain as
Combining the dimensionless initial conditions of Equations (A.13) and (A.14), the governing equations in the Laplace domain are
The boundary conditions in the Laplace domain can be expressed as
It is clear that Equation (A.18) is a second-order homogeneous linear differential equation. Its general solution is
Substituting Equation (A.22) into (A.21) yields
Substituting Equation (A.22) into (A.20) yields
Solving Equations (A.23) and (A.24) yields
With Equations (A.25) and (A.26), taking partial differentiation on both sides of Equation (A.22) gets where
Substituting Equation (A.27) into (A.17) yields
The general solution of Equation (A.29) is
Substituting Equation (A.30) into (A.19) yields
Hence, the analytical solutions in the Laplace domain are obtained as
Thus, the analytical solutions of temperature distribution are derived as Equations (A.32) and (A.33) in the Laplace domain. These analytical solutions should be inversely transformed back to the time domain to obtain their physical meaning. However, Equations (A.32) and (A.33) are too complex to perform Laplace inverse transform. In this study, a numerical inverse transform is performed to derive the solution in the time domain.
A.3. Stehfest Method of Numerical Inverse Transform
The Stehfest method is a numerical approximation for inverse transform to derive the original function in time domain. If the function has an image function in the Laplace domain, the Stehfest method expresses the transform of back to as where
Generally, is taken as an even number between 10 and 18.
Data Availability
Data is available on request.
Conflicts of Interest
None of the authors have any conflicts of interest.
Acknowledgments
The authors are grateful to the financial support from the National Natural Science Foundation of China (Grant No. 51674246), the Postgraduate Research & Practice Innovation Program of Jiangsu Province (Grant No. SJKY19_1859), and the Postgraduate Research & Practice Innovation Program of China University of Mining and Technology (Grant No. KYCX19_2170).