Geothermal-Related Thermo-Elastic Fracture Analysis by Numerical Manifold Method

One significant factor influencing geothermal energy exploitation is the variation of the mechanical properties of rock in high temperature environments. Since rock is typically a heterogeneous granular material, thermal fracturing frequently occurs in the rock when the ambient temperature changes, which can greatly influence the geothermal energy exploitation. A numerical method based on the numerical manifold method (NMM) is developed in this study to simulate the thermo-elastic fracturing of rocklike granular materials. The Voronoi tessellation is incorporated into the pre-processor of NMM to represent the grain structure. A contact-based heat transfer model is developed to reflect heat interaction among grains. Based on the model, the transient thermal conduction algorithm for granular materials is established. To simulate the cohesion effects among grains and the fracturing process between grains, a damage-based contact fracture model is developed to improve the contact algorithm of NMM. In the developed numerical method, the heat interaction among grains as well as the heat transfer inside each solid grain are both simulated. Additionally, as damage evolution and fracturing at grain interfaces are also considered, the developed numerical method is applicable to simulate the geothermal-related thermal fracturing process.


Introduction
Geothermal energy is one of the most abundant renewable energies on the Earth.Its environmental friendliness and reliable features compared with the conventional fossil fuels greatly encourage further exploration and exploitation of geothermal energy.One significant factor that influences the exploitation of geothermal energy is the variation of the mechanical properties of rock in high temperature environments [1][2][3].Therefore, investigating the thermo-mechanical mechanism of rock behavior under different thermal circumstances is attracting more and more attention nowadays.As a powerful and fundamental research approach, a great deal of laboratory tests have been conducted to investigate the thermo-mechanical mechanisms of rock.Based on laboratory tests, scholars have found that rocks are typical heterogeneous granular materials and the heterogeneity of rocks are mainly attributed to their random shaped mineral grains, cemented interfaces as well as microdefects (Figure 1a) [4][5][6].The heterogeneity of rock can dramatically influence the thermal and mechanical properties of rocks.Under thermal circumstances, the random shaped grain structure of rock can bring about high temperature gradients, uncoordinated thermal stresses and strains [7].Consequently, thermal fracturing of rocks occurs frequently (Figure 1b), which is an important factor that should be considered during geothermal energy exploitation.Thus, investigating the thermal fracturing process of rocklike granular material is of great importance for revealing the thermo-mechanical mechanism of the rock and improving the geothermal exploitation technologies.Indeed, considerable achievements have been attained based on laboratory tests.However, carrying out large numbers of laboratory tests is time-wasting and costly.Additionally, it is still difficult to completely monitor the micro failure process of materials in laboratory tests.Since rock are typical heterogeneous granular materials, developing numerical methods to investigate the thermal fracture effects of the rocklike granular materials is of great significance.that should be considered during geothermal energy exploitation.Thus, investigating the thermal fracturing process of rocklike granular material is of great importance for revealing the thermomechanical mechanism of the rock and improving the geothermal exploitation technologies.Indeed, considerable achievements have been attained based on laboratory tests.However, carrying out large numbers of laboratory tests is time-wasting and costly.Additionally, it is still difficult to completely monitor the micro failure process of materials in laboratory tests.Since rock are typical heterogeneous granular materials, developing numerical methods to investigate the thermal fracture effects of the rocklike granular materials is of great significance.

Basic Formulas
Basically, the thermal conduction of granular materials can be divided into two separate processes, i.e., the heat conduction of solid grains and the heat interactives among bonded grains [7].The authors consider, without loss of generality, the simple example shown in Figure 8, in which the schematic of the thermal conduction of a granular material consists of two separate grains, Ω and Ω , are presented.Figure 8 shows that except for the heat conduction in Ω and Ω , respectively, the heat transferred through the cemented interface Γ can also dramatically influence the temperature field of the granular material.Thus, apart from the prescribed temperature boundary condition on Γ and the prescribed heat flux boundary condition on Γ , the interface heat interactive condition at Γ should also be satisfied.Considering this, the basic formulas for thermal conduction of granular materials are given as follows:  The first law of thermodynamics To date, numerous kinds of numerical simulations have been conducted to investigate the thermal fracture effects of different kinds of materials.Via incorporating fracture mechanics, many numerical simulations focus on evaluating the thermal stress intensity factors of originally existing fractures, which is essential for predicting the thermal induced fracture propagation.During these numerical investigations, the continuum-based numerical methods, i.e., the finite element method (FEM) [8], the extended finite element method (XFEM) [9], the meshless methods (MMs) [10] and the boundary element methods (BEMs) [11], are generally adopted.However, due to the limited capacities of the continuum-based numerical methods in simulating complex discontinuous geometries as well as the fracture mechanics on treating complex fracture patterns, only simple macroscale thermal fracture problems are studied in these numerical simulations.
Since the thermal fracturing of rocklike materials are greatly affected by their microstructure, continuum-based microscale thermo-mechanical coupling numerical algorithms are also developed to investigate the thermal fracturing process of granular materials.One typical example is the RFPA-thermal software developed by the Tang's group [12], which is also an extension of the FEM.The RFPA-thermal represents the micro structure of rock with square finite elements.Through the incorporation of the statistically distributed material parameters and the thermo-elastic damage constitutive law, complex thermal fracturing from microscale to macroscale can be simulated by the RFPA-thermal method.However, since the displacements at the square grain boundaries are continuously simulated, the fracture paths are not explicitly represented in the RFPA-thermal.
Since the continuous numerical methods run into difficulties when simulating thermal fracturing of granular materials, the discontinuous numerical methods are also widely adopted to deal with such problems [13][14][15][16].One widely used discontinuous numerical method for thermal fracturing is the particle flow code (PFC), in which the micrograins are represented by circular or spherical elements [13].Incorporating the grain-based model and cluster model, which bond the same types of circular or spherical elements together, the PFC can realistically represent the grain structures of the rock [4].Thus, complex thermal fracturing processes that consider the micro structure of rock are simulated [3,14,17].As alternatives to the circle or sphere assembly-based numerical methods described above, polygon assembly-based simulations, such as the Voronoi polygon assembly-based simulations of Universal Distinct Element Code (UDEC) [18] and the triangular block assembly-based simulations of Discontinuous Deformation Analysis (DDA), are also extensively adopted in thermo-mechanical coupling problems [15].By incorporating a contact-based fracture criterion, complex thermal fracturing of rock can also be successfully simulated by the polygon assembly-based numerical methods.In the above discontinuous particle assembly numerical methods, one temperature degree of freedom is assigned to each grain to simulate the heat transfer process.The heat transfer of the material are then approximated by the heat transfer among blocks and the heat transfer inside each polygonal grain is generally ignored [15,18].Therefore, to obtain accurate temperature results, fine grain size, which means more temperature degrees of freedom, is generally necessary in these numerical methods.Additionally, the discontinuous numerical methods also run into difficulties in capturing the complex deformation field inside solid grains.
Since both the continuous and discontinuous numerical methods suffer limitations, combined continuous-discontinuous numerical methods have been developed during the past two decades.One of the most popular combined numerical methods is the numerical manifold method (NMM) [19].The NMM inherits the contact theory of the DDA and adopts dual cover systems, i.e., the mathematical cover (MC) and the physical cover (PC), thus providing unique advantages for modelling continuous and discontinuous problems on a unified numerical platform.Therefore, the NMM has been widely applied to simulate both continuous and discontinuous problems [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37].Nowadays, the NMM has also been extended to model thermo-mechanical fracture problems by Zhang [38,39].However, only thermal fracture stress intensity factors were investigated in Zhang's study.
Considering the advantages of the NMM in simulating both continuous and discontinuous problems, the NMM is extended to analyze the thermal fracturing of the rocklike granular materials.Since laboratory experiments have revealed that the microstructure of rocklike granular materials appear more like random polygons [40], in this study random Voronoi polygon assemblies are constructed to approximate the microstructure of the granular materials.Basically, two separate heat transfer processes should be considered, i.e., the heat transfer in solid grains and the heat interactions among grains.To more realistically simulate the heat interactions among grains, a contact-based heat transfer model is developed based on the contact algorithm of NMM.Then the transient heat transfer algorithm is discretized and incorporated into the NMM platform.Additionally, a damage-based contact fracture model is developed to reflect the bonding-cracking mechanism of grain boundaries.Since both the heat transfer in the solid grain and the heat interactions among grains are considered, the developed numerical method can better simulate the heterogeneous heat transfer process of rocklike granular materials.Additionally, as damage evolution and fracturing at grain interfaces is also considered, the developed numerical method is applicable to simulate the geothermal-related thermal fracturing process.Two simple examples are simulated to validate the developed NMM algorithm.Then, the extended NMM algorithm is applied to investigate the thermal fracturing processes induced by elevated temperature and cooling respectively.

Cover System
The most distinctive feature of the NMM is its dual coverage system, i.e., the mathematical cover and the physical cover, which have already been elaborated in many references [22].Attributed to the dual cover system, the NMM is capable of handling problems with complex physical meshes, such as cracks, holes, material interfaces, etc.The mathematical cover of NMM is a collection of small patches which can be of arbitrary shape and in theory it does not need to be consistent with the physical meshes.Thus, it is convenient to construct a numerical model with complex physical meshes in the NMM.Generally, the mathematical cover patches may partly overlap with each other and should cover the whole problem domain.In addition, the mathematical cover might intersect with the physical meshes.The subdivisions of the mathematical cover by the physical meshes are termed as the physical cover in NMM.The intersection of different physical covers is the manifold elements (ME) of NMM.
To clarify, a simple example is presented to illustrate the process of constructing the NMM cover system.Figure 2a shows a problem domain Ω is subdivided into three separate blocks, i.e., Ω 1 , Ω 2 and Ω 3 , by its physical meshes.To discretize Ω, two MCs, i.e., M 1 and M 2 in Figure 2b, are generated to completely cover Ω.Then, M 1 is subdivided into three PCs, termed as P 1  1 , P 2 1 and P 3 1 (Figure 3a), whereas M 2 is subdivided into three PCs, termed as P 1 2 , P 2 2 and P 3 2 (Figure 3b), by the physical meshes.The six PCs partly overlap each other and cut Ω into seven MEs, as shown in Figure 4.
Energies 2018, 11, x FOR PEER REVIEW 4 of 21 the physical meshes.The subdivisions of the mathematical cover by the physical meshes are termed as the physical cover in NMM.The intersection of different physical covers is the manifold elements (ME) of NMM.
To clarify, a simple example is presented to illustrate the process of constructing the NMM cover system.Figure 2a shows a problem domain Ω is subdivided into three separate blocks, i.e., Ω , Ω and Ω , by its physical meshes.To discretize Ω , two MCs, i.e., and in Figure 2b, are generated to completely cover Ω.Then, is subdivided into three PCs, termed as , and (Figure 3a), whereas is subdivided into three PCs, termed as , and (Figure 3b), by the physical meshes.The six PCs partly overlap each other and cut Ω into seven MEs, as shown in Figure 4. Following generating the cover system, a partition of unity function (weight function) ( ) and a polynomial local field function are constructed for each PC [41].Then, the global approximation of the field function in consideration is established by weight averaging the local field functions: the physical meshes.The subdivisions of the mathematical cover by the physical meshes are termed as the physical cover in NMM.The intersection of different physical covers is the manifold elements (ME) of NMM.
To clarify, a simple example is presented to illustrate the process of constructing the NMM cover system.Figure 2a shows a problem domain Ω is subdivided into three separate blocks, i.e., Ω , Ω and Ω , by its physical meshes.To discretize Ω , two MCs, i.e., and in Figure 2b, are generated to completely cover Ω.Then, is subdivided into three PCs, termed as , and (Figure 3a), whereas is subdivided into three PCs, termed as , and (Figure 3b), by the physical meshes.The six PCs partly overlap each other and cut Ω into seven MEs, as shown in Figure 4.   Following generating the cover system, a partition of unity function (weight function) ( ) and a polynomial local field function are constructed for each PC [41].Then, the global approximation of the field function in consideration is established by weight averaging the local field functions: Energies 2018, 11, x FOR PEER REVIEW 4 of 21 the physical meshes.The subdivisions of the mathematical cover by the physical meshes are termed as the physical cover in NMM.The intersection of different physical covers is the manifold elements (ME) of NMM.
To clarify, a simple example is presented to illustrate the process of constructing the NMM cover system.Figure 2a shows a problem domain Ω is subdivided into three separate blocks, i.e., Ω , Ω and Ω , by its physical meshes.To discretize Ω , two MCs, i.e., and in Figure 2b, are generated to completely cover Ω.Then, is subdivided into three PCs, termed as , and (Figure 3a), whereas is subdivided into three PCs, termed as , and (Figure 3b), by the physical meshes.The six PCs partly overlap each other and cut Ω into seven MEs, as shown in Figure 4.   Following generating the cover system, a partition of unity function (weight function) ( ) and a polynomial local field function are constructed for each PC [41].Then, the global approximation of the field function in consideration is established by weight averaging the local field functions: Following generating the cover system, a partition of unity function (weight function) ω i (x) and a polynomial local field function are constructed for each PC [41].Then, the global approximation of the field function in consideration is established by weight averaging the local field functions: where x is the position vector, N is the total number of the PC; P T represents the polynomial basis vector and a is the corresponding degree of freedom vector of the PC.In this study, the 0th-order polynomial local function is adopted to approximate the local field function.Thus, the global displacement field is written as: where d is the displacement degree of freedom vector and N u is the displacement shape function.
Similarly, the global temperature field can be written as: where θ is the temperature degree of freedom vector and N θ is the temperature shape function.

Contact Theory
The NMM inherits the contact theory from DDA to deal with block movements and block interactions [19].To avoid large displacements in a single time step, a maximum permissible displacement d 0 is prescribed in the NMM.Generally, the contact detection process of the NMM can be divided into two separate steps.The first step is detecting the possible contact pairs.The block pair is considered as a possible contact pair if the minimum distance between two blocks is no more than 2d 0 .Basically, three types of contact are defined in the NMM after possible contact detection, i.e., the angle-to-edge contact, the edge-to-edge contact and the angle-to-angle contact, as shown in Figure 5a-c where, x is the position vector, N is the total number of the PC; represents the polynomial basis vector and a is the corresponding degree of freedom vector of the PC.In this study, the 0th-order polynomial local function is adopted to approximate the local field function.Thus, the global displacement field is written as: ( ) ( ) where d is the displacement degree of freedom vector and is the displacement shape function.Similarly, the global temperature field can be written as: where θ is the temperature degree of freedom vector and is the temperature shape function.

Contact Theory
The NMM inherits the contact theory from DDA to deal with block movements and block interactions [19].To avoid large displacements in a single time step, a maximum permissible displacement is prescribed in the NMM.Generally, the contact detection process of the NMM can be divided into two separate steps.The first step is detecting the possible contact pairs.The block pair is considered as a possible contact pair if the minimum distance between two blocks is no more than 2 .Basically, three types of contact are defined in the NMM after possible contact detection, i.e., the angle-to-edge contact, the edge-to-edge contact and the angle-to-angle contact, as shown in Figures 5a-c, respectively.The second step is the invasion judgement of the possible contact pairs.Figure 6 shows when invasion exists in a possible contact pair, one or two point-edge pair are used to completely represent the invasion status of the possible contact pair.To be specific, for angle-to-edge contact and angle-toangle contact, one point-edge pair can be found, e.g., -in Figure 6a,c.However, for edge-toedge contact, two point-edge pairs should be specified to completely represent invasion status, e.g., the -and -in Figure 6b.After determining the invasion status of each possible contact pair, the normal invasion displacement Δ and the tangential invasion displacement Δ can be calculated according to the relative position of the corresponding point-edge pair.In NMM, linear elastic contact force−invasion displacement relationships are adopted.Thus the normal contact force and tangential contact force are calculated according to following equations [19]: The second step is the invasion judgement of the possible contact pairs.Figure 6 shows when invasion exists in a possible contact pair, one or two point-edge pair are used to completely represent the invasion status of the possible contact pair.To be specific, for angle-to-edge contact and angle-to-angle contact, one point-edge pair can be found, e.g., V 1 -V 2 V 3 in Figure 6a,c.However, for edge-to-edge contact, two point-edge pairs should be specified to completely represent invasion status, e.g., the V 1 -V 3 V 4 and V 2 -V 3 V 4 in Figure 6b.After determining the invasion status of each possible contact pair, the normal invasion displacement ∆o and the tangential invasion displacement ∆s can be calculated according to the relative position of the corresponding point-edge pair.In NMM, linear elastic contact force-invasion displacement relationships are adopted.Thus the normal contact force and tangential contact force are calculated according to following equations [19]: where F n and F s are the normal contact force and the tangential contact force, respectively.k n and k s are the corresponding normal contact stiffness and tangential contact stiffness, respectively.The selection of normal and tangential stiffness can greatly influence the output of the numerical model.According to Shi's suggestion, the ratio of the contact stiffness to the elastic modulus of the block is preferably between 20 and 100 [19].
Energies 2018, 11, x FOR PEER REVIEW 6 of 21 where and are the normal contact force and the tangential contact force, respectively.and are the corresponding normal contact stiffness and tangential contact stiffness, respectively.The selection of normal and tangential stiffness can greatly influence the output of the numerical model.According to Shi's suggestion, the ratio of the contact stiffness to the elastic modulus of the block is preferably between 20 and 100 [19].

Representation of the Micro Structure of Granular Materials
To more realistically model the thermo-elastic fracturing of the rocklike granular materials, the microstructure of these materials should be considered.The microstructure of the rocklike granular materials appear much like random polygons, as shown in Figure 7a [40].The Voronoi tessellation is introduced into the pre-processor of the NMM to represent the randomness of the microstructure of granular materials [42].The Voronoi tessellation begins with generating a set of random points.These random points can be moved step-by-step through an iteration procedure to make the space between the points more uniform.Then triangulation is carried out based on these random points.Following that, the Voronoi polygons are generated by constructing the perpendicular bisectors of all the triangle edges.Finally, by truncating the Voronoi polygons at the boundaries, the random microstructure of granular material is established (Figure 7b).The random microstructure is then treated as the physical meshes of the NMM pre-processor, and the NMM model for granular material finally is constructed by the pre-processor of the NMM (Figure 7c).

Basic Formulas
Basically, the thermal conduction of granular materials can be divided into two separate processes, i.e., the heat conduction of solid grains and the heat interactives among bonded grains [7].The authors consider, without loss of generality, the simple example shown in Figure 8, in which the schematic of the thermal conduction of a granular material consists of two separate grains, Ω and Ω , are presented.Figure 8 shows that except for the heat conduction in Ω and Ω , respectively, the

Representation of the Micro Structure of Granular Materials
To more realistically model the thermo-elastic fracturing of the rocklike granular materials, the microstructure of these materials should be considered.The microstructure of the rocklike granular materials appear much like random polygons, as shown in Figure 7a [40].The Voronoi tessellation is introduced into the pre-processor of the NMM to represent the randomness of the microstructure of granular materials [42].The Voronoi tessellation begins with generating a set of random points.These random points can be moved step-by-step through an iteration procedure to make the space between the points more uniform.Then triangulation is carried out based on these random points.Following that, the Voronoi polygons are generated by constructing the perpendicular bisectors of all the triangle edges.Finally, by truncating the Voronoi polygons at the boundaries, the random microstructure of granular material is established (Figure 7b).The random microstructure is then treated as the physical meshes of the NMM pre-processor, and the NMM model for granular material finally is constructed by the pre-processor of the NMM (Figure 7c).

Basic Formulas
Basically, the thermal conduction of granular materials can be divided into two separate processes, i.e., the heat conduction of solid grains and the heat interactives among bonded grains [7].The authors consider, without loss of generality, the simple example shown in Figure 8, in which the schematic of the thermal conduction of a granular material consists of two separate grains, Ω and Ω , are presented.Figure 8 shows that except for the heat conduction in Ω and Ω , respectively, the heat transferred through the cemented interface Γ can also dramatically influence the temperature field of the granular material.Thus, apart from the prescribed temperature boundary condition on Γ and the prescribed heat flux boundary condition on Γ , the interface heat interactive condition at

Basic Formulas
Basically, the thermal conduction of granular materials can be divided into two separate processes, i.e., the heat conduction of solid grains and the heat interactives among bonded grains [7].The authors consider, without loss of generality, the simple example shown in Figure 8, in which the schematic of the thermal conduction of a granular material consists of two separate grains, Ω 1 and Ω 2 , are presented.Figure 8 shows that except for the heat conduction in Ω 1 and Ω 2 , respectively, the heat transferred through the cemented interface Γ c can also dramatically influence the temperature field of the granular material.Thus, apart from the prescribed temperature boundary condition on Γ θ and the prescribed heat flux boundary condition on Γ q , the interface heat interactive condition at Γ c should also be satisfied.Considering this, the basic formulas for thermal conduction of granular materials are given as follows: • The first law of thermodynamics ∇q(x, t) where ρ is the mass density, q(x, t) is the heat flux, c θ is the specific heat, Q is the heat source per unit area, k is the thermal conductivity of the solid grains, ∇ is the gradient operator.
Correspondently, the boundary and initial conditions for thermal conduction are given as follows: • Initial conditions where n is the outside unit normal vector.Γ θ , Γ q are the boundaries with prescribed temperature θ and prescribed heat flux q, respectively.θ 0 represents initial temperature of the problem domain.Moreover, the solution of temperature field should also satisfy the interface heat interactive conditions: q(x, t) where n + and n − denote the corresponding outside unit normal vector of the upper boundary Γ + c and lower boundary Γ − c , respectively.q + c and q − c represent the amount of heat flux across the boundary Γ + c and Γ − c , respectively.
Energies 2018, 11, x FOR PEER REVIEW 7 of 21 heat transferred through the cemented interface Γ can also dramatically influence the temperature field of the granular material.Thus, apart from the prescribed temperature boundary condition on Γ and the prescribed heat flux boundary condition on Γ , the interface heat interactive condition at Γ should also be satisfied.Considering this, the basic formulas for thermal conduction of granular materials are given as follows:  The first law of thermodynamics where is the mass density, ( , ) is the heat flux, is the specific heat, Q is the heat source per unit area, k is the thermal conductivity of the solid grains, ∇ is the gradient operator.Correspondently, the boundary and initial conditions for thermal conduction are given as follows:  Initial conditions where n is the outside unit normal vector.Γ , Γ are the boundaries with prescribed temperature and prescribed heat flux , respectively.represents initial temperature of the problem domain.Moreover, the solution of temperature field should also satisfy the interface heat interactive conditions: ( ) where and denote the corresponding outside unit normal vector of the upper boundary Γ and lower boundary Γ , respectively.and represent the amount of heat flux across the boundary Γ and Γ , respectively.n Cemented Interface ( ) . Schematic of thermal conduction of granular materials.

Contact Based Heat Transfer Model
To obtain the temperature field of the granular material from the basic formulas presented above, the amount of the heat flux transferred through the cemented interfaces should be prescribed first.Inspired from the fact that the NMM adopts contact algorithm to deal with the interactions

Contact Based Heat Transfer Model
To obtain the temperature field of the granular material from the basic formulas presented above, the amount of the heat flux transferred through the cemented interfaces should be prescribed first.
Inspired from the fact that the NMM adopts contact algorithm to deal with the interactions between blocks, a contact based heat transfer model is developed to simulate the heat interactive among grains in this study.The heat interactive between the two grains are considered only when the two grains are cemented together or contact each other.Since the contact status of a possible contact can be completely represented by the geometric relation of corresponding point-edge relations, the contact-based heat transfer model in this study is also based on the point-edge relations.Additionally, according to Fourier's law, the amount of the heat flux transferred across the interface depends on the heat conduction capacity of the interface.Thus, in the rest of this article, a contact thermal conductivity k c is defined to represent the heat conduction capacity of the interface.The physical meaning of k c is the amount of heat flux transferred across unit length of interface per unit time in case of unit temperature difference.The heat conduction capacity of the cemented interfaces might be related to many factors, such as the contact force, the roughness of the interface and the hardness of the material [43,44].However, this is not the main research objective of this work.
To explain the contact-based heat transfer model in detail, a representative example is shown in Figure 9, in which block B (1) contacts block B (2) at point V 1 of B (1) and point V 0 of B (2) respectively.The contact status of this contact pair is represented with the point-edge relation V 1 -V 4 V 5 .The heat interactive between B (1) and B (2) is then approximated by the heat interactive between boundary V 4 V 5 and point V 1 .
Energies 2018, 11, x FOR PEER REVIEW 8 of 21 between blocks, a contact based heat transfer model is developed to simulate the heat interactive among grains in this study.The heat interactive between the two grains are considered only when the two grains are cemented together or contact each other.Since the contact status of a possible contact can be completely represented by the geometric relation of corresponding point-edge relations, the contact-based heat transfer model in this study is also based on the point-edge relations.Additionally, according to Fourier's law, the amount of the heat flux transferred across the interface depends on the heat conduction capacity of the interface.Thus, in the rest of this article, a contact thermal conductivity is defined to represent the heat conduction capacity of the interface.The physical meaning of is the amount of heat flux transferred across unit length of interface per unit time in case of unit temperature difference.The heat conduction capacity of the cemented interfaces might be related to many factors, such as the contact force, the roughness of the interface and the hardness of the material [43,44].However, this is not the main research objective of this work.
To explain the contact-based heat transfer model in detail, a representative example is shown in Figure 9, in which block ( ) contacts block ( ) at point of ( ) and point of ( ) respectively.The contact status of this contact pair is represented with the point-edge relation -.The heat interactive between ( ) and ( ) is then approximated by the heat interactive between boundary and point .
Contact area Contact based heat transfer model.
As the thickness and the volume of the interfaces are not explicitly represented in the NMM.The cemented interfaces are assumed to be with zero thickness and zero volume.Based on this kind of treatment, it is supposed that the cemented interface between blocks have no heat storage capacity and ignoring the frictional heat.Thus, the following balance equation of the heat flux at the interface is thus obtained: According to Fourier's law, the amount of heat flux transferred through the interface is related to the contact thermal conductivity, the interface temperature difference as well as the contact area.Thus, the amount of interface heat flux can be written as: with l representing the contact area as shown in Figure 9: ( ) where , and are the position vectors of , and respectively.∆ in Equation ( 14) represents the temperature difference between and , which is approximated by the temperature difference between and the projection point as shown in Figure 9: where and represent the temperatures at point and , respectively, which can be approximated by the linear interpolation of the corresponding PC temperatures: As the thickness and the volume of the interfaces are not explicitly represented in the NMM.The cemented interfaces are assumed to be with zero thickness and zero volume.Based on this kind of treatment, it is supposed that the cemented interface between blocks have no heat storage capacity and ignoring the frictional heat.Thus, the following balance equation of the heat flux at the interface is thus obtained: According to Fourier's law, the amount of heat flux transferred through the interface is related to the contact thermal conductivity, the interface temperature difference as well as the contact area.Thus, the amount of interface heat flux can be written as: with l representing the contact area as shown in Figure 9: where x 1 , x 2 and x 3 are the position vectors of V 1 , V 2 and V 3 respectively.∆θ c in Equation ( 14) represents the temperature difference between V 4 V 5 and V 1 , which is approximated by the temperature difference between V 1 and the projection point V 0 as shown in Figure 9: Energies 2018, 11, 1380 9 of 21 where θ + c and θ − c represent the temperatures at point V 1 and V 0 , respectively, which can be approximated by the linear interpolation of the corresponding PC temperatures: Thus, Equation ( 16) can be shortened as:

Discretization and Solution
To obtain the temperature field of the granular material, the governing equation should be established by discretizing the basic formulas described in Section 2.2.2 and the contact-based heat transfer model in Section 2.2.3.In this section, the weight residual method is adopted to establish the governing equation [41].The weak form of the governing equation can be expressed as: where δθ h is the first order variation of θ h (x, t) and λ T is the penalty number to enforce the prescribed temperature boundary condition.The last term in Equation ( 19) is the contact heat transfer terms which have to satisfy the heat flux continuity condition given in Equation ( 13) and the temperature condition given in Equation (16).By substituting Equation ( 14) and (18) into Equation ( 19), the discretized governing equation for transient heat conduction with contact heat transfer among cemented grains can be established: Equation ( 20) is a first-order time dependent differential equation, which can be solved by the time integration scheme [45].Generally, the time integration scheme for transient heat conduction can be uniformly written as: in which, k is the serial number of the time step, ∆t k+1 = t k+1 − t k is the time interval of the k + 1 time step, θ k and θ k+1 are the temperature degree of freedom vectors of the k time step and the k + 1 time step respectively.The parameter β is adopted to choose the specific time integration scheme.Regarding 1/2 ≤ β ≤ 1, the time integration scheme is unconditionally stable [45].In this study, β = 1 is adopted to solve the governing equation as it is unconditionally stable and relatively easy to be implemented into NMM.Thus, Equation ( 32) can be simplified as: 2.3.Thermo-Elastic Fracture Analysis

Thermo-Mechanical Coupling
The variation of temperature can lead to thermal expansion or thermal shrinkage of materials.As materials are generally under different external constrains, the thermal expansion or shrinkage of the materials cannot occur freely in all directions.Subsequently, thermal stresses and even thermal fractures usually appears in the materials.The thermal induced stress can be written as: where σ θ and ε θ are the thermal induced stress and strain vector respectively, ∆θ h is the amount of temperature variation, α represents the thermal expansion coefficient, D is the linear elastic constitutive matrix and I is the second order identity tensor.The thermal induced stress can be easily incorporated into the NMM by treating σ θ as the initial stress of the NMM [38].

A Damage Based Contact Fracture Model
The micrograins of granular materials in Nature are bonded together by the cemented interface.However, the cohesion effects among blocks are neglected in the original contact algorithm of the NMM.Additionally, under thermal circumstances, the thermal induced stresses could result in fracture evolution in the granular materials.Thus, the contact algorithm of the NMM should be improved to successfully simulate the cohesion effects and the thermal fracture evolution of the granular materials.To achieve this goal, a damage-based contact fracture model is incorporated into the NMM in this study [46].
As mentioned in Section 2.1.2,after the invasion judgement, the invasion status of contact block pairs are obtained.Generally, the interactions among grains are divided into normal behavior and tangential behavior.In NMM, normal contact springs and tangential contact springs are adopted to simulate the normal behaviors and tangential behaviors among blocks as shown in Figure 10.To be consistent with this and to simulate the cohesion effects among blocks as well, the normal force-aperture relationship shown in Figure 11a and the tangential force-aperture relationship shown in Figure 11b are adopted and assigned to the normal spring and tangential spring respectively.Correspondently, to reflect the damage evolution process of the cemented interface, two damage factors, i.e., the normal damage factor D n and the tangential damage factor D s , are defined for the normal contact spring and tangential contact spring respectively.The normal behavior and tangential behavior of the damage based contact fracture model are described in detail as follows.
As shown in Figure 11a, the normal behavior of the contact is divided into three parts according to the normal damage factor.The normal damage factor can be calculated by the following formula: where o is normal aperture, o u and o f represent the normal aperture corresponding to the initial tensile strength of the contact F t0 and the tensile failure of the contact, respectively: where k n0 is the initial stiffness of the normal spring and G I represents the mode I critical energy release rate.The tensile strength of the contact F t0 can be determined by F t0 = l f t0 , in which l is the contact area defined in Figure 9 and f t0 is the initial tensile strength of the interface per unit length.
granular materials.To achieve this goal, a damage-based contact fracture model is incorporated into the NMM in this study [46].As mentioned in Section 2.1.2,after the invasion judgement, the invasion status of contact block pairs are obtained.Generally, the interactions among grains are divided into normal behavior and tangential behavior.In NMM, normal contact springs and tangential contact springs are adopted to simulate the normal behaviors and tangential behaviors among blocks as shown in Figure 10.To be consistent with this and to simulate the cohesion effects among blocks as well, the normal forceaperture relationship shown in Figure 11a and the tangential force-aperture relationship shown in Figure 11b are adopted and assigned to the normal spring and tangential spring respectively.Correspondently, to reflect the damage evolution process of the cemented interface, two damage factors, i.e., the normal damage factor and the tangential damage factor , are defined for the normal contact spring and tangential contact spring respectively.The normal behavior and tangential behavior of the damage based contact fracture model are described in detail as follows.As shown in Figure 11a, the normal behavior of the contact is divided into three parts according to the normal damage factor.The normal damage factor can be calculated by the following formula: where o is normal aperture, and represent the normal aperture corresponding to the initial tensile strength of the contact and the tensile failure of the contact, respectively: where is the initial stiffness of the normal spring and represents the mode I critical energy release rate.The tensile strength of the contact can be determined by 0 0 t t F lf = , in which l is the contact area defined in Figure 9 and is the initial tensile strength of the interface per unit length.Based on the above descriptions, the normal force of the damage based contact fracture model can be expressed as follows: ( ) Based on the above descriptions, the normal force of the damage based contact fracture model can be expressed as follows: Similar to the normal behavior, the tangential behavior of the contact is also divided into three parts according to the tangential damage factor (Figure 11b).The tangential damage factor can also be calculated according to the tangential aperture: where s represents the tangential aperture of the contact.s u and s f are the tangential apertures corresponding to the initial shear strength of the contact F s0 and the shear failure of the contact respectively, which can be calculated as follows: where G I I denotes the mode II critical energy release rate; c represents the initial cohesion of the contact; k s0 is the initial stiffness of the tangential spring.According to Coulomb's friction law, the initial shear strength F s0 of the contact can be written as: where ϕ is the internal friction angle of the contact.As the cohesion of the contact decreases linearly with the increase of tangential damage factor as shown in Figure 11b, the tangential contact force-aperture relationship can thus be expressed as: Furthermore, to distinguish the mode I fracture, the mode II fracture and the mixed mode fracture as well as judge the failure of the cemented interface, a mixed mode damage factor is defined.The mixed mode damage factor is calculated by the following equation: It is supposed that whenever the mixed damage factor exceeds 1, the contact is set as "failure".

Validation of Contact Heat Transfer Model
The thermal conduction in a plane consists of two square blocks is simulated in this example to verify the effectiveness of the developed contact heat transfer model.The edge length of the square blocks is L = 1.5 mm.The plane is subjected to prescribed temperature conditions on the top boundary and bottom boundary, as shown in Figure 12a.The two blocks are bonded together by the center cemented interface.The heat transfer across the cemented interface should influence the temperatures on the two sides of the cemented interface, which is termed as θ 1c and θ 2c in this example.The steady state theoretical solutions for θ 1c and θ 2c can be calculated according to the following equations [43]: where η indicates the ratio of the contact thermal conductivity to the material thermal conductivity: The NMM cover system shown in Figure 12b is constructed to implement the numerical simulation.The numerical results of θ 1c and θ 2c for cases with different η are compared with the corresponding theoretical solutions in Figure 13a, which indicates good agreements.Moreover, it is illustrated in Figure 13 that the temperature difference at the cemented interface decreases with the increase of η.Therefore, it is inferred that when η is large enough, the interface temperature difference would diminish to 0. To validate this point of view, Figure 13b shows the temperature on the center vertical line of the plane for case with η = 1 × 10 6 .Obviously, the temperature agrees well with the exact temperature solution of a continuous plane with the same size and boundary conditions.would diminish to 0. To validate this point of view, Figure 13b shows the temperature on the center vertical line of the plane for case with η = 1 × 10 6 .Obviously, the temperature agrees well with the exact temperature solution of a continuous plane with the same size and boundary conditions.
Cemented interface

Validation of Thermo-Mechanical Coupling
A simple thermo-mechanical coupling problem is analyzed to validate effectiveness of the extended NMM on simulating thermo-mechanical coupling problems.As shown in Figure 14a,   would diminish to 0. To validate this point of view, Figure 13b shows the temperature on the center vertical line of the plane for case with η = 1 × 10 6 .Obviously, the temperature agrees well with the exact temperature solution of a continuous plane with the same size and boundary conditions.
Cemented interface

Validation of Thermo-Mechanical Coupling
A simple thermo-mechanical coupling problem is analyzed to validate effectiveness of the extended NMM on simulating thermo-mechanical coupling problems.As shown in Figure 14a,

Validation of Thermo-Mechanical Coupling
A simple thermo-mechanical coupling problem is analyzed to validate effectiveness of the extended NMM on simulating thermo-mechanical coupling problems.As shown in Figure 14a, a square plane with edge length L = 1 m and initial temperature θ 0 = 0 • C is subjected to thermal loading θ 1 = 1 • C on the top boundary.The vertical displacement of the bottom boundary as well as the horizontal displacements of the right boundary and left boundary are fixed.The plane is subjected to plane strain conditions.According to thermodynamics, the heat flux would flow gradually from the top boundary to the bottom boundary.The resulting temperature variation could cause thermal expansion and thermal stresses in the plane.To investigate these phenomenon, this transient heat transfer process is numerically simulated by the extended NMM.The variation of the temperatures and lateral stresses at point A (0.5, 0) and point B (0.5, 0.5) with time for cases with different time step sizes are extracted and investigated as shown in Figure 15a,b.To compare, the corresponding theoretical solutions for the temperatures and stresses are obtained according to the following equations [47]: It is clearly indicated in the figure that good agreements between the numerical results and theoretical solutions are obtained, which confirms the effectiveness of the extended NMM on simulating thermo-mechanical coupling problems.Moreover, the figures illustrate that the numerical results are more accurate for smaller time step size, e.g., ∆ = 0.002 s in Figure 15.According to thermodynamics, the heat flux would flow gradually from the top boundary to the bottom boundary.The resulting temperature variation could cause thermal expansion and thermal stresses in the plane.To investigate these phenomenon, this transient heat transfer process is numerically simulated by the extended NMM.The variation of the temperatures and lateral stresses at point A (0.5, 0) and point B (0.5, 0.5) with time for cases with different time step sizes are extracted and investigated as shown in Figure 15a,b.To compare, the corresponding theoretical solutions for the temperatures and stresses are obtained according to the following equations [47]: It is clearly indicated in the figure that good agreements between the numerical results and theoretical solutions are obtained, which confirms the effectiveness of the extended NMM on simulating thermo-mechanical coupling problems.Moreover, the figures illustrate that the numerical results are more accurate for smaller time step size, e.g., ∆t = 0.002 s in Figure 15.According to thermodynamics, the heat flux would flow gradually from the top boundary to the bottom boundary.The resulting temperature variation could cause thermal expansion and thermal stresses in the plane.To investigate these phenomenon, this transient heat transfer process is numerically simulated by the extended NMM.The variation of the temperatures and lateral stresses at point A (0.5, 0) and point B (0.5, 0.5) with time for cases with different time step sizes are extracted and investigated as shown in Figure 15a,b.To compare, the corresponding theoretical solutions for the temperatures and stresses are obtained according to the following equations [47]: It is clearly indicated in the figure that good agreements between the numerical results and theoretical solutions are obtained, which confirms the effectiveness of the extended NMM on simulating thermo-mechanical coupling problems.Moreover, the figures illustrate that the numerical results are more accurate for smaller time step size, e.g., ∆ = 0.002 s in Figure 15.

Thermo-Elastic Fracturing Induced by Elevated Temperature
The thermo-elastic fracturing of a model consists of an outer cylinder and an inner disc is numerically investigated in this example.The radii of the cylinder and the disc are shown in Figure 16a.The initial temperature of the model is 0 • C. The cylinder and disc are made of two different materials respectively, termed as Material 1 and Material 2 as shown in Figure 16a, where Material 1 is supposed to be a typical granular material.The temperature of the model is supposed to be increased uniformly, in which condition only the thermal fracturing process is considered and the heat transfer across the model can be ignored.The material parameters of Material 1 and Material 2 are presented in Table 1.Since the elastic modulus and the thermal expansion coefficient of Material 2 are larger than those of Material 1, thermal fracturing would occur in Material 1 when the temperature of the model is uniformly elevated.To investigate this thermal fracturing process, the cylinder is subdivided into 1241 Voronoi polygons as shown in Figure 16b.The contact-based heat transfer model and the damage-based contact fracture model are used to simulate the heat interactive and the bonding-cracking effects between the micro grains.The parameters of the interface are presented in Table 1 as well.

Thermo-Elastic Fracturing Induced by Elevated Temperature
The thermo-elastic fracturing of a model consists of an outer cylinder and an inner disc is numerically investigated in this example.The radii of the cylinder and the disc are shown in Figure 16a.The initial temperature of the model is 0 °C.The cylinder and disc are made of two different materials respectively, termed as Material 1 and Material 2 as shown in Figure 16a, where Material 1 is supposed to be a typical granular material.The temperature of the model is supposed to be increased uniformly, in which condition only the thermal fracturing process is considered and the heat transfer across the model can be ignored.The material parameters of Material 1 and Material 2 are presented in Table 1.Since the elastic modulus and the thermal expansion coefficient of Material 2 are larger than those of Material 1, thermal fracturing would occur in Material 1 when the temperature of the model is uniformly elevated.To investigate this thermal fracturing process, the cylinder is subdivided into 1241 Voronoi polygons as shown in Figure 16b.The contact-based heat transfer model and the damage-based contact fracture model are used to simulate the heat interactive and the bonding-cracking effects between the micro grains.The parameters of the interface are presented in Table 1 as well.The NMM cover system with 11,187 MEs and 10,955 PCs are generated to simulate this problem (Figure 16b).According to the analysis of Example 2, the time step size adopted for the numerical simulation is ∆ = 0.001 s .During the numerical simulation, the temperature of the model is increased uniformly by 0.3 °C per time step.The numerical results of the thermal fracturing  The NMM cover system with 11,187 MEs and 10,955 PCs are generated to simulate this problem (Figure 16b).According to the analysis of Example 2, the time step size adopted for the numerical simulation is ∆t = 0.001 s.During the numerical simulation, the temperature of the model is increased uniformly by 0.3 • C per time step.The numerical results of the thermal fracturing processes are shown in Figure 17.A fracture initiates at the lower right corner of the cylinder first (Figure 17a).This fracture propagates rapidly until it completely cut-through the cylinder at 480 • C (Figure 17b).Subsequently, two secondary fractures initiate at the upper left and lower left corners of the cylinder at 651 • C (Figure 17c) and 750 • C (Figure 17d), respectively.Similar numerical simulation of this example has been conducted by the RFPA-thermal [12].Figure 18a shows the final thermal fracture morphology of the RFPA-thermal, in which one main radial fracture that cut-through the outer cylinder as well as some shorter secondary radial fractures are also obtained.A similar experiment test has also conducted by Abdalla [49].In Abdalla's experiment test, a concrete prism symmetrically reinforced by FRP bars is subjected to uniform rise of temperature.The thermal expansion coefficient of FRP in the test is much larger than that of the concrete.As a result, transverse thermal stresses gradually increase in the concrete, which further leads to the radial fracturing in the concrete surrounding the FRP bars. Figure 18b shows the thermal fracture morphology of the test.Similar to the numerical results, one main radial fracture is initiated and finally cut-through the surrounding concrete.Along with the main fracture, there are also some shorter radial fractures in the surrounding concrete.processes are shown in Figure 17.A fracture initiates at the lower right corner of the cylinder first (Figure 17a).This fracture propagates rapidly until it completely cut-through the cylinder at 480 °C (Figure 17b).Subsequently, two secondary fractures initiate at the upper left and lower left corners of the cylinder at 651 °C (Figure 17c) and 750 °C (Figure 17d), respectively.Similar numerical simulation of this example has been conducted by the RFPA-thermal [12].Figure 18a shows the final thermal fracture morphology of the RFPA-thermal, in which one main radial fracture that cutthrough the outer cylinder as well as some shorter secondary radial fractures are also obtained.A similar experiment test has also conducted by Abdalla [49].In Abdalla's experiment test, a concrete prism symmetrically reinforced by FRP bars is subjected to uniform rise of temperature.The thermal expansion coefficient of FRP in the test is much larger than that of the concrete.As a result, transverse thermal stresses gradually increase in the concrete, which further leads to the radial fracturing in the concrete surrounding the FRP bars. Figure 18b shows the thermal fracture morphology of the test.Similar to the numerical results, one main radial fracture is initiated and finally cut-through the surrounding concrete.Along with the main fracture, there are also some shorter radial fractures in the surrounding concrete.processes are shown in Figure 17.A fracture initiates at the lower right corner of the cylinder first (Figure 17a).This fracture propagates rapidly until it completely cut-through the cylinder at 480 °C (Figure 17b).Subsequently, two secondary fractures initiate at the upper left and lower left corners of the cylinder at 651 °C (Figure 17c) and 750 °C (Figure 17d), respectively.Similar numerical simulation of this example has been conducted by the RFPA-thermal [12].Figure 18a shows the final thermal fracture morphology of the RFPA-thermal, in which one main radial fracture that cutthrough the outer cylinder as well as some shorter secondary radial fractures are also obtained.A similar experiment test has also conducted by Abdalla [49].In Abdalla's experiment test, a concrete prism symmetrically reinforced by FRP bars is subjected to uniform rise of temperature.The thermal expansion coefficient of FRP in the test is much larger than that of the concrete.As a result, transverse thermal stresses gradually increase in the concrete, which further leads to the radial fracturing in the concrete surrounding the FRP bars. Figure 18b shows the thermal fracture morphology of the test.Similar to the numerical results, one main radial fracture is initiated and finally cut-through the surrounding concrete.Along with the main fracture, there are also some shorter radial fractures in the surrounding concrete.(Figure 20a).These small fractures propagate downwards as the plane cools gradually from the top boundary to the bottom boundary.Additionally, the number of fractures that continue to grow downwards decrease monotonically.Generally, the ultimate thermal fracture patterns are mainly in the vertical direction, which is parallel to the average temperature gradient.However, due to the heterogeneity of the grain structure, the thermal fracture patterns are slightly deviated from the vertical direction.This numerical result is consistent with that obtained by Huang [50].
Energies 2018, 11, x FOR PEER REVIEW 18 of 21 mainly in the vertical direction, which is parallel to the average temperature gradient.However, due to the heterogeneity of the grain structure, the thermal fracture patterns are slightly deviated from the vertical direction.This numerical result is consistent with that obtained by Huang [50].

Conclusions
Because of the environmental friendliness and reliable features of geothermal energy, exploration and exploitation of geothermal energy is attracting more and more attention.However, since rocks are typical heterogeneous granular materials, the randomly shaped microstructure of the rock can bring about uncoordinated thermal stresses and even result in thermal fracturing when the deep high temperature circumstance is disturbed during geothermal energy exploitation.The geothermal-related thermal fracturing of rock can significantly influence the geothermal energy exploitation.
In this study, a thermo-mechanical coupling algorithm based on the numerical manifold method was developed to analyze the thermo-elastic fracturing process of granular materials.The developed NMM algorithm represented the randomly shaped microstructure of granular materials using Voronoi polygon assemblies.Based on the Voronoi polygon assembly, a transient thermal conduction algorithm for granular materials was established, in which the heat interactions among Voronoi grains were reflected by a newly developed contact-based heat transfer model.Then, a damage-based contact fracture model was incorporated in the NMM algorithm to simulate the bonding-cracking effects between grains.The newly developed NMM algorithm was validated by two simple examples, the results of which agreed well with that of theoretical solutions.At last, the thermoelastic fracturing processes induced by elevated temperature and cooling were simulated by the developed NMM algorithm.The obtained thermal fracture morphologies were consistent with those obtained by experiment test and available numerical simulations.
The developed numerical method is a coupled continuous-discontinuous numerical method.Since both the heat transfer in the solid grain and the heat interactions among grains are considered, the developed numerical method can better simulate the heterogeneous heat transfer process of rocklike granular materials.Additionally, as damage evolution and fracturing at grain interfaces is also considered, the developed numerical method is applicable to simulate the geothermal-related thermal fracturing process.However, in the numerical method developed, the effects of seepage, which is another significant influence of the geothermal energy exploitation, are not considered.In the future, the coupled thermal-hydraulic-mechanical numerical method will be developed and applied to analyze geothermal related problems.

Conclusions
Because of the environmental friendliness and reliable features of geothermal energy, exploration and exploitation of geothermal energy is attracting more and more attention.However, since rocks are typical heterogeneous granular materials, the randomly shaped microstructure of the rock can bring about uncoordinated thermal stresses and even result in thermal fracturing when the deep high temperature circumstance is disturbed during geothermal energy exploitation.The geothermal-related thermal fracturing of rock can significantly influence the geothermal energy exploitation.
In this study, a thermo-mechanical coupling algorithm based on the numerical manifold method was developed to analyze the thermo-elastic fracturing process of granular materials.The developed NMM algorithm represented the randomly shaped microstructure of granular materials using Voronoi polygon assemblies.Based on the Voronoi polygon assembly, a transient thermal conduction algorithm for granular materials was established, in which the heat interactions among Voronoi grains were reflected by a newly developed contact-based heat transfer model.Then, a damage-based contact fracture model was incorporated in the NMM algorithm to simulate the bonding-cracking effects between grains.The newly developed NMM algorithm was validated by two simple examples, the results of which agreed well with that of theoretical solutions.At last, the thermo-elastic fracturing processes induced by elevated temperature and cooling were simulated by the developed NMM algorithm.The obtained thermal fracture morphologies were consistent with those obtained by experiment test and available numerical simulations.
The developed numerical method is a coupled continuous-discontinuous numerical method.Since both the heat transfer in the solid grain and the heat interactions among grains are considered, the developed numerical method can better simulate the heterogeneous heat transfer process of rocklike granular materials.Additionally, as damage evolution and fracturing at grain interfaces is also considered, the developed numerical method is applicable to simulate the geothermal-related thermal fracturing process.However, in the numerical method developed, the effects of seepage, which is another significant influence of the geothermal energy exploitation, are not considered.In the future, the coupled thermal-hydraulic-mechanical numerical method will be developed and applied to analyze geothermal related problems.

Figure 2 .
Figure 2. Problem domain and MC (a) problem domain (b) MC.

Figure 2 .
Figure 2. Problem domain and MC (a) problem domain (b) MC.

Figure 2 .
Figure 2. Problem domain and MC (a) problem domain (b) MC.

Figure 5 .
Figure 5. Three types of possible contact pairs (a) angle to edge contact (b) edge to edge contact (c) angle to angle contact.

Figure 5 .
Figure 5. Three types of possible contact pairs (a) angle to edge contact (b) edge to edge contact (c) angle to angle contact.

Figure 6 .
Figure 6.Invasion of the three types of contact (a) angle to edge contact (b) edge to edge contact (c) angle to angle contact.

Figure 7 .
Figure 7. Representation of microstructure of granular materials (a) micro structure of rock like granular material [40] (b) Voronoi polygons (c) Numerical manifold method (NMM) model for granular material.

Figure 6 .
Figure 6.Invasion of the three types of contact (a) angle to edge contact (b) edge to edge contact (c) angle to angle contact.

Figure 7 .
Figure 7. Representation of microstructure of granular materials (a) micro structure of rock like granular material [40] (b) Voronoi polygons (c) Numerical manifold method (NMM) model for granular material.

Figure 7 .
Figure 7. Representation of microstructure of granular materials (a) micro structure of rock like granular material [40] (b) Voronoi polygons (c) Numerical manifold method (NMM) model for granular material.

Figure 8 .
Figure 8. Schematic of thermal conduction of granular materials.

Figure 9 .
Figure 9. Contact based heat transfer model.

Figure 10 .
Figure 10.Normal contact spring and tangential contact spring.Figure 10.Normal contact spring and tangential contact spring.

Figure 10 .
Figure 10.Normal contact spring and tangential contact spring.Figure 10.Normal contact spring and tangential contact spring.Energies 2018, 11, x FOR PEER REVIEW 11 of 21

Figure 11 .
Figure 11.A damage based contact fracture model (a) normal force-aperture relationship (b) tangential force-aperture relationship.

Figure 11 .
Figure 11.A damage based contact fracture model (a) normal force-aperture relationship (b) tangential force-aperture relationship.

Figure 12 .Figure 13 .
Figure 12.Simple contact heat transfer problem (a) geometry and boundary conditions (b) NMM cover system.
a square plane with edge length L = 1 m and initial temperature = 0 °C is subjected to thermal loading ̅ = 1 °C on the top boundary.The vertical displacement of the bottom boundary as well as the horizontal displacements of the right boundary and left boundary are fixed.The plane is subjected to plane strain conditions.

Figure 12 .
Figure 12.Simple contact heat transfer problem (a) geometry and boundary conditions (b) NMM cover system.

Figure 12 .Figure 13 .
Figure 12.Simple contact heat transfer problem (a) geometry and boundary conditions (b) NMM cover system.
a square plane with edge length L = 1 m and initial temperature = 0 °C is subjected to thermal loading ̅ = 1 °C on the top boundary.The vertical displacement of the bottom boundary as well as the horizontal displacements of the right boundary and left boundary are fixed.The plane is subjected to plane strain conditions.

Figure 13 .
Figure 13.Effects of contact thermal conductivity on the temperature (a) evolution of interface temperatures with η (b) temperature on center vertical line for case with very large η.

Figure 14 .
Figure 14.Simple thermo-mechanical coupling example (a) geometry and boundary conditions (b) NMM cover system.

Figure 15 .
Figure 15.Numerical results of thermo-mechanical coupling (a) variation of temperature with time for A and B (b) variation of lateral stresses with time for A and B.

Figure 14 .
Figure 14.Simple thermo-mechanical coupling example (a) geometry and boundary conditions (b) NMM cover system.

Figure 14 .
Figure 14.Simple thermo-mechanical coupling example (a) geometry and boundary conditions (b) NMM cover system.

Figure 15 .
Figure 15.Numerical results of thermo-mechanical coupling (a) variation of temperature with time for A and B (b) variation of lateral stresses with time for A and B.Figure 15.Numerical results of thermo-mechanical coupling (a) variation of temperature with time for A and B (b) variation of lateral stresses with time for A and B.

Figure 15 .
Figure 15.Numerical results of thermo-mechanical coupling (a) variation of temperature with time for A and B (b) variation of lateral stresses with time for A and B.Figure 15.Numerical results of thermo-mechanical coupling (a) variation of temperature with time for A and B (b) variation of lateral stresses with time for A and B.

Figure 16 .
Figure 16.Example of thermo-elastic fracturing induced by elevated temperature (a) geometry (b) NMM cover system.

Figure 16 .
Figure 16.Example of thermo-elastic fracturing induced by elevated temperature (a) geometry (b) NMM cover system.