Numerical Investigation on the Evolution of Mechanical Properties of Rock Affected by Micro-Parameters

Abstract

Investigating the micro-parameters of rock is vital for understanding the macro-properties of rock, such as the uniaxial compressive strength (UCS), Young’s modulus, failure patterns, etc. In this paper, based on the experimental results of rock material, a parallel-bond model in three-dimensional particle flow code (PFC3D) was applied to investigate the effects of the joint action of bond stiffness ratio and bond stress ratio on macro-properties of rock. The uniaxial compressive strength, stress–strain relationships, and failure characteristics, as well as underlying compression and failure mechanisms, in the process of parameter calibration, were systematically studied. The results indicated that the interaction of several micro-parameters would obviously change the response characteristics of the macro-properties of the model. The mechanism of the effects of various micro-parameters on the macro-properties of the model was further revealed. The change of the micro-parameters would change the strength and stress state of the bond between particles. The research results could promote the understanding of the failure mechanism of rock and improve the efficiency of micro-parameter calibration and the accuracy of calibration results.

1. Introduction

The discrete element method (DEM) was first introduced by Cundall in 1971 [1]. Since this time, the DEM has been widely applied in fields, such as mining and geotechnical engineering [2,3,4,5,6,7,8,9]. Particle flow code (PFC) is one of the widely used particle DEMs. The PFC is convenient for dealing with problems in non-continuous media and can effectively simulate non-continuous phenomena (e.g., cracking, media separation), which has greatly aided research into the damage, fracture evolution, and failure mechanism of rocks [10,11,12,13,14].

The first step in a PFC simulation is to select rational micro-parameters according to the material macro-responses. The accuracy of the simulation results depends on the selection of micro-parameters [15,16]. Many studies have explored the relationship between micro- and macro-parameters and achieved substantial success. Potyondy and Cundall [17] analyzed the sensitivity of the results of biaxial, triaxial, and Brazilian splitting tests of the PFC model to micro-parameters. Yoon et al. [18,19] proposed the calibration and optimization method of PFC micro-parameters based on the design of experiments method. Yang et al. [20,21] analyzed the quantitative relationship between the macro-properties and micro-parameters of the calculation model by using a two-dimensional particle flow code (PFC2D) parallel-bond model. Cong [22] combined the triaxial loading and unloading test of marble to study the micro-parameters applicable to the mechanical analysis of rock materials. Shi et al. [23,24,25] established a two-dimensional discrete element model, which could reflect the mineral components of rocks, and studied the influence of mineral components and microdamage distribution on the macro-properties of rocks. Ding et al. [26,27] studied the influence of model scale and particle size distribution on PFC simulation results.

However, most of these studies are based on PFC2D, and few are based on a three-dimensional particle flow code (PFC3D). Substantial differences have been shown to exist between PFC2D and PFC3D simulation results [17,20,27,28,29]. Additionally, plane strain and plane stress conditions can be applied to a wide variety of applications [30,31,32,33], but some engineering and experiment problems in the field cannot be simplified into plane stress or plane strain problems [19,26,34,35]. The research on the relationship between macro-properties and micro-parameters generally only considers the relationship between macro-properties and a single micro-parameter, while the research on the influence of multi-parameter coupling on macro-properties is insufficient. However, limited research shows that the coupling of multiple parameters has a more significant impact on the macro characteristics of the model [21,36].

In the PFC bond model, the contact-bond model (CBM) and parallel-bond model (PBM) are widely used. CBM can be regarded as a spring with certain normal and shear stiffness, which binds the particle element at the contact. PBM can be regarded as a series of springs with normal and shear stiffness, which are evenly distributed on a certain width of the contact surface and can resist the action of force and moment at the same time [19,37]. In PBM, the bond failure will immediately reduce the macro-stiffness, which is more in line with the brittle failure characteristics of rock [10,19,37]. Therefore, PBM is selected to simulate the mechanical behavior of hard rock in the process of failure.

In this study, the parallel-bond model in PFC was applied to investigate the effects and effect mechanisms of bond normal strength, shear strength, and stiffness ratio on the macro-properties of the PFC3D model of rocks. Section 2 briefly reviews the relationships between micro-parameters and macro-properties in the parallel-bond model. Section 3 describes the test scheme, Section 4 presents the test results, and the key problems found in this study are discussed in Section 5. The main conclusions are given in Section 6.

2. Parallel-Bond Model and Its Micro-Parameters

PFC3D simulates the mechanical behavior of continuous medium materials by defining the contact model to bond mutually independent particle units. Several studies have addressed the calibration of micro-parameters for various contact models in PFC and also analyzed the relationships between these micro-parameters and macro-properties [18,38,39,40,41,42].

The following parameters in PBM have been shown to exert major effects on the macro-properties [22,26,28,37,38]: $R_{\max }/R_{\min }$ ($R_{\min }$ is minimum particle radius in the entire population of particles, $R_{\max }$ is particle radius in the entire population of particles), $R_{\min }$, $L/d$ ($L$ is the smallest characteristic model length, such as the diameter of a cylindrical specimen, or the smallest width of a parallelepiped specimen, $d$ is the median particle diameter), $k^{\ast }$ (bond normal-to-shear stiffness ratio, $k^{\ast }=k_n/k_s$, $k_n$ is normal stiffness, $k_s$ is shear stiffness), $E_c$ (effective modulus, $E_c=k^{\ast }{L}_R/A$, $L_R$ is the sum of particle radius on both sides of the bond, $A=2r$ (2D) or $A=\pi r^2$ (3D), $r$ is the minimum radius of particles on both sides of the bond), μ (grain contact friction coefficient), ${\overline{k}}^{\ast }$ (bond normal-to-shear stiffness ratio, ${\overline{k}}^{\ast }={\overline{k}}_n/{\overline{k}}_s$, ${\overline{k}}_n$ is bond normal stiffness, ${\overline{k}}_s$ is bond shear stiffness), ${\overline{E}}_c$ (bond effective modulus, ${\overline{E}}_c={\overline{k}}^{\ast }{L}_R$), ${\overline{\sigma }}_n$ (bond normal stress), ${\overline{\tau }}_n$ (bond shear stress), and λ (bond width multiplier). $\overline{R}=\lambda r$, $\overline{R}$ is the bond radius. By setting λ, the bond radius is enlarged or reduced.

Among these parameters, ${\overline{k}}^{\ast }$, ${\overline{\sigma }}_n$, and ${\overline{\tau }}_n$ have significant effects on model macro-properties. The ${\overline{k}}^{\ast }$ value has a substantial effect on the model Poisson’s ratio and Young’s modulus [17,19,22]. With increasing ${\overline{k}}^{\ast }$, the Poisson’s ratio increases, while Young’s modulus decreases [19,20,21,22]. In contrast, the existing studies on the effect of ${\overline{k}}^{\ast }$ on uniaxial compressive strength (UCS) are insufficient. Most researchers have not considered this effect or only considered it on the condition that other micro-parameters are unchanged. They have not considered coupling effects between ${\overline{k}}^{\ast }$ and other micro-parameters on UCS, as well as the model deformation and failure characteristics [18,19,20,21,22,38,43].

Table 1. Summary of the PFC study on the relationship between micro-parameters and UCS.

References	Value Range of Bond Normal-to-Shear Stiffness Ratio	Influence Degree of Bond Normal-to-Shear Stiffness Ratio on UCS	Note
Deng, S.X. et al. [19]	0.6–3.6	Not significant	PFC3D
Wang, Y.N. et al. [38]	0–0.8	Not considered	PFC3D
Zhao, G.Y. et al. [21]	0.1–6	Not considered	PFC2D
Yang, B.D. et al. [20]	0.6–2	Not significant	PFC2D
Cong, Y. et al. [22]	1–5	Not significant	PFC2D
Yoon, J. [18]	0.5–10	Not significant	PFC2D
Xu, X.M. et al. [43]	1–100	Not considered	PFC3D

Uniaxial compressive strength (UCS), three-dimensional particle flow code (PFC3D), two-dimensional particle flow code (PFC2D)

shows the results of the relationship between ${\overline{k}}^{\ast }$ and UCS in different researchers’ studies.

Previous studies have shown that UCS of the model is also highly dependent on ${\overline{\sigma }}_n$ and ${\overline{\tau }}_n$, both of which have a positive effect [18,19]. Zhao [21] reported that when ${\overline{\sigma }}_n/{\overline{\tau }}_n$ > 0.5, the strength of the model was mainly affected by ${\overline{\tau }}_n$, and when ${\overline{\sigma }}_n/{\overline{\tau }}_n$ < 0.5, it was mainly affected by ${\overline{\sigma }}_n$. However, changes in ${\overline{\sigma }}_n$ and ${\overline{\tau }}_n$ generally do not affect Young’s modulus or Poisson’s ratio. Earlier studies have primarily investigated the effect of a particular single micro-parameter on the macro-properties of the model, whereas the combined effects of multiple micro-parameters on the macro-properties of the model have not been taken into account. Additionally, some studies [18,19,20,21,23,24,25,26,27] set $k^{\ast }$ = ${\overline{k}}^{\ast }$ as a default to facilitate analyses, even though these values have been shown to have considerable effects on the model macro-properties under different conditions [7,22].

In summary, the effects of ${\overline{k}}^{\ast }$, ${\overline{\sigma }}_n$, and ${\overline{\tau }}_n$ in PFC3D on model macro-properties require a more thorough investigation. It is of great significance to improve the efficiency of micro-parameter calibration and the accuracy of calibration results.

3. Experiment and Calibration

In order to make the simulation results conform to the physical properties of rock materials, the rock was first tested under uniaxial compression, and then the micro-parameters calibrated according to the test results were used as the benchmark for further study.

The rock material used in the test was sandy mudstone (Figure 1a), which was taken from −437.62 m of buertai mine, Shandong mining area, Inner Mongolia, China. The main mineral compositions of rock materials were analyzed by D/max-2500 X-ray diffractometer in the School of materials science and engineering, Shandong University. The content of clay mineral was 34.8%, quartz was 51.9%, K-feldspar was 5.6%, and plagioclase was 7.7%. The rock was made into a standard rock sample with a height of 100 mm and a radius of 25 mm. The test loading rate was 1 mm/min.

A cylindrical model with a height of 100 mm and a radius of 25 mm was constructed using PFC3D, and the damping coefficient used in the experiment was 0.7. Specimens were loaded using the top and bottom walls during uniaxial compression. A loading rate of 0.1 m/s was used in the simulation process to ensure the comparability and prevent interference of the loading rate on the results [26,44] (Quasi-static modeling with PFC3D usually made usage of a high local damping coefficient to efficiently remove kinetic energy from the system. Therefore, typical quasi-static deformation could be performed at a much higher rate than in real experiments [26]. Because of the high damping coefficient used, the loading rate in PFC was different from that used in the actual test [42]. The loading rate of 0.1 m/s could meet the requirements of PFC3D uniaxial compression simulation).

The numerical simulation software was PFC3D5.0, purchased from Itasca Consulting Co., Ltd. (Wuhan, China). Based on the experimental data of sandy mudstone, a numerical model simulating sandy mudstone was established, and a set of micro-parameters to match the numerical results with the measured mechanical properties was determined. Calibrated micro-parameters were used to systematically explore the effects of ${\overline{\sigma }}_n$, ${\overline{\tau }}_n$, and ${\overline{k}}^{\ast }$ in PFC3D on the macro-properties of rocks and the underpinning failure mechanisms by varying ${\overline{\sigma }}_n/{\overline{\tau }}_n$ and ${\overline{k}}^{\ast }$. The micro parameters were adjusted by a trial-and-error method using the following specific steps: (1) According to the uniaxial compression test, the UCS, E, μ of sandy mudstone were obtained; (2) The numerical model was established by PFC3D software, and a group of micro-parameters was set up according to the experience, and the uniaxial compression test was simulated; (3) The UCS, E, and μ obtained from the experiment and simulation were compared; (4) The reasonable micro-parameters could be obtained by adjusting the micro-parameters until the simulation results were close enough to the experimental results. The calibrated micro-parameters of the PFC3D model are presented in

Table 2. Calibrated micro-parameters of the PFC3D model.

Micro-Parameter	Description	Calibrated Value
$R_{\min }$/mm	Minimum particle radius	1
$R_{\max }/R_{\min }$	Particle radius ratio	1.66
$E_c$/GPa	Effective modulus	3
$k^{\ast }$	Normal-to-shear stiffness ratio	3.3
μ	Particle friction coefficient	0.5
${\overline{E}}_c$/MPa	Bond effective modulus	3
${\overline{k}}^{\ast }$	Bond normal-to-shear stiffness ratio	3.3
${\overline{\sigma }}_n$/MPa	Normal bond strength (Mean ± standard deviation*)	9 ± 3
${\overline{\tau }}_n$/MPa	Shear bond strength (Mean ± standard deviation)	15 ± 3
λ	Bond width multiplier	1

* According to the test results, it was set by PFC code

.

Results of PFC simulations and uniaxial compression test are presented in

Table 3. Experimental data and simulation results.

Test Method	Uniaxial Compression Strength (σc/MPa)	Young’s Modulus (E/GPa)	Poisson’s Ratio (ν)	Sample Size (R × H/mm)
Laboratory test	31.91	4.61	0.241	25 × 100
PFC simulation	33.07	4.57	0.248	25 × 100

. As shown in Figure 1b, the stress–strain characteristics of the DEM model were compared with the experimental data obtained from sandy mudstone under uniaxial compression. The simulation results were basically consistent with the test data. The difference in stress–strain curve was caused by the loading mechanism of numerical simulation [42]. In the process of sample formation, the particle had achieved self-balance under the effect of its gravity acceleration, and the particle had reached an ideal state through effective contact, while the microfracture of the actual rock sample led to initial compression phase.

To study the effects of ${\overline{\sigma }}_n$, ${\overline{\tau }}_n$, and ${\overline{k}}^{\ast }$ on the macro-properties of the model, a total of 10 different bond normal-to-shear stress ratios, ${\overline{\sigma }}_n/{\overline{\tau }}_n$ < 1 (0.1, 0.2, 0.4, 0.6, 0.8) and ${\overline{\sigma }}_n/{\overline{\tau }}_n$ > 1 (1, 2, 3, 4, 6), were set based on the micro-parameters presented in Table 2. When simulating rock material, the reasonable range of values of ${\overline{k}}^{\ast }$ was 1–11 [45,46]. Thus, each group ${\overline{\sigma }}_n/{\overline{\tau }}_n$ corresponded to ten groups ${\overline{k}}^{\ast }$. The value of ${\overline{k}}^{\ast }$ was set to 1.11, 2, 3, 4, 5, 6, 7, 8, 9, 10.11. The influence on the rock macro-properties of the PFC3D model was monitored by varying ${\overline{\sigma }}_n/{\overline{\tau }}_n$ and ${\overline{k}}^{\ast }$. Additionally, the stiffness ratio factor (KK = $k^{\ast }/{\overline{k}}^{\ast }$) was taken as KK = 0.1, 1, 2, 4, 6, 8, or 10 to study the effects of $k^{\ast }$ and ${\overline{k}}^{\ast }$ on the model macro-properties under different conditions.

4. The Results of Micro-Parameters on Mechanical Properties of Rock

4.1. The Influence of Micro-Parameters on UCS

4.1.1. Combined Effects of ${\overline{\sigma }}_n/{\overline{\tau }}_n$ and ${\overline{k}}^{\ast }$ on UCS

The relationships between UCS and ${\overline{\sigma }}_n/{\overline{\tau }}_n$ and ${\overline{k}}^{\ast }$ are shown in Figure 2. It could be seen from Figure 2a that when ${\overline{\sigma }}_n/{\overline{\tau }}_n$ < 0.8, the UCS of the model decreased continuously with increasing ${\overline{k}}^{\ast }$. The rate of the decrease declined with increasing ${\overline{k}}^{\ast }$, whereas the extent of the decrease grew with increasing ${\overline{\sigma }}_n/{\overline{\tau }}_n$. When 0.8 ≤ ${\overline{\sigma }}_n/{\overline{\tau }}_n$ < 4, the UCS initially increased and then decreased with increasing ${\overline{k}}^{\ast }$. The extent of the increase grew consistently with increasing ${\overline{k}}^{\ast }$, and the curve maximum continuously shifted toward the right. As shown in Figure 2b, when ${\overline{\sigma }}_n/{\overline{\tau }}_n$ ≥ 4, the UCS initially increased with increasing ${\overline{k}}^{\ast }$ and then leveled off. Additionally, as ${\overline{\sigma }}_n/{\overline{\tau }}_n$ increased, the change in UCS with ${\overline{k}}^{\ast }$ became consistently larger.

The results indicated that when ${\overline{k}}^{\ast }$ was constant, the UCS increased with increasing ${\overline{\sigma }}_n/{\overline{\tau }}_n$. When ${\overline{k}}^{\ast }$ was relatively low, there was a minor increase in the UCS with increasing ${\overline{\sigma }}_n/{\overline{\tau }}_n$. For ${\overline{k}}^{\ast }$ = 1.11, when ${\overline{\sigma }}_n/{\overline{\tau }}_n$ < 0.8, the UCS increased by 1.77-fold. When 0.8 ≤ ${\overline{\sigma }}_n/{\overline{\tau }}_n$ < 1, the increase in the UCS was 20.58%, and when ${\overline{\sigma }}_n/{\overline{\tau }}_n$ ≥ 1, there was an increase of 8.94%. With increasing ${\overline{k}}^{\ast }$, the increase in UCS with ${\overline{\sigma }}_n/{\overline{\tau }}_n$ became gradually larger. For ${\overline{k}}^{\ast }$ = 10.11, when ${\overline{\sigma }}_n/{\overline{\tau }}_n$ < 1, the UCS increased by 6.24-fold; when ${\overline{\sigma }}_n/{\overline{\tau }}_n$ ≥ 1, the increase in UCS was 3.22-fold.

The response surface plots of UCS are shown in Figure 3. When ${\overline{\sigma }}_n/{\overline{\tau }}_n$ < 0.8, the UCS exhibited a negative response to ${\overline{k}}^{\ast }$ that became more evident with increasing ${\overline{\sigma }}_n/{\overline{\tau }}_n$. When ${\overline{\sigma }}_n/{\overline{\tau }}_n$ ≥ 0.8, the response of UCS to ${\overline{k}}^{\ast }$ shifted from negative to positive. This positive response also became more evident with increasing ${\overline{\sigma }}_n/{\overline{\tau }}_n$, and UCS exhibited a positive response to ${\overline{\sigma }}_n/{\overline{\tau }}_n$. When ${\overline{\sigma }}_n/{\overline{\tau }}_n$ < 1, the positive response of UCS was attenuated with increasing ${\overline{k}}^{\ast }$; when ${\overline{\sigma }}_n/{\overline{\tau }}_n$ > 1, the positive response was enhanced by increased ${\overline{k}}^{\ast }$.

4.1.2. Combined Effects of KK and ${\overline{\sigma }}_n/{\overline{\tau }}_n$ on UCS

The relationships between UCS versus ${\overline{\sigma }}_n/{\overline{\tau }}_n$ and KK are shown in Figure 4. The UCS curve exhibited a consistent relationship when ${\overline{\sigma }}_n/{\overline{\tau }}_n$ < 1 and > 1. With increasing KK, the UCS decreased continuously at a gradually decreasing rate. When KK > 6, the UCS tended to level off. Additionally, with increasing ${\overline{\sigma }}_n/{\overline{\tau }}_n$, the overall decrease in the UCS became gradually larger.

4.1.3. Combined Effects of KK and ${\overline{k}}^{\ast }$ on UCS

The relationships between UCS versus ${\overline{k}}^{\ast }$ and KK at ${\overline{\sigma }}_n/{\overline{\tau }}_n$ = 0.8 and 4 are shown in Figure 5. The UCS consistently decreased with increasing KK. When ${\overline{k}}^{\ast }$ was relatively low, there was an approximately linear relationship between KK and UCS; the curve exhibited a very gradual decrease, indicating a minor effect of KK on UCS for relatively low ${\overline{k}}^{\ast }$. However, with increasing ${\overline{k}}^{\ast }$, there existed a nonlinear relationship between KK and UCS; the curve exhibited a larger decrease, and the effect of KK on the UCS became increasingly more evident.

4.2. The Influence of KK and ${\overline{k}}^{\ast }$ on Young’s Modulus and Poisson’s Ratio

The relationships between Young’s modulus versus ${\overline{k}}^{\ast }$ and KK at ${\overline{\sigma }}_n/{\overline{\tau }}_n$ = 0.8 and 4 are shown in Figure 6. When KK < 4, Young’s modulus decreased gradually with increasing KK, and there was a minor change in Young’s modulus when ${\overline{k}}^{\ast }$ was relatively low. However, the change in Young’s modulus became gradually larger with increasing ${\overline{k}}^{\ast }$ at a decreasing rate. When KK > 4, Young’s modulus generally remained constant, and KK had almost no effect. Additionally, as compared with ${\overline{k}}^{\ast }$, the change in KK had less effect on Young’s modulus, while the change in ${\overline{\sigma }}_n/{\overline{\tau }}_n$ had no effect.

The relationships between Poisson’s ratio versus ${\overline{k}}^{\ast }$ and KK at ${\overline{\sigma }}_n/{\overline{\tau }}_n$ = 0.8 and 4 are shown in Figure 7. When KK < 4, the Poisson’s ratio increased gradually with increasing KK. For relatively low ${\overline{k}}^{\ast }$, the Poisson’s ratio showed a minor change with KK. However, with increasing ${\overline{k}}^{\ast }$, the change in Poisson’s ratio became more evident, whereas the rate of increase diminished gradually. When KK > 4, the Poisson’s ratio was generally constant, and KK almost had no effect on the Poisson’s ratio. As compared with ${\overline{k}}^{\ast }$, the change in KK had less effect on the Poisson’s ratio, while the change in ${\overline{\sigma }}_n/{\overline{\tau }}_n$ had no effect.

4.3. The Effects of ${\overline{\sigma }}_n/{\overline{\tau }}_n$ and ${\overline{k}}^{\ast }$ on Stress–Strain Relationships and Failure Characteristics

The stress–strain relationships with different ${\overline{k}}^{\ast }$ at ${\overline{\sigma }}_n/{\overline{\tau }}_n$ = 0.8 and 4 are shown in Figure 8. For ${\overline{\sigma }}_n/{\overline{\tau }}_n$ = 0.8, when ${\overline{k}}^{\ast }$ < 3, the UCS and strain of specimens increased with increasing ${\overline{k}}^{\ast }$. A clear yield stage was not observed, and the curve rapidly dropped after reaching the peak value. There was a minor strain change in the post-peak stage. When ${\overline{k}}^{\ast }$ > 3, the sample UCS decreased with increasing ${\overline{k}}^{\ast }$, whereas little change occurred in the strain. Plastic flow increased considerably in the yield and post-peak stages, and sample failure gradually changed from brittle to ductile failure.

For ${\overline{\sigma }}_n/{\overline{\tau }}_n$ = 4, when ${\overline{k}}^{\ast }$ = 10.11, the samples showed no clear yield stage, and the stress dropped rapidly after reaching the peak value. With decreasing ${\overline{k}}^{\ast }$, the sample UCS and strain decreased continuously, whereas the strain in the yield stage increased gradually, and the plastic flow was evident. This trend was distinct from the case for ${\overline{\sigma }}_n/{\overline{\tau }}_n$ = 0.8.

As shown in Figure 9, PFC3D used the break of contact bond between particles to characterize the generation of micro-fractures. The propagation of micro-fractures in the test specimens with different ${\overline{k}}^{\ast }$ at ${\overline{\sigma }}_n/{\overline{\tau }}_n$ < 1 is illustrated in Figure 10. When ${\overline{\sigma }}_n/{\overline{\tau }}_n$ < 1, the formation of micro-fractures within the samples mainly resulted from tension. Under constant ${\overline{k}}^{\ast }$, the number of micro-fractures gradually decreased with increasing ${\overline{\sigma }}_n/{\overline{\tau }}_n$. When ${\overline{\sigma }}_n/{\overline{\tau }}_n$ = 0.8, the number of micro-fractures decreased from 44% to 63% compared to ${\overline{\sigma }}_n/{\overline{\tau }}_n$ = 0.1. When ${\overline{\sigma }}_n/{\overline{\tau }}_n$ was relatively low, the micro-fractures formed mainly by tensile failure. As ${\overline{\sigma }}_n/{\overline{\tau }}_n$ increased, the number of micro-fractures formed by shear failure gradually increased. Under constant ${\overline{\sigma }}_n/{\overline{\tau }}_n$, the number of micro-fractures within the samples increased with increasing ${\overline{k}}^{\ast }$, whereas the number of micro-fractures formed by shear decreased gradually. When ${\overline{k}}^{\ast }$ = 10.11, the number of micro-fractures increased from 28% to 96% compared to ${\overline{k}}^{\ast }$ = 1.11. Meanwhile, the bond fracture between particles gradually evolved from combined tension-shear to tension alone.

The propagation of micro-fractures within specimens with different ${\overline{k}}^{\ast }$ for ${\overline{\sigma }}_n/{\overline{\tau }}_n$ > 1 is illustrated in Figure 11. With increasing ${\overline{\sigma }}_n/{\overline{\tau }}_n$, shear gradually played a dominant role in the formation of micro-fractures within the samples. When ${\overline{k}}^{\ast }$ < 5, the number of micro-fractures gradually increased with increasing ${\overline{\sigma }}_n/{\overline{\tau }}_n$. When ${\overline{\sigma }}_n/{\overline{\tau }}_n$ = 6, the number of micro-fractures increased from 80% to 118% compared to ${\overline{\sigma }}_n/{\overline{\tau }}_n$ = 1. However, when ${\overline{k}}^{\ast }$ ≥ 5, the number of micro-fractures exhibited a reverse trend. When ${\overline{\sigma }}_n/{\overline{\tau }}_n$ = 6, the number of micro-fractures decreased from 0.08% to 55% compared to ${\overline{\sigma }}_n/{\overline{\tau }}_n$ = 1. With increasing ${\overline{k}}^{\ast }$, gradually more micro-fractures formed within the samples. The number of micro-fractures formed by shear decreased, whereas the number of micro-fractures formed by tension increased. When ${\overline{k}}^{\ast }$ = 10.11, the proportion of shear fractures in total micro-fractures decreased from 20 to 62%. Moreover, with increasing ${\overline{\sigma }}_n/{\overline{\tau }}_n$, the change in the micro-fracture number with increasing ${\overline{k}}^{\ast }$ reversed from an increase to a decrease.

5. Discussion

5.1. Effects of Different ${\overline{k}}^{\ast }$ Assignment Methods on the Macro-Properties

In previous studies, the assignment of model ${\overline{k}}^{\ast }$ is mostly implemented by using the term contact method pb_deformability. Using this method, the model parallel-bond stiffness is calculated as follows [37]:

$${\overline{k}}_n={\overline{E}}_c/L_R$$

(1)

$${\overline{k}}_s={\overline{k}}_n/{\overline{k}}^{\ast }$$

(2)

where ${\overline{k}}_n$ is bond normal stiffness, ${\overline{k}}_s$ is bond shear stiffness, and $L_R$ is the contact length. If the contact is particle-particle, then $L_R=R^1+R^2$, where $R^1$ and $R^2$ are the radii of the two particles, respectively; if the contact is particle-wall, then $L_R=R^1$.

According to the above formula, this assignment method actually achieves the target ${\overline{k}}^{\ast }$ by maintaining constant ${\overline{k}}_n$ and adjusting ${\overline{k}}_s$. Increasing ${\overline{k}}^{\ast }$ essentially leads to a continuous decrease of the overall bond stiffness of the specimens.

In this study, the term contact property was, therefore, used to assign values to ${\overline{k}}_n$ and ${\overline{k}}_s$. The target ${\overline{k}}^{\ast }$ was set by maintaining constant while varying ${\overline{k}}_n$ to evaluate the effects of different ${\overline{k}}^{\ast }$ assignment methods on the model macro-properties.

The effects of different ${\overline{k}}^{\ast }$ assignment methods on the macro-properties for ${\overline{\sigma }}_n/{\overline{\tau }}_n$ = 0.8 and 4 are shown in Figure 12. In case Ι, the assignment was implemented using the contact method pb_deformability; in case II, the assignment was implemented using contact property. In these two schemes, the UCS and Poisson’s ratio generally exhibited a consistent relationship with ${\overline{k}}^{\ast }$, whereas Young’s modulus showed distinct relationships with ${\overline{k}}^{\ast }$. In case Ι, Young’s modulus was inversely proportional to ${\overline{k}}^{\ast }$ and decreased continuously with increasing ${\overline{k}}^{\ast }$. In case II, Young’s modulus was positively proportional to ${\overline{k}}^{\ast }$ and increased approximately linearly with increasing ${\overline{k}}^{\ast }$.

The stress–strain relationships with different ${\overline{k}}^{\ast }$ for ${\overline{\sigma }}_n/{\overline{\tau }}_n$ = 0.8 and 4 in case II are shown in Figure 13. In the two schemes, the UCS changed with ${\overline{k}}^{\ast }$ following consistent trends; however, the trends of strain change were distinct between schemes. In case II, owing to the constant ${\overline{k}}_s$ and varying ${\overline{k}}_n$, the overall stiffness of specimens increased with increasing ${\overline{k}}^{\ast }$, while the deformability decreased continuously. Thus, the strain decreased constantly with increasing ${\overline{k}}^{\ast }$.

5.2. Mechanisms Underpinning the Effects of Micro-Parameters on Macro-Properties

The failure envelope for the PBM of PFC3D is shown in Figure 14. PBM determines the bond fracture between particles, mainly by the following method [37]. First, the tensile strength limit was determined. If the tensile strength limit was exceeded ($\overline{\sigma }>{\overline{\sigma }}_n$), the bond underwent tensile fracture. If the bond did not undergo tensile fracture, the shear strength limit was determined. The shear strength was expressed as:

$${\overline{\tau }}_c={\overline{\tau }}_n-\sigma \tan \overline{\phi }$$

(3)

$$\sigma ={\overline{F}}_n/\overline{A}$$

(4)

where ${\overline{F}}_n$ is the normal component of the parallel-bond contact force, and $\overline{A}=\pi R^2$ is the cross-sectional area of the contact bond between particles. If the shear strength limit was exceeded ($\overline{\tau }>{\overline{\tau }}_c$), then the bond underwent shear fracture.

Combined with

Table 4. Effect characteristics of bond normal-to-shear stress ratio (${\overline{\sigma }}_n/{\overline{\tau }}_n$) and parallel-bond stiffness ratio (${\overline{k}}^{\ast }$) on the macro-properties.

Type	Law of UCS Change	${\overline{\mathit{\sigma }}}_{\mathit{n}}/{\overline{\mathit{\tau }}}_{\mathit{n}}$	Curve Characteristics	Main Failure Characteristics of Samples
1		${\overline{\sigma }}_n/{\overline{\tau }}_n$ < 0.8	The UCS decreased with increasing ${\overline{k}}^{\ast }$
2		0.8≤${\overline{\sigma }}_n/{\overline{\tau }}_n$ < 4	The UCS initially increased and then decreased with increasing ${\overline{k}}^{\ast }$. The curve maximum continuously shifted toward the right
3		${\overline{\sigma }}_n/{\overline{\tau }}_n$ ≥ 4	The UCS initially increased with increasing ${\overline{k}}^{\ast }$ and then leveled off.

, the results showed that under the conditions of different ${\overline{\sigma }}_n/{\overline{\tau }}_n$, the effects of ${\overline{k}}^{\ast }$ on the macro-properties could be divided into three types. ${\overline{\sigma }}_n/{\overline{\tau }}_n$ < 0.8: In this case, the strength of specimens was mainly restricted by ${\overline{\sigma }}_n$, and the inter-particle bond failure was dominated by tensile fracture.

By selecting the contact force chain, which is an option in the PFC software, the stress state of the bond between particles could be observed. Additionally, Figure 15 shows a growing number of contact force chains under tension within the specimens. Owing to the relatively low ${\overline{\sigma }}_n$, the bond within the specimens first underwent tensile fracture because the external load exceeded ${\overline{\sigma }}_n$, and the sample strength decreased continuously. 0.8 ≤ ${\overline{\sigma }}_n/{\overline{\tau }}_n$ < 4: When ${\overline{k}}^{\ast }$ was relatively low, the sample strength was mainly restricted by ${\overline{\tau }}_n$ and the inter-particle bond failure was dominated by shear fracture. With increasing ${\overline{k}}^{\ast }$, the effect of ${\overline{\sigma }}_n$ on the sample strength increased. Accordingly, the sample strength changed from being restricted by ${\overline{\tau }}_n$ to being restricted by ${\overline{\sigma }}_n$, and the inter-particle bond failure was dominated by shear fracture. ${\overline{\sigma }}_n/{\overline{\tau }}_n$ ≥ 4: In this case, the sample strength was mainly restricted by ${\overline{\tau }}_n$, and the inter-particle bond failure was dominated by shear fracture. With increasing ${\overline{k}}^{\ast }$, there was a decreasing number of contact force chains under tension within the specimens, whereas the strength increased continuously.

As shown in Figure 2, the sample strength increased with ${\overline{\sigma }}_n/{\overline{\tau }}_n$. However, Figure 16 shows that for ${\overline{\sigma }}_n/{\overline{\tau }}_n$ > 10, the sample strength ceased to increase with increasing ${\overline{\sigma }}_n/{\overline{\tau }}_n$. Furthermore, Figure 16 shows that an increase in ${\overline{\sigma }}_n$ or ${\overline{\tau }}_n$ resulted in higher strength. In the research scheme of this paper, increasing ${\overline{\sigma }}_n/{\overline{\tau }}_n$ referred to maintaining a constant ${\overline{\tau }}_n$ while increasing ${\overline{\sigma }}_n$. When ${\overline{\sigma }}_n/{\overline{\tau }}_n$ reached a certain value (that we would call N), the ${\overline{\sigma }}_n$ was relatively high, and the bond did not undergo tensile fracture. The shear strength limit was then determined. Because ${\overline{\tau }}_n$ did not increase, the sample strength was restricted by shear strength and, thus, ceased to increase.

As demonstrated in Figure 14 and Figure 17, the value of n was related to ${\overline{k}}^{\ast }$ and the internal friction angle ($\overline{\varphi }$). The N decreased with increasing $\overline{\varphi }$ and increased with increasing ${\overline{k}}^{\ast }$ (Figure 18). This indicated that the range of the effect of ${\overline{\sigma }}_n/{\overline{\tau }}_n$ (0 < ${\overline{\sigma }}_n/{\overline{\tau }}_n$ < N) was affected by a combination of $\overline{\varphi }$ and ${\overline{k}}^{\ast }$.

According to Figure 19, the parallel-bond model degenerated into a linear contact model after bond failure [10,37]. In this case, the overall sample stiffness decreased rapidly, whereas the KK increased; that is, the $k^{\ast }$ of the linear contact portion increased. There was a larger decrease in the overall sample stiffness, while Young’s modulus decreased, and the Poisson’s ratio increased. The deformation parameter of the specimens deteriorated, and their strength declined.

6. Conclusions

Based on PFC3D simulation, this study analyzed the effects of different ${\overline{\sigma }}_n/{\overline{\tau }}_n$ and ${\overline{k}}^{\ast }$ on the macro-properties of the model using PFC3D simulations. Ten different bond normal-to-shear stress ratios (${\overline{\sigma }}_n/{\overline{\tau }}_n$ = 0.1, 0.2, 0.4, 0.6, 0.8, 1, 2 3, 4, and 6) were considered with 10 different bond normal-to-shear stiffness ratios (${\overline{k}}^{\ast }$ = 1.11, 2, 3, 4, 5, 6, 7, 8, 9, and 10.11) set for each stress ratio. The effects of different KK on the macro-properties of the model were also considered. The simulation results were analyzed, and the following results were obtained:

• Changes in ${\overline{k}}^{\ast }$ altered the bond stress state within specimens and their overall stiffness, thereby affecting their strength and stress–strain characteristics. Changes in bond stress state with ${\overline{k}}^{\ast }$ altered the extent of influence that ${\overline{\sigma }}_n/{\overline{\tau }}_n$ exerted on sample strength. The change of ${\overline{k}}^{\ast }$ assignment mode also had a significant effect on the deformation ability of the model.
• ${\overline{\sigma }}_n/{\overline{\tau }}_n$ altered the inter-particle bond strength. An increase in ${\overline{\sigma }}_n/{\overline{\tau }}_n$ improved the effect level of ${\overline{k}}^{\ast }$ on the sample strength. The sample strength increased with increasing ${\overline{\sigma }}_n/{\overline{\tau }}_n$. However, when ${\overline{\sigma }}_n/{\overline{\tau }}_n$ exceeded a certain value N, the sample strength was restricted by ${\overline{\tau }}_n$ and did not increase further. The N value was determined by a combination of $\overline{\varphi }$ and ${\overline{k}}^{\ast }$.
• An increase in KK resulted in a larger decrease in the overall stiffness of specimens after parallel-bond failure. The sample strength decreased as their deformability increased.
• Considerable changes in the sample brittleness and ductility were observed under the combined effects of ${\overline{\sigma }}_n/{\overline{\tau }}_n$ and ${\overline{k}}^{\ast }$. When ${\overline{\sigma }}_n/{\overline{\tau }}_n$ < 1, specimens changed from brittle to ductile with increasing ${\overline{k}}^{\ast }$, and when ${\overline{\sigma }}_n/{\overline{\tau }}_n$ > 1, they exhibited a reverse trend.
• According to the research results, the appropriate adjustment of the micro-parameters of PFC could make the mechanical properties of the simulation model more consistent with the actual mechanical properties of the rock material, improve the accuracy of the simulation results, and provide some help for the research of the rock mechanical behavior. Considering the variation range of micro-parameters in this paper, the results of this paper could only be used to simulate the mechanical properties of rock materials.