# Crack nucleation in the adhesive wear of an elastic-plastic half-space

Lucas Frérot, Guillaume Anciaux, Jean-François Molinari

► To cite this version:

Lucas Frérot, Guillaume Anciaux, Jean-François Molinari. Crack nucleation in the adhesive wear of an elastic-plastic half-space. Journal of the Mechanics and Physics of Solids, 2020, 145, pp.104100. 10.1016/j.jmps.2020.104100. hal-04760509

HAL Id: hal-04760509

https://hal.science/hal-04760509v1

Submitted on 30 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Crack Nucleation in the Adhesive Wear of an Elastic-Plastic Half-Space

Lucas Frérot<sup>a,*</sup>, Guillaume Anciaux<sup>a</sup>, Jean-François Molinari<sup>a</sup>

<sup>a</sup>Institute of Civil Engineering, cole Polytechnique Fdrale de Lausanne (EPFL), CH 1015 Lausanne, Switzerland

## Abstract

The detachment of material in an adhesive wear process is driven by a fracture mechanism which is controlled by a critical length-scale. Previous efforts in multi-asperity wear modeling have applied this microscopic process to rough elastic contact. However, experimental data shows that the assumption of purely elastic deformation at rough contact interfaces is unrealistic, and that asperities in contact must deform plastically to accommodate the large contact stresses. We therefore investigate the consequences of plastic deformation on the macro-scale wear response using novel elastoplastic contact simulations. The crack nucleation process at a rough contact interface is analyzed in a comparative study with a classical $J_2$ plasticity approach and a saturation plasticity model. We show that plastic residual deformations in the $J_2$ model heighten the surface tensile stresses, leading to a higher crack nucleation likelihood for contacts. This effect is shown to be stronger when the material is more ductile. We also show that elastic interactions between contacts can increase the likelihood of individual contacts nucleating cracks, irrespective of the contact constitutive model. This is supported by a statistical approach we develop based on a Greenwood-Williamson model modified to take into account the elastic interactions between contacts and the shear strength of the contact junction.

**Keywords:** wear, crack nucleation, plasticity, saturation, roughness, contact

Adhesive wear, as defined by Burwell (1957), happens when two asperities, belonging to different surfaces being rubbed together, form a junction (which can be adhesive, chemical, mechanical, etc.) that supports the shear caused by the relative sliding of the surface, such that cracks propagate in the bulk of the materials, causing the detachment of a third-body. In this setting, the removal of a wear debris particle from the surfaces is mainly driven by a fracture process, and as such obeys a balance between the energy release rate (i.e. the energy released by the crack front advancing) and the fracture toughness (i.e. the energy required to create new surfaces). This Griffith (1921) energy balance has been verified in atomistic simulations (Aghababaei et al., 2018, 2016; Brink and Molinari, 2019) and reduced to a critical length-scale $d^* \propto G\Delta w/\tau_j^2$, with G being the shear modulus, $\Delta w$ the fracture energy and $\tau_j$ the junction shear strength. This gives a simple geometric criterion for the formation of hemispherical wear particles: if the contact diameter between two hemispherical asperities is larger than $d^*$ then a wear particle detaches from the surface upon shearing of the system. The issue of transposing asperity-scale wear mechanisms to multi-asperity contact is key in the goal of formulating predictive wear models (Meng and Ludema, 1995; Vakis et al., 2018). Popov and Pohrt (2018) and Pham-Ba et al. (2019) have recently proposed energy-based models for the formation of wear particles in multi-asperity settings. The former investigates the formation of hemispherical wear particles in an elastic rough surface contact by computing an energy-favored particle diameter based on the elastic deformation energy of the contact solution. The latter formulates the energetic competition between the formation of a single vs. multiple wear particles (for 2D line contacts), thus giving an energy approach to the crack shielding mechanism that leads to *disjoint but sufficiently close* contacts forming a single wear particle (Aghababaei et al., 2018).

In a previous work (Frérot et al., 2018), we have applied the critical length-scale concept to rough elastic contact by defining a critical cluster area $A^* \propto (d^*)^2$ above which micro-contacts should form a wear particle. This work however suffers from two model inadequacies: (a) the contact solution is given by

<footnote>*Corresponding Author
Email address: lfrerot1 <at> jhu <dot> edu (Lucas Frérot)</footnote>
an elastic contact model, (b) it assumes that $A^*$ exists and is proportional to the square of $d^*$. The latter is related to the topography and shape of contacts. Contacts resulting from interfaces with rough surfaces are not disk-shaped and the crack is not expected to produce hemispherical wear particles. Moreover, this does not account for disjoint contacts that may form a single particle (Aghababaei et al., 2018; Pham-Ba et al., 2019). Thence, it is unclear if the Griffith balance can be characterized with a comparison as “naive” as $A \stackrel{?}{\ge} A^*$ with A being the area of a single contact cluster.

The former shortcoming (a) provides to the wear models developed in (Frérot et al., 2018) an unrealistic contact solution. Since both $A^*$ and the contact solution indirectly depend on $\sigma_y^1$, the outcome of an elastic contact problem (which is independent of $\sigma_y$) leads to a paradox: more ductile materials (with lower $\sigma_y$) have a higher $A^*$ and thus wear less than more brittle materials (with higher $\sigma_y$). This is due to the contact solution being insensitive to changes in $\sigma_y$, but also to the lack of surface roughness evolution in sliding. When the contact of two asperities does not create a wear particle, plastic smoothing of the asperities occurs, thus creating larger contacts. In any case, a contact model incorporating plastic effects is needed.

Experimental data clearly shows that some form of plasticity must occur at rough contact interfaces (Bowden and Tabor, 1939; Greenwood and Williamson, 1966; Weber et al., 2018; Zhang et al., 2019). Modeling these interfaces with a non-linear constitutive behavior is however a challenge because of the multi-scale nature of rough surfaces. Pei et al. (2005) were the first to use the finite-element method to study elastic-plastic rough contact with a classical von Mises formulation (Simo and Hughes, 1998). Jacq et al. (2002) have developed a volume integral method that we have refined with a Fourier approach to be able to handle the large discretization requirements of multi-scale rough surfaces (Frérot et al., 2019b). The majority of published works on elastic-plastic contact does not rely on classical formulations of plastic flow, but rather on the concept of surface flow pressure, which is associated to the hardness of a material (Archard, 1953; Greenwood and Williamson, 1966; Majumdar and Bhushan, 1991; Persson, 2001a). The surface flow pressure is usually taken as the maximum value of the mean pressure caused by an indenter of a given shape (it is therefore shape-dependent). Tabor (1951) has shown that for a spherical indenter, the mean pressure saturates at a value close to $3\sigma_y$ (with $\sigma_y$ being the yield stress). The models previously mentioned are thereafter referred to as “saturation models”, in the sense that they apply this concept of a maximum average pressure to a multi-asperity contact model (see e.g Tabor, 1951, chap. 9) and assume that a given contact cannot have a pressure exceeding the saturation pressure noted $p_m^2$. They have been used in conjunction with boundary integral approaches (Almqvist et al., 2007) to study friction (Weber et al., 2018), but have to our knowledge never been compared to classical plasticity formulations, and the relevance of the choice between the two plasticity models has never been studied.

Akchurin et al. (2016) and Li et al. (2019) have used a saturation plasticity model to compute the contact solution and applied a stress based criterion for the removal of debris particles: from the contact pressure profile, they computed the resulting von Mises stress caused in a *purely elastic* medium. Then the zones of the material where the von Mises stress exceeds the yield stress are removed, changing the surface profile. This has the advantage of foregoing any geometrical consideration, at the expense of providing an ad-hoc removal process that is not derived from the fracture energy balance, as well as using a stress distribution that does not account for plastic deformations.

In this work, we wish to investigate the multi-asperity wear process from a fracture-mechanics perspective and understand the influence of plasticity in the contact model on the global wear response. To this end, we focus on the crack nucleation process in the contact of a rigid self-affine rough surface with an elastic-plastic flat half-space (without surface energy effects). The contact is resolved using full-order simulations (Frérot et al., 2019b). One measure of particular importance is the crack spatial density. While it is not a measure of wear itself, crack nucleation is a necessary process of wear, and understanding what are the roles of the normal load, the critical nucleation stress, the junction resistance, and plastic behavior in crack nucleation is a fundamental step towards predictive wear models. We first highlight the importance of the choice of a plasticity model and the implications it may have on the contact response (section 1). We then show how

<footnote n="1">¹We have $A^* \propto \sigma_y^{-4}$ if one assumes $\tau_j \propto \sigma_y$ and the total contact area $A_c \propto \sigma_y^{-1}$ according to saturation models</footnote>

<footnote n="2">²More often than not, the saturation pressure is referred to as “hardness”. As Burwell and Strang (1952) discuss, the saturation (or flow) pressure cannot be absolutely known but is of the same order of magnitude as the value given by usual hardness tests. We therefore keep separate notations for clarity.</footnote>
the crack nucleation density in a rough surface elastic-plastic contact depends on the fracture mechanics properties of the material, as well as the applied load and the junction shear strength (section 2). To rationalize the differences between the elastic, the saturation and the von Mises plasticity approaches, we study the contact behavior of a single asperity to understand under which conditions a crack can nucleate and what is the influence of residual plastic deformations on this process (section 3). These findings are applied to a simple multi-asperity contact model (Greenwood and Williamson, 1966) in order to obtain qualitative analytical predictions of the scaling of the crack density with respect to system parameters like the applied normal load (section 4). These predictions are confronted to the elastic and elastoplastic simulations results which are able to reproduce the contact shielding effect under shear loading, as seen in molecular dynamics simulations (Aghababaei et al., 2018; Pham-Ba et al., 2019). The elastic-plastic results show that ductile materials in contact with rough surfaces produce more crack nucleation sites than brittle materials due to the residual stresses caused by plastic deformations. This effect is not captured by the elastic contact model nor the saturation plasticity model, indicating that the resolution of the aforementioned wear paradox should include the full plastic contact response. This further implies that the true contact area is not the only key quantity in wear modeling.

# 1. Elastic-plastic contact

At our disposal are (at least) two formulations of the elastic-plastic contact of solids, the choice of which may have an impact on the subsequent results we wish to obtain. The first formulation, which has been used in the finite-element studies of Pei et al. (2005), follows the classical modeling hypothesis of metal plasticity (Simo and Hughes, 1998), which have both experimental (Bui, 1969) and theoretical backgrounds (Reddy and Martin, 1994), and additionally are valid in other context than contact. The second, developed by Bowden and Tabor (1939) and extended by Almqvist et al. (2007) in conjunction with a boundary integral approach, postulates that the surface contact pressure should nowhere exceed a maximum value $p_m$. This is based on observations that for spherical indentation the mean contact pressure does not exceed a value around $3\sigma_y$ (Tabor, 1951). Recent finite-element simulations (Ghaednia et al., 2017; Krithivasan and Jackson, 2007; Song and Komvopoulos, 2013) show that $p_m/\sigma_y$ may depend on the ratio $\sigma_y/E^*$ (with $E^* := E/(1 - \nu^2)$ being the contact modulus) as well as the wavenumber in the case of sinusoidal contact surfaces. Despite these reports, saturation models are often used in computational tribology (Akchurin et al., 2016; Li et al., 2019; Weber et al., 2018) due to their simplicity and ease of implementation. Besides increasing the magnitude of the true contact area compared to elastic contact, plasticity influences other aspects of the contact interface (such as contact pressures). These additional aspects may be key ingredients in wear modeling. For this reason we wish to provide a comprehensive comparison between the von Mises associated plasticity and the saturation plasticity with $p_m = 3\sigma_y$ in a rough contact situation, and determine the consequences of the choice of one model over the other. We start by giving the full mechanical formulation for both models, then proceed to the comparison.

*Definitions.* In this work, we consider a deformable three-dimensional solid $\mathcal{B}$ spanning a half-space, with its (flat) boundary noted $\partial\mathcal{B}$. Moreover, we suppose a horizontal periodicity in the cell $\mathcal{B}_p = [0, L]^2 \times \mathbb{R}^+$. We note $\sigma$ the Cauchy stress tensor, which is related to the small-strain tensor $\varepsilon$ and the plastic strain tensor $\varepsilon_p$ by the relation $\sigma = C : (\varepsilon - \varepsilon_p)$ where $C$ is the usual isotropic linear elasticity tensor. The strain tensor is given by kinematic compatibility as a function of the displacement field: $\varepsilon = \nabla^{\text{sym}}u$. Finally, $\sigma$ is expected to be divergence-free to satisfy conservation of momentum without volume forces.

We additionally define some surface quantities: $t$ and $p := t \cdot e_3$ are respectively tractions and normal pressures applied on $\partial\mathcal{B}$. Other surface quantities are noted with an over-bar $\bar{\cdot}$ when not explicitly defined on $\partial\mathcal{B}$, e.g. $\bar{u}$ is the surface displacement.

*Saturation: perfect plasticity*

The simplest form of saturation model, conceptually close to the notion of “perfect plasticity”, is given as (Almqvist et al., 2007)

$$ \min_{p} \left\{ \frac{1}{2} \int_{\partial \mathcal{B}_p} p \overline{\mathcal{M}}[p] \, dS - \int_{\partial \mathcal{B}_p} ph \, dS \right\}, \qquad (1a) $$
which is a problem of finding the surface pressures p minimizing the complementary energy of the system under the constraints

$$p \ge 0, \qquad (1b)$$

$$p \le p_m, \qquad (1c)$$

$$\int_{\partial \mathcal{B}_p} p \, dS = W. \qquad (1d)$$

The linear operator $\overline{\mathcal{M}}$ gives the normal surface displacement due to the applied pressure p if $\mathcal{B}$ is assumed elastic; h is a continuous function representing the rough surface brought in contact with $\partial \mathcal{B}$ and W is the total applied normal load in the periodic cell boundary $\partial \mathcal{B}_p$. The gap is defined as $g := \overline{\mathcal{M}}[p] - h$ and should satisfy weak Hertz-Signorini-Moreau conditions (Weber et al., 2018):

$$g \ge 0 \quad \text{where } p < p_m, \qquad (2a)$$

$$p \ge 0, \qquad (2b)$$

$$p g = 0 \quad \text{where } p < p_m. \qquad (2c)$$

The solution to the above constrained optimization problem yields a negative gap where $p = p_m$. The magnitude of the negative gap is often assumed to be the magnitude of the residual plastic displacements. Since the weak optimality conditions do not represent a physical system (the gap should be non-negative everywhere to avoid body interpenetration), it is necessary to replace h in eq. (1a) by $h_{\text{mod}}^s := h + h_{\text{pl}}$, with $h_{\text{pl}}$ being in principle a correction due to residual plastic displacements, therefore $h_{\text{pl}} := -(\overline{\mathcal{M}}[p] - h)$ where $p = p_m$. Weber et al. (2018) propose an iterative scheme to solve for $h_{\text{pl}}$ which we have implemented and made available in the open-source contact library Tamaas (Frérot et al., 2019a) (https://c4science.ch/tag/tamaas/).

The “perfect plasticity” aspect of the model comes from the fact that $p_m$ is homogeneous on $\partial \mathcal{B}$ and constant. Weber et al. (2018) have amended this hypothesis to include a form of hardening. The saturation pressure is simply expressed as a linear function of $h_{\text{pl}}$ (i.e. the initial saturation stress is zero, and rises in proportion with $h_{\text{pl}}$). We will however not discuss this particular model here.

### $J_2$ von Mises plasticity

For the Cauchy stress tensor $\sigma$, the von Mises yield function $f_y$ is defined as

$$f_y(\boldsymbol{\sigma}) = \sqrt{\frac{3}{2}} \|\boldsymbol{s}\|, \quad \text{where } \boldsymbol{s} := \boldsymbol{\sigma} - \frac{1}{3} \text{Tr}(\boldsymbol{\sigma}) \mathbf{I}. \qquad (3)$$

The equivalent cumulated plastic strain is expressed as the integral of the plastic strain rate $\dot{\boldsymbol{\epsilon}}^p$ from some reference time $t_0$:

$$e^p := \sqrt{\frac{2}{3}} \int_{t_0}^{t} \|\dot{\boldsymbol{\epsilon}}^p\| \, dt. \qquad (4)$$

The admissibility and consistency conditions are written as:

$$f_y(\boldsymbol{\sigma}) - f_h(e^p) \le 0, \qquad (5a)$$

$$(f_y(\boldsymbol{\sigma}) - f_h(e^p))\dot{\boldsymbol{\epsilon}}^p = 0, \qquad (5b)$$

where $f_h$ is the hardening function. In this work, we will only consider functions of the form $f_h(e^p) = \sigma_y + E_h e^p$, with $\sigma_y$ the initial yield stress and $E_h$ the hardening modulus$^3$. The associated flow rule that determines $\dot{\boldsymbol{\epsilon}}^p$ is given by (Johnson, 1985):

$$\dot{\boldsymbol{\epsilon}}^p = \frac{3\dot{e}^p}{2f_y(\boldsymbol{\sigma})} s(\boldsymbol{\sigma}). \qquad (6)$$

The numerical integration of the relations above is typically done with a backwards Euler scheme and is classical to the solid mechanics literature (Simo and Hughes, 1998). Its coupling with the equilibrium and contact conditions is however non-trivial.

<footnote>$^3$This corresponds to linear isotropic hardening.</footnote>
**Solution strategy.** Jacq et al. (2002) established a numerical method for the solution of the elastic-plastic rough contact problem, which we summarize here. The method consists in solving the contact and the plasticity problems separately. The contact problem is solved for fixed plastic deformations: it is effectively an elastic contact problem with a rough surface $h_{\text{mod}} := h - \bar{u}_3^p$, with $\bar{u}_3^p$ the vertical component of the actual<fnref n="4" /> plastic residual displacement. Various solution strategies for the elastic rough contact problem are available in the literature (Bemporad and Paggi, 2015), and we use here the modified conjugate gradient algorithm of Polonsky and Keer (1999) coupled with the spectral approach of Stanley and Kato (1997) for the gradient computation involving the operator $\overline{M}$.

The plastic problem is solved with fixed boundary tractions, meaning that the contact area does not evolve during the resolution of the plastic strain increment. The procedure we employ, fully detailed in (Frérot et al., 2019b), relies on an implicit incremental volume integral equation formulation proposed by Telles and Carrer (1991). The total strain increment is shown to be expressed as:

$$ \Delta\varepsilon = \nabla^{\text{sym}}\mathcal{M}[\Delta t] + \nabla^{\text{sym}}\mathcal{N}[\mathcal{C} : \Delta\varepsilon^p(\Delta\varepsilon; S)], \quad (7) $$

where $S := (e^p, \varepsilon^p)$ is the current plastic state, $\Delta t$ is the increment of surface tractions (in our case $\Delta t = \Delta p e_3$ as we are in a normal contact situation). The function $\Delta\varepsilon^p(\Delta\varepsilon; S)$ represents the radial-return algorithm classically used in incremental plastic analysis (Simo and Hughes, 1998). Equation (7) is a non-linear equation that can be solved with the DF-SANE algorithm (La Cruz et al., 2006) which has the advantage of being jacobian-free.

The operators $\mathcal{M}$ and $\mathcal{N}$, which are at the heart of the method developed in (Frérot et al., 2019b), are linear integral operators which compute in $\mathcal{B}$ the displacement due to periodic distributions of surface traction and volume eigenstress respectively<fnref n="5" />. Their complete formulation and application in a discretized setting is extensively discussed in (Frérot et al., 2019b). The coupling between the elastic contact problem and the plasticity problem is done with a relaxed fixed point strategy (Frérot et al., 2019b; Jacq et al., 2002). The full implementation of the described solution method is also freely available in *Tamaas*.

## Comparison: rough surface

While both plasticity models are phenomenological, associated plasticity is soundly grounded in experimental observations (Bui, 1969) as well as thermodynamic principles (Reddy and Martin, 1994; Simo and Hughes, 1998), and expresses a macroscopic picture of dislocation systems at the micro-scale. This is not the case for the saturation models: they depend on the observation that the mean contact pressure saturates for spherical indentation (Tabor, 1951), which has been challenged by recent finite-element simulations (Krithivasan and Jackson, 2007; Song and Komvopoulos, 2013).

We aim here to provide a direct comparison for a rough surface between a perfectly plastic $J_2$ model and the saturation model. The rough surfaces we use throughout this work are self-affine random surfaces. Their power-spectrum density is defined as

$$ \phi(\mathbf{q}) = \begin{cases} C \left( \frac{q_l}{|\mathbf{q}|} \right)^{-2(H+1)} & q_l \le |\mathbf{q}| \le q_s, \\ 0 & \text{otherwise} \end{cases} \quad (8) $$

where $q_l, q_s$ are the spatial frequencies associated to the long cutoff wavelength $\lambda_l$ and short cutoff wavelength $\lambda_s$ respectively, while $H$ is the Hurst exponent. We define two surfaces with Hurst exponent $H = 0.8$: $S_1$ has a rather narrow spectrum with $L/\lambda_l = 3$ and $\lambda_l/\lambda_s = 3$; $S_2$ has a broader spectrum and a larger system size, with $L/\lambda_l = 16$ and $\lambda_l/\lambda_s = 8$. The root-mean-square of slopes $h'_{\text{rms}}$ is 0.07 for $S_1$ and 0.1 for $S_2$. Since $S_1$ is used for visual comparison, the discretization size $\Delta l$ is chosen such that $\lambda_s/\Delta l = 3$, while a depth of $L/5$ is represented with 64 points (totaling 243 × 243 × 64 points). $S_2$ is used to study crack nucleation at the interface, and therefore has a finer discretization $\lambda_s/\Delta l \approx 5.7$, with a depth $L/5$

<footnote n="4">⁴In this approach the residual displacement is directly computed from $\varepsilon^p$, whereas in the saturation plasticity model it is merely *assumed* equal to the negative gap.</footnote>

<footnote n="5">⁵For reference, we can express with $\mathcal{M}$ and $\mathcal{N}$ both the surface vertical displacement due to an applied pressure $\overline{M}[p] = \mathcal{M}[p \cdot e_3]|_{\partial\mathcal{B}} \cdot e_3$ and the residual vertical displacement $\bar{u}_3^p = \mathcal{N}[\mathcal{C} : \varepsilon^p]|_{\partial\mathcal{B}} \cdot e_3$.</footnote>
represented with 32 points (totaling 729 × 729 × 32). Naturally, the saturated simulations do not account for volumetric behavior, so the number of points are 243 × 243 and 729 × 729 for S₁ and S₂ respectively.

The normal loads applied in both models are adimensionalized by $W_0 = E^*L^2h'_{rms}$. In elasticity, this normalization collapses load (W) vs. true contact area ($A_c$) for different values of $h'_{rms}$ (but the same spectrum parameters) (Bush et al., 1975; Hyun et al., 2004). We do not intend here to modify the spectrum parameters but merely scale $h'_{rms}$, which therefore becomes a non-dimensional measure of surface summit amplitude. This is convenient to compare the results of the $J_2$ and saturated models to an elastic reference, as we expect the contact behavior to depend on the surface peak amplitude because of plasticity.

<figure><img src="image_2.png" /><figcaption>Figure 1: Pressure profiles for $J_2$ plasticity (a), saturation (b) and elasticity (c) models. Although fig. (a) shows that the local pressure can exceed values of $3\sigma_y$, the average pressure in contacts $\langle p \rangle = W/A_c$ is closer to $2\sigma_y$, whereas the saturated model gives an average of $2.5\sigma_y$ with large saturated portions of the micro-contacts. In this case, the normal load is $W/W_0$ is $6.5 \cdot 10^{-2}$, and the saturated model predicts a contact area 20% smaller than the $J_2$ prediction. As a result the connectivity of micro-contacts is different between the two models.</figcaption></figure>

Figure 1 shows the contact pressures for $J_2$ (fig. 1a) and saturated (fig. 1b) plasticity, as well as elasticity (fig. 1c), at the load $W/W_0 = 6.5 \cdot 10^{-2}$. The yield stress is $\sigma_y = 10^{-2} \cdot E$ and the surface used is $S_1$. The total contact ratio is 25% for $J_2$ plasticity, 20% for saturation and 15% for elasticity, resulting in about 20% error in the contact area of the saturated model. Moreover, while the maximum pressure in the $J_2$ model exceeds $3\sigma_y$ (cf. fig. 1a), the average pressure on micro-contacts is closer to $2\sigma_y$, which the saturation model fails to capture with an average of $2.5\sigma_y$. Local features of the contact patches also differ due to the three models being in different contact stages.

Figure 2 shows the probability density function of the surface pressures for the three models, with the addition of reference data from Pei et al. (2005) (fig. 10a) which results from a $J_2$ plasticity criterion used in a finite-element approach. The features of the curve corresponding to the $J_2$ plasticity models are not qualitatively reproduced by the saturation model: the peak at $p = 2.5\sigma_y$ is non-existent and as expected the distribution of the saturated pressures tends to a Dirac at $p = p_m$ whereas the $J_2$ distribution tends to zero. The difference between our results and those of Pei et al. (2005) can be explained by the coarseness of the mesh they used, as well as different spectrum parameters (e.g. $H = 0.5$ in their case).

Figure 3 shows the contact area ratio and the secant slope of the load/contact area relationship, which is a quantity that has been extensively studied in elastic contact (Bush et al., 1975; Campañá and Müser, 2007; Hyun et al., 2004; Pastewka and Robbins, 2014; Persson, 2001b; Yastrebov et al., 2015, 2017). Figure 3a shows the result for surface $S_1$: as in fig. 1, the saturated model underestimates the true contact area compared to the $J_2$ plasticity model. All three model nonetheless exhibit a similar behavior, which is not the case in fig. 3b, where we show the same data for surface $S_2$. Due to computational issues, the $J_2$ approach in this case includes hardening with $E_h/E = 0.05$. This locally increases the yield stress, but in our simulations the cumulated plastic strain averaged over plastic zones never exceeds 1%, which means that the yield stress only increases by 5% on average. One can note that this time the saturated model overestimates the true contact area, meaning that the average pressure on contacts in the $J_2$ case is larger than $3\sigma_y$. This is consistent with the findings of Pei et al. (2005) and Krithivasan and Jackson (2007) with elastic perfectly-plastic materials<fnref n="6" />. Unlike for the saturated and elastic models, the secant slope of the $J_2$ model seems to stay constant, as observed by Pei et al. (2005). This observation, and in general the

<footnote n="6">⁶We have $\langle p \rangle \approx 3.6\sigma_y$, which is 20% larger than the common $3\sigma_y$. Since hardening only increases $\sigma_y$ by 5% on average, this difference in $\langle p \rangle$ cannot be explained by hardening alone.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2: **Probability density function of surface pressures.** Neither the elastic nor the saturated models qualitatively reproduce the pressure distribution of the elastic-plastic model. As expected, the pressure distribution of the saturated model tends to a Dirac distribution for $p = p_m$, whereas the distribution for the elastic-plastic model tends to zero. Note that the results of Pei et al. (2005) have been renormalized (cf. companion notebook (Frérot, 2019)).</figcaption></figure>

dependency of the contact area to load relationship on surface and plasticity parameters warrants further investigations, similar to Pei et al. (2005), with refined simulations, which unfortunately go out of the scope of this work.

To study macroscopic quantities (such as the true contact area), the saturation model is only suitable if one has reference data to fit $p_m$. However, it fails to qualitatively reproduce local surface quantities (as can be seen in figs. 1 and 2), as well as global quantities in certain cases (e.g. ??b), in addition to providing no information on the complete stress state at and below the contact surface. The saturation model may be useful in applications where quantitative errors in the contact area magnitude or its topography may be accepted, only to obtain approximate qualitative relations. It can however give unequivocally wrong results when local stress based quantities drive the phenomenon one wishes to study. We shall see in the next section that this can be the case in wear modeling when we consider the crack nucleation process in an elastic-plastic rough contact.

## 2. Crack nucleation in rough surface contact

Most of the atomistic investigations of adhesive wear processes use geometries that contain stress concentrators (Aghababaei et al., 2016; Brink and Molinari, 2019; Milanese et al., 2019), such that in their model system the debris formation is only controlled by the Griffith energy balance. However, without a defect/stress concentration, a crack described by linear elastic fracture mechanics cannot nucleate. In a half-space geometry, with plasticity constitutive behavior, no such concentration can exist. For that reason, it is necessary to introduce a critical nucleation tensile stress $\sigma_c$.

In linear elasticity the stresses are unbounded and depend linearly on the applied load. The picture however changes in plasticity, as one can expect the tensile stress to saturate, possibly below $\sigma_c$, preventing crack nucleation altogether. Brink and Molinari (2019) have also shown that the resistance to shear of the contact junction plays a fundamental role in the wear particle formation. From a stress perspective, there is a competition between $\sigma_c$ and $\tau_j$ (the junction shear resistance) for the formation of a crack: if the junction is strong enough, the maximum tensile stress may, under conditions depending on $\sigma_y$, reach $\sigma_c$ and nucleate a crack. Conversely, if the junction is weaker, its slip may prevent the tensile stress from reaching $\sigma_c$.

We investigate the interplay of these effects in the contact of the rigid self-affine rough surface $S_2$ (defined previously) with a solid that we have discretized with 729 × 729 points for the elastic and saturation models (with $\sigma_y/E = 10^{-2}$), and 729 × 729 × 32 for the $J_2$ model ($\sigma_y/E = 10^{-2}$, $E_h/E = 5 \cdot 10^{-2}$ as in the previous section). The applied normal mean pressure varies between $10^{-2}h'_{\text{rms}}E$ and $8 \cdot 10^{-2}h'_{\text{rms}}E$ for the elastic case and between $10^{-2}h'_{\text{rms}}E$ and $4 \cdot 10^{-2}h'_{\text{rms}}E$ for the saturation and $J_2$ models. Since
<figure><img src="image_4.png" /><figcaption>Figure 3: **Contact ratio and secant slope of the load/contact area relationship.** Figure (a) shows results for surface S₁ and elastic perfectly-plastic materials. Figure (b) shows results for surface S₂ (which has a broader spectrum), and 5% hardening for the J₂ model. In fig. (a), the saturated model overestimates the true contact area, whereas it is underestimated in fig. (b). This suggests that adjusting the saturated model to quantitatively reproduce the true contact area may only be done on a case-by-case basis.</figcaption></figure>

we investigate the effect of the junction shear strength, we also apply a shear stress on the contacts at constant normal load. There is a simple correspondence between the applied shear stress and the junction strength $\tau_j$: any shear loading larger than $\tau_j$ should not modify the stress state of the system since all contacts should be slipping and the stresses should not increase. We therefore interchangeably refer to the applied shear and the junction strength as $\tau_j$. In linear elasticity, the application of a constant shear stress on a patch of the surface creates a stress singularity at the edge of the patch because the derivative of the surface tangential displacement diverges (Menga and Carbone, 2019). In a physical system, a small amount of slip and rearrangement of the solids would occur at the edge of the contact junction so that the shear stress carried should be reduced on this zone. We therefore regularize<fnref n="7" /> the constant shear distribution over a transition zone of width $\epsilon_\tau$ small compared to the smallest surface wavelength ($\epsilon_\tau \approx \lambda_s/3$), removing any numerical discrepancy due to the singularity.

We call a crack nucleation site a connected zone of the surface where the largest eigenvalue of the Cauchy stress tensor $\sigma_I$ is larger than $\sigma_c$. Although the wear particle formation process is deterministic, the inherent randomness of the rough surfaces makes the process epistemically random (Frérot et al., 2018). We apply a similar concept here and study the probability that a given contact nucleates a crack. Although we focus here on crack nucleation, and do not model crack propagation, we assume that the propagation of a crack has a minimal influence on the nucleation of other cracks. This assumption holds for small true contact areas, where the micro-contacts are dilute, which would correspond to a mild wear regime<fnref n="8" />. We shall see that despite this assumption, our model still includes elastic interactions between contacts and plastic zones, which have an effect on the crack nucleation sites in close proximity, particularly when shear is applied.

Because of the complex topography of the micro-contacts, there is no one-to-one correspondence between crack nucleation sites and contacts and the probability that a given contact nucleates a crack cannot be explicitly evaluated. This can be seen in fig. 4, where we show in grey the true contact area, in red the plastic zones and in black the crack nucleation sites (i.e. the zones where $\sigma_I > \sigma_c$) for $\sigma_c = 0.1h'_{\text{rms}}E^*$ and

<footnote n="7">⁷Regularization is done by convolution with a function of the form $\phi_{\varepsilon}(r) = \exp(-1/(1 - (r/\varepsilon)))\varepsilon^2$ (David and Gosselet, 2015).</footnote>

<footnote n="8">⁸It is possible that this assumption holds for large contact area ratios (Sevostianov and Kachanov, 2012), but current numerical methods do not allow for detailed contact simulations with explicit crack propagation modeling to verify if this is the case.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 4: Crack nucleation sites for the saturation pressure (a) and $J_2$ plasticity (b) models. An “upwards” shear is applied on each contact. The true contact area is shown in light grey, the plastic zones in red and the crack nucleation sites in black. We can see that the $J_2$ model has more crack nucleation sites than the saturation model. Since both models give approximately the same true contact area, this discrepancy must be due to plastic residual deformations which are not represented in the saturation approach. Figure (b) shows that there can be multiple crack nucleation sites per contact.</figcaption></figure>

upwards applied shear stress $\tau_j = 2 \cdot 10^{-3} E$. Figure 4a shows the result for the saturation plasticity model and fig. 4b shows the $J_2$ plasticity model. One can easily see on the latter that a single micro-contact may have several crack sites. We also recognize the expected crescent shape crack nucleation on some contacts. Finally, the number of cracks is larger in the $J_2$ model than in the saturation plasticity model. Since the contact areas predicted by both models are essentially the same, this difference can only be explained by the residual plastic deformations, which are not modeled in the saturation plasticity approach.

Because we still want to investigate a quantity akin to the probability that a contact nucleates a crack, we propose an adimensional measure called the crack nucleation likelihood (CNL) given by $A_0 \cdot n_{\text{crack}}/A_c$, where $A_0 = L^2$ is the apparent contact area and $n_{\text{crack}}$ is the number of crack nucleation sites (i.e. the number of connected black zones in fig. 4). The CNL is conceptually a normalization of the number of crack nucleation sites by the density of contacts and relates to the probability of crack nucleation at a contact (this relationship will be detailed later on in this article).

Figures 5 and 6 show the CNL as a function of $\sigma_c$ when the normal load and the applied shear stress/junction strength are respectively varied, for an elastic (a), a saturated and a $J_2$ contact (b). On fig. 5a, the CNL curves are shifted to the right with larger normal loads. This means that for a fixed value of $\sigma_c$, the CNL increases exponentially when the load is increased. The vertical lines indicate a quantity $\sigma_{i/c}$, which is the stress for which the probability density function of $\sigma_I$ on the whole surface is maximum: in other words it is the most frequent stress value, and is typically found between contacts (hence the term “inter-contact stress”). The horizontal shift in the CNL curves corresponds to $\sigma_{i/c}$. Since the latter depends on the spatial proximity of contacts, the CNL must depend on elastic interactions between contacts. Figure 5b shows that the two plasticity models have widely different behavior: the crack nucleation likelihood is much higher in the case of $J_2$ plasticity, but also decays faster for values of $\sigma_c > \sigma_y$. It is surprising that despite being plastic with some hardening, the $J_2$ model is more likely to lead to surface cracks than the saturation model. This may be a first step towards resolving the paradox highlighted in introduction.

Figure 6 shows the CNL when the junction strength increases (for the last normal load of fig. 5). Unlike previously, the elastic results are not shifted to the right when $\tau_j$ increases but are instead scaled rightward. The same happens with the saturation plasticity model, whereas the $J_2$ CNL seems relatively unaffected by $\tau_j$. This is due to plastic deformations having reached the contact surface (cf. fig. 4) and preventing an increase of $\sigma_I$ as fast as the elastic and saturation models, although hardening still allows some increase at a lower rate. In order to rationalize these results and provide evidence for the conclusions we have induced from our multi-asperity simulations, we now study the crack nucleation process for a single asperity.
<figure><img src="image_6.png" /><figcaption>Figure 5: Crack nucleation likelihood (CNL) as a function of $\sigma_c$ and normal load. Figure (a) shows the results for an elastic contact, while (b) shows the saturation and $J_2$ plasticity results. In (a) the CNL curves are uniformly shifted to the right when the load is increased, indicating an exponential increase in the CNL. The normalization of $n_{\text{crack}}$ by the true contact area makes explicit that this increase is due to stronger elastic interactions between contacts. The magnitude of the shift is given by the most frequent value of $\sigma_I$ called $\sigma_{I/c}$. The two plasticity models in (b) have drastically different behavior: the crack nucleation is much more likely in the $J_2$ approach because of plastic residual stresses, and the CNL experiences a faster decay for values of $\sigma_c > \sigma_y$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Crack nucleation likelihood (CNL) as a function of $\sigma_c$ and junction strength. As for fig. 5, (a) shows the elastic model and (b) the plastic models. Unlike fig. 5a, the CNL for the elastic model is *scaled* to the right when the junction strength increases. The same can be said of the curves corresponding to the saturation plasticity in (b), but not of the $J_2$ curves, which are relatively insensitive to changes in $\tau_j$. This is due to plasticity preventing increases in $\sigma_I$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: **Maximum tensile stress as a function of applied shear stress across single asperity contact with $J_2$ plasticity.** The loading curves consist of two parts: an initial linear loading and a non-linear saturation of $\sigma_I$. Although the initial slope is linear, the loading is not elastic as plastic deformations still evolve in the system. Their influence on the maximum tensile stress at the surface is however minimal. A sufficient shear stress may eventually create plastic strain at the surface, which causes a transition in the loading curve. Moreover, given a low enough yield stress, the initial indentation may cause surface yield, as seen for the curve where $\sigma_y/E = 2 \cdot 10^{-2}$.</figcaption></figure>

### 3. Single asperity crack nucleation

To investigate the effect of plasticity in the competition between $\sigma_c$ and $\tau_j$, we simulate a spherical indenter of radius $R$ pushed onto an elastic-perfectly-plastic solid. The resulting contact junction is then subjected to a shear distribution (with $\varepsilon_\tau = R/64$), and the principal tensile stress $\sigma_I$ at the surface is recorded.

Figure 7 shows the maximum tensile stress $\sigma_I$ as a function of the applied shear stress $\tau_j$ across the contact. The different curves correspond to different yield stresses, with the dashed lines indicating the value of $\sigma_y$ for reference. Stresses here are normalized by the maximum tensile stress in Hertz contact $\sigma_I^{\text{Hertz}} = (1 - 2\nu)p_0/3$, with $p_0$ being the maximum hertzian contact pressure (Johnson, 1985). We can observe that the initial tensile stress (without applied shear) depends on the amount of plastic deformation: if $\sigma_y$ is decreased (or conversely the applied load increases), the initial tensile stress at the edge of the contact is higher. This is due to the residual stresses created by the plastic deformations that accommodate the indentation: the localized nature of the plastic strains causes the unloaded equilibrium position to not be stress-free. The additional stresses are tensile and add to the stress on the contact rim. In the case of $\sigma_y/E = 2 \cdot 10^{-2}$, the plastic zone has reached the surface and the von Mises stress at the edge of contact has reached $\sigma_y$ by indentation alone (not shown here).

If the von Mises stress at the surface is below $\sigma_y$, the application of a shear stress will cause an elastic loading phase, offset by the initial $\sigma_I$ value, as seen for the higher values of $\sigma_y/E$. The loading continues until $\sigma_I$ reaches values close to $\sigma_y$, as indicated by the dashed lines. Since the stress state at the edge of contact is triaxial, the maximum value $\sigma_I$ can reach is not $\sigma_y$, as is seen for the most plastic case. After a certain point, eq. (7) becomes unsolvable because a plastic failure mechanism develops (Drucker and Prager, 1952): we supposed that further loading will not increase the value of $\sigma_I$. Of course, for a hardening material $\sigma_I$ should not saturate and instead increase further at a lower rate.

Figure 8 shows the competition between the junction strength $\tau_j$ and the tensile strength $\sigma_c$ for different values of $\sigma_y$ (first row no hardening, i.e. $E_h/E = 0$) and $E_h$ (second row with $\sigma_y/E = 2 \cdot 10^{-2}$) the hardening modulus. The zones below the white curve show for which values of $(\sigma_c, \tau_j)$ a crack may nucleate, and the zones above show when the interface breaks (i.e. slips) before crack nucleation. The dashed lines indicate the values of $\sigma_y$. We can see that for materials with high yield there is an affine boundary between the two mechanisms, which shows their competition. Of course materials with $\sigma_c < \sigma_y$ are not plastic, so such transition can only happen for brittle materials. It does however exist for hardening materials, as
<figure><img src="image_9.png" /><figcaption>Figure 8: Failure regimes for a sheared spherical indentation. For the first row of graphs, the ratio $\sigma_y/E$ is varied for a perfectly plastic material. For the second row, the hardening ratio $E_h/E$ is varied for a yield ratio of $\sigma_y/E = 2 \cdot 10^{-2}$. On each graph the dashed line shows the yield stress. The white curve marks the transition between failure driven by slip rupture of the junction and failure by crack nucleation. The competition between the junction strength $\tau_j$ and the critical stress $\sigma_c$ is influenced by $\sigma_y$ because of the saturation effect shown in fig. 7. While plasticity gives a failure mechanism independent of $\tau_j$, hardening allows the tensile stress to grow past the initial yield limit, giving a linear transition between failure mechanisms.</figcaption></figure>

opposed to perfectly plastic ones. With no hardening, the failure mechanism is purely determined by the value of $\sigma_c/\sigma_I^{\text{Hertz}}$, which depends on the applied load.

With this single-asperity analysis, we have explained why the crack nucleation likelihood is higher in the $J_2$ model for no applied shear stress: the plastic residual deformations cause tensile stresses which combine with the contact stresses and increase $\sigma_I$, thus increasing the CNL. This does not occur in the saturated model because it ignores plastic residual deformations in the stress computation. We have also explained why the $J_2$ CNL is relatively insensitive to the applied shear/junction strength: plasticity has a saturation effect on $\sigma_I$: when the system is sheared, the rim of contacts is in the plastic regime, and the increase in $\sigma_I$ is purely driven by hardening, which in the case of fig. 6 is only 5% of the Young's modulus. In order to rationalize the other aspects of the CNL highlighted by figs. 5 and 6, we resort to a statistical model for multi-asperity contact.

## 4. Multi-asperities

One can now apply a statistical approach to estimate the proportion of contacts that nucleate cracks in a multi-asperity setting. We thereafter use a Greenwood-Williamson (GW) model (Greenwood and Williamson, 1966) with an exponential distribution of asperity heights to obtain simple, qualitative, analytical results to help rationalize the data of figs. 5 and 6. Since the asperities are randomly distributed (all with the same radius $R$), $\sigma_I^{\text{Hertz}}$ becomes in turn a random variable. There is however a significant difference between the single asperity case we have studied and the multi-asperity setting. Because of elastic interactions, the tensile stress at the edge of a contact depends on the proximity and magnitude of the neighboring contacts. In a traditional GW approach, contacts are independent of each other. We assume this is the case, but that the stress state is determined by the local contact with an additional contribution $\sigma_{i/c}$, the inter-contact stress, determined from the neighboring contacts. The radial stress outside the area of a single contact of radius $a$ is given by (Johnson, 1985):

$$ \sigma_r(r) = \frac{1-2\nu}{3} \cdot \frac{a^2}{r^2} p_0 = \kappa \frac{sR}{r^2} (z^* - h)^{\frac{3}{2}}, \qquad (9) $$
with

$$ \kappa := \frac{2(1 - 2\nu)E^*}{3\pi} \cdot \sqrt{\frac{s}{R}}, \qquad (10) $$

where $r$ is the euclidean distance from the contact center, $z^* = z/s$ is the asperity height random variable normalized by the standard deviation of heights $s$ and $h$ is the normalized surface approach. We assume a spatial asperity density $\eta$ and a contact density $\eta_c = \eta e^{-h} = A_c/(A_0\pi sR)$ (Greenwood and Williamson, 1966), with $A_c$ being the true contact area. To compute the inter-contact stress $\sigma_{i/c}$, we assume the stress state outside each contact is given by the mean of eq. (9), averaged over contacting asperities, and we compute the largest stress eigenvalue at the center of a series of circles whose diameters are multiples of $d_c = 1/\sqrt{\eta_c}$ which is the characteristic distance between contacts. The calculation process is detailed in appendix Appendix A and leads to the following expression:

$$ \sigma_{i/c} = 3\kappa \frac{\xi}{\sqrt{\pi}} \cdot \frac{A_c}{A_0}, \qquad (11) $$

with $\xi \le \sqrt{3\zeta(3)} \approx 1.9$ and $\zeta$ is the Riemann Zeta function. We use this upper bound for $\sigma_{i/c}$ in the rest of this work. We can see that $\sigma_{i/c}$ depends linearly on the contact ratio, and therefore is linear with the load. We suppose that $\sigma_{i/c}$ acts as a “background” stress, and that the maximum tensile stress is the sum of the local contact tensile stress $\sigma_I^{\text{Hertz}}$ and the inter-contact stress. We can therefore quantify the probability that a contact nucleates a crack:

$$ \begin{aligned} P_{\text{crack}} &= P \left( \frac{\sigma_c}{\sigma_I^{\text{Hertz}} + \sigma_{i/c}} < \omega(\tau_j) \mid z^* - h \ge 0 \right) \\ &= \exp \left( - \left( \frac{\sigma_c/\omega(\tau_j) - \sigma_{i/c}}{\kappa} \right)^2 \right), \end{aligned} \qquad (12) $$

where $\omega$ is the function describing the failure mechanism transition (white line in fig. 8). The calculation details of eq. (12) can be found in appendix Appendix B.

*Remark 1.* The crack nucleation likelihood $A_0 \cdot n_{\text{crack}}/A_c$ and $P_{\text{crack}}$ are related in our GW approach: indeed $n_{\text{crack}} = P_{\text{crack}}\eta_c A_0 = P_{\text{crack}}A_c/\pi sR$, cf. (Greenwood and Williamson, 1966, p. 303).

*Remark 2.* The quantity $\sigma_c/\omega(\tau_j) - \sigma_{i/c}$ strikingly explains the features of figs. 5 and 6 for the elastic model. In fig. 5, the curves are shifted to the right as the load increases, which is apparent in eqs. (11) and (12): $\sigma_{i/c}$ increases linearly with the load, and thus causes a rightward shift in the graph of the CNL. Similarly, as $\omega$ is linear in $\tau_j$, $P_{\text{crack}}$ is scaled horizontally, which can also be seen in the CNL on fig. 6.

*Remark 3.* When $\sigma_c/\omega(\tau_j) = \sigma_{i/c}$ the probability is one, meaning that all contacts, regardless of size, nucleate cracks; in other words the whole surface should be cracking in a catastrophic breakdown. This does not happen in practice, as the normal loading process should nucleate and propagate cracks at single asperities before the breakdown is reached, thus relaxing the tensile stresses in the system.

## Comparison to a rough surface

We wish to assess the validity of the above developments with simulations of self-affine rough surface contact. Because of the simplifying assumptions of a GW model, we do not hope to establish a quantitative agreement, especially since the asperity curvature is not unequivocally defined on a self-affine rough surface (Nayak, 1971). Instead, we will focus on the qualitative relations between $P_{\text{crack}}$, $\sigma_c$, $\sigma_{i/c}$ and $\tau_j$ highlighted above.

**Elasticity Results.** We first consider $\sigma_{i/c}$ for a rough surface. Recall that the inter-contact stress in a rough contact is the most frequent value of $\sigma_I$, i.e. the value for which the probability density function $p_{\sigma_I}$ of the surface tensile stress is maximum. This is illustrated in the inset of fig. 9. In the latter, we plot the evolution of $\sigma_{i/c}$ for the elastic rough contact defined previously and for eq. (11). We can observe that both curves behave linearly with the contact area ratio, with different slopes. To compute the value of $\kappa$, which depends on $\sqrt{s/R}$, cf. eq. (10), we have used Nayak’s approach (Nayak, 1971) to estimate the mean curvature radius of the zones of the rough surfaces in contact, i.e. $R = \sqrt{3/m_4}/I(z_c^*)$, where $m_4$ is the
<figure><img src="image_10.png" /><figcaption>Figure 9: **Inter-contact stress as a function of contact ratio.** Inset shows the probability density function of the largest stress eigenvalue at the surface. The stress value corresponding to the peak in the probability density is defined as the inter-contact stress $\sigma_{i/c}$. We see that in the rough surface simulation and the Greenwood-Williamson model $\sigma_{i/c}$ evolves linearly with the contact ratio. The value of $\kappa$ is estimated from the mean curvature of contacting summits in the rough surface.</figcaption></figure>

fourth moment of the surface spectrum, $z_c^*$ is the normalized height of the surface in contact and $I$ is a function defined in (Nayak, 1971). As can be seen in fig. 9, the slope of the GW curve is approximately constant, showing the weak dependency on $z_c^*$. Note that although the values of $R$ and $A_c/A_0$ for the GW model are informed from the rough surface simulation, there is no fit parameter to the GW prediction. The lack of quantitative agreement between the theoretical approach and the rough contact simulation shows the prediction limit of asperity-based models.

As previously mentioned, $P_{\text{crack}}$ is not directly measurable on a rough contact interface (cf. fig. 4). However, the crack nucleation likelihood acts as an alternative measure for $P_{\text{crack}}$, cf. remark 1. Figure 10a, confirms that this is indeed the case and we find the squared exponential dependency predicted by eq. (12), with all curves collapsed due to the shift caused by $\sigma_{i/c}$ (recall that $\omega(0) = 1$). Figure 10b, on the other hand, shows that the CNL does not follow eq. (12) for non-zero $\tau_j$. While each curve remains close to a straight line, they do not overlap, but seem to converge to a master curve. More strikingly, the CNL decreases as $\tau_j$ increases, meaning that $\omega(\tau_j)$, which was computed directly from the data of fig. 7, over-normalizes the data. This is again due to interactions between asperities. For our “single-asperity” analysis, because of periodicity, we in fact consider many interacting asperities on a square lattice, each separated by a distance $L$. When shear is applied, a positive $\sigma_1$ is created at the trailing edge of the contact and a negative $\sigma_1$ appears at the leading edge. Because the periodic images are equidistant and far apart, they weakly affect the stress distribution in the vicinity of the contact. However, when two contacts are close to each other, creating local anisotropy, the inter-contact stress distribution of each asperity is compensated by the other. The positive peak in $\sigma_1$ at the trailing edge of one contact is then reduced, thereby reducing $P_{\text{crack}}$ as seen in fig. 10b. This phenomena is akin to the crack shielding mechanism uncovered by Aghababaei et al. (2018).

*Plasticity Results.* As for the elastic model, we compare the inter-contact stress computed from the plastic rough contact simulations to our GW approach. Figure 11 shows the results for both the saturation and $J_2$ plasticity. Compared to the elastic results, the slope of the plastic models is smaller. The $J_2$ plasticity model has the smallest slope, showing that residual plastic deformations play a role in the inter-contact stress. One should note that the plastic model in fig. 11 includes hardening, hence the reduced contact ratio. It seems both models still give a linear dependency of $\sigma_{i/c}$ on the contact ratio, although some more data may be required to draw an affirmative conclusion in this regard.

Finally, fig. 12 shows the data of fig. 5b normalized to compare the results to eq. (12). While the saturation model seems to follow our GW prediction (which is based on elasticity assumption), it is clear that the $J_2$ model does not conform to our scaling predictions for $P_{\text{crack}}$. However, in light of fig. 11, it is
<figure><img src="image_11.png" /><figcaption>Figure 10: Crack nucleation likelihood as a function of re-normalized critical stress, normal load and junction shear strength. The material in contact behaves elastically. We normalize by $h'_{\text{rms}}E^*$ instead of $\kappa$ because of the difficulty of defining asperity curvature on a rough surface. One can see on figure (a) that the curves corresponding to different normal loads are collapsed on a straight line, showing that the CNL does indeed follow the scaling established in eq. (12). When the strength of the junction is taken into account (figure (b)), or, equivalently, if a shear stress is applied, we observe qualitative deviations from eq. (12). There is a decrease in crack event density due to the interference of close contacts, which tends to unload the tensile stresses at the trailing edge of leading contacts.</figcaption></figure>

interesting to note that although $\sigma_{i/c}$ is lowest for the $J_2$ approach (indicating less interactions between contacts), the latter has the largest crack nucleation likelihood, because of the local effect of plastic residual deformations. As shown in fig. 7, this local effect of plastic deformations is stronger the more ductile a material is, as expected from experimental data which shows that softer materials wear more.

## Conclusion

We have investigated in this work the nucleation of cracks at an elastic-plastic rough contact interface. This was motivated by the necessity for an accurate description of the process of crack nucleation for adhesive wear that includes plasticity. By comparing a classical $J_2$ plasticity model and a saturation plasticity approach commonly used in tribology, we have concluded that the latter can only qualitatively reproduce the true contact (in some cases) area and fails to give satisfactory results on local quantities. This can be seen in the crack nucleation likelihood, which is much higher in the $J_2$ approach. This is caused by plastic residual deformations which increase tensile stresses in the vicinity of contacts. In this regard, the saturation model is not applicable to study crack nucleation because it does not capture plastic deformations. We show with a single asperity analysis that more ductile materials can have larger surface tensile stresses and nucleate more cracks at the interface.

We have also showed that elastic interactions play a role in the crack nucleation likelihood of a single contact. They may increase the latter through proximity of contacts, or decrease it in shearing by elastic shielding. This was further supported by an analytical approach based on a Greenwood-Williamson model modified to take interactions into account.

Future improvements of this work could include adhesive contact at the interface due to surface energy. Attractive forces usually produce large stresses at the surface, which may also contribute to the development of plastic zones in the vicinity of contact edges.

*Acknowledgments.* All authors acknowledge the insightful discussions with T. Brink and E. Milanese, as well as the financial support of the Swiss National Science Foundation (grant #162569, “Contact mechanics of rough surfaces”).

*Supplementary data.* All codes used in this work are available on Zenodo (Frérot, 2019; Frérot et al., 2019a).
<figure><img src="image_12.png" /><figcaption>Figure 11: **Inter-contact stress as a function of contact ratio for the saturation and $J_2$ plasticity models.** While the curves do not match the analytical GW approach, their slopes are smaller than in the elastic case. The $J_2$ model shows the smallest slope, indicating that the stresses due to plastic residual deformations have an influence on the inter-contact stress and actually reduce it compared to the underlying elastic stresses of the saturation pressure model.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: **Crack nucleation likelihood as a function of re-normalized critical stress and normal load for the saturation and $J_2$ plasticity models.** The pressure saturation model reproduces a scaling similar to the elastic case in fig. 10a, while the $J_2$ model shows a fundamentally different behavior. The crack density is higher in the plastic case because of the additional tensile stresses caused at the edge of contacts by residual plastic deformations.</figcaption></figure>
Appendix A. Inter-contact stress computation

We first compute $\overline{\sigma_r}(r)$, which is the average radial stress outside contacts:

$$ 
\begin{aligned}
\overline{\sigma}_r(r) &= \frac{1}{P(z^* \ge h)} \int_h^\infty \kappa \frac{sR}{r^2} (z^* - h)^{\frac{3}{2}} e^{-z^*} dz^* \\
&= \kappa \frac{3\sqrt{\pi}}{4} \cdot \frac{sR}{r^2}.
\end{aligned}
$$

We then suppose for simplicity that all contacts have the same radial stress $\overline{\sigma}_r(r)$. Accordingly, their hoop stress is $\overline{\sigma}_\theta(r) = -\overline{\sigma}_r(r)$ (Johnson, 1985). We assume that all contacts are spatially uniformly distributed with density $\eta_c$, so that the characteristic distance between contacts is $d_c = 1/\sqrt{\eta_c}$. We divide the infinite surface into concentric rings of width $d_c$ and diameters $d_i \in \{d_c, 2d_c, 3d_c, \dots\}$. Each ring can be reduced to a circle of contacts with linear density $\sqrt{\eta_c}$. We now wish to compute for the sum of all circles of diameter $d_{1,2,\dots}$ the largest stress eigenvalue at the center. For a single contact positioned at an angle $\theta$ on a circle of diameter $d_i$, the stress state in Cartesian coordinates is:

$$ \sigma = \overline{\sigma}_r(d_i/2) \begin{pmatrix} \cos(2\theta) & -\sin(2\theta) \\ -\sin(2\theta) & -\cos(2\theta) \end{pmatrix}. $$

For the circle number $i$, the expected number of contacts is $\pi i d_c \sqrt{\eta_c} = i \cdot \pi$, but we simplify by assuming the expected number to be $n_i := 3i$. The total stress state at the center, summing all contacts per circle and all circles:

$$ \sigma = \sum_{i=1}^{\infty} \sigma_r^i \sum_{k=1}^{n_i} \begin{pmatrix} \cos(2\theta_k) & -\sin(2\theta_k) \\ -\sin(2\theta_k) & -\cos(2\theta_k) \end{pmatrix}, $$

with $\sigma_r^i := \overline{\sigma}_r(id_c/2)$. The largest eigenvalue of $\sigma$ is given by:

$$ 
\begin{aligned}
\lambda^2 &= \left( \sum_{i=1}^{\infty} \sigma_r^i \sum_{k=1}^{n_i} \cos(2\theta_k) \right)^2 + \left( \sum_{i=1}^{\infty} \sigma_r^i \sum_{k=1}^{n_i} \sin(2\theta_k) \right)^2 \\
&= \sum_{i=1}^{\infty} (\sigma_r^i)^2 (c_i^2 + s_i^2) + \sum_{i<j}^{\infty} \sigma_r^i \sigma_r^j (c_i c_j + s_i s_j),
\end{aligned}
$$

where:

$$ c_i := \sum_{k=1}^{n_i} \cos(2\theta_k), \quad s_i := \sum_{k=1}^{n_i} \sin(2\theta_k) $$

The angular position of each contact is assumed to be uniformly distributed in $[0, 2\pi]$. The expected value of $\lambda$, which gives the inter-contact stress, is given by:

$$ \sigma_{i/c} = \mathbb{E}[\lambda] \le \sqrt{\mathbb{E}[\lambda^2]} = \sqrt{\sum_{i=1}^{\infty} (\sigma_r^i)^2 \mathbb{E}[c_i^2 + s_i^2] + \sum_{i<j}^{\infty} \sigma_r^i \sigma_r^j \mathbb{E}[c_i c_j + s_i s_j]}. $$
where we have used Jensen's inequality for a simple estimation. We note that:

$$ 
\begin{aligned}
c_i^2 + s_i^2 &= \left( \sum_{k=1}^{n_i} \cos(2\theta_k) \right)^2 + \left( \sum_{k=1}^{n_i} \sin(2\theta_k) \right)^2 \\
&= \sum_{k=1}^{n_i} \cos^2(2\theta_k) + 2 \sum_{k<l}^{n_i} \cos(2\theta_k) \cos(2\theta_l) + \sum_{k=1}^{n_i} \sin^2(2\theta_k) + 2 \sum_{k<l}^{n_i} \sin(2\theta_k) \sin(2\theta_l) \\
&= n_i + 2 \sum_{k<l}^{n_i} \cos(2(\theta_k - \theta_l))
\end{aligned}
$$

$$ 
\begin{aligned}
c_i c_j + s_i s_j &= \sum_{k=1}^{n_i} \cos(2\theta_k) \sum_{l=1}^{n_j} \cos(2\theta_l) + \sum_{k=1}^{n_i} \sin(2\theta_k) \sum_{l=1}^{n_j} \sin(2\theta_l) \\
&= \sum_{k=1}^{n_j} \sum_{l=1}^{n_j} \cos(2(\theta_k - \theta_l))
\end{aligned}
$$

Computing the expected value of the above expressions gives integrals of the form:

$$ \int_{0}^{2\pi} \int_{0}^{2\pi} \cos(2(\theta - \gamma)) \, d\theta d\gamma = 0, $$

and we simply obtain $E[c_i^2 + s_i^2] = n_i = 3 \cdot i$ and $E[c_i c_j + s_i s_j] = 0$. Therefore:

$$ 
\begin{aligned}
\sigma_{i/c} &\le \sqrt{\sum_{i=1}^{\infty} (\sigma_r^i)^2 n_i} \\
&\le 3\kappa \sqrt{\pi} \frac{sR}{d_c^2} \sqrt{\sum_{i=1}^{\infty} \frac{3}{i^3}} \\
&\le 3\kappa \sqrt{\pi s R \eta_c} \sqrt{3\zeta(3)}
\end{aligned}
$$

where $\zeta$ is the Riemann Zeta function. We can now use the GW contact model to replace $sR\eta_c = sR\eta e^{-h} = A_c/(\pi A_0)$:

$$ \sigma_{i/c} \le 3\kappa \sqrt{\frac{3\zeta(3)}{\pi}} \cdot \frac{A_c}{A_0}. $$

Note that only the $\sqrt{3\zeta(3)}$ term depends on the estimation from Jensen's inequality, so $\sigma_{i/c}$ is indeed linear with respect to the contact ratio.

## Appendix B. Nucleation probability

$P_{crack}$ as defined in eq. (12) is a conditional probability. It expresses the question “knowing an asperity is in contact, what is the probability that a crack nucleates at the contact edge?” The final expression for this probability is obtained by manipulating the inequality. We evaluate $\sigma_r$ at the contact radius $a = R \cdot s\sqrt{z^* - \bar{h}}$:

$$ \sigma_I^{\text{Hertz}} = \sigma_r(R \cdot s\sqrt{z^* - \bar{h}}) = \kappa(z^* - h)^{\frac{1}{2}}, $$
which is replaced in the inequality:

$$ 
\begin{aligned}
& \frac{\sigma_c}{\sigma_I^{\text{Hertz}} + \sigma_{i/c}} \le \omega(\tau_j) \\
\Leftrightarrow & \frac{\sigma_I^{\text{Hertz}} + \sigma_{i/c}}{\sigma_c} \ge \frac{1}{\omega(\tau_j)} \\
\Leftrightarrow & \sigma_I^{\text{Hertz}} \ge \frac{\sigma_c}{\omega(\tau_j)} - \sigma_{i/c} \\
\Leftrightarrow & z^* \ge \left( \frac{\sigma_c/\omega(\tau_j) - \sigma_{i/c}}{\kappa} \right)^2 + h.
\end{aligned}
$$

Let us call X the event corresponding to the above inequality. We have:

$$ 
\begin{aligned}
P_{\text{crack}} &= P(X \mid z^* \ge h) \\
&= \frac{P(X \text{ and } z^* \ge h)}{P(z^* \ge h)} \\
&= \frac{P(X)}{P(z^* \ge h)} \\
&= \exp\left(-\left(\frac{\sigma_c/\omega(\tau_j) - \sigma_{i/c}}{\kappa}\right)^2\right),
\end{aligned}
$$

since $z^*$ follows the canonical exponential distribution.

## References

Aghababaei, R., Brink, T., Molinari, J.F., 2018. Asperity-Level Origins of Transition from Mild to Severe Wear. Physical Review Letters 120, 186105. doi:10.1103/PhysRevLett.120.186105.

Aghababaei, R., Warner, D.H., Molinari, J.F., 2016. Critical length scale controls adhesive wear mechanisms. Nature Communications 7, 11816. doi:10.1038/ncomms11816.

Akchurin, A., Bosman, R., Lugt, P.M., 2016. A Stress-Criterion-Based Model for the Prediction of the Size of Wear Particles in Boundary Lubricated Contacts. Tribology Letters 64, 35. doi:10.1007/s11249-016-0772-x.

Almqvist, A., Sahlin, F., Larsson, R., Glavatskih, S., 2007. On the dry elasto-plastic contact of nominally flat surfaces. Tribology International 40, 574–579. doi:10.1016/j.triboint.2005.11.008.

Archard, J.F., 1953. Contact and Rubbing of Flat Surfaces. Journal of Applied Physics 24, 981–988. doi:10.1063/1.1721448.

Bemporad, A., Paggi, M., 2015. Optimization algorithms for the solution of the frictionless normal contact between rough surfaces. International Journal of Solids and Structures 69–70, 94–105. doi:10.1016/j.ijsolstr.2015.06.005.

Bowden, F.P., Tabor, D., 1939. The Area of Contact between Stationary and between Moving Surfaces. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 169, 391–413.

Brink, T., Molinari, J.F., 2019. Adhesive wear mechanisms in the presence of weak interfaces: Insights from an amorphous model system. Physical Review Materials 3, 053604. doi:10.1103/PhysRevMaterials.3.053604.

Bui, H.D., 1969. Étude de l'évolution de la frontière du domaine élastique avec l'écrouissage et relations de comportement élasto-plastique des métaux cubiques. Ph.D. thesis. Faculté des Sciences de l'Université de Paris. Paris.
Burwell, J.T., 1957. Survey of possible wear mechanisms. Wear 1, 119-141. doi:10.1016/0043-1648(57)90005-4.

Burwell, J.T., Strang, C.D., 1952. On the Empirical Law of Adhesive Wear. Journal of Applied Physics 23, 18-28. doi:10.1063/1.1701970.

Bush, A.W., Gibson, R.D., Thomas, T.R., 1975. The elastic contact of a rough surface. Wear 35, 87-111. doi:10.1016/0043-1648(75)90145-3.

Campañá, C., Müser, M.H., 2007. Contact mechanics of real vs. randomly rough surfaces: A Green's function molecular dynamics study. EPL (Europhysics Letters) 77, 38005. doi:10.1209/0295-5075/77/38005.

David, C., Gosselet, P., 2015. Équations aux dérivées partielles: cours et exercices corrigés. Dunod, Paris.

Drucker, D.C., Prager, W., 1952. Soil mechanics and plastic analysis or limit design. Quarterly of Applied Mathematics 10, 157-165. doi:10.1090/qam/48291.

Frérot, L., 2019. Supplementary codes and data to "Crack Nucleation in the Adhesive Wear of an Elastic-Plastic Half-Space". Zenodo. doi:10.5281/zenodo.3479539.

Frérot, L., Aghababaei, R., Molinari, J.F., 2018. A mechanistic understanding of the wear coefficient: From single to multiple asperities contact. Journal of the Mechanics and Physics of Solids 114, 172-184. doi:10.1016/j.jmps.2018.02.015.

Frérot, L., Anciaux, G., Rey, V., Pham-Ba, S., Molinari, J.F., 2019a. Tamaas, a high-performance library for periodic rough surface contact. Zenodo. doi:10.5281/zenodo.3479237.

Frérot, L., Bonnet, M., Molinari, J.F., Anciaux, G., 2019b. A Fourier-accelerated volume integral method for elastoplastic contact. Computer Methods in Applied Mechanics and Engineering 351, 951-976. doi:10.1016/j.cma.2019.04.006.

Ghaednia, H., Wang, X., Saha, S., Xu, Y., Sharma, A., Jackson, R.L., 2017. A Review of Elastic-Plastic Contact Mechanics. Applied Mechanics Reviews 69, 060804-060804-30. doi:10.1115/1.4038187.

Greenwood, J.A., Williamson, J.B.P., 1966. Contact of Nominally Flat Surfaces. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 295, 300-319. doi:10.1098/rspa.1966.0242.

Griffith, A.A., 1921. VI. The phenomena of rupture and flow in solids. Phil. Trans. R. Soc. Lond. A 221, 163-198. doi:10.1098/rsta.1921.0006.

Hyun, S., Pei, L., Molinari, J.F., Robbins, M.O., 2004. Finite-element analysis of contact between elastic self-affine surfaces. Physical Review E 70, 026117. doi:10.1103/PhysRevE.70.026117.

Jacq, C., Nélias, D., Lormand, G., Girodin, D., 2002. Development of a Three-Dimensional Semi-Analytical Elastic-Plastic Contact Code. Journal of Tribology 124, 653. doi:10.1115/1.1467920.

Johnson, K.L., 1985. Contact Mechanics. Cambridge University Press, Cambridge. doi:10.1017/CB09781139171731.

Krithivasan, V., Jackson, R.L., 2007. An analysis of three-dimensional elasto-plastic sinusoidal contact. Tribology Letters 27, 31-43. doi:10.1007/s11249-007-9200-6.

La Cruz, W., Martínez, J., Raydan, M., 2006. Spectral residual method without gradient information for solving large-scale nonlinear systems of equations. Mathematics of Computation 75, 1429-1448. doi:10.1090/S0025-5718-06-01840-0.

Li, T., Shi, J., Wang, S., Zio, E., Ma, Z., 2019. Mesoscale Numerical Modeling for Predicting Wear Debris Generation. Tribology Letters 67, 38. doi:10.1007/s11249-019-1150-2.

Majumdar, A., Bhushan, B., 1991. Fractal Model of Elastic-Plastic Contact Between Rough Surfaces. Journal of Tribology 113, 1-11. doi:10.1115/1.2920588.
Meng, H.C., Ludema, K.C., 1995. Wear models and predictive equations: Their form and content. Wear 181, 443-457. doi:10.1016/0043-1648(95)90158-2.

Menga, N., Carbone, G., 2019. The surface displacements of an elastic half-space subjected to uniform tangential tractions applied on a circular area. European Journal of Mechanics - A/Solids 73, 137-143. doi:10.1016/j.euromechsol.2018.07.011.

Milanese, E., Brink, T., Aghababaei, R., Molinari, J.F., 2019. Emergence of self-affine surfaces during adhesive wear. Nature Communications 10, 1116. doi:10.1038/s41467-019-09127-8.

Nayak, P.R., 1971. Random Process Model of Rough Surfaces. Journal of Lubrication Technology 93, 398-407. doi:10.1115/1.3451608.

Pastewka, L., Robbins, M.O., 2014. Contact between rough surfaces and a criterion for macroscopic adhesion. Proceedings of the National Academy of Sciences 111, 3298-3303. doi:10.1073/pnas.1320846111.

Pei, L., Hyun, S., Molinari, J.F., Robbins, M.O., 2005. Finite element modeling of elasto-plastic contact between rough surfaces. Journal of the Mechanics and Physics of Solids 53, 2385-2409. doi:10.1016/j.jmps.2005.06.008.

Persson, B.N.J., 2001a. Elastoplastic Contact between Randomly Rough Surfaces. Physical Review Letters 87, 116101. doi:10.1103/PhysRevLett.87.116101.

Persson, B.N.J., 2001b. Theory of rubber friction and contact mechanics. The Journal of Chemical Physics 115, 3840-3861. doi:10.1063/1.1388626.

Pham-Ba, S., Brink, T., Molinari, J.F., 2019. Adhesive wear and interaction of tangentially loaded micro-contacts. International Journal of Solids and Structures doi:10.1016/j.ijsolstr.2019.10.023.

Polonsky, I.A., Keer, L.M., 1999. A numerical method for solving rough contact problems based on the multi-level multi-summation and conjugate gradient techniques. Wear 231, 206-219. doi:10.1016/S0043-1648(99)00113-1.

Popov, V.L., Pohrt, R., 2018. Adhesive wear and particle emission: Numerical approach based on asperity-free formulation of Rabinowicz criterion. Friction doi:10.1007/s40544-018-0236-4.

Reddy, B.D., Martin, J.B., 1994. Internal Variable Formulations of Problems in Elastoplasticity: Constitutive and Algorithmic Aspects. Applied Mechanics Reviews 47, 429-456. doi:10.1115/1.3111086.

Sevostianov, I., Kachanov, M., 2012. Effective properties of heterogeneous materials: Proper application of the non-interaction and the "dilute limit" approximations. International Journal of Engineering Science 58, 124-128. doi:10.1016/j.ijengsci.2012.03.031.

Simo, J.C., Hughes, T.J.R., 1998. Computational Inelasticity. Number v. 7 in Interdisciplinary Applied Mathematics, Springer, New York.

Song, Z., Komvopoulos, K., 2013. Elastic-plastic spherical indentation: Deformation regimes, evolution of plasticity, and hardening effect. Mechanics of Materials 61, 91-100. doi:10.1016/j.mechmat.2013.01.003.

Stanley, H.M., Kato, T., 1997. An FFT-Based Method for Rough Surface Contact. Journal of Tribology 119, 481-485. doi:10.1115/1.2833523.

Tabor, D., 1951. The Hardness of Metals. Monographs on the Physics and Chemistry of Materials, Clarendon Press, Oxford.

Telles, J.C.F., Carrer, J.A.M., 1991. Implicit procedures for the solution of elastoplastic problems by the boundary element method. Mathematical and Computer Modelling 15, 303-311. doi:10.1016/0895-7177(91)90075-I.
Vakis, A.I., Yastrebov, V.A., Scheibert, J., Nicola, L., Dini, D., Minfray, C., Almqvist, A., Paggi, M., Lee, S., Limbert, G., Molinari, J.F., Anciaux, G., Aghababaei, R., Echeverri Restrepo, S., Papangelo, A., Cammarata, A., Nicolini, P., Putignano, C., Carbone, G., Stupkiewicz, S., Lengiewicz, J., Costagliola, G., Bosia, F., Guarino, R., Pugno, N.M., Müser, M.H., Ciavarella, M., 2018. Modeling and simulation in tribology across scales: An overview. Tribology International 125, 169–199. doi:10.1016/j.triboint. 2018.02.005.

Weber, B., Suhina, T., Junge, T., Pastewka, L., Brouwer, A.M., Bonn, D., 2018. Molecular probes reveal deviations from Amontons’ law in multi-asperity frictional contacts. Nature Communications 9, 888. doi:10.1038/s41467-018-02981-y.

Yastrebov, V.A., Anciaux, G., Molinari, J.F., 2015. From infinitesimal to full contact between rough surfaces: Evolution of the contact area. International Journal of Solids and Structures 52, 83–102. doi:10.1016/j.ijsolstr.2014.09.019.

Yastrebov, V.A., Anciaux, G., Molinari, J.F., 2017. The role of the roughness spectral breadth in elastic contact of rough surfaces. Journal of the Mechanics and Physics of Solids 107, 469–493. doi:10.1016/j. jmps.2017.07.016.

Zhang, F., Liu, J., Ding, X., Wang, R., 2019. Experimental and finite element analyses of contact behaviors between non-transparent rough surfaces. Journal of the Mechanics and Physics of Solids 126, 87–100. doi:10.1016/j.jmps.2019.02.004.