Microparameters Calibration for Discrete Element Method Based on Gaussian Processes Response Surface Methodology

Abstract

Microparameter calibration is an important problem that must be solved in the discrete element method. The Gaussian process (GP) response surface methodology was proposed to calibrate the microparameters based on the Bayesian principle in machine-learning methods, which addresses the problems of uncertainty, blindness, and repeatability in microparameter calibration methods. Using the particle flow code (PFC) as an example, the effects of the microparameters on the macroparameters were evaluated using the control-variable method, and the range of the microparameters was determined based on the macroparameters. The uniform design (UD) method and numerical calculation were used to obtain training samples, and a GP response surface methodology suitable for multifactor, multilevel, and nonlinear processes was used to establish the response surface relationships for macro–micro parameters of rock-like materials in discrete element method. According to the macroparameters obtained from the uniaxial experiments conducted on rock specimens, the microparameters were calibrated using the GP response surfaces. Numerical calculations of uniaxial compression and Brazilian splitting were performed using microparameters, and the results were compared with laboratory experiments for verification. The results showed that the relative errors of the GP response surface and laboratory test values were 5.3% for the modulus of elasticity, −7.8% for compressive strength, and −2.6% for tensile strength. The nonlinear GP response surface considered the characteristics of multiple interacting factors, and the established nonlinear response surface relationship between the microparameters and macroparameters can be used for the calibration of microparameters. The accuracy of the microparameters was verified according to the stress–strain curve and failure morphology of the rock specimens. The method of using the GP response surface to establish the macro–micro parameter relationship in the discrete element method can also be extended to other numerical simulation methods and can provide a basis for accurately analysing the microdamage mechanism of rock materials under complex loading conditions.

1. Introduction

Numerical simulations have become an effective tool for studying the microdamage mechanisms of rocks [1]. Among them, the discrete element method can more realistically reflect the microscopic characteristics of materials using a simple contact model, which is convenient for dealing with deformation problems and the associated damage related to continuous and discontinuous media [2,3,4,5]. It is widely used in civil engineering, transportation, water conservation, and other engineering fields. Among these discrete element methods, the representative particle flow code (PFC) is used to study the mechanical properties by simulating the mutual motion of particles after the material is discretised into rigid particles and using simple contact models, including the contact-bond model (CBM), parallel-bond model (BPM), and flat-joint model (FJM). However, the prerequisite is to select the appropriate microparameters [6,7,8,9,10,11,12,13]. Because these macroparameters, such as the elastic modulus, Poisson’s ratio, cohesion, and friction angle, no longer have an explicit correspondence with the microparameters in PFC in conventional numerical calculation methods based on continuous media, a parameter calibration needs to be carried out. Researchers have used trial and error [14,15,16], simulation experiments [17,18,19], and experimental designs [20,21,22,23], such as the Plackett–Burman design, central group and design, and orthogonal design method, to conduct relevant studies on parameter calibration. However, multifactor interactions, poor generality of parameters for multiple contact models, simplicity, and efficiency require further exploration. An improved method of establishing the relationship between macro–micro parameters in discrete element methods must be sought to improve accuracy, reduce the number of tests, and simplify the calibration method.

Because of the interaction of several microparameters, manual calibration methods are uncertain, blind, repetitive, and pose several other problems. Therefore, researchers must seek more effective methods for the quantitative calibration of microparameters. Intelligent algorithms based on machine learning are more adept at dealing with multiparameter interactions. Among them, the Gaussian process (GP), which is one of the most important Bayesian machine-learning approaches based on a particularly effective method for placing a prior distribution over the space of functions, is a generalization of the Gaussian probability distribution [24], which has been rapidly developed in recent times as a new machine-learning methodology and has been successfully applied in many fields [25]. The GP model is determined by the mean and covariance functions, which are not restricted by a specific form of mathematical structure and can effectively model the correlations among complex physical systems. In addition, it has better adaptability to deal with small samples, highly nonlinear, and high-dimensional numbers, which can provide a variance corresponding to the calculated results and a confidence interval of the given probability. The factorial design of a GP response surface is based on the statistical variance analysis method and eliminates unimportant correction parameters, which can reduce the size (dimensionality) of the optimisation problem. Additionally, it has an advantage over traditional correction methods, neural networks, support-vector machines, and other methods.

The aim of this study is to propose a novel method (response surface method) for microparameter calibration in the discrete element method, which is based on the GP in machine learning. This method has the advantages of small samples, multifactor, multilevel, and high fitting accuracy, and is suitable for establishing a nonlinear response surface between two datasets with high generality [26]. The basic theory is based on the GP response surface, and the general steps of microparameter calibration are introduced in this study. The effect of the microparameters on the macroparameters was analysed using the control-variable method to obtain the range of microparameters. Training datasets were generated using the uniform design (UD) method and numerical calculations to establish a GP response surface for macro–micro parameters of rock, which used the discrete element method. The microparameters were calibrated from the macroparameters obtained from laboratory uniaxial compression and Brazilian splitting experiments using GP response surfaces. The stress–strain curves of rock specimens under uniaxial compression and Brazilian splitting were calculated and compared with laboratory experimental results to verify the accuracy of the macro–micro parameter response relationships established by the GP response surfaces. The method, as an effective method to solve the parametric inverse analysis problem, can also be extended to other numerical models for parameter calibration. The research is helpful to those who are interested in analysing the multicrack extension process and reveals its microscopic fracture mechanism under complex loading conditions using the discrete element method in geotechnical engineering.

2. Gaussian Process Theory

2.1. Gaussian Process Model

A Gaussian process is a collection of random variables, any finite number of which have a joint Gaussian distribution [27]. There is a training dataset $D$ of $n$ observations, $D=\left\{({\mathbf{x}}_i,y_i)\left|i=1,2,\dots ,n\right.\right\}$, where $\mathbf{x}=(x_1,x_2,\dots ,x_n)$ denotes an input vector (covariates) of dimension $d$, that is, the number of variables, $x$ denotes any random variable, and $y$ denotes a scalar output or target (dependent variable); the column vector inputs for all $n$ cases are aggregated in the $d\times n$ design matrix $X=({\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_3)$, and the targets are collected in the vector $\mathbf{y}=(y_1,y_2,\dots ,y_n)$, so we can write $D=(X,\mathbf{y})$. Given this training dataset $D$ and outputs ${\mathbf{y}}^{\ast }$, we wish to make predictions for new inputs ${\mathbf{x}}^{\ast }$ that we have not seen in the training dataset [25].

A Gaussian process is completely specified by its mean function and covariance function. Define the mean function $m(\mathbf{x})$ and the covariance function $k(\mathbf{x},\mathbf{x}’)$ of a real process $f(\mathbf{x})$ as

$$\left.\begin{array}{l}m(\mathbf{x})=E[f(\mathbf{x})]\\ k(\mathbf{x},\mathbf{x}’)=E[(f(\mathbf{x})-m(\mathbf{x}))(f(\mathbf{x}’)-m(\mathbf{x}’))]\end{array}\right\}$$

(1)

where $E[f(\mathbf{x})]$ denotes the mean of the function $f(\mathbf{x})$, and write the Gaussian process as

$$f(\mathbf{x})\sim GP(m(\mathbf{x}),k(\mathbf{x},\mathbf{x}’))$$

(2)

$f(\mathbf{x})$ and $f(\mathbf{x}’)$ are jointly Gaussian with zero mean and covariance, and the function values $f({\mathbf{x}}_1),\dots ,f({\mathbf{x}}_n)$ corresponding to any number of input points $n$ are jointly Gaussian. Given a set of input points $\left\{{\mathbf{x}}_i\left|i=1,\dots ,n\right.\right\}$ we can compute the Gram matrix $\mathbf{K}$ whose entries are ${\mathbf{K}}_{ij}=k({\mathbf{x}}_i,{\mathbf{x}}_j)$. If $k$ is a covariance function, we call the matrix $\mathbf{K}$ the covariance matrix. A real $n\times n$ matrix $\mathbf{K}$ is called positive semidefinite.

2.2. Selection of Covariance Function

From Equations (1) and (2), the GP was determined using the mean and covariance functions. The covariance function (kernel function) plays a decisive role in the prediction results of the GP model and measures the closeness between the training and prediction dataset. Covariance functions are typically used in the form of square exponential, linear, and periodic functions. A covariance function can also be generated by summing several functions, products, and other operations of the new covariance function, reflecting the flexibility of the GP kernel function.

In this study, we used the squared exponential covariance function, which has the advantage of being a smooth and infinitely derivable function and has the characteristic that neighbouring inputs produce similar outputs. making it the most used covariance function with a form shown in Equation (3) [28]:

$$k({\mathbf{x}}_i,{\mathbf{x}}_j)={\psi }^2\exp [-\frac{1}{2}{\displaystyle \sum _{m=1}^d{(\frac{{\mathbf{x}}_{i,m}-{\mathbf{x}}_{j,m}^{}}{l_m})}^2}],$$

(3)

where $x_{i,m}$ is the m-th variable of the matrix ${\mathbf{x}}_i=\left(x_{i,1}^{},x_{i,2}^{},\dots ,x_{i,d}^{}\right)$; $d$ is the dimensionality of the input space; ${\psi }^2$ is the signal variance, which is used to control the degree of local correlation; and $(l_1,l_2,\dots ,l_m,\psi )$ is the hyperparameter inside the covariance function and set $\theta =(l_1,l_2,\dots ,l_m,\psi )$ to simplify the expression.

2.3. Solution of the Hyperparameter

The hyperparameters determine the expression of the GP, that is, the GP response surface; thus, an optimal choice of the covariance function and calculation method for the hyperparameters are required. In practical applications, the mean function is set as a constant to simplify the calculation process.

To optimise the prediction ${\mathbf{x}}^{\ast }$, the hyperparameter $\theta$ is often determined using the great likelihood method, which is generally performed by taking the negative logarithm of the marginal-likelihood function. The negative logarithm is the likelihood function. The partial derivatives of the hyperparameters are obtained, which are further optimised by the Newton, conjugate gradient, and other methods to find the optimal solution of the hyperparameters [29]. The optimal hyperparameter $\theta$ can be expressed as Equation (4):

$$\widehat{\theta }={\arg }_{\theta }\min L(\theta ),$$

(4)

where

$$L(\theta )=-\lg [p(\mathbf{y}|X,\theta )]=\frac{1}{2}{\mathbf{y}}^T{\mathbf{K}}^{-1}\mathbf{y}+\frac{1}{2}\lg (|\mathbf{K}|)+\frac{n}{2}\lg (2\pi ).$$

(5)

The partial derivatives of the hyperparameters are as follows:

$$\frac{\partial L(\theta )}{\partial {\theta }_n}=\frac{1}{2}tr({\mathbf{K}}^{-1}\frac{\partial \mathbf{K}}{\partial {\theta }_n})-\frac{1}{2}{\mathbf{y}}^T{\mathbf{K}}^{-1}\frac{\partial \mathbf{K}}{\partial {\theta }_n}{\mathbf{K}}^{-1}\mathbf{y},$$

(6)

where ${\arg }_{\theta }\min L(\theta )$ denotes the value of the variable when the function $L(\theta )$ is minimized, tr$(•)$ denotes the trace, and ${\theta }_n$ denotes the n-th optimal hyperparameter.

The iteration for the training-sample set optimisation is an ideal method for solving the hyperparameters, which is also a significant advantage in Gaussian algorithms [28,30]. Its specific operation is as follows: 100 sets of initial points are randomly selected to calculate the corresponding $L(\theta )$. Thereafter, ten sets of initial points corresponding to the smallest $L(\theta )$ are selected, and, among the ten sets of results, the smallest $L(\theta )$ is the optimal hyperparameter $\widehat{\theta }$.

3. Selection of Parameters

The process of microparameter calibration first requires the establishment of a Gaussian response surface using the microparameter set ($X$) and the macroparameter set ($\mathbf{y}$), which are generated by the UD method and numerical simulation, respectively, as shown in Figure 1.

3.1. Selection of Microparameters

PFC simulates the mechanical behaviour of continuous-media materials by bonding particle elements that are independent of each other. The primary bonding models commonly used to simulate rock materials are CBM and PBM. The difference between CBM and PBM is shown in Figure 2. A CBM can be regarded as a spring with normal and tangential stiffness, which can resist normal and tangential forces. The PBM can be considered as a series of springs with normal and tangential stiffnesses uniformly distributed over a certain width of the contact surface, which resists normal forces, tangential forces, and moments simultaneously. In the PBM, the breakdown of the PBM immediately leads to a reduction in macroscopic stiffness, which is more consistent with the brittle failure of rock and discontinuous mediums after failure. Therefore, the PBM was chosen for this study to simulate the mechanical behaviour of rocks during failure.

The parameters of the particles and parallel bonds are included in the PBM (see

Table 1. Microparameters of the PBM.

Parameters of the Particles	Parameters of the Parallel Bonds
Modulus of the particles: $E_c$	Modulus of parallel bonds: ${\overline{E}}_c$
Ratio of normal to tangential stiffness of the particles: $k_n/k_s$	Ratio of normal to tangential stiffness of parallel bonds: ${\overline{k}}_n/{\overline{k}}_s$
Friction coefficient of the particles: $u$	Average value of normal strength of parallel bonds: ${\overline{\sigma }}_c$
	Average value of tangential strength of parallel bonds: ${\overline{\tau }}_c$
	Radius increase factor of parallel bonds: $\overline{\lambda }$

). To simplify the calculation, it was necessary to assume that the particles and parallel bonds have the same deformation properties, that is, $E_c={\overline{E}}_c$ and $k_n/k_s={\overline{k}}_n/{\overline{k}}_s$. It was also assumed that the radius of the parallel bonds is the same as the radius of the smaller particles in contact with each other, that is, $\overline{\lambda }=1.0$. Therefore, the microparameters ${\overline{E}}_c$, ${\overline{k}}_n/{\overline{k}}_s$, ${\overline{\sigma }}_c$, ${\overline{\tau }}_c$ and $u$ were selected as test factors.

3.2. Range Determination of Microparameters

The training set ($D$) needs to be established by substituting the microparameter set ($X$) into the discrete element numerical model to generate the macroparameter set ($\mathbf{y}$). Therefore, it is necessary to first determine the range of the microparameters and subsequently generate the training set ($X$) within this range. The control-variable method was used to determine the range of microparameters by calculating the linear relationship between the microparameters and the macroparameters, whereafter the range of the microparameters based on the empirical values of the macroparameters is deduced.

The Initial microparameters were ”et a’: $E_c={\overline{E}}_c=50\mathrm{GPa}$, $k_n/k_s={\overline{k}}_n/{\overline{k}}_s=1.0$, $u=0.5$, $R_{\min }=0.3\mathrm{mm}$, and $R_{\max }/R_{\min }=1.66$. The single-factor analysis method, where only one microparameter (${\overline{\sigma }}_c$) is changed repeatedly while other microparameters are kept constant, was used to calculate the linear relationship. The range ${\overline{\sigma }}_c={\overline{\tau }}_c$ was set to 20–300 MPa. Numerical calculation models were established to simulate uniaxial compression experiments to obtain macroparameters and study the influence of microparameters on macroparameters. The relationship between the parallel-bond and uniaxial compression strengths were obtained, as shown in Figure 3. Linear fitting of the numerical test results yielded $y=-0.4+1.65x$. According to the petrophysical tests results of the rock mass [31], the uniaxial compression strengths of marble range from 40 to 120 MPa. Substituting $x=40$ and $x=120$ into $y=-0.4+1.65x$ yields $y=24.1$ and $y=72.3$, respectively. Therefore, the range of values for the cohesive strengths of ${\overline{\sigma }}_c$ and ${\overline{\tau }}_c$ was set to 24.1–72.3 MPa. Following similar calculations, the range of the parallel-bond modulus (${\overline{E}}_c$) was set to 34.35–103.05 GPa, and that of the stiffness ratio (${\overline{k}}_n/{\overline{k}}_s$) to 0.455–1.365, with their fitting lines shown in Figure 4 and Figure 5. As suggested by Potyonndy and Cundall [17], for rocky materials, $u$ can be considered to have a nonzero value of 0.5. To investigate the effect of $u$ on the macroscopic response, the variation range of $u$ was set to 0.25–0.75.

3.3. Selection of Macroparameters

In numerical calculations, the mechanical parameters typically used in rock materials include compressive strength ${\sigma }_c$, elastic modulus $E$, and tensile strength ${\sigma }_t$. Therefore, compressive strength ${\sigma }_c$, elastic modulus $E$, and Brazilian splitting strength ${\sigma }_t$ were selected as macroparameters in this study, which were also the physical–mechanical parameters commonly used in finite element numerical calculations.

3.4. Discussion

A Gaussian process is completely specified by its mean function and covariance function, while there are many kinds of covariance functions, selecting different covariance functions will get a different response surface, which will affect the calculation accuracy of the results; so, the user needs to select the appropriate covariance function according to the existing experience. Meanwhile, sampling methods, such as Uniform Design (UD), Latin Hypercube Sampling (LHS), Hammersley Sequence Sampling (HSS), etc., are needed in the process of selecting the training samples, but different sampling methods will get different sample sets, which will also have some impact on the results.

4. Parameter Calibration Process

According to the principle of the GP response surface methodology, a training dataset, including inputs $X$ and outputs $\mathbf{y}$, are required to construct the response surface. The generation of an appropriate training set is the key to constructing the response surface. Therefore, it is necessary to design numerical and physical experiments to generate a training dataset. The purpose of the experimental design was to obtain as much information as possible with as few trials as possible. The choice of the experimental scheme directly affects the fitting accuracy of the response surface and the cost of the test. If the number of tests is too small, it is difficult to reflect the response relationship between the input parameters and output results. If the number of tests is too large, it improves the accuracy but also increases the test cost. Moreover, an excessive number of tests is not necessary [32]. The commonly used experimental design methods are UD, Latin Hypercube Sampling, and Hammersley Sequence Sampling [29]. In this study, the UD method was selected to construct a GP response surface based on the balance between the prediction accuracy and the experimental cost.

4.1. Generation of Training Dataset for Microparameters

The selection of input and output parameters needs to be based on experience to select the most sensitive parameters to the calculation results and discard the insensitive parameters as much as possible to reduce the calculation dimension and improve the calculation efficiency. Meanwhile, we need as many sample sets as possible, and the more sample sets we have, the higher the accuracy of the Gaussian response surface.

The training set $X$ needs to be generated first according to the range of microparameters and subsequently substituted into the numerical model (PFC) to generate the training set $\mathbf{y}$. In the experimental design, the five parameters, ${\overline{E}}_c$, ${\overline{k}}_n/{\overline{k}}_s$, ${\overline{\sigma }}_c$, ${\overline{\tau }}_c$, and $u$ were used as input parameters for $X$. The three parameters, ${\sigma }_c$, $E$ and ${\sigma }_t$, obtained by numerical calculations were used as the output response results $\mathbf{y}$.

According to the range of microparameters, 30 sets of training samples were generated using the UD method. The results (input parameters: $X$) are presented in

Table 2. Training dataset.

No.	Input Parameters: $\mathit{X}$ (UD Results)					Output Parameters: $\mathbf{y}$ (Numerical Calculation Results)
	${\overline{\mathit{E}}}_{\mathit{c}}$/GPa	${\overline{\mathit{k}}}_{\mathit{n}}/{\overline{\mathit{k}}}_{\mathit{s}}$	${\overline{\mathit{\sigma }}}_{\mathit{c}}$/MPa	${\overline{\mathit{\tau }}}_{\mathit{c}}$/MPa	$\mathit{u}$	$\mathit{E}$/GPa	${\mathit{\sigma }}_{\mathit{c}}$/MPa	${\mathit{\sigma }}_{\mathit{t}}$/MPa
1	48.09	1.001	33.74	33.74	0.55	70.24	55.47	4.91
2	48.09	1.183	24.10	43.38	0.75	71.82	56.83	4.13
3	75.57	1.001	33.74	62.66	0.55	126.61	81.15	4.51
4	103.05	1.183	72.30	43.38	0.75	201.37	79.33	4.82
5	48.09	0.637	33.74	53.02	0.35	68.36	68.59	5.72
6	89.31	0.455	72.30	62.66	0.55	149.98	100.11	4.91
7	103.05	1.183	33.74	72.30	0.25	201.37	79.29	4.81
8	61.83	0.637	53.02	33.74	0.55	61.44	56.06	4.62
9	61.83	1.183	43.38	24.10	0.35	74.93	45.05	4.13
10	48.09	0.637	53.02	72.30	0.75	68.36	107.37	5.71
11	75.57	0.455	33.74	24.10	0.75	67.51	38.31	4.73
12	89.31	1.365	53.02	53.02	0.55	131.96	90.60	5.31
13	103.05	0.819	43.38	43.38	0.45	153.37	71.44	4.87
14	34.35	0.455	62.66	43.38	0.45	43.28	70.68	5.01
15	103.05	0.637	24.10	53.02	0.65	159.84	61.31	8.02
16	34.35	0.819	72.30	24.10	0.65	48.75	43.35	3.81
17	34.35	1.365	43.38	72.30	0.65	66.07	88.60	3.72
18	89.31	1.001	53.02	33.74	0.25	131.20	59.50	4.71
19	48.09	1.365	72.30	53.02	0.25	82.30	96.79	4.11
20	34.35	1.001	53.02	62.66	0.35	54.79	100.42	5.21
21	75.57	0.455	43.38	43.38	0.25	67.51	66.88	4.72
22	61.83	1.183	62.66	62.66	0.45	74.93	110.14	4.14
23	89.31	0.819	43.38	62.66	0.75	106.06	93.50	4.18
24	75.57	0.819	72.30	72.30	0.35	133.63	112.27	4.33
25	61.83	1.001	62.66	53.02	0.65	77.88	91.54	4.45
26	75.57	1.365	62.66	33.74	0.65	116.83	64.21	4.61
27	89.31	1.365	24.10	24.10	0.45	131.96	43.19	5.31
28	61.83	0.455	24.10	72.30	0.45	67.50	69.92	4.72
29	34.35	0.819	24.10	33.74	0.25	48.75	51.10	3.82
30	103.05	0.637	62.66	24.10	0.35	158.48	40.55	11.41

. Thirty sets of numerical calculations were required to generate the training set $\mathbf{y}$.

4.2. Generation of Training Dataset for Macroparameters

Numerical calculations, including uniaxial compression and Brazilian splitting, were performed using the microparameters above to obtain the training set (output parameters: $\mathbf{y}$). The sizes of the rock samples were 50 × 100 mm (diameter × height) for the uniaxial compression simulation and 50 mm (diameter) for the Brazilian splitting simulation. The minimum radius of the particle was $R_{\min }=2\mathrm{mm}$, and the radius ratio of the particle was $R_{\max }/R_{\min }=1.66$. The numerical model is shown in Figure 6.

Based on the above model, 30 sets of numerical simulations, including uniaxial compression and Brazilian splitting, were performed using the input parameters ($X$) in Table 2. The stress and strain values in each direction during loading were recorded, stress–strain curves of the numerical calculations were obtained, and uniaxial compression strengths ${\sigma }_c$, elastic modulus $E$, and Brazilian splitting strength ${\sigma }_t$ were calculated as output results ($\mathbf{y}$) (see Table 2).

4.3. Establishment of GP Response Surface and Calibration Process of Microparameters

The calculation for the GP model in Equations (1)–(6) was prepared as a toolkit embedded in MATLAB, and the calculation process was completed in MATLAB. The training set must be transformed into a matrix format recognisable by the MATLAB programme, as shown in Equations (7) and (8).

Training set (inputs):

$$\mathbf{x}=\left[\begin{array}{ccccc}{\overline{E}}_c& {\overline{k}}_n/{\overline{k}}_s& {\overline{\sigma }}_c& {\overline{\tau }}_c& u\end{array}\right],$$

(7)

where the first, second, third, fourth, and fifth columns are the modulus of parallel bonds, the ratio of normal to tangential stiffness of parallel bonds, the average value of normal strength of parallel bonds, the average value of tangential strength of parallel bonds, and the friction coefficient of the particles, respectively. There are 30 rows in total.

Training set (outputs):

$$\mathbf{y}=\left[\begin{array}{ccc}E& {\sigma }_c& {\sigma }_t\end{array}\right],$$

(8)

where the first, second, and third columns are the elasticity modulus, uniaxial compression strengths, and Brazilian splitting strength, respectively.

Laboratory uniaxial compression and Brazilian splitting experiments of marble were carried out, and the corresponding macroparameters were obtained, as shown in

Table 3. Macroparameters of marble.

$\mathit{E}$/GPa	${\mathit{\sigma }}_{\mathit{c}}$/MPa	${\mathit{\sigma }}_{\mathit{t}}$/MPa
66.1	63.9	3.8

.

The macroparameters were used as outputs ${\mathbf{y}}^{\ast }$ (Equation (9)) and were put into the GP response surface model to obtain the predictions ${\mathbf{x}}^{\ast }$ (Equation (10)). The calibrated microparameters are listed in

Table 4. Calibration values of microparameters.

${\overline{\mathit{E}}}_{\mathit{c}}$/GPa	${\overline{\mathit{k}}}_{\mathit{n}}/{\overline{\mathit{k}}}_{\mathit{s}}$	${\overline{\mathit{\sigma }}}_{\mathit{c}}$/MPa	${\overline{\mathit{\tau }}}_{\mathit{c}}$/MPa	$\mathit{u}$
45.72	0.94	44.83	43.37	0.497

.

Outputs (macroparameters):

$${\mathbf{y}}^{\ast }=\left[\begin{array}{ccc}66.1& 63.9& 3.8\end{array}\right].$$

(9)

Inputs (microparameters):

$${\mathbf{x}}^{\ast }=\left[45.72\hspace{1em}\hspace{1em}0.94\hspace{1em}\hspace{1em}44.83\hspace{1em}\hspace{1em}43.37\hspace{1em}\hspace{1em}0.497\right].$$

(10)

4.4. Experimental Validation of Calibrated Microparameters

To verify the accuracy of the calibrated microparameters, the values in Table 4 were used to establish a numerical model for uniaxial compression and Brazilian splitting numerical calculations, and the calculated results were compared with those of laboratory experiments.

The stress–strain curves obtained from the numerical simulations and laboratory experiments are shown in Figure 7, and the calculated macroparameters are listed in

Table 5. Calculated and experimental values of macroparameters for marble.

	$\mathit{E}$/GPa	${\mathit{\sigma }}_{\mathit{c}}$/MPa	${\mathit{\sigma }}_{\mathit{t}}$/MPa
Experimental values	66.1	63.9	3.8
Calculated values	62.8	68.9	3.9
Relative error	5.3%	−7.8%	−2.6%

. In the laboratory experimental results, the stress–strain curves showed a downward trend (nonlinear growth trend), indicating that the rock samples underwent a compression-density process. However, the stress–strain curves in the numerical calculations could not characterise this phenomenon, and the calculated values were slightly higher than the experimental values. In the numerical simulations and laboratory experiments, the relative errors of the GP response surface calculated values and laboratory test values were 5.3% for the modulus of elasticity, −7.8% for uniaxial compression strengths, and −2.6% for tensile strength. The specimens exhibited splitting failure characteristics in the simulation of uniaxial compression and cracking failure along the middle of the rock samples in the Brazilian splitting simulation, which was maintained at the same level as the failure characteristics of the laboratory experiments.

The uniaxial compressive strength from the experiment is 63.9 MPa while that of the numerical calculation is 68.9 MPa. The calculated values were slightly higher than the experimental values, and the relative error between them is −7.8%, which shows that the calibrated microparameters can reflect the macro failure strength. But, the calculation accuracy needs to be further improved. The numerical calculation only considers a limited number of factors, and the calculated constitutive model fails to consider the nonlinear characteristics of the rock. But, in the actual stress environment, there may be a combination of multiple factors. Therefore, more factors may need to be considered in the numerical calculation.

5. Conclusions

In this study, a novel method for the calibration of microparameters based on the discrete element method and the GP response surface methodology was proposed. A GP model was developed to establish the GP response surface using a training dataset generated by the UD method and numerical simulations. The microparameters were calibrated using the macroparameters obtained from uniaxial compression and Brazilian splitting tests, and the accuracy of the numerical model was verified using the laboratory test data.

(1) The sensitivity parameters (i.e., ${\sigma }_c$, $E$, $v$, and ${\sigma }_t$) were screened, and the range of the microparameters was initially determined according to their effect on the macroparameters using the control-variable method;

(2) The training set $X$ of the microparameters was generated according to the UD method, and a numerical model was established to generate the training set $\mathbf{y}$ of the macroparameters. The response relationship of rock macro–micro parameters in the PFC was established using the GP response surface methodology, and the microparameters (${\mathbf{x}}^{\ast }$) were calibrated according to the macroparameters (${\mathbf{y}}^{\ast }$) determined by laboratory uniaxial and Brazilian splitting experiments;

(3) Numerical models were established based on the calibrated microparameters, and numerical calculations were performed to compare the stress–strain curves and failure morphology of the rock specimens with laboratory tests under uniaxial compression and Brazilian splitting. The results showed that the rock specimens had the same failure strength and morphology in laboratory experiments and numerical calculations. The relative errors of the GP response surface and laboratory test values were 5.3% for the modulus of elasticity, −7.8% for compressive strength, and −2.6% for tensile strength, which verified the accuracy of the macro–micro parameter relationship established by the GP response surface methodology in the PFC for rock materials;

(4) GP response surfaces can consider the interaction of multiple factors, which is suitable for establishing nonlinear response surfaces between the macro–micro parameters. This method can be extended to other contact and numerical models.