## 2 Numerical Scheme
### Wear Load Calculation Model  
The DEM can be defined as a process of iteratively calculating the particle behavior and impact load. The collision load can be classified into a normal force and a tangential force of a contact pair, and each load is defined as a function of the stiffness and overlap of the contact pair. The normal force was defined using the hysterical linear spring (HLS) model. The HLS model is an elastic-plastic model that is used by separately defining the loading stiffness, in which the direction of motion is maintained in the direction of the initial collision, and the unloading stiffness, where the direction of motion is reversed after the initial collision. This relationship can be schematized as shown in Figure 1 where \(F_{n}\), \(S_{n}\), and \(S_{nres}\) denote the normal force, overlap, and residual overlap, respectively. The slope of \(\overset{=}{AB}\) is the loading stiffness, and that of \(\overset{=}{BC}\) is the unloading stiffness. \(\overset{=}{CA}\), which is an overlap that remains after being completely unloaded numerically, is plastic deformation, and the area of the triangle is energy dissipation due to noise or heat.  
The linear spring Coulomb limit (LSCL) model is used to define the tangential force, the most important load that causes wear. The Coulomb friction model, which defines the tangential force as the product of the friction coefficient (\(\mu\)) and the normal force (\(F_{n}\)), has the disadvantage that the contact area and sliding stiffness cannot be considered. The LSCL  
Figure 1: Overlap–normal force relationship in HLS model.  
model has the advantage that the contact area and sliding stiffness can be considered under the tangential force defined by the Coulomb friction model. The normal force model and tangential force model are defined by the following equations and shown in Figure 2.  
\[F_{n}^{t}=\begin{cases}\min\bigl{(}K_{nl}\,s_{n}^{t},\;F_{n}^{t-\Delta t}+K_{nu} \Delta s_{n}\bigr{)}\;\;if\;\Delta s_{n}\geq 0\\ \max\bigl{(}F_{n}^{t-\Delta t}+K_{nu}\Delta s_{n},\;\lambda K_{nl}\;s_{n}^{t} \bigr{)}if\;\Delta s_{n}<0\end{cases} \tag{1}\]  
\[\Delta s_{n}=s_{n}^{t}-s_{n}^{t-\Delta t} \tag{2}\]  
\[F_{\tau}^{t}=\min\bigl{(}\left|\,F_{\tau,\varepsilon}^{t}\right|,\;\mu F_{n}^ {t}\bigr{)}\frac{F_{\tau,\varepsilon}^{t}}{\left|F_{\tau,\varepsilon}^{t}\right|} \tag{3}\]  
\[F_{\tau,\varepsilon}^{t}=F_{\tau}^{t-\Delta t}-K_{\tau}s_{\tau} \tag{4}\]  
where \(K_{nl}\) and \(K_{nu}\) denote the contact loading stiffness and the contact unloading stiffness, respectively. \(S_{n}^{t}\) is the normal overlap at time \(t\), \(\lambda\) is the dimensionless coefficient, \(\mu\) is the friction coefficient, \(K_{t}\) is the tangential stiffness, and \(S_{t}\) is the tangential overlap. The loading stiffness \(K_{nl}\) and the unloading stiffness \(K_{nu}\) are respectively defined as follows:  
\[\frac{1}{K_{nl}}=\frac{1}{K_{nl,p1}}+\frac{1}{K_{nl,p2}} \tag{5}\]  
\[K_{nu}=\frac{K_{nl}}{\varepsilon^{2}} \tag{6}\]  
\[K_{nl,p}=E_{p}L \tag{7}\]  
where \(K_{nl,p1}\) and \(K_{nl,p2}\) are the stiffness values of particles 1 and 2, respectively, L is the size of the particle, and \(\varepsilon\) is the restitution coefficient, which is not an individual property of particles but a relative property of a collision pair. \(E_{p}\) is the bulk Young's modulus of the particles. \(K_{\tau}\) is the sliding stiffness and is defined as the product of \(r_{K}\) (tangential stiffness ratio), and \(K_{nl}\) as expressed in the following equation:  
\[K_{T}=r_{K}K_{nl} \tag{8}\]