# A Multicomponent Particle Model and Linear Fitting Calibration Method for Heterogeneous Rocks

Authors: Qifeng Guo; Ying Zhang; Xun Xi; Xu Wu; Meifeng Cai
Journal: Advances in Civil Engineering
 (2019)

---

## Abstract

Deformation and failure of rock are very important to underground space engineering. In this paper, a new three-dimensional multicomponent bonded-particle model is developed for heterogeneous rocks. Granite samples from the Sanshandao Gold Mine are first analyzed using a microscope. Cylindrical and disc models consisting of four minerals (i.e., plagioclase, potash feldspar, quartz, and biotite) are built to simulate the behavior of the granite under compressive and tensile tests. To improve the calibration efficiency, a new method for determining the microparameters of minerals is proposed. Uniaxial compression, triaxial compression, and Brazilian tests are carried out in the lab. The failure pattern and stress-strain curves are obtained from numerical simulations and verified with those observed in the experiments. Furthermore, the Mohr–Coulomb and Hoek–Brown strength parameters are obtained and compared with the experimental results. The multicomponent particle model is shown to well reproduce the behavior of granite under the compressive and tensile tests. The multicomponent bonded-particle model significantly improves the accuracy of three parameters: the ratio of tensile strength to uniaxial compressive strength, the friction angle, and the parametermi of the Hoek–Brown criterion. The linear fitting calibration method provides a fast and efficient way to determine the microparameters of particles in heterogeneous rocks.


---

## Body

## 1. Introduction

Due to the rapidly increasing demands for city space, more and more underground constructions and infrastructures are designed and built in the world. The deformation and failure of rock are important to underground space engineering. Rock is normally considered as a heterogeneous geomaterial consisting of multiple minerals. The failure of rock is characterized by complicated micromechanical processes, including the nucleation, growth, and coalescence of microcracks. The behavior of rock deformation and failure has a large influence on the stability of geotechnical engineering projects. Landslides, tunnel collapse, and rock burst greatly threaten the lives of people and engineering safety [1, 2]. Therefore, a thorough investigation of the mechanisms of rock deformation and failure behavior is well justified.The deformation, damage accumulation, and crack patterns of rocks are commonly examined by tensile and compressive tests [3, 4, 5, 6]. Over the past two decades, numerical methods have been developed to provide increasing insight into the mechanical behaviors of rock samples. Lan et al. [7] developed a grain-based model considering grain-scale heterogeneity to simulate the microstructure of brittle rock by UDEC and found that crack initiation stress is controlled by the microscale geometric heterogeneity. Particle-based models, which were originally developed to simulate the behavior of cohesive frictional materials (e.g., soils and sands), are now often used for rock modeling. Compared with element- or mesh-based models (e.g., FEM, BEM, lattice model, and FDEM), particle-based models can reproduce the micromechanical behaviors of rock, such as microcrack nucleation, growth, and coalescence. Moreover, complex constitutive relationships are replaced with simple particle contact logic in particle-based models. Therefore, many researchers have resorted to particle-based models in simulating the deformation and failure behavior of rocks under tensile and compressive tests at both meso- and microscales. Potyondy et al. simulated the macromechanical parameters of Lac du Bonnet granite by using the bonded-particle model (BPM) [8]. Fakhimi et al. [9] employed the BPM to model biaxial compressive tests of notched sandstone and found that the model could reproduce the damage zone of samples as experimental results. Zhang and Wong [10–12] captured microcrack initiation, growth, and coalescence by using the BPM, with the simulation results being in good agreement with the experimental results. Zhang et al. [13] employed parallel bonded-particle models to simulate crack coalescence between two nonparallel flaws. Ji et al. [14] developed a new method to model the crack closure behavior of rock under compression by introducing preexisting opened microcracks with notional surfaces and gap in PFC. Kulatilake et al. [15] used the BPM to investigate the effect of joints on the uniaxial compressive strength (UCS) of samples. Park and Song [16] carried out direct shear tests by using the BPM and found that the simulated shear failure behavior is consistent with that observed in experiments. Poulsen and Adhikary [17] employed the BPM to simulate compressive tests of coal samples and obtained the stiffness and strength of the coal samples. More studies have proved that the BPM can reproduce the behavior of samples under mechanical tests and yield Young’s modulus, Poisson’s ratio, and UCS of rock. However, these models do not consider the multiple minerals in rock, and the existing models face three difficulties: (1) the ratio of tensile strength (TS) to UCS of rock, i.e., TS/UCS, yielded is larger than the real value; (2) the friction angle is smaller than the real value; and (3) the parameter mi of the Hoek–Brown (HB) failure criterion is smaller than the real value.For the particle-based modeling methodology, an extensive calibration based on experimentally measured macroscale properties is required to determine the contact parameters, which are very cumbersome and complex. Most studies commonly used the trial-and-error method to determine the contact parameters. The trial-and-error method is as follows: by keeping other parameters constant, investigating the effect of a single parameter or a combination of parameters on simulation results, comparing the macromechanical parameters with the experimental results, and then repeating the above steps again and again, the parameters that best predict the observed macroscale response can be obtained. To improve the efficiency of parameter determination, some researchers have resorted to using some new algorithms, e.g., the global optimization package SNOBFIT, orthogonal design or BP neural network. Zhou et al. [18] built a nonlinear BP neural network model for connecting macromechanical parameters to contact parameters. Liu et al. [19] used the orthogonal design method to determine the contact parameters and discussed the relationship between macromechanical parameters and contact parameters by a multivariate variance analysis. Wu and Xu [20] believed that any parameters in any calibration step will affect each other in terms of the macroresponse and hence developed a program to iteratively calibrate the contact parameters. Vallejos et al. [21, 22] developed an interactive calibration program that can change the contact parameters in real time until the parameters best predicting the macromechanical behavior are obtained. These new methods combining computation and mathematics have improved the efficiency of parameter determination. However, no universally applicable method with which to achieve this determination with high efficiency within a short amount of time exists.This paper attempts to develop a new multicomponent bonded-particle model (MBPM) for simulating rock micromechanical behavior. Granite samples from the Sanshandao Gold Mine are chosen as the work example. The mineral composition and content of the granite are first scanned and analyzed by using a microscope. Next, the MBPM is built to simulate the behavior of granite under tensile and compressive tests. A new calibration method for determining the microparameters of particles representing different minerals is then proposed. The effect of multiple minerals on the macromechanical response is considered. Compressive and Brazilian tests are carried out in the lab. Then, compressive and Brazilian tests based on the numerical model are carried out. The failure pattern and stress-strain curves are subsequently obtained from the numerical simulations and verified with the experimental results. Finally, the Mohr–Coulomb and HB strength parameters are obtained and compared with the experimental results.


## 2. A Multicomponent Bonded-Particle Model

In the traditional BPM model with a single mineral component, the mechanical parameters of uniaxial compression and tensile tests can easily be calibrated, at least in part. However, no matter the combination of normal and tangential bond strengths, the ratio of TS/UCS is always approximately 0.25, which is typically larger than the experimental results. Moreover, the internal friction angle and HB strength parametersmi are often smaller. This outcome likely occurs because all the particles have the same bond and stiffness parameters; thus, the difference in the strength of various minerals and the coordinated deformation of different minerals under loading cannot be reflected. In this paper, the minerals constituting the rock under observation are considered to build a multiphase BPM.Some granite samples were taken from the Sanshandao Gold Mine, the mining depth of which exceeds 1000 m. The samples were first scanned by using a microscope to analyze their mineral components. As shown in Figure1, the mineral composition and content can be described as follows: plagioclase (44%), potash feldspar (20.0%), quartz (32%), and biotite (4%). Based on the content of the four minerals, a MBPM is built by PFC and an in-house script is written in FISH. Figure 2 illustrates the multicomponent particle model, where blue, green, red, and yellow particles represent the plagioclase, potash feldspar, quartz, and biotite, respectively.Figure 1
Minerals of the Sanshandao granite under the microscope.Figure 2
Multicomponent bonded-particle model (MBPM) [10].As shown in Figure3(a), a cylindrical MBPM with a 50 mm diameter and 100 mm height is generated to simulate the uniaxial compression tests, which is in agreement with the recommendation on rock sample size by ISRM [23]. The model consists of 3858 spherical particles with varying diameters. The axial load is applied by 2D upper and lower rigid plates with a dimension of 75 mm × 75 mm (length and width, respectively). Figure 3(b) illustrates the MBPM for a triaxial compression test, which is achieved by applying confining pressure to the same model used in the uniaxial compression test. The confining stress is applied by using a cylindrical wall using a servo. Figure 3(c) shows a disc MBPM with a 50 mm diameter and 25 mm height. The model, which is generated to simulate a Brazilian disc split test, contains 3741 spherical particles with varying diameters. The particles are compacted, the internal structure is stable, and the particles are in full contact with each other after the equilibrium state is reached. The models are subjected to a constant displacement load rate until failure. Zhang et al. [11–13, 24] studied the loading rate effect on the simulation results by the particle model and found that the maximum loading rates for the uniaxial compression and Brazilian tests are 0.02 m/s and 0.01 m/s, respectively. To reduce the simulation time, the loading plate speeds used in this paper are the maximum values as the references.Figure 3
Numerical models for mechanical tests: (a) uniaxial compression; (b) triaxial compression; (c) Brazilian disc test.


## 3. Calibration of Micromechanical Parameters

When calibrating the micromechanical parameters for particle flow simulations, most previous studies considered only Young’s modulus, Poisson’s ratio, and UCS and ignored the damage threshold parameters (e.g., crack initiation stress and crack damage stress) [25, 26, 27, 28]. In this paper, we calibrated the micromechanical parameters to match Young’s modulus, Poisson’s ratio, UCS, TS, crack initiation stress, crack damage stress, cohesion, and friction angle of the model with respect to those of the rock sample.Table1 shows the microparameters in the particle flow simulation, which affect the macrobehaviors of the numerical sample and should be determined via the trial-and-error method. According to previous studies [20, 21, 22], five microparameters, i.e., parallel bond Young’s modulus, linearly contacted bonded Young’s modulus, ratio of normal to shear stiffness, magnification factor of the strength ratio, and friction angle, have significant effects on the behavior of rock samples. In addition, other microparameters, e.g., radius multiplier, maximum-to-minimum grain diameter ratio, and density, can be determined empirically using the references. As illustrated in Figure 4, the tensile Young’s modulus is calibrated by using the Brazilian tests, while the other five microparameters are calibrated by using the uniaxial compressive tests. The final microparameters are used for a triaxial compression test and verified against the experimental results. For the MBPM, the microparameters of particles representing different minerals must be determined, which is much more time-consuming compared to the traditional BPM model. Therefore, we propose a new calibration method with linear fitting that is explained by considering the calibration of the tensile Young’s modulus as an example. As shown in Figure 4, the initial values of the parallel bonded modulus of four minerals are determined empirically and then multiplied by n coefficients as the input for the Brazilian tests. The contact parameters between two particles representing two different minerals are simplified as the average values of the contact parameters of them. For example, the contact Young’s modulus between the quartz particle and the plagioclase particle is the average value of those in quartz particles and plagioclase particles. According to the n simulation results, the function relationship between the macrotensile Young’s modulus and the parallel bond Young’s modulus of particles representing individual mineral components can be obtained via the linear fitting method. Finally, the parallel bond Young’s modulus can be calculated by substituting the experimental tensile Young’s modulus into the function. Note that the other microparameters are constant when determining the parallel bond modulus. The other microparameters are also determined by using this method one by one. This method has the following characteristics: procedural flow, strong pertinence, stability, reliability, and ease of operation. Compared with the trial-and-error calibration method and iterative calibration method, the efficiency of the proposed method is greatly improved.Table 1
Microparameters of the contact bond model and parallel bond model.
Microparameter
Definition
Unit
E
¯
∗
Parallel bond Young’s modulus
GPa
E
∗
Particle contact Young’s modulus
GPa
k
n
/
k
s
Ratio of particle normal-to-shear stiffness
—
k
¯
n
/
k
¯
s
Ratio of parallel bond normal-to-shear stiffness
—
σ
¯
c
Parallel bond tensile strength
MPa
c
¯
Parallel bond cohesion
MPa
ϕ
¯
Friction angle
°
μ
Friction coefficient
—
λ
¯
Radius multiplier
—
R
min
Minimum particle radius
Mm
R
max
/
R
min
Particle size ratio
—
ρ
Particle bulk density
kg·m−3
n
Porosity
—Figure 4
The calibration procedure and microparameters determination.


## 4. Worked Examples

### 4.1. Mechanical Experiments on the Sanshandao Granite

The size of the experimental sample and that of the numerical model are the same. No crack is produced when cutting the specimen. Figure5(a) shows the WEP-600 universal testing machine at the University of Science and Technology Beijing, where the uniaxial compression and Brazilian tests are conducted. The loading rates are 0.5–1.0 MPa/s for the uniaxial compression test and 0.1–0.3 MPa/s for the Brazilian test. The triaxial compression tests are carried out by using the MTS-810 testing machine at the University of Minnesota, USA. The tests are controlled by a ring extensometer with a loading rate of 0.0005 mm/s. The confining pressures are 5 MPa, 10 MPa, 20 MPa, 30 MPa, 40 MPa, and 50 MPa. Eighteen specimens are divided into 6 groups for the uniaxial compression tests. The groupings and numbers of specimens for the Brazilian tests and triaxial compression tests are the same as those for the uniaxial compression tests.Figure 5
Mechanical test equipment: (a) uniaxial compression test and Brazilian disc test (WEP-600); (b) triaxial compression test (MTS-810).
(a)
(b)The crack initiation stress and crack damage stress are two important parameters indicating the strength of rock based on the stress-strain response curve [29]. The crack initiation stress of granite is defined as the stress when the sample starts to have a crack volume strain. The crack damage stress is defined as the stress when the volume strain changes from increasing to decreasing. According to the crack strain model [30, 31], the crack volume strain εvc is given as follows:(1)εvc=εv−εve,where εv and εve are the total volume strain and elastic volume strain, respectively.The elastic volume strain could be calculated as follows:(2)εve=1−2υEσ1−σ3,where σ1 is the axial stress, σ3 is the lateral stress, E is the elastic modulus, and υ is Poisson’s ratio in the elastic stage.The mechanical parameters of the granite obtained from the experimental tests are given in Table2. The granite parameters are calibrated to make the results of the numerical tests and those of the experimental tests consistent.Table 2
Mechanical properties of Sanshandao granite from experiments.
Property
Density (kg·m−3)
Tensile strength (MPa)
Tensile Young’s modulus (GPa)
Uniaxial compressive strength (MPa)
Comprehensive Young’s modulus (GPa)
TS/UCS
Crack initiation (MPa)
Crack damage (MPa)
Poisson’s ratio
Friction angle (°)
Cohesion (MPa)
Value
2686
12.40
40.99
94.37
43.98
0.13
37.75
81.16
0.20
49.96
17.19


### 4.2. Calibration of the Tensile Young’s ModulusEt

According to previous studies on the calibration of Young’s modulus [32, 33, 34, 35] and experiments on minerals [1], the initial microparameters of the four minerals in the granite are as shown in Table 3. The initial parallel bond Young’s modulus of the plagioclase (E¯x∗), potash feldspar (E¯j∗), quartz (E¯s∗), and biotite (E¯h∗) are 88.1 GPa, 96.8 GPa, 94.5 GPa, and 30 GPa, respectively. It should be mentioned that the value of bond Young’s modulus of every mineral represents the contact modulus between two particles representing the same mineral. The bond Young’s modulus between two different mineral particles is equal to the average value of them. The values of E¯x∗, E¯j∗, E¯s∗, and E¯h∗ are multiplied by five coefficients, i.e., 1.0, 0.8, 0.6, 0.4, and 0.2, while the other microparameters are kept constant; these are used as the input parameters for the Brazilian disc test. Figure 6 illustrates the axial stress-strain curves of the five Brazilian tests. The numerical tensile Young’s modulus can be measured based on the stress-strain curve. Furthermore, the relationship between the tensile Young’s model and parallel bond Young’s modulus of mineral particles can be obtained by using the linear fitting method, as illustrated in Figure 7 and equations (3)–(6):(3)plagioclase:E¯x∗=3.38746Et−1.21061,(4)potash feldspar:E¯j∗=3.72197Et−1.33016,(5)quartz:E¯s∗=3.63354Et−1.29856,(6)biotite:E¯h∗=1.1535Et−0.41224.Table 3
Initial values of the microparameters.
Microparameter
Plagioclase
Potash feldspar
Quartz
Biotite
E
¯
∗
E
¯
x
∗ = 88.1 GPa
E
¯
j
∗ = 96.8 GPa
E
¯
s
∗ = 94.5 GPa
E
¯
h
∗ = 30 GPa
E
∗
E
x
∗ = 88.1 GPa
E
j
∗ = 96.8 GPa
E
s
∗ = 94.5 GPa
E
h
∗ = 30 GPa
k
n
/
k
s
1
1
1
1
k
¯
n
/
k
¯
s
1
1
1
1
σ
¯
c
400 MPa
450 MPa
420 MPa
350 MPa
c
¯
800 MPa
900 MPa
840 MPa
700 MPa
ϕ
¯
25°
30°
27°
20°
μ
0.5
0.5
0.5
0.5
λ
¯
1
1
1
1Figure 6
Axial stress-strain curves from Brazilian disc tests for the different coefficients ofE¯∗.Figure 7
Relationships between the tensile Young’s modulus and parallel bond Young’s modulus of different minerals.The experimental tensile Young’s modulus of the granite is 40.99 GPa. Substituting that value into equations (3)–(6), the parallel bond Young’s modulus of the plagioclase, potash feldspar, quartz, and biotite are calculated to be 137.64 GPa, 151.23 GPa, 147.64 GPa, and 46.86 GPa, respectively.


### 4.3. Calibration of the Compressive Young’s ModulusEc

The parallel bond Young’s modulus of minerals is replaced with the values determined above. Then, the particle contact Young’s modulus for different minerals is used to calibrate the compressive Young’s modulus of the granite. Similar to the calibration of the tensile Young’s modulus, the initial values of the particle contact Young’s modulus of the four minerals are multiplied by five coefficients, i.e., 1.0, 0.8, 0.6, 0.4, and 0.2, while the other microparameters are held constant. Figure8 shows the axial stress-strain curves from the five uniaxial compression tests. The relationships between the compressive Young’s modulus and particle contact Young’s modulus of every mineral are shown in Figure 9 and expressed as follows:(7)plagioclase:Ex∗=2.25526Ec−2.77281,(8)potash feldspar:Ej∗=2.47798Ec−3.05077,(9)quartz:Es∗=2.4191Ec−2.97424,(10)biotite:Eh∗=0.76797Ec−0.9442.Figure 8
Axial stress-strain curves from uniaxial compression tests for the different coefficients ofE∗.Figure 9
Relationships between the compressive Young’s modulus and particle contact Young’s modulus of different minerals.By substituting the experimental compressive Young’s modulus of 43.98 GPa into equations (7)–(10), the following values of the particle contact Young’s modulus of the plagioclase, potassium feldspar, quartz, and biotite can be obtained: 96.41 GPa, 105.93 GPa, 103.42 GPa, and 32.83 GPa, respectively.


### 4.4. Calibration of Poisson’s Ratioυ

The ratio of particle normal-to-shear stiffness and the ratio of parallel bond normal-to-shear stiffness affect Poisson’s ratio of the model. Referring to previous research [33], the two ratios are considered to have the same value (i.e., kr=kn/ks=k¯n/k¯s). In addition, the difference in the stiffness ratio for different minerals is ignored. The value of the stiffness ratio is multiplied by five coefficients (0.5, 1.0, 1.5, 2.0, and 2.5), while the other microparameters are held constant. The lateral and axial stress-strain curves from the uniaxial compression simulations are plotted in Figure 10, from which Poisson’s ratio can be calculated. Furthermore, the relationship between the stiffness ratio and Poisson’s ratio is obtained by linear fitting, as shown in Figure 11. The regression formula is expressed as follows:(11)υ=0.07643+0.09282kr.Figure 10
Stress-strain curves from uniaxial compression tests for the different stiffness ratios.Figure 11
The relationship between Poisson’s ratio and the stiffness ratio.Poisson’s ratio obtained from the uniaxial compressive tests is 0.20. By substituting it into equation (11), the particle stiffness ratio is determined to be 1.33.


### 4.5. Calibration of the Failure Mode

The failure mode of numerical samples is known to be controlled by the ratio of bond TS to bond cohesive strength [20, 36, 37]. If the value of σ¯c/c¯ is small, the failure mode will be brittle. If the value of σ¯c/c¯ is large, ductile failure will occur in the simulations. The ratio of bond TS to bond cohesive strength is considered to be the same for different mineral components. Based on the microparameters determined above, the initial value of σ¯c/c¯ is multiplied by five coefficients (0.5, 1.0, 1.2, 1.5, and 2.0). The number of microcracks and the failure pattern under different values of σ¯c/c¯ are illustrated in Figure 12. Clearly, when the ratio of bond TS to bond cohesive strength is 0.5, the microcracks are almost all tensile cracks. When the ratio is larger than 1.0, shear microcracks appear first. Moreover, the larger the ratio is, the more the shear cracks dominate. The experimental results indicate that the failure of the granite under the uniaxial compression test involves both tensile and shear fractures but is dominated by tensile cracks. Therefore, the ratio of bond TS to bond cohesive strength is determined to be 1.0.Figure 12
Microcrack number developments and failure pattern under different values ofσ¯c/c¯: (a) 0.5, (b) 1.0, (c) 1.2, (d) 1.5, and (e) 2.0.
(a)
(b)
(c)
(d)
(e)


### 4.6. Calibration of the Uniaxial Compressive Strength

The initial parallel TS of the plagioclase, potash feldspar, quartz, and biotite is 400 MPa, 450 MPa, 420 MPa, and 350 MPa, respectively. A factorξ is introduced as the coefficient of parallel TS, which is set to 0.25, 0.5, 1.0, 1.5, and 2.0. Five uniaxial compression tests are carried out based on the determined microparameters and the five different parallel TS values. Figure 13 shows the axial stress-strain curves under different coefficients ξ. Clearly, the larger the parallel strength is, the larger the UCS is. The relationship between UCS and the coefficient ξ is illustrated in Figure 14. By linear fitting, the function is expressed as follows:(12)σc=7.65488+160.24488ξ.Figure 13
Axial stress-strain curves from uniaxial compression tests for the different coefficients ofξ.Figure 14
The relationship between the compressive strength and the strength coefficientξ.By substituting the experimental UCS of 94.37 MPa into equation (12), ξ is calculated to be 0.54. Furthermore, the parallel TS of each mineral can be obtained by multiplying the initial value by 0.54.


### 4.7. Calibration of the Friction Angle

Based on previous studies [35]; the initial friction angles of the plagioclase, potash feldspar, quartz, and biotite are determined to be 25°, 30°, 27°, and 20°, respectively. The friction angles of the minerals are multiplied by the coefficients of 0.25, 0.5, 1.0, 1.5, and 2.0, while the other microparameters are held constant. Figure 15 shows the axial stress-strain curves from the uniaxial compression tests for different friction angles. Clearly, the friction angle of minerals has no effect on the mechanical behavior of rock before reaching peak strength. In addition, the effect of the particle friction angle on postpeak behavior is little. Therefore, the friction angle coefficient is set to 0.5.Figure 15
Axial stress-strain curves from uniaxial compression tests for the different coefficients of the friction angle.


### 4.8. Results and Verification

The final microparameters obtained after performing calibration are given in Table4. The microparameters are used for the Brazilian test and compression test. The tensile strength in the simulations is determined by the stress-strain curve. The tensile strength of the bonds σ¯c was calibrated by matching the tensile strength of the rock σt obtained from direct tensile tests. The mineral of strength from high to low is potash feldspar, quartz, plagioclase, and biotite. The parameters for the biotite are significantly smaller than those for other minerals. It can be found in the simulations that the crack initiation is prone to occur at the points between the biotite particles or contact biotite particles. The crack initiation and postpeak response of the granite are mainly controlled by the heterogeneous microstructure [7]. Therefore, the multicomponent particle model could reproduce a more realistic mechanical behavior of the granite. Figure 16 shows the failure modes and axial stress-strain curves from the Brazilian experiment and numerical simulation. Clearly, the failure behaviors of the granite sample from the numerical simulation and the experiment are in good agreement. Notably, a brittle failure occurs in the experiments; thus, the stress-strain curve after the peak cannot be captured. Figure 17 illustrates the failure modes and stress-strain curves from the uniaxial compressive experiment and numerical simulation. According to the numerical results, the axial stress-strain curves before the peak are consistent with those of the experiments. Moreover, the failure patterns from the experiment and the numerical simulation are quite similar. Figure 18 shows the axial stress-strain curve and volumetric-axial strain curve from the simulations, which can be used to calculate the crack initial stress and crack damage stress. According to the crack strain model and equations (1) and (2), the crack initial stress and crack damage stress from the simulation (Figure 18) are approximately 36% of the UCS (34.66 MPa) and 89% of the UCS (85.72 MPa), respectively.Table 4
Microparameters of the multicomponent bonded-particle model after performing calibration.
Mineral
R
min (mm)
R
max
/
R
min
ρ (kg·m−3)
μ
E
∗ (GPa)
k
r
E
¯
∗ (GPa)
λ
¯
n
σ
¯
c (MPa)
c
¯ (MPa)
ϕ
¯ (°)
Plagioclase
1.5
1.66
2560
0.5
96.41
1.33
137.64
1
0.36
224
224
12.5
Potash feldspar
1.5
1.66
2630
0.5
105.93
1.33
151.23
1
0.36
252
252
15.0
Quartz
1.5
1.66
2650
0.5
103.42
1.33
147.64
1
0.36
235.2
235.2
13.5
Biotite
1.5
1.66
3050
0.5
32.83
1.33
46.86
1
0.36
196
196
10Figure 16
Comparison of experimental results and numerical results on the Brazilian disc test.Figure 17
Comparison of experimental results and numerical results on the uniaxial compression test.Figure 18
Axial stress-strain curve and volumetric-axial strain curve from the numerical uniaxial compression test.The comparison between the simulation results and experimental results on the basic macromechanical parameters is shown in Table5. As shown, the errors of the tensile and compressive Young’s moduli are 1.61% and 2.83%, respectively. The value of TS/UCS from the simulations is 0.12, which is 7.7% higher than that observed in the experiments. The numerical friction angle is just 0.5% larger than the experimental value. Clearly, the properties obtained from the simulations are in very good agreement with those from the experiments.Table 5
Comparison of mechanical properties of the Sanshandao granite from the experimental and numerical results.
Mechanical properties
Experimental result
Numerical result
Error (%)
E
c (GPa)
43.98
43.27
1.61
E
t (GPa)
40.99
42.15
2.83
υ
0.2
0.18
10.00
σ
ci (MPa)
37.75
34.66
8.20
σ
cd (MPa)
81.16
85.72
5.62
σ
c (MPa)
94.37
96.16
1.90
σ
t (MPa)
12.40
11.40
8.06
TS/UCS
0.13
0.12
7.7
C (MPa)
17.19
17.40
1.22
φ (°)
49.96
50.21
0.50


## 4.1. Mechanical Experiments on the Sanshandao Granite

The size of the experimental sample and that of the numerical model are the same. No crack is produced when cutting the specimen. Figure5(a) shows the WEP-600 universal testing machine at the University of Science and Technology Beijing, where the uniaxial compression and Brazilian tests are conducted. The loading rates are 0.5–1.0 MPa/s for the uniaxial compression test and 0.1–0.3 MPa/s for the Brazilian test. The triaxial compression tests are carried out by using the MTS-810 testing machine at the University of Minnesota, USA. The tests are controlled by a ring extensometer with a loading rate of 0.0005 mm/s. The confining pressures are 5 MPa, 10 MPa, 20 MPa, 30 MPa, 40 MPa, and 50 MPa. Eighteen specimens are divided into 6 groups for the uniaxial compression tests. The groupings and numbers of specimens for the Brazilian tests and triaxial compression tests are the same as those for the uniaxial compression tests.Figure 5
Mechanical test equipment: (a) uniaxial compression test and Brazilian disc test (WEP-600); (b) triaxial compression test (MTS-810).
(a)
(b)The crack initiation stress and crack damage stress are two important parameters indicating the strength of rock based on the stress-strain response curve [29]. The crack initiation stress of granite is defined as the stress when the sample starts to have a crack volume strain. The crack damage stress is defined as the stress when the volume strain changes from increasing to decreasing. According to the crack strain model [30, 31], the crack volume strain εvc is given as follows:(1)εvc=εv−εve,where εv and εve are the total volume strain and elastic volume strain, respectively.The elastic volume strain could be calculated as follows:(2)εve=1−2υEσ1−σ3,where σ1 is the axial stress, σ3 is the lateral stress, E is the elastic modulus, and υ is Poisson’s ratio in the elastic stage.The mechanical parameters of the granite obtained from the experimental tests are given in Table2. The granite parameters are calibrated to make the results of the numerical tests and those of the experimental tests consistent.Table 2
Mechanical properties of Sanshandao granite from experiments.
Property
Density (kg·m−3)
Tensile strength (MPa)
Tensile Young’s modulus (GPa)
Uniaxial compressive strength (MPa)
Comprehensive Young’s modulus (GPa)
TS/UCS
Crack initiation (MPa)
Crack damage (MPa)
Poisson’s ratio
Friction angle (°)
Cohesion (MPa)
Value
2686
12.40
40.99
94.37
43.98
0.13
37.75
81.16
0.20
49.96
17.19


## 4.2. Calibration of the Tensile Young’s ModulusEt

According to previous studies on the calibration of Young’s modulus [32, 33, 34, 35] and experiments on minerals [1], the initial microparameters of the four minerals in the granite are as shown in Table 3. The initial parallel bond Young’s modulus of the plagioclase (E¯x∗), potash feldspar (E¯j∗), quartz (E¯s∗), and biotite (E¯h∗) are 88.1 GPa, 96.8 GPa, 94.5 GPa, and 30 GPa, respectively. It should be mentioned that the value of bond Young’s modulus of every mineral represents the contact modulus between two particles representing the same mineral. The bond Young’s modulus between two different mineral particles is equal to the average value of them. The values of E¯x∗, E¯j∗, E¯s∗, and E¯h∗ are multiplied by five coefficients, i.e., 1.0, 0.8, 0.6, 0.4, and 0.2, while the other microparameters are kept constant; these are used as the input parameters for the Brazilian disc test. Figure 6 illustrates the axial stress-strain curves of the five Brazilian tests. The numerical tensile Young’s modulus can be measured based on the stress-strain curve. Furthermore, the relationship between the tensile Young’s model and parallel bond Young’s modulus of mineral particles can be obtained by using the linear fitting method, as illustrated in Figure 7 and equations (3)–(6):(3)plagioclase:E¯x∗=3.38746Et−1.21061,(4)potash feldspar:E¯j∗=3.72197Et−1.33016,(5)quartz:E¯s∗=3.63354Et−1.29856,(6)biotite:E¯h∗=1.1535Et−0.41224.Table 3
Initial values of the microparameters.
Microparameter
Plagioclase
Potash feldspar
Quartz
Biotite
E
¯
∗
E
¯
x
∗ = 88.1 GPa
E
¯
j
∗ = 96.8 GPa
E
¯
s
∗ = 94.5 GPa
E
¯
h
∗ = 30 GPa
E
∗
E
x
∗ = 88.1 GPa
E
j
∗ = 96.8 GPa
E
s
∗ = 94.5 GPa
E
h
∗ = 30 GPa
k
n
/
k
s
1
1
1
1
k
¯
n
/
k
¯
s
1
1
1
1
σ
¯
c
400 MPa
450 MPa
420 MPa
350 MPa
c
¯
800 MPa
900 MPa
840 MPa
700 MPa
ϕ
¯
25°
30°
27°
20°
μ
0.5
0.5
0.5
0.5
λ
¯
1
1
1
1Figure 6
Axial stress-strain curves from Brazilian disc tests for the different coefficients ofE¯∗.Figure 7
Relationships between the tensile Young’s modulus and parallel bond Young’s modulus of different minerals.The experimental tensile Young’s modulus of the granite is 40.99 GPa. Substituting that value into equations (3)–(6), the parallel bond Young’s modulus of the plagioclase, potash feldspar, quartz, and biotite are calculated to be 137.64 GPa, 151.23 GPa, 147.64 GPa, and 46.86 GPa, respectively.


## 4.3. Calibration of the Compressive Young’s ModulusEc

The parallel bond Young’s modulus of minerals is replaced with the values determined above. Then, the particle contact Young’s modulus for different minerals is used to calibrate the compressive Young’s modulus of the granite. Similar to the calibration of the tensile Young’s modulus, the initial values of the particle contact Young’s modulus of the four minerals are multiplied by five coefficients, i.e., 1.0, 0.8, 0.6, 0.4, and 0.2, while the other microparameters are held constant. Figure8 shows the axial stress-strain curves from the five uniaxial compression tests. The relationships between the compressive Young’s modulus and particle contact Young’s modulus of every mineral are shown in Figure 9 and expressed as follows:(7)plagioclase:Ex∗=2.25526Ec−2.77281,(8)potash feldspar:Ej∗=2.47798Ec−3.05077,(9)quartz:Es∗=2.4191Ec−2.97424,(10)biotite:Eh∗=0.76797Ec−0.9442.Figure 8
Axial stress-strain curves from uniaxial compression tests for the different coefficients ofE∗.Figure 9
Relationships between the compressive Young’s modulus and particle contact Young’s modulus of different minerals.By substituting the experimental compressive Young’s modulus of 43.98 GPa into equations (7)–(10), the following values of the particle contact Young’s modulus of the plagioclase, potassium feldspar, quartz, and biotite can be obtained: 96.41 GPa, 105.93 GPa, 103.42 GPa, and 32.83 GPa, respectively.


## 4.4. Calibration of Poisson’s Ratioυ

The ratio of particle normal-to-shear stiffness and the ratio of parallel bond normal-to-shear stiffness affect Poisson’s ratio of the model. Referring to previous research [33], the two ratios are considered to have the same value (i.e., kr=kn/ks=k¯n/k¯s). In addition, the difference in the stiffness ratio for different minerals is ignored. The value of the stiffness ratio is multiplied by five coefficients (0.5, 1.0, 1.5, 2.0, and 2.5), while the other microparameters are held constant. The lateral and axial stress-strain curves from the uniaxial compression simulations are plotted in Figure 10, from which Poisson’s ratio can be calculated. Furthermore, the relationship between the stiffness ratio and Poisson’s ratio is obtained by linear fitting, as shown in Figure 11. The regression formula is expressed as follows:(11)υ=0.07643+0.09282kr.Figure 10
Stress-strain curves from uniaxial compression tests for the different stiffness ratios.Figure 11
The relationship between Poisson’s ratio and the stiffness ratio.Poisson’s ratio obtained from the uniaxial compressive tests is 0.20. By substituting it into equation (11), the particle stiffness ratio is determined to be 1.33.


## 4.5. Calibration of the Failure Mode

The failure mode of numerical samples is known to be controlled by the ratio of bond TS to bond cohesive strength [20, 36, 37]. If the value of σ¯c/c¯ is small, the failure mode will be brittle. If the value of σ¯c/c¯ is large, ductile failure will occur in the simulations. The ratio of bond TS to bond cohesive strength is considered to be the same for different mineral components. Based on the microparameters determined above, the initial value of σ¯c/c¯ is multiplied by five coefficients (0.5, 1.0, 1.2, 1.5, and 2.0). The number of microcracks and the failure pattern under different values of σ¯c/c¯ are illustrated in Figure 12. Clearly, when the ratio of bond TS to bond cohesive strength is 0.5, the microcracks are almost all tensile cracks. When the ratio is larger than 1.0, shear microcracks appear first. Moreover, the larger the ratio is, the more the shear cracks dominate. The experimental results indicate that the failure of the granite under the uniaxial compression test involves both tensile and shear fractures but is dominated by tensile cracks. Therefore, the ratio of bond TS to bond cohesive strength is determined to be 1.0.Figure 12
Microcrack number developments and failure pattern under different values ofσ¯c/c¯: (a) 0.5, (b) 1.0, (c) 1.2, (d) 1.5, and (e) 2.0.
(a)
(b)
(c)
(d)
(e)


## 4.6. Calibration of the Uniaxial Compressive Strength

The initial parallel TS of the plagioclase, potash feldspar, quartz, and biotite is 400 MPa, 450 MPa, 420 MPa, and 350 MPa, respectively. A factorξ is introduced as the coefficient of parallel TS, which is set to 0.25, 0.5, 1.0, 1.5, and 2.0. Five uniaxial compression tests are carried out based on the determined microparameters and the five different parallel TS values. Figure 13 shows the axial stress-strain curves under different coefficients ξ. Clearly, the larger the parallel strength is, the larger the UCS is. The relationship between UCS and the coefficient ξ is illustrated in Figure 14. By linear fitting, the function is expressed as follows:(12)σc=7.65488+160.24488ξ.Figure 13
Axial stress-strain curves from uniaxial compression tests for the different coefficients ofξ.Figure 14
The relationship between the compressive strength and the strength coefficientξ.By substituting the experimental UCS of 94.37 MPa into equation (12), ξ is calculated to be 0.54. Furthermore, the parallel TS of each mineral can be obtained by multiplying the initial value by 0.54.


## 4.7. Calibration of the Friction Angle

Based on previous studies [35]; the initial friction angles of the plagioclase, potash feldspar, quartz, and biotite are determined to be 25°, 30°, 27°, and 20°, respectively. The friction angles of the minerals are multiplied by the coefficients of 0.25, 0.5, 1.0, 1.5, and 2.0, while the other microparameters are held constant. Figure 15 shows the axial stress-strain curves from the uniaxial compression tests for different friction angles. Clearly, the friction angle of minerals has no effect on the mechanical behavior of rock before reaching peak strength. In addition, the effect of the particle friction angle on postpeak behavior is little. Therefore, the friction angle coefficient is set to 0.5.Figure 15
Axial stress-strain curves from uniaxial compression tests for the different coefficients of the friction angle.


## 4.8. Results and Verification

The final microparameters obtained after performing calibration are given in Table4. The microparameters are used for the Brazilian test and compression test. The tensile strength in the simulations is determined by the stress-strain curve. The tensile strength of the bonds σ¯c was calibrated by matching the tensile strength of the rock σt obtained from direct tensile tests. The mineral of strength from high to low is potash feldspar, quartz, plagioclase, and biotite. The parameters for the biotite are significantly smaller than those for other minerals. It can be found in the simulations that the crack initiation is prone to occur at the points between the biotite particles or contact biotite particles. The crack initiation and postpeak response of the granite are mainly controlled by the heterogeneous microstructure [7]. Therefore, the multicomponent particle model could reproduce a more realistic mechanical behavior of the granite. Figure 16 shows the failure modes and axial stress-strain curves from the Brazilian experiment and numerical simulation. Clearly, the failure behaviors of the granite sample from the numerical simulation and the experiment are in good agreement. Notably, a brittle failure occurs in the experiments; thus, the stress-strain curve after the peak cannot be captured. Figure 17 illustrates the failure modes and stress-strain curves from the uniaxial compressive experiment and numerical simulation. According to the numerical results, the axial stress-strain curves before the peak are consistent with those of the experiments. Moreover, the failure patterns from the experiment and the numerical simulation are quite similar. Figure 18 shows the axial stress-strain curve and volumetric-axial strain curve from the simulations, which can be used to calculate the crack initial stress and crack damage stress. According to the crack strain model and equations (1) and (2), the crack initial stress and crack damage stress from the simulation (Figure 18) are approximately 36% of the UCS (34.66 MPa) and 89% of the UCS (85.72 MPa), respectively.Table 4
Microparameters of the multicomponent bonded-particle model after performing calibration.
Mineral
R
min (mm)
R
max
/
R
min
ρ (kg·m−3)
μ
E
∗ (GPa)
k
r
E
¯
∗ (GPa)
λ
¯
n
σ
¯
c (MPa)
c
¯ (MPa)
ϕ
¯ (°)
Plagioclase
1.5
1.66
2560
0.5
96.41
1.33
137.64
1
0.36
224
224
12.5
Potash feldspar
1.5
1.66
2630
0.5
105.93
1.33
151.23
1
0.36
252
252
15.0
Quartz
1.5
1.66
2650
0.5
103.42
1.33
147.64
1
0.36
235.2
235.2
13.5
Biotite
1.5
1.66
3050
0.5
32.83
1.33
46.86
1
0.36
196
196
10Figure 16
Comparison of experimental results and numerical results on the Brazilian disc test.Figure 17
Comparison of experimental results and numerical results on the uniaxial compression test.Figure 18
Axial stress-strain curve and volumetric-axial strain curve from the numerical uniaxial compression test.The comparison between the simulation results and experimental results on the basic macromechanical parameters is shown in Table5. As shown, the errors of the tensile and compressive Young’s moduli are 1.61% and 2.83%, respectively. The value of TS/UCS from the simulations is 0.12, which is 7.7% higher than that observed in the experiments. The numerical friction angle is just 0.5% larger than the experimental value. Clearly, the properties obtained from the simulations are in very good agreement with those from the experiments.Table 5
Comparison of mechanical properties of the Sanshandao granite from the experimental and numerical results.
Mechanical properties
Experimental result
Numerical result
Error (%)
E
c (GPa)
43.98
43.27
1.61
E
t (GPa)
40.99
42.15
2.83
υ
0.2
0.18
10.00
σ
ci (MPa)
37.75
34.66
8.20
σ
cd (MPa)
81.16
85.72
5.62
σ
c (MPa)
94.37
96.16
1.90
σ
t (MPa)
12.40
11.40
8.06
TS/UCS
0.13
0.12
7.7
C (MPa)
17.19
17.40
1.22
φ (°)
49.96
50.21
0.50


## 5. Triaxial Compression Tests and Discussion

The microparameters determined are used to simulate the mechanical behavior of samples under triaxial compressive tests. The confining pressures are set to 5 MPa, 10 MPa, 20 MPa, 30 MPa, 40 MPa, and 50 MPa for the simulations and experiments. The axial stress-strain curves from the experiments and simulations are illustrated in Figure19. As shown, the larger the confining pressure is, the larger the peak strength is. The stress-strain curves from the simulations before the peak stress is reached are in good agreement with those from the experiments. Because of the brittle or quasibrittle failure and the loading procedure and data capture procedure used in the experiments, the experimental curves after reaching the peak strength are more volatile than those of the simulations.Figure 19
Comparisons of deviatoric stress-strain curves from the experimental results and numerical results of triaxial compression tests.Using the Mohr–Coulomb failure criterion, the strength envelope curves of the experiments and simulations are obtained, as shown in Figure20. Clearly, when the confining pressure is in the range of 0–50 MPa, the strength envelope curves are almost linear, and the regression equations are as follows:(13)experiment:σ1=7.53116σ3+112.30462,simulation:σ1=7.63534σ3+113.04038.Figure 20
Mohr–Coulomb strength envelopes from experimental results and numerical results.It is important to compare the simulations results with the experimental results based on the parameters of the HB failure criterion. Most previous studies could not calibrate the strength parameter of the HB failure criterion. The HB failure criterion is expressed as follows [38, 39]:(14)σ1=σ3+σcmiσ3σc+sα,where mi is the parameter accounting for the rock hardness; s is the measure of how fractured the rock is; and α reflects the system’s bias.The fitting equations obtained from the experiments and simulations are given in equations (15) and (16), respectively. Figure 21 presents the strength envelope curves obtained from the experiments and simulations, which are in good agreement. Clearly, the mi parameter from the simulations is 17% smaller than that from the experiments:(15)σ1=σ3+74.0380810.91621σ374.03808+1.327080.79835,(16)σ1=σ3+96.76979.02805σ396.7697+1.142640.86144.Figure 21
Hoek–Brown strength envelopes from experimental results and numerical results.


## 6. Conclusions

In this paper, we proposed a new three-dimensional MBPM for heterogeneous rocks. The models, consisting of four mineral components, were built based on a microscopic analysis of Sanshandao granite. The cylindrical and disc models were built to simulate the behavior of Sanshandao granite under uniaxial and triaxial compression tests and Brazilian disc tests. To improve the efficiency of calibration, a new method for determining the microparameters of the minerals was proposed. The uniaxial compression, triaxial compression, and Brazilian tests were carried out in the lab. By using the linear fitting method, the microparameters were determined quickly and accurately. After the calibrations, the multicomponent particle model well reproduced the macrobehavior of the rock. The findings of our study are as follows:(1)
The MBPM becomes a linear contact model after cementation failure. The parallel bond parameters affect the behaviors under both tension and compression, while the particle contact parameters affect only the compressive behavior. Therefore, a positive correlation exists between the particle contact Young’s modulus and compressive Young’s modulus of rocks.(2)
The ratio of bond TS to bond cohesive strength greatly affects the failure mode and crack evolution under compression. When the ratio is larger than 1.0, shear microcracks appear first. Moreover, the larger the ratio is, the more dominant the shear cracks are. The friction angle of minerals has no effect on the mechanical behavior of rock before the peak strength is reached. In addition, the effect of the particle friction angle on postpeak behavior is little.(3)
The Mohr–Coulomb and HB strength envelopes from the numerical results are in good agreement with those from the experiments. The value of TS/UCS from the simulations is 0.12, which is 7.7% higher than that from the experiments. The numerical friction angle is just 0.5% larger than the experimental value. In addition, themi parameter from the simulations is 17% smaller than that from the experiments. The MBPM significantly improves the accuracy of the above three parameters, and the linear fitting calibration method provides a fast and efficient way to determine the microparameters of particles for heterogeneous rocks.



---

*Source: 1509719-2019-04-01.xml*