Peridynamic simulations of damage in indentation and scratching of 3C-SiC

The cubic silicon carbide (3C-SiC) has broad application prospects due to its excellent properties. The modeling of damage in 3C-SiC due to contact loads is important yet challenging. In this paper, simulations based on peridynamics theory are proposed to model damage of 3C-SiC in indentation and scratching. In indentation, the pop-in phenomenon of load-depth curve and initiation and propagation of cracks are observed. During scratching, the specific cutting energy (energy consumed per unit volume of material removed) increases nonlinearly with the decrease of scratching depth. Initiation and propagation of radical cracks on both sides of the groove are found. In addition, for scratching with double indenters, the volume of material removal is more than twice that of scratching with an indenter when indenter interval is less than the threshold. This paper demonstrates that peridynamics is a powerful tool to investigate damage in indentation and scratching of brittle materials.


Introduction
The cubic silicon carbide (3C-SiC) has broad application prospects in many fields such as microelectromechanical systems (MEMS), laser, detection, and other optoelectronic devices due to its excellent mechanical, chemical, electrical, and thermal properties [1][2][3]. In most applications, SiC materials have to undergo contact loads inevitably, which will lead to damage such as cracks in SiC. Since damage usually deteriorates the mechanical properties of SiC devices [4], it is fundamentally important to understand the damage characteristics in the indentation and scratching process of 3C-SiC.
At present, a great number of experiments [5][6][7][8] and simulations [9][10][11][12][13] have been carried out to investigate the deformation mechanisms of 3C-SiC. For instance, Mishra et al. performed large-scale molecular dynamics (MD) simulations of nanoindentation and found that the primary deformation mechanisms of 3C-SiC are dislocation nucleation and propagation in the low pressure (zincblende) phase [11]. Goel et al. investigated the influence of crystal anisotropy of 3C-SiC during its cutting behavior by MD simulations [12]. Liao et al. employed MD simulations to analyze the elastic-plastic deformation mechanism of 3C-SiC ceramic parts with hemispherical and regular pyramid indenters under ultimate loads [13]. Although great insights into the deformation and damage of 3C-SiC have been gained using MD simulations, MD simulations usually fail to capture the spatial and temporal scales in experiments due to the high computational cost. Therefore, many simulation methods based on classical continuum mechanics (CCM) such as finite element method (FEM) [14] and smooth particle hydrodynamics (SPH) have been adopted to study the indentation and scratching of SiC [15]. For example, Zhao et al. employed FEM and experiments to reveal the coexistence of microscopic plastic deformation and brittle fracture of 3C-SiC at different indentation depths [16]. The scratching process of SiC was also investigated using SPH [17,18]. However, the FEM and SPH are based on continuum approaches and cannot properly deal with the discontinuity problems such as crack initiation and propagation occurring in the indentation and scratching of SiC [19].
Given the limitations of existing methods for simulating discontinuity problems such as crack initiation and propagation, peridynamics (PD) as a new nonlocal theory developed by Silling et al. [20] provides a more refined approach for materials damage analysis [21]. The PD employs the integral equations of motion instead of the conventional spatial differential equations [20][21][22][23][24][25], thus overcoming the challenges of modeling the discontinuity problems with spatial differential equations used in CCM. In fact, PD theory has been widely used in various fields to deal with discontinuity problems, including quasi-static fracture problems [26], dynamic fracture of brittle materials [27,28], composites cutting [29,30] and impact [21,31], fatigue cracking [32], heat conduction phenomena [33,34], etc. For example, Cao et al. used PD simulations to describe the nanoindentation of an archetypical soda-lime silicate window glass [35]. Sayna et al. [36] adopted two-dimensional state-based PD model to study the nanoscale friction and wear behaviors of thin amorphous carbon films. However, so far little attention has been paid to the scratching process, especially for SiC materials. The ability of PD to properly describe the response to indentation and scratching of 3C-SiC has been seldom explored. Therefore, in this paper, 3C-SiC indentation and scratching simulations are carried out using PD method. In indentation, the observed crack initiation and propagation are consistent with the crack morphology in the indentation experiment [37]. In scratching, the variation trends of specific cutting energy with the decrease of scratching depth and the coupling effect of scratching with double indenters on the volume of material removal are in good agreement with the experimental results [38]. The work demonstrates that PD method is a powerful tool to investigate the indentation and scratching process of 3C-SiC. It can also provide a reference for the modeling of damage evolution of other brittle materials.

Indentation simulations
The load-depth curves of the PD indentation simulation at different indentation speeds are shown in Fig. 1. For the indentation speeds of 5 m/s, 10 m/s, and 20 m/s, the load-depth curves in the loading process almost completely overlap. However, the load exhibits a certain degree of decrease when the indentation speed decreases to 1 m/s. In the unloading process, the load decreases faster for a higher indentation speed. The "pop-in" phenomenon is observed at almost the same indentation depth. This indicates that the material reaches its yield limit at a certain indentation depth, and the burst of plastic deformation in MD simulations [39,40] is also reproduced in the PD simulation.
The damage of substrate during indentation is almost the same for speeds of 5 m/s, 10 m/s, and 20 m/s. Therefore, the simulation results with a speed of 5 m/s are taken as the representative at these three speeds and compared with those at 1 m/s. The top view and sectional view of the 3C-SiC are shown in Figs. 2 and 3, respectively. It can be seen that for a small indentation depth, the damage area is small and the deformation is mainly elastic. As the indentation depth becomes larger, the deformation and damage of SiC material increase. When the deformation reaches the elastic limit of the SiC, plastic deformation occurs. And the damage area of SiC material increases significantly, accompanied by initiation and propagation of cracks. When the indentation speed is 1 m/s, initiation of cracks is observed and the damage is relatively small at the indentation depth of 15 μm. When the indentation speed reaches 5 m/s, radical cracks, lateral cracks and median cracks are observed and propagate in the substrate at the indentation depth of 15 μm. Four radical cracks with different lengths are found in Fig. 3(c) and (f). The indentation morphology and crack expansion observed in PD simulations agree well with those from experiments [37], which indicates the effectiveness of the PD indentation model.

Scratching with a single indenter
The scratching process is first simulated with a single indenter. The friction force and normal force in the scratching process at different scratching depths are shown in Fig. 4. It can be found that the variations of friction force and normal force have a similar trend. The maximum friction force and normal force increase with the increase of scratching depth. At the scratching depths of 1 μm and 2 μm, the friction force and normal force show slight variations in the scratching process. However, at the scratching depth of 3 µm and 4 µm, the friction force and normal force showed cyclic fluctuations during the scratching process. The big fluctuations of friction and normal force during scratching [19] are similar to the variations of the cutting forces in the cutting simulations using PD. In addition, the wear volumes during the scratching process at different scratching depths are also calculated, as shown in Fig. S-2. The wear volume increases with the scratching depth. By calculating the total damage value of substrate, it can be found that the occurrence of wear is consistent with the variation of the friction force curves, as shown The damage contours of the substrate at different scratching depths are also shown in Figs. 5 and 6. With the increase of scratching depth, the number of severely damaged particles gradually becomes greater and the damage range expands. It can be seen that at a scratching depth of 1 μm, the width and depth of the scratch are small and insignificant, see Fig. 5(a, c). When the scratching depth increases to 2 μm, an obvious groove appears, see Fig. 5(b, d). No cracks are observed around the scratch groove for the scratching depths of 1 μm and 2 μm, indicating the ductile mode of material removal. When the scratching depth increases to 3 μm, the size of the groove continues to increase and initiation of radical cracks is found, as shown in Fig. 6(a, c). For the scratching depth of 4 μm, the scratch width further increases and radial cracks are distributed on the edge of the scratch groove, as shown in Fig. 6(b, d). The formation of cracks on both sides of the scratch groove during the scratching process at the scratching depth of 4 μm is shown in Fig. 7. As the indenter moves along the -x direction, a damage area appears in front of the indenter [ Fig. 7 To further study the effect of scratching depth, the specific cutting energy defined as the energy consumed per unit volume of material removed is obtained, as shown in Fig. S-3. It can be seen that the specific cutting energy increases nonlinearly as the scratching depths decrease. When the scratching depth is 1 μm, the specific cutting energy fluctuates around 22 GPa, which is much larger than that for a larger scratching depth. In the 3C-SiC scratching experiment, size effect is a main reason for the nonlinear increase of specific cutting energy with the decrease of scratching depth, which is consistent with simulation results [6].
As scratching speed plays a vital role in the scratching process, we further investigate the effect of scratching speed. According to the study above, the material removal behaviors of 3C-SiC in scratching are different at different scratching depths. Therefore, to study the effect of scratching speed, the scratch-  the scratching depth of 1 μm, the scratching speed has almost no effect on the scratching phenomenon. In the case of scratching depth of 3 μm, the variations of friction and normal forces exhibit little difference at the scratching depth of 1 μm in the range of speeds studied. It can be seen that when the scratching speed is reduced to a certain level (1 m/s), the amplitudes of the friction and normal forces decrease, but the average values of forces keep almost fixed. At the same time, the amount of wear and specific cutting energy maintain the same levels in the range of scratching speeds studied.

Scratching with double indenters
As real contact surfaces are always rough, the contact happens at multiple asperities. The multi-asperity interactions play an important role in the tribological process. Moreover, in the actual grinding process, multiple abrasives are present on the grinding wheel. Different abrasives produce different scratches, which may interact with each other. Thus, in this part, PD simulations of scratching with double indenters are performed. It has been found that the interaction between adjacent scratches depends on the scratching depth and separation distance [38]. In this paper, a typical scratching depth of 3 μm is chosen to better compare the results of a single indenter with those of double indenters. The force curves of the PD scratching simulations for double indenters with different indenter intervals of 40 μm, 50 μm, and 60 μm are shown in Fig. S-8. The friction force and normal force for scratching with double indenters show no significant difference in the magnitude of the average values throughout the scratching process at different intervals. Figure 8 shows the damage contours for different indenter intervals. It can be found from Fig. 8(c) that when the interval is 60 μm, two scratches are present and the width of each scratch is equal to that of the single indenter with the same parameters. For the interval of 50 μm, the scratch also shows a clear periodicity and the scratch width is about 92 μm, which is more than twice the width of the single indenter, as shown in Fig. 8(b). When the indenter interval continues to decrease to 40 μm, the scratch width decreases, and the scratch exhibits irregular undulation instead of periodicity, see Fig. 8(a). In other words, for the same scratch depth, the region of interaction between adjacent scratches depends on the indenter interval. The area of interaction between adjacent scratches increases first and then decreases as the indenter interval increases.      The total wear volume at different indenter intervals is also counted, as shown in Fig. S-9(a). It can be seen that when the interval is 60 μm, the total wear volume is almost twice the wear volume of the single indenter. When the interval is 50 μm, there is a significant increase in wear volume compared to the 60 μm interval. When the interval is 40 μm, the wear volume is slightly larger than that of the 60 μm interval and less than that of the 50 μm interval. We further calculate the specific cutting energy at different intervals, as shown in Fig. S-9(b). It can be found that the maximum cutting specific energy at the interval of 60 μm is around 0.6 GPa, which is consistent with that for the single indenter with the same parameters. The scratching for 50 μm interval has the smallest cutting specific energy, which facilitates the material removal and leads to the maximum wear volume. The scratching for 40 μm interval has a slightly lower cutting specific energy compared to that of 60 µm interval, which also results in a slightly larger wear volume. Figure 9 shows the sectional view of the substrate in scratching. It can be found that when the indenter interval is 60 μm, two separate scratches appear, see Fig. 9(c). As shown in Fig. 9(b), when the indenter interval is reduced to 50 μm, the superposition of the displacement field makes the scratches overlap and the damage area increases significantly. As the indenter interval continues to decrease to 40 μm, the scratches still overlap, but the damage area is reduced compared to the case with indenter interval of 50 μm, see Fig. 9(a).
In general, the interaction between adjacent scratches depends on the indenter interval for the certain scratching depth. When the scratching depth is 3 μm and the indenter interval is less than 50 μm, there is an interaction between the two indenters due to the superposition of the displacement field generated in the scratching. The volume of material removed by two scratches is more than twice that of a single scratch. This is attributed to the superposition of the displacement field that can lead to interaction-induced material removal caused by cracks between two scratches and the reduction of specific cutting energy. Moreover, there is a critical interval at which the volume of material removed reaches its peak. This is because the total influence region of the displacement field decreases as the indenter interval becomes rather small. Meanwhile, the interaction of damage between adjacent scratches becomes weaker as the indenter interval is large enough due to the limited diffusion length of the damage [38]. Therefore, when the indenter interval is large enough, the interaction between two scratches becomes insignificant. The results are consistent with the existing experimental results [38].

Conclusions
In this paper, an ordinary state-based PD theory is employed to model the indentation and scratching process of 3C-SiC. The constitutive parameters of 3C-SiC for PD simulation are obtained by performing MD simulations of nanoindentation. The PD simulation results are in good agreement with experimental results. The main conclusions are drawn as follows: (1) The PD indentation simulations reproduce the "pop-in" phenomenon of load-depth curve in MD simulations as well as the crack initiation and propagation in the experiments. (2) During the scratching process, the variation of friction is consistent with that of material wear. We observe no cracks in the 3C-SiC for a small scratching depth. However, initiation and propagation of cracks on both sides of the groove are found in scratching. Moreover, the specific cutting energy increases nonlinearly as the scratching depth decreases due to size effect. (3) For the scratching with double indenters, the material removal is greatly affected by the indenter interval for the certain scratching depth. The volume of material removal by two indenters in the scratching is more than twice that of one indenter due to the coupling effect of the two indenters when the indenter interval is less than the threshold.
In this paper, the implementation of PD theory for the indentation and scratching simulations of a typical brittle solid material (3C-SiC) has been carried out. This study provides detailed parameters for PD indentation and scratching models. The methods for constructing PD models of indentation and scratching developed in this paper can be extended to other brittle materials. Overall, this work demonstrates that PD is a powerful tool to investigate the damage in indentation and scratching process of brittle materials.

Methods
PD theory employs a finite number of particles to discretize a solid. The equations of motion in PD theory are given by Eq. (1) [23]: where ρ is the material mass density, H x is the domain of the spherical horizon with a radius δ, T is the force vector state field, b is the body force density field, t is the time, and dV x ′ is the volume of a particle.
The scalar force state t for ordinary isotropic elastoplastic model is given as Eq. (2) [25]: where K is the bulk modulus, θ is the dilatation, ω is the influence function, e d is the deviatoric extension, and e dp is the plastic deviatoric extension. α is a positive constant and m is the weighted volume which is defined as Eq. (3): where x represents the scalar reference state and · is the dot product of two states.
The value of plastic deviatoric extension can be obtained by plastic flow rule which is defined by Eq. (4): where ψ is the scalar-valued plastic potential function, is the proportionality constant, and ∇ d ψ is the constrained Frechet derivative of ψ in the space of co-deviatoric force states.
A PD elastoplastic material model is proposed by Silling et al. [23,25] for perfect plasticity model and constant influence function ω = 1 . According to the proposed Silling's elastic-plastic material model, the plastic potential function can be shown as a function of the norm of the deviatoric force state as Eq. (5): In addition, yield function f is defined as Eq. (6): where ψ 0 can be determined by calculating ψ for simple shear loading and equating ψ = ψ 0 where yielding occurs. For this condition, ψ 0 can be calculated by Eq. (7): (4) e dp = ∇ d ψ, .
where σ y is the tensile yield stress in the case of uniaxial tension. A short-range force model has been introduced to describe the interaction between indenter and substrate [24]. The contact force is independent of the initial distance between material particles, only depending on the current relative positions of particles. The short-range force density between material particles p and i is defined by Eq. (8): where c sh = 15c, and c = 12E/ πδ 4 is a peridynamic modulus. d pi is the critical distance to determine the contact.
where x p and r p are the position and radius of particle p in the vicinity of particle i, respectively. r i is the radius of particle i, which is set equal to one-half of the grid size.
Considering the Coulomb law, the constitutive model for friction is defined by Eq. (10) [39,40]: where ε fric is a regularization parameter with dimensions of speed (indicating the maximum slip speed that is allowed in the static friction regime), and c fric is the coefficient of friction. During the simulation, the friction force can be obtained by considering the contact force and friction coefficient.
Local damage at a point is defined as the weighted ratio of the number of eliminated interactions to the total number of initial interactions of a material particle with its family members. To describe the degree of the local damage of material particles, by using statistical methods, the damage D at a point x is defined by Eq. (11) [24,41]: where µ is a scalar-valued function and also a discontinuous function to characterize whether the bond between two material particles is broken. When the fracture occurs, µ is 0; otherwise, µ equals 1. The local damage at a point can be quantified as a value which ranges from zero to one [24]. When the local damage is one, all the interactions initially associated with the point have been eliminated. While a local damage of zero means that all interactions are intact, the measure of local damage is an indicator of possible crack formation within a body.
In this paper, the open-source software Peridigm is employed for PD simulations [42].

PD indentation
The indentation model consists of a diamond indenter and a 3C-SiC substrate. An ordinary state-based PD model of 3C-SiC substrate with dimensions of 297 × 297 × 297 μm 3 is considered, as shown in Fig. 1(a). The interparticle distance is 3.3 μm. This study intends to observe the damage evolution in the indentation and scratching process. Brittle materials such as 3C-SiC experience elastic-plastic deformation when subjected to stress [8,13]. Therefore, the effect of plastic deformation should be considered in addition to brittle fracture during 3C-SiC indentation and scratching. The elastoplastic constitutive model of 3C-SiC was parameterized based on existing experimental data. Specifically, the density of 3C-SiC substrate is set to 3200 kg/m 3 and Poisson's ratio is set to 0.168 [43]. Young's modulus used herein is 320 GPa, which is obtained by performing MD simulations, as shown in Fig. S [44]. Similar to the previous studies [45], the horizon δ is 3.05 times the interparticle distance and the critical stretch is set to be 5 × 10 -4 . The particles in the bottom are fixed to serve as the boundary to anchor the 3C-SiC substrate. The timestep is set as 50 ps, which is small enough to ensure the stability of the simulations based on stability analysis [24].
At the beginning of the indentation simulation, the indenter is above the substrate. Then, the indenter indents with a fixed speed (loading) until the maximum indentation depth (15 μm) is reached. Subsequently, the indenter moves upward with the same fixed speed (unloading). The load is computed as the resultant force applied by the indenter on the substrate. Specifically, all the individual components of force applied on the indenter at a given moment are recorded and then added up to calculate the resultant force. Different indentation speeds of 1 m/s, 5 m/s, 10 m/s, and 20 m/s are also employed to further study the effect of indentation speed on the deformation mechanism of materials.

PD scratching
The material and key parameters of the PD scratching model are the same as those of the indentation model. The 3C-SiC scratching model with a size of 594 × 198 × 198 μm 3 is established. The indenter parameters remain the same as those in the indentation, as shown in Fig. 10(b). The indenter indents to a certain depth of the substrate and then scratches the substrate along the negative x direction with a fixed speed. The coefficient of friction between the substrate and the indenter is set as 0.15 [46]. The total scratching distance is 350 μm.
To quantitatively analyze the scratching process, the friction force, wear volume, and specific cutting energy during scratching are calculated. The friction force f is obtained by summing up the individual forces on the indenter particles along the x direction during the simulation. The wear volume is defined as ∑D·V, where D and V are the damage and volume of the particles, respectively. The power during scratching is calculated by W = f·v, and the trapezoidal integral method is used to integrate the time to obtain the cutting energy. Combined with the wear volume obtained previously, the specific cutting energy can be obtained by Eq. (12): To explore the influence of scratching depth, the scratching depths of 1 μm, 2 μm, and 3 μm are selected for the simulations.  To investigate the effect of scratching speed, different scratching speeds of 1 m/s, 10 m/s, 20 m/s, and 40 m/s are chosen for two typical scratching depths of 1 μm and 3 μm. The speeds used in the simulations are constrained by the computational cost during the simulations, which are typically larger than the real speeds in the experiments. It is noted that even though higher scratching speeds are used in the simulations compared to the experiments, the scratching speed in this range has little influence on the simulation results, as shown in Figs. S-4-S-7.

Author contributions
YX contributed toward methodology, software, data curation, formal analysis, writing-original draft, and visualization. PZ contributed toward conceptualization, validation, resources, writing-review and editing, supervision, project administration, and funding acquisition.

Data availability
The raw/processed data required to reproduce these findings cannot be shared at this time as the data also form part of an ongoing study.

Conflict of interest
The authors declare that they have no conflict of interest.