Trans. Nonferrous Met. Soc. China 25(2015) 2367-2379

Thermo-mechanical coupled particle model for rock

Ming XIA

Computational Geosciences Research Centre, Central South University, Changsha 410083, China

Received 10 March 2014; accepted 21 July 2014

Abstract: A thermo-mechanical coupled particle model for simulation of thermally-induced rock damage based on the particle simulation method was proposed. The simulation results of three verification examples, for which the analytical solutions are available, demonstrate the correctness and usefulness of the thermo-mechanical coupled particle model. This model is applied to simulating an application example with two cases: one is temperature-independent elastic modulus and strength, while the other is temperature-dependent elastic modulus and strength. The related simulation results demonstrate that microscopic crack initiation and propagation process with consideration of temperature-independent and temperature-dependent elastic modulus and strength are different and therefore, the corresponding macroscopic failure patterns of rock are also different. On the contrary, considering the temperature-dependent elastic modulus and strength has no or little effect on the heating conduction behavior. Numerical results, which are obtained by using the proposed model with temperature-dependent elastic modulus and strength, agree well with the experimental results. This also reveals that the rock subjected to heating experiences much more cracking than the rock subjected to cooling.

Key words: particle simulation method; micromechanics; rock fracture; thermo-mechanical coupled model

1 Introduction

With the characteristics of high strength, rock has been widely used in the engineering structure, such as tunnels, rock slopes, dams and petroleum reservoirs [1]. To ensure the safety and stability of these engineering structures, many researchers have studied its behaviors under mechanical condition during the past few decades. However, it is often under the thermo-mechanical coupled condition in many scientific and engineering fields. For example, in environment protection engineering, geological sequestration of CO₂ into deep saline aquifers could become a promising solution for reducing the CO₂ concentration in the atmosphere. In a deep aquifer, the injected CO₂ flow in the near-well injection zone can create a temperature gradient that leads to a macroscopic failure as a result of microscopic crack initiation, propagation, and coalescence in the surrounding rock mass. This implies that thermally- induced rock damage should be considered in simulating the behavior of near-well injection zone. Similarly, the geological isolation has been commonly adopted to dispose high-level radioactive wastes (HLW) produced from nuclear power stations. The radioactive heat generated by the decay of nuclear elements in the buried HLW can create a high temperature gradient around the surrounding rock mass. Due to the long half-life features of radioactive elements, the maximum temperature around surrounding rock mass can reach 200-300 °C [2]. To ensure the safety and stability of HLW repositories, it is necessary to consider thermally-induced rock damage on the 10⁴-10⁶ year time-scale. Although all the above mentioned problems may involve two or more processes, such as thermal, hydrological, mechanical and chemical processes (THMC), the focus of this work is mainly on the thermally-induced rock damage associated with thermal and mechanical coupling processes. Thus, a thermo-mechanical coupled model for rock, which can reproduce its mechanical and thermal behavior accurately, is of great significance for simulating the rock mechanical and thermal behavior under thermo- mechanical condition.

From the macromechanical point of view, the acoustic emission technology can be used to locate the microscopic acoustic emission source in rocks for monitoring thermally-induced rock damage. In particular, DONG and LI [3,4] optimized and simplified the sensor location coordinates to find the analytical solution of the acoustic emission source location coordinates. On the other hand, from the micromechanical point of view, the failure of rock under thermo-mechanical condition is governed by microscopic crack initiation and propagation. However, the continuum-based numerical methods, such as the finite element method and boundary element method, are difficult, if not impossible, to simulate the predominant mechanism of rock failure processes associated with microscopic crack initiation, propagation and coalescence [5]. In this way, a micromechanical model of rock based on the particle simulation method [6,7], which can track the initiation and propagation of mechanically-induced microscopic crack, has been proposed. XIA and ZHOU [8] used the model to simulate the rock failure process. ZHAO et al [9] used the model to study the effects of shear on solute retardation coefficient in rock fractures. ZHAO [10] studied the gouge particle evolution in a rock fracture undergoing shear. YOON et al [11] simulated the fracture and friction of Aue granite under compression using clumped particle model. Recently, XIA and ZHAO [12] simulated rock deformation and mechanical characteristics using clump parallel-bond models. Although these studies enhance the understanding of rock failure processes, they are concentrated on the mechanical field. As stated previously, rock is often under thermo-mechanical condition in some engineering problems. While in this work, the particle simulation method is extended from the mechanical field to thermo-mechanical coupling fields. It should be noted that a thermo-mechanical coupled particle model for rock has been proposed previously [13]. However, it has following two disadvantages: 1) the microscopic parameters of this model cannot be directly determined from the related macroscopic ones; 2) the temperature-dependent elastic modulus and strength are not considered. Furthermore, WANNE and YOUNG [14] used this model to simulate thermally fractured granite, but the microscopic crack initiation and propagation processes cannot be simulated at the cooling stage.

To overcome the above mentioned two disadvantages, a new thermo-mechanical coupled particle model for rock in the thermo-mechanical coupling system was proposed. For the purpose of verifying the proposed model for simulating the thermo-mechanical coupling problem, three examples, for which analytical solutions are available, were used to show the correctness and usefulness of the proposed model. Finally, an application example was simulated to investigate the thermally-induced rock damage using the present model with two different cases: one is temperature-independent elastic modulus and strength, while the other is temperature-dependent elastic modulus and strength.

2 Thermo-mechanical coupled particle model

Even though the details of the model and the related computational algorithms can be found in Ref. [15], for the sake of completeness of this work, only some key mathematical formulations of the thermo-mechanical coupling particle model based on the particle simulation method are briefly given below. In this model, rock material was simulated as an assembly of particles, which were connected to each other through their bonds in the case of simulating mechanical deformation, but connected to each other through thermal pipes in the case of simulating heat conduction.

Figure 1 shows the mathematical formulations used in the thermo-mechanical coupling simulation. In Fig. 1, the simulation of mechanical deformation is shown on the left side, while the simulation of heat conduction is shown on the right side. For the equation of motion, F_i is the resultant force exerted on the mass center of the particle; m is the mass of the particle;  is the acceleration; g_i is the gravity acceleration; M_i is the rotation moment exerted on the mass center of the particle; I is the moment of inertia of the particle;  is the angular acceleration. For the force-displacement equation,  is the normal component of the temperature-dependent contact force;  is the temperature-dependent secant normal stiffness;  is the normal displacement;  is the incremental shear component of the temperature- dependent contact force;  is the temperature- dependent tangential shear stiffness;  is the incremental shear displacement. For the temperature calculation equation,  and  are the temperature at time t and , respectively;  is the thermal time step; c_V is the specific heat of the particle;  is the out-of-balance heat power of the heat reservoir. For the heat flow calculation equation, Q ^p is the heat power in pipe that is flowing out of the heat reservoir;  is the temperature difference between the two heat reservoirs at the each end of pipe; and  is the length of pipe.

Fig. 1 Thermo-mechanical coupling system

In the present model, a thermal pipe is associated with a specific contact. Once the thermal microscopic properties are assigned, subsequent loading or damage (in the form of bond breakages) will modify the number of active pipes, and thereby change the ability of the material to conduct heat. For example, as a bonded material is compressed and new pipes form, the macroscopic thermal conductivity will be increased, but as bonds break and associated pipes become inactive, the macroscopic thermal conductivity will be reduced. On the other hand, the thermally-induced force will affect the total force. The detailed formulations are given below.

When a contact bond is present at the contact associated with a pipe, it is assumed that only the normal component of the force carried by the contact bond, , will be affected by the expansion/contraction of the particle material due to the temperature change. Based on this assumption, a coupling relationship between the thermally-induced force increment () and the temperature increment () at the contact bond can be established. For an isotropic thermal expansion/contraction of the particle material, which has an effective length (), the corresponding change in the normal component of the force () carried by the contact bond can be expressed as follows:

    (1)

where  is the thermal volumetric expansion coefficient of the bond material, and the value of which is equal to the average value of the thermal volumetric expansion coefficients of the particles at two ends of the pipe associated with the bond;  is the change in the normal component of the temperature- dependent force () carried by the contact bond;  is the temperature increment, and the value of which is equal to the average temperature change of the two particles at the end of the pipe associated with the bond;  is the temperature-dependent contact normal stiffness; the value of the effective length () of the particle material is equal to the distance between the centroids of the two particles at the end of the pipe.

Note that in the thermo-mechanical coupling system, the contact force exerted on a particle contains two parts (i.e., the mechanically-induced force and the thermally- induced force increment), which can be evaluated as follows:

    (2)

           (3)

It needs to be pointed out that only the normal component of the force is affected by the temperature difference at two ends of the pipe associated with the contact bond.

Unlike the calculation of the normal component of the temperature-dependent contact force, the tangential component of the temperature-dependent contact force is computed in an incremental fashion. For any given time instant, t, the tangential component of the temperature- dependent contact force can be calculated by adding the temperature-dependent contact force increment expressed in Eq. (3) into the tangential component at the previous time, t-Δt_mechanical. Meanwhile, the tangential component of the contact force has to satisfy the Coulomb law of friction, which is determined by

        (4)

where  and  are the tangential components of the temperature-dependent contact force at t and t-Δt_mechanical, respectively; is the normal component of the temperature-dependent contact force at t;  is the corresponding mechanical incremental tangential contact displacement at t; Δt_mechanical is the mechanical time step in the numerical simulation and μ is the contact friction coefficient. Equation (4) holds true only when the normal component of the temperature- dependent contact force is greater than zero.

For the thermo-mechanical coupling model, crack initiation and generation are determined with the following crack criteria:

                              (5)

                             (6)

where  and  are the temperature- dependent normal and shear contact forces at a contact between any two particles; b_n(T) and b_s(T) are the temperature-dependent normal and tangential strengths of a contact bond when the bond average temperature is T.

3 Verification

It should be noted that any numerical models must be validated before they are put into scientific research and engineering application. For the verification purpose, three examples, which are also used as verification cases in the particle flow code in two dimensions (PFC^2D) [13], are used to conduct the validation for the present model. The first example (the transient thermal response of a sheet with constant temperature boundaries) is used to verify the particle simulation of heat conduction. The second example (the free expansion of a heated specimen) is used to verify the particle simulation of thermally- induced stress (strain). The third example (the transient response of a semi-infinite solid with applied heat flux) is used to verify the particle simulation of thermally- induced stress and heat conduction simultaneously. It should be noted that the macroscopic parameters, such as elastic modulus and strength, are treated as constant within these verification examples, which means that the microscopic parameters in the particle model are not changed with the temperature. The reason for this is that the corresponding analytical solution, which is used to verify the numerical results, is also derived with the temperature-independent case. Furthermore, it is very difficult to derive an analytical solution to the thermo-mechanical coupled problems when the elastic modulus and strength of a material is a function of temperature.

3.1 Transient thermal response of sheet with constant temperature boundaries

A planar sheet of width (L) is initially at a uniform temperature of 0 °C. The left side of the sheet is exposed to a constant temperature of 100 °C, while the right side is kept at a constant temperature of 0 °C. The sheet eventually reaches an equilibrium state at a constant heat flux and unchanging temperature distribution. The analytical solution for temperature within the sheet as a function of distance from the left side, x, and time, t, can be expressed as follows [16]:

       (7)

where T₁ is the temperature at the left side of the sheet (equal to 100 °C in this example);  is the thermal diffusivity that can be expressed as follows:

                                  (8)

where k is the thermal conductivity; ρ is the material density; c_V is the specific heat.

This heat conduction problem is one-dimensional, and can be modeled as a thin slice of material with constant temperatures applied to its left and right boundaries, and adiabatic conditions on the upper and lower boundaries to represent thermal symmetry planes. Figure 2 shows two different particle simulation models consisting of either a cubic packing or a hexagonal packing. The analytical solution is computed using the following macroscopic parameters: the density of the particle material is 1000 kg/m³; the specific heat is 0.2 J/(kg·°C); the thermal conductivity is 1.6 W/(m·°C). In order to match these macroscopic parameters, the thermal microscopic parameters are set as follows: thermal resistance coefficients are 0.3125 °C/(W·m) and 0.5413 °C/(W·m) for the cubic pack and hexagonal pack, respectively. The thermal time step used in the simulation is 62.5 s. Each model is run for 3125, 15625, 564937.5 s (for cubic packing) and 607437.5 s (for hexagonal packing), at which the thermal equilibrium state is reached in the particle simulation model. Figure 3 shows the normalized temperature (T/T₁) versus normalized distance (x/L) distribution for both the cubic-and hexagonal-packed specimens at different time. It is obvious that there is a good agreement between the numerical simulation results obtained from the particle simulation model and the analytical solutions obtained from Eq. (7). This indicates that the proposed model can be used to simulate heat conduction problems correctly.

Fig. 2 Two particle simulation models of different packing

Fig. 3 Comparison of simulated results with corresponding analytical ones

3.2 Free expansion of heated specimen

When a uniform temperature change of ΔT is applied to a square specimen with unconstrained boundaries, the thermally-induced strains can be measured and compared with analytical solutions for the square specimen consisting of an isotropic elastic continuum. The analytical solution for the thermally- induced strains can be expressed as follows:

     (9)

where α_t is the linear thermal expansion coefficient;  is Poisson ratio. The analytical solution is computed under the plane-stress condition with α_t=3×10^-6 °C^-1 in this study.

A particle simulation model of 2 m × 2 m is constructed and shown in Fig. 4(a). In this model, the average particle radius is 49.1 mm, resulting in a resolution of 20 particles across the specimen width. After the simulated specimen is constructed, the surrounding walls are removed, and the particle simulation model is run to allow free expansion to release the locked-in stresses until these stresses are less than 20 kPa in the simulated specimen. The strains were measured using a measurement circle, which has a radius of 0.95 m and is located at the center of the simulated specimen. The thermal expansion coefficient of the particle material is 3×10^{-6 °C-1}. To investigate the thermal expansion and contraction effects, both ΔT=  100 °C and ΔT=-100 °C are uniformly applied to all the particles within the particle simulation model. In this case, the mechanical time step used in the simulation is 5.542×10^-6 s. The mechanical parameters are as follows: the density of the particle material is 3155 kg/m³; the friction coefficient of the particle material is 0.5; the particle stiffness (in both the normal and the tangential directions) is 138 GN/m; both the macroscopic tensile strength and shear strength of the particle material are 100 MPa. The tensile strength and shear strength are artificially set large enough to prevent the failure of the particle material, because the analytical solution expressed in Eq. (9), which is used to verify the numerical results, is only valid in the elastic range. The model was then run in mechanical mode until it reaches a mechanically equilibrium state.

Fig. 4 Particle simulation model and simulated displacement distribution

In the case of ΔT=100 °C, the simulated numerical results for thermally-induced strains are , and , while in the case of ΔT=-100 °C, the simulated numerical results for thermally- induced strains are: ,  and . As shown in Figs. 4(b) and (c), the displacement distribution obtained from the particle simulation in the case of ΔT=100 °C displays an isotropic expansion, while the displacement distribution obtained from the particle simulation in the case of ΔT=-100 °C displays an isotropic contraction. As expected, the stresses in the particle simulation model after the application of each thermal increment remain near their initial values (i.e., the initial stresses just before the application of each thermal increment) during the whole particle simulation process. On the other hand, the corresponding analytical solutions are:  in the case of ΔT=100 °C and  in the case of ΔT=-100 °C. Compared with the analytical solutions, the relative errors of the simulated strains obtained from the particle simulation model are less than 0.1%. These small differences may be caused by the following two factors: 1) the particle simulation model before the application of the temperature increment is not completely stress free because the locked-in stresses resulted from the construction of the particle simulation model cannot be totally released; 2) the particle simulation model is not isotropic because the arbitrary packing of particles during the construction of the particle simulation model can produce a heterogeneity in the microstructure of material.

3.3 Transient response of semi-infinite solid with applied heat flux

The analytical solution for temperature within the semi-infinite solid as a function of distance from the surface, x, and time, t, can be expressed as follows [16]:

   (10)

where q is the applied heat flux on the surface of the semi-infinite solid (x=0); erfc is the complementary error function.

Under the plane strain condition, the analytical solution for the thermally-induced stress within the semi-infinite solid can be expressed as follows [17]:

                             (11)

where E is the elastic modulus of the solid material.

Since the transient response of semi-infinite solid with heat flux applied on its surface can be treated as a one-dimensional problem, the semi-infinite solid can be simulated as a thin slice of the solid material, as shown in Fig. 5. In this thin slab, a constant heat flux needs to be applied at its left boundary, and adiabatic conditions on the upper and lower boundaries need to be used to represent thermal symmetrical planes. In addition, the right boundary should be extended far enough to simulate the infinite extension of the semi-infinite solid in the x direction. The upper and lower boundaries of the thin slab need to be mechanically fixed only in the vertical direction to represent shear-free symmetry planes.

Fig. 5 Conceptual model for simulating thermal problem

Based on the conceptual model for simulating the transient response of semi-infinite solid with heat flux applied on its surface, a particle simulation model of 2 m × 25 m is constructed. In this model, the average particle size is 49.1 mm, resulting in a resolution of 20 particles in the height direction (y direction) of the thin slab. After the simulated specimen is constructed, the surrounding walls are removed, and the particle simulation model is run to allow free expansion to release the locked-in stresses until these stresses are less than 60 kPa in the simulated thin slab. The top and bottom of the slab are mechanically fixed in the y direction, while the right side of the slab is mechanically fixed in the x direction (by specifying the appropriate fixity condition for a thin layer of particles along these boundaries). The stress distribution is measured using a set of measurement circles, which have a radius of 0.9 m and are located with equal distance along the slab axis. The temperature distribution is measured by dividing the slab into 12 equal slices along the x axis, and computing the average temperature of the particles within each slice. In order to compare the simulated numerical results with the corresponding analytical ones, the distance, x, from the slab left boundary, is measured relatively to the original location of the left boundary (i.e., 12.5 m from the center of the particle simulation model).

The analytical solutions are computed using the following macroscopic parameters: the specific heat of the solid slab is 1015 J/(kg·°C); the linear thermal expansion coefficient of the solid slab is 3×10^{-6 °C-1}; the thermal conductivity of the solid slab is 3.5 W/(m·°C); the elastic modulus and Poisson ratio of the solid slab are 60.2 GPa and 0.0928, respectively; the tensile strength of the solid slab is 20 MPa; the shear strength of the solid slab is 200 MPa. To match these macroscopic parameters of the solid slab on the basis of method of ZHAO et al [18,19], the following thermal and mechanical microscopic parameters are used in the particle simulation model: the specific heat of the particle material is 1015 J/(kg·°C); the linear thermal expansion coefficient of the particle material is 3.711×10^{-6 °C-1}; the thermal resistance coefficient of the particle material is 3.126 °C/(W·m); the particle stiffness (in both the normal and the tangential directions) is 120.4 GN/m; the bond strength of the particle material is 20 MPa and 200 MPa in the normal and tangential directions, respectively. Note that the linear thermal expansion coefficient of the particle material is determined under the plane strain condition.

After the initial temperatures of all the particles in the particle simulation model are set to be 0 °C, a constant heat flux of 100 W/m² is applied to a thin layer of particles at the left end of the particle simulation model. The time step used in the simulation of heat conduction is 358.3 s, while it is 5.517×10^-6 s in the simulation of mechanical deformation. When this heat flux is lasted during a period of 3×10⁷ s, the particle simulation model is continuously running until it reaches a mechanical equilibrium state, at which the stress distribution can be obtained. Since no bond breakages occur during the particle simulation, the particle simulated slab behaves elastically during this period of the simulation time. Figure 6 shows both the temperature and stress distributions at the time of 3×10⁷ s. It is obvious that there is a good agreement between the numerical simulation results obtained from the particle simulation model and the analytical solutions obtained from Eqs. (10) and (11). Compared with the analytical solutions, the relative errors of the particle simulated numerical results are less than 5% at the left-most measurement circle in the particle simulated slab. This indicates that the proposed thermo-mechanical coupled particle model can be used to simulate heat conduction and thermally-induced deformation problems.

Fig. 6 Comparison of simulated results with corresponding analytical solutions at t=3×10⁷ s

4 Application

4.1 Description of numerical model and simulating process

After the proposed model is verified for simulating the thermo-mechanical coupling problem, it can be applied to the simulation of thermally-induced rock damage problems. As an application example, rock damage processes in the granite specimen with a heater placed in the center of a borehole, which was investigated through a laboratory experiment [20,21], are simulated to further demonstrate the applicability and usefulness of the proposed thermo-mechanical coupled particle model. Although WANNE and YOUNG [14] have simulated the same experiment using the PFC^2D, their model cannot reproduce the microscopic cracks initiation in the cooling phase. Meanwhile, the microscopic parameters in their model are all treated as a constant. As mentioned in the introduction, this treatment is not accurate and it is temperature-dependent. On the other hand, the microscopic parameters in their model cannot be determined directly from the corresponding macroscopic parameters. For these reasons, the proposed model, in which microscopic parameters are expressed as a temperature-dependent function, may predict more accurately the thermally- induced damage in a thermo-mechanical system.

As shown in Fig. 7, the particle simulation model with a diameter of 30 cm representing the outer boundary and a diameter of 3 cm representing a borehole is simulated using 3148 randomly distributed particles. To match macroscopic parameters and microscopic ones directly [18,19], the corresponding microscopic parameters adopted in this study are shown in Table 1.

Fig. 7 Geometry of computational model

Table 1 Microscopic parameters of model

After the particle simulation model reaches a static equilibrium state, a set of white particles which lie around the central borehole are identified as the heater particles (Fig. 7). Then the temperature of these particles is controlled at both the heating and steady stages. The temperature of outer boundary white particles is fixed to 20 °C for simulating room temperature. The initial temperate of the whole model is 20 °C. As stated in the experiment conducted by JANSER et al [20], an electrical resistance cartridge heater, which is placed in the center of the vertical borehole, is used to progressively heat the sample until it reaches the desired peak temperature. To realistically simulate the heating stage in the particle simulation model, the temperature applied to the heater particles is gradually increased until it reaches the assigned temperature, so that the same heating process in the laboratory experiment [20] can be simulated in the particle simulation model. The temperature increment is set to be 10 °C, while the assigned ultimate temperature is 200 °C. The simulation processes can be divided into three stages: 1) a heating stage from 0 to 2160 s, at which the temperature applied to the heater particles is gradually increased from 20 °C to 200 °C, so that the heater particles can reach the assigned temperature; 2) a steady stage from 2161 to 6870 s, at which the temperature of the heater particles is always 200 °C, so as to ensure that the heater particles can keep a thermal steady state at 200 °C; 3) a cooling stage from 6871 to 18000 s, at which the temperature of the heater particles can drop freely, but the temperature of the outer boundary particles is always 20 °C. It should be noted that three temperature monitoring particles are chosen to monitor the temperature change during the simulation process. The coordinates of these particles are A (25 mm, 0), B (50 mm, 0) and C (100 mm, 0), respectively.

In this research, two different cases, namely Case I (temperature-independent elastic modulus and strength) and Case II (temperature-dependent elastic modulus and strength), are considered in the numerical analysis to examine the effect of temperature-dependent elastic modulus and strength on the evolution of both temperature and microscopic cracks in a thermo- mechanical coupled system. In Case I, the elastic modulus and strength in all three stages are treated as a constant, which can be expressed as E(T)=E(T₀),  . While in Case II, from the experimental results of Lac du Bonnet granite [21], the elastic modulus and strength at temperatures between the normal temperature (20 °C) and 205 °C in the heating stage and steady stage can be fitted and described as follows:

      (12)

  (13)

  (14)

where E(T) is the macroscopic elastic modulus at temperature T;  is the macroscopic tensile strength at temperature T;  is the macroscopic shear strength at temperature T.

It should be noted that the thermally-induced damage is an irreversible process. For this reason, the value of both elastic modulus and strength in the cooling stage should be treated as a constant in Case II, which is the same as that at the end of the steady stage.

4.2 Computational simulation results

Figure 8 shows the time-history variation of particle temperature and the number of microscopic cracks with Case I and Case II. It is observed that although microscopic cracks are generated in the heating stage and cooling stage in both cases, the number of microscopic crack in Case II is less than that in Case I. For instance, in Case II, 87 microscopic cracks are generated in the heating stage, while in Case I, 126 microscopic cracks are generated in the heating stage. This indicates that the proposed model with the temperature-dependent elastic modulus and strength has a great effect on the microscopic crack initiation and propagation processes. On the other hand, 22 microscopic cracks are generated in the cooling stage in both Case I and Case II, indicating that the rock subjected to heating can experience much more cracking than that subjected to cooling. As shown in Fig. 8(a), it also can be observed that most of the microscopic cracks are generated in the heating stage in both cases. The reason for this is that the thermally-induced stress increases in the heating stage with the increase of time until its value is greater than the corresponding bond strength, so that most of the microscopic cracks are generated in the heating stage. There is almost no further damage in the steady stage because only few microscopic cracks are generated in this stage. However, in the cooling stage, microscopic cracks continue to generate only during a very short period of time, starting at the beginning of the cooling stage. This is because the temperature around the borehole decreases quickly during a very short period of time at the beginning of the cooling stage, so that the thermally-induced stress will be greater than the bond strength, and then microscopic cracks generate. Compared with the results by using the PFC^2D, which cannot be used to simulate the microscopic crack initiation and propagation in the cooling stage [14], the main advantage by using the present model is that it can be used to simulate the thermally-induced rock damage in the cooling stage, no matter whether the temperature-dependent elastic modulus and strength are considered or not. This demonstrated the correctness and usefulness of the proposed model for simulating thermally-induced rock damage problems.

Fig. 8 Time-history variation of temperature and microscopic cracks in particle simulation model

As shown by the temperature evolution of particles A, B and C in Fig. 8, the particle temperature evolution is similar, indicating that considering the temperature- dependent elastic modulus and strength or not has little effect on the heating conduction behavior. The thermally-induced damage is only formed at the local place in this computational model. Therefore, there is no prominent change in the topological structure of the active pipe networks in both Case I and Case II. Furthermore, the macroscopic thermal parameters in both Case I and Case II, which were used in the heat conduction and heat flow calculation in Fig. 1, are exactly the same, then the heating conducting behavior is very similar.

Fig. 9 Evolution processes of thermally-induced microcracks distribution

Figure 9 shows the evolution processes of the thermally-induced microscopic cracks during the whole particle simulation process. In the heating stage, when the temperature of the heater is 90 °C, microscopic cracks generate in both cases. Then, more microscopic cracks are randomly generated when the temperature of heater is 100 °C. As time goes on, these microscopic cracks tend to propagate and coalesce, leading to the formation of macroscopic fracture, which grow in a plane along the radial direction in the particle simulation model. It is noted that in the heating stage, the number of microscopic cracks increase slowly until the temperature of the heater reaches a critical value (100 °C), after which the number of microscopic cracks shows a dramatic increase and finally the thermally-induced macroscopic fracture is formed, and then the specimen is failed. At this time, the temperature of the heater reaches the thermal cracking threshold temperature for the Lac du Bonnet granite, so that a large amount of microscopic cracks are generated and failure occurs. As shown in Fig. 9, the final failure is initiated from the outer surface of the particle simulation model and propagates toward the borehole in both cases. This phenomenon is similar to what was observed from the laboratory experiment and is attributed to a temperature- gradient cracking mechanism within a brittle, anisotropic medium [20]. However, the macroscopic fracture patterns in these two cases are different. For instance, three macroscopic fractures are formed in Case I in the heating stage, while in Case II, only one macroscopic fracture is formed in the heating stage. After the macroscopic fracture is suddenly formed, few microscopic cracks are generated in the steady stage in both cases (Fig. 9). In the cooling stage, most of secondary thermally-induced microscopic cracks are generated and propagated around the borehole normal to the macroscopic fracture in Case I, while in Case II, they are mainly initiated from the inner surface of the particle simulation model and propagated along the radial direction. Then these microscopic cracks stop growing, which implies that macroscopic fracture reaches a steady state. This yields a characteristic macroscopic failure pattern as presented in Fig. 9. The related simulation results demonstrate that the present model with temperature-dependent elastic modulus and strength affects the microscopic crack initiation and propagation and therefore, may affect the macroscopic failure pattern of rock in the thermo-mechanical coupled system.

Fig. 10 Post-failure AE events of cylindrical sample (a) [22], and macroscopic failure pattern of numerical model in Case I (b) and Case II (c)

Compared the numerical results of final macroscopic failure pattern in both Case I and Case II with the experimental result [22], it can be found that the numerical result from Case II agrees well with the experimental result (Fig. 10). This demonstrated the correctness and usefulness of the proposed model with temperature-dependent elastic modulus and strength for simulating thermally-induced rock damage problems.

In order to illustrate the advantage of Case II over Case I, the numerical simulation in Case II was carried out by using two other different models, which have the same initial geometry, mechanical and thermal properties but different total number of particles. One model contains 6424 particles, while the other model contains 9521 particles. Figure 11 shows the microscopic crack distribution at the end of the cooling stage with these two different models. It can been seen that the macroscopic fracture in the model with 9521 particles has the highest resolution, within that the macroscopic fracture is formed with one main fracture plane in the heating stage. While in the cooling stage, there are three shorter fractures along the radial direction. However, as shown in Fig. 9(a), the macroscopic pattern in Case I mainly contains three fractures in the heating stage. Meanwhile, there is only one fracture formed along the circumferential direction in the cooling stage. Combined Fig.11 with Fig. 10, it can further demonstrate the correctness and usefulness of the proposed model with temperature-dependent elastic modulus and strength for simulating thermally-induced rock damage problems.

Fig. 11 Microscopic crack distribution in Case II

5 Conclusions

1) A thermo-mechanical coupled particle model for rock based on the particle simulation method is presented. Compared with the PFC^2D, the use of the present model has following three advantages: input microscopic (particle-level) parameters can be directly determined from the related macroscopic (sample-level) ones; temperature-dependent elastic modulus and strength are considered; microscopic crack initiation and propagation processes can be reasonably simulated in the cooling stage.

2) Three benchmark examples, for which analytical solutions are available, have been considered to verify the proposed model. Comparison of the simulated numerical results with the corresponding analytical solutions demonstrates that the proposed thermo- mechanical coupled particle model is suitable for simulating thermally-induced rock damage problems.

3) This model is applied to simulating an application example with two cases: one is temperature- independent elastic modulus and strength, while the other is temperature-dependent elastic modulus and strength. Meanwhile, numerical results with these two cases were compared. The related simulation results demonstrate that microscopic crack initiation and propagation process with consideration of temperature- independent and temperature-dependent elastic modulus and strength are different and therefore, the corresponding macroscopic failure patterns of rock are also different. On the contrary, considering the temperature-dependent elastic modulus and strength has no or little effect on the heating conduction behavior. Numerical results, which are obtained using the proposed model with temperature-dependent elastic modulus and strength, agree well with the experimental results. This also reveals that the rock subjected to heating experiences much more cracking than the rock subjected to cooling.

References

[1] HARRISON J P, HUDSON J A. Engineering rock mechanics. Part 2: Illustrative workable examples [M]. Oxford: Pergamon, 2000.

[2] HUDSON J A. Lesson learned from 20 years of UK rock mechanics research for radioactive waste disposal [J]. International Society for Rock Mechanics News Journal, 1999, 6(1): 27-41.

[3] DONG Long-jun, LI Xi-bing. Three-dimensional analytical solution of acoustic emission or microseismic source location under cube monitoring network [J]. Transactions of Nonferrous Metals Society of China, 2012, 22(12): 3087-3094.

[4] DONG Long-jun, LI Xi-bing, ZHOU Zi-long, CHEN Guang-hui, MA Ju. Three-dimensional analytical solution of acoustic emission source location for cuboid monitoring network without pre-measured wave velocity [J]. Transactions of Nonferrous Metals Society of China, 2015, 25(1): 293-302.

[5] JING L. A review of techniques, advances and outstanding issues in numerical modeling for rock mechanics and rock engineering [J]. International Journal of Rock Mechanics and Mining Sciences, 2003, 40(3): 283-353.

[6] CUNDALL P A, STRACK O D L. A discrete numerical model for granular assemblies [J]. Geotechnique, 1979, 29 (1): 47-65.

[7] ZHAO C, NISHIYAMA T, MURAKAMI A. Numerical modeling of spontaneous crack generation in brittle materials using the particle simulation method [J]. Engineering Computations, 2006, 23(5-6): 566-584.

[8] XIA M, ZHOU K. Particle simulation of the failure process of brittle rock under triaxial compression [J]. International Journal of Minerals, Metallurgy and Materials, 2010, 17(5): 507-513.

[9] ZHAO Zhi-hong, JING Lan-ru, NERETNIEKS I. Particle mechanics model for the effects of shear on solute retardation coefficient in rock fractures [J]. International Journal of Rock Mechanics and Mining Sciences, 2012, 52: 92-102.

[10] ZHAO Zhi-hong. Gouge particle evolution in a rock fracture undergoing shear: A microscopic DEM study [J]. Rock Mechanics and Rock Engineering, 2013, 46(6): 1461-1479.

[11] YOON J S, ZANG A, STEPHANSSON O. Simulating fracture and friction of Aue granite under confined asymmetric compressive test using clumped particle model [J]. International Journal of Rock Mechanics Mining Sciences, 2012, 49: 68-83.

[12] XIA M, ZHAO C. Simulation of rock deformation and mechanical characteristics using clump parallel-bond models [J]. Journal of Central South University, 2014, 21(7): 2885-2893.

[13] Itasca Consulting Group, Inc. Particle flow code in two dimensions (PFC^2D) [M]. Minneapolis: Minnesota, 2004.

[14] WANNE T S, YOUNG R P. Bonded-particle modeling of thermally fractured granite [J]. International Journal of Rock Mechanics and Mining Sciences, 2008, 45(5): 789-799.

[15] XIA M, ZHAO C, HOBBS B E. Particle simulation of thermally- induced rock damage with consideration of temperature-dependent elastic modulus and strength [J]. Computers and Geotechnics, 2014, 55: 461-473.

[16] CARSLAW H S, JAEGER J C. Conduction of heat in solids [M]. 2nd ed. London: Oxford Press, 1959.

[17] TIMOSHENKO S P, GOODIER J N. Theory of elasticity [M]. 3rd ed. New York: McGraw-Hill, 1970.

[18] ZHAO C, HOBBS B E, ORD A, PENG S, LIU L. An upscale theory of particle simulation for two-dimensional quasi-static problems [J]. International Journal for Numerical Methods in Engineering, 2007, 72(4): 397-421.

[19] ZHAO C, HOBBS B E, ORD A. Fundamentals of computational geoscience [M]. Berlin: Springer, 2009.

[20] JANSEN D P, CARLSON S R, YOUNG R P, HUTCHINS D A. Ultrasonic imaging and acoustic emission monitoring of thermally induced microcracks in Lac du Bonnet granite [J]. Journal of Geophysical Research, 1993, 98(B12): 22231-22243.

[21] CARLSON S R, JANSEN D P, YOUNG R P. Thermally induced fracturing of lac du bonnet granite [R]. Kingston: Queen’s University, 1993: 1-13.

[22] CHERNIS P J, ROBERTSON P B. Thermal cracking in Lac du Bonnet granite during slow heating to 205 °C [R]. Kingston: Atomic Energy of Canada Limited, 1993: 1-15.

一种岩石热-力耦合颗粒模型

夏 明

中南大学 计算地球科学研究中心，长沙 410083

摘  要：为了模拟岩石热损伤问题，基于颗粒法提出一个新的热-力耦合颗粒模型。利用3个存在解析解的标准算例验证该模型的正确性和有用性。然后，利用该模型模拟一个考虑2种不同情况的应用例子：一是考虑温度无关的弹性模量和强度，二是考虑温度相关的弹性模量和强度。计算结果表明：考虑温度相关和温度无关弹性模量和强度这两种不同情况时，微裂纹的萌生和扩展过程不同，最后导致所形成的宏观破裂形态不同；相反，考虑温度相关或无关弹性模量和强度对热传导行为影响不大。所提出的模型考虑温度相关弹性模量和强度时，相应的计算结果与试验结果吻合较好；与在冷却条件下相比，岩石在加热条件下能形成更多的微裂纹。

关键词：颗粒法；微观力学；岩石破裂；热-力耦合模型

(Edited by Xiang-qun LI)

Foundation item: Project (41372338) supported by the National Natural Science Foundation of China

Corresponding author: Ming XIA; Tel: +86-731-88830020; E mail: xiaming105@126.com

DOI: 10.1016/S1003-6326(15)63852-3