A Stress-State-Dependent Sliding Wear Model for Micro-Scale Contacts

Sliding wear is a common tribological phenomenon for many engineering applications, and it is one of the main causes of the failure of mechanical systems. The underlying physics behind the phenomenon is very complex and, to this day, the topic is still not very well understood. Sliding wear occurs as a result of the interaction of asperities of mating surfaces in relative motion. The real area of contact is significantly smaller than the nominal area of contact. This causes the contact to take place over some asperity summits, which results in plastic deformations. The onset of sliding then causes shearing, which eventually leads to wear. Plastic deformation, fracture of material, thermal softening, and strain hardening are some of the many processes that affect sliding wear. It is known that wear is the result of fracture growth due to large stresses at the contacting asperities between two sliding surfaces [1–3]. The fracture growth strongly depends on the stress state, that is, compressive, tensile, or shear stresses. It is, therefore, important to incorporate this when simulating sliding wear. While there have been many attempts at modeling wear mechanisms, see e.g., [4–20], none of them account for stress-state dependency, and there is a need for a model that considers the material detachment and large deformations that occur during sliding wear.

In our previous work [21], we developed a numerical model based on the finite element method (FEM) and continuum damage mechanics (CDM) to simulate wear of sliding asperities. The main disadvantage of that approach is that the elements were eroded (deleted) when a failure criterion was reached, thus violating the law of conservation of mass. To address the issues in our previous study, we are now proposing a mass-conserving, meshfree, particle method for the prediction of wear. This model also takes the thermal response, resulting from the plastic work, into account. Meshfree particle methods are particularly suited to deal with the large deformations and material detachment that occur during fracture growth and wear. The most common method is to use molecular dynamics (MD), see, e.g., Refs. [22–30]. The main disadvantage of MD is that material properties such as yield stress, strain hardening, and plastic flow, which are known to be of importance when studying sliding wear and the main mechanisms driving it, are difficult to take into account. Another disadvantage of MD is that it is only possible to study very small material volumes during extremely short sliding times. The short sliding times simulated with MD are nowhere near the actual time it takes for two micro-scale asperities to collide.

In a recent article by Dimaki et al. [31], the discrete element method was utilized to model brittle and ductile fracture of sliding asperities. It is known that fracture is dependent on the stress state, stress triaxiality, as well as the Lode parameter [32]. The stress-state dependency of the fracture criterion was, however, not considered in the model presented in Ref. [31], and it is therefore unclear if it can be used to accurately capture the material detachment that occurs during wear.

In the present work, we show that the true meshfree momentum-consistent smoothed particle Galerkin (MC-SPG) method, developed by Pan et al. [33], can be combined with the state-of-the-art generalized incremental stress-state-dependent damage model (GISSMO) [34,35] and a fully coupled thermal solver to simulate the fracture growth due to large deformations in colliding asperities. This approach will thus allow us to simulate the wear and the thermal response of colliding asperities. The method uses an advanced bond failure method to simulate crack propagation in both ductile and brittle materials, and it is suitable for problems where extremely large deformations occur. The main advantage of this method over the coupled finite element and CDM method presented in our previous work [21] is that it is mass conserving and that fracture is modeled as a stress-state-dependent phenomenon.

The damage parameters that are used as input to GISSMO are based on fracture tests for various stress triaxialities and Lode parameters. The GISSMO damage model is particularly well suited for simulating wear and the thermal response of surface asperities because it considers the stress triaxiality, Lode parameter, temperature dependency, as well as the postnecking behavior of the material in question. The main goal of the present work is to establish a numerical model that can be utilized to increase the understanding of the wear process, occurring due to multiphysics phenomena such as plastic deformation, fracture mechanics, material nonlinearity, and thermal behavior. The results obtained by using this model can then, for example, be used to calibrate much simpler and more empirical wear laws, such as the well-known Archard model, as shown in Ref. [21].

In this section, the main idea behind meshfree methods is explained, and a short theoretical background is established.

In areas where large deformations, as well as material failure, are present, the traditional mesh-based FEM can experience difficulties in correctly capturing these processes and is more prone to volumetric and shear locking [36–41]. For this reason, alternative, meshfree methods were developed. Meshfree methods are numerical methods to solve partial differential equations, which nowadays constitute an important branch in the fields of computational fluid and solid mechanics. The main idea behind meshfree methods is to represent a body with a set of discrete nodes. In contrast to the traditional finite element method, the meshfree shape functions are not defined in an isoparametric space and are thus not based on elements. The construction of the meshfree shape functions remains one of the major challenges because it requires a predefined relationship of the nodes. Another challenge is to construct shape functions that satisfy the Kronecker delta property, that is, the value of the shape function of a particle is unity at the particle itself. This property is important because it makes applying Dirichlet boundary conditions straightforward and reduces the computational cost in the traditional FEM [36–41]. If the meshfree shape functions satisfy the Kronecker delta property, the treatment of Dirichlet boundary conditions can be applied in exactly the same way as in traditional FEM. The resulting linear system of equations derived with meshfree methods is similar to that derived with traditional FEM, with the only difference being that the stiffness and mass matrices are calculated with the meshfree shape functions instead of, for example, isoparametric shape functions.

Although meshfree methods do not require a mesh of elements for the field variable interpolation, some methods do still require a background mesh for integrating the system of matrices over the domain. The background mesh is only used for integration purposes, and thus, any form of mesh can be used, provided it can give sufficient integration accuracy.

No matter if the basis function is convex (i.e., nonnegative and reproducing exactly affine function) or nonconvex, the shape function declared in Eq. (8) guarantees the unique solution inside a convex hull with a minimum distributed data set and always inhibits a weak Kronecker delta property. The weak Kronecker delta property means that essential (Dirichlet) boundary conditions can be treated the same way as in traditional FEM [43]. This is because the generalized meshfree convex method can achieve local convexity at the boundary nodes, which is a necessary condition for satisfying the Kronecker delta property [43]. Another method of constructing the meshfree shape functions, which satisfies the reproduction property of Eq. (9), is the moving least squares method with linear basis functions γ^T(x) = [1, x, y, z] [44]. This type of shape function, however, does not inhibit the Kronecker delta property, and special treatments, such as Lagrange multipliers, for enforcing essential boundary conditions must be used. This method, however, is not always preferred because the system of equations to solve becomes larger.

The updated unsmoothed nodal velocities are then inserted into Eq. (16), and the whole procedure is repeated until the desired simulated time is reached. The stabilization in this method is implicitly included in the smoothed velocity term of Eq. (16). If the velocity term was calculated with the standard method (i.e., unsmoothed velocity was used), then nodal integration would severely underestimate the contribution of the strain energy to the global stiffness matrix. In such cases, an extra stabilization term could be added to Eq. (19) at the expense of extra computational cost.

Fracture occurs for D = 1. For the SPG particles, when the damage variable reaches 1, the bond is broken by modifying the shape functions of the particle pair as shown in Eq. (27). Figure 2 shows the difference between a coupled and uncoupled damage model.

Finally, the general steps in the meshfree solution method can be shortly summarized as follows:

For each time-step:

1. Find the influencing nodes, construct the shape functions, and calculate its derivatives.

2. Compute the nodal mass from Eq. (15) and smooth the nodal velocity using Eq. (16).

3. Calculate strain using the smoothed nodal velocity, and calculate the stresses as well as the stiffness matrix, contact, and internal forces using Eqs. (20) and (21). Lastly, update the nodal density using Eq. (22).

4. Solve the discretized system of equations and obtain the accelerations from Eq. (19). Update the nodal velocity and displacement using Eqs. (23) and (24) as well as the unsmoothed velocity from Eq. (25).

5. Insert the unsmoothed velocity into Eq. (16) and repeat the whole procedure.

During the solution of the mechanical problem, the Von-Mises stresses are always checked and compared with the yield stress criterion. If the stress is in the elastic regime, Hook’s law of elasticity can be used. If the material has yielded, the stresses are scaled back accordingly. This is done using a special algorithm that uses the materials (experimentally measured) stress–strain data as input and calculates the correct stresses for a particular plastic strain, which is done in step 3. During steps 1–4, the damage variable D is calculated, and when it reaches D_crit, stresses are once again scaled using Eq. (31). For a given particle pair, when the plastic strain reaches the fracture strain e_crit, the shape functions for the particles are set to zero according to Eq. (27). This ensures that a discontinuity is created in the displacement field and thus simulates the fracture.

In the present work, we simulate a collision between two hemispherical asperities, which are vertically separated by the overlap denoted with δ. The simulation setup can be seen in Fig. 3, and it consists of two hemispheres sliding relative to each other with speed v_l and v_u.

Both hemispheres were discretized by SPG particles with an average nodal spacing of 0.5 µm, and each hemisphere consisted of about 90,000 particles. Nodal points in the eight-noded hexahedral elements were used in the original mesh to generate the meshfree discretization. The simulated time is 40 µs, which is enough time to ensure that both hemispheres pass each other. Initially, the asperities were deliberately positioned at a considerable distance from each other to ensure their dynamics would have maximum impact prior to the collision. It is worth mentioning that, in the context of simulating small micro-scaled asperity collisions, shock waves may not be the dominant phenomenon. At micro-scales, the impact energy levels are typically expected to be much lower compared to macro-scale impact events [21]. Therefore, the formation of shock waves may be less pronounced or negligible.

The material properties as well as all material models in the meshfree MC-SPG method are handled in the same way as in traditional FEM. The material properties such as density, Young’s modulus, and Poisson's ratio are used to calculate the mass matrix as well as the internal force vector of Eqs. (19) and (21). These are then used to calculate the stresses in the elastic regime. If, however, the loading crosses over to the plastic regime, the elastic stresses need to be scaled back at each new time iteration. This is done with a material model that uses known material stress–plastic strain relationships (input by the user) to correctly scale the elastic stresses.

For simplicity, all entities are given the same material, which is based on high-strength transformation-induced plasticity (TRIP) steel. The stress–strain curve is obtained from the study by Bleck et al. [46] and is shown in Fig. 4. The effective plastic strain refers to the true strain after elastic unloading, and the stress is simply the true effective stress.

In order to account for thermal effects, the stress–strain relationship for two different temperatures is used and input to the ls-dyna *MAT_106 material card. The material model *MAT_106 is an elastic-viscoplastic material model with Von-Mises yield criterion that can account for temperature-dependent material properties. Thermal expansion is also included in the model through the *MAT_ADD_EXPANSION card. All relevant material properties are presented in Table 1. In order to obtain reasonable simulation times, the density for the mechanical problem was increased by 1000 in order to decrease the dilatational wave speed c_d in the material and thus increase the critical time-step of Eq. (34). This artificial added mass does not affect the results due to the extremely small volumes that give almost no contribution to the total energy of the system, that is, kinetic energy ≪ internal strain energy [21].

Lastly, the input to the GISSMO damage model is based on fracture tests of high-strength TRIP steel obtained by Bai and Wierzbicki [45] and is summarized in Table 2.

The fracture data from Table 2 are used to calibrate the MMC fracture model, that is, by finding the coefficients C₁, C₂, and C₃ of Eq. (28) that minimized the residual of Eq. (30). Table 3 shows the values of the optimized parameters. An example 2D MMC curve is shown in Fig. 5. As can be seen, the fracture strain is not constant and strongly dependent on the stress state.

The damage exponent N and fade exponent M from Eq. (32) are important parameters as they describe the postnecking behavior of the material and determine whether the fracture is brittle or ductile. These parameters are usually material dependent and thus calibrated based on tensile fracture tests. The chosen values of these are presented in Table 4.

Finally, the general methodology is summarized as follows [21]:

• Calibrate the MMC model for fracture strains as a function of different stress states as well as damage parameters to be used as input to GISSMO.

• Generate background mesh for the spheres and use the nodal points for the meshfree particle discretization.

• Set up the contact problem in the multiphysics finite element software ls-dyna.

• Run the model and analyze the wear results.

This section presents the numerical results. Before showing the results for the three cases described in the previous section, a verification will be done with the Hertz solution for the linear elastic contact between the two spheres. All simulations were performed with the MPP R13.0.0 version of ls-dyna on 64 Xeon Gold 6248R 3 GHz cores computer cluster, which took roughly 5 h in real time to solve for each simulation. Depending on the situation, the words “failed,” “fracture,” and “material separation” will be used interchangeably as deemed best, but all of them are referring to the same phenomenon, that is, a discontinuity in the displacement field that occurs when the effective plastic strain reached the failure strain $εf$⁠, and the bond between two particles are broken according to Eq. (27).

In order to verify the accuracy of the present model, a comparison with the Hertz solution, with the same discretization size of 5 µm and 90,000 particles per sphere used for the three cases and material properties of Table 2, has been made for the perfectly elastic contact between two spheres. A normal load of 10 N has been applied on spheres with a radius of 0.1 mm, as shown in Fig. 6.

The comparison of the calculated pressure profile between the Hertz solution and the present work is shown in Fig. 7. As can be seen, there is a good agreement between the analytical and numerical solutions. The biggest deviation occurs at the edge of the contact where the pressure is very low. To handle that, a very fine mesh is required at the edge. However, such an approach would not yield a substantial enhancement in the precision of the pressure distribution within the contact area. The result shows, therefore, that the method can accurately predict the contact pressure for the perfectly linear elastic case and that the discretization size used in the present work is adequate.

A few snapshots from the simulations are shown in Figs. 8 and 9, along with the coefficient of friction (COF) for cases 2, 3, and 9 shown in Fig. 8. Figure 8 displays the Von-Mises stress through a color map, while Fig. 9 portrays the plastic strain using the same visual representation. For cases 2 and 3, there is mostly plastic deformation with minor number of detached particles. For case 9, the behavior is changed, and a transition to brittle fracture can be seen. As expected, both spheres deform and wear equally. This is due to the symmetrical nature of the problem as both spheres are identical and are moving toward each other with the same kinetic energy. The maximum stress for case 9 is highest around time 10 µs, which happens to be just before the fracture occurs. As the asperities start to fracture, the COF starts to decline. As the bonds between the particles fail, the stress is reset, and the failed particles act as wear debris. It is worth mentioning that some of the debris is still left on the surface of the spheres and hence may be mistaken for attached particles. However, since stresses are reset to zero when failure occurs, the stress fringe plot shown in Fig. 8 is used to find where the maximum damage occurs, and it is an effective tool to expose the detached particles and identify the fracture path. In cases 2 and 3, the smaller overlaps (4 and 6 µm) are mostly leading to plastic deformation and almost no material separation occurs. It is, however, evident that the amount of lost material (de-bonded particles) is larger in case 3 when the overlap is higher. For an overlap of 8 µm (case 9), a significantly larger portion of the asperities is being torn off due to the large stresses and strains, indicating that a transition to brittle fracture has occurred. To further understand the fracture behavior, the damage variable D for case 9 is fringe plotted and shown in Fig. 10. The area around the fracture zone is where the damage variable D has reached unity and is also where the fracture occurs.

A study of the stress states at the fracture zone revealed that the asperities are undergoing both shear and tension in those areas. The particles are being dragged in the sliding directions, and the areas where the damage variable reaches unity are also where shear and tension stresses are dominant. One of the main strengths of the present method lies in the fact that the fracture criteria are stress-state dependent, and thus, a more accurate modeling of the wear mechanism is obtained. This is because, for most metals, the fracture strains are not constant and depend on the triaxiality factor as well as the Lode parameter, which are both taken into consideration here.

The fracture behavior is observed to be a combination of both ductile and brittle in nature, but this behavior is purely based on the selection of the damage and fade exponent N and M. While studying the deformations, it was noted that the standard finite element method adopted from our previous work [17] has difficulties in modeling the large strains as well as correctly predicting the fracture path. This is precisely due to the difficulties in handling large deformations as well as element erosion (deletion) to simulate fracture (i.e., violates conservation of mass). The strength of the particle method is that it is not dependent on elements, which would otherwise become much distorted and cause negative Jacobian volumes in the stiffness matrix. The extreme element distortions would also require finer and finer time-steps to be taken, which would put further limitations on the method. These problems are avoided in the present solution, as the MC-SPG is a true meshfree method, and the discretization is not based on elements.

To better understand the effect of the stress-state dependency of the fracture criteria on the wear formation, simulations with constant fracture criteria were performed for all three cases. This implies that there is no stress-state dependency of the fracture criteria, and fracture occurs at the same stress level regardless of the type of deformation. Because the colliding asperities undergo a significant amount of shearing, the shear fracture strain criterion from Table 3 was used as a constant criterion for the GISSMO damage model. A comparison for case 9 is shown in Fig. 11, and the results clearly show the importance to account for the stress-state dependency to predict wear. The formation of wear particles occurs in different ways for the two cases. If the overlap had been lower, the wear formation would have been smaller and so also the stress-state dependency. As the wear increases due to increasing overlap values, the stress-state dependency becomes increasingly important. At higher overlap values, both the shear and tensile stresses become important, and the asperities are more prone to fracture as compared to compressive stresses. As shown in Fig. 10(c), when the deformations are low, in the first half of the collision, the COF behaves similarly for both cases. As the asperities are getting more and more deformed, thus reaching fracture, the COF starts to increase for the stress-state-dependent model. This behavior agrees well with results obtained by Ying and Hsu [47] for single asperity collision, thus verifying the importance of accounting GISSMO.

Next, the thermal response of the asperities is fringe plotted and shown in Fig. 12 along with the average transient temperature rise for each spherical asperity. As expected, the temperature development is the same for both upper and lower spheres in all simulated cases. As the problem is symmetrical, the same amount of plastic work is converted to heat for both spheres. The temperature is highest at the point of impact, which is also where the largest plastic deformation is observed. The temperature quickly reaches steady state as the fracture propagates, and the spheres cannot undergo more plastic deformation. For the third case, the temperature rise for the spheres is higher due to more plastic deformations as compared to the second case. The temperature development of the ninth case displays the same kind of behavior as seen for cases 2 and 3 except that the maximum temperature rise is much higher. A temperature study of this kind can help to explain if wear particles undergo a kind of quenching when they are rapidly cooled when they leave the wear zone. The wear particles will then become harder than the original material and potentially cause abrasive wear at a later stage.

During material failure, it is common for ductile materials to lose the ability to carry load and thus experience reduced stresses just before fracture, as shown in Fig. 2. This softening behavior is important because it determines the correct amount of plastic work converted to heat, and the GISSMO model in the present method accounts for this. For brittle materials, the fracture is sudden, and there is no softening, meaning that there is less plastic work and heat that is generated in the material. It is worth mentioning that the model can also take thermal effects into account for the fracture criteria. The experiments for fracture strains shown in Table 4 can be performed at different temperatures and then used as input to the GISSMO model.

Finally, the wear volume, coefficient of friction, and average normal and tangential forces are shown in Fig. 13. The COF tends to increase exponentially at around the critical overlap. The observed behavior can be attributed to the reduction in the required normal force caused by a significant increase in the overlap; meanwhile, the tangential force reaches a steady state. At around 7 µm, the wear mechanism changes from mostly plastic deformation to brittle fracture. This means that the asperities are no longer only deforming, but also a large portion is torn off due to the propagation of cracks inside the asperities, as shown in Fig. 10. This change in mechanism clearly shows a rapid increase in wear.

The average temperature rise is also shown for all cases in Fig. 13(c). The temperature shows a linear relationship with smaller amounts of overlap. As the transition from plastic deformation to brittle fracture is passed, the average temperature variation rate decreases, and the temperature reaches a constant value at the highest overlap. The temperature undergoes its most pronounced rate of change (steep slope) from 4 µm up to the critical overlap, which aligns with the occurrence of substantial plastic deformations. Once the critical overlap is exceeded, ductility diminishes, resulting in brittle fracture. Importantly, the brittle fracture phenomenon has minimal to no impact on heat generation, thus verifying the anticipated thermal behavior.

A meshfree method based on the momentum-consistent smoothed particle Galerkin has been adopted to numerically predict the sliding wear between two hemispherical asperities. The traditional FEM is known for having difficulties in handling large deformations as well as violating the conservation of mass when simulating fracture. With the method presented in this article, these issues are eliminated, and the fracture behavior can be simulated more correctly as compared to that by FEM. The presented model allows the user to better model and understand wear processes. From the presented results, the following conclusions can be drawn:

• Because the MC-SPG method does not rely on elements, it can more accurately simulate the large deformations and wear that occur in the asperities, as compared to FEM.

• Numerical issues such as negative Jacobian volumes and instabilities that occur with FEM are avoided with the presented method.

• The average temperature displays a linear behavior up until the critical overlap and converges to a constant value afterwards.

• For small overlap, the wear formation is small, and a constant fracture criterion may be sufficient to predict wear. However, as the overlap increases, it becomes increasingly important to account for the stress-state dependency when predicting wear.

• The average coefficient of friction increases exponentially with increasing overlap.

It must be pointed out that adhesive effects are disregarded in the present model. Adhesive effects will present additional complexity and will be treated in upcoming investigations.

• The Swedish Research Council (VR) (Grant No. 2020-03635).

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.