1. Introduction
With the rapid development of information technology and computational mechanics theory, more and more commercial software based on finite element method (FEM), discrete element method (DEM), boundary element method (BEM), or finite difference method (FDM) has been developed [
1,
2,
3,
4,
5,
6], such as Abaqus [
7], Ansys [
8,
9], FLAC (Fast Lagrangian Analysis of Continua) [
10,
11], Rocscience [
12], UDEC (Universal Distinct Element Code) [
13], and PFC (Particle Flow Code) [
14]. Since these pieces of software possess significant merits like low cost, high efficiency, and broad range of applications [
15,
16,
17], they have been widely applied to investigate both the macroscopic and microscopic mechanical behavior of large-scale rock engineering or small-size rock specimens.
DEM is a numerical solution devoted to describing the mechanical behavior of discontinuous media, particularly the rock and rock-like materials. Inspired by molecular dynamics method, Cundall [
18] first proposed a DEM model to simulate the gradual motion of blocky rock mass. Afterwards, it was employed to study the mechanical characteristics of granular assemblies, and the results indicated that the obtained force vector plots by a computer program show good agreement with those from photoelastic disc tests [
19]. Later on, a series of versions of DEM-based PFC software was released by Itasca Consulting Group, Inc. This piece of promising software has a wide range of applications, e.g., slope stability, rock fracture, underground mining, hydraulic fracturing, and flow simulation.
In PFC model, the medium is idealized as an assembly of various sized rigid particles and the bonds between them. The motion and interaction of the particles affect the macromechanical behavior. The shape of the particles is circular in PFC2D (Two-dimensional PFC) and spherical in PFC3D (Three-dimensional PFC). There are two types of bonds: contact bond and parallel bond [
20] (
Figure 1). For the contact-bond model, the adjacent particles are connected by point-formed bonds; that is, the particles can only move in translation and not in rotation. As long as the particles are in contact, the contact stiffness remains the same regardless of whether the bond is broken [
21]. Clearly, this is not realistic in rock mechanics. With regard to the parallel-bond model, particles are contacted via a bond with a certain thickness and radius. When particles translate or rotate, both the force and moment can be transmitted through the bond. Once the normal or tangential stress on the bond exceeds its normal or tangential strength, fracture will occur as tensile or shear cracks. This behavior is consistent with the actual fracture of the rock. Hence the parallel-bond model was selected for PFC modeling in this study. During modeling, the macroscopic physical and mechanical properties as well as the constitutive model of the material are not required to predefine. The reason is that the mechanical response is controlled by the motion and interaction of the particles. Thus, determining a set of micro-scaled physical and mechanical parameters of the particle and contact to match the macro-response of the material is essential. Obviously, the accuracy of PFC simulation is strongly associated with these parameters.
Currently, the trial-and-error method is commonly used to calibrate the micro-scale parameters of the PFC model based on uniaxial compression test results of rock specimens; that is, they need to repeat performing uniaxial compression simulation by adjusting these parameters until the simulated uniaxial compressive strength, elastic modulus, Poisson’s ratio, and failure mode of the model are consistent with the laboratory test results. Indeed, the number of the micro-scale parameters is quite large, and determining a proper set of micro-scale parameter values is extremely time consuming and troublesome. Therefore, it is meaningful to find an approach to address this issue quickly and reliably. To explore the relationship between micro-scale parameters and macroscopic mechanical properties, a few attempts have been made by some researchers. Yoon [
22] applied the Plackett–Burman design method to analyze the sensitivities of micro-scale parameters to macro-response of contact-bond model, and singled out two most significant micro-scale parameters for each macroscopic property. By using Taguchi method, Hanley et al. [
23] performed a calibration on micro-structural parameters of parallel-bonded agglomerates in both 2D and 3D models. It is found that the normal strength of bond, bond stiffness, and particle friction are the dominant parameters influencing the strain at failure. Cheung et al. [
24] conducted a parametric study to identify the effect of parallel bond-contact stiffness ratio, bond area, and bond strength on stress–strain response of PFC3D model. Results show that the overall stiffness increases with the increase of the first two parameters, but is not affected by the bond strength. Shi et al. [
25] summarized that the macroscopic Young’s modulus is controlled by the modulus of the bond and particle, the Poisson’s ratio depends on the ratio of particle normal to shear stiffness, and the uniaxial compressive strength is related to the bond shear strength. Also, the linear equations between them are fitted. Castro-Filgueira et al. [
26] argued that bond cohesion has a remarkable impact on the model strength. Ding et al. [
27] further examined the effects of maximum to minimum particle size ratio and particle size distribution on the macroscopic mechanical properties of PFC3D, which show good agreement with the results of Wang et al. [
28] using similitude rules. Besides, based on numerous uniaxial compression simulations, Tawadrous et al. [
29] utilized artificial neural networks to predict the micro-properties required for PFC3D modeling under three given macro-properties. By means of gray relational analysis and factorial screening, Zou and Lin [
30] obtained the dominant micro- parameters influencing the uniaxial compressive strength, Young’s modulus, and Poisson’s ratio of coal. Moreover, Wang and Cao [
31] used an improved simulated annealing algorithm to find the approximate values of these micro-parameters, which does not require establishing the relationship between microscopic and macroscopic parameters. Sun et al. [
32] further developed FFD-ANN (Full Factorial Designs and Artificial Neural Networks) method to predict the micro-parameters of PFC3D model for triaxial compression tests.
From the literature review, it can be concluded that relationship between partial micro- and macro- parameters has been revealed. However, so far, it is still unclear that how these parameters interact with each other and collectively affect each macro-scaled mechanical property. Moreover, information regarding effective methods for finding the micro-scale parameters quickly is relatively rare. Therefore, in this research, a novel experimental design method incorporating Box–Behnken response surface methodology and desirability function optimization (BB-DFO) approach was proposed to conduct sensitivity analysis and optimization of micro-scale parameters of parallel-bond model. After that, the fracture development around circular, square, and D-shaped openings in rock specimens loaded in uniaxial compression was simulated using the calibrated PFC model and then was compared with the results of experimental investigation.
4. Experimental Design and Analysis
4.1. Experimental Scheme and Results
As stated above, three influencing factors and three response variables were considered in this study. Hence, a three-level-three-factor model was established using Box–Behnken design. According to the laboratory test results, the defined levels of each factor are given in
Table 2. To make statistical calculation simple, dimensionless processing was conducted on these factors according to
where
Xi and
xi (
i = 1, 2, and 3) denote the dimensionless coded value and actual value of influencing factor, respectively,
x0 is the actual value of the factor at the center point, and Δ
xi means the step change value of the influencing factor.
According to the Box–Behnken design scheme, a total of 17 groups of uniaxial compression simulations were carried out by the PFC software.
Table 3 shows the simulation results. Accordingly, different models, including the linear, 2FI (Two-Factor Interaction), quadratic, and cubic models were employed to fit the second-order response surface equation. Lack-of-fit test results indicate that the quadratic model fits well the nonlinear relationship between the independent and dependent variables (see
Table 4). It is noted that the cubic model and higher are aliased. The fitted equations of the response variables with respect to the actual factors are given as
Equations (11)−(13) show that all the correlation coefficients are greater than 0.99, suggesting that the fitted response surface functions are reliable. Also, the predicted values of the response variables are found to be extremely agreeable with the simulated results (
Figure 3). The average prediction errors of the uniaxial compressive strength, Young’s modulus, and Poisson’s ratio reach 0.70%, 0.85%, and 0.21%, respectively.
4.2. Analysis of Variance
To further identify the significance of A, B, and C as well as their interactions on y1, y2, and y3, analysis of variance was performed on each response variable model. In general, the significance of the model terms can be evaluated on the basis of F-value and p-value. It is calculated that the F-values of the three models are 1480.09, 1372.12, and 3553.51, respectively, all of which are greater than F0.95(8,8) = 3.438. For the p-values, they are all less than 0.0001. This implies that these fitted models are extremely significant. Note that the larger the F-value or the smaller the p-value is, the more significant the model will be. With regard to the probability of the model term, p < 0.0001, p < 0.05, and p > 0.05 means it is extremely significant, significant, and insignificant, respectively. Also, the signal-to-noise ratios of the three models are 113.49, 131.69, and 163.63, respectively, which are much larger than 4. Clearly, these models with adequate signals can be applied to navigate the design space. Moreover, their coefficients of variation are 1.23%, 1.38%, and 0.42%, respectively, indicating that the models have high accuracy.
4.2.1. Influence of Micro-Scale Parameters on Uniaxial Compressive Strength
Table 5 shows that the effects of
A and
C on uniaxial compressive strength are extremely significant, while that of
B is not significant. As the value of
A increases, the uniaxial compressive strength increases linearly and significantly. By contrast, the uniaxial compressive strength increases slowly with the increase of
C, but the factor
B has little effect on the uniaxial compressive strength (see
Figure 4a). When the coded values of
B and
C are both zero, the uniaxial compressive strength (192 MPa) when
A is 1 is found to be 2.87 times that (67 MPa) when it is −1. Within the range of
B, the strength changes from 120 to 130 MPa. It can be observed that the effect of
A on uniaxial compressive strength is more significant than that of
C.
According to the
p-value obtained from the analysis of variance, it is further found that the interactions between
A and
C and between
B and
C have a significant effect on the uniaxial compressive strength, which are illustrated in
Figure 4b,c.
Figure 4b shows that the uniaxial compressive strength increases as
A and
C increase simultaneously. The response surface plot also verifies that
A has a greater influence on intensity than
C. Further on, the strength is also significantly affected by the interaction between
B and
C (see
Figure 4c). When the coded values of
B and
C are near zero, the uniaxial compressive strength has a maximum value.
4.2.2. Influence of Micro-Scale Parameters on Young’s Modulus
Similarly, we can obtain the impact of these three micro-scale parameters on the macroscopic Young’s modulus.
Figure 5 clearly demonstrates the one-factor and multiple factor interaction effects. The influence of
B and
C on Young’s modulus are extremely significant, while
A does not affect it significantly. The Young’s modulus linearly changes from 22 to 70 GPa as the coded value of
B increases from −1 to 1. When
C varies from −1 to 1, the Young’s modulus decreases nonlinearly from 59 to 40 MPa, the reduction rate gradually decreases. Besides, it is found that the interaction between
B and
C has an extremely significant effect on the Young’s modulus, whereas the interactions of
A–
B and
A–
C are not significant. With the rise of
B and the decrease of
C, the Young’s modulus gradually increases.
4.2.3. Influence of Micro-Scale Parameters on Poisson’s Ratio
The effects of
A,
B, and
C and their interactions on the Poisson’s ratio can also be derived based on the results of the analysis of variable.
Figure 6a shows that the influence of
B and
C on the Poisson’s ratio are significant and extremely significant, respectively, while that of
A is not remarkable. The Poisson’s ratio increases slowly with the increase of
B, but grows first sharply and then slowly as
C rises. Within the range of
C, the minimum and maximum Poisson’s ratios are 0.16 and 0.29, respectively, with an increase rate of 81.25%. By contrast, with the increase of
B, the increase of Young’s modulus is very small, and basically remains unchanged at 0.25. For the effect of the interaction, it is observed that there are no significant interactions between
A and
B and between
A and
C. The
p-value of the
B–
C interaction is 0.3506, implying that it is also not significant.
Figure 6b displays that the contours of the Poisson’s ratio are approximately horizontal.
4.3. Optimization of Micro-Scale Parameters
According to the principle of the desirability function method, first, we need to determine the optimization interval for each factor. The purpose of this study is to find the micro-scale parameters to fit the macroscopic mechanical response. Therefore, Equation (8) was utilized to calculate the individual desirability, and then the overall desirability could be obtained according to Equation (9). Assuming these three factors are equally important, their weights were all set to 5 (
r1 =
r2 =
r3 = 5). The optimization intervals of
A,
B, and
C were defined as (50,150), (20,60) and (1,3), respectively, and their target values were set to be equal to those of the uniaxial compression test results, i.e.,
y1 = 102.61 MPa,
y2 = 20.78 GPa, and
y3 = 0.258. The effect of each factor on the overall desirability is plotted in
Figure 7. It is shown that the desirability is the maximum when
A is 75~85 MPa,
B is 20~23 GPa, and
C is 1.6~2.7. In these intervals, the desirability increases first and then decrease with the increasing
A and
C, while it declines as
B increases. Moreover,
Figure 8 illustrates that the interactions between these factors have a significant effect on the overall desirability. It can be observed clearly that the overall desirability is the largest under the condition of
A is about 82 MPa,
B is about 20 GPa, and
C is about 2.2.
Through optimization, we found that the overall desirability of the solution that A = 82.07 MPa, B = 20.00 GPa, and C = 2.198 is the maximum, namely, D = 0.914. In such a case, the simulated uniaxial compressive strength, Young’s modulus, and Poisson’s ratio are 102.61 MPa, 21.77 GPa, and 0.258, respectively, which are exceedingly agreeable with the laboratory results.
In addition to the simulated responses being consistent with laboratory test results, the failure mode of the model should also coincide with that of the rock specimen. Thus, minor adjustments are required on these parameters and the standard deviations of normal and shear strength (
and
). The final micro-scale parameters obtained are shown in
Table 6, and the comparison between experimental and numerical stress–strain curves is presented in
Figure 9. The peak stresses and slopes of the two curves are basically equal, with relative errors of 1.19% and 5.92%, respectively. This indicates that the uniaxial compressive strength and Young’s modulus are consistent. Since the built models compressed by two discs are strictly two dimensional, there is no initial pores and cracks compaction stage during the loading process. Consequently, the peak stain of the model is slightly smaller than that of the real rock specimen.