Abstract
This study aims to identify discrete element model parameters of rock-like materials. An inverse procedure is developed to determine the discrete element model parameters from experimental measurements. This involves the solution of an inverse problem through minimizing the misfit function which describes the error between numerical computation and experiment by an optimization procedure. In this procedure, the discrete element method is adopted as the numerical calculation method of the forward problem. The orthogonal experimental design is used for parameter sensitivity analysis. Besides, the approximation model with radial basis function is adopted instead of the actual calculation model to reduce the time of forward calculation. The ant-colony optimization algorithm is employed as the inverse operator. Therefore, the parameters of the discrete element model are optimized by this procedure. The three-point bending experiment with discrete element simulation is provided to verify the validity and accuracy of the inversion results. The results indicate that it can rapidly obtain the available and reliable model parameters just through a few sets of experimental data. As a result, this inverse procedure can be applied more widely to parameter identification of the discrete element model for brittle materials.
1. Introduction
With the rapid development of rock engineering, especially the tunnel construction for high-speed rails [1], shale gas exploitation [2], and construction of protective engineering [3], the working efficiency and life of large rock breaking machinery such as the tunnel boring machine (TBM) are widely concerned. In the process of rock breaking of the TBM, hard rock is very prone to brittle failure, which is easy to cause disasters such as low pressure and rock burst. Therefore, it is most significant and essential to understand the mechanical performance and breakage behaviors of rock materials when they are suffering from different loading. Traditionally, mechanical properties of rock materials can be studied through experimental measurements [4–7]. However, the laboratory experiments or tests are the costly and time-consuming process. Thus, numerical computation has provided a new way for investigating the mechanical properties of materials. Among these numerical methods, the discrete element method (DEM) is one of the best tools to study the mechanical property of rock materials [8].
The DEM considers many kinds of discontinuities, material failure characteristics, and fracture modes. It is commonly used to numerically compute the mechanical properties of rock materials by considering them as a collection system with a set of units such as springs, beams, or separate particles bonded together with contact. Many researchers have studied the mechanical properties of rock-like materials through the DEM [9–14]. These research results showed that the DEM had a great advantage in micromechanical representation, which provided a good avenue to understand material failure characteristics.
For numerical calculation using the DEM, an accurate set of material input parameter values is needed, which determines the accuracy of numerical simulation results. Some microscopic input parameters can be determined directly by experiments, while some parameters describing the mechanical behavior of materials at the microscopic level can only be obtained indirectly. Generally, an initial value is used in the DEM and gradually adjusted until the error between the experiment and the simulation reaches the allowable value [15]. This way of parameter identification is called the “trial and error” method. Yang et al. studied the relationship between the microscopic parameters and the macroscopic parameters of the parallel bond model and found that these quantitative relationships between microscopic and macroscopic properties are empirical [16]. Tan et al. determined the DEM parameters with a set of tests, comprising the uniaxial compression test, Brazilian test, and fracture toughness test [17]. However, the “trial and error” method for identifying the microscopic parameters depends on previous experience, with randomness and blindness. Simultaneously, this method is very expensive and time-consuming. In order to overcome these obstacles, it is necessary to seek an effective way to determine these microscopic input parameters. By now, the inverse method for parameter identification is more suitable.
The inverse method for parameter identification is aimed at searching the optimal input parameters such that the model responses best match the experimental data. In recent years, the computational inverse method has made good progress in the field of model parameter identification [18–23]. The inverse method for identifying the DEM parameters uses the complex relationship between the macroscopic responses from experimental measurements and the microscopic input constants. This complex relationship is generally reconstructed by a given mathematical model that describes forward problems. Hence, as a number of accurate experimental responses are obtained, the unknown microscopic input parameters of the DEM can be determined through solving the inverse problem, which is properly established.
The aim of this study, therefore, is to propose an effective inverse technique for determining the model parameters of the DEM for rock-like materials. The remainder of this study is organized as follows: In Section 2, the outline of the inverse procedure is presented and uniaxial compression and three-point bending tests of granite were carried out to provide input data and validation data. Moreover, the microscopic parameters were successfully determined by the inverse analysis, which then combines with forward calculation, sensitivity analysis, approximate model, and optimization algorithm. In Section 3, the identified results are given and discussions on the three-point bending tests are performed to verify the reliability of the inverse procedure.
2. Inverse Analysis: An Objective Identification Process
The inverse process for the DEM parameter identification for rock-like brittle materials is shown in Figure 1. In the process, firstly, it is necessary to analyze the mechanical behaviors of the material and perform the physical mechanics experiments and then clarify the corresponding forward and inverse problem forms. According to the mechanics experiment, the corresponding forward problem model is built. Secondly, in order to ensure the existence and solvability of the inverse problem, the model parameters need to be analyzed by sensitivity analysis and then the inversed parameters are determined. Then, the forward solver needs to be repeatedly called in the parameters’ inverse process. Therefore, the approximate model is usually adopted instead of the numerical calculation model to enhance the computational efficiency. Furthermore, the calculation response obtained by solving the forward problem is compared with the experimental measurement response, and the inverse objective function is established according to the practical problems and the requirements of the solution. Furthermore, it is necessary to use the efficient and reliable optimization algorithm to obtain the inverse results. In this work, the objective function is defined to estimate the gap between experimental measures and numerical data. Objective function used here iswhere is the vector of inversed parameters, is the response measured from the experiment, is the computational response, and is the number of data points.
Finally, if the convergence criteria can be satisfied, the inversed parameters are obtained. On the contrary, if not, it is needed to add new samples and calculate the forward problem again to reconstruct the approximate model. Furthermore, the previous steps are repeated until the criteria of convergence are satisfied, as shown in Figure 1.
2.1. Mechanical Property Tests for Granite
2.1.1. Sample Preparation and Experimental Equipment
Combining the inverse technique, the following two basic physical mechanics experiments are carried out to obtain a prepared and reliable source of experimental data. They are the uniaxial compression test and three-point bending test.
As shown in Figure 2, the rock sample selected for the test is granite. The rock samples used in the test are prepared according to the relevant experimental rules of rock mechanics. All granite samples were taken from the same large rock mass with good homogeneity. During the coring process, the drilling directions are parallel to each other to ensure the verticality and parallelism of the experimental sample. The end face of each sample was carefully ground using a stone grinder. Among them, the samples used for the uniaxial compression tests are cylindrical samples with a diameter of 50 mm and a height of 100 mm. The samples used for the three-point bending tests are cylindrical specimens having a diameter of 50 mm, a height of 220 mm, a slit length of 18 mm, and a slit width of 2 mm. The sample size meets the requirements of the International Rock Mechanics Institute for determining the type I static fracture toughness method. In the tests, AX and CX were used to number the test samples in three different test modes, where A represents a uniaxial compression test, C represents a three-point bending test, and X is the serial number of the sample under each test mode. Among them, the uniaxial compression test uses the INSTRON 1346 electrohydraulic servo test machine, and the three-point bending test uses the INSTRON 1342 electrohydraulic servo test machine. The installation position of the test device and samples are shown in Figures 3 and 4.
(a)
(b)
(a)
(b)
2.1.2. Experimental Results
The force-displacement curves and stress-strain curves of the uniaxial compression test are shown in Figure 5. The uniaxial compression strength (UCS) of the granite sample has a maximum value of 140.25 MPa, a minimum value of 137.19 MPa, and an average value of 138.8 MPa. Young’s modulus is up to 44.29 GPa, the lowest Young’s modulus is 39.95 GPa, and the average Young’s modulus is 41.5 GPa. Poisson’s ratio is usually expressed as the ratio of the average of the lateral strain along the X-axis to the average of the lateral strain along the Y-axis. So based on the experimental measurement, the value of Poisson’s ratio is at most 0.24, the minimum is 0.22, and the average is 0.23. Through the analysis of the above data, it has been found that the UCS, Young’s modulus, and Poisson’s ratio between the samples are not much different, the error of which is in the range of 10%. These indicate that the homogeneity of the granite sample is better. The influence of material heterogeneity on the results of uniaxial compression tests is within the error range. Meanwhile, the peak strain corresponding to the peak stress of sample A1 was found to be 0.425%, the peak strain of A2 was 0.451%, and the peak strain of A3 was 0.482%. The axial strain of the sample before failure was between 0.425% and 0.482%, both less than 1%. According to the principle of brittleness and toughness of brittle materials, the axial strain before material failure is less than 3% that of the brittle material. Therefore, the granite samples have strong brittleness.
(a)
(b)
The force and load point displacement curves of the three-point bending test are shown in Figure 6. The maximum force of the rock sample is 1.509 kN, the minimum force is 1.396 kN, and the average force is 1.453 kN, and the fracture toughness can be obtained by the following formula:where is the straight incision depth, is the diameter of the test piece, is the maximum load of the fracture failure, and is the distance between the two support points. In the experiment of this paper, is 160 mm.
By calculation based on the above formulas, the maximum fracture toughness is 1.062 MPa·m1/2, the minimum is 1.048 MPa·m1/2, and the average is 1.057 MPa·m1/2. When the maximum load occurs, the displacement of the load point is basically concentrated at 0.15-0.16 mm.
The above test data were collated, and the obtained mechanical properties of granite are summarized in Table 1.
2.2. Forward Calculation
In this work, the DEM was used as the numerical calculator of the forward problem. The DEM should be corresponded with the actual material test. In this study, the uniaxial compression test was used to identify the DEM parameters. Based on the tests, the DEM could be built. As the scale of the specimen and the applied load were plane stress problems, a two-dimensional axisymmetric DEM was built, as shown in Figure 7. In this numerical model, for the uniaxial compression test, the length W, height L, and density ρ of this model were 50 mm, 100 mm, and 2622 kg/m3, respectively.
In the DEM, the number of particles is very important for the numerical model. Usually, the more the number of particles is, the higher the accuracy of the DEM is. Some studies have shown that as the number of particles is more than 6000, the model could be more accurate [24]. In this work, the minimum radius of particles was 0.25 mm, the average particle diameter was 0.33 mm, the maximum-to-minimum radius ratio was 1.65, and the number of particles was 12,952.
A suitable numerical model for describing the mechanical behaviors of the material is of great importance for the computational accuracy of the DEM. There are many numerical models for the material in the DEM, such as the contact stiffness model, contact sliding model, and bond particle model (BPM). The BPM is one of the most commonly used models, which is a combination of particles, in which the particles are joined together by parallel bond. This model can represent many features of material behaviors, such as elasticity, fracturing, and failure. Hence, the BPM was used to describe the mechanical behaviors of granite in this work.
2.3. BPM and Its Parameter Analysis
Potyondy and Cundall proposed the BPM to reproduce many features of rock behavior [15]. The BPM contains two forms, which are contact bond and parallel bond. As shown in Figure 8, the parallel bond is not only able to transmit the tension and shear force between the particles but also able to transmit the moment and torque between the particles.
The contact force between the particles would be divided into normal force and shear force aswhere the contact force vector denotes the contact behavior of particle A on B, as shown in Figure 9, is the normal force component, is the shear force component, is the normal unit vector of the contact plane, and is the tangential one.
The magnitude of the tangential force is calculated in increments. When the two particles are in contact, is initially zero, and then it is incremented with the corresponding tangential displacement increment. The normal force and the increment of shear force are calculated bywhere is the normal stiffness value at the contact, is the shear stiffness value, denotes the overlap of the couple of particles at contact, denotes the increment of , and denotes the relative shear displacement.
Once the tension or shear force exceeds the limit of normal strength or tangential strength between the particles, the bond between the two particles should break. At the same time, the corresponding force and torque can be removed. Figure 10 shows the action of the bond at contact [25, 26].
Hence, the DEM constructed with the BPM contains eight microscopic parameters, which are listed as follows: the friction coefficient of the particles , the radius multiplier of the parallel bond which is used to set the parallel-bond radii, the ball effective modulus , the bond effective modulus , the normal-to-shear stiffness ratio of the ball , the normal-to-shear stiffness ratio of the parallel bond , the tensile strength of the parallel bond , and the shear strength of the parallel bond . Some previous studies indicate that, for brittle materials, the radius multiplier of the parallel bond is usually set to be 1. The ball effective modulus and the bond effective modulus are set to the same value, and the stiffness ratio of the ball and the stiffness ratio of the parallel bond are set to the same value [27]. However, it is difficult to determine the other five parameters (, , , , and ) by the existing theory or test method. Therefore, this paper proposes the inverse procedure to identify these microscopic input parameters.
2.4. Sensitivity Analysis
Generally, the input variables should be highly sensitive to the output for reducing the ill-posed problem during parameter identification. In other words, the input and output should have a great causal-effect relationship. As the orthogonal experimental design (OED) method is an available tool to investigate the problem with many factors and levels, in this work, it was adopted for sensitivity analysis of unknown parameters. Then, the friction coefficient of the particles , the effective modulus of the particles , the normal-to-shear stiffness ratio of the particles , the tensile strength of the parallel bond , and the shear strength of the parallel bond were selected as variables and an orthogonal array was adopted, as listed in Table 2. Through forward numerical computations, the three macroscopic response variables (Young’s modulus , Poisson’s ratio , and uniaxial compressive strength ) could be obtained with different sets of unknown parameters. So the DEM calculations of 16 groups of different model microscopic parameters were carried out through the orthogonal experimental table, and the results are listed in Table 3.
The range analysis is one of the commonly used methods in OED. It can show the difference in the computational results was different with different levels. Usually, the greater the difference is, the greater the influence of the factor is. According to the range analysis, the orthogonal experimental data of DEM computation were analyzed. In this study, it is noted that is the summation of the “i” level (i = 1, 2, 3, and 4) in the A factor, , and . The greater the value of R is, the greater the impact of this factor on the result is [28]. Hence, the analysis results are listed in Tables 4–6.
It can be seen from the value of the range analysis that the effective modulus of the particles has the greatest effect on Young’s modulus and its corresponding value is 29.07. The stiffness ratio of the particles has the greatest effect on Poisson’s ratio and its corresponding value is 0.195, while the tensile strength of the parallel bond has the greatest influence on UCS and its is 74.74. In addition, in the range analysis of Poisson’s ratio, the friction coefficient has the second largest impact on Poisson’s ratio. However, the value of stiffness ratio is still more than 6 times the value of the friction coefficient. Therefore, the remaining two parameters friction coefficient and shear strength have little effect on the calculation results. Then, it can be found that the three most sensitive parameters were , , and . Simultaneously, based on the results of parameter sensitivity analysis and experimental measures, the value range of these three parameters could be initially designed. The value ranges of inversed parameters were as follows: (20 GPa and 30 GPa), [1, 3], and (30 MPa and 50 MPa). Moreover, according to other research [29, 30], the friction coefficient of the particles and the shear strength of the parallel bond can be set to 0.5 and 90 MPa, respectively.
2.5. RBF Approximate Model
The approximate model technique was used to replace the actual calculation model for promoting the calculation efficiency. The construction of the approximate model can be divided into three parts, namely, experimental design, model selection, and model evaluation. The Latin hypercube design (LHD) method was adopted for multifactorial experimental design. It has a great sample-recording ability and can effectively avoid repeated sampling. More importantly, it has good equality and can extract samples distributed in the boundary position. In this work, 30 samples were produced by the LHD method, as shown in Table 7. In this paper, the radical basis function (RBF) was adopted, which has simple principles and high accuracy for high-dimensional nonlinear problems, and the Gaussian function was used for RBF. The approximate model can be described aswhere is the weight coefficient of the linear combination of the function , is the Euclidean distance between the unknown point and the sample point, and represents the radial basis function, as shown in Table 8.
Thirty samples were used to construct the RBF approximate model, and other three random samples were used to verify whether the model was accurate. The actual sample output and the expected output of test samples are shown in Table 9, and it can be found that the basic prediction error was within 5%. Therefore, it indicates that this RBF approximation model is reliable.
Hence, the objective function can be expressed aswhere is the vector of inversed parameters (, , and ) and is the vector of responses (E, υ, and σ).
2.6. Ant-Colony Optimization Algorithm
To promote the computational inverse accuracy, the ant-colony optimization (ACO) method was adopted for the optimization algorithm. It is a very effective metaheuristic algorithm to solve the complicated engineering problems. In this algorithm, the search way is according to population, and it has very strong parallel computing capability. The ACO algorithm has several advantages, such as distributed computation, constructive greedy heuristics, and positive feedback [31]. The greatest advantage of this method is that it can tune a lot of parameters simultaneously and it can be used to solve optimization problems.
Because of these advantages, the ACO algorithm was employed as the inverse operator to determine the inversed parameters of the DEM. When using this algorithm, the number of ant colonies m is 50, the maximum number of iterations is 200, the information heuristic factor α is set to 1, the information residual factor ρ is 0.5, and the expected heuristic factor β is 5. The objective function information is applied to guiding and judging the search process. The relative errors between the measured data and the numerical computation results based on the inversed parameters are below 5%. If these above conditions are satisfied, these inversed parameters can be considered reliable.
3. Identification Results and Discussion
3.1. Results
Based on the above-mentioned inverse procedure, parameters of the DEM for granite were determined and are shown in Table 10. The convergence curves of ACO for determination of objective function are shown in Figure 11(a). It can be seen that the ACO has a fine convergence performance in this study. The computational efficiency curves of ACO for three inversed parameters are shown in Figures 11(b)–11(d). They show that the computational efficiency is fine, and the largest number of generations is about 80, which occurs in the process of inversely determining the stiffness ratio.
(a)
(b)
(c)
(d)
In addition, the value of objective function was 0.0133, which satisfies the evaluation criteria. The parameters were substituted into DEM numerical calculation, and the results and the experimental curves are shown in Figure 12.
For estimating the goodness of fit between the calculated results and the experimental data, the relative error RE and the correlation coefficient CC were used. The functions of the expression are stated as follows:where denotes the number of sampling points and denotes the average of the responses.
According to the evaluation principle, the smaller the RE and the larger the CC, the higher the coincidence degree of the two curves. If the RE is 0 and the CC is 1, the two curves are completely coincident. Through the calculation, the two curves’ RE is 0.035 and the correlation coefficient CC is 0.996. According to the evaluation principle, it can be shown that the numerical results based on the inversed parameter values are in good agreement with the measurement curves from the experiment.
At the same time, the effective and reliable DEM parameters determined by the inverse procedure are substituted into the DEM positive problem calculation. The simulated sample damage form is compared with the experimental sample damage form, as shown in Figure 13. Because the loading surface is wide and the damage form is greatly affected by the material homogeneity, the comparison results will be somewhat fuzzy but overall consistent.
(a)
(b)
3.2. Discussion
The three-point bending test is used to verify the model parameters. In the corresponding DE modeling process, the three-point bending test is simulated as a plane stress problem, which is treated as a two-dimensional DEM. The basic dimensions are consistent with those of the test samples, and the detailed dimensions are shown in Figure 14. In this numerical model, for the three-point bending test, the length W, height L, and density ρ were 220 mm, 50 mm, and 2622 kg/m3, respectively. In this model, the minimum radius of particles was 0.25 mm, the average particle diameter was 0.33 mm, the maximum-to-minimum radius ratio was 1.65, and the number of particles was 28458.
The DEM parameters identified in the previous work are applied to the three-point bending model, and the loading conditions are set to be consistent with those of the test. The comparison between the obtained macroscopic mechanical indexes and the test results is listed in Table 11. It can be found that the error is about 10%, within the accuracy requirements. As shown in Figure 15, it can be clearly seen that the simulated sample damage results almost coincide with the damage results from the experiment, and it indicates the availability of the inversed parameters. The above work also indicates that this inverse procedure can be effectively applied to the parameter identification problem of DEMs.
(a)
(b)
4. Conclusions
In this work, the microscopic parameters of the DEM for granite have been determined by using the inverse procedure. To identify the microscopic parameters, the error between the synthetic measurements and the calculated data has been minimized by the inverse technique. A few of uniaxial compressive experimental data are used as the input for determining the unknown parameters. Then, the parameters of the DEM are identified through the above work, and the agreement between the computational curves and the experimental results is good. Subsequently, it verifies the accuracy of the inverse results with the three-point bending tests. And the results demonstrate the availability of this inverse technique. As a consequence, the inverse procedure for parameter identification can be applied to quickly determine the unknown model parameters with high efficiency and accuracy, which provides an effective understanding of the mechanical properties of rock-like brittle materials.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (11602212, 51605409, 11802258, and 51775468) and the Natural Science Foundation of Hunan Province of China (2018JJ3509).