# Calibrating the Micromechanical Parameters of the PFC2D(3D) Models Using the Improved Simulated Annealing Algorithm

Authors: Min Wang; Ping Cao
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

PFC2D(3D) is commercial software, which is commonly used to model the crack initiation of rock and rock-like materials. For the PFC2D(3D) numerical simulation, a proper set of microparameters need to be determined before the numerical simulation. To obtain a proper set of microparameters for PFC2D(3D) model based on the macroparameters obtained from physical experiments, a novel technique has been carried out in this paper. The improved simulated annealing algorithm was employed to calibrate the microparameters of the numerical simulation model of PFC2D(3D). A Python script completely controls the calibration process, which can terminate automatically based on a termination criterion. The microparameter calibration process is not based on establishing the relationship between microparameters and macroparameters; instead, the microparameters are calibrated according to the improved simulated annealing algorithm. By using the proposed approach, the microparameters of both the contact-bond model and parallel-bond model in PFC2D(3D) can be determined. To verify the validity of calibrating the microparameters of PFC2D(3D) via the improved simulated annealing algorithm, some examples were selected from the literature. The corresponding numerical simulations were performed, and the numerical simulation results indicated that the proposed method is reliable for calibrating the microparameters of PFC2D(3D) model.


---

## Body

## 1. Introduction

The discrete element method (DEM) was firstly proposed by Cundall in 1971 [1]. The DEM was investigated further over the following years [2–7]. It has been widely employed in modeling the damage and nonlinear behaviors of materials. The particle flow code in 2 or 3 dimensions (PFC2D(3D)) models the movement and interaction of circular (2D) or spherical (3D) particles using the DEM. It has many advantages [8–10]: It is efficient, as contact detection between circular objects is much simpler than contact detection between angular objects, it is possible for the blocks to break (because they are composed of bonded particles), among others. Thanks to these advantages of the PFC2D(3D), PFC2D(3D) is extensively utilized to solve rock mechanics and rock engineering problems [11–23]. However, the software also has its disadvantages: it requires calibration. In other words, some microparameters must be specified to result in a material with desired macroparameters such as the uniaxial compressive strength (UCS), Young’s modulus, Poisson’s ratio, and tensile strength.The relationship between the microparameters and the macroparameters is difficult to quantify and the microparameters cannot be directly determined according to the macroparameters obtained from the physical experiments. In practice, however, the microparameters of a numerical simulation model in PFC2D(3D) can be calibrated based on the macroparameters determined by the physical experiments, for example, UCS, Poisson’s ratio, Young’s modulus, and tensile strength. According to the difference between the macroparameters obtained from the physical experiments and the numerical simulation, the microparameters are calibrated until the macroparameters obtained from the numerical simulation are sufficiently closed to those from physical experiments. This calibration procedure is called the “trial and error” method [24]. While the drawbacks of the “trial and error” method are obvious, on the one hand, the calibration procedure is subjective, as it depends on the experiments and the experimenter; if the experimenter is not experienced at calibrating the microparameters, then the calibration process will take a very long time. On the other hand, calibrating the microparameters means modifying the microparameters in the command flow  .txt file of the numerical simulations for each step of calibration, which would be hard work. In summary, calibrating the microparameters via the “trial and error” method does not illustrate how to calibrate the microparameters of PFC2D(3D) specifically, as it depends on the experiments or the experimenter. Despite these drawbacks of the “trial and error” approach, the method has been commonly adopted by many researchers [25–33] primarily because there is no other better way to determine the microparameters of PFC2D(3D) models.To avoid the subjectivity in the process of calibrating the microparameters, Yoon [34] carried out a new approach for calibrating the microparameters of contact-bond models in PFC2D. The relationships between microparameters and UCS, Young’s modulus, and Poisson’s ratio were constructed. By combining the numerical simulation results, Plackett-Burman designed a central composite method. The optimum set of microparameters were determined, and the macroparameters obtained from the numerical simulation results were in good agreement with the laboratory results, whereas the interaction of different microparameters is not considered in the method presented by Yoon [34], and the approach can only be used in the contact-bond model of PFC2D. Additionally, the method can only determine the microparameters of rock materials with their physical properties falling within the following ranges: UCS (40–170 MPa), Young’s modulus (20–50 GPa), and Poisson’s ratio (0.19–0.25). In summary, the approach proposed by Yoon [34] has a limited application. Tawadrous et al. [35] conducted a large number of PFC3D numerical simulations. By combining the numerical simulation results and artificial neural networks, the microparameters of parallel-bond of PFC3D can be predicted. However, the numbers of determined microparameters are limited; namely, the parallel-bond and particle elastic modulus, normal-to-shear stiffness ratio, and parallel-bond strength can be determined by utilizing this method. Moreover, the key for the success of artificial neural networks is sufficient training data [36–38], which implies that a large number of numerical simulations should be conducted using this method. Its applicability to the other models of PFC2D(3D) also requires further investigation.According to the references given above, the difficulty of screening out a proper set of microparameters can be classified into three categories:(1) The subjectivity of calibrating the microparameters: the calibration process should be more objective and should not depend on the experiments or the experimenter; (2) the hard work of calibrating the microparameters: during the calibration process, the microparameters must be changed by hand for each calibration step, which is time-consuming and tedious; (3) the limited use of the microparameters calibration method: the calibration method should be applied to different kinds of numerical models in both PFC2D and PFC3D software.For convenience of singling out a proper set of microparameters of PFC2D(3D) based on some basic experimental macroparameters (UCS, Young’s modules, Poisson’s ratio, tensile strength, etc.), a new approach for calibrating microparameters is proposed in this paper. The method is based on the improved simulated annealing algorithm. In addition, Python scripts were developed to accomplish the calibration process automatically. The main merit of the proposed method is decreasing the difficulty of calibrating microparameters in calibration process. Additionally, it avoids the subjectivity in calibrating microparameters, and it can be applied to calibrate the microparameters of contact-bond materials and parallel-bond models in both PFC2D and PFC3D. Additionally, the numbers of microparameters and macroparameters can be increased or decreased according to the specific circumstances in the presented method, which is quite flexible in practical use.


## 2. Calibration Process via the Improved Simulated Annealing Algorithm

### 2.1. Introduction of the Simulated Annealing Algorithm

The simulated annealing algorithm was a stochastic search method that was first carried out by Metropolis et al. [46], and, then, the simulated annealing algorithm was successfully applied to solving the optimization problems by Kirkpatrick et al. [47]. The simulated annealing algorithm is analogous to the annealing process of materials. Boltzmann [48] reckoned that if a system was in thermal equilibrium at a temperatureT, the probability PT(s) of the system being in a given state s could be expressed as follows:(1)PTs=exp⁡-Es/kT∑s∈Sexp⁡-Es/kT,where E(s) is the energy of the state s, k is the Boltzmann constant (cooling coefficient), and S is the set of all the possible states. However, (1) does not give any information on how the material reaches thermal equilibrium at a given temperature. Metropolis et al. [46] proposed an algorithm that simulated the process described by Boltzmann: when the system is in the original state sold with original energy E(sold), a random neighborhood state snew is selected, which leads to a new energy E(snew). Based on the Metropolis criterion, if E(snew)≤E(sold), then the new state is accepted. If E(snew)>E(sold), then the probability of accepting the new state can be written as follows:(2)Paccept=exp⁡-Esold-EsnewkT.To reach the thermal equilibrium completely, the process will be repeatedMarkov times at each temperature, andMarkov is thus called the Markov chain [49]. For better understanding of the simulated annealing algorithm, its flowchart is illustrated in Figure 1.Figure 1
Flowchart of the simulated annealing algorithm.In recent years, the simulated annealing algorithm has been widely applied in the field of optimization [50–58]. The process of calibrating microparameters can be described as minimizing the difference between the macroparameters obtained from the numerical simulations and the physical experiments, which is a typical optimization problem: the calibration process is to find the minimum value of the difference between the macroparameters obtained from the numerical simulation and the physical experiments, and the corresponding set of microparameters are the proper one. Accordingly, the calibration process is transformed into an optimization problem. However, the simulated annealing algorithm should be improved to minimize the total iterations, which will be given in the following section.


### 2.2. Disadvantages of the Simulated Annealing Algorithm

The simulated annealing algorithm has a drawback: the computation times are large [59, 60], which is mainly caused by the characteristics of the simulated annealing algorithm itself. There are four conditions for the success of the simulated annealing algorithm which are as follows [61–65]: (1) The initial temperature must be high enough; (2) the cooling speed for the temperature must be slow enough; (3) for each temperature, the length of the Markov chain must be long enough; (4) the final temperature must be low enough. These four conditions increase the total computation times directly or indirectly, which results in the slow convergence speed. Thereafter, it is necessary to improve the simulated annealing algorithm so as to decrease the total iterations. To improve the simulated annealing algorithm, some procedures of the simulated annealing algorithm were changed, which will be explained explicitly in the following part.


### 2.3. Improved Simulated Annealing Algorithm

In this paper, the simulated annealing algorithm would be improved from two aspects: the disturbance method and the cooling method.


#### 2.3.1. The Change of the Cooling Temperature Method

For the original simulated annealing algorithm, the cooling process can be expressed as follows:(3)T=kT0,where T0 is the initial temperature, k is the cooling efficient, and T is the cooling temperature.According to Ingber [66–68], the convergence speed of the simulated annealing algorithm will increase when the following cooling strategy is adopted:(4)T=T0exp⁡-Ci1/N,where i is the iteration time, C is a constant, and N is the number of input parameters in the simulated annealing algorithm.In practice use, (4) can be simplified as follows:(5)T=T0βi1/N,where β is a constant, which is falling within 0.7<β<1. In this paper, β is set to 0.99.


#### 2.3.2. The Change of the Disturbance Method for the Variables

In the traditional simulated annealing algorithm, the disturbance for the variable can be denoted as follows:(6)vj-new=vj-old+Bvj-upper-Bvj-lowerunormal-0.5,where vj-new is the new generated value for the variable vj, vj-old is the original value of vj, Bvj-upper is the upper bound of the variable vj, Bvj-lower is the lower bound of the variable vj, and unormal is a random number between 0 and 1, which is subjected to the normal distribution.Based on the study by Ingber [66, 69, 70], the Cauchy random distribution function for the disturbance is conducive for the convergence of the simulated annealing algorithm:(7)vj-new=vj-old+Bvj-upper-Bvj-loweruCauchy-0.5,where uCauchy is the random number between 0 and 1 and the generated value uCauchy is subjected to the Cauchy distribution.


### 2.4. Calibrating Microparameters via the Improved Simulated Annealing Algorithm

Before using the improved simulated annealing algorithm for calibrating the microparameters, some changes are required. After the PFC2D(3D) numerical simulations based on the set of microparameters, the corresponding macroparameters can be obtained. In the calibration process, the macroparameters determined by numerical simulations are defined as the state of the system. To implement the microparameters calibration process via the improved simulated annealing algorithm, the energy function (objective function)E(MA1-PFC,MA2-PFC,…,MAi-PFC) calculated by the macroparameters of the numerical simulation is denoted as follows:(8)EMA1-PFC,MA2-PFC,…,MAi-PFC=max⁡absMA1-experiment-MA1-PFCMA1-experiment,absMA2-experiment-MA2-PFCMA2-experiment,…,absMAi-experiment-MAi-PFCMAi-experiment.where MAi-experiment is the set of macroparameters obtained from the physical experiments and MAi-PFC is the set of macroparameters determined by the numerical simulations. The objective function in (8) indicates that the value of EMA1-PFC,MA2-PFC,…,MAi-PFC is the maximum difference between the macroparameters obtained from the numerical simulation and the physical experiments. The macroparameters MA1-PFC,MA2-PFC,…,MAi-PFC of the numerical simulation are determined by performing the numerical simulations by using the microparameters (MI1,MI2,…,MIj). It should be noted that the numbers of macroparameters and microparameters can be increased or decreased according to the specific circumstances, which is one of the advantages of the proposed technique.To search for the microparameters completely within the determined range, the length of the Markov chain is utilized. The length of the Markov chain represents the iteration times at each temperature. With increasing length of the Markov chain, the corresponding iteration time will increase as well; in this paper, the length of the Markov chain is set to 100.Another important factor is the bounding of the microparameters. The new microparameters should be generated within a certain range during the calibration, or the numerical simulations of PFC2D(3D) cannot proceed successfully. For example, all the microparameters in PFC2D(3D) should be larger than 0, or the numerical simulations cannot run successfully.In this paper, the termination criterion of the calibration process isEMA1-PFC,MA2-PFC,…,MAi-PFC<10%, meaning that the maximum error between the macroparameters obtained from the numerical simulations and the physical experiments is less than 10%.By combining the improved simulated annealing algorithm and some basic parameters for the calibration process, the procedure of calibrating microparameters via the improved simulated annealing algorithm can be described step-by-step as follows.Step 1.
The initial temperatureT, the cooling efficient k, the length of Markov chainMarkov, the numbers of the microparametersN, and the range of microparameters are initialized, and the initial iteration timei is set to 1. A set of microparameters MI1,MI2,…,MIjold is randomly generated; the microparameters should be within the determined range. By using the microparameters MI1,MI2,…,MIjold, PFC2D or PFC3D numerical simulations are conducted to obtain the macroparameters MA1-PFC,MA2-PFC,…,MAi-PFCold. EMA1-PFC,MA2-PFC,…,MAi-PFCold is estimated based on (8).Step 2.
IfEMA1-PFC,MA2-PFC,…,MAi-PFCold<10%, then go to Step 10.Step 3.
M = 0.Step 4.
A new set of random neighborhood microparametersMI1,MI2,…,MIjnew near the microparameters MI1,MI2,…,MIjold is generated; the new generated microparameters should be within the determined range. Taking the microparameter MIj as example, the new generated microparameter MIj-new can be expressed as follows:(9)MIj-new=MIj-old+Bj-upper-Bj-lower∗uCauchy-0.5,Bj-lower<MIj-new<Bj-upper,where MIj-new is the newly generated microparameter, MIj-old is the original microparameter, Bj-upper is the upper bound of the microparameterMIj, Bj-lower is the lower bound of the microparameter MIj, and uCauchy represents a random value between 0 and 1; besides, the generated value is subjected to the Cauchy distribution. The macroparameters MA1-PFC,MA2-PFC,…,MAi-PFCnew are determined by the PFC2D or PFC3D numerical simulations by using the microparameters MI1,MI2,…,MIjnew. EMA1-PFC,MA2-PFC,…,MAi-PFCnew is calculated based on (9).Step 5.
IfEMA1-PFC,MA2-PFC,…,MAi-PFCnew<EMA1-PFC,MA2-PFC,…,MAi-PFCold, then the new microparameters MI1,MI2,…,MIjnew and EMA1-PFC,MA2-PFC,…,MAi-PFCnew are accepted, while if EMA1-PFC,MA2-PFC,…,MAi-PFCnew≥EMA1-PFC,MA2-PFC,…,MAi-PFCold, then the probability of accepting the new microparameters MI1,MI2,…,MIjnew and EMI1,MI2,…,MIjnew can be expressed as follows:(10)Paccept=exp⁡-EMA1-PFC,MA2-PFC,…,MAi-PFCold-EMA1-PFC,MA2-PFC,…,MAi-PFCnewT.Step 6.
IfEMA1-PFC,MA2-PFC,…,MAi-PFCnew<10%, then go to Step 10.Step 7.
i
=
i
+
1.Step 8.
M
=
M
+
1; if M≤Markov, then go to Step 4.Step 9.
T
=
0.9
9
i
1
/
N
T; go to Step 3.Step 10.
End.Because the microparameters need to be changed and the PFC2D(3D) software needs to be manipulated many times for each step of calibration, the process of calibrating microparameters is accomplished using Python script. Python (https://www.python.org) is a scripting language that is very easy to use and widely applied in many fields [71–76].The numerical simulation mentioned in this section includes four parts:(1) The command flow  .txt file based on the microparameters is written; (2) the particle assembly is carried out; (3) the numerical simulation is performed to obtain the macroparameters; (4) the output of the macroparameters is produced.


#### 2.4.1. Writing the Command Flow  .txt File Based on the Microparameters

In practice, the command flow is frequently complied in a  .txt file, which will be convenient for the PFC2D(3D) numerical simulations. Due to the change in the microparameters for each step of calibration, the command flow  .txt file is written by the Python scripts, and the microparameters for each calibration are changed according to the improved simulated annealing algorithm.


#### 2.4.2. The Particle Assembly

Before the numerical simulation, a PFC2D(3D) assembly is generated, and the process involves particle generation, packing the particles, isotropic stress installation (stress initialization), floating particle elimination, and bond installation, which is described in detail by Itasca [9, 77] and will not be repeated in this paper.


#### 2.4.3. The Numerical Simulation to Obtain the Macroparameters

To obtain the macroparameters, the corresponding numerical simulations should be carried out. The UCS, Young’s modulus, and Poisson’s ratio can be obtained by the uniaxial compression test. Meanwhile, the tensile strength can be determined via the Brazilian disc test. Other numerical simulations can also be utilized to obtain the corresponding macroparameters. In practice, the numbers of macroparameters and microparameters can be increased or decreased according to the specific circumstances. The corresponding command flow  .txt file should be changed accordingly.


#### 2.4.4. The Output of the Macroparameters

The aim of the PFC2D(3D) numerical simulations is to obtain the numerical macroparameters. Thereafter, the macroparameters will be output in the log file when the numerical simulations terminate [9, 10].


## 2.1. Introduction of the Simulated Annealing Algorithm

The simulated annealing algorithm was a stochastic search method that was first carried out by Metropolis et al. [46], and, then, the simulated annealing algorithm was successfully applied to solving the optimization problems by Kirkpatrick et al. [47]. The simulated annealing algorithm is analogous to the annealing process of materials. Boltzmann [48] reckoned that if a system was in thermal equilibrium at a temperatureT, the probability PT(s) of the system being in a given state s could be expressed as follows:(1)PTs=exp⁡-Es/kT∑s∈Sexp⁡-Es/kT,where E(s) is the energy of the state s, k is the Boltzmann constant (cooling coefficient), and S is the set of all the possible states. However, (1) does not give any information on how the material reaches thermal equilibrium at a given temperature. Metropolis et al. [46] proposed an algorithm that simulated the process described by Boltzmann: when the system is in the original state sold with original energy E(sold), a random neighborhood state snew is selected, which leads to a new energy E(snew). Based on the Metropolis criterion, if E(snew)≤E(sold), then the new state is accepted. If E(snew)>E(sold), then the probability of accepting the new state can be written as follows:(2)Paccept=exp⁡-Esold-EsnewkT.To reach the thermal equilibrium completely, the process will be repeatedMarkov times at each temperature, andMarkov is thus called the Markov chain [49]. For better understanding of the simulated annealing algorithm, its flowchart is illustrated in Figure 1.Figure 1
Flowchart of the simulated annealing algorithm.In recent years, the simulated annealing algorithm has been widely applied in the field of optimization [50–58]. The process of calibrating microparameters can be described as minimizing the difference between the macroparameters obtained from the numerical simulations and the physical experiments, which is a typical optimization problem: the calibration process is to find the minimum value of the difference between the macroparameters obtained from the numerical simulation and the physical experiments, and the corresponding set of microparameters are the proper one. Accordingly, the calibration process is transformed into an optimization problem. However, the simulated annealing algorithm should be improved to minimize the total iterations, which will be given in the following section.


## 2.2. Disadvantages of the Simulated Annealing Algorithm

The simulated annealing algorithm has a drawback: the computation times are large [59, 60], which is mainly caused by the characteristics of the simulated annealing algorithm itself. There are four conditions for the success of the simulated annealing algorithm which are as follows [61–65]: (1) The initial temperature must be high enough; (2) the cooling speed for the temperature must be slow enough; (3) for each temperature, the length of the Markov chain must be long enough; (4) the final temperature must be low enough. These four conditions increase the total computation times directly or indirectly, which results in the slow convergence speed. Thereafter, it is necessary to improve the simulated annealing algorithm so as to decrease the total iterations. To improve the simulated annealing algorithm, some procedures of the simulated annealing algorithm were changed, which will be explained explicitly in the following part.


## 2.3. Improved Simulated Annealing Algorithm

In this paper, the simulated annealing algorithm would be improved from two aspects: the disturbance method and the cooling method.


### 2.3.1. The Change of the Cooling Temperature Method

For the original simulated annealing algorithm, the cooling process can be expressed as follows:(3)T=kT0,where T0 is the initial temperature, k is the cooling efficient, and T is the cooling temperature.According to Ingber [66–68], the convergence speed of the simulated annealing algorithm will increase when the following cooling strategy is adopted:(4)T=T0exp⁡-Ci1/N,where i is the iteration time, C is a constant, and N is the number of input parameters in the simulated annealing algorithm.In practice use, (4) can be simplified as follows:(5)T=T0βi1/N,where β is a constant, which is falling within 0.7<β<1. In this paper, β is set to 0.99.


### 2.3.2. The Change of the Disturbance Method for the Variables

In the traditional simulated annealing algorithm, the disturbance for the variable can be denoted as follows:(6)vj-new=vj-old+Bvj-upper-Bvj-lowerunormal-0.5,where vj-new is the new generated value for the variable vj, vj-old is the original value of vj, Bvj-upper is the upper bound of the variable vj, Bvj-lower is the lower bound of the variable vj, and unormal is a random number between 0 and 1, which is subjected to the normal distribution.Based on the study by Ingber [66, 69, 70], the Cauchy random distribution function for the disturbance is conducive for the convergence of the simulated annealing algorithm:(7)vj-new=vj-old+Bvj-upper-Bvj-loweruCauchy-0.5,where uCauchy is the random number between 0 and 1 and the generated value uCauchy is subjected to the Cauchy distribution.


## 2.3.1. The Change of the Cooling Temperature Method

For the original simulated annealing algorithm, the cooling process can be expressed as follows:(3)T=kT0,where T0 is the initial temperature, k is the cooling efficient, and T is the cooling temperature.According to Ingber [66–68], the convergence speed of the simulated annealing algorithm will increase when the following cooling strategy is adopted:(4)T=T0exp⁡-Ci1/N,where i is the iteration time, C is a constant, and N is the number of input parameters in the simulated annealing algorithm.In practice use, (4) can be simplified as follows:(5)T=T0βi1/N,where β is a constant, which is falling within 0.7<β<1. In this paper, β is set to 0.99.


## 2.3.2. The Change of the Disturbance Method for the Variables

In the traditional simulated annealing algorithm, the disturbance for the variable can be denoted as follows:(6)vj-new=vj-old+Bvj-upper-Bvj-lowerunormal-0.5,where vj-new is the new generated value for the variable vj, vj-old is the original value of vj, Bvj-upper is the upper bound of the variable vj, Bvj-lower is the lower bound of the variable vj, and unormal is a random number between 0 and 1, which is subjected to the normal distribution.Based on the study by Ingber [66, 69, 70], the Cauchy random distribution function for the disturbance is conducive for the convergence of the simulated annealing algorithm:(7)vj-new=vj-old+Bvj-upper-Bvj-loweruCauchy-0.5,where uCauchy is the random number between 0 and 1 and the generated value uCauchy is subjected to the Cauchy distribution.


## 2.4. Calibrating Microparameters via the Improved Simulated Annealing Algorithm

Before using the improved simulated annealing algorithm for calibrating the microparameters, some changes are required. After the PFC2D(3D) numerical simulations based on the set of microparameters, the corresponding macroparameters can be obtained. In the calibration process, the macroparameters determined by numerical simulations are defined as the state of the system. To implement the microparameters calibration process via the improved simulated annealing algorithm, the energy function (objective function)E(MA1-PFC,MA2-PFC,…,MAi-PFC) calculated by the macroparameters of the numerical simulation is denoted as follows:(8)EMA1-PFC,MA2-PFC,…,MAi-PFC=max⁡absMA1-experiment-MA1-PFCMA1-experiment,absMA2-experiment-MA2-PFCMA2-experiment,…,absMAi-experiment-MAi-PFCMAi-experiment.where MAi-experiment is the set of macroparameters obtained from the physical experiments and MAi-PFC is the set of macroparameters determined by the numerical simulations. The objective function in (8) indicates that the value of EMA1-PFC,MA2-PFC,…,MAi-PFC is the maximum difference between the macroparameters obtained from the numerical simulation and the physical experiments. The macroparameters MA1-PFC,MA2-PFC,…,MAi-PFC of the numerical simulation are determined by performing the numerical simulations by using the microparameters (MI1,MI2,…,MIj). It should be noted that the numbers of macroparameters and microparameters can be increased or decreased according to the specific circumstances, which is one of the advantages of the proposed technique.To search for the microparameters completely within the determined range, the length of the Markov chain is utilized. The length of the Markov chain represents the iteration times at each temperature. With increasing length of the Markov chain, the corresponding iteration time will increase as well; in this paper, the length of the Markov chain is set to 100.Another important factor is the bounding of the microparameters. The new microparameters should be generated within a certain range during the calibration, or the numerical simulations of PFC2D(3D) cannot proceed successfully. For example, all the microparameters in PFC2D(3D) should be larger than 0, or the numerical simulations cannot run successfully.In this paper, the termination criterion of the calibration process isEMA1-PFC,MA2-PFC,…,MAi-PFC<10%, meaning that the maximum error between the macroparameters obtained from the numerical simulations and the physical experiments is less than 10%.By combining the improved simulated annealing algorithm and some basic parameters for the calibration process, the procedure of calibrating microparameters via the improved simulated annealing algorithm can be described step-by-step as follows.Step 1.
The initial temperatureT, the cooling efficient k, the length of Markov chainMarkov, the numbers of the microparametersN, and the range of microparameters are initialized, and the initial iteration timei is set to 1. A set of microparameters MI1,MI2,…,MIjold is randomly generated; the microparameters should be within the determined range. By using the microparameters MI1,MI2,…,MIjold, PFC2D or PFC3D numerical simulations are conducted to obtain the macroparameters MA1-PFC,MA2-PFC,…,MAi-PFCold. EMA1-PFC,MA2-PFC,…,MAi-PFCold is estimated based on (8).Step 2.
IfEMA1-PFC,MA2-PFC,…,MAi-PFCold<10%, then go to Step 10.Step 3.
M = 0.Step 4.
A new set of random neighborhood microparametersMI1,MI2,…,MIjnew near the microparameters MI1,MI2,…,MIjold is generated; the new generated microparameters should be within the determined range. Taking the microparameter MIj as example, the new generated microparameter MIj-new can be expressed as follows:(9)MIj-new=MIj-old+Bj-upper-Bj-lower∗uCauchy-0.5,Bj-lower<MIj-new<Bj-upper,where MIj-new is the newly generated microparameter, MIj-old is the original microparameter, Bj-upper is the upper bound of the microparameterMIj, Bj-lower is the lower bound of the microparameter MIj, and uCauchy represents a random value between 0 and 1; besides, the generated value is subjected to the Cauchy distribution. The macroparameters MA1-PFC,MA2-PFC,…,MAi-PFCnew are determined by the PFC2D or PFC3D numerical simulations by using the microparameters MI1,MI2,…,MIjnew. EMA1-PFC,MA2-PFC,…,MAi-PFCnew is calculated based on (9).Step 5.
IfEMA1-PFC,MA2-PFC,…,MAi-PFCnew<EMA1-PFC,MA2-PFC,…,MAi-PFCold, then the new microparameters MI1,MI2,…,MIjnew and EMA1-PFC,MA2-PFC,…,MAi-PFCnew are accepted, while if EMA1-PFC,MA2-PFC,…,MAi-PFCnew≥EMA1-PFC,MA2-PFC,…,MAi-PFCold, then the probability of accepting the new microparameters MI1,MI2,…,MIjnew and EMI1,MI2,…,MIjnew can be expressed as follows:(10)Paccept=exp⁡-EMA1-PFC,MA2-PFC,…,MAi-PFCold-EMA1-PFC,MA2-PFC,…,MAi-PFCnewT.Step 6.
IfEMA1-PFC,MA2-PFC,…,MAi-PFCnew<10%, then go to Step 10.Step 7.
i
=
i
+
1.Step 8.
M
=
M
+
1; if M≤Markov, then go to Step 4.Step 9.
T
=
0.9
9
i
1
/
N
T; go to Step 3.Step 10.
End.Because the microparameters need to be changed and the PFC2D(3D) software needs to be manipulated many times for each step of calibration, the process of calibrating microparameters is accomplished using Python script. Python (https://www.python.org) is a scripting language that is very easy to use and widely applied in many fields [71–76].The numerical simulation mentioned in this section includes four parts:(1) The command flow  .txt file based on the microparameters is written; (2) the particle assembly is carried out; (3) the numerical simulation is performed to obtain the macroparameters; (4) the output of the macroparameters is produced.


### 2.4.1. Writing the Command Flow  .txt File Based on the Microparameters

In practice, the command flow is frequently complied in a  .txt file, which will be convenient for the PFC2D(3D) numerical simulations. Due to the change in the microparameters for each step of calibration, the command flow  .txt file is written by the Python scripts, and the microparameters for each calibration are changed according to the improved simulated annealing algorithm.


### 2.4.2. The Particle Assembly

Before the numerical simulation, a PFC2D(3D) assembly is generated, and the process involves particle generation, packing the particles, isotropic stress installation (stress initialization), floating particle elimination, and bond installation, which is described in detail by Itasca [9, 77] and will not be repeated in this paper.


### 2.4.3. The Numerical Simulation to Obtain the Macroparameters

To obtain the macroparameters, the corresponding numerical simulations should be carried out. The UCS, Young’s modulus, and Poisson’s ratio can be obtained by the uniaxial compression test. Meanwhile, the tensile strength can be determined via the Brazilian disc test. Other numerical simulations can also be utilized to obtain the corresponding macroparameters. In practice, the numbers of macroparameters and microparameters can be increased or decreased according to the specific circumstances. The corresponding command flow  .txt file should be changed accordingly.


### 2.4.4. The Output of the Macroparameters

The aim of the PFC2D(3D) numerical simulations is to obtain the numerical macroparameters. Thereafter, the macroparameters will be output in the log file when the numerical simulations terminate [9, 10].


## 2.4.1. Writing the Command Flow  .txt File Based on the Microparameters

In practice, the command flow is frequently complied in a  .txt file, which will be convenient for the PFC2D(3D) numerical simulations. Due to the change in the microparameters for each step of calibration, the command flow  .txt file is written by the Python scripts, and the microparameters for each calibration are changed according to the improved simulated annealing algorithm.


## 2.4.2. The Particle Assembly

Before the numerical simulation, a PFC2D(3D) assembly is generated, and the process involves particle generation, packing the particles, isotropic stress installation (stress initialization), floating particle elimination, and bond installation, which is described in detail by Itasca [9, 77] and will not be repeated in this paper.


## 2.4.3. The Numerical Simulation to Obtain the Macroparameters

To obtain the macroparameters, the corresponding numerical simulations should be carried out. The UCS, Young’s modulus, and Poisson’s ratio can be obtained by the uniaxial compression test. Meanwhile, the tensile strength can be determined via the Brazilian disc test. Other numerical simulations can also be utilized to obtain the corresponding macroparameters. In practice, the numbers of macroparameters and microparameters can be increased or decreased according to the specific circumstances. The corresponding command flow  .txt file should be changed accordingly.


## 2.4.4. The Output of the Macroparameters

The aim of the PFC2D(3D) numerical simulations is to obtain the numerical macroparameters. Thereafter, the macroparameters will be output in the log file when the numerical simulations terminate [9, 10].


## 3. Examples of Calibrating Microparameters via the Improved Simulated Annealing Algorithm

### 3.1. A Detailed Example

#### 3.1.1. Basic Information about the Numerical Simulations

To illustrate the calibration process of microparameters more specifically, an example was selected. In [78], sandstone was utilized to investigate the Kaiser effect. Moreover, the PFC2D was employed to reproduce the Kaiser effect, the deformation, and the rate of deformation. Meanwhile, the contact-bond model in PFC2D was utilized, and the size of the numerical simulation for the macroparameters was 61.14mm×153.32mm. According to the physical experiments, the UCS and Young’s modulus of the sandstone were 83.1 MPa and 14.3 GPa, respectively.In total, 10 microparameters determine the contact-bonded model in PFC2D. The microparameters and their corresponding descriptions are listed in Table1.Table 1
Microparameters of a contact-bonded material.
Parameter
Description
R
m
i
n
Minimum ball radius [m]
R
m
a
x
/
R
m
i
n
Ball size ratio, uniform distribution []
ρ
Ball density [kg/m3]
E
c
Ball-ball contact modulus [Pa]
k
n
/
k
s
Ball stiffness ratio []
μ
Ball friction coefficient
σ
c (mean)
Contact-bond normal strength, mean [Pa]
σ
c (std. dev.)
Contact-bond normal strength, std. dev. [Pa]
τ
c (mean)
Contact-bond shear strength, mean [Pa]
τ
c (std. dev.)
Contact-bond shear strength, std. dev. [Pa]As listed in Table1, the density of the ball can be treated as the density of the sandstone, which is 2397 kg/m3 [78]. The effect of gravity does not need to be considered, as the specimens are small in the physical experiments and the gravity-induced stress has a negligible effect on the macroparameters [77]. In addition, the microparameters should be in a certain range, which guarantees the success of the numerical simulations. For example, Rmax/Rmin, the ball size ratio, should be greater than 1, or the numerical simulation cannot run successfully. Therefore, it is necessary to determine a reasonable range of microparameters. Table 2 gives the range of the microparameters of the contact-bonded model in the numerical simulations.Table 2
The range of the microparameters of the contact-bonded material.
Parameter
Description
Range
Rmin
Minimum ball radius [m]
0.1e − 3~1e-3
Rmax/Rmin
Ball size ratio, uniform distribution []
1~20
Ec
Ball-ball contact modulus [Pa]
0~100e9
kn/ks
Ball stiffness ratio []
0~10
μ
Ball friction coefficient
0~1
σ
c (mean)
Contact-bond normal strength, mean [Pa]
0~1000e6
σ
c (std. dev.)
Contact-bond normal strength, std. dev. [Pa]
0~100e6
τ
c (mean)
Contact-bond shear strength, mean [Pa]
0~1000e6
τ
c (std. dev.)
Contact-bond shear strength, std. dev. [Pa]
0~100e6It should be noted that the minimum ball radius should be larger than0.1e-3 m, or the error “Error: Fewer balls generated” will occur during the numerical simulations, which will stop the numerical simulation from running. The upper bound of the microparameters in Table 2 can be increased slightly, but this will increase the computation times.


#### 3.1.2. Basic Parameters of the Improved Simulated Annealing Algorithm

Some basic parameters of the improved simulated annealing algorithm need to be determined before the numerical simulations; these are listed in Table3.Table 3
Basic parameters of the improved simulated annealing algorithm.
Parameter
Value
The length of the Markov chain
100
Cooling efficientk
0.95
The initial temperatureT
100
The number of microparametersN
9


#### 3.1.3. Numerical Simulation Platform

In the numerical simulation, the operation system is Linux Mint 18 Sarah (on a personal computer with an Intel Core i5, 3.30 GHz CPU, and 8 GB RAM); the operation system is open source. The Python version is 3.5. Additionally, to implement the improved simulated annealing algorithm and manipulate the PFC2D(3D) software, two extra Python packages were installed,numpy (to accomplish the improved simulated annealing algorithm) andsubprocess (to manipulate the PFC2D(3D) software). The PFC2D version 3.1 in this paper was adapted to run on the Linux system.Wine, an open source software for running the Windows system application, was used for this purpose.


#### 3.1.4. Numerical Simulations Results

Based on the improved simulated annealing algorithm, after 2739 iterations, the termination criterion is satisfied and the numerical simulation is terminated. The UCS and Young’s modulus determined by the numerical simulation are 74.356 MPa and 14.824 GPa, respectively, which is quite close to the physical experiment results (UCS 82.1 MPa and Young’s modulus 14.3 GPa.). The corresponding microparameters are obtained and are listed in Table4.Table 4
Microparameters obtained by the improved simulated annealing algorithm.
Parameter
Value
R
m
i
n
0.746e-3 m
R
m
a
x
/
R
m
i
n
1.567
E
c
20.688 GPa
k
n
/
k
s
1.833
μ
0.590
σ
c (mean)
88.916 MPa
σ
c (std. dev.)
36.950 MPa
τ
c (mean)
83.768 MPa
τ
c (std. dev.)
28.932 MPaThe microparameters can successfully be determined by the improved simulated annealing algorithm. This saves a tremendous amount of hard work, as it is not necessary to change the microparameters and manipulate PFC2D by hand for each step of calibration, which is convenient for use.


### 3.2. More Examples

From the analysis of the detailed example in Section3.1, it is found that the calibration method is independent of the numerical simulation model; that is, once the command flow for numerical simulations is fixed, the calibration process can change the microparameters in the command flow  .txt file for each step of calibration. Thus, it can be used for contact bond and parallel bond in both PFC2D and PFC3D. Recall that the numbers of microparameters and macroparameters can be increased or decreased according to the specific circumstances, in which case the corresponding Python scripts should also be changed correspondingly. To further validate the proposed method, more numerical simulation tests were performed. The numerical simulation results are listed in Table 5.Table 5
Macroparameters obtained by the improved simulated annealing algorithm.
Author
Material species
Numerical model
Macroparameters determined by the numerical simulations (error)
Chang et al.2001 [39]
Kimachi sandstone
PFC2DContact-bond
Experimental results: σc: 48 MPa υ: 0.27E: 6.5 GPa
Numerical results: σc: 44.372 MPa (7.558%)υ: 0.275 (1.851%) E: 6.477 GPa (0.353%)
Chang 2002[40]
Hwangdeung granite
PFC2DContact-bond
Experimental results: σc: 162 MPa υ: 0.28E: 50.7 GPaNumerical results: σc: 171.357 MPa (5.775%)υ: 0.283 (1.071%) E: 48.492 GPa (4.355%)
Hong 2004[41]
Daejeon granite
PFC2DContact-bond
Experimental results: σc: 161.7 MPa υ: 0.22E: 50.9 GPaNumerical results: σc: 159.172 MPa (1.563%)υ: 0.239 (8.636%) E: 46.943 GPa (7.774%)
Wang et al.2016 [12]
Rock cored from the rock mass in the Heishan open pit mine in Chengde city, China
PFC2DContact-bond
Experimental results: σc: 132.3 MPa E: 12.0 GPaNumerical results: σc: 130.329 MPa (1.489%)E: 11.924 GPa (0.633%)
Yang et al.2014 [42]
Red sandstone
PFC2DParallel-bond
Experimental results: σc: 52.98 MPa E: 14.42 GPaε1c: 5.113Numerical results: σc: 57.832 MPa (9.158%)E: 15.490 GPa (7.420%) ε1c: 4.650 (9.055%)
Zhou et al.2016 [43]
Concrete
PFC2DParallel-bond
Experimental results: σc: 23.3 MPa υ: 0.17E: 37 GPaNumerical results: σc: 23.525 MPa (0.965%)υ: 0.159 (6.470%) E: 35.822 GPa (3.183%)
Yang et al.2015 [44]
Sandstone
PFC2DParallel-bond
Experimental results: σc: 190.80 MPa E: 35.64 GPaNumerical results: σc: 173.583 MPa (9.023%)E: 39.035 GPa (9.525%)
Fan et al.2015 [25]
Rock-like materials
PFC3DContact-bond
Experimental results: σc: 8.27 MPa E: 4.04 GPaNumerical results: σc: 8.947 MPa (8.186%)E: 4.228 GPa (4.653%)
Duan et al2015 [45]
Beishan granite
PFC3DParallel-bond
Experimental results: σc: 136.4 MPa υ: 0.266E: 43.167 GPaNumerical results: σc: 127.017 MPa (6.879%)υ: 0.265 (0.375%) E: 42.652 GPa (1.193%)
σ
c: Uniaxial compressive strength; υ: Poisson’s ratio E: Young’s modulus; σt: tensile strength; ε1c: peak axial strain.As listed in Table5, calibrating the microparameters via the improved simulated annealing algorithm is a viable strategy for both parallel bond and contact bond in PFC2D(3D). It has wide practical use in PFC2D(3D) numerical simulations.


## 3.1. A Detailed Example

### 3.1.1. Basic Information about the Numerical Simulations

To illustrate the calibration process of microparameters more specifically, an example was selected. In [78], sandstone was utilized to investigate the Kaiser effect. Moreover, the PFC2D was employed to reproduce the Kaiser effect, the deformation, and the rate of deformation. Meanwhile, the contact-bond model in PFC2D was utilized, and the size of the numerical simulation for the macroparameters was 61.14mm×153.32mm. According to the physical experiments, the UCS and Young’s modulus of the sandstone were 83.1 MPa and 14.3 GPa, respectively.In total, 10 microparameters determine the contact-bonded model in PFC2D. The microparameters and their corresponding descriptions are listed in Table1.Table 1
Microparameters of a contact-bonded material.
Parameter
Description
R
m
i
n
Minimum ball radius [m]
R
m
a
x
/
R
m
i
n
Ball size ratio, uniform distribution []
ρ
Ball density [kg/m3]
E
c
Ball-ball contact modulus [Pa]
k
n
/
k
s
Ball stiffness ratio []
μ
Ball friction coefficient
σ
c (mean)
Contact-bond normal strength, mean [Pa]
σ
c (std. dev.)
Contact-bond normal strength, std. dev. [Pa]
τ
c (mean)
Contact-bond shear strength, mean [Pa]
τ
c (std. dev.)
Contact-bond shear strength, std. dev. [Pa]As listed in Table1, the density of the ball can be treated as the density of the sandstone, which is 2397 kg/m3 [78]. The effect of gravity does not need to be considered, as the specimens are small in the physical experiments and the gravity-induced stress has a negligible effect on the macroparameters [77]. In addition, the microparameters should be in a certain range, which guarantees the success of the numerical simulations. For example, Rmax/Rmin, the ball size ratio, should be greater than 1, or the numerical simulation cannot run successfully. Therefore, it is necessary to determine a reasonable range of microparameters. Table 2 gives the range of the microparameters of the contact-bonded model in the numerical simulations.Table 2
The range of the microparameters of the contact-bonded material.
Parameter
Description
Range
Rmin
Minimum ball radius [m]
0.1e − 3~1e-3
Rmax/Rmin
Ball size ratio, uniform distribution []
1~20
Ec
Ball-ball contact modulus [Pa]
0~100e9
kn/ks
Ball stiffness ratio []
0~10
μ
Ball friction coefficient
0~1
σ
c (mean)
Contact-bond normal strength, mean [Pa]
0~1000e6
σ
c (std. dev.)
Contact-bond normal strength, std. dev. [Pa]
0~100e6
τ
c (mean)
Contact-bond shear strength, mean [Pa]
0~1000e6
τ
c (std. dev.)
Contact-bond shear strength, std. dev. [Pa]
0~100e6It should be noted that the minimum ball radius should be larger than0.1e-3 m, or the error “Error: Fewer balls generated” will occur during the numerical simulations, which will stop the numerical simulation from running. The upper bound of the microparameters in Table 2 can be increased slightly, but this will increase the computation times.


### 3.1.2. Basic Parameters of the Improved Simulated Annealing Algorithm

Some basic parameters of the improved simulated annealing algorithm need to be determined before the numerical simulations; these are listed in Table3.Table 3
Basic parameters of the improved simulated annealing algorithm.
Parameter
Value
The length of the Markov chain
100
Cooling efficientk
0.95
The initial temperatureT
100
The number of microparametersN
9


### 3.1.3. Numerical Simulation Platform

In the numerical simulation, the operation system is Linux Mint 18 Sarah (on a personal computer with an Intel Core i5, 3.30 GHz CPU, and 8 GB RAM); the operation system is open source. The Python version is 3.5. Additionally, to implement the improved simulated annealing algorithm and manipulate the PFC2D(3D) software, two extra Python packages were installed,numpy (to accomplish the improved simulated annealing algorithm) andsubprocess (to manipulate the PFC2D(3D) software). The PFC2D version 3.1 in this paper was adapted to run on the Linux system.Wine, an open source software for running the Windows system application, was used for this purpose.


### 3.1.4. Numerical Simulations Results

Based on the improved simulated annealing algorithm, after 2739 iterations, the termination criterion is satisfied and the numerical simulation is terminated. The UCS and Young’s modulus determined by the numerical simulation are 74.356 MPa and 14.824 GPa, respectively, which is quite close to the physical experiment results (UCS 82.1 MPa and Young’s modulus 14.3 GPa.). The corresponding microparameters are obtained and are listed in Table4.Table 4
Microparameters obtained by the improved simulated annealing algorithm.
Parameter
Value
R
m
i
n
0.746e-3 m
R
m
a
x
/
R
m
i
n
1.567
E
c
20.688 GPa
k
n
/
k
s
1.833
μ
0.590
σ
c (mean)
88.916 MPa
σ
c (std. dev.)
36.950 MPa
τ
c (mean)
83.768 MPa
τ
c (std. dev.)
28.932 MPaThe microparameters can successfully be determined by the improved simulated annealing algorithm. This saves a tremendous amount of hard work, as it is not necessary to change the microparameters and manipulate PFC2D by hand for each step of calibration, which is convenient for use.


## 3.1.1. Basic Information about the Numerical Simulations

To illustrate the calibration process of microparameters more specifically, an example was selected. In [78], sandstone was utilized to investigate the Kaiser effect. Moreover, the PFC2D was employed to reproduce the Kaiser effect, the deformation, and the rate of deformation. Meanwhile, the contact-bond model in PFC2D was utilized, and the size of the numerical simulation for the macroparameters was 61.14mm×153.32mm. According to the physical experiments, the UCS and Young’s modulus of the sandstone were 83.1 MPa and 14.3 GPa, respectively.In total, 10 microparameters determine the contact-bonded model in PFC2D. The microparameters and their corresponding descriptions are listed in Table1.Table 1
Microparameters of a contact-bonded material.
Parameter
Description
R
m
i
n
Minimum ball radius [m]
R
m
a
x
/
R
m
i
n
Ball size ratio, uniform distribution []
ρ
Ball density [kg/m3]
E
c
Ball-ball contact modulus [Pa]
k
n
/
k
s
Ball stiffness ratio []
μ
Ball friction coefficient
σ
c (mean)
Contact-bond normal strength, mean [Pa]
σ
c (std. dev.)
Contact-bond normal strength, std. dev. [Pa]
τ
c (mean)
Contact-bond shear strength, mean [Pa]
τ
c (std. dev.)
Contact-bond shear strength, std. dev. [Pa]As listed in Table1, the density of the ball can be treated as the density of the sandstone, which is 2397 kg/m3 [78]. The effect of gravity does not need to be considered, as the specimens are small in the physical experiments and the gravity-induced stress has a negligible effect on the macroparameters [77]. In addition, the microparameters should be in a certain range, which guarantees the success of the numerical simulations. For example, Rmax/Rmin, the ball size ratio, should be greater than 1, or the numerical simulation cannot run successfully. Therefore, it is necessary to determine a reasonable range of microparameters. Table 2 gives the range of the microparameters of the contact-bonded model in the numerical simulations.Table 2
The range of the microparameters of the contact-bonded material.
Parameter
Description
Range
Rmin
Minimum ball radius [m]
0.1e − 3~1e-3
Rmax/Rmin
Ball size ratio, uniform distribution []
1~20
Ec
Ball-ball contact modulus [Pa]
0~100e9
kn/ks
Ball stiffness ratio []
0~10
μ
Ball friction coefficient
0~1
σ
c (mean)
Contact-bond normal strength, mean [Pa]
0~1000e6
σ
c (std. dev.)
Contact-bond normal strength, std. dev. [Pa]
0~100e6
τ
c (mean)
Contact-bond shear strength, mean [Pa]
0~1000e6
τ
c (std. dev.)
Contact-bond shear strength, std. dev. [Pa]
0~100e6It should be noted that the minimum ball radius should be larger than0.1e-3 m, or the error “Error: Fewer balls generated” will occur during the numerical simulations, which will stop the numerical simulation from running. The upper bound of the microparameters in Table 2 can be increased slightly, but this will increase the computation times.


## 3.1.2. Basic Parameters of the Improved Simulated Annealing Algorithm

Some basic parameters of the improved simulated annealing algorithm need to be determined before the numerical simulations; these are listed in Table3.Table 3
Basic parameters of the improved simulated annealing algorithm.
Parameter
Value
The length of the Markov chain
100
Cooling efficientk
0.95
The initial temperatureT
100
The number of microparametersN
9


## 3.1.3. Numerical Simulation Platform

In the numerical simulation, the operation system is Linux Mint 18 Sarah (on a personal computer with an Intel Core i5, 3.30 GHz CPU, and 8 GB RAM); the operation system is open source. The Python version is 3.5. Additionally, to implement the improved simulated annealing algorithm and manipulate the PFC2D(3D) software, two extra Python packages were installed,numpy (to accomplish the improved simulated annealing algorithm) andsubprocess (to manipulate the PFC2D(3D) software). The PFC2D version 3.1 in this paper was adapted to run on the Linux system.Wine, an open source software for running the Windows system application, was used for this purpose.


## 3.1.4. Numerical Simulations Results

Based on the improved simulated annealing algorithm, after 2739 iterations, the termination criterion is satisfied and the numerical simulation is terminated. The UCS and Young’s modulus determined by the numerical simulation are 74.356 MPa and 14.824 GPa, respectively, which is quite close to the physical experiment results (UCS 82.1 MPa and Young’s modulus 14.3 GPa.). The corresponding microparameters are obtained and are listed in Table4.Table 4
Microparameters obtained by the improved simulated annealing algorithm.
Parameter
Value
R
m
i
n
0.746e-3 m
R
m
a
x
/
R
m
i
n
1.567
E
c
20.688 GPa
k
n
/
k
s
1.833
μ
0.590
σ
c (mean)
88.916 MPa
σ
c (std. dev.)
36.950 MPa
τ
c (mean)
83.768 MPa
τ
c (std. dev.)
28.932 MPaThe microparameters can successfully be determined by the improved simulated annealing algorithm. This saves a tremendous amount of hard work, as it is not necessary to change the microparameters and manipulate PFC2D by hand for each step of calibration, which is convenient for use.


## 3.2. More Examples

From the analysis of the detailed example in Section3.1, it is found that the calibration method is independent of the numerical simulation model; that is, once the command flow for numerical simulations is fixed, the calibration process can change the microparameters in the command flow  .txt file for each step of calibration. Thus, it can be used for contact bond and parallel bond in both PFC2D and PFC3D. Recall that the numbers of microparameters and macroparameters can be increased or decreased according to the specific circumstances, in which case the corresponding Python scripts should also be changed correspondingly. To further validate the proposed method, more numerical simulation tests were performed. The numerical simulation results are listed in Table 5.Table 5
Macroparameters obtained by the improved simulated annealing algorithm.
Author
Material species
Numerical model
Macroparameters determined by the numerical simulations (error)
Chang et al.2001 [39]
Kimachi sandstone
PFC2DContact-bond
Experimental results: σc: 48 MPa υ: 0.27E: 6.5 GPa
Numerical results: σc: 44.372 MPa (7.558%)υ: 0.275 (1.851%) E: 6.477 GPa (0.353%)
Chang 2002[40]
Hwangdeung granite
PFC2DContact-bond
Experimental results: σc: 162 MPa υ: 0.28E: 50.7 GPaNumerical results: σc: 171.357 MPa (5.775%)υ: 0.283 (1.071%) E: 48.492 GPa (4.355%)
Hong 2004[41]
Daejeon granite
PFC2DContact-bond
Experimental results: σc: 161.7 MPa υ: 0.22E: 50.9 GPaNumerical results: σc: 159.172 MPa (1.563%)υ: 0.239 (8.636%) E: 46.943 GPa (7.774%)
Wang et al.2016 [12]
Rock cored from the rock mass in the Heishan open pit mine in Chengde city, China
PFC2DContact-bond
Experimental results: σc: 132.3 MPa E: 12.0 GPaNumerical results: σc: 130.329 MPa (1.489%)E: 11.924 GPa (0.633%)
Yang et al.2014 [42]
Red sandstone
PFC2DParallel-bond
Experimental results: σc: 52.98 MPa E: 14.42 GPaε1c: 5.113Numerical results: σc: 57.832 MPa (9.158%)E: 15.490 GPa (7.420%) ε1c: 4.650 (9.055%)
Zhou et al.2016 [43]
Concrete
PFC2DParallel-bond
Experimental results: σc: 23.3 MPa υ: 0.17E: 37 GPaNumerical results: σc: 23.525 MPa (0.965%)υ: 0.159 (6.470%) E: 35.822 GPa (3.183%)
Yang et al.2015 [44]
Sandstone
PFC2DParallel-bond
Experimental results: σc: 190.80 MPa E: 35.64 GPaNumerical results: σc: 173.583 MPa (9.023%)E: 39.035 GPa (9.525%)
Fan et al.2015 [25]
Rock-like materials
PFC3DContact-bond
Experimental results: σc: 8.27 MPa E: 4.04 GPaNumerical results: σc: 8.947 MPa (8.186%)E: 4.228 GPa (4.653%)
Duan et al2015 [45]
Beishan granite
PFC3DParallel-bond
Experimental results: σc: 136.4 MPa υ: 0.266E: 43.167 GPaNumerical results: σc: 127.017 MPa (6.879%)υ: 0.265 (0.375%) E: 42.652 GPa (1.193%)
σ
c: Uniaxial compressive strength; υ: Poisson’s ratio E: Young’s modulus; σt: tensile strength; ε1c: peak axial strain.As listed in Table5, calibrating the microparameters via the improved simulated annealing algorithm is a viable strategy for both parallel bond and contact bond in PFC2D(3D). It has wide practical use in PFC2D(3D) numerical simulations.


## 4. Discussion

PFC2D(3D) is a typical DEM, which is widely applied to model the damage and nonlinear behaviors of rock materials. However, the microparameters in PFC2D(3D) need to be determined before the numerical simulations. In theory, the microparameters can be calibrated on the basis of macroparameters determined by the physical experiments. In practice, determining the proper set of microparameters based on the macroparameters obtained from the physical experiments is very difficult but is critically important for the success of the PFC2D(3D) numerical simulations.The microparameters are mainly calibrated based on the macroparameters determined by the physical experiments. In recent years, the most commonly used method for calibrating microparameters in PFC2D(3D) is called the “trial and error” method [24]. Its main process for calibrating the microparameters is comparing the macroparameters determined by the numerical simulation and the physical experiments; then, the microparameters are changed according to the difference between the macroparameters of the numerical simulation and the physical experiments. However, the calibration process is quite subjective; it is influenced by the experiment and experimenter significantly. In principle, this means that if a user is familiar with the influence of microparameters on the macroparameters, it will not take very long. In contrast, if a user is inexperienced, the calibration process will be lengthy. Another significant issue is that this calibration process is accomplished by hand. In the calibration process, each step of calibrating microparameters requires a series of repeated work: changing the microparameters in the command flow  .txt file, opening the PFC2D(3D) software, and calling the command flow  .txt file. With increasing iterations, it takes much time and work. Therefore, accomplishing the calibration process automatically has great practical utility in PFC2D(3D) numerical simulations. As described above, there exist two unsolved problems in calibrating microparameters: how to avoid the subjectivity in the calibration process and how to accomplish the calibration process automatically.In this paper, the improved simulated annealing algorithm was applied to calibrate the microparameters based on the macroparameters determined from the physical experiments, and the calibration process was completely accomplished by Python scripts (the corresponding Python scripts in this paper will be published on GitHub (https://github.com/)). The improved simulated annealing algorithm was mainly utilized to obtain the proper set of microparameters. By using this method, the microparameters are changed on the basis of the improved simulated annealing algorithm; when the macroparameters of the numerical simulations are quite close to the macroparameters obtained from the physical experiments, the calibration process terminated. The calibration process was accomplished by Python scripts, eliminating the necessity of human input at each step. Additionally, this method can be applied to different kinds of the bond models in PFC2D and PFC3D. Consequently, the method proposed in this paper has a practical use in PFC2D(3D) numerical simulations.However, the calibration process is realized in the Linux operation system in this paper; meanwhile, the PFC2D(3D) are mainly used in the Windows operation system. Thereafter, whether the complying Python scripts are suitable in the Windows operation system needs to be further investigated and will be our next work.


## 5. Conclusions

To determine a proper set of microparameters based on the macroparameters obtained by physical experiments, the improved simulated annealing algorithm was utilized to calibrate the microparameters of the PFC2D(3D) models. Moreover, the Python scripts were developed to accomplish the calibration process successfully. The main conclusions of this paper are as follows.(
1
) The improved simulated algorithm was utilized to calibrate the microparameters based on the macroparameters obtained from the physical experiments, avoiding the subjectivity in the calibration process. Moreover, the Python scripts were developed to calibrate the microparameters by changing the microparameters in the command  .txt file for each step of calibration based on the improved simulated annealing algorithm. This method could be implemented automatically, and the calibration program terminated when a proper set of microparameters were obtained.(
2
) By using the improved simulated algorithm, a set of numerical simulations were performed. It is indicated that the proposed approach can be used for the calibration of microparameters. This method can be used for the contact-bond or parallel-bond models in PFC2D(3D). The method has a practical value for the PFC2D(3D) numerical simulations. However, the corresponding Python scripts were developed on the Linux operation system platform, and it is necessary to realize the calibration process on the Windows operation system; it requires further investigation.



---

*Source: 6401835-2017-04-26.xml*