1 Introduction

Atomic-scale studies of the phenomenon of adhesive wear [1,2,3,4] suggest that the size of the smallest wear particles is of the order of a critical length scaleFootnote 1 of the material

$$d^{*} \sim \frac{\gamma E}{\sigma _{\text {m}}^{2}},$$
(1)

where E is the Young’s modulus, \(\gamma\) the surface energy, and \(\sigma _{\text {m}}\) the strength of the interfaces. While the smallest wear particles contribute to the formation of a third-body layer, TBL [5,6,7,8,9], there is no mention in the literature of a direct link between the material’s critical length scale and the TBL characteristics (thickness, rheology, etc.).

Experimental and numerical studies typically treat TBL and gouges macroscopically as (granular) fluids [10,11,12], whereas numerical studies at the atomic scale on comparatively very short durations clearly showcase the formation of distinct wear particles of various sizes and their aggregation [1, 13]. Investigating the transition between scales and behaviors is a challenging task, both for experimentalists and numericians. In experiments, the transition happens at durations that are short and challenging to observe [5]. Numerical simulations of wear, typically performed using molecular dynamics (MD), are more suited for assessing small and short events. However, the computational cost becomes too high for larger systems and longer times (e.g., Brink et al. [13]). Thankfully, a new formulation of inter-particular forces was developed to tackle this particular problem using the discrete element method, DEM [14], while staying close to the physics achievable with MD. This model is capable of reproducing the single-asperity wear behaviors captured by MD while bringing the computational cost down by 5 orders of magnitude. This significant gain is possible thanks to a coarse-graining procedure, by modeling deformable material with discrete particles 10 times larger than atoms.Footnote 2 This new method opens the door to studying TBL formation and evolution using adhesive wear physics.

While the aforementioned method was not yet used in a context involving more than the creation of a single wear particle, other numerical ways of investigating the evolution of third body were developed. For example, Zhang et al. [15] and Mollon [16] use a discrete element (DE) model with deformable grains to probe the role of discrete particle properties (e.g., stiffness and friction coefficient) on the overall physical behavior of the sliding interface. Mollon’s work [17] showcases the emergence of several types of microstructures, linked to different flow regimes. As it is commonly accepted with DEM simulations, the parameters given to the particles are not directly linked to material or loading parameters. The choice of the size of the particles is also not necessarily physically motivated.

The coarse-grained DE model of Pham-Ba and Molinari [14] allows performing similar simulations while directly linking the physics of the system to material properties. Systems studied using this model are shown to produce wear particles having a size independent of the size of the DEs (in a given range) and controlled instead by material properties and geometrical parameters. This feature removes the concern about the impact of the size of the DE (which has not always physical meaning) on the behavior of the whole system.

In this exploratory study, we perform numerical simulations of adhesive wear under sliding motion using the coarse-grained DE method [14]. Several material and loading parameters are evaluated, and the effects on the properties of the formed TBL are assessed, putting an emphasis on the role of the material’s critical length scale \(d^{*}\) on the rheology of the TBL.

2 Method

The chosen numerical method allows us to simulate a dynamical tribological system where detachment and reattachment of matter can occur, using spherical DE with cohesive forces. The inter-particular force parameters can be chosen to match several material parameters, notably the Young’s modulus, the strength, and the surface energy.

The force \({\varvec{F}}\) acting between a given pair of particles acts along the line going through the center of the particles (the normal direction, relative to particles’ surface). It is the sum of a stiffness or cohesive component \(F_{\text {N}}\) and a velocity damping force:

$${\varvec{F}} = -(F_{\text {N}} + c_{\text {N}}v_{\text {N}}) {\varvec{n}}_{\text {N}},$$
(2)

where \({\varvec{n}}_{\text {N}}\) is the unit vector pointing in the normal direction, \(v_{\text {N}}\) is the corresponding relative velocity, and \(c_{\text {N}}\) is a damping factor. The formulation used here is simplified compared to the full formulation established by Pham-Ba and Molinari [14]. The tangential forces between the particles are not modeled, reducing the number of model parameters, the trade-off being that the target Poisson’s ratio of the modeled material is fixed at \(\nu = 0.25\). Macroscopic apparent shear forces still exist since the normal inter-particular forces are sufficiently long-ranged. The expression of the normal force is

$$F_{\text {N}} = \left\{ \begin{array}{ll} k_{\text {N}} \delta _{\text {N}} &{} {\text {if }} \delta _{\text {N}} \leqslant \delta _{\text {e}}, \\ \displaystyle -\frac{k_{\text {N}} \delta _{{\text {e}}}}{\delta _{\text {f}} - \delta _{\text {e}}} (\delta _{\text {N}} - \delta _{\text {f}}) &{} {\text {if }} \delta _{\text {e}} < \delta _{\text {N}} \leqslant \delta _{\text {f}} , \\ 0 &{} {\text {otherwise,}} \end{array}\right.,$$
(3)

where \(\delta _{\text {N}}\) is the distance between the particles’ surfaces, \(\delta _{\text {e}}\) is the maximum elastic distance, \(\delta _{\text {f}}\) is the fracture distance, and \(k_{\text {N}}\) is the Hookean stiffness. The profile of the force is shown in Fig. 1.

Fig. 1
figure 1

Normal force between two particles as a function of inter-particular distance \(\delta _{\text {N}}\). There is interpenetration when \(\delta _{\text {N}} < 0\). The force has a cohesive part when \(\delta _{\text {N}} > 0\)

The stiffness, distance, and damping parameters are linked to the desired material properties through the following equations:

$$\begin{aligned}&k_{\text {N}} = \frac{A_{\text {eff}}E}{r_{i} + r_{j}} , \end{aligned}$$
(4)
$$\begin{aligned}&\delta _{\text {e}} = \frac{(r_{i} + r_{j}) \sigma _{\text {m}}}{E} , \end{aligned}$$
(5)
$$\begin{aligned}&\delta _{\text {f}} = \displaystyle \frac{4\gamma }{\sigma _{\text {m}}}, \end{aligned}$$
(6)
$$\begin{aligned}&c_{\text {N}} = \frac{2(1 - \eta )}{\pi } \sqrt{k_{\text {N}}m_{\text {eff}}} , \end{aligned}$$
(7)
$$\begin{aligned}&r_{\text {eff}} = \min (r_{i}, r_{j}), \end{aligned}$$
(8)
$$\begin{aligned}&A_{\text {eff}} = 2\sqrt{2} \, r_{\text {eff}}^{2} , \end{aligned}$$
(9)
$$\begin{aligned}&m_{\text {eff}} = \frac{m_{i} m_{j}}{m_{i} + m_{j}} , \end{aligned}$$
(10)

where \(r_{i}\) and \(r_{j}\) are the radii of the interacting particles, \(m_{i}\) and \(m_{j}\) their respective masses, E is the target Young’s modulus, \(\sigma _{\text {m}}\) the target strength, \(\gamma\) the target surface energy, and \(\eta\) the target restitution coefficient. These equations are valid as long as the radii of the particles satisfy

$$r_{i} + r_{j} \leqslant d_{\text {m}} ,$$
(11)

where

$$d_{\text {m}} = 4 \frac{\gamma E}{\sigma _{\text {m}}^{2}}$$
(12)

can be seen as a maximum discrete particle diameter [specific to this formulation but still closely related to the material’s critical length scale \(d^{*}\) (1)]. If (11) is not satisfied, the actual values of \(\sigma _{\text {m}}\) and \(\gamma\) (of the DEM system) may deviate from their imposed targets.

The Young’s modulus E, the density \(\rho\), and the size of the system L are all arbitrarily set to 1 (a coherent system of units can also be chosen arbitrarily). Their actual value is not important, as all the other dimensional quantities can be expressed relative to those three. The restitution coefficient between particles (dimensionless) is set to \(\eta = 0.95\).Footnote 3 Thus, the only remaining free parameters needed to define the material are the strength \(\sigma _{\text {m}}\) and the surface energy \(\gamma\).

Since our goal is to relate the material parameters to the properties of the TBL, such as its thickness, it is sensible to combine the material parameters into a quantity that has units of length. We can use the critical length of the material (1), which was found numericallyFootnote 4 in our situation to be roughly

$$d^{*} = 32 \frac{\gamma E}{\sigma _{\text {m}}^{2}},$$
(13)

following the scaling of (1), with the interface strength taken as the material bulk tensile strength \(\sigma _{\text {m}}\). The last parameter of interest is the normal pressure \(p_{\text {N}}\) applied on the system. Our simulations differ from those in Mollon [17], where the gap between the surfaces is kept constant instead of imposing a normal load.

Fig. 2
figure 2

Schematic of the simulated DEM systems. On the right part of the figure, some relevant sizes are represented. The whole system is filled with particles whose mean diameter is \(d_{0} = 0.01\ L\). The thickness of the system (quasi two-dimensional) is \(B = 3\ d_{0}\). The material properties are chosen such that the resulting \(d^{*}\) is equal to either \(0.05\ L\), \(0.1\ L\), or \(0.2\ L\), which are represented on the right. The thin darker areas at the top and bottom of the system are modeled as rigid bodies. The particles making them cannot move relative to each other. The bottom rigid body is fixed, while the top one receive a constant compressive normal pressure \(p_{\text {N}}\) and is dragged horizontally at a constant velocity v. At the initial state of the simulations, the top and bottom part of the system are disjoined, i.e., there is no particle that is part of both the orange body and the blue body. The area enclosed by the dashed line is the region represented in the subsequent figures (Color figure online)

All studied systems fit in a quasi-2D square of side length L (see Fig. 2). The thickness B of the system is adjusted to the average size \(d_{0} = 0.01 \, L\) of the particles to be \(B = 3 \, d_{0}\). Using such a system instead of a fully two-dimensional one allows us to use spherical particles with sizes distributed around a mean value to create a disordered system without weak planes, while ensuring a tight packing of the particles. The sizes of the particles follow a log-normal distributionFootnote 5 with a mode (most frequent value) of \(d_{0}\) and a standard deviation of \(0.1 \, d_{0}\). The distribution is truncated between \(0.75 \, d_{0}\) and \(1.25 \, d_{0}\). It is shown in Pham-Ba and Molinari [14] that \(d_{0}\) and the whole distribution can be varied without significantly changing the properties of the resulting medium or its adhesive wear behavior, as long as \(d_{0}\) is smaller than \(d^{*}\) (the minimum wear particle size for the target material). An amorphous arrangement of particles is generated by first inserting the particles at random positions into the system, then performing a critically damped dynamical simulation, which makes the particles rearrange into a configuration that minimizes overlaps. The boundaries of the system are allowed to move such that the average internal stresses drop to 0. More details regarding the exact parameters of the relaxation procedure are given by Pham-Ba and Molinari [14].

The simulated systems consist of two halves filling the space of size \(L \times L \times B\) and having a non-matching interface (in the sense of the discrete particles’ arrangement). The bottom half is connected to a fixed rigid body (made of non-moving particles) of height \(1.5 \, d_{0}\), and the top half is dragged by a rigid body of the same size at a fixed velocity of \(0.1 \sqrt{E / \rho }\) over a total distance of \(50 \, L\). The other directions have periodic boundary conditions. A normal pressure \(p_{\text {N}}\) is applied on the top rigid body.

The normal pressure is varied between \(0.01 \, E\), \(0.02 \, E\), and \(0.04 \, E\). For the other two parameters, one could boldly expect the thickness of the TBL to be only dependent on the value of \(d^{*}\) (linked to the minimum size of wear particles), but not on the individual values of \(\gamma\) and \(\sigma _{\text {m}}\). So, at first, to assess the direct influence of \(d^{*}\),  the strength is set to \(\sigma _{\text {m}} = 0.1 \, E\) and the surface energy \(\gamma\) is adjusted using (13) to reach \(d^{*}\) values of \(0.05 \, L\), \(0.1 \, L\), and \(0.2 \, L\). Then, in a second stage, to assess the individual effects of \(\sigma _{\text {m}}\) and \(\gamma\), the critical size is fixed to \(d^{*} = 0.1 \, L\) and the strength is varied between \(0.1 \, E\), \(0.14 \, E\), and \(0.2 \, E\), while the surface energy \(\gamma\) is still adjusted (increased when \(\sigma _{\text {m}}\) increases) to match \(d^{*}\). In summary, for each normal pressure \(p_{\text {N}}\), the tested sets of \(\sigma _{\text {m}}\) and \(\gamma\) are shown in Table 1. To get a sense of how material parameters affect the interaction between particles, we can also compute the maximum interaction distance \(\delta _{\text {f}}\) given by (6), which is related to ductility.

Table 1 Sets of tested parameters (apart from \(p_{\text {N}}\)) and particles’ maximum interaction distance

3 Results and Discussion

All the performed simulations can be qualitatively categorized into one of the three aspects shown in Fig. 3 and in Supplemental Videos. In the first case (a), wear particles of various sizes are formed, continuously agglomerated together, and reattached to the sliding bodies. The size of the gap remains visually stationary. In the second aspect (b), comparable to a fluid flow, the two bodies are fully connected and some mixing (vertical diffusion) occurs over time. In some instances, the simulations transition into another state (c), featuring large wear particles, growing in size, and reaching up to the full size of the system. In any case, there is always a first stationary phase, (a) or (b), lasting for an arbitrary amount of time. In the quantitative analysis of the simulations, reported thereafter, we dropped the ending parts involving large growing wear particles (c). It is unknown whether all systems would eventually reach this third regime, as our simulations are currently relatively limited in time (sliding distance of \(50 \, L\)). Also, the limited size of the system hinders the analysis, as the TBL can travel vertically and reach the system’s boundary, triggering a transition toward an unphysical regime, not represented here. Therefore, there is still room for investigation about the mechanisms triggering the third regime and for inspecting the structure of the resulting TBL. We suspect that there might be specific geometric and loading requirements for the transition to happen, that are randomly met in some simulations. For further studies, one way of simulating taller systems over longer periods of time without absurd computational costs could be to use coupling methods between finite elements and DEM [18].

Fig. 3
figure 3

Examples of qualitative regimes. The images are cropped according to the region shown in Fig. 2. Videos for each case are available as Supplementary Material. a Example of mixed regime, featuring the formation of wear particles, their aggregation, and reattachment to the sliding bodies. The gap remains constant in size. The colors represent the initial vertical position of the particles. The parameters of the simulation are \(p_{\text {N}} = 0.01 \, E\), \(d^{*} = 0.1 \, L\), and \(\sigma _{\text {m}} = 0.1 \, E\). b Example of shear band. There is no visible gap between the two bodies, and the particles are migrating vertically. The parameters of the simulation are \(p_{\text {N}} = 0.01 \, E\) and \(d^{*} = 0.1 \, L\), as in a, and \(\sigma _{\text {m}} = 0.2 \, E\). c Example of formation of large wear particles, that grow in size until reaching the size of the system. The parameters of the simulation are \(p_{\text {N}} = 0.01 \, E\) and \(d^{*} = 0.1 \, L\), as in a and b, and \(\sigma _{\text {m}} = 0.14 \, E\) (Color figure online)

The emergence of a stationary state featuring wear particles and reattachment (Fig. 3a) is worth noting, as it is not something that could be observed using MD simulations limited to a smaller scale, or larger scale DEM simulations missing the details of wear particle formation.

The geometric properties of the TBL can also be measured quantitatively. The boundaries of the TBL (and thus its thickness) are computed by looking at the horizontal velocity profile along the vertical axis (see Fig. 4a). This measurement is not very accurate because the velocity profile is oscillating, meaning that arbitrary thresholds must be used to determine the boundaries of the TBL. Still, we can confirm the stationary nature of the main simulation phases by using the evolution of the TBL thickness over time (see Fig. 5). We notice that in the case of the shear banding (Fig. 3b), even if vertical mixing of the particles is observed over the whole height of the body, the TBL thickness (related to the vertically confined shear strain) keeps a finite size. A similar mixing behavior was observed in MD simulations [19].

Fig. 4
figure 4

Examples of velocity and density profiles. The vertical coordinate is normalized by the current height H of the system (initially equal to L). a Example of averaged horizontal velocity profile as a function of the vertical position. The vertical axis is cut into 400 bins and the average velocity of each section is computed. The particles have an average diameter of \(d_{0} = 0.01 \, L\). The boundaries of the third-body layer (TBL, shaded area) are determined by the velocity profile and velocity thresholds. The bottom boundary of the system is fixed, and the top has an imposed velocity of \(v_{0}\). Here, the thresholds are \(v_{\text {th}} = 0.1 \, v_{0}\) and \(v'_{\text {th}} = 0.9 \, v_{0}\). The TBL thickness \(h_{\text {TBL}}\) can be deduced. b Example of estimated density profile. As for the velocity profile, the density is computed by binning on 400 sectors (‘filtered’ curve). The density profile is filtered using a Savitzky–Golay filter of order 3 on 7 samples. A Gaussian curve (chosen arbitrarily) is also least-square fitted to the unfiltered data, leading to the minimum density inside the TBL \(\rho _{\text {TBL}}\). TBL boundaries estimated from the velocity profile are shown for comparison (gray region). It is possible to have a homogeneous density with \(\rho _{\text {TBL}} = \rho _{0}\), while still having well-defined TBL boundaries (given by the velocity profile)

Fig. 5
figure 5

Examples of TBL thickness evolution with \(p_{\text {N}} = 0.04\ E\) and \(d^{*} = 0.1\ L\). Raw measurements are shown as transparent lines, whereas the solid lines are moving averages. At the very beginning of sliding, before any amount of damage, the shear deformation is uniform, which is seen here as a peak starting from almost \(h_{\text {TBL}} / H = 1\), H being the current height of the system. A steady regime is quickly reached, being either a ‘mixed’ regime or a shear band. After some time, some simulations transition into the regime of large wear particles formation, identifiable here by erratic changes of TBL thickness. In other cases, the TBL reaches a boundary of the system (which is not physical), and the TBL thickness drops

The average density of particles inside the TBL can also be measured (see Fig. 4b). This measurement is more reliable and accurate than the thickness measurement. When there is a visible gap between the sliding surfaces (populated by wear particles, see Fig. 3a), the TBL density is lower than the density of the bodies.

The dependence of the TBL density \(\rho _{\text {TBL}}\) on \(d^{*}\) and the normal pressure is shown in Fig. 6a (the strength is kept constant at \(\sigma _{\text {m}} = 0.1 \, E\)). Surprisingly, \(d^{*}\) has no significant influence on \(\rho _{\text {TBL}}\), compared to the effect of \(p_{\text {N}}\). Imposing a larger normal load on the bodies has the expected effect of increasing the TBL density, which can be interpreted as increased wear volume production. The positive correlation between load and wear rate is compatible with the experimentally supported Archard’s wear model [20]. With these parameters (strength fixed), neither \(d^{*}\) nor \(p_{\text {N}}\) has a significant effect on the TBL thickness (not shown here).

Fig. 6
figure 6

Dependence of TBL minimum density on material properties. a Dependence on pressure and material’s characteristic size (critical length scale). The material strength is kept constant at \(\sigma _{\text {m}} = 0.1 \, E\) and the surface energy \(\gamma\) is accommodated to obtain different values of \(d^{*}\). The density depends mostly on the normal pressure \(p_{\text {N}}\). b Dependence on pressure and material strength at a constant characteristic size. The surface energy \(\gamma\) is increased conjointly with the strength \(\sigma _{\text {m}}\) in order to keep the characteristic size constant at \(d^{*} = 0.1 \, L\). Even if \(d^{*}\) is constant, the TBL density is strongly dependent on the other parameters

Figure 6b shows the individual effects of the material parameters, where \(d^{*}\) is kept constant. There is a clear effect of increasing the strength on the TBL density. A larger strength (and surface energy) keeps the TBL density high, meaning that the gap between the bodies remains more easily closed. The maximum density is a bit over the initial density of the material due to the normal pressure applied on the system. Increasing the strength and the surface energy (at constant \(d^{*}\)) also has the overall effect of increasing the TBL thickness, intuitively comparable to a more viscous flow.

Regarding friction, at the scale we are working, tangential forces are dominated by adhesive forces. Therefore, rather than looking at the friction coefficient \(p_{\text {T}} / p_{\text {N}}\), we look directly at the average value of \(p_{\text {T}}\) in the different steady-state regimes. \(p_{\text {T}}\) is computed from the total horizontal force acting on the top or bottom boundary. Being mainly controlled by adhesive forces, \(p_{\text {T}}\) is mostly dependent on the material’s (tensile) strength \(\sigma _{\text {m}}\) and on the contact area between the top and bottom bodies. When there is full contact (shear band regime shown in Fig. 3b, maximum TBL density), \(p_{\text {T}}\) is equal to the shear strength of the material, which is measured to be roughly equal to \(\tau _{\text {m}} = 0.22 \, \sigma _{\text {m}}\). In the mixed regime (Fig. 3a), the contact area is lower, and it depends on the normal pressure \(p_{\text {N}}\). In fact, the contact area can be estimated to be proportional to the TBL density, and the measured \(p_{\text {T}}\) is indeed near \(\tau _{\text {m}} \rho _{\text {TBL}} / \rho _{0}\). As such, \(p_{\text {T}}\) is not directly dependent on \(d^{*}\) (like \(\rho _{\text {TBL}}\)). Finally, in the regime with large wear particles (Fig. 3c), the contact area between the sliding bodies and the rolling particles is very small (not correlated with the TBL density anymore). This leads to a lubricated behavior, with tangential loads as low as \(0.001 \, E\). This lubricative wear behavior is also observed experimentally [21]. Understanding how it arises and how it is controlled could be very beneficial for the tribology community.

4 Conclusion

In this short exploratory numerical study, we demonstrated the capabilities of the coarse-grained DE model of Pham-Ba and Molinari [14] in representing adhesive wear situations. We showcased the emergence of a mixed regime involving wear particle creation and their aggregation into a TBL. The material’s critical length scale, known to play a major role in the sizing of the first formed wear particles, was shown to have a limited effect on the properties of the TBL. Instead, the strength and the surface energy of the material must be considered directly when treating third body macroscopically. This works opens the door to the exploration of more ambitious setups involving adhesive wear and the formation and evolution of a TBL. The next steps could involve larger systems, fully three-dimensional systems, and other loading conditions such as gap-controlled sliding.

5 Supplementary Material

Videos for the simulations shown in Fig. 3 are available along the online version of this article.