Discrete element modeling of the machining processes of brittle materials: recent development and future prospective

In recent years, the discrete element method (DEM) has been used to model bulk material, especially brittle materials (such as rocks, ceramics, concrete, ice, etc.) with various mechanical properties or responses by setting serials of contact properties (such as bonds) in the particle assembly. These bonds can withstand a certain amount of force and/or moment, so that the stresses executed in the bond can be used for determining the initiation and propagation of micro-cracks. There are increasing evidences over the last 20 years that the DEM is becoming an effective numerical method to simulate the cracking, crushing, and deformation of continuous media under external loads. The DEM has now been widely used in the field of processing and machining of rock, ceramic, concrete, and other brittle materials. In this paper, the theoretical principles, formulations, and contact models as well as the numerical solving processes of DEM are introduced. The applications of DEM for the machining processes of brittle and rigid materials such as ceramics are described and reviewed in detail, and the future development trend is also discussed.


Introduction
With the rapid development of modern technologies in industry, hard-brittle materials (such as engineering ceramics, semiconductors, single crystal silicon, sapphire, rock, etc.) are more and more widely used in various fields. For example, engineering ceramics (such as silicon nitride, silicon carbide, zirconia, etc.) due to their excellent properties have been widely used in the aerospace, optical, semiconductor, mechanical, and military fields [1,2]. Semiconductor materials such as silicon carbide as power devices will have important applications and broad development prospects in 5G communications, electric vehicles, intelligent Internet of Things, and other fields. Due to the inherent hardness and brittleness of the hardbrittle materials, micro-cracks and surface and sub-surface damages are usually induced during the machining operations [3,4]. This process is not only related to the structure, shape, and characteristics of the tool but also closely related to the microstructure and mechanical properties of the workpiece, and the machining parameters are also one of the main influencing factors. How to improve the processing efficiency of brittle materials, clarify the machining damage mechanism, and explore the factors affecting the strength of parts has become an emerging challenge.
In fact, many phenomena that are difficult to observe in the physical experiments can be investigated using validated computational simulations. All of these simulation techniques are very useful to investigate the damage mechanisms in machining of brittle materials. Among them, according to how the numerical domain is discretized, they can be roughly categorized into continuum-based methods (CBMs) (such as FEM, XFEM, and BEM) and discontinuum-based methods (DBMs) (such as SPH, MD, and DEM). The most important simulation method among the CBMs is the FEM, which is mainly based on the continuum mechanics to define the stress-strain relationship. There is no doubt that FEM is the most popular and powerful simulation method for engineering and physical problems with multiphysic fields and complex load boundary conditions. However, due to the fact that associated governing equations are based on the continuum mechanics, FEM is difficult to describe the transformation process of continuous body to discontinuous body during the processing of brittle materials [25]. Therefore, it is difficult for FEM to take into consideration the influence of the material microstructures, such as defects, dislocations, etc., and these microscopic properties are likely to be closely related to the macroscopic failure of the brittle materials. In the fracture simulation process of FEM, the form of embedded crack source is usually needed and the area of interest needs to be re-meshed, which is timeconsuming and tedious. To overcome the shortcomings of FEM, Moës et al. [8] developed XFEM to model the fracture problems. Compared with FEM, XFEM can simulate the crack initiation and crack propagation process along an arbitrary or the related path without re-meshing. BEM as another method of CBMs is also used to simulate the crack propagation problems in the brittle materials [22,23]. BEM is advantageous in reducing the problem dimension and representing the arbitrary curve and boundaries much easier. However, it is usually an issue to obtain the approximate Green functions for the solution of the boundaries.
In parallel with the development of CBMs, the mesh-free methods have also attracted extensive attention from research scholars, and gradually developed and formed the DBMs. SPH was proposed by Gingold [10] and Lucy [26] in 1977, and it has been well developed in the field of fluid mechanics [9]. In the field of brittle material processing, SPH is often used to study the crack formation mechanism and material removal process during the machining process [11,24,27]. As one of the DBMs, MD used Newton's motion theorem as its basic principle, which can well study the motion, interaction force, stress, strain, heat, temperature, and temperature field distribution of particle motion in the system, and changes in the internal microstructure of the solid material and the resulting lattice defects. At present, MD has become a powerful tool for analyzing micro-machining mechanisms [28]. However, MD simulation technology has not yet formed a complete theoretical system, which is not perfect in the establishment of potential functions of many brittle materials, and has certain limitations in the simulation scale.
The discrete element method (DEM) was originated from MD and it is a numerical method specifically designed to solve the problem of continuous or discontinuous media. The obvious advantages are applicable to the simulation of the deformation and failure process of discrete particle assemblies under quasi-static or dynamic conditions, as well as the analysis of large dispersion, large deformation, and significant flow. When the DEM is used to solve the continuum problem, it regards the material as a collection of discrete particles, which are connected by bonds, and the damage of the bond is used to describe the generation of micro-cracks in the material, so it can effectively study the damage and crack in the machining process of brittle materials. The hard-brittle materials are usually treated as a collection of bonded particles, so the generation and the propagation of micro-cracks in the materials are described as the damage of bonds. In this way, the damage and crack that happened in the hard-brittle material during the machining process can be effectively reflected [21,29]. Due to the advantages of DEM different from the above numerical methods, it is widely used in the fields of rock [17,[30][31][32], soil mechanics [33][34][35], brittle material processing [19,29,36], etc. In the simulation process, the external temperature field and the internal heat generated by friction, grain deformation, or breakage of materials can be considered and controlled, which makes the DEM become a new method to simulate the high-temperature mechanical properties [37] and high-temperature damage and destruction process of brittle materials [38]. Meanwhile, since the hard-brittle materials normally contain a lot of micro-cracks with different sizes, shapes, and orientations, these internal defects will affect the mechanical properties [39][40][41] and processing quality [27] of the hard-brittle materials (such as ceramics). By deleting particles to generate random defects inside the DEM model, researchers can simulate the physical structure of the real material and explore the crack propagation and mechanical properties of the brittle material containing the defect under external load [21,42]. Since the DEM was first applied to the simulation of rock materials, researchers in various countries have gradually used DEM as a useful simulation method for the motion of discontinuous materials and the damage of continuous materials. Figure 2 shows the scientific papers published around the keywords DEM or discrete element method in the past 24 years, which the statistical data were obtained from the Web of Science. It can be seen that the published paper based on the DEM increase exponentially. It also indirectly proves that the DEM becomes an effective numerical simulation method for engineering applications.
Aiming to enrich and develop hard-brittle material machining simulation methods, this paper provides a complete review around the basic principles, modeling methods, and application examples based on the discrete element method. It is divided into 6 sections. Following the introduction, Section 2 gives an introduction on the theory of DEM from the aspects of particle motion, contact searching, particle packing, constitutive equations, and determination of timestep. Section 3 explores the relations between microscopic parameters and macroscopic properties, and discusses the common methods of model calibration and validation. Section 4 lists a few of application examples of DEMbased modeling of the machining process of brittle materials, such as cutting, grinding, milling, polishing, and other machining methods. Section 5 discusses the advantages and disadvantages of DEM with other numerical methods, and the challenges and opportunities of DEM are also discussed. And Section 6 presents a short summary of this review paper.

Basic equation of particle motion
The DEM is firstly introduced by Cundall [43,44] in the 1970s for the analysis of rock mechanism problems and then applied to the study of soil materials [45]. It has been widely used to analyze the complex mechanical model of discontinuous media, such as granular materials. However, DEM has several advantages over the CBMs. For example, DEM can trace the detail motion of each individual particle unit and can also simulate the generation, propagation, and coalescence of cracks as the bond breakage between the particles. At present, the types of discrete units are mainly divided into two categories: non-spherical units and spherical-based units. Nonspherical units include polyhedron, dilated polyhedron, super-quadric particles and a lot of irregular elements, etc., which are often used for modeling rock materials. Sphericalbased units include circular particles, spherical particles, ellipsoid particles, and combined spherical particles. The combined spherical particles can be divided into cluster and clump, which can be used to simulate super-particles or particles with complex structural shapes in the simulation domain. If you want to consider the crushing of super-particles, you can use the cluster model; otherwise, use the clump model. Spherical particles are one of the earlier discrete units. Due to the convenient contact search and simple equations of motion, spherical particles are widely used in the modeling of hardbrittle materials. In the following description, unless otherwise specified, all particles refer to spherical particles.
In the DEM, each divided unit has its own independent motion; the overall motion law can be obtained by studying the motion law of each unit, which means the overall motion law can be observed by analyzing and calculating each independent motion unit. The formula describing the movement of the basic unit is as follows: where m i is the quality of particle i; θ ! is the angular vector of particle i; r ! i is the distance from the centriod of the particle i to the contact point, J i is the moment of inertia of particle i; q ! ij is the arm of force from the acting force to the center of the particle i, F ! ij is the force applied to particle i by particle j, N i is the neighbor number of particle i, F ! t e is the force applied to particle i, and K ! t e is the external torque applied to particle i. The motion equation of each element is solved by the explicit finite-difference method for continuum analysis, and then, the overall movement form of the research object can be obtained. In the numerical simulation process, the dynamic behavior of elements is represented by a timestep algorithm in which it is assumed that the velocities and accelerations are constant within each time step. The DEM uses the iterative method such as the dynamic relaxation method or the static relaxation method to perform loop iterative calculation and to determine the force and displacement of all contacts in each time step, and then updates the displacement changes of all contacts at any time.

Contact searching for spherical particles
In the DEM simulation, it is necessary to determine the contact relationship of the units before calculating the contact force between them. Since the DEM is usually used to deal with the mechanical problems of discontinuous media such as granular particles, and the mechanical problems of continuous media to discontinuous media, a large number of particles often come into contact or separation over time. Therefore, the contact searching process is computationally intensive and time-consuming. The process of judging the contact relationship of discrete units is called the contact searching process, and its algorithm is called the contact searching algorithm. Contact searching has become one of the bottlenecks restricting the DEM for large-scale calculations. In order to meet the needs of large-scale computing and expand the application range of DEM, a large number of scholars have conducted extensive and in-depth research on the DEM contact searching algorithm.
The simplest and direct contact searching algorithm is to judge each discrete unit in the system one by one with the other discrete units. If there are N discrete units in the system, the calculation of this algorithm is O (N 2 ). It can be seen that, if the number of system units is large, the algorithm will consume a lot of calculation time. In recent years, scholars have proposed a series of efficient contact searching algorithms suitable for large-scale DEM systems, such as NBS (No Binary Search) algorithm [46], C-Grid (D-Cell) algorithm [47], ADT (Alternating Digital Tree) algorithm [48], ASDT (Augmented Spatial Digital Tree) algorithm [49], etc. These contact searching algorithms are mainly divided into two categories: tree-based contact searching algorithms and latticebased searching algorithms. Among them, the ADT algorithm and ASDT algorithm belong to the tree-based contact searching algorithm. If there is a good tree structure, the average calculation amount is O (NlnN), while the NBS algorithm and C-Grid (D-Cell) algorithm belong to the grid-based searching algorithm, and its average calculation is O (N) in  April 15, 2020) theory. Han, Feng, and Owen [50] made a systematic comparison of several commonly used contact searching algorithms for NBS, C-Grid, ADT, and ASDT, and concluded through numerical experiments. Their research results show that the computational efficiency of the grid-based search algorithm is affected by the grid size, and the tree-based contact searching algorithm is affected by the tree structure. For larger-scale contact judgment, the computational efficiency of the NBS algorithm and the C-Grid algorithm is much higher than that of the ADT algorithm and ASDT algorithm [50]. If readers need to know more about various contact algorithms, please refer to related literatures [51][52][53][54][55].

Particle packing
In order to establish a DEM model similar to the microstructure of hard-brittle materials, it is necessary to closely arrange particles with different sizes to achieve the specified porosity and initial stress state. It usually includes the following 5 steps: (1) Randomly generate a certain number of particles in a specified area, and give a radius change range of the particles, as shown in Fig. 3a.
(2) Calculate the current porosity n and compare it with the target porosity n 0 to calculate the radius enlargement factor m. The relevant calculation equations are as follows: (3) Enlarge the particle radius and let it move under the action of unbalanced force, and finally reach equilibrium with desired isotropic stress, as shown in Fig. 3b. (4) Find particles with contact number less than n c (n c ≥ 3) and mark them as floating particles (such as blue particles in Fig. 3c), and use the method of enlarging the particle radius to eliminate the floating particles in the simulation domain. (5) Set the microscopic parameters between the particles (see Section 3.1), delete the walls, and relax the model to obtain a block model of the hard-brittle material, as shown in Fig. 3d.

Contact model and its failure criterion
In order to make the continuous model of hard-brittle material to bear a load, the adjacent particles are usually connected by "cement." Common particle connection models include bonded particle model (BPM) and cohesive beam model (CBM). When the active load between the particles is greater than the strength set by the BPM or CBM, then the bond between the particles will break, resulting in the initiation of a micro-crack. When the BPM or CBM is broken between the particles, the contact action between the particles is generally described by a linear elastic model or Hertz model etc., and the frictional sliding between particles can be characterized by a sliding model [31].

Bonded particle model
By using the bonds at the contact between particles, the DEM can reproduce the complex mechanical behavior of solid materials. Bonds conceptually represent the materials can be used to carry a lot of load. If stresses are introduced below the given value of the bonds, the simulated materials may act similar to a continuum approximation of a solid. However, according to whether the bonds can withstand a torque, it can be divided into linear contact bond model (LCBM) and liner parallel bond model (LPBM). LCBM provides the behavior of an infinitesimal and linear bonded at the contact and it can only use to carry a force for tensile or shear. That is to say if a LCBM was set between two particles, and one of them was fixed, then the other particle can rotate around the first particle under the combined action of gravity and bonded force. LPBM refers to the characterization of a stacking particlebonded material in a finite size. This bond creates an elastic interaction between the particles, which does not prevent the occurrence of slip, and can transmit both force and torque, but the force can only be transmitted at the point of contact, as shown in Fig. 4. Due to the existence of parallel bond stiffness, the relative motion on the contact (when the parallel bond is created) will cause forces and moments to act on the bond material. Therefore, compared with LCBM, LBPM is more suitable for the modeling of hard-brittle materials. The calculations for the maximum tensile stress σ max and maximum shear stress τ max acting at the parallel bond are: where the F n i and F s i are represented as vectors of normal and tangential parts, respectively; M n and M s are the normal and tangential moment parts associated with the parallel bond; J is the polar moment of inertia of the parallel bond cross section; A is the area of the bond cross section, and I is the moment of inertia about the axis through the contact point and in the angular direction of Δθ; and R L is the radius of the parallel bond.
If the maximum tensile stress exceeds the normal strength of the bond (σ max ≥σ t ) or the maximum shear stress exceeds the shear strength of the bond (τ max ≥ τ t ), then the parallel bond will break [51], as shown in Fig. 5.

Cohesive beam model
The cohesive beam model was first introduced by Stanley [53]; this model was initially used in lattice networks and then has been developed in DEM filed. And this model is able to simulate the behavior of brittle and elastic media, which is defined by Poisson's ratio, Young's modulus, and tensile failure strength. Figure 6 shows the particle model boned by the  cohesive beam model and the two independent parameters' length L beam and radius r beam which can describe the cylindrical geometry of the beam. However, the cohesive beams are massless, and the L beam depends on the distance between the discrete element's center. In addition, there are two mechanical property parameters of cohesive beam: Young's modulus and Poisson's ratio. Similar to the parallel bond model, the micro properties of the cohesive beam model do not directly reflect the macro properties of the material model, and researchers usually need to calibrate its parameters when using the cohesive beam model [54].
The calculations for the maximum normal stress σ n μmax and maximum shear stress τ μmax acting at the cohesive beam are [55]: where the σ nt μmax and σ nb μmax are the maximal normal stresses due to the tensile loading and bending loading, respectively; τ μmax is the maximal shear stress; r beam is the beam radius; I μ is the moment of inertia along y and z; M μx , M μy , and M μz are the torsional moment along x, the bending moment along y, and the bending moment along z in the cohesive beam; Io μ is the polar moment of inertia; N μ is the normal force in the cohesive beam; and the S μ is the section of the cohesive beam.
If the maximum normal stress exceeds the normal strength of the beam (σ μmax ≥ σ μ ) or the maximum shear stress exceeds the shear strength of the beam (τ μmax ≥ τ μ ), then the cohesive beam will break.

Determination of time step
In order to reliably update the model state, one must determine a suitable time step for the integration of Newton's laws. The determination of the time step in the DEM preferably does not exceed the critical time step associated with the minimum natural period of the entire system, because a small time step will make the disturbance of the unit only propagate to the nearest unit and the calculation efficiency extremely slow. On the other hand, a large time step will increase the calculation amount and diverge the calculation result. Therefore, it is very important to choose the appropriate time step. Consider a onedimensional mass-spring simple harmonic motion system described by a point mass m, and spring with stiffness k, the critical time step corresponding to a second-order finite-difference scheme can be calculated as [56]: where i and j designate two particles in contact; m ij is the reduced mass which is defined as m ij = 2m i m j /(m i + m j ), where m i and m j are the masses respectively associated to i and j particles; J ij = min(J i , J j ), where J i and J j are the moments of inertia respectively associated to i and j spherical particles; K n and K t are the classical normal and tangential stiffness and R ij = min(R i , R j ), where R i and R j are the radius respectively associated to i and j particles.

Microscopic parameters and macroscopic properties
When using DEM to simulate the fracture or machining of bulk materials, the primary task is to calibrate the microscopic parameters in the model. This calibration process is usually a typical inverse problem and it is also one of the difficulties in the DEM simulation process. Unlike the continuum simulation method, the DEM usually cannot directly input the mechanical properties of the material through laboratory tests, Fig. 6 Cohesive beam model between the particles because the material parameters are set from the microscopic scale to the contact behavior between particles in the simulation model.
However, if the target model is simulated by regularly arranged particles, it is possible to directly derive the correlation between the microscopic parameters and the macroscopic properties through the analytical solution method. For instance, Thornton [57] provides expressions for the strength of face-centered cubic arrays of spheres in 2D. For a columnar array of disc-shaped particles of thickness t (in which each particle has four neighbors), expressions for modulus and tensile strength are easy to obtain. The elastic modulus E and tensile strength σ t of materials can be calculated as follows: Here, k n is the normal stiffness of the particle, t is the thickness of the disc-shaped particles, S n is the normal strength of the bond between the particles, and R is the radius of the particle.
For a 3D regularly arranged particles model, Wang and Mora [58] have developed the relationships between microscopic and macroscopic properties by the principle of energy conservation. Their research results show that the normal stiffness of the particle is related to the macroscopic Young's modulus E, Poisson's ratio v, and particle size, and Poisson's ratio v of the material will be affected by the ratio of the particle's normal stiffness k n to the tangential stiffness ks . The relevant calculation formulas are as follows: In addition, other scholars have also studied and deduced the relationship between the microscopic parameters and macroscopic properties which focus on the estimation of the rotational stiffness or spring stiffness [59][60][61]. Although there may be some differences in the final relationship, it still provides theoretical support for us in the initial selection of material microscopic parameters.
However, it should be pointed out that the above research is derived on the basis of the equal-diameter and regular arrangement of particles. In order to represent the random defects in those brittle materials, random packing of particles is preferred, and as a result a calibration process is needed. As mentioned in Section 2.4, most researchers used the bonded particle model or cohesive beam model to cement the particles together at the contact position for simulating the block materials. As shown in Table 1, the microscopic parameters of the parallel bonded model and the cohesive beam model are listed.

Calibration methods
In the BPM model of brittle materials, the characterization of these material parameters is achieved by setting the bond model between particles [21]. In the BPM model, the mechanical properties of particles and parallel bonds are defined on a microscopic scale, such as the particle stiffness, the coefficient of friction between the particles and the strength of the parallel bond and so on. These microscopic parameters are difficult to obtain by test methods, and there is no method or theory that directly matches the macroscopic mechanical properties of brittle materials with the parameters in the corresponding BPM model. Therefore, the relationship between the mechanical properties of the actual material and the properties of discrete particles and parallel bonds must be established macroscopically by means of test calibration. In order to determine whether the BPM model of the established brittle material is reasonable, scholars mostly calibrate the BPM model by the uniaxial compression test (UCT), Brazilian test (BT), threepoint bending test (TPBT), single edge notched beam test (SENBT), and uniaxial tensile test (UTT) [19,62], as shown in Fig. 7.
Furthermore, if the CBM model was used to simulate the behavior of brittle materials, some microscopic parameters, such as Young's modulus, Poisson's ratio, and the length and radius for the beam, were needed to be defined [63]. With the same like BPM, how to build the transition laws between the microscopic parameters of the cohesive beam and the macroscopic parameters of the materials is the main difficulty. For example, André et al. [54] used the numerical quasistatic uniaxial tensile test to calibration these microscopic parameters, and investigated the relations between the microscopic parameters and the macroscopic mechanical properties.
For the microscopic parameters of the above two particle contact models, trial-and-error method was usually used to calibrate the material model. However, this method is inefficient, especially for beginners. In addition, there are other calibration methods that have also attracted widespread attention. Detailed discussions are as follows.

Trial-and-error method
The trial-and-error method is an empirical learning method that achieves the goal through continuous experimentation and elimination of errors. Since the attribute microscopic parameters entered in the DEM cannot directly reflect the macroscopic physical properties of the materials, and there is currently no universally applicable parameter correspondence relationship for researchers to use, many researchers choose the trial-and-error method to calibrate their material microscopic parameters for different brittle materials [63,64]. In order to develop more efficient method for DEM calibration, the method of design of experiments (DOE) was used to improve the trial-and-error process in recent years. Yoon et al. [65] used a Plackett-Burman design and response surface analysis to search the suitable microscopic parameters for uniaxial compression of rock materials. In order to optimize the DEM calibration, Hanley et al. [66] applied the Taguchi method to analyze the relation between the input parameters and the bonded agglomerate, and the mechanical response of each agglomerate was measured in a uniaxial compressive test simulation. It can be seen that using the DOE method to calibrate the microscopic parameters of particles should find the relationship between the model inputs (microscopic parameters) and outputs (macroscopic parameters), which requires many computational models to be run. Nevertheless, for a large number of micro-parameters that need to be calibrated, the trial-and-error method is cumbersome and stupid, because that there is no scientific basis for adjustment to give a guide for how to modify the multiple parameters. Although the trialand-error method can be used to obtain the relationship between the microscopic parameters and the macroscopic  attributes in the specific DEM model, the obtained conclusions require a large amount of computing resources.

Optimization-based method
Unlike the trial-and-error method, the most commonly reported optimization-based method is a computational alternative model. The optimization calculation process is used to replace the tedious numerical simulation process, and the main idea is to use iteration to approximate the optimal solution of the problem as much as possible. The procedure includes several steps as follows: (1)  Chen et al. [68] developed an inverse procedure for the calibration of BPM model used in the rock-like materials. The schematic chart for the parameter calibration process is shown in Fig. 8. The uniaxial compressive tests were used to determine the unknown parameters for input, and the threepoint bending test was used to verify the accuracy of the inverse results. The research results show that after about 80 iterations, the calculated results tend to be stable. Feng et al. [67] used the optimization method to carry out an automated calibration procedure for bonded SiC ceramic materials. The UCT and TPBT models were performed with Adam algorithm in a few iterative steps to match the bond parameters. The schematic chart of the calibration process is shown in Fig. 9. The results show that the method proposed by Feng has higher calibration accuracy with less iteration. It can be seen that this method transforms the inverse method for parameter identification to an optimization problem, which it can quickly and effectively determine model parameter by few material experiments and simulations. However, it is not advisable to reverse too much of the parameters; otherwise, the calculation accuracy will decrease, especially when there is a correlation between the parameters, the convergence will be reduced.

Validation tests
When using the DEM to build a model of hard-brittle materials, the macroscopic mechanical properties of the material are closely related to the microscopic parameters in the model, which has been mentioned in Section 3.1. Meanwhile, one of the inherent advantages of using DEM simulation is the analysis of the crack's propagation, so the failure criterion for the bond is very important. Section 2.4 has mentioned the failure judgment methods of BPM and CBM. According to the fracture mechanics of hard-brittle materials, it can be known that the failure behavior of hard-brittle materials under complex load conditions is often caused by tensile failure, shear failure, or mixed failure. In order to simulate the shear failure behavior of materials, researchers usually use a compression test [21,42,69] to study the mechanical properties, and use this to calibrate the shear parameters in the micro-parameters. Concurrently, the tensile test [62], Brazilian test [16,70], or bending test [71,72] is usually used to study the tensile failure of materials and to calibrate the tensile parameters in the micro-parameters. In the calibration and validation process, the numerical simulation method is used to simulate the conventional mechanical properties of the material, and the macromechanical properties obtained by simulation are compared with the macro-mechanical properties obtained by experiment. If the error between them does not exceed the set allowable value (usually set to 5%), then the microscopic parameters of the material are calibrated. For various validation tests, detailed discussions are as follows.

Compression test
Compression test is a common mechanical property test method used to evaluate the shear failure of brittle materials [64,[66][67][68]73]. Its schematic diagram and DEM simulation failure diagram are shown in Fig. 10 a and b. During the compression test, the shape of the specimen is usually a cylinder or a rectangular parallelepiped. Place the specimen in the center of the lower platen of the testing machine, and let the upper platen move towards the lower platen at a slow and uniform speed, and then record the load and displacement curve of the upper platen (Fig. 10c). Record the maximum load during the process of compression failure, and the compressive strength of the material can be calculated according to Eq. (14). At the same time, by recording the deformation of the specimen in the radial and axial process during compression, Young's modulus and Poisson's ratio of the material can be calculated, as shown in Eq. (15) and Eq. (16).
Compressive strength σ mc : Young's modulus E c : where |F c | max is the peak axial force, d T is the diameter of the specimen, σ a is the axial stress, ε a is the axial stress, and ε r is the radial strain.

Tensile test
Unlike compression testing, tensile testing is mainly used to evaluate the elastic-plastic properties of metal materials and the tensile strength of brittle materials [63,74]. As shown in Fig. 11, clamp the specimen by the upper and lower chucks to perform low-speed stretching until the specimen breaks, and record the force-displacement curve during the tensile test of the specimen. According to Eq. (17), the tensile strength can be calculated. However, as the tensile strength of the brittle material is much smaller than the compressive strength, it is easy to cause the specimen to break when clamping. Therefore, tensile testing is rarely used directly to measure the tensile strength of the brittle materials.
Tensile strength σ b : Fig. 9 Schematic chart of the calibration process for bond parameters (reproduced from ref. [67]; Copyright 2020, with permission from Elsevier) where |F b | max is the peak force and d c0 is the diameter at which the specimen breaks.

Brazilian test
The Brazilian test is also called indirect tensile test, which is mainly used to measure the tensile strength of brittle materials. It is a common method for testing the tensile strength of rocklike materials [63,75]. As shown in Fig. 12, a linear load is applied in the diameter direction of the cylindrical specimen to cause it to fail along the diameter direction of the specimen. Record the maximum load during the process of Brazilian failure, and the tensile strength of the material can be calculated according to Eq. (18). Because the Brazilian test process is relatively easier, the measured tensile strength is close to the tensile strength measured by the direct tensile method, so this method is commonly used to determine the tensile strength of hard-brittle materials. Tensile strength σ B : where |F a | max is the peak force and d B and t B are the diameter and thickness of the Brazilian disk, respectively [76].

Three-point bending test
The three-point bending test is a common test method for testing the flexural strength of materials [67,68], and it is also used to calibrate the tensile strength of microscopic parameters. The schematic diagram and DEM simulation failure diagram are shown in Fig. 13. In this test, the tested specimen is placed on two support points at a certain distance, and a downward load is applied above the midpoint of the two support points. Let the indenter above the specimen start to move vertically downward slowly, and the displacement-load curve of the material begins to be recorded (Fig. 13c). Record the maximum load during the failure process of TPBT, and the flexural strength of the material can be calculated according to Eq. (19). Flexural strength σ bb : where L P is the span of the specimen, |F P | max is the peak force, b P is the width of the specimen, and h P is the height of the specimen.

Single edge notched beam test
In order to study the sensitivity of brittle materials to cracks or defects, researchers usually use an artificially prefabricated defect to study their mechanical properties, such as the single edge notched beam test [19,74]. This test method has a notch on one side of the center of the test specimen and pre-made a sharp crack, and applied pressure using the three-point bending method, as shown in Fig. 14. For instance, Tan et al. [19] used the single edge notched beam test when validating the fracture toughness of the polycrystalline SiC model. The fracture toughness of the material can be calculated according to the peak load the specimen during the failure process, as shown in Eqs. (20) and (21). Fracture toughness K Ic : where F i is the peak force, L N is the span of the specimen, h N is the height of the specimen, b N is the width of the specimen, and a is the length of the crack in the specimen. Table 2 lists the basic contents of the above five validation methods. In addition to the above commonly used verification methods, there are other test methods (e.g., torsional test [63], cracked chevron notched Brazilian discs [79]) that can also be used to validate the macroscopic properties of the material model.

Application examples of DEM-based for modeling the machining of hard-brittle materials
The above description shows that DEM has a unique advantage in dealing with the initiation and propagation of microcracks in brittle materials without the requirement on mesh generation, which has attracted many scholars to carry out studies in machining process [21,36,[80][81][82][83], machining damages [19,84,85], and surface quality [86][87][88] of brittle materials under various processing conditions.

Construction of machining model
When using DEM to simulate the machining of hard-brittle materials, how to establish a reasonable tool and workpiece model is the first problem that researchers need to solve. In addition, when paying attention to the impact of chips on tool life, the construction of chip models is also very important. Generally, there are two methods for modeling the tools: wallbased method [18,36,80,83,86,[89][90][91][92][93] (Fig. 15) and  (Fig. 16). If the tool has been set with an ideal stiffness, and the influence of deformation and damage on the workpiece is ignored, you can usually use a folding wall to build the tool (represented as crosssection in two dimensions). If the morphology of the tool involved is more complicated, the geometric model of the tool can usually be generated by the 3D software, and then imported (such as stl format) into the DEM model [97,98]. The tool generated by this method is simplified into a folding wall with ideal stiffness, the deformation will not occur in the simulation process, and the researchers can simulate different processing conditions by changing the location or velocity of the wall. For example, Jiang et al. [19,82] used the wall-based method to establish a cutter, and a series of machining simulations of ceramic materials were carried out. However, if the tool damage and machining quality of the workpiece are the research objects, the tool models can be made up of a lot of particles. For example, Jiang et al. [99] designed a wheel model by particle-based model to study the process of crushing or shedding of abrasive particles on the surface of the grinding wheel, as shown in Fig. 17. However, when using particles to build a tool model, researchers need to calibrate the microscopic parameter of the model reasonable. The details about calibration and verification methods are shown in Section 3.
In the DEM simulation of the machining of hard-brittle materials, the generation process of the workpiece model can be found in Section 2.3. Particle-based method can also be used to build chips to simulate the impact of chip flow behavior during machining for studying the tool damages. Jiang and Peng et al. [74,84] have studied the crack propagation and damage of coated tool during machining and predicted the wear damage of ceramic tool, as shown in Fig. 18. The establishment of the chip model is similar to that of the workpiece model. As the tool enters the steady state of cutting, the contact length between the tool rake face and the chip is  [95]; Copyright 2020, with permission from Elsevier) basically unchanged. In order to reduce the particle scale of the simulation model, the ends of the chip can be set as periodic boundaries [74,84,100]. Similarly, some scholars set the boundary condition of the abrasive part as a periodic boundary [96] when studying the process of polishing, which greatly simplified the calculation model and improved the simulation calculation efficiency. Table 3 lists some research results of hard-brittle material machining simulation based on DEM simulation. It can be seen that the workpiece model is usually constructed by particles and connected by bonds. However, the tool model is usually constructed by walls. At the same time, in order to study the wear and damage of the tool, the tool model can also be constructed by the particles.

Machining simulation
Due to the advantages in dealing with the micro-scale problems, the nonlinear deformation and the failure problems, DEM is widely used in the simulation of the machining process of brittle materials [30,36,92,104]. The processing of brittle materials such as ceramic, rock, and glass are generally processed in precision and ultra-precision, especially in grinding [105][106][107][108] and polishing [109][110][111]. In addition to traditional machining methods such as cutting, milling, and grinding, there are also some emerging auxiliary machining methods, such as abrasive water-jet machining [112,113], laser-assisted machining [13,92,105,107,[114][115][116][117][118][119][120], ultrasonic machining [121][122][123][124][125][126][127], electron beam machining [128][129][130][131][132], chemical machining [111,133], electro-chemical machining [133], and so on. In the current literatures about the machining simulation of hard-brittle materials through DEM and its coupling method, a lot of works focus on studying the material removal mechanism, crack propagation, machining quality, and the optimization of the machining parameters. Some of the works of machining methods performed on hard-brittle materials are shown in Table 4. This review paper presents detailed discussions from the application of DEM and its coupling method in the machining simulation of hard-brittle materials.

Cutting simulation
Cutting is one of the most fundamental ways to achieve material removal. Because of its few influencing factors, in the process of researching grinding, milling, and polishing, researchers usually simplify the effect of cutters or abrasive particles on the workpiece to single abrasive particle cutting or micro cutting. Due to the unique mechanical properties of hard-brittle materials, it is easy to cause cracks and pits on the workpiece surface during the cutting process. At present, numerical simulations of hard-brittle materials mainly focus on the cutting mechanism and the optimized of machining parameters.
Peng and Tan et al. [78,148] used the commercial software PFC to simulate the dynamical cutting process or scratching process of alumina and silicon ceramics. In their simulation, the cutting tool is made up of folding walls with ideal rigidity. Different cutting depth and speed of the tool were carried out to simulate the initiation and propagation of cracks within the materials, and the effects of cutting speeds, cutting depths, and tool rake angles on the number of cracks, the depth of maximum crack, and chip formation after machining were also investigated. Additionally, by analyzing the residual stress  [82], as shown in Fig. 19. Roostai et al. [101] used the flat-joint contact model to construct the quartz ceramics, and carried out a serials of DEM simulation of the dynamic cutting process of the quartz ceramics to study the damages of the surface/subsurface of the workpiece, chip formation, and cutting force by different processing parameters. The research results reveal the process of brittle propagation of micro-cracks in the ceramic during cutting.
In addition, in the research of rock material cutting, many scholars also try to use DEM for simulation. Gong et al. [138] tried to use the discontinuous code UDEC to simulate the cutting process of rock mass in 2005. Two types of crack initiation and propagation of jointed rock mass during the TBM process were found in their research. In order to explore the optimal rock cutting conditions of the disc cutter, Moon et al. [83] conducted a series of rock indentation test simulations, and obtained the optimal rate of spacing to penetration. Su et al. [135] attempted to conduct a three-dimensional DEM simulation of rock cutting experiments. As shown in Fig. 20, a conical tool model was used to cut the rock model consisting of graded particles in the simulation. By comparing the experimental and simulation data, Su found that although there is a significant difference, there is also a significant correlation. In 2013, Huang et al. [89] found that the workpiece has a transition from ductile failure to brittle failure in the DEM  Cemented carbide Zhao et al. [98] simulation of rock cutting when the cutting depth increases. Aiming at the transition of the ductile-brittle failure mode of the rock, Liu et al. [136] used PFC2D to simulate the rock cutting process. And they reproduced the transition from ductile to brittle failure during the rock cutting process, and the mechanical specific energy (MSE) is used as the index to establish the critical cutting depth calculation model. After then, a new method by using MSE as the index to identify the transition from ductile to brittle failure mode of the rock is proposed. For studying the cutting performance of rock, Li et al. [80,90] firstly studied the effect of rock brittleness on rock crushing and cutting performance, by changing the ratio of the bond shear strength and tensile strength in the DEM model of rock. Through a series of simulations (Fig. 21), Li found that as the cutting depth increases, the failure mode of brittle rock will change from the toughness/ductile mode to the brittleness as the macro cracks grow. Meanwhile, they found that the formation of a large chip model is similar to the conclusion reached by Huang et al. By using a scale factor to reduce the effective modulus of bonding, tensile strength, and shear strength, the rock model with various degrees of damage can be obtained, and a linear relationship between the scale factor and the damage factor can be also established.
In addition, the DEM simulations of rock cutting with different damage levels have shown that as the damage factor increases, the failure mode of rock cutting will change from the ductile mode to the brittle model. As a new technology, torsional impact cutting (TIC) shows the great rock crushing efficiency, and rate of penetration (ROP). However, few studies have paid attention to the advantages of TIC tools for ROP improvement. Zhu et al. [16] used PFC2D to conduct the simulation of TIC and discuss the rock breakage under steady state and torsional impact cutting. By discussing the relationship of the rock crushing energy, the microscopic crushing process, and the chip characteristic at different impact frequencies and amplitudes, they also found the improvement mechanism of ROP under TIC technology. Qiu et al. [36] conducted 2D and 3D simulations of glass cutting by DEM, and found that the tool rake angle is a significant factor affecting the cutting deformation and cutting force. Alkotami et al. [140] built the DEM simulations of orthogonal cutting of cracked soda lime glass. Meanwhile, they built an algorithm to simulate the surface roughness of the material. In their simulations, different types of seed cracks were prefabricated on the glass model, and a series of orthogonal cutting simulations were carried out. The cutting simulations have shown that the prefabricated cracks in the workpieces can reduce the cutting force, and the cutting force is minimum when the crack angle is 45°or 135°. Moreover, the surface roughness can be minimized when the crack angle is 135°.
In addition, some researchers have focused their research on cutting tools and studied the crack propagation and wear on the surface of cutting tool during the cutting process. In 2015, Jiang et al. [84] modeled carbide-coated tools based on the DEM, and established a tool-chip contact model based on the Merchant cutting model [149]. By applying boundary conditions to the chip to simulate the actual cutting process, the dynamic process of crack propagation of the coating part during the simulation process was recorded, and the influence of cutting parameters on the damage of the coated tool was predicted. Subsequently, in order to study the crack propagation and wear of ceramic tools, Peng et al. [74] established a toolchip model, and used the boundary conditions to simplify the simulation model, as shown in Fig. 22. A thermal-mechanical coupling model was introduced to their model for studying the effect of mechanical heat on tool life. These two examples demonstrate the reliability and superiority of the DEM in predicting tool wear, and more relative research is expected in the future. Table 5 shows some DEM simulation of cutting made by researchers, which mainly based on the simulation software PFC2D/3D.

Grinding simulation
Due to the unique mechanical properties of brittle materials, grinding is one of the most effective methods for precision and high-precision machining of brittle materials [150]. At present, the DEM research on the grinding process of hard-brittle materials is relatively scattered, which can be roughly divided into the following parts: grinding mechanism of hard-brittle materials [18,94,102,144,145], modeling analysis of grinding wheel surface [143,151,152], wear and breakage of grinding wheel abrasives [95,153], etc.
The traditional grinding simulation of hard-brittle materials mostly focuses on the grinding mechanism of the workpiece and the grinding damage of the workpiece. For example, Jiang et al. [18] proposed a DEM model to simulate the grinding process of polycrystalline SiC, as shown in Fig. 23. The relationship between material removal, crack initiation and propagation on the surface of the material, and the changes of grinding force during grinding was investigated. In order to reduce the subsurface damage caused by diamond grinding on the fused silica surface, Blaineau et al. [94] studied the relationship between the subsurface damage of fused silica and the grinding parameters and the mechanical force applied by the grinding wheel. They also carried out a series of numerical researches on the grinding interface, which used the simulation platform GranOO to build a DEM model of the grinding interface. In their model, the abrasives of the grinding wheel were composed of spherical particles bonded by cohesive beams, and the quartz part is considered to be a quasi-rigid body of a rectangular parallelepiped. The results of simulation reflected that reducing the dispersion of the abrasive grains radius of the grinding wheel with large grain size and high abrasive particle concentration can reduce the depth of subsurface defects of fused silica. In order to get a deeper understanding of the removal mechanism of sawn granite, Xu et al. [144] proposed an improved overlapping rigid cluster (ORC) technique to construct irregular grain properties, which ensures the consistency of simulated abrasive grains and the actual abrasive grain shape. The dynamic DEM simulation of granite grinding was carried out in their studies, and the correctness and effectiveness of the numerical calculation method of the surface discrete element were verified. In order to analyze the influence of different tooth tip angles and grinding parameters on the residual stress distribution of the workpiece, Liang et al. [145] dynamically simulated the grinding process of granite by using the single diamond grain grinding model, and proved that DEM is an effective method to analyze residual stress.
Some of the above studies are basically based on the linear scratch test of abrasive particles on the surface of the grinding wheel to study the surface quality and crack propagation of the workpiece, which ignores the impact of the pendulum scratching of the abrasive on the quality of the workpiece. In 2018, Tan et al. [102] explored the impact of abrasive particles randomly distributed on the grinding wheel periodically pounding the ceramic workpiece during the grinding process. Their results showed that the cracks produced by the pendulum scratching test on the workpiece are more obvious than those generated by the linear scraping test. Therefore, in exploring the grinding mechanism, it is necessary to pay due attention to the pendulum scratching of the grinding wheel abrasive grains on the workpiece surface.
From the existing literatures, we can find that the grinding wheel is mostly treated as a single abrasive grain model or a two-dimensional model with a simple shape in DEM simulation. However, the actual distribution of abrasive grains on the grinding wheel surface and the shape of the abrasive grains are often more complex. Constructing a more realistic grinding wheel model is more conducive for researchers to study the grinding mechanism. Osa et al. [151] discussed the potential of using DEM to establish grinding wheel models in his paper published in 2018, and separately introduced the grinding wheel model established by Li [143] and Osa [152]. In order to fill the gap in the problem of numerical simulation research   Fig. 24. Meanwhile, these beams are also served to connect the abrasive grains to the grinding wheel. Similarly, in a numerical model proposed by Osa, the contact length between the grinding wheel and the workpiece during surface grinding has been mentioned; the model of the grinding wheel shares the truncated octahedrons used by Li, and simplify it to a spherical shape. However, a single beam is used to bond each truncated octahedron to build the grinding wheel instead of the network of beams, as is shown in Fig. 25. The purpose is to model the stiffness of the wheel, so they ignored the adhesive fraction and its shape. In the aforementioned study of Blaineau et al. [94], the abrasive grains in the grinding wheel model were replaced by spherical particles and connected by cohesive beams, which is similar to the grinding wheel modeling method proposed by Osa. In addition to the several grinding wheel modeling methods proposed above, some researchers have also established DEM models of abrasive grain [95], as shown in Fig. 16b. But the main purpose of such modeling is to explore the wear of abrasive grain during grinding.
In addition to the two research priorities above, the problems related to the wear and breakage of grinding wheel abrasive grain in the grinding process began to be considered by researchers, because the life of the grinding wheel is largely related to the wear of the abrasive grains. Godino et al. [95] established a single abrasive model and revealed the entire wear tendency and evolution of abrasive during grinding contact based on a workbench-GranOO developed by André et al. [154]. However, this study focuses on the effects of tribochemical reactions and the formation of a third body on the wear evolution of sol-gel alumina and white fused alumina abrasive grains under specific contact conditions, which ignored the possible effects of dynamic mechanical thermal behavior, and this factor deserves our attention. As is shown in Table 6, some detailed setting parameters of various numerical simulations of the grinding process around brittle materials are listed.

Milling simulation
It is shown that micro-milling can also meet the requirements for precision machining of hard-brittle materials and machining complex surfaces [155]. Therefore, the research idea based on "milling instead of grinding" has been adopted by many researchers. In 2010, Shen et al. [92] established a DEM model of laser-assisted milling by converting curved cutting edges in milling to straight cutting edges in two-dimensional cutting, and conducted a laser-assisted milling simulation of silicon nitride ceramics at a specified temperature, as is shown in Fig. 26. By observing the cutting forces and subsurface Fig. 23 Simulation of the grinding process of SiC ceramic (reproduced from ref. [18]; Copyright 2015, with permission from Elsevier) damage at different depths of cutting, they have successfully constructed a simulation model that can predict subsurface damage of materials under laser-assisted milling. And the collaborative experiment proves that the basic mechanism of ceramic material removal in LAM is brittle fracture. In 2013, Wu et al. [93] conducted the low-speed milling dynamic process simulation of zirconia ceramics. The authors have analyzed the formation of cracks on the workpiece surface by different processing parameters, and put forward the feasibility of high-speed machining on zirconia ceramics. In the paper published by Qiu et al. [36] in 2015, the authors also established a DEM model of glass and simulated the indentation, scratching, and micro-milling process to dynamically observe the change of cutting force under different processing parameters and the details of crack initiation and propagation in the glass. Subsequently, Du et al. [104] used PFC3D to conduct a 3D simulation of the milling process of ceramics and simulated the crack propagation in ceramics. And they obtained the changing rules between different milling parameters and milling force, and the number of ceramic material cracks. By comparing with the experimental data, they have verified the effectiveness of the DEM simulation. Nowadays, the DEM has also been applied to the study of asphalt mixture milling. Wu et al. [97] have established a DEM model of asphalt mixture to simulate its milling process and analyzed its mechanical response. In this model, the aggregate in the asphalt mixture is expressed by particles of different diameters, and the viscoelasticity of the asphalt mortar is defined by the parallel bonds among the particles. The milling tool is an imported 3D model, which realizes the authenticity of the milling tool model. Meanwhile, the influences of cutting speed, cutting angle, and cutting depth on cutting force were also studied. However, this research is still in the early stages of exploration. As is shown in Table 7, we list some detailed setting parameters in numerical simulations of the milling process around the brittle materials.

Polishing simulation
The existing polishing techniques for brittle materials are mainly divided into mechanical polishing (MP) and chemical mechanical polishing (CMP) [156]. In the DEM simulation study on the polishing of hard-brittle materials, Iordanoff et al. [157] firstly used DEM to systematically simulate the wear mechanism of the workpiece. Although a simple model was used in this study and the simulation results cannot truly reflect the complex mechanism of polishing, the result can reflect the superiority of the DEM model in the study of abrasive grains in the polishing process. At present, the numerical simulation based on DEM and its coupling simulation with other numerical methods has been able to establish a more complete polishing simulation model after continuous exploration and optimization by scholars [81,147,158]. In these DEM-based simulation researches on the polishing of hard-brittle materials, it can be roughly divided into single abrasive polishing and multiple abrasive polishing simulations.  (a) Single abrasive polishing Single abrasive polishing simulation simplifies the polishing process by simulating the cutting effect of a single abrasive grain on the workpiece, and the researchers can obtain the effects of different processing parameters of an abrasive grain on the material removal and polishing quality of workpiece. Han et al. [86] conducted a single abrasive grain mechanical polishing simulation of ceramic workpieces by PFC2D, as is shown in Fig. 27. In this study, they have studied the formation mechanism of engineering ceramic surface during mechanical polishing, and confirmed the reasonable trend of the correlation between the micro-damage of the workpiece surface and the polishing conditions. Meanwhile, they evaluated the ductile transition of hard-brittle materials under polishing. In order to simulate the consolidating abrasive grain polishing process of optical hardbrittle materials, Wang et al. [103,159] have established an average cutting depth model for consolidating abrasive grain grinding based on the principle of contact mechanics. And they used the angle polishing method to measure the depth of the subsurface damage layer of optical hard-brittle materials with different particle sizes. The results show that the simulation predictions are basically consistent with the experimental results.
(b) Multiple abrasives polishing However, in order to more intuitively simulate the polishing process of hard-brittle materials, many researchers have focused on the polishing simulation of multiple abrasive polishing [81,96,147]. For example, some of them focus on the trajectory of the abrasive flow during the polishing process, or the interactions between the particle and particle or the particle and workpiece [160]. In the work of Xiu et al. [161], the polishing process was simplified to the movement of particles in the parallel plate shear friction system; the complex characteristics of the internal force chain in the abrasive system were observed. They put forward the conclusion to optimize the polishing efficiency. Meanwhile, there are also some studies focusing on the subsurface damage of hard-brittle materials during the surface polishing. As shown in Fig. 28, in order to study the mechanism of subsurface damage in the polishing process of silica more deeply, Iordanoff et al. [96] established a three-dimensional DEM model of silica polishing, and proposed a simplified model for studying the polishing process. In 2020, Zhao et al. [98] established the simulation of the tool edge preparation process based on DEM and hertz contact theory when studying the edge machining mechanism of cemented carbide milling cutter. And they have studied the effects of speed, preparation time, abrasive particle ratio, abrasive particle speed, and direction of rotation on the edge of the milling cutter, which provided a basis for optimizing the trimming effect and efficiency.
In addition to the content described above, in the current multi-abrasive polishing simulation, many researchers couple DEM with other numerical methods, especially in the study of simulating the CMP. CMP is a combination of mechanical polishing and chemical etching, which involves the combination of solid phase and fluid phase; simple DEM simulation cannot fully reflect the actual situation. Tan et al. [81] used the CCFD module in the PFC3D software to conduct a fluid calculation based on the coupling of computational fluid mechanics and computational bulk mechanics, and simulated the flow behavior of solid-liquid two phases in the composite abrasive polishing fluid. By comparison with the results of experiments, the feasibility of using PFC3D software to simulate the two nano-phase flow problems was verified. In the study of Ji et al. [147], to overcome the shortcomings of traditional modeling methods for soft abrasive flow processing of hard-brittle materials, an abrasive flow modeling method based on the coupling of CFD and DEM was used to analyze the distribution of abrasive-wall collision and the material removal of workpiece surface. And they have also studied the uniformity of surface constrained soft abrasive particle flow processing. The coupling of DEM and other numerical methods shows us the broader application prospect of DEM. Table 8 shows some researches of MP and CMP based on DEM simulation.

Other machining methods
In recent years, a series of micro processing methods for brittle materials have appeared, such as laser processing [115], ultrasonic assisted processing [127], electrical discharge machining (EDM) [131], etc. At present, the advantages of non-contact, non-polluting, low-noise, high-precision, and easy-toimplement digital control of laser micro-machining technology make it occupy an important position in the field of fine processing [115]. In modern industry, EDM has become the greatest choice for machining conductive ceramic materials [132]. In the EDM process, there is no contact between the tool and the workpiece, so the surface layer of the workpiece can be quickly melted and removed at a high temperature. EDM uses the electric spark to corrode and produce the brittle materials that efficiently produce the desired shape and size as well as finer surface features [162]. Although laser processing and EDM can greatly improve the efficiency of machining, it depends on the conductivity of the material, which will produce thermophysical effects, and will destroy the surface structure of the workpiece after processing, and even burns and other phenomena. Ultrasound-assisted machining (UAM) is widely used in the processing of brittle materials such as engineering ceramics [163,164], but the traditional ultrasonic processing method is complicated in operation and high in cost.

DEM-FEM
The FEM can be used to simulate structural response [165], elastic wave propagation [166], heat conduction [167], and large deformation and large displacement of slopes on a macro scale [168], but it is impossible to analyze the discontinuous surface or to simulate the progressive destruction process of materials [169]. The DEM-FEM coupling analysis method can simulate the machining process [170,171] and thermal damage [172,173] of intact brittle materials and brittle materials with discontinuous surfaces, and simulate the generation and propagation of cracks after using the relevant criteria of fracture mechanics [170]. Yu et al. [174] has introduced a comprehensive model of single crack and dispersion crack as the stressstrain relationship and failure criterion of the joint element based on the DEM-FEM coupling analysis method, and simulated the initiation and propagation of the crack inside the rock by the initiation, expansion and failure of the joint unit, and then the simulation of the initiation and expansion was validated by the Brazilian split test. Xu et al. [175] established a multi-scale method combining DEM with FEM. The two-dimensional FEM is used to combine the plane stress-strain element with the three-dimensional DEM to deeply study the impact response and damage of complex structures and heterogeneous materials. The process of the mechanism can use the DEM to simulate the local part of interest, and the FEM for macroscopic simulation. And a special transition layer was usually used to connect the discrete element area and the finite element area, which can be observed the pre-stressed damage response under the laser spoke.

DEM-CFD
At present, the description models for two-phase gas-solid flow mainly include the Euler-Lagrange (EL) model and the Euler-Euler (EE) two-fluid model [176]. The EL model [176][177][178] is mainly applied to the dilute phase flow, and the EE model is used for the dense phase flow [176]. However, due to the discontinuity of the particle field itself, the simulation results obtained by the above model methods have some deviation from the actual situation. Discrete element method-computational fluid dynamics (DEM-CFD) is a new method to simulate two-phase flow in recent years, and firstly proposed by Tsuji [179]. The method mainly uses the DEM model to establish a parametric model of the solid particle system, describes the characteristics of collision and agglomeration between particles, and combines the advantages of CFD in the treatment of the gas phase flow field, which can effectively improve the calculation efficiency and precision of numerical solution [180].
Al-Arkawazi et al. [181] coupled the DEM and CFD to construct a simple model for calculating the interaction between fluid and particles, and the effect of porosity on the hydrodynamic behavior of fluids was studied. In this coupling process, the motion of the fluid is mainly described by the Navier-Stokes equation. The motion of the particle is mainly obtained by Newton's second law. The coupling between the fluid-solid phases is realized by Newton's third law.
Based on the above theory, scholars have applied DEM-CFD coupling technology to the field of mechanical processing, especially in the field of abrasive flow processing technology. Li et al. [182] used the combination of CFD and DEM to numerically analyze the processing of abrasive flow during the study of micropore structure precision machining, revealing the micro-cutting behavior of abrasive flow. Ji et al. [147] obtained the abrasive-wall collision distribution and the material removal rate distribution on the surface of the workpiece by DEM-CFD coupling method. Based on this, the processing characteristics and the processing uniformity of the surfaceconstrained soft abrasive flow were studied. The DEM-CFD coupling simulation study can be used in the analysis of the fluid-solid two-phase processing method such as the abrasive flow and the water jet of the brittle material, which can fully exploit the advantages of the DEM.

DEM-SPH
In dealing with the problem of fluid-solid two-phase flow, in addition to the DEM-CFD coupling technology previously mentioned, many scholars have been working on the coupling technique of the smooth particle hydrodynamics method (SPH) and the DEM in recent years, to numerically analyze the coupling problem of solid-liquid two phases [183][184][185][186][187][188][189].
Both SPH and DEM are meshless methods, which have certain advantages in solving problems such as large deformation. SPH is a pure Lagrangian meshless numerical method and suitable for large deformation, free surface flow, and motion interface. Combining with the DEM, it can be used to analyze the deformation and the fracture during brittle material processing [11]. For now, the practical application of SPH-DEM coupling technology is mostly dealing with landslide problems [190]. In terms of simulating the machining of hard-brittle materials, Wu et al. [191] have studied the rock-crushing behavior under the impact of a water jet based on the coupling of SPH and FEM/DEM, and systematically studied the influence of rock microstructure and microscopic performance on the rockcrushing performance under the water jet impact. This example has shown the advantage of the DEM-SPH coupling in dealing with the problems of complex fluid-solid coupling process.

Comparison of DEM with other numerical methods
Compared with other numerical methods, the biggest advantage of the DEM is that it can naturally simulate the initiation and propagation of cracks in hard-brittle materials, without the need for other additional assumptions. For a detailed description and discussion, please refer to Section 1. However, the DEM still has some shortcomings, such as the difficulty of calibrating the microscopic parameters of the material, the lack of sufficient contact constitutive models, and the limitation of the simulation system due to the calculation capabilities. In order to facilitate readers to choose different simulation methods reasonably, different numerical simulation methods are compared in detail, as shown in Table 9.

DEM simulation with multi-field coupling
In reality engineering, there are many physical fields, such as stress field, temperature field, humidity field, etc. In the current research of ceramic material machining, it has not only involved a single physical field, such as the coupling problem of flow field and stress field has become a common concern in the engineering community. The interaction between the stress field and the seepage field constitutes the coupling relationship between the stress field and the seepage field. Liu et al. [199] established a dual medium seepage stress coupling based on continuum discrete elements. At the same time, in the LAM and EDM machining of ceramic materials, the problem of coupling the temperature field and the stress field is encountered. In the research of the grinding process, the grinding heat is also used to remove the material and the surface quality, so it involves the coupling problem between temperature field and stress field [173,200].

DEM simulation with multiphase medium
When simulating the processing of ceramic materials, multiphase coupling problems are often involved, such as gas-solid coupling and fluid-structure coupling. For the fluid-solid coupling problem, the overlapping range of solid-phase and liquid-phase interactions can be divided into two categories. The first major category is that the two phases overlap partially or completely; the second major category is that the coupling of two phases occurs at the interface of the two phases. For the second major problem, the DEM no longer has an advantage because its solid phase medium is a continuous whole. In fluid-solid systems, solids are modeled by discrete circular (or spherical) particles, and the fluid is described by the continuity equation and the Navier-Stokes equation. The coupling method is often applied to the pneumatic transportation process of bulk materials, the powder and polymer dispersion mixing process, the geotechnical engineering, and the agricultural mechanical engineering. In the machining simulation of ceramic materials, the processing involving fluidsolid coupling includes chemical mechanical polishing (CMP) and particle water jet (PWJ). The simulation of this multi-phase coupling problem cannot be completed simply by discrete element simulation technology but, usually, coupling with other numerical simulation software, such as DEM-SPH can simulate the gas-solid coupling problem, DEM-CFD can simulate the solid-liquid coupling problem.

Challenges and opportunities
After several decades of development, the numerical methods on researching brittle materials have evolved from FEM and/ or DEM to the coupling of numerical methods such as FEM-DEM, CED-DEM, SPH-DEM, BEM-DEM, etc. The problem is also developed from the simulation of static problems to the simulation of dynamic problems, and the complexity of the simulation has also evolved from pure mechanical simulation to simulation of multi-phase and multi-field coupling problems. The successful application of the above numerical methods provides useful help for the study of material removal mechanism, crack propagation mechanism, and optimization of processing parameters during the processing of ceramic materials. However, compared with the actual processing, there are still many computing challenges and opportunities for further development. The authors believe that the method based on particle simulation for the machining process of brittle materials can be improved or optimized from the following aspects.

Selection and improvement of material constitutive model
For hard-brittle materials, the mechanical properties are much more complicated than general engineering materials. In most cases, we regard it as a continuous medium with pure brittleness, but defects in the material, grain boundary orientation, or slip surface of the rock often affect its mechanical properties. Using of the BPM model and CBM model described in Section 2.4 can effectively simulate the crack propagation and failure of hard-brittle materials under complex conditions, such as withstand the external loads or heat effects. However, the parameters describing the constitutive relationship between the particles are very limited, and the particles are usually assumed to be ideally rigid. How to invert complex and rich material mechanical properties through limited contact constitutive parameters is a huge challenge for researchers. In addition, the current contact constitutive model in DEM simulation basically only considers the mechanical response of the material under static or quasi-static (including the calibration process). But in the processing of hard-brittle materials, the strain rate of the material is usually up to 10 2 -10 7 s −1 [201]. Therefore, how to improve the constitutive model of the material and calibrate its dynamic response through experiments such as Split Hopkins Pressure Bar (SHPB) [202] may be a direction worth studying in the future of DEM numerical simulation processing of hard-brittle materials.

Calibration of microscopic parameters of the DEM model
In the particle-based method, it is usually needed to build a simulation model of the workpiece or the abrasive based on the accumulation and bonding of the particles. However, the attributed parameters among the particles belong to the micro level, how to construct the correlation between the micro parameters of the particle and the macroscopic mechanical properties of the workpiece or the abrasive tool is one of the questions that researchers should be worth considering. Currently, this is still a black box issue. In the existing research, the mechanical properties of the workpiece or the abrasive tool were mainly inferred by the trial and error method. However, this method has problems such as cumbersome matching process, lack of theoretical support, and non-unique combination of micro parameters, etc. Nowadays, the intervention of artificial intelligence, big data, and computational reversal will  [7] Menezes et al. [192] Hu et al. [193] MD • Describe the interaction between materials from the viewpoint of molecular or atomic • Suitable for ultra-precision machining simulation • Small simulation scale • Difficult to establish the potential functions Meng et al. [13] Shockly et al. [14] Abdulkadir et al. [15] SPH • Not limited by scale • Can simulate large deformation problems • Particle defects at free boundaries are difficult to guarantee accuracy • Tension is unstable Lv et al. [11] Mi et al. [194] Cao et al. [195] BEM • Effectively reduce the spatial dimension of solving problems Bencheikh et al. [196] Zhai et al. [197] Dong et al. [198] help to quickly and efficiently reverse the microscopic parameters among particles.

Size effect
For brittle materials, the size effect is their basic mechanical properties. The mechanical properties of brittle materials will change with the changing of the geometric dimensions. Generally, the strength and deformation characteristics of the brittle material specimens obtained by the numerical or experimental tests under certain dimensions will be different from the strength and deformation characteristics under actual service conditions. Therefore, the research on the size effect of brittle materials is of great significance to guide the engineering practice.
In the study of the size effect, generally the structure is similar and the size of the sample or model is constantly changed to study the change of its mechanical properties. Bazant et al. [203] put forward a lot of theories about the strength and size effect of brittle materials, and pointed out that there is a size effect when using the DEM to study the mechanical properties of brittle materials. Han et al. [204] found that the elastic modulus of single-crystal silicon has no obvious size effect through the micro-bridge experiment of single-crystal silicon, while the bending strength has the theory of obvious size effect. The authors [205] used DEM to study the mechanical properties of single crystal silicon. The results show that the size effects of Poisson's ratio, compressive strength, elastic modulus and fracture toughness are not obvious, while the bending strength decreases with the increase of the model size small. Di and Xue et al. [206] found that the macro-mechanical properties decreased with the increase of particle size and sample size D/L when they studied the sea-ice sample.
In order to clarify the effect of size effect, some scholars explained it from the perspective of energy absorption and release [207,208]. They believe that for brittle materials such as rocks and ceramics, there are more or less original defects such as cracks or holes in the interior of materials, which leads to anisotropy in mechanical properties and non-uniformity in material properties. When an external force is applied, a very complex internal stress field will be formed inside the sample, and the internal stress field will be released by the boundary surface. If the sample is large, this release effect can be ignored. However, if the sample is small to a certain extent, this release effect is very obvious. This results in a situation where the strength of the material decreases as the volume of the test piece increases, and the elastic modulus increases with the volume of the test piece. In other words, in order to reduce the influence of sample size and particle size on the accuracy of simulation calculation, it is recommended to study the size effect of the material before using DEM numerical simulation, and then choose the sample size and particle size reasonably according to the simulation results to reduce or eliminate the effects of size effects.

High-performance computing
Due to the limit of computational ability, most particle-based method simulations currently have a particle size of no more than one million. The simulation system constructed is much smaller than the actual physical system, or the time span of the simulation is much shorter than the actual processing time. How to improve the scale and system of calculation is another issue that researchers should pay attention to. At present, although parallel processing technology, GPU, and other methods have improved the computational scale and computational efficiency to a certain extent, there is still a big gap with reality. In the future, methods that rely on computational hardware to improve computational efficiency will be more widely used. At the same time, the scaling model simulation based on similarity theory will also provide powerful theoretical support for solving the large computational system problems.

Coupling with other numerical methods
Although DEM has the inherent advantage of simulating crack propagation in simulating ceramic material processing, it is often impossible to satisfy the material removal mechanism under various complex processing conditions by relying on DEM alone. With the development of processing technology, the contact between the abrasive grains and the workpiece is not only related to mechanical problems, but also involves the fluid phase (chip liquid, grinding fluid, polishing fluid, etc.), gas phase (air, multi-phase coupling such as protection of inert gas and electric field (ELID), magnetic field (magnet or heological polishing), and sound field (ultrasonic cutting, ultrasonic grinding), etc. At this condition, the DEM needs to be coupled with other numerical methods, especially with numerical methods such as CFD, SPH, and FEM. For example, when analyzing the chemical polishing and abrasive flow processing of brittle materials, it involves the problem of fluid-solid coupling. At this time, it is suitable to use the coupling method such as DEM-CFD or DEM-SPH; when the problem is equal, the DEM-FEM coupling method can be used [173].
It is hard to satisfactorily satisfy the processing requirements of ceramic materials by a single machining method or special processing method. With the development of ceramic processing theory and equipment, the comprehensive use of mechanical processing and advanced laser, EDM, ultrasonic, microwave, and other composite processing technologies to achieve high-quality, high-efficiency, and low-damage processing of engineering ceramics is the inevitable trend of engineering ceramic processing technology development. In the process of composite processing of ceramic materials, multifield coupling, such as thermo-mechanical coupling, is needed to further study the dynamic behavior of ceramic materials under transient, local high stress, high strain and the deepening of DEM research, the DEM will also exert its greater advantages in simulating the machining process of brittle materials. By means of DEM and other particle-based numerical simulation techniques, it is useful to reveal the scientific essence of the ceramics material removal process under multifield, and it also can provide a theoretical basis for the development of processing techniques and equipment for hard and brittle materials. With the deepening of DEM research, the DEM will also exert its greater advantages in simulating the machining process of brittle materials.

Conclusions
In this work, we first introduce the advantages and disadvantages of several numerical simulation methods used commonly, and highlight the advantages of DEM in the study of discontinuous media problems. And then the basic principles of DEM, contact models, the coupling of other numerical methods were introduced, and how to construct and calibrate the DEM model of brittle materials was also discussed in this paper. Finally, the applications of the DEM in the machining processes (such as cutting, milling, grinding, etc.) of brittle materials were listed. Some of the conclusions of this review are summarized as follows: (1) The particle-bonded model is commonly used in the construction of hard-brittle materials, which can easily simulate the crack propagation process, but the microscopic contact parameter model of the material needs to be calibrated by trial and error or parameter optimization. (2) During the calibration process, the compression tests are usually used to calibrate the shear microscopic parameters. Concurrently, the tensile test, Brazilian test, and bending test can be used to calibrate the tensile microscopic parameters. (3) In the process of machining simulation, the tool model can use the wall-based method or the particle-based method. If the damage and deformation of the tool are not considered, the wall-based method is more efficient. Meanwhile, if the tool wear or damage needs to be considered, the particle-based method is recommended. (4) Although the DEM method can be used to simulate the damage of tools or workpieces in the processing of hardbrittle materials, due to the limitations of the single numerical method, the DEM method can be considered to be coupled with other methods. And the influence of various working conditions can be considered in more detail to improve the calculation accuracy and efficiency of simulation.