Rock Macro–Meso Parameter Calibration and Optimization Based on Improved BP Algorithm and Response Surface Method in PFC 3D

Abstract

In order to obtain the calibration law of rock macro and meso parameters under three-dimensional conditions, based on the parallel bond model, starting with the basic theory of PFC and the qualitative relationship between macro and meso parameters, an orthogonal experimental scheme is designed. An improved BP algorithm is proposed, which has a function with gradient factor, adaptive Nesterov momentum method, and adaptive learning rate for the lightweight analysis of meso parameters. The sensitivity between macro and meso parameters is quantified, and the key meso parameters are screened out. Based on the lightweight model, the correlation and influence mechanisms between macro and meso parameters are analyzed. It was found that the elastic modulus increases linearly with the increase in equivalent modulus. The parallel bond stiffness ratio can inhibit the growth of the elastic modulus, and the elastic modulus decreases greatly when the stiffness is relatively high. There is a linear relationship between Poisson’s ratio and stiffness ratio, and the increase in the friction coefficient will cause the nonlinear decrease in it. The strength parameters have an incentive effect on the peak strength. When the tensile and shear strengths of the parallel bond are at a high level, the combination has the most significant effect on the increase in the peak strength. The internal friction angle mainly has a certain influence on the postpeak strength of the rock, because it has a control effect on the particle sliding on both sides of the shear zone when the sample is loaded after the peak. Based on the central composite experimental design and response surface method, a nonlinear model of macro–meso parameters described by multiple subresponse surfaces is obtained. Finally, the mathematical model of parameter calibration is established, and the optimal solution of rock meso parameters is obtained by using optimization techniques. Through the example verification, it was found that the numerical experiment and laboratory test results are close in the stress characteristics, stress evolution, and failure mode of the sample, which proves the effectiveness and reliability of the calibration method. The research results have a certain reference value for PFC parameter calibration.

1. Introduction

Rock is a complex mineral aggregate material. Under the long geological process, the rock generally contains primary structural planes, such as joints and fissures, and even faults pass through. Therefore, it shows complex nonlinear characteristics, such as discontinuity, heterogeneity, and anisotropy. The mechanical properties and damage evolution mechanism of rock have important theoretical significance and engineering value for the research of damage, disaster prevention, and reduction of geotechnical engineering.

In 1971, Cundall [1] introduced the principle of molecular dynamics; then, he proposed the particle flow theory and developed the particle flow (PFC) program. Since then, he joined with Strack [2] and established the granular discrete element calculation theory suitable for geotechnical engineering. In 1991, Wang et al. [3] introduced the discrete element theory to China and used it to solve the rock mechanics problems of underground engineering. In recent years, with the continuous development of the particle flow theory, it has become one of the most important numerical methods for researching geotechnical problems.

Choosing reasonable and reliable input parameters is the premise of accurate numerical simulation. Therefore, the calibration of PFC parameters has always been a challenging aspect of geotechnical engineering numerical simulation. The trial-and-error method generally has no specific rules and depends on the experience of researchers. At the same time, it has encountered several limitations in practical application, including low efficiency, high calculation costs, and difficulty using the gradient algorithm. Therefore, scholars around the world have explored this issue extensively.

The correlation and intrinsic connection between macro and meso parameters are the basis of PFC calibration research [4]. Potyondy et al. [5] systematically expounded the principle of the granular discrete element method and denoted that the geometric parameters of particles not only control the accuracy of the model but also have a certain impact on the mechanical properties and damage characteristics of materials. Zhao et al. and Chen et al. [6,7] applied theoretical and numerical methods to study the correlation between macro and meso parameters of different contact models and established a semi-quantitative function model between two scale parameters based on PFC2D. The calibration efficiency of the trial-and-error method improved. Cho et al., Zhang et al., and Shi et al. [8,9,10] researched the influence of PFC cluster particle geometry and crack distribution on the macroscopic mechanical behavior of rock. Cong et al. and Abi et al. [11,12] explored the correlation between macro- and mesoscopic parameters of geotechnical materials under plane conditions comprehensively and considered the influence of meso parameters on the shear mechanical behavior of rock.

In respect of parameter calibration methods, there are currently methods based on statistics, optimization techniques, neural networks, data mining, machine learning, and other technologies. Yang et al and Cui et al. [13,14] proposed a calibration method for gravel materials. Yoon et al. and Deng et al. [15,16,17] obtained the quantitative relationship between macro and meso parameters of rock by combining experimental design and optimization technology. This method has a good calibration effect on various rocks. For unsaturated soil, De pue et al. [18] simulated a great quantity of the uniaxial compression test and calibrated the meso parameters by the Kriging method in PFC2D. Zou et al. [19] took coal as the research object and proposed a combined optimization calibration method based on the Grey–Taguchi method, response surface method, and Mahalanobis distance measurement method. Xu et al. [20,21] introduced the iterative idea into parameter calibration, which significantly improved the efficiency and accuracy of calibration. In addition, Zhou et al. [22] used data mining and neural network to preliminarily predict the meso parameters of materials.

To sum up, the response law and calibration method between macro and meso parameters of geotechnical materials are important research directions in the particle discrete element method. At present, there have been some research achievements in this field using theoretical analysis, laboratory tests, and numerical simulation. However, the following aspects are still worthy of further study and discussion: (1) Most of the existing studies tend to explore the response law between macro and meso parameters through experimental design and sensitivity analysis to obtain the empirical relationship or fitting equations between a few parameters. It is difficult to obtain a complete set of meso parameters for numerical calculations. (2) Most of the research results are based on PFC2D and biaxial tests. However, rocks in practical engineering are often in a three-dimensional complex stress state, and not all rock materials can be simplified as plane problems. Therefore, this method has certain limitations. (3) The current calibration method is still in the semiempirical stage, which can only calibrate some parameters. There are still a few parameters that need to be selected according to experience or the literature, which makes it difficult to fully reflect the real physical and mechanical behaviors of materials. The theories and methods need to be further improved and explored.

Therefore, based on the research results of many scholars, this paper aims at the calibration of the meso parameters in the PFC parallel bond model. The research mainly contains the following three parts: (1) screening out the key meso parameters by the improved BP algorithm; (2) analyzing the influence mechanism and mathematical model between macro and meso parameters using the response surface method; (3) calibration and solution through the optimization technique based on the achievements in (1) and (2). First, some research samples were obtained through orthogonal design and triaxial numerical tests. Then, based on the improved BP algorithm, lightweight processing was performed to accurately quantify the influence of each meso parameter on macro parameters, and the key meso parameters that have a significant impact on macro parameters were screened out. Next, according to test results, the linear model between macro and meso parameters was regressed. The interaction law and influence mechanism between them were analyzed. Then, the nonlinear mathematical model between the key meso and macro parameters was obtained by using the central composite experimental design and response surface method. Finally, based on previous research results, the mathematical model for parameter calibration was established. The meso parameters of the parallel bond model were obtained by solving the calibration model through optimization techniques. The calibration results were compared with the laboratory test results, which verifies the rationality of the method. This method was based on the basic theory of PFC and solved by machine learning and optimization techniques, which can obtain the meso parameters in PFC3D reasonably and accurately, overcoming the defects of traditional methods, such as low efficiency and incomplete parameter calibration. We hoped to provide some reference for parameter matching and a solution for PFC numerical calculation.

2. PFC Contact Model Selection

PFC provides three contact models [23] to describe the contact constitutive relationships of materials, including the linear contact model, bond model, and Hertz–Mindlin contact model. The bond model includes the contact bond model and parallel bond model. Rock is a complex mixture composed of mineral particles and cement. The parallel bond model is usually used to characterize its contact behavior. The model is divided into linear and bond parts. The linear part is used to describe the contact behavior of mineral particles, and the bond part describes cement. The contact schematic and contact model are shown in Figure 1.

It can be seen from the figure that the linear and the parallel bond elements of the model act in parallel and establish an elastic interaction between two contact parts. The model has axial, tangential, and rotational stiffness. In addition, the damping term is introduced to complete the energy dissipation process of the model.

The linear contact force is generated by two linear springs with constant normal and tangential stiffness in the linear element, which can only resist the external force. The parallel bond contact force and moment are generated by two springs with constant normal and tangential stiffness in the parallel bond element. Due to the bonding effect, they can resist the relative rotation between components. The damping force is generated by the damper, and the viscosity of the damper is related to the critical damping ratio.

The force and failure mechanism of the particles in the parallel bond model is shown in Figure 2. It can be seen that for the linear part when the shear contact force between the two particles exceeds the allowable maximum, the particles will slip. The strength criterion of the bonding part can be described according to the beam bending theory. When the maximum stress between particles exceeds its allowable stress, the bonding will be broken. A tensile or shear crack will occur.

3. Meso Parameters and Experimental Design

3.1. Qualitative Relationship Description of Macro–Meso Parameters

There are many physical and mechanical parameters of rock, and its macro mechanical properties are mainly characterized by deformation and strength parameters. In this paper, the deformation parameters elastic modulus and Poisson’s ratio (E, $\nu$) and the strength parameter peak strength (${\sigma }_{\mathrm{c}}$) were selected to study. The meso parameters of rock mainly include three categories: geometric parameters (L, W, L/R, $R_{\max }/R_{\min }$, and n), material and deformation parameters ($\rho$, $E^{*}$, ${\overline{E}}^{*}$, $k_{\mathrm{n}}$, $k_{\mathrm{s}}$, $\overline{d}$, and $\mu$), and particle strength parameters (${\overline{\sigma }}_{\mathrm{c}}$, ${\overline{\tau }}_{\mathrm{c}}$, and $\phi$).

According to the research of Potyondy et al. [5], the qualitative relationship between macro and meso parameters can be expressed as follows:

$$Y_i={\Phi }_i\left[\left(L,W,n,\frac{L}{R},\frac{R_{\max }}{R_{\min }}\right),\left({\overline{\sigma }}_c,{\overline{\tau }}_c,\phi \right),\left(\rho ,E^{*},{\overline{E}}^{*},\frac{k_n}{k_s},\frac{{\overline{k}}_n}{{\overline{k}}_s},\overline{d},\mu \right)\right]$$

(1)

where L and W are the height and width of model. n is the particle porosity. L/R is the ratio of height to average particle size. $R_{\max }/R_{\min }$ are the maximum and minimum particle size ratio. $\rho$ is the particle density. $E^{*}$ and ${\overline{E}}^{*}$ are the equivalent modulus and parallel bond equivalent modulus, respectively. $k_n/k_s$ is the particle stiffness ratio. ${\overline{k}}_n/{\overline{k}}_s$ is the parallel bond stiffness ratio. $\overline{d}$ is the control gap of the parallel bond. $\mu$ is the friction coefficient. ${\overline{\sigma }}_{\mathrm{c}}$ and ${\overline{\tau }}_{\mathrm{c}}$ are the tensile and shear strengths of the parallel bond, respectively. $\phi$ is the internal friction angle.

3.2. Scheme of Orthogonal Test

This paper mainly studies the calibration of the meso parameters in the parallel bond model. Therefore, the geometric size, particle density, and porosity of the model are fixed firstly. According to the provisions of the rock test code for water resources and hydropower for sample size, the cylindrical sample with a bottom radius of 50 mm and a height of 100 mm was used in the numerical test. The research of Nohut et al. [24] shows that the particle density has little effect on the calibration, so, the density is 2700 kg/m³ with reference to the macro density of rock. The porosity is recommended to be 0.16 according to the literature [6,12,25].

The research scope of the remaining parameters was appropriately expanded based on the literature [11,24,25] to expand the applicability and reliability of the results. At the same time, to reduce variables, this paper assumed that $E^{*}={\overline{E}}^{*}$. The research scope of each meso parameter finally was determined as follows: L/R were set to 50–250; $R_{\max }/R_{\min }$ were set to 1–5; $E^{*}$(${\overline{E}}^{*}$) was set to 10–150 GPa; $k_n/k_s$ and ${\overline{k}}_n/{\overline{k}}_s$ were set to 0.8–4.0; $\overline{d}$ was set to 0–1; ${\overline{\sigma }}_{\mathrm{c}}$ and ${\overline{\tau }}_{\mathrm{c}}$ were set to 20–200 MPa; $\phi$ was set to 40–60°, and $\mu$ was set to 0.1–1.5.

The traditional test scheme requires a large number of tests and calculation resources because of the large number of parameters involved in the research and the large cross and uncertainty. The orthogonal test method can comprehensively consider the cross effects among various factors and effectively reduce the number of tests on the premise of reliable experimental results. According to the principle of orthogonal design, a L₅₀5¹⁰ (50 trials in total with 10 factors and 5 levels) orthogonal table was designed to carry out the triaxial compression numerical test in this paper. The specific scheme is shown in

Table 1. Scheme of orthogonal test.

No.	L/R	${\mathit{R}}_{\mathbf{max}}/{\mathit{R}}_{\mathbf{min}}$	${\mathit{E}}^{*}({\overline{\mathit{E}}}^{*})/\mathbf{GPa}$	${\mathit{k}}_{\mathit{n}}/{\mathit{k}}_{\mathit{s}}$	${\overline{\mathit{k}}}_{\mathit{n}}/{\overline{\mathit{k}}}_{\mathit{s}}$	$\overline{\mathit{d}}$	${\overline{\mathit{\sigma }}}_{\mathbf{c}}/\mathbf{MPa}$	${\overline{\mathit{\tau }}}_{\mathbf{c}}/\mathbf{MPa}$	$\mathit{\phi }$/°	$\mathit{\mu }$
1	50	1	10	0.8	0.8	0	20	20	40	0.1
2	50	2	10	1	1	0.2	50	20	40	0.1
3	50	3	10	2	2	0.5	100	20	40	0.1
4	50	4	10	3	3	0.8	180	20	40	0.1
5	50	5	10	4	4	1	200	20	40	0.1
…	…	…	…	…	…	…	…	…	…	…
47	250	2	100	0.8	2	0.8	200	20	45	0.5
48	250	3	100	1	3	1	20	20	45	0.5
49	250	4	100	2	4	0	50	20	45	0.5
50	250	5	100	3	0.8	0.2	100	20	45	0.5

Notes: The five levels of L/R are 50, 100, 150, 200, and 250, respectively. The values of $R_{\max }/R_{\min }$ are 1, 2, 3, 4, and 5, respectively. The values of $E^{*}$(${\overline{E}}^{*}$) are 10, 20, 50, 100, and 150, respectively. The values of $k_n/k_s$ and ${\overline{k}}_n/{\overline{k}}_s$ are 0.8, 1, 2, 3, and 4, respectively. The values of $\overline{d}$ are 0, 0.2, 0.5, 0.8, and 1.0, respectively. The values of ${\overline{\sigma }}_{\mathrm{c}}$ and ${\overline{\tau }}_{\mathrm{c}}$ are 20, 50, 100, 180, and 200, respectively. The values of $\phi$ are 40, 45, 50, 55, and 60, respectively. The values of $\mu$ are 0.1, 0.3, 0.5, 1.0, and 1.5, respectively. Due to the limited space, tests 6–46 were omitted.

.

In order to eliminate the adverse effect on the calibration because of the magnitude difference between the meso parameters and improve the accuracy and reliability of calibration, the meso parameters were standardized before calibration, and this can be transformed as follows:

$$y_i=\frac{x_i-\overline{x}}{s}$$

(2)

where $x_i$ and $y_i$ are the sequence of parameters before and after the transformation, respectively. $\overline{x}$ and s are mean and standard deviation, respectively.

After transformation, the new parameter sequence obeys the standard normal distribution $y_i\sim N\left(0,1\right)$, which is convenient for analysis and comparison among parameters.

3.3. Numerical Test and Servo Loading

In this paper, PFC 3D was used to simulate the triaxial compression test of rock. The model was a cylindrical sample with a bottom radius of 50 mm and a height of 100 mm, which is shown in Figure 3. In the numerical test, the basic meso parameters were selected according to

Table 2. Basic meso parameters in model.

$\mathit{\rho }$/(kg/m3) 1	n/% 2	${\mathit{k}}_{\mathit{n}}/\mathbf{GPa}{ }^3$	c 4	${\mathit{k}}_{\mathit{n}}{}^{\mathit{w}}/\mathbf{GPa}{ }^5$	${\overline{\mathit{k}}}_{\mathit{n}}/\mathbf{GPa}{ }^6$	$\overline{\mathit{\lambda }}{ }^7$
2700	16	2.0	0.7	0.3	50	1.0

¹ $\rho$ is particle density. ² n is porosity. ³ $k_n$ is normal stiffness of particles. ⁴ c is particle local damping coefficient. ⁵ $k_n{}^w$ is normal stiffness of wall. ⁶ ${\overline{k}}_n$ is normal stiffness of parallel bond. ⁷ $\overline{\lambda }$ is parallel bond radius multiplier.

, and the confining pressure was set to 10 MPa.

During the triaxial test in PFC, the loading of the model was achieved by applying a fixed velocity to the wall. The stress state of each particle was updated every moment. If the static loading mode is used, it will be difficult to maintain the constant stress in three directions due to the interaction between particles and walls. In order to provide a stable confining pressure, it was necessary to servo control the wall [26]. The velocity of the boundary wall during the servo process can be expressed as

$${\dot{u}}^w=G({\sigma }^{mea}-{\sigma }^{req})=G\Delta {\sigma }^w$$

(3)

where G is the servo coefficient. ${\sigma }^{mea}$ and ${\sigma }^{req}$ are current stress and objective stress, respectively.

$$G=\frac{\alpha A}{n{k}^w\Delta t}$$

(4)

where n is number of contacts. $k^w$ is the average contact stiffness of walls. $\alpha$ is the control coefficient to ensure iterative stability and takes a value between 0 and 1.

4. Lightweight Modeling Method Based on Improved BP Algorithm

The research results of most scholars show that there is a close relationship between macro and meso parameters. The meso parameters not only affect the macro parameters independently but also have a cross and correlation between the parameters, which makes the analysis and calibration in PFC a typical high-dimensional problem. If all parameters are taken into account without treatment, although the calibration accuracy is increased to a certain extent, the calibration difficulty and cost will be greatly increased and will even affect the interpretability of the model. Therefore, lightweight processing of the model, which can reasonably screen out key information and eliminate redundant data, is very important for accurate and efficient parameter calibration.

The BP neural network is a favorable tool to deal with nonlinear problems, and it is also one of the most popular intelligent algorithms in machine learning. Through its strong learning, parallelism, and fault tolerance capabilities, it can accurately and efficiently realize nonlinear mapping from the input to the output layers. It is widely used in nonlinear systems, function approximation, system identification, etc. It can solve problems in geotechnical engineering fields, such as parameter calibration, prediction, inversion, modeling, and other complex engineering technology problems.

In this paper, an improved BP algorithm is proposed, which can flexibly change the form of function and adjust the momentum coefficient and learning rate adaptively during the calculation process. It makes up for the shortcomings of traditional algorithms, such as poor stability and slow convergence through training the orthogonal test results. The degree of interaction between the macro and meso parameters can be sorted. Then, key meso parameters can be screened to achieve the purpose of lightweight.

4.1. Basic Principle of BP Algorithm

Unlike the statistical model, the neural network does not require test samples to satisfy the conditions of normal distribution and homogeneity of variance, and the number of samples tends to infinity. Through its strong learning, parallelism, and fault tolerance, it can accurately and efficiently realize the nonlinear mapping from the input to the output layers [27,28].

The multilayer perceptron neural network based on the BP algorithm is a typical feedforward neural network. It is generally composed of three parts shown in Figure 4. A set of input layers consisting of source nodes and several hidden layers and output layers consisting of computing nodes. In the figure, x is the input; w is the weight, and u is the output. The basic principle of the BP algorithm can be summarized as two processes: the forward propagation of the input signal and the back propagation of the error. The model error is minimized through repeated alternating iterations.

The sigmoid function is used as the function of neurons to establish the mapping between input (meso parameters) and output (macro parameters). It can be expressed as:

$$f(x)=\frac{1}{1+e^{-x}}$$

(5)

When the actual output of the network is not equal to the objective output (macro parameters measured in laboratory test), there is an error, and the error function is defined as:

$$J_p=\frac{1}{2}{\displaystyle \sum _{j=1}^{N_L}{\left(u_{L,j}-d_j\right)}^2}$$

(6)

$$J={\displaystyle \sum _{p=1}^P{J}_p}$$

(7)

where $u_{L,j}$ is the output of node j in layer L. $d_j$ is the objective output. $J_p$ is the error function of sample p, and J is total error of the system. L and P are the number of network layers and samples, respectively.

Taking the error function as the objective function, the weight of each layer is iterated and corrected using the gradient descent method, giving,

$$\Delta w_{l,i,j}(n)=-\eta \frac{\partial J}{\partial w_{l,i,j}}$$

(8)

where $\eta$ is the learning rate, and $\eta \in (0,1)$. $w_{l,i,j}$ is the weight between node i in layer l-1 and node j in layer l.

For the output layer, its weight can be corrected according to Equation (9).

$$\Delta w_{L,i,j}(n)=-\eta {\displaystyle \sum _{p=1}^P\frac{\partial J_p}{\partial w_{L,i,j}}}=-\eta {\displaystyle \sum _{p=1}^P\frac{\partial J_p}{\partial u_{L,j}}}·\frac{\partial u_{L,j}}{\partial w_{L,i,j}}=-\eta {\displaystyle \sum _{p=1}^P{\left[(u_{L,j}-d_j)u_{L,j}(1-u_{L,j})u_{L-1,i}\right]}_p}$$

(9)

For the hidden layer, there is a recurrence formula,

$$\Delta w_{l,i,j}(n)=-\eta {\displaystyle \sum _{p=1}^P\frac{\partial J_p}{\partial u_{l,j}}}·\frac{\partial u_{l,j}}{\partial w_{l,i,j}}=-\eta {\displaystyle \sum _{p=1}^P{\displaystyle \sum _{k=1}^{N_{L+1}}\frac{\partial J_p}{\partial u_{l+1,k}}}{u}_{l+1,j}(1-u_{l+1,j})·w_{l+1,k,j}{u}_{l,j}(1-u_{l,j})u_{l-1,i}}$$

(10)

According to Equations (9) and (10), the weight correction process is the back propagation of the error signal from the output layer to the input layer. When the error reaches the given accuracy, the training is terminated.

The traditional BP algorithm has a simple structure and is easy to use, with strong nonlinear mapping and generalization ability. However, it also has inherent defects in practical application, such as ease to fall into a local minimum, slow convergence speed, nonadjustable learning rate, short learning and memory time, etc. [29]. Therefore, aiming at the above defects, this paper appropriately optimizes and improves the traditional BP algorithm to increase the computational efficiency and stability of the algorithm.

4.2. Improved BP Algorithm

4.2.1. Introducing Gradient Factor into Sigmoid Function

For the sigmoid curve, there is a saturation region at each end of the definition domain. In the saturation region, the node output is not sensitive to the change of the activity value. If the error of the network is still large, it will be difficult to reduce the error through the correction of the weight in a short time. At this time, it can be determined that it has entered the flat area of the error surface [30], and the convergence speed of the network will be greatly affected. Therefore, a gradient factor $\delta$ can be introduced into the sigmoid function:

$$f(x)=\frac{1}{1+e^{-x/\delta }}$$

(11)

The gradient factor is generally initialized to 1 before training. If the network error variation is very small, but the value is large during training, it can be considered that it has entered the flat area of the error surface. It should appropriately increase $\delta$ to change the shape of the sigmoid function and increase the range of the unsaturated region. After exiting the flat area, $\delta$ will be restored to the initial value.

4.2.2. Introducing Adaptive Momentum

It can be seen from Equation (8) that the traditional BP algorithm only considers the current gradient direction of the error surface during weight iteration and ignores the previous gradient direction. It will result in oscillation effects in the training process, slow convergence, and will easily fall into a local optimum. [31]. Therefore, the momentum term is introduced on the basis of Equation (8). After improvement, the correction amount of the weight is determined by the modification direction of weight in the last time and the current gradient direction of the error surface. The momentum reflects the experience accumulated in the previous weight iteration and acts as a damping effect on the current correction. Based on the classical momentum method, this paper adopted the Nesterov momentum method [32] to iterate the weights. The classical momentum method can be expressed as:

$$\Delta w_{l,i,j}(n)=\alpha \Delta w_{l,i,j}(n-1)-\eta \frac{\partial J(w)}{\partial w_{l,i,j}}$$

(12)

where

$\alpha$ is the momentum coefficient, and $\alpha \in [0,1)$.

The Nesterov momentum method can be expressed as:

$$\Delta w_{l,i,j}(n)=\alpha \Delta w_{l,i,j}(n-1)-\eta \frac{\partial J(w+\alpha \Delta w_{l,i,j}(n-1))}{\partial w_{l,i,j}}$$

(13)

Comparing Equations (12) and (13), it can be found that compared with the classical method, the gradient was calculated after the current velocity in the Nesterov momentum method, so, it can be regarded as adding a correction factor to the classical momentum.

Then, the value of momentum coefficient was discussed. Generally speaking, it is often selected according to the experience of researchers. For example, Zhao et al. [33] considered that the convergence speed is faster when the momentum coefficient is between 0.9 and 0.98. The selection strategies of different scholars often lead to different training effects. In practical application, if the momentum coefficient is too small, the convergence speed will be reduced, and if it is too large, it will lead to divergence, and the fixed momentum coefficient is difficult to adapt to the shape change of the complex error surface. In this paper, the adaptive momentum coefficient was calculated according to the gradient of the error surface and the adjustment direction of the error function, which improved the training speed and stability of the network.

$$\alpha (n+1)=\left\{\begin{array}{cc}{e}^{-‖\frac{\partial J}{\partial w}‖}& J(n+1)<J(n)\\ 0& J(n+1)>1.04J(n)\\ \alpha (n)& \mathrm{other}\end{array}\right.$$

(14)

where $‖\frac{\partial J}{\partial w}‖$ is 2-norm of the error function gradient.

4.2.3. Adaptive Learning Rate

The learning rate is a crucial hyperparameter in the BP algorithm, and it directly affects the training performance of the network. It is generally taken as a fixed constant according to experience, but in applications, it is difficult to give an optimal learning rate from beginning to end [34]. In the initial stage or the flat area of the error surface, a larger learning rate should be selected to accelerate the training speed and avoid falling into a local optimum. At the later stage or in the steep area, an excessive learning rate will lead to oscillation effects and be difficult to converge. On this issue, this paper proposes a calculation method that adaptively adjusts the learning rate according to the gradient of the error surface and the number of iterations. The learning rate can be expressed as:

$$\eta =\sqrt{1-{\left(\frac{d-1}{d_{\max }}\right)}^2}\frac{2}{1+e^{\lambda \left|\frac{\partial J}{\partial w}\right|}}$$

(15)

where d is iterations. $d_{\max }$ is maximum iterations, and $\lambda$ is an empirical coefficient.

It can be seen from Equation (15) that the adaptive learning rate can flexibly adjust the iterative step size. Compared with the fixed learning rate algorithm, the training speed, accuracy, and stability all improved.

The meso parameter lightweight modeling process based on the improved BP algorithm is shown in Figure 5. For every independent variable, the variation of the predicted value of the network model caused by the change of its value is measured, in order to calculate and screen out the importance of each variable. The meso parameters with less importance were reasonably eliminated, making the analysis model simple and lightweight.

4.3. Lightweight Analysis Results

The single hidden layer neural network structure was used to perform macro–meso parameter lightweight analysis using traditional and improved BP algorithms, respectively. During training, the initial value of the momentum coefficient was 0.9, and the learning rate was 0.05. Figure 6 shows the lightweight analysis results obtained from the training. It can be seen from the figures that the conclusions drawn by the two analysis methods are basically the same, which mutually corroborate the correctness of the analysis results and indicate that the BP algorithm with global search capability is reasonably effective in lightweight analysis.

From the comparison of the two methods, it is clear from the results of the lightweight analysis of the elastic modulus (Figure 6a) that there are some differences between the traditional algorithm and the method in this paper in the calculation of the parallel bond stiffness ratio. The importance calculated by the traditional method is not prominent compared with other meso parameters (except $E^{*}$), while the lightweight result using the improved BP algorithm obviously reflects the importance of ${\overline{k}}_n/{\overline{k}}_s$ to the elastic modulus.

For Poisson’s ratio (Figure 6b), the two methods are mainly different in the lightweight results of the particle and parallel bond stiffness ratios. The improved BP algorithm shows that ${\overline{k}}_n/{\overline{k}}_s$ has the greatest impact on Poisson’s ratio, and it is also found that increasing ${\overline{k}}_n/{\overline{k}}_s$ can significantly improve Poisson’s ratio of rocks according to the research of this paper and the literature [11,12]. In terms of peak strength (Figure 6c), there is little difference between the two methods. In comparison, the calculation results of the improved BP algorithm are more distinct among the variables.

Figure 7 shows the iterative convergence curve of the two algorithms for the lightweight analysis of the elastic modulus, and the training process of Poisson’s ratio and peak strength are not shown due to space limitations. As can be seen from the figure, the improved BP algorithm converges quickly; the training process is stable, and it converges to an error value near 0.01 after about 50 iterations. On the other hand, the traditional algorithm decreases the error faster in the first 15 iterations, but the convergence of the network becomes significantly slower when it is iterated to about 15–60 times. It falls into the local optimum twice and takes several iterations to jump out. Once the number of iterations is greater than 60, the error starts to decrease slowly and tends to converge, but the network oscillation amplitude is large in the convergence process, and the oscillation amplitude can reach 22% when the iterations reach 80–120. It was accompanied by a small oscillation and convergence near the error value of 0.019 after the number of iterations reached 160. Compared with the method in this paper, the error increased by about 47%. In summary, it can be seen that the improved BP algorithm is significantly higher than the traditional method in terms of convergence speed, stability, and accuracy, and has significant advantages in the lightweight analysis of macro and meso parameters.

Based on the analysis results of the above two methods, it can be concluded that the meso parameters that have a significant impact on the elastic modulus are the equivalent modulus ($E^{*}$, ${\overline{E}}^{*}$), parallel bond stiffness ratio (${\overline{k}}_n/{\overline{k}}_s$), and control gap of the parallel bond ($\overline{d}$). The meso parameters that have a great influence on Poisson’s ratio are the parallel bond stiffness ratio (${\overline{k}}_n/{\overline{k}}_s$), particle stiffness ratio ($k_n/k_s$), and friction coefficient ($\mu$). The meso parameters that control the peak strength are the tensile and shear strengths of the parallel bond (${\overline{\sigma }}_{\mathrm{c}}$, ${\overline{\tau }}_{\mathrm{c}}$) and the internal friction angle ($\phi$).

To sum up, the stability and convergence speed of the traditional algorithm are increased by applying the improved BP algorithm introducing the sigmoid function with the gradient factor, the adaptive Nesterov momentum method, and the adaptive learning rate. It reasonably and effectively achieves the objective of lightweight and lays a foundation for an accurate analysis of macro and meso parameter influence characteristics and parameter calibration.

5. Analysis of Macro–Meso Parameter Correlation and Influence Characteristic

According to the orthogonal test results, the influence coefficients between macro and meso parameters as shown in

Table 3. Analysis of orthogonal test results.

Parameters	E/GPa		$\mathit{\nu }$		${\mathit{\sigma }}_{\mathit{c}}/\mathbf{MPa}$
	Influence Coefficient	t	Influence Coefficient	t	Influence Coefficient	t
constant	96.412	0.000	0.214	0.000	151.686	0.000
L/R	0.605	0.754	−0.011	0.038	−14.465	0.082
$R_{\max }/R_{\min }$	2.456	0.454	0.011	0.034	−12.155	0.161
$k_n/k_s$	−0.010	0.103	0.007	0.001	−19.434	0.058
${\overline{k}}_n/{\overline{k}}_s$	−28.877	0.000	0.047	0.000	−10.657	0.256
$\overline{d}$	14.583	0.001	0.006	0.862	−2.708	0.452
${\overline{\sigma }}_{\mathrm{c}}$	5.424	0.565	0.003	0.568	67.016	0.000
${\overline{\tau }}_{\mathrm{c}}$	1.199	0.621	−0.013	0.063	62.579	0.000
$\phi$	−0.359	0.802	0.001	0.058	3.045	0.016
$\mu$	−0.961	0.671	−0.044	0.000	5.310	0.224
$E^{*}$$({\overline{E}}^{*}$)	74.853	0.000	0.001	0.303	1.156	0.196

are such that the linear regression model between macro and meso parameters can be expressed as Equations (16)–(18).

$$\begin{array}{l}{E}^L=0.605L/R+2.456R_{\max }/R_{\min }-0.010k_n/k_s-28.877{\overline{k}}_n/{\overline{k}}_s+\\ 14.583\overline{d}+5.424{\overline{\sigma }}_c+1.199{\overline{\tau }}_c-0.359\phi -0.961\mu +74.583E^{*}+96.412\end{array}$$

(16)

$$\begin{array}{l}{\nu }^L=-0.011L/R+0.011R_{\max }/R_{\min }+0.007k_n/k_s+0.047{\overline{k}}_n/{\overline{k}}_s+\\ 0.006\overline{d}+0.003{\overline{\sigma }}_c-0.013{\overline{\tau }}_c+0.001\phi -0.044\mu +0.001E^{*}+0.214\end{array}$$

(17)

$$\begin{array}{l}{\sigma }_c{}^L=-14.465L/R-12.155R_{\max }/R_{\min }-19.434k_n/k_s-10.657{\overline{k}}_n/{\overline{k}}_s-\\ 2.708\overline{d}+67.016{\overline{\sigma }}_c+62.579{\overline{\tau }}_c+3.045\phi +5.310\mu +1.156E^{*}+151.686\end{array}$$

(18)

An analysis of variance was performed on the test results, and t in Table 3 was used to test whether the influence of meso parameters on macro parameters is significant, and 0.05 was used as the critical significance level. If t is less than 0.05, it is judged that the influence on macro parameters is significant. It can be seen from the table that the meso parameters that have a significant impact on the elastic modulus are: $E^{*}$(${\overline{E}}^{*}$), ${\overline{k}}_n/{\overline{k}}_s$, and $\overline{d}$. The meso parameters that have a significant impact on Poisson’s ratio are: ${\overline{k}}_n/{\overline{k}}_s$, $\mu$, $k_n/k_s$, $R_{\max }/R_{\min }$, and L/R. The meso parameters that have a significant impact on the peak strength are: ${\overline{\sigma }}_{\mathrm{c}}$, ${\overline{\tau }}_{\mathrm{c}}$, and $\phi$. In a word, this is basically consistent with the conclusion of the lightweight analysis.

5.1. Influence of Single Meso Parameters on Macro Parameters

According to the above lightweight analysis results, the selected meso parameters with significant impact were further studied, and the boxplot between macro and meso parameters was obtained, as shown in Figure 8, Figure 9 and Figure 10. The upper and lower edges of each box line in the figure, respectively, represent the maximum and minimum values at the same level, except for outliers. The bottom of the intermediate box is the 25% quantile, also known as the lower quartile. The top is the 75% quantile, also known as the upper quartile. The middle red line is the median, and the small circles represent outliers. Outliers are judged to be higher than the upper limit or lower than the lower limit. The upper and lower limits are calculated from the corresponding quartile and interquartile ranges (the difference between the upper and lower quartiles). The lower limit is the difference between the lower quartile and 1.5 times the interquartile range, and the upper limit is the sum of the upper quartile and 1.5 times the interquartile range.

(1)

Influence of Meso Parameters on Elastic Modulus

It can be seen from Figure 8a that the elastic modulus increases overall with the increase in the equivalent modulus. It can be divided into two stages. When the equivalent modulus is small, the increase is not obvious. When it is greater than 20 GPa (−0.87 in Figure 8a), the elastic modulus starts to increase significantly, and the median basically appears as linear growth. From each box, it can be seen that the macro elastic modulus at the same level is approximately normal distribution. When the equivalent modulus is small, the interquartile range is small. As the equivalent modulus increases, the elastic modulus distribution is also gradually discretized.

In the aspect of the parallel bond stiffness ratio (Figure 8b), the E shows a nonlinear decreasing trend with the increase in the stiffness ratio. When ${\overline{k}}_n/{\overline{k}}_s$≥ 2 (−0.13 in Figure 8b), this trend can be basically linearized. For each box, the distribution of the elastic modulus at the same level is relatively discrete, showing a right-skewed distribution. In general, the rate of the elastic modulus decreases with the stiffness ratio less than the increasing speed with the equivalent modulus. So, the equivalent modulus is the most important factor affecting E. This is consistent with the conclusion of the lightweight analysis using the improved BP algorithm.

The elastic modulus increases first and then decreases with the change of the control gap of the parallel bond (Figure 8c) and reaches the maximum when $\overline{d}$ = 0.5 (0 in Figure 8c). However, the lower quartile and lower limit were not significantly affected by $\overline{d}$. It can be considered that hard rock with a higher elastic modulus is more susceptible to the influence of the control gap of the parallel bond. At the same level, the elastic modulus has a right-skewed distribution, and if the elastic modulus is larger, the distribution will be more discrete.

(2)

Influence of Meso Parameters on Poisson’s Ratio

The particle stiffness ratio, parallel bond stiffness ratio, and friction coefficient have a great influence on Poisson’s ratio, as shown in Figure 9. Specifically, the variation trend of Poisson’s ratio with the particle and parallel bond stiffness ratios is basically the same, and both increase approximately linearly (Figure 9a,b). The difference is that under the influence of the particle stiffness ratio, Poisson’s ratio is roughly left-skewed, while under the influence of the parallel bond stiffness ratio, it is a right-skewed distribution, indicating that increasing ${\overline{k}}_n/{\overline{k}}_s$ can improve Poisson’s ratio more significantly.

With the increase in the friction coefficient, Poisson’s ratio basically shows a decreasing law of the quadratic curve (Figure 9c). The results of this paper also show that it has a positive excitation effect on the elastic modulus and peak strength. When the friction coefficient is small, Poisson’s ratio decreases greatly; when $\mu$ > 0.5 (−0.35 in Figure 9c), the decreasing trend is gradually gentle, and the value fluctuates around 0.2–0.3. This range is also Poisson’s ratio for common rock in engineering. So, it is suggested that 0.5 should be taken as the initial value when calibrating the friction coefficient.

(3)

Influence of Meso Parameters on Peak Strength

It is obvious that the peak strength is mainly affected by meso strength parameters, all of which have a positive excitation effect, as shown in Figure 10. It can be seen from Figure 10a that the peak strength increases with ${\overline{\sigma }}_{\mathrm{c}}$ and can be divided into two stages. When ${\overline{\sigma }}_{\mathrm{c}}$ < 50 MPa (−0.84 in Figure 10a), the increase in the peak strength is not obvious, and the data distribution is relatively concentrated. After 50 MPa, the growth rate increases significantly, and the data fluctuation range also increases significantly, indicating that the tensile strength of the parallel bond is not the only factor affecting the peak strength. Under the same level of ${\overline{\sigma }}_{\mathrm{c}}$, the changes and cross effects of other strength parameters also have an important impact on the peak strength.

In terms of the shear strength of the parallel bond, the peak strength shows a parabolic law with a gradually decreasing growth rate (Figure 10b). It increases most obviously when ${\overline{\tau }}_{\mathrm{c}}$ is in the range of 50–100 MPa (−0.84~−0.14 in Figure 10b). Then, the curve of the median and the lower quartile gradually flattened. The upper quartiles and maximum increased continuously. In the aspect of data distribution, it is basically similar to ${\overline{\sigma }}_{\mathrm{c}}$ and will not be repeated here.

The peak strength increases linearly with the increase in the internal friction angle (Figure 10c). At the same level, it shows a right-skewed distribution. The larger internal friction angle, the larger fluctuation range of the data. At the same time, the internal friction angle has a certain influence on the residual strength of rock, because the internal friction angle has a control effect on the particle sliding of the shear band when the sample is loaded after peak strength.

5.2. Influence of Multiple Meso Parameters on Macro Parameters

The preceding sections have provided a preliminary analysis of the variation law between macro and meso parameters. Generally speaking, macro parameters are often affected by multiple meso parameters. It is difficult for a single parameter to play a decisive role in a macro parameter. This also explains the reason for the large interquartile range in the boxplot. Therefore, this section will cross-combine the meso parameters screened above and study the variation law of macro parameters under the influence of multiple meso parameters.

(1)

Influence of Multiple Meso Parameters on Elastic Modulus

It can be seen from Figure 11a that when ${\overline{k}}_n/{\overline{k}}_s$ is fixed, the elastic modulus increases linearly with the increase in $E^{*}$. When ${\overline{k}}_n/{\overline{k}}_s$ is 0.8, the growth rate is the largest, and as ${\overline{k}}_n/{\overline{k}}_s$ increases, the growth rate gradually decreases. However, when the equivalent modulus is low, E is not significantly affected by the change of ${\overline{k}}_n/{\overline{k}}_s$; when it is greater than 100 GPa (0.64 in Figure 11a), the reduction range is relatively large. At the same level of $E^{*}$, ${\overline{k}}_n/{\overline{k}}_s$ can inhibit the growth of elastic modulus. From Figure 11b, we can see that the trends of each curve are basically the same, and E will be greatly reduced when the stiffness ratio is relatively high. On the other hand, the increase in the elastic modulus is relatively uniform due to the linear relationship between $E^{*}$ and E.

As seen in Figure 11c, when fixed $E^{*}$, the elastic modulus first increases and then decreases with the increase in $\overline{d}$, which is basically consistent with the rules of the boxplot in the previous section. In general, when $\overline{d}$ > 0.2 (−0.80 in Figure 11c), the impact on the elastic modulus is gradually reduced.

(2)

Influence of Multiple Meso Parameters on Poisson’s Ratio

According to previous analysis results, there is a linear relationship between Poisson’s ratio and $k_n/k_s$, ${\overline{k}}_n/{\overline{k}}_s$. As shown in Figure 12a,b, when $\mu$ is fixed, Poisson’s ratio increases linearly with the increase in the stiffness ratio, and curves are basically parallel. On the other hand, the increase in the friction coefficient will cause a nonlinear decrease in Poisson’s ratio (Figure 12c), and the smaller friction coefficient will result in a greater reduction range.

(3)

Influence of Multiple Meso Parameters on Peak Strength

Figure 13a shows the influence of ${\overline{\sigma }}_{\mathrm{c}}$ and ${\overline{\tau }}_{\mathrm{c}}$ on the peak strength. It can be seen from the figure that curves grow nonlinearly. When ${\overline{\tau }}_{\mathrm{c}}$ ≤ 50 MPa, the growth rate changes from large to small and tends to be flat after ${\overline{\sigma }}_{\mathrm{c}}$ > 50 MPa (−0.85 in Figure 13a). With the increase in ${\overline{\tau }}_{\mathrm{c}}$, the curve shape and growth rate also change. When ${\overline{\sigma }}_{\mathrm{c}}$ > 100 MPa (−0.14 in Figure 13a), the growth rate begins to decrease, and there is a sudden increase in peak strength when ${\overline{\sigma }}_{\mathrm{c}}$ in the range of 180–200 MPa (0.98–1.27 in Figure 13a). The reason for this phenomenon is that the rock has a strong hardening effect when the meso strength parameters are high.

It can be seen from Figure 13b that with fixed ${\overline{\sigma }}_{\mathrm{c}}$ and increased ${\overline{\tau }}_{\mathrm{c}}$, the peak strength increases cumulatively, and the slope of the curve is also from large to small. The peak strength increased most significantly when ${\overline{\tau }}_{\mathrm{c}}$ was in the range of 20–100 MPa (−1.27~−0.14 in Figure 13b). The trend of each curve is basically the same. We can see from the figure that the combination of ${\overline{\sigma }}_{\mathrm{c}}$ and ${\overline{\tau }}_{\mathrm{c}}$ at a higher level has a very significant effect on the improvement of peak strength.

Compared with ${\overline{\sigma }}_{\mathrm{c}}$ and ${\overline{\tau }}_{\mathrm{c}}$, the influence of the internal friction angle on peak strength is slightly weaker, which is also reflected in the lightweight analysis. Specifically, when ${\overline{\sigma }}_{\mathrm{c}}$ is fixed, the peak strength increases nonlinearly with $\phi$, and the curve is roughly “S”-shaped, with the maximum growth rate in the middle. When ${\overline{\sigma }}_{\mathrm{c}}$ ≤ 100 MPa, the growth of the curves is gentle. Only when the value of ${\overline{\sigma }}_{\mathrm{c}}$ is high, can changing the internal friction angle significantly improve the peak strength.

6. Response Surface Analysis Based on Central Composite Design Method

Through orthogonal test, lightweight analysis, and regression fitting, the linear relationship between macro and meso parameters was obtained; the meso parameters with significant influence were screened out, and the influence law between two scale parameters was analyzed. To achieve continuous, quantitative, and rapid matching, the accuracy of the linear regression model is not enough. The response surface method is a comprehensive optimization technique based on mathematical and statistical methods. Its basic idea is to obtain an approximate model between the research objective and its control factors through reasonable experimental schemes and iterative methods. The model can replace the complex implicit function between variables with an easy analyze function. Generally speaking, the second-order function can reflect the nonlinear relationship between input and output variables, consider the cross influence between input variables, and meet the accuracy requirements. Its general form can be expressed as:

$$y={\beta }_0+{\displaystyle \sum _{i=1}^k{\beta }_i{x}_i+{\displaystyle \sum _{i<j}{\displaystyle \sum _{j=2}^k{\beta }_{ij}{x}_i{x}_j+}}{\displaystyle \sum _{i=1}^k{\beta }_{ii}{x}_i{}^2}}+\epsilon$$

(19)

where x and y are the input and output variables, respectively. k is number of input variables. $\beta$ is the coefficient to be solved, and $\epsilon$ is the error term.

In order to obtain the second-order response surface function between macro and meso parameters, the central composite experimental design method (CCD) was used to select sample points. Based on the results of the lightweight analysis, the two most important meso parameters were selected for each macro parameter to perform a two-factor response surface analysis. For the elastic modulus, the selected meso parameters were the equivalent modulus ($E^{*}$) and parallel bond stiffness ratio (${\overline{k}}_n/{\overline{k}}_s$). For Poisson’s ratio, it was the particle ($k_n/k_s$) and parallel bond stiffness ratio (${\overline{k}}_n/{\overline{k}}_s$). The most significant influence on the peak strength are the tensile (${\overline{\sigma }}_{\mathrm{c}}$) and shear (${\overline{\tau }}_{\mathrm{c}}$) strengths of the parallel bond.

As shown in Figure 14, the two-factor CCD has 13 test points, including 5 center points, 4 factorial points, and 4 star points. The center points are used to provide error estimates; the factorial points are two-level factorial design, which are used to estimate the linear ${\beta }_i$ and interactive terms ${\beta }_{ij}$; the star points are used to expand the experimental area and estimate the squared terms ${\beta }_{ii}$. In order to make the experimental design satisfy the rotatability, the factor level at star points is determined according to Equation (20).

$$\alpha =2^{k/4}$$

(20)

Therefore, in this paper $\alpha =\sqrt{2}$.

In order to describe the quantitative relationships between macro and meso parameters more accurately, the research scope of the meso parameters was further subdivided into three equal subintervals, and the division of intervals is shown in

Table 4. Division of subintervals for response surface.

Term		S-1 *	S-2	S-3	S-4	S-5	S-6	S-7	S-8	S-9
E	$E^{*}$	[10, 56.67)	[56.67, 103.34)	[103.34, 150)	[10, 56.67)	[56.67, 103.34)	[103.34, 150)	[10, 56.67)	[56.67, 103.34)	[103.34, 150)
	${\overline{k}}_n/{\overline{k}}_s$	[0.8, 1.87)	[0.8, 1.87)	[0.8, 1.87)	[1.87, 2.94)	[1.87, 2.94)	[1.87, 2.94)	[2.94, 4.0)	[2.94, 4.0)	[2.94, 4.0)
$\nu$	$k_n/k_s$	[0.8, 1.87)	[1.87, 2.94)	[2.94, 4.0)	[0.8, 1.87)	[1.87, 2.94)	[2.94, 4.0)	[0.8, 1.87)	[1.87, 2.94)	[2.94, 4.0)
	${\overline{k}}_n/{\overline{k}}_s$	[0.8, 1.87)	[0.8, 1.87)	[0.8, 1.87)	[1.87, 2.94)	[1.87, 2.94)	[1.87, 2.94)	[2.94, 4.0)	[2.94, 4.0)	[2.94, 4.0)
${\sigma }_c$	${\overline{\sigma }}_{\mathrm{c}}$	[20, 80)	[80, 140)	[140, 200)	[20, 80)	[80, 140)	[140, 200)	[20, 80)	[80, 140)	[140, 200)
	${\overline{\tau }}_{\mathrm{c}}$	[20, 80)	[20, 80)	[20, 80)	[80, 140)	[80, 140)	[80, 140)	[140, 200)	[140, 200)	[140, 200)

* S-X is the abbreviation of Subinterval X.

. The corresponding second-order response surfaces were obtained in each subinterval, so that for each macro parameter, nine independent subresponse surfaces were used to describe it. The experimental design schemes and numerical calculation results of the elastic modulus, Poisson’s ratio, and peak strength are shown in

Table 5. CCD test scheme and calculation results of elastic modulus.

Type	${\mathit{E}}^{*}({\overline{\mathit{E}}}^{*})$	${\overline{\mathit{k}}}_{\mathit{n}}/{\overline{\mathit{k}}}_{\mathit{s}}$	Subinterval	E/GPa	Subinterval	E/GPa	Subinterval	E/GPa
Center points	0	0	$E^{*}$$({\overline{E}}^{*}$): 10–56.67 ${\overline{k}}_n/{\overline{k}}_s$: 0.8–1.87	55.23	$E^{*}$$({\overline{E}}^{*}$): 56.67–103.34 ${\overline{k}}_n/{\overline{k}}_s$: 0.8–1.87	132.62	$E^{*}$$({\overline{E}}^{*}$): 103.34–150 ${\overline{k}}_n/{\overline{k}}_s$: 0.8–1.87	209.57
	0	0		55.23		132.62		209.57
	0	0		55.23		132.62		209.57
	0	0		55.23		132.62		209.57
	0	0		55.23		132.62		209.57
Factorial points	1	1		86.56		157.88		228.86
	−1	1		15.15		86.55		157.88
	−1	−1		19.20		109.27		195.68
	1	−1		109.27		195.68		288.88
Star points	a	0		110.39		187.16		254.76
	0	a		49.43		118.87		187.66
	−a	0		0.47		78.09		155.54
	0	−a		70.59		169.44		267.26

,

Table 6. CCD test scheme and calculation results of Poisson’s ratio.

Type	${\mathit{k}}_{\mathit{n}}/{\mathit{k}}_{\mathit{s}}$	${\overline{\mathit{k}}}_{\mathit{n}}/{\overline{\mathit{k}}}_{\mathit{s}}$	Subinterval	$\mathit{\nu }$	Subinterval	$\mathit{\nu }$	Subinterval	$\mathit{\nu }$
Center points	0	0	$k_n/k_s$: 0.8–1.87 ${\overline{k}}_n/{\overline{k}}_s$: 1.87–2.94	0.157	$k_n/k_s$: 1.87–2.94 ${\overline{k}}_n/{\overline{k}}_s$: 1.87–2.94	0.189	$k_n/k_s$: 2.94–4.0 ${\overline{k}}_n/{\overline{k}}_s$: 1.87–2.94	0.194
	0	0		0.157		0.189		0.194
	0	0		0.157		0.189		0.194
	0	0		0.157		0.189		0.194
	0	0		0.157		0.189		0.194
Factorial points	1	1		0.177		0.223		0.211
	−1	1		0.113		0.177		0.223
	−1	−1		0.089		0.143		0.184
	1	−1		0.143		0.184		0.175
Star points	a	0		0.164		0.193		0.241
	0	a		0.175		0.208		0.213
	−a	0		0.171		0.165		0.176
	0	−a		0.127		0.156		0.162

and

Table 7. CCD test scheme and calculation results of peak strength.

Type	${\overline{\mathit{\sigma }}}_{\mathbf{c}}$	${\overline{\mathit{\tau }}}_{\mathbf{c}}$	Subinterval	${\mathit{\sigma }}_{\mathit{c}}/\mathbf{MPa}$	Subinterval	${\mathit{\sigma }}_{\mathit{c}}/\mathbf{MPa}$	Subinterval	${\mathit{\sigma }}_{\mathit{c}}/\mathbf{MPa}$
Center points	0	0	${\overline{\sigma }}_{\mathrm{c}}$: 20–80 ${\overline{\tau }}_{\mathrm{c}}$: 140–200	145.57	${\overline{\sigma }}_{\mathrm{c}}$: 80–140 ${\overline{\tau }}_{\mathrm{c}}$: 140–200	288.07	${\overline{\sigma }}_{\mathrm{c}}$: 140–200 ${\overline{\tau }}_{\mathrm{c}}$: 140–200	377.80
	0	0		145.57		288.07		377.80
	0	0		145.57		288.07		377.80
	0	0		145.57		288.07		377.80
	0	0		145.57		288.07		377.80
Factorial points	1	1		228.77		401.39		439.80
	−1	1		66.37		228.77		401.39
	−1	−1		63.73		226.46		319.73
	1	−1		226.46		319.73		344.94
Star points	a	0		271.85		360.32		379.53
	0	a		153.51		321.96		421.63
	−a	0		28.64		182.33		353.93
	0	−a		132.91		274.11		310.06

, respectively. Due to space limitations, for each macro parameter, only three representative subintervals were selected to show the calculation results, and the rest of the calculation results were omitted.

(1)

Response Surface Analysis of Elastic Modulus

According to the CCD test results,

Table 8. Response surface equation of elastic modulus.

Terms	S-1	S-2	S-3	S-4	S-5	S-6	S-7	S-8	S-9
$E^{*}{}^2$	0.05	−0.27	−1.85	0.09	0.11	−0.10	0.09	0.02	0.42
${({\overline{k}}_n/{\overline{k}}_s)}^2$	2.35	5.53	9.38	0.45	1.19	2.28	0.17	0.46	0.89
$E^{*}·({\overline{k}}_n/{\overline{k}}_s)$	−4.67	−3.77	−5.55	−1.67	−1.68	−1.66	−0.87	−0.93	−0.89
$E^{*}$	39.68	39.06	38.12	33.93	33.82	32.05	31.48	31.32	31.22
${\overline{k}}_n/{\overline{k}}_s$	−7.10	−16.53	−26.34	−2.43	−5.80	−9.23	−1.29	−3.08	−4.92
Constant	55.23	132.62	209.56	47.90	115.03	181.77	44.52	107.06	168.76
$R^2$	1.0	1.0	0.99	1.0	1.0	1.0	1.0	1.0	1.0

lists the coefficients of each subresponse surface equation for the elastic modulus. The R² of each response surface equation is above 0.99, and the fitting is relatively accurate. These equations will be used later in the parameter calibration solution. It can be seen from the Figure 15 that the elastic modulus is positively correlated with $E^{*}$ and is negatively correlated with ${\overline{k}}_n/{\overline{k}}_s$. On the whole, it is greatly affected by $E^{*}$. When the equivalent modulus is low, the effect of ${\overline{k}}_n/{\overline{k}}_s$ on the elastic modulus is not obvious (Figure 15a,d,g), and the gradient direction of the response surface is gradually deflected with the increase in $E^{*}$. As shown in Figure 15h,f,i, the elastic modulus is affected by two parameters. So, it is necessary to use multiple subresponses to describe the relationship between the macro and meso parameters. It is suggested that the differences in the influence of meso parameters on macro parameters in different regions should be considered when performing parameter calibration.

(2)

Response Surface Analysis of Poisson’s Ratio

Similarly, the coefficients of Poisson’s ratio response surface equation (

Table 9. Response surface equation of Poisson’s ratio.

Terms	S-1	S-2	S-3	S-4	S-5	S-6	S-7	S-8	S-9
${(k_n/k_s)}^2$	−0.001	−0.004	−0.001	−0.002	−0.004	0.008	−0.002	−0.005	0.009
${({\overline{k}}_n/{\overline{k}}_s)}^2$	−0.018	−0.015	−0.015	−0.010	−0.003	−0.003	−0.010	−0.001	−0.001
$(k_n/k_s)·({\overline{k}}_n/{\overline{k}}_s)$	0.007	0.004	−0.002	0.002	0.001	−0.001	0.001	0.001	0.001
$k_n/k_s$	0.010	0.013	0.003	0.013	0.016	0.009	0.014	0.017	0.009
${\overline{k}}_n/{\overline{k}}_s$	0.043	0.050	0.052	0.016	0.018	0.018	0.009	0.010	0.010
Constant	0.105	0.134	0.141	0.157	0.188	0.193	0.181	0.215	0.219
$R^2$	0.98	0.99	0.99	0.94	0.99	0.98	0.94	0.98	0.95

) and the contour plot of the response surface were obtained based on the experimental results. As seen in Figure 16, the gradient direction of Poisson’s ratio response surface is also deflected with the change of the two meso parameters. When ${\overline{k}}_n/{\overline{k}}_s$ is small (Figure 16g,h,i), Poisson’s ratio is less affected by $k_n/k_s$; as ${\overline{k}}_n/{\overline{k}}_s$ increases, Poisson’s ratio also increases, and the effect of $k_n/k_s$ on Poisson’s ratio starts to be significant. As shown in Figure 16b,d,e, the gradient of the contour line is basically along 45°, indicating that the gradient of Poisson’s ratio response surface along the two coordinate axes is equivalent, and ${\overline{k}}_n/{\overline{k}}_s$ and $k_n/k_s$ have equally significant effects on Poisson’s ratio. It is worth noting that, as shown in Figure 16a, Poisson’s ratio in the local range has a quadratic curve law that increases first and then decreases slightly with the increase in ${\overline{k}}_n/{\overline{k}}_s$. It may be due to the joint influence of the particle stiffness ratio and the parallel bond stiffness ratio. This is of concern in the parameter calibration.

(3)

Response Surface Analysis of Peak Strength

On the other hand, the response surface equation of peak strength was obtained by regression (

Table 10. Response surface equation of peak strength.

Terms	S-1	S-2	S-3	S-4	S-5	S-6	S-7	S-8	S-9
${\overline{\sigma }}_c{}^2$	−16.014	−0.855	2.285	0.369	−13.987	−6.038	2.224	−5.991	−3.007
${\overline{\tau }}_c{}^2$	−12.968	3.448	12.176	−4.608	−12.350	−9.330	−1.265	7.254	−3.446
${\overline{\sigma }}_c·{\overline{\tau }}_c$	26.172	13.147	2.097	11.630	10.110	4.209	−0.082	19.837	3.301
${\overline{\sigma }}_c$	31.604	7.847	1.230	75.071	37.051	9.170	83.469	64.564	12.445
${\overline{\tau }}_c$	35.424	72.499	85.166	10.367	42.507	52.815	4.258	18.912	41.697
Constant	107.855	121.826	122.321	135.680	260.799	293.830	145.569	288.096	377.836
$R^2$	0.99	0.99	1.0	0.99	0.99	1.0	1.0	0.99	0.98

), and it can be seen that for each subresponse surface equation, the R² > 0.98 and most of them are above 0.99. The regression effect of the equation is relatively significant. As can be seen from Figure 17, the peak strength is controlled by ${\overline{\sigma }}_{\mathrm{c}}$ and ${\overline{\tau }}_{\mathrm{c}}$ together, both of which are positively correlated. When the magnitude of the two is small (Figure 17d,e,g), ${\sigma }_{\mathrm{c}}$ is jointly affected by the two meso parameters. With the increase in a certain meso parameter, the gradient direction of the response surface will deflect, showing a “saturation” effect, that is, the peak strength is mainly controlled by ${\overline{\tau }}_{\mathrm{c}}$ when ${\overline{\sigma }}_{\mathrm{c}}$ is large and vice versa (Figure 17a,f,i). This phenomenon is obvious in the process of ${\overline{\sigma }}_{\mathrm{c}}$ changed.

7. Calibration and Solution Method of Macro and Meso Parameters

To calibrate the meso parameters of rock, the optimization technique was used to solve the meso parameters, and the objective of the solution was to minimize the error between the numerical and laboratory test that is expressed as a function:

$$f(x)=\min (y_s-y_l)$$

(21)

where $y_s$ is the numerical test result, and here the subresponse surface equations were used. $y_l$ is the macroscopic parameter of the rock obtained by laboratory tests.

In order to solve the objective function, the following boundary conditions were defined for the problem.

1.

Linear equality constraint.

The macro parameters obtained by the linear regression model should be equal to the laboratory test.

$$\left\{\begin{array}{l}{E}^L-E^{lab}=0\\ {\nu }^L-{\nu }^{lab}=0\\ {\sigma }_c{}^L-{\sigma }_c{}^{lab}=0\end{array}\right.$$

(22)

2.

Linear inequality constraints.

According to the literature [15,25,35] and the research of this paper, it was found that there is a certain empirical relationship between the elastic modulus and the peak strength of rock, with the ratio of the two being about 0.25~1.0, and the ratio of ${\overline{\sigma }}_{\mathrm{c}}$ and ${\overline{\tau }}_{\mathrm{c}}$ is 1~3. Therefore, the linear regression model has the following constraints

$$\left\{\begin{array}{l}0.25{\sigma }_c{}^L-E^L\le 0\\ E^L-{\sigma }_c{}^L\le 0\end{array}\right.$$

(23)

$$\left\{\begin{array}{l}{\overline{\sigma }}_{\mathrm{c}}-{\overline{\tau }}_{\mathrm{c}}\le 0\\ -3{\overline{\sigma }}_{\mathrm{c}}+{\overline{\tau }}_{\mathrm{c}}\le 0\end{array}\right.$$

(24)

3.

Nonlinear inequality constraint.

Applying the inequality constraints above to the nonlinear models in Table 8, Table 9 and Table 10 constitutes the nonlinear inequality constraint.

$$\left\{\begin{array}{l}0.25{\sigma }_c-E\le 0\\ E-{\sigma }_c\le 0\end{array}\right.$$

(25)

4.

Side constraints.

As mentioned above, each meso parameter has its research range. This constraint was determined by the research range in Section 3.2.

The above constitutes the mathematical model for the parameter calibration solution, which contains three sets of objective functions, three equality constraints, six inequality constraints, and three side constraints. It is a typical multiobjective, nonlinear optimization problem. The above mathematical model was solved by using the multiobjective goal attainment solution function (fgoalattain) in the Matlab optimization toolbox. The optimal meso parameters can be iteratively solved by setting the solution objective value (0, which means the minimum value of the objective function), the initial value (all set to 1), and the weight (each objective function has the equal weight) and specifying the corresponding objective functions, constraints, and boundary conditions according to the above-calibrated solution model.

8. Case Study

8.1. Calibration of Rock Meso Parameters

The griotte researched in the literature [12] was selected for calibration, and its macro parameters are: E = 39.5 GPa, $\nu$ = 0.12, and ${\sigma }_{\mathrm{c}}$ = 111 MPa. The meso parameters were obtained after optimal solution using the mathematical model established above, as shown in

Table 11. Mesoscopic parameters of griotte.

L/R	${\mathit{R}}_{\mathbf{max}}/{\mathit{R}}_{\mathbf{min}}$	${\mathit{E}}^{*}({\overline{\mathit{E}}}^{*})/\mathbf{GPa}$	${\mathit{k}}_{\mathit{n}}/{\mathit{k}}_{\mathit{s}}$	${\overline{\mathit{k}}}_{\mathit{n}}/{\overline{\mathit{k}}}_{\mathit{s}}$	$\overline{\mathit{d}}$	${\overline{\mathit{\sigma }}}_{\mathbf{c}}/\mathbf{MPa}$	${\overline{\mathit{\tau }}}_{\mathbf{c}}/\mathbf{MPa}$	$\mathit{\phi }$/°	$\mathit{\mu }$
50	3.67	25.6	1.86	3.0	0.29	140	106.7	47	0.57

.

8.2. Analysis of Calibration Results

The calculated meso parameters were input into the PFC for a triaxial compression test, and the confining pressure was set to 10 MPa. The deviatoric stress–strain curve, contact force chain, and crack distribution obtained from the test are shown in Figure 18. From the figure, it can be seen that

(1)

The sample is ductile failure as a whole, and the compression contact force chain between the particles is dominant. Some tensile contact force chains appear in the middle of the sample, which reflects the shear expansion characteristics during compression.

(2)

According to numerical calculation results, the macro parameters of the sample are: E = 42.9 GPa, $\nu$ = 0.11, and ${\sigma }_{\mathrm{c}}$ = 112.9 MPa. Compared with laboratory test results, the absolute errors are: 3.4 GPa, 0.01, and 1.99 MPa; The relative errors are 8.6%, 8.3%, and 1.7%, respectively. The main reason for the high elastic modulus and strength and low Poisson’s ratio in the numerical test is that the numerical test is easier to obtain uniform and complete samples, so, its physical and mechanical parameters are also slightly better than laboratory test. In general, the numerical test results can accurately reflect the stress evolution law of rock.

(3)

The cracks of the sample mainly developed along the bottom surface in the direction of $63°$. In addition, there are a few cracks in the upper left corner of the sample, and finally an obvious main fracture surface was formed as the numerical experiments were conducted. The damage pattern and crack distribution of the sample are relatively close to the laboratory tests.

In summary, the above method can quickly and accurately obtain the corresponding meso parameters of the parallel bond model based on the macroscopic parameters of the rock. The comparison of the numerical tests with the laboratory tests also reveals that the stress characteristics, stress evolution, and failure mode of the sample are relatively close, which proves that this method is reasonable and effective.

9. Conclusions

In this paper, starting from the basic theory of PFC3D and qualitative relationship between macro and meso parameters, based on the parallel bond model and reasonable experimental design methods, a lightweight analysis method based on the improved BP algorithm is proposed, and the sensitivity between macro and meso parameters is quantified. The influence law of meso parameters on macro parameters was studied, and the linear and nonlinear models between them were obtained. The mathematical model of rock meso parameter calibration was established, and the optimal solution was obtained by using optimization technology, which can provide some reference for PFC parameter calibration. The main conclusions of this paper are as follows:

(1)

An improved BP algorithm is proposed that has a function with gradient factor, Nesterov momentum method with adaptive momentum coefficients, and adaptive learning rate and applied to the PFC meso parameter lightweight analysis. In the comparison with the traditional BP algorithm, it was found that the convergence speed is increased by 140%; the oscillation effect in the iteration is avoided, and the error is reduced by about 47%. The results are more accurate and reasonable and effectively achieve the purpose of lightweighting the macro and meso parameter analysis model. At the same time, the algorithm proposed in this paper is also applicable to similar problems in the field of water conservancy and geotechnical engineering.

(2)

The improved BP algorithm was used for the lightweight analysis, and it was found that the equivalent modulus, parallel bond stiffness ratio, and parallel bond control gap have a significant impact on the elastic modulus, in which the equivalent modulus plays a leading role. On the other hand, the improved algorithm obviously reflects the importance of the stiffness ratio to the elastic modulus compared with the traditional method. The meso parameters that have a great influence on Poisson’s ratio are the parallel bond stiffness ratio, particle stiffness ratio, and friction coefficient. The algorithm in this paper is slightly more accurate in describing the importance of the two stiffness ratios. The compressive strength, shear strength, and internal friction angle of the parallel bond can control the peak strength.

(3)

The interaction law and influence mechanism between macro and meso parameters were analyzed. The elastic modulus increases linearly with the increase in $E^{*}$; ${\overline{k}}_n/{\overline{k}}_s$ has an inhibitory effect on the growth of the elastic modulus, and E decreases greatly when the stiffness is relatively high. There is a linear relationship between Poisson’s ratio and $k_n/k_s$ and ${\overline{k}}_n/{\overline{k}}_s$. The increase in the friction coefficient will cause the nonlinear decrease in Poisson’s ratio, and the smaller friction coefficient, the greater reduction. All the strength parameters have an incentive effect on the peak strength, and the combination of ${\overline{\sigma }}_{\mathrm{c}}$ and ${\overline{\tau }}_{\mathrm{c}}$ at higher levels has the most significant effect on the enhancement of peak strength. On the other hand, the internal friction angle has a certain effect on the residual strength of the rock, which is due to the control effect of the internal friction angle on the particle sliding on both sides of the shear zone when the sample is loaded after the peak.

(4)

The research interval of the meso parameters with significant effects was further subdivided, and based on the central composite experimental design and response surface method, a nonlinear mathematical model of macro–meso parameters described by multiple subresponse surfaces was obtained to achieve a refined description between macro and meso parameters.

(5)

The mathematical model of the meso parameter calibration is established, solved by optimization technology, and verified by an example. Through the comparison of numerical experiment and laboratory test results, it was also found that the mechanical characteristics, stress evolution, and failure mode of the sample are relatively close, which proves the validity and reliability of the calibration method.

(6)

Limitations and Prospects. ① The research results of this paper are mainly based on the parallel bond model, and the applicability in other contact models of PFC remains to be verified. ② Less consideration was given to the damage evolution and failure morphology matching of rocks during loading, which will be improved in future research. ③ The influence of clump properties was not considered, and the role of complex particles in the calibration of rock macro and meso parameters needs further study in the future.