Numerical simulation of rock cutting using YADE

Laboratory rock cutting tests are commonly used for the study of the cutting process and the evaluation of the expected performance of excavation machines. The cutting process is dominated by a great number of parameters, most of them relating with the intrinsic mechanical characteristics of the rock. In this study the open software Yade is used for simulating laboratory cutting tests on sandstone samples using a drag pick as the cutting tool. The process is simulated based on a conventional formulation, combined with an enhancement of the sample microstructure through the manipulation of the interaction range coefficient, which provides for the sample a very realistic initial strength ratio UCS/UTS. On cylindrical samples, four different cuts in four different paths were carried out. The mean cutting force for each cut was calculated and was chosen to represent the macroscopic response of the numerical model. The optimum set of microparameters is obtained through an experimental design with the Placket-Burman and Central Composite Design methods, and then optimized, in regard to the microparameters’ values, so that the rock cutting simulation is in close accordance with the observations from the actual laboratory cutting tests.


Introduction
In rock excavation, and mining in particular, machines are commonly used for directly attacking the rock. Drag tools, which are attached on these machines, are responsible for the actual excavation, that is for the creation of rock chips. When examining the action of a drag tool, the force is typically resolved into three forces, namely the cutting force, the thrust force and the sideway force. Although these three forces are acting simultaneously during cutting, the force responsible for the removal of rock is the cutting force. Estimation of the cutting force magnitude is possible through either laboratory cutting tests, theoretical and empirical models, or numerical modeling. In the absence of appropriate equipment for the conduction of rock cutting tests, and due to the limited scope of theoretical and empirical criteria, the use of numerical modeling is recommended.
In this particular field, the discrete element method (DEM) is realistically reproducing the mechanical behavior of rock by representing the system as an assembly of bonded spherical particles. Due to its discontinuous nature and the simple constitutive laws that govern DEM, the method provides a tool to study the initiation and propagation of cracks in brittle materials. However, DEM fails in achieving high UCS/UTS ratios for the assembly, which are characteristic of brittle rocks regardless of the micromechanical parameters [1]. Clustered models, which form large complex-shaped particles, were used as a first solution to study the influence of the cluster size on the UCS/UTS ratio and on the strength envelope [1,2]. Although the results from the clustered models were very encouraging, cracking was observed only along the cluster boundaries, due to the very high strength assigned on intracluster bonds. Rock cutting simulations were conducted in 2D and the relation between the critical depth of cut, which is responsible for the transition from ductile to brittle failure, the material toughness KIC, and the unconfined strength σ C was found [3]. During rock cutting simulations the grain crushing phenomenon was taken into account in order to achieve a realistic strength ratio [4]. The cluster logic [2] was adopted in the sensitivity analysis in combination with the inter cluster strength. As a result, the cluster can break, thereby simulating the case where big rocks can be crushed into fine particles during rock cutting, which is commonly observed in experiments.
In order to realistically simulate the mechanical behavior of rocks by using the discrete element method, a calibration procedure is needed for identifying suitable values for the microparameters. This is frequently accomplished by 'trial-error' [1]. Statistical methods are also used, so as to reduce calibration time and provide acceptable results [5]. Moreover, the experimental design and the optimization method of enumeration was used in rock cutting simulations with PFC3D in [6].
In this paper, the open-software Yade [7], which has been used before for solving rock engineering and geomechanics problems [8], provides the computational framework for simulating the rock cutting process. In order to achieve realistic results, firstly, the interaction range coefficient [9], as defined for spherical elements, was used to enhance the numerical texture by yielding a macroscopic strength ratio UCS/UTS close to that recorded in actual laboratory experiments. Then, an experimental design of several stages was applied on the rock cutting simulations: A sensitivity analysis of the effect the microparameters have on the mean cutting force was conducted. Then, the Placket-Burman method was used for identifying the two microparameters with the highest impact on the mean cutting force and a linear relation between them was derived. In addition, by applying the Central Composite Design method, a quadratic relation between the microparameters and the mean cutting force was also obtained. The Sequential Least Squares Programming (SLSQP) optimization method was implemented for finding the best set of microparameters, so that the difference between the rock cutting force recorded in laboratory tests and the mean cutting force from the simulation was minimized. Finally, the performance of the finetuned rock cutting model was evaluated by executing a series of simulations and comparing the results to the measurements from actual laboratory tests.

Laboratory rock cutting tests
At NTUA's Laboratory of Excavation Engineering rock cutting tests were executed on sandstone samples [10], with mechanical characteristics as listed in Table 1. The geometrical specifications of the test setup and the cutting tool are illustrated in Figure 1. As shown, each cylindrical sample was cut in four different paths, parallel to its longitudinal axis, and the cutting force was continuously recorded. In Figure 2 representative cutting force measurements can be seen. The mean cutting force for each cut was calculated and is listed in Table 2.

Numerical simulation of rock cutting tests
In simulations of the bonded particle model with Yade [7], a rock material can be described as a collection of spherical elements clasped together at their contact points and confined by walls. The evaluation of the particles' kinematic behavior (i.e. acceleration, velocity, displacement, etc.) is based on the fundamental principle of dynamics [11]. A basic contact model supported in Yade is the contact bond model, which approximates the physical behavior of a vanishingly small cement-like substance lying between and joining the two bonded particles. Specific values in terms of tensile and shear strength are defined for the contacts on the bonds. Because of the nature of the contact bond, it is not able to resist to bending. The microstructure of the assembly is enhanced by the use of the interaction range coefficient, which increases the density of bonds between elements [9].
For the creation of a bonded particle model in Yade, the following microparameters have to be defined: the Elasticity Modulus (E) between spheres, the stiffness ratio (KNKS), the friction coefficient between spheres (phi), the normal (t) and shear bond strength (c), and the upper and lower boundaries of particle radii of the spherical elements. A numerical specimen with more than 10000 discrete elements, provides the best compromise in regards to both computational efficiency and realism [9].
For the sandstone described in Section 2, a Uniaxial Tension Strength (UTS) of 8.6 MPa has been estimated from the Brazilian Tensile Test (BTS) by using the correlation provided in [12]. Thus, the macroscopic strength ratio σ σ ⁄ is 15.8. By changing the values of interaction range coefficient that correspond to different bond densities, as seen in Table 3, a series of uniaxial tension test simulations were carried out on the same numerical assembly, in order to calibrate the microparameters. The macroscopic tensile behavior was matched for the values of the interaction range coefficient , interparticle stiffness E , and local tensile strength t that are shown in Figure 3. The stress-strain curves for uniaxial compression strength test simulations corresponding to microparameters in Table 3 are shown in Figure  4. The simulated macroscopic strength ratio obtained from UTS and UCS is listed in the last column of Table 4. It is concluded that an interaction range coefficient equal to 1.25 provides the more realistic macroscopic strength ratio.   Table 3). To constrain the movement of the numerical assembly due to the action of the cutting tool, only the top wall of the specimen was deleted, thus allowing the recording of the cutting force and chip formation. The cutting tool was modelled as a rigid body by using facet elements to resemble the actual cutting tool's geometrical characteristics and was moved at a fixed speed in the direction of the z-axis. In order to achieve computational efficiency and a good dynamic profile, a cutting speed of 2.0 m/s was chosen. The cutting force was continuously recorded during the simulation process. A local damping coefficient, a, of 0.8 was specified, which provided good results in uniaxial tension and compression tests.

Plackett-Burman
The microparameters having the highest significance were determined with a Plackett-Burman design. By varying each microparameter independently, its effect on the mean cutting force was examined separately [14]. The microparameters that have been chosen for the sensitivity analysis were the stiffness ratio (KNKS), the interparticle elasticity modulus (E), the strength ratio (c/t), and the interparticle friction coefficient (phi) [15]. The variation of the mean cutting force with the change of each microparameter is shown in Figure 6. The factors that result in positive correlation to the mean cutting force are the interparticle friction phi and the bond strength ratio c/t, while the influence of E and KNKS is not obvious. Based on the sensitivity analysis, the plus and minus level for each microparameter is presented in Table 6, where the equations that relate the coded and uncoded values are listed in the last column. The mean cutting force was chosen as the macroscopic response in the calibration process. The complete Placket-Burman design matrix is shown in Table 7 with the simulation results in the last column.   The results of the Plackett-Burman screening design are presented in Table 8. The microparameter having the highest impact on each variable is underscored. The impact of each microparameter is determined by the non-negative value of the effect, meaning that higher values correspond to greater significance of the effect on the response. The negative or positive correlation of the effect is determined by its sign. The significance of the microparameter on the response is determined from the p-value. An acceptable level of significance is usually 0.05. If a microparameter is having a p-value less than or equal to   Table 9 provides the order of each microparameter from the strongest to the weakest on the macroscopic response. The mean cutting force is mostly affected by the tensile strength t and the interparticle friction coefficient phi. With respect to the microparameter t, the higher its value the higher the force required to break the bonds.

Central composite design
From the Plackett-Burman method the microparameters having the highest impact on the mean cutting force were identified. In order to calculate the non-linear relation of these two microparameters with the mean cutting force, the Response Surface Method was adopted.
The Central Composite Design used in this study was based on an imbedded factorial with center points that have been increased by a group of star points for the estimation of curvature. In the design space, when the distance from center to a factorial point is 1 unit for each factor, then the distance from the center to a star point will be a. In this study, two elements with two levels exist, therefore the number of elementary runs is four, resulting √2 [16]. The complete design matrix and the result of the CCD analysis along with the coefficients are presented in Tables 10 and 11 respectively. The interaction between the two factors, along with the non-linear relation to the response, is given in Table 11, from which the quadratic equation (2) was constructed. The contour line for the response is plotted in Figure 7.

Implementation and optimization validation
The Sequential Least Squares Programming (SLSQP) method was used to minimize the absolute difference of the simulated and measured value of the mean cutting force: Objective function: Minimize|MCF Laboratory test result| 3 The microparameters were also required to meet equality constraints of the linear fit: Inequality constraints: |MCF Laboratory test result| 1 4 1 x 1 5 Last step was to check the adequacy of the optimization result. The coded values for a set of four cuts, shown in Table 12, were uncoded by using the transformation equations from the last column of Table 6. These were used to rebuild the numerical model. In Figure 8, the cutting force diagrams obtained from the numerical simulations are presented.