Investigation of parameters influencing erosive wear using DEM

The effects of several parameters on the erosive wear were studied using the discrete element method (DEM). The Finnie model was implemented using an open-source code. Regarding the time integration, the Gear algorithm was used, and to ensure the accuracy of the DEM results, a time-step sensitivity analysis was performed. The problem was modeled in two parts: first, the impact of a single particle on a surface was modeled, and then a more general model was prepared to examine the wear of surfaces under the flow of particles. The effects of the surface area, impact angle, speed, particle size, particle density, Young’s modulus, Poisson’s ratio, and restitution coefficient on the wear were studied numerically, and the results are discussed herein.


Introduction
Wear is among the important issues from the theoretical, practical, economic, and environmental viewpoints. For instance, approximately 3% of the world's total energy is used to remanufacture parts that fail owing to wear [1]. Therefore, there are sufficient motivations for scientific and research investment in this field. Because of the complexity of the issue, it is difficult to be modeled theoretically and predicted practically. No consensus exists in the literature regarding the classification of wear. For instance, in Ref. [2], the wear is categorized into four main types-abrasive, adhesive, corrosive, and surface fatigue-and the others, such as erosion, are classified as minor types. As another tactic in Ref. [3], the surface state, interaction mechanism, and relative motion were considered as three main indices for categorizing wear, and a unified wear classification was proposed [3]. A literature review indicates that several studies have been conducted in this field, especially using experimental methods. The experimental study performed by Archard [4] is fundamental research that serves as a cornerstone to wear analysis. Lynn et al. [5] employed an experimental approach for understanding the effects of the particle size inside a slurry on the flow erosive wear. They concluded that the wear rate is related to the kinetic energy of the particles. They proved that if the kinetic energy exceeds a threshold, it causes surface failure and wear [5]. Chacon-Nava et al. [6] used a rubber wheel tester to study the effect of the surface hardness and particle size on the abrasive wear. In the particle-size range that they tested, smaller particles led to larger wear. Such behavior is linked to the cutting effect of small particles, while large particles cause plastic deformation [6]. According to the results of Ojala et al. [7], the particle size affects the stress level, but no absolute conclusion can be drawn regarding the effects of the particle size. In this regard, they found that when the particle size was less than 1 mm, a low-stress condition occurred, and the tested elastomers showed better wear resistance than the tested steels [7]. According to their results, the dependency of the wear on the particle size can exhibit both ascending and descending trends [7]. Woldman et al. [8] studied the effects of particle properties on the wear of parts that operate in sandy environments. Although they showed the importance of the particle size and shape, no implicit relation was observed [8]. Many other studies have employed experimental methods; however, theoretical approaches have also been developed to study wear.
Finnie proposed a mathematical model to predict the wear under the impact of particles [9]. The mass, speed, and impact angle of the particles play roles in this model. Other analytical models, such as those of Archard [10], Wellinger et al. [11], Hutchings [12], and Rabinowicz [13], have also been proposed. The finite-element method has been employed as a conventional numerical method to model wear; however, this method has generally been applied to study the abrasive and adhesive wear mechanisms [14−17]. On the other hand, the discrete element method (DEM) is becoming favorable in this field, as it is able to apply the introduced mathematical models under complex conditions. The study of Cleary [18] is among the first published reference in the field of using DEM for wear analysis and focused on mill liners. Kalala et al. [19,20] studied the wear of mill liners using the DEM. They analyzed the effects of the liner profile on the wear regime. Ashrafizadeh et al. [21] used the DEM to study the effects of the impact angle and particle energy on the wear of a surface. They showed that the maximum shear impact energy occurs at an impact angle of approximately 30° and concluded that this angle causes the highest wear. However, we believe  [22] and Tan et al. [23] used the DEM to study the erosive wear in pipes carrying a slurry. They showed that the effect of the speed on the wear rate in the elbow region was larger than that in the straight parts and determined the critical region of the elbow. On the same topic, Uzi et al. [24] explored variations of the wear in different regions of the elbow in conveying pipelines. Powell et al. [25] used the DEM to determine the impact energy of particles in a mill. They estimated the profile of the liners undergoing wear gradually. Varga et al. [26] employed both an experimental approach and the DEM to evaluate the wear in pipes carrying a slurry. As another case study, Jafari et al. employed the DEM to examine the effects of the vibration screen characteristics on its mesh wear [27]. Chu et al. [28] used the DEM to predict the wear in dense medium cyclones (DMCs) and proposed this approach for studying effects of the wear on the DMC performance. Forsström et al. [29] employed the DEM with Archard's wear model, finding that the critical region of tipper bodies ruptured owing to abrasive wear. They concluded that such numerical simulations are useful for optimizing the tipper body geometry in order to avoid local failure and improve its service life. According to the foregoing literature review, the DEM has the potential to model wear; nonetheless, this field is not fully developed, and further investigations can be valuable for enhancing the modeling of the wear under realistic conditions. Moreover, this can assist researchers to reduce the requirements of experimental setups, which are generally expensive and take considerable time. In this regard, the present study employs the DEM to investigate the parameters that affect the wear following two procedures. In the first part, a single-particle model is used to observe the basics of the wear due to the impact of a particle, and in the second part, simulations inspired by a mixture-type experimental setup are performed.

DEM
According to the DEM, motion equations are solved to determine the position, velocity, and acceleration of all the particles in the system. For a system of N particles, the equilibrium equations of the force and moment for the i th particle are as follows:   is the angular velocity, i m is the mass, and i J is the moment of inertia of the i th particle.
are the resultant force and moment on the i th particle, respectively, and are determined as follows, 1, Here, ij F  and ij M  are the interacting force and moment, respectively, acting from particle or element j on particle i. The interaction force arises from direct contact between particles and contact between particles and surface elements. On the other hand, the moment is due to the action of the tangential force on the particle. The equation (s) for the estimation of the tangential force will be provided subsequently. Several models have been proposed for estimating the interaction forces between contacting particles [27]. The viscoelastic Kelvin-Hertz model employed here relates the normal interaction force to the normal deformation and deformation rate, as follows [30]. According to this model, the interaction force depends on normal deformation and its rate, and no explicit relationship with the orientation is considered in Eq. (5), Here, ij  is the total normal deformation of collided particles i and j at their contact point. * ij R and * ij E are the equivalent radius and Young's modulus, respectively, which are defined as follows,  is the internal viscosity related to the restitution coefficient, ij  , via the following Eq. (7) [30, 31], Here, * ij m is the equivalent mass of the interacting particles i and j, and n ij K is the normal equivalent spring between these particles, Furthermore, the tangential interaction force between two bodies is obtained as follows, where ij t v , ij t  , and ij  are the relative tangential velocity and the viscous and dry friction coefficients between the interacting bodies, respectively. The governing equations in the system are time-dependent, and numerical integration with respect to time is performed in process of obtaining the numerical solution. Here, the Gear integration method [32] is employed for this purpose, and sensitivity analysis is performed to ensure adequate accuracy.
To estimate the wear using the DEM, the Finnie model was implemented [33]. For this purpose, each surface was meshed into small triangular elements, and the wear of each element is evaluated as a summation of the wear due to the impact of all the particles on that element. During the process of solving the DEM, the position and velocity of all the particles and surface elements are known at each time step. Therefore, the following equation can be employed to evaluate the wear due to each impact at any time step [9,33]: where W is the removed volume from a point due to the impact of a particle with mass m, speed ν, and impact angle  (Fig. 1). K F is the Finnie wear constant, and ( ) Z  represents the dependency on the impact angle, according to the following Eqs. (12a) and (12b) [9,33], Here, k is a conditional constant depending on the material characteristics and is related to the impact angle that causes the maximum wear. A few experimental tests are needed to determine the constants k and F K in Eqs. (11) and (12). The experimental data should match the mathematical model. However, in theoretical studies, including the numerical simulations of the present study, values are assigned to these constants. In this regard, assumed values and relevant explanations are presented in Section 3. In the DEM, the removed volume is divided by the affected surface to obtain the depth of the wear, and the unit of the wear becomes the depth of the wear. To clarify the process of wear estimation via the DEM, the contact time can be divided into the penetration and repulsion periods. During each of these two periods, the magnitude and direction of the particle velocity change gradually. Thus, the total wear is determined via integration of the wear rate over the contact time. For this purpose, the wear rate is calculated during each impact by using the following Eq. (13), It is a logical assumption that during the penetration  | https://mc03.manuscriptcentral.com/friction period, the particle transfers its energy to the surface, causing wear, and negligible wear occurs during the repulsion. Therefore, integration of the wear rate over the penetration period determines the wear. On the other hand, the direction of the particle velocity is assumed to remain almost unchanged during the penetration period; thus, the last term in Eq. (13) can be neglected. Hence, by substituting Eq. (11) into Eq. (13) and employing the well-known equilibrium Here, F is the contact force. Finally, the wear due to a collision is obtained via integration of Eq. (14) over the penetration time, as follows, Here, ( , ) H v q   is the Heaviside function defined below, which applies the rule that the integration kernel is non-zero only during the penetration period, Here, q  is a unit vector pointing from the particle center to the contact point.
To conduct the simulations, adequate codes were developed by using the open-source DEM code LIGGGHTS [33], which is a special version of LAMMPS for modeling granular media. A set of parallel computers with a total of 24 cores and a Linux operation system was used to run the simulations. Graphical results were obtained using ParaView [34].

Numerical results for single-particle model
As mentioned previously, the Finnie model is employed here to evaluate the erosive wear. According to this model, the particle mass, impact angle, and Finnie coefficient participate in the mathematical formulation implicitly. However, in addition to these parameters, others, including the particle size, surface mesh size, Young's modulus, Poisson's ratio, and restitution coefficient, are studied using a single-particle model. In this part, a spherical particle is dropped at the center of an equilateral triangular surface. The orientation of the surface is changed to provide different impact angles. Table 1 presents the ranges of the parameters investigated in this part. The effect of each parameter is studied individually, meaning that the others are fixed at the reference values presented in Table 1. In the reported results, including the graphs, when the value of any parameter is not given, the reference value should be assumed.
As mentioned previously, the accuracy and convergence of the results depend on the time step. Therefore, sensitivity analysis is performed to obtain a proper time step leading to reliable results. For this purpose, several simulations are conducted by selecting different values for the time step in the range of 0.01−1.84 μs, and the results are plotted in Fig. 2. To provide a clear indication of the sensitivity to the time step, the results in this figure are normalized by dividing them by the smallest one. The dependence of the wear on the time step disappears when the time step is reduced. When the time step is <0.075 μs, the graph almost becomes plateau, which means that acceptable accuracy is obtained regardless of the time step.

Effect of surface mesh size
In the implementation of the DEM, each surface must be meshed into small triangular elements via discretization. To reveal the effects of the meshing on the results, a spherical particle is dropped at the center of an equilateral triangular surface, as shown in Fig. 3. The same problem is simulated under three different meshing cases consisting of 1, 4, and 16 elements. By dividing the data by the smallest one, the normalized wear results are obtained. As shown in Fig. 3, the wear values of the affected element are 1, 4, and 16 in the cases where 1, 4, and 16 elements are used, respectively. In the case where 16 elements are used, the evaluated wear is 4 and 16 times larger than those for meshing with 4 elements and 1 element, respectively. Moreover, the wear occurs only in the element that is impacted by the particle. Thus, the particle transfers its kinetic energy to only the area of the impacted element, rather than the entire area of the surface. However, multiplying the evaluated wear by the area of the affected element gives an identical value for all three meshing cases. Hence, the reported value for the wear is the depth of the wear of the affected element. It can be concluded that a particle impact causes a certain volume to be removed from the target body; however, a smaller affected element yields greater depth of the local wear. That is, the local wear depends on the meshing, such that refinement of the meshing leads to the evaluation of more localized wear. Nevertheless, as explained previously, the overall wear of the entire surface is identical for all meshing numbers.

Effect of impact angle
In the aforementioned Finnie analytical model, the wear depends on the impact angle according to certain functions. However, there is the question of how the DEM code captures this model and the question of their coincidence. As explained in the previous section, the constant k in the Finnie model is a material characteristic; however, as an example case, we assumed that the maximum wear occurs at the impact angle   c 38°, leading to 12 k  . To compare the simulation results with the analytical Finnie model, several simulations were conducted with different values of the impact angle, while the other parameters were fixed at the reference values. The normalized wear (obtained by dividing by the largest value) versus the impact angle is plotted in Fig. 4. According to these graphs, in most ranges of the impact angle, both approaches are in good agreement; however, a deviation between the DEM and the analytical Finnie model is observed. The experimental results indicate that the impact angle that causes the largest wear depends on the material behavior; for ductile materials, this angle is in the range of 20°−60°, whereas for brittle materials, it is approximately 90° [13,35]. Therefore, the current  model is suitable for ductile materials. The deviations between the analytical and numerical results can be due to differences between the implementation of the wear model in the DEM and the analytical approach. In the analytical model, the particle speed at the beginning of the collision is considered, whereas in the DEM, the more complicated scenario illustrated in Section 2.1 is implemented.

Effects of impact velocity
To examine the effect of the impact speed on the wear, simulations were conducted for different values of the impact speed in the range of 2−10 m/s. To properly indicate the dependency, the results were divided by the largest value, and the obtained normalized results are plotted in Fig. 5. As shown in Fig. 5, the DEM results imply that the dependency of the wear on the speed has the form of a power function, such as 1.87 .
v The same trend in the form of n v was observed for other simulation conditions, and the value of the power factor, n, depends on the impact angle. This agrees with previous experimental results, such as those reported in Ref. [36]. Further discussion regarding this matter is omitted here, and adequate explanations will be presented in the next section.

Effect of particle size
According to previously reported experimental studies, increasing the particle size leads to a higher wear [5, 37−39]. Furthermore, it has been proven that there is no unique dependency for all ranges of the particle size [40,41]. No consensus exists in the literature regarding the mathematical model for the dependence of the wear on the particle size. For instance, a polynomial function was presented in [42], and a power function was proposed in Refs. [43−45]. In the present study, for different values of the particle diameter in the range of 0.2−40 mm, DEM simulations were conducted, and the normalized results are shown in Fig. 6. Normalization was performed by dividing the results by the largest one. In Fig. 6, a fitted curve for the DEM results is also plotted. The wear is related to the particle diameter almost in the form of 2.81 d . Thus, for the range of particles studied, the power function is suggested for the dependence of the wear on the particle diameter.

Effect of particle density
According to the Finnie model, the wear is linearly related to the particle mass; however, the DEM model was verified via simulation of several problems with different particle densities in range of 500−5,000 kg/s. The normalized results depicted in Fig. 7 show that the wear has a linear relationship with the particle density. This result demonstrates the idea that the wear of a surface is due to the kinetic energy of the impinging particle. The kinetic energy varies linearly with respect to the particle mass; thus, in the case where the particle volume is constant, this energy has a linear relationship with the particle density. However, other wear estimator models have been proposed for brittle materials that suggest a nonlinear dependence of the wear on the particle density [46,47]. The present results are suitable for ductile materials.

Effects of restitution coefficient, Young's modulus, and Poisson's ratio
According to the Finnie model, there is no implicit relationship between the wear and the restitution coefficient, Young's modulus, and Poisson's ratio. However, the wear depends on the velocity and the impact force, which is related to these characteristics. Thus, several simulations were conducted with different values for each of the introduced properties. Here, the impact angle was 50°, and the other characteristics were fixed at their reference values. For the ranges of the introduced characteristics, various wear results were obtained. They were normalized by dividing them by the largest value for each case study and are plotted in Figs. 8−10. As shown in Fig. 8, the wear is    with the increase of this coefficient, the wear increases slightly. This dependency can be due to the fact that by increasing the restitution coefficient, the dissipated energy is decreased, and the higher remaining kinetic energy can contribute to the increase of the wear. On the other hand, the collision time varies with respect to the restitution coefficient [48] and can affect the integration domain introduced in Eq. (15). This can be another illustration of the dependency of the wear on the restitution coefficient; however, it is not a substantial relationship.
Regarding the dependency of the wear on the elastic constants-the Young's modulus and Poisson's ratio-the results are plotted in Figs. 9 and 10. Similar to the case of the dependency on the restitution coefficient, small variations are observed in these graphs. In the Finnie model, there is no apparent relationship between the wear and these constants. Nevertheless, owing to Eqs. (5) and (9) and as proved in [48], these parameters participate in the overall elastic stiffness of the contact zone and thus affect the collision time and interacting force. Therefore, weak relationships between the wear and elastic constants are expected, which agrees with the findings in Ref. [49]. The graphs show that with the increase of the surface Young's modulus and Poisson's ratio, the wear is increased slightly. To discuss such dependencies, Eqs. (5) and (9) should be reloaded as indicate that by increasing these properties, the contact region becomes stiffer leading to a creation a higher energy for making the wear.

Flow erosion
Up to now, the wear was studied using the singleparticle model to reveal the basic correlations between the wear and influencing parameters. As described in this section, another set of simulations was conducted to study the erosion of a surface under a flow of particles. The flow erosion was studied using the model depicted in Fig. 11, in which two bodies rotate around a fixed axis inside a container filled with spherical particles. This model is similar but not identical to the experimental rigs reported in Refs. [5,45,50].
The physical characteristics of the particles and fixed parameters are presented in Table 2. To examine the effects of the sample orientation and speed on the wear, different values were assigned to these parameters in the ranges shown in Table 3. To clarify the details of the model and the definition of the surface orientation (angle), different parts and a schematic diagram were drawn, as presented in Fig. 12. As shown, the surface angle  , is the angle between the surface and the circular path in the container.
As an example, a contour plot of the flow erosion for a typical case is presented in Fig. 13. It is observed that the wear on each sample surface is not uniform. As proved previously, the impact direction is an important parameter affecting the wear and can cause    such wear regimes. Each surface of the sample is meshed with 128 elements, each of which experiences its own flow state, leading to deviation of the local wear between different points on the surface. To facilitate comparison and numerical study, the total erosion of all the elements for each flat surface was determined as the overall erosive wear of this surface, as presented in the following subsections.

Erosion versus time
Variations of the erosion versus time for different speeds and different sample angles are plotted in Fig. 14. The distance can be determined when the rotation speed, radius of rotation, and elapsed time are known. For all speeds and all angles, the wear varies linearly with respect to time (distance). This trend agrees with experimental results obtained for a steel material [51]. However, zooming in on the initial steps presented in Fig. 15 indicates that initial regions of the graphs are nonlinear. This can be due to the fact that the system started in a static condition and took some time to reach steady-state conditions. Thus, the wear varies non-linearly in transient conditions; however, it gradually approaches a linear trend.

Erosion versus sample angle
The experimental studies proved that the angle that causes maximum wear depends on the material behavior. For ductile materials, the maximum wear occurs at small angles in the range of 20°-60°, whereas for brittle materials, it occurs around 90° [13]. In the present study, the role of the sample angle, which is introduced in Fig. 12, is studied by adjusting the angle  . For a constant traveled distance (14.72 revolutions = 0.000555 m), the graphs of the wear versus the surface orientation are plotted in Fig. 16. For all speeds of the surface inside the particulate medium, the maximum wear occurs at 45°  . Therefore, it can be concluded that the constants assigned to the employed Finnie model are suitable for ductile materials. At angles of 0° and 90°, the wear is too small. This can be due to the fact that when the angle is 0°, the particles slip on the surface, and when the angle is 90°, the particles penetrate the surface and return without considerable wear. The simulation conditions used in this study are similar to the  experimental setup in Ref. [45]. Although the numerical results cannot be compared directly, the trends of the dependency of the wear on the sample orientation in this study are similar to those reported in Ref. [45]. Moreover, we observed that the maximum wear occurred for the sample orientation of 45°  , which is in agreement with the experimental results in Ref. [45].

Erosion versus speed
In the previous section, the effects of the particle speed on the wear were studied using a single-particle model. It was proven that the wear is related to the particle speed via a power function; however, there is a question about this correlation when the sample is under the impact of a particle flow. To answer this question for more general conditions, several problems were simulated for different values of the sample orientation and various speeds. For different sample angles, graphs of the erosion versus the speed are presented in Fig. 17. For identical elapsed time, different speeds lead to different traveled distances; thus, these graphs are plotted for an identical traveled distance. The values of the wear change with respect to the sample angle; however, similar trends are observed with regard to the speed dependency. Furthermore, we attempted to extract a function to represent the speed dependency, as shown in Fig. 18, for a few cases and found that at  all sample angles, the erosion depends on the speed via a power function n Cv . Interestingly, these graphs are similar to the one previously proposed for the single-particle model. Regarding the correlation function, agreement is observed between this form and the results in the literature, especially the experimental studies in Refs. [11, 12, 36, 52−54]. In fact, the constant C and the power n depend on the working conditions and are obtained by fitting experimental/numerical data to the power function. For instance, values in the range of 1.5−3.5 have been reported for the power, n [11, 12, 36, 52−54]. These constants were quantified via interpolation of the DEM results, and their values for different sample angles are presented in Table 4. The power constant, n, is approximately 2.3 for angles of 0°  and 90°  and approximately 1.9 for the other angles.

Conclusions
The DEM was employed successfully to study the surface wear under the impact of granular materials. First, a single-particle model was used to calibrate the wear parameters and identify the basic factors that affect the wear process. It was proven that the impact angle, impact speed, particle size, and particle density significantly affect the wear. Furthermore, it was observed that physical characteristics-the Young's modulus, Poisson's ratio, and restitution coefficientaffect the wear to some extent. As a more general case, simulations were conducted to identify the parameters influencing the surface erosion under the flow of particles. The results showed that in steady-state conditions, the flow erosion increases linearly with respect to time. The orientation of the surface inside the particulate medium is a key factor affecting the wear. It was observed that the erosion is related to the speed