DDA Simulation Study on Mechanical Failure of Heterogenous Rock

Heterogeneity is an important characteristic that affects the mechanical behavior of rock. In the present work, a statistical rock mesoheterogeneity model based on the Weibull distribution function is introduced into the discontinuous deformation analysis (DDA) method to simulate the mechanical failure of heterogeneous rock, in which the general heterogeneity degree is controlled by a heterogeneity index and the mechanical property of each subblock element is randomly assigned. Brazilian disc and uniaxial compressive rectangular specimens are simulated as examples. Results show that it is more reasonable to consider the heterogeneity of elasticity properties (the elastic modulus and Poisson’s ratio) and strength properties (the tensile strength, cohesion, and friction angle) simultaneously in the heterogeneity model. It is also shown that with a larger heterogeneity index, which means a lower degree of heterogeneity, the reproducibility of the macroscopic response curves of a specimen gets better, while the exact cracking always differs but with less scattered cracks, and the global fracturing failure pattern and mode are weakly influenced by the heterogeneity. Moreover, with the increase in the heterogeneity index, the macroscopic equivalent modulus and strength get larger and approach those of a homogeneous specimen. This work indicates the importance of heterogeneity for rock mechanical behaviors including the macroscopic equivalent response and the fracturing failure. By the subblock DDA method to simulate fracturing realistically, the fracturing failure process of heterogeneous rock can be successfully reproduced, which builds good foundation for the simulation study of heterogeneous rock fracturing in practical problems, e.g., coal and rock fracturing in fluidization mining in the future.


Introduction
Rock is generally regarded as a kind of brittle material with nonlinearity, heterogeneity, and anisotropy. The mechanical behavior of rock including the macroscopic equivalent response and the fracturing failure is strongly influenced by the characteristic of heterogeneity. Thus, the heterogeneity of rock should be carefully considered in practice, e.g., in the fluidization mining [1,2], in which the mechanical behavior of coal and rock is a key concern [3][4][5]. In the mesoscale, heterogeneity means that inner rock has different properties (elasticity properties, strength properties, etc.) at different locations. In view of its good economy, abundant information, and easy manipulation, numerical simulation is a favorable approach for the mechanical behavior study of heterogeneous rock.
To describe heterogeneous rock mechanical behavior in numerical simulations, the heterogeneity model is required to be introduced in the numerical method first. In statistical modeling, the Weibull distribution function is considered a suitable function to reflect the mesoheterogeneity of rock and is usually used to describe heterogeneous rock. In previous studies, some simulation methods can already simulate the rock heterogeneity, such as the RFPA, DEM, and DDA method [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21]. Based on the assumption that the rock mesoelement parameters meet the Weibull distribution, Tang et al. [6,7] introduced the heterogeneity model in the RFPA method, which was proven to be feasible. With the RFPA method, Zhu et al. [8] simulated the concrete and rock fracture by using the maximum tensile strain criterion. Liu et al. [9] further analyzed the effect of the heterogeneity index on the rock cracking process and cracking strength by the RFPA method. Besides, Chen and Konietzky [10] introduced the heterogeneity model into the DEM method and analyzed the rock mechanical behavior effectively. Molladavoodi and Rahimirezaei [11] introduced a constitutive model based on the Weibull distribution into the DEM to study the influence of heterogeneity on rock mechanical behavior. Tang et al. [12][13][14][15] introduced a spatial correlation length factor into the traditional Weibull distribution and studied the heterogeneity character of concrete failure using an equivalent probabilistic model. Jiao et al. [16] introduced the rock heterogeneity model into the DDA method and investigated the influence of heterogeneity on rock fracture but only considered the heterogeneity of the elasticity properties in the heterogeneity model without considering the effects of the strength parameters' heterogeneity on rock fracturing. Moreover, based on the SPH method, Sun et al. [17] established a numerical model to simulate the fracture process of heterogeneous rock. Chen et al. [18] studied the heterogeneous mechanical response of tight conglomerate under triaxial compression through experiment and FEM analysis and found that the deformation of rock samples before reaching peak strength is linear-elastic. Liu et al. [19] proposed a thermomechanical coupling simulation method to reflect the heterogeneity of rock mass under high temperature environment. Zhang et al. [20] conducted a series of uniaxial compression simulation using the COM-SOL software and discussed the influence of rock heterogeneity. Zhang et al. [21] investigated the damage evolution and mechanical properties of heterogeneous rock based on a grain-based finite-discrete element model.
What is more, due to the stochastic character of the heterogeneity model when the Weibull distribution function is involved, under the same degree of heterogeneity, namely, with the same heterogeneity index, the simulated results of a problem may be different. So, it may be not convincing when some conclusion is derived from only one time of simulation. However, there were only a few studies that considered this issue. Ma et al. [22] aimed to model the axial loading test and checked the reproducibility of simulation results, but they only carried out two simulations for each example, which seems not to be persuasive and systematic enough. Tang et al. [12] carried out a series of numerical uniaxial compression tests on rock specimens. Thereinto, five samples were generated in each group based on the Monte Carlo method, which were used to reflect the randomness of mesoproperties in rock. But they did not give a specific discussion on the reproducibility.
Since the DDA method is good at solving discontinuous large deformation and large displacement problems [23][24][25], this method is adopted to study the mechanical failure of heterogeneous rock in the present study. Firstly, based on the Weibull distribution function, the heterogeneity model is introduced into the DDA method, in which the heterogeneity of elasticity properties as well as strength properties is all considered so that the rock heterogeneity characteristics can be described more properly. Thereafter, a number of Brazilian disc and uniaxial compressive specimens are simulated to check the reproducibility of simulation results with the same heterogeneity index. The influences of the heterogeneity index on the macroscopic equivalent response and the fracturing of the specimens are also investigated through the simulation examples. DDA is a numerical method that computes the mechanical response of discrete deformable block systems. The subblock DDA fracturing modeling approach can be adopted for rock fracturing simulation, in which subblock elements are glued by artificial joints of high strength to represent continuum, and tensile/shear fracturing can take place along artificial joint surfaces. In the improved subblock DDA fracturing modeling method [26] adopted here, the fracturing along artificial joints is judged based on the adjacent subblock stress levels rather than the contact stresses between subblocks, which greatly reduces the influence of artificial joint orientations and subblock element distributions on the fracturing path and failure strength modeling results.

Introduction of the Heterogeneity Model.
According to previous studies by other researchers, the Weibull distribution is suitable for describing rock mesoheterogeneity characterization. Here, an element-level mesoheterogeneity model is introduced into the DDA based on the Weibull distribution function. Because of the linear-elastic constitutive relation and first-order displacement function used in the DDA method, the stress-strain relationship for a DDA block is linear. However, like that in the RFPA [6,7], by introducing the heterogeneity model, the nonlinear characteristics of rock at the macrolevel could be reproduced. From the Weibull distribution formula, 10º Figure 1: Geometrical model of the disc.
2 Geofluids the cumulative distribution function can be obtained by integration as Let Q ðxÞ = y and simplify; it is derived that where y is the random number between 0 and 1, m is the heterogeneous index, x 0 is the homogeneous material property parameter, and x is the heterogeneous material property parameter. Substitute the homogeneous elastic modulus E, Poisson's ratio μ, friction angle θ, cohesion c, and tensile strength σ t into x 0 , respectively; a corresponding heterogeneous material property parameter x related to the random number y and heterogeneous index m can be obtained. By giving the heterogeneous index value m and inputting a homogeneous material property value x 0 , a statistical heterogeneous material property value x can be generated for each DDA subblock element by a random number created by the computer. In this way, the heterogeneity of each elasticity or strength parameter can be realized for a DDA model.

Considering Both Elasticity and Strength Heterogeneity
The previous work done by other researchers generally only considered the heterogeneity of the elasticity (the modulus and Poisson's ratio) and ignored the strength (the tensile strength, cohesion, and friction angle) heterogeneity. However, in reality, the strength parameter values vary with locations like that of the elasticity parameters in real rock, which means that the heterogeneity of strength parameters should be considered. Therefore, the strength parameters are processed to be heterogeneous in this paper to better simulate real rocks. Brazilian disc examples are simulated to verify this issue below. The Brazilian disc is assumed to be linear-elastic with a radius of 50 mm, an elastic modulus of E = 20 GPa, and Poisson's ratio of μ = 0:2. There is a platform of a disc center angle of 10 degrees at each loading end of the disc which can reduce the stress concentration and prevent the disc rotation, as illustrated in Figure 1. The strength parameters of the rock are taken as a tensile strength of σ t = 5 MPa, a cohesion of c = 20 MPa, and a friction angle of θ = 30°. The friction angle between the rigid loading plates and the disc is 0°. The maximum displacement ratio, the time step, and the spring stiffness in the DDA simulation are 5 × 10 −4 , 1 × 10 −6 s, and 100 GPa, respectively. The disc is subjected to linear velocity loading of 2 mm/s at its two ends by the rigid  3 Geofluids loading plates. In order to measure the macroscopic equivalent response of the disc, one more plate is set between the upper loading plate and the disc, and two measuring points are located in this plate. The disc is discretized into triangular elements with an element number of 4633. Figure 2 shows the DDA simulation results of Brazilian discs when the crack just goes through the specimen with different degrees of heterogeneity, namely, with different heterogeneity index values. Two groups of numerical examples are carried out: the heterogeneity of elasticity and strength parameters is all taken into account in the first group (condition 1: a 1 -d 1 ), and only the heterogeneity of elasticity parameters is considered in the second group (condition 2: a 2 -d 2 ). Comparing the failure modes under the two conditions, when the heterogeneity index is the same, more scattered cracks appear in the discs under condition 1, and the cracks propagate through the disc faster under condition 1, as shown in Table 1. This phenomenon is more obvious when the value of m is smaller. That is to say, with the same heterogeneity index, the distribution of the crack in the disc is more dispersed and the disc is more prone to destruction when considering the heterogeneity of elasticity and strength parameters simultaneously. Figure 3 depicts the macroscopic response curves (the equivalent stress measured by the measuring plate) of the Brazilian disc when index m = 1:5 and m = 200 under different conditions. It can be found that when the index m is the same, compared with condition 2, in condition 1, cracks in the disc are generated earlier, the macroscopic equivalent strength of the disc is lower, and the macroscopic nonlinear characteristics are more obvious. Meanwhile, when the index m increases, the macroscopic equivalent strength under the two conditions gets closer, indicating a weaker effect of the heterogeneity.
Apparently, the above results show that it is more reasonable to consider the heterogeneity of elasticity properties and strength properties simultaneously in the heterogeneity model, which is more in line with the reality of real rock. So, in the following examples, heterogeneity of both the elasticity and strength parameters is taken into consideration.

Reproducibility Investigation with the Same Heterogeneity Index
According to the heterogeneity generation method described above, the elastic and strength parameter values of each subblock element in a DDA model are calculated based on a random number; that is, the strong and weak elements are randomly distributed in the model, which means that even with the same heterogeneity index, the created DDA model will be different every time, and of course the simulation results will show a certain degree of divergence. Therefore, it is very important to study the reproducibility of simulation results when considering the heterogeneity. In this section, a series of numerical examples including the Brazilian disc and compressive rectangular specimens are simulated to reveal this issue, as shown in Table 2.  where the geometrical parameters and mechanical parameters of the disc are the same as those in Section 3. It can be seen clearly that with the same heterogeneity index, the exact cracking in every simulation is different. However, compared with the random distribution of the scattered cracks, the reproducibility of the main crack path is quite obvious.
The main crack is formed in the center of the disc as a crack band in every model, and all the main cracks initiate from the disc center and propagate along the loading direction. With a larger heterogeneity index, the main crack band is narrower and clearer, and less scattered cracks are generated throughout the model. In general, the heterogeneity does lead to the discreteness of the exact cracking results, but the global failure pattern and mode are similar. Moreover, it can be seen that the macroscopic response curves are linear-elastic at the early stage, and before the quick fall at the peak strength, the curves show a little extent of nonlinearity, especially under the small heterogeneity index condition. When the heterogeneity index is the same, compared with the obvious discreteness of the curve after reaching peak strength, the reproducibility of the curve in the elastic stage is higher.
By comparing the macroscopic response curves under different heterogeneity index m values, it can be found that (1) (2) (3) (4)   5 Geofluids the curves match better as m increases. To reflect the influence of heterogeneity on the reproducibility of macroscopic response results quantitatively, the peak strength discreteness under the same m value is defined by where ξ is the discreteness of peak strength, Δσ t is the difference between the highest peak strength and the lowest peak strength, and σ t is the mean value of the highest peak strength and the lowest peak strength. The peak strength discreteness with different heterogeneity index values is listed in Table 3. When index m = 1:5, 2.5, and 5, the discreteness of peak strength is about 1/3, 1/4, and 1/6, respectively. Obviously, as the heterogeneity index increases, the reproducibility of macroscopic response gets higher.

Uniaxial Compressive Rectangular Specimen Examples.
In this part, the uniaxial compressive failure of rectangular rock specimens with different heterogeneous degrees is       7 Geofluids and a friction angle of θ = 30°. The friction angle between rigid loading plates and the specimen is 0°, and a measuring plate is also placed between the upper loading plate and the specimen. The specimen is discretized into 11486 triangular elements. Other material and DDA control parameters are the same as those of the disc examples above. Figures 8 and 9 show the simulation results when index m = 5 and 20, respectively. It can be seen clearly that the global failure pattern and mode are similar under all conditions; namely, the specimen finally fails due to the formation of inclined main cracks. However, with the increase in the heterogeneity index, i.e., the decrease in heterogeneity, the inclination angle of the main crack path increases, and there are fewer scattered cracks generated throughout the specimen.
For the macroscopic response, it can be seen that with the same heterogeneity index, the curves at the linearelastic stage almost duplicate, but with the increase in loading displacement, the dispersion of the curves becomes obvious at the softening stage, which leads to different residual strength. Generally, the degree of the coincidence of the curves when m = 20 is better than that when m = 5. As listed in Table 4, when m = 20, the discreteness of the peak and residual strength is 2.9% and 58.7%, respectively, which is smaller than that of 5.0% and 73.9% when m = 5, respectively.
From the above disc and rectangular specimen simulation examples, it can be concluded that the increase in the heterogeneity index, which means the decrease in the degree of heterogeneity, will lead to better reproducibility of the simulation results regarding the macroscopic equivalent response, while the global fracturing failure pattern and mode are weakly influenced by the heterogeneity, but a lower degree of heterogeneity will bring less scattered cracks in the model. It is also found that with the same heterogeneity index value, the reproducibility of the macroscopic equivalent response for the rectangular specimen is better than that for the disc specimen.

Influence of Heterogeneity on Mechanical Failure of Rock
A series of Brazilian discs are simulated to further investigate the influence of heterogeneity on the mechanical   Figure 10.
As can be seen from the figure, obviously, with the increase in the value of heterogeneity index m, the failure pattern of the heterogeneous discs gradually gets close to that of the homogeneous disc. Specifically, when m = 1:5, many scattered cracks are generated in the disc and a wide crack band is formed in the center of the disc, which indicates the high degree of heterogeneity. As the value of heterogeneity index m increases, the randomly distributed scattered cracks decrease obviously, and the main crack band becomes narrower and clearer along the loading direction. When m = 200, the fracturing failure pattern is very close to that of the homogeneous disc, indicating that a heterogeneity index of 200 corresponds to a very low degree of heterogeneity.
For the macroscopic response curves, when the value of heterogeneity index m is small, the curve reaches the peak strength through a relatively short linear-elastic stage, which means that cracks are generated in the disc at an early loading stage. It is also found that a small value of the heterogeneity index leads to a small peak strength as well as a small equivalent elastic modulus of the model. When the value of index m is 200, the macroscopic response curve almost duplicates that of the homogeneous disc, which also proves that a heterogeneity index of 200 corresponds to a very low degree of heterogeneity.

Conclusions
In the present work, a heterogeneity model based on the Weibull function is introduced into the subblock DDA method, in which the general heterogeneity degree is controlled by a heterogeneity index and the mechanical property of each subblock element is randomly assigned. A large number of Brazilian disc and compressive rectangular rock specimen examples are simulated to study the effect of heterogeneity on rock mechanical failure.
First, as compared with the condition in which only the elastic parameters' heterogeneity is considered, when the heterogeneity of elasticity and strength is both taken into consideration, a lower macroscopic equivalent strength of the specimen will be derived, and more scattered cracks will be generated throughout the specimen. Second, the reproducibility of the macroscopic response curve increases with the increase in the heterogeneity index, while the exact fracturing routes always differ but with less scattered cracks. Under the same heterogeneity index value, the linearelastic stage of the macroscopic response curves is less affected by the heterogeneity as compared with that of the nonlinear stage, the global failure pattern and mode are not obviously influenced by the heterogeneity, and the reproducibility of the macroscopic equivalent response for the rectangular specimen is better than that for the disc specimen. Third, a smaller value of the heterogeneity index means a decrease in the macroscopic equivalent strength and more random scattered cracks. As the heterogeneity index value increases, the macroscopic equivalent modulus and strength of rock will be higher, and the mechanical response gradually gets close to that of homogeneous rock.
In the present study, the numerical simulation reproduces some interesting and important phenomena that can be found in rock specimen experiments in the aspects of macroscopic response and fracturing failure patterns. Future research, e.g., by correlating the simulation results with experimental results, could be carried out to ascertain the heterogeneity index value for different rocks and the threshold heterogeneity index value by which the rock can be regarded homogeneous.

Data Availability
Data is available upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.  10 Geofluids