Elastic shakedown and roughness evolution in
repeated elastic-plastic contact

Lucas Frérot,* Lars Pastewka†

Department of Microsystems Engineering, University of Freiburg,
Georges-Köhler-Allee 103, 79110 Freiburg, Germany

December 4, 2023

Surface roughness emerges naturally during mechanical removal of material, fracture, chemical deposition, plastic deformation, indentation, and other processes. Here, we use continuum simulations to show how roughness which is neither Gaussian nor self-affine emerges from repeated elastic-plastic contact of rough and rigid surfaces on a flat elastic-plastic substrate. Roughness profiles change with each contact cycle, but appear to approach a steady-state long before the substrate stops deforming plastically and has hence “shaken-down” elastically. We propose a simple dynamic collapse for the emerging power-spectral density, which shows that the multi-scale nature of the roughness is encoded in the first few indentations. In contrast to macroscopic roughness parameters, roughness at small scales and the skewness of the height distribution of the resulting roughness do not show a steady-state, with the latter vanishing asymptotically with contact cycle.

**Keywords:** surface; roughness; contact; elastoplasticity; tribology; elastic shakedown

# 1 Introduction

Friction, wear, thermal transfer and other tribological phenomena depend on the roughness of surfaces in contact: Its presence means that contact only occurs on a limited subset of the apparent area of contact, the *true contact area* A, whose magnitude and geometry are of prime importance for tribology [1]. Understanding the origins of roughness, with its multiscale nature [2, 3], and making predictions on the true contact area are therefore objectives of the study of rough interfaces between solids. They are also both intimately linked to material properties.

Early models for contact of rough surfaces assumed that deformation is purely plastic [4–8]. Assuming a surface hardness $p_m$, this means that for a normal force $F_n$ the true area of contact is given by $A = F_n/p_m$, because the mean pressure within contacts equals $p_m$. However, researchers in the 1950s and 60s discussed the idea of an “elastic shakedown” [9]: It seemed unlikely, that a contact would continue to behave plastically during repeated loading. Eventually, the contact

should “shake down” to a point where it only responds elastically. This eventually lead to the development of elastic contact models, such as the famous 1966 model of Greenwood and Williamson [10]. The idea of elastic shakedown is probably one of the oldest ideas in the rough contact literature. Here, we present numerical calculations of the contact of rough surfaces with an elastic-plastic substrate to gain insights into the deformation processes during shakedown.

It seems reasonable to assume that a signature of shakedown is hidden in surface roughness: Self-affine surface roughness is typically associated with irreversible deformation, such as fracture [11, 12] or wear [2, 13–17], but plastic deformation can also contribute to the emergence of roughness [18–20]. In this work, we focus on the roughness created through quasi-static repeated indentation (i.e. without sliding), at a small contact area ratio, of an initially flat elastic-plastic half-space with a number of rough, rigid, periodic, isotropic counter-faces, each different. Our work stands apart from prior theoretical work on elastic-plastic contact [21, 22], which has investigated the con-
tact of a rough deformable half-space with a rigid but flat plane. While those calculations show smoothening of asperities, we find that the initially flat counterbody continuously roughens during our calculations. However, the statistical properties of the deformed body differ significantly from the indenting surface: Profiles are non-Gaussian and do not appear to be self-affine, even if the indenting surface has ideal self-affine properties. Even though simple macroscopic properties, such as the RMS height, appear to saturate at a steady-state value after a few contact cycles, our calculations reveal that surface topography and subsurface plastic deformation continues to evolve beyond this apparent “shake down” of the contact.

# 2 Methods

## 2.1 Elastoplastic Contact

Figure 1 illustrates the repeated indentation procedure described in the introduction. We simulate the contact between an elastic-plastic half-space B and a periodic rigid rough surface using a volume-integral equation approach [23–25]. In these repeated contact simulations, without horizontal movement, we denote $h_k$ the rigid surface of the k-th indentation. Unlike conventional boundary-integral equation methods used in elastic and elastic-plastic rough contacts [26–28, 21, 29], our approach fully considers plastic residual strains in the subsurface region of the half-space. The subsurface (volumetric) plastic strain field contributes to displacements and stress [30]. Equilibrium equations are solved with Green’s functions applied in the Fourier domain [25].

The plasticity model assumes an additive strain decomposition and a von Mises ($J_2$) yield criterion, $f_y$:

$$ \sigma = C : (\varepsilon - \varepsilon^p), \qquad (1) $$

$$ f_y(\sigma) = \sqrt{\frac{3}{2}} \|s\|, \quad \text{with } s = \sigma - \frac{1}{3}\text{tr}(\sigma)I. \qquad (2) $$

We have defined $\sigma$ the Cauchy stress tensor, $\varepsilon$ the total strain, $C$ the small-strain elasticity tensor, $\varepsilon^p$ the plastic strain and $s$ the deviatoric stress. The colon operator here introduces contraction on the two right-most indices of a tensor. This type of model, also referred to as $J_2$ plasticity [31], is the most widely-used canonical plasticity model. It assumes plastic deformation progresses as volume-conserving laminar flow beyond

a yield shear (von Mises) stress, but not under hydrostatic conditions. The model is empirical and material-unspecific, but captures the basic phenomenology of plastic flow at macroscopic scales.

Since our simulations are quasi-static, and the pseudo-time variable is given by the number of different surfaces used in indentation, $k$, we describe plastic strain rates as increments [31]. The equivalent cumulated plastic strain is defined as:

$$ e^p = \sqrt{\frac{2}{3}} \sum_{i=1}^{k} \|\Delta\varepsilon_i^p\|. \qquad (3) $$

The plasticity conditions can be written as, in terms of $e^p$ and $\sigma$:

$$ f_y(\sigma) - f_h(e^p) \le 0, \qquad (4) $$

$$ \Delta e^p \ge 0, \qquad (5) $$

$$ (f_y(\sigma) - f_h(e^p)) \Delta e^p = 0, \qquad (6) $$

where $f_h(e^p) = \sigma_0 + E_h e^p$ is the linear hardening function with initial yield stress $\sigma_0$ and hardening modulus $E_h$. These conditions express that the von Mises stress should never exceed the yield stress and that plastic deformation is irreversible. We assume an associated flow rule, so the plastic strain increment is expressed as:

$$ \Delta\varepsilon^p = \frac{3\Delta e^p}{2f_y(\sigma + C : \Delta\varepsilon)} s(\sigma + C : \Delta\varepsilon). \qquad (7) $$

The usual return-mapping algorithm [31] is used to compute $\Delta e^p$. Solving the coupled elastic-plastic contact problem is done with a fixed-point iterative approach [24, 25] accelerated with an Anderson mixing procedure [32, 33]. The code, available as a supplementary material, is built on top of the open-source high-performance contact library TAMAAS [34].

For reference, we also used the surface plasticity approach of Almqvist et al. [21] and a rigid-plastic, as known as a bearing area model. The former solves an elastic contact problem with the constraint that the normal surface pressure be less than or equal to the indentation hardness $p_m$ [35], correcting the counterface so that the constraint is satisfied. We use this correction as the residual plastic displacement [21, 29]. The choice of the hardness $p_m$ is a debated question: While $p_m \approx 3\sigma_y$ is a common choice stemming from spherical indentation [35], finite element simulations of sinusoidal [36, 37] and rough surfaces [38] report values in the order of $6\sigma_y$, and a quantitative agreement
<figure><img src="image_1.png" /><figcaption>Fig. 1:</figcaption></figure>

Two dimensional schematic view of the repeated indentation procedure we employ in this work. Indentation is made with a different surface at each iteration to reflect that in practice, contacts would be unlikely to occur twice with the same alignment of surfaces. B is a semi-infinite elastic-plastic solid obeying a $J_2$ isotropic linear hardening plasticity model. The rigid, rough, isotropic counterface has a period $L$.

with the $J_2$ model requires a surface-specific adjust-ment of $p_m$ [30]. Here we nonetheless set $p_m = 3\sigma_y$. The model is hereafter referred to as the saturated plasticity model. Finally, the rigid-plastic, or bearing area approach, consists in solving the contact entirely without elastic deformation, by finding a plane intersecting the rough counterface such that the intersection area equals the normal force, $F_n$, divided by $p_m$. Note that these two models corresponds roughly to historic plastic contact models by Holm [8].

## 2.2 Roughness Statistics

Gaussian, self-affine rough surfaces with a period $L$, are generated using a random-phase algorithm [39, 40, 3, 41] for each indentation step. They follow the isotropic power-spectral density (PSD)

$$|\mathcal{F}[h_k]|^2 = \begin{cases} C & \text{if } q_0 \le q \le q_1, \\ C \left(\frac{q}{q_1}\right)^{-2(H+1)} & \text{if } q_1 \le q \le q_2, \\ 0 & \text{otherwise,} \end{cases} \quad (8)$$

where $\mathcal{F}$ is the Fourier transform, $q_i$ are angular wavenumbers for the cutoffs in the surface spectrum, $H$ is the Hurst exponent, and $C$ an unspecified constant, adjusted such that the root-mean-square of surface slopes has the desired expected value.

The residual displacement resulting from plastic deformation at step $k$ is noted $u_k^P$ and is computed using a Green's function [42], which accounts for “spring-back” effects,

$$u_k^P(x) = \int_B \nabla_y U(x, y) : C : \epsilon_k^P dV_y, \quad (9)$$

where $U$ is the Mindlin tensor [43, 44, 25]. We note $h_k^P$ the trace on $\partial B$, the boundary of $B$, of the component of $u_k^P$ normal to the surface: It is the emerging surface roughness due to successive indentations.

Be it from plastic indentation or from synthetic generation, we derive from a rough surface $h$ the following statistical properties [41],

$$(h_{\text{rms}}^{(\alpha)})^2 = \frac{1}{L^2} \int_{A_L} |\nabla^{\alpha} h|^2 \, dA, \quad (10)$$

$$h_{\text{rms}} = h_{\text{rms}}^{(0)}, \quad h'_{\text{rms}} = h_{\text{rms}}^{(1)}, \quad h''_{\text{rms}} = \frac{1}{2} h_{\text{rms}}^{(2)}, \quad (11)$$

$$s = \frac{1}{h_{\text{rms}}^3 L^2} \int_{A_L} (h - \langle h \rangle)^3 \, dA, \quad (12)$$

where $A_L = [0, L]^2$, $\nabla^\alpha$ is the fractional Riesz derivative, computed in the Fourier domain, which links $h_{\text{rms}}^{(\alpha)}$ to Nayak's definition of PSD moments [45]. The parameters $h_{\text{rms}}$, $h'_{\text{rms}}$ and $h''_{\text{rms}}$ are the root-mean-square (RMS) of heights, RMS of slopes and RMS of curvatures, respectively. The parameter $s$ is the skewness of the surface height distribution.

Rather than reporting the curvature directly, we report it in form of Nayak's parameter [45],

$$\alpha = 4(h_{\text{rms}} h''_{\text{rms}} / h'^2_{\text{rms}})^2. \quad (13)$$
The advantage of using Nayak's parameter is that it removes the influence of changes in roughness amplitude.

Since rough surfaces are multi-scale objects, we use a multi-scale quantitative roughness analysis in addition to the scalar measures defined above. As one of the most common statistical descriptors, we report the isotropic PSD $\phi(q) = |F[h]|^2$. We also apply a variable bandwidth method (VBM) in form of detrended fluctuation analysis [46–48, 41]. Specifically, the surface with side length $L$ is divided into square patches of side length $\ell$ (see inset in fig. 7 for an illustration). Each patch is then independently detrended by subtracting the plane that minimizes $h_{\text{rms}}$ in that patch. Evaluating the above roughness parameters on the detrended patches for decreasing $\ell$ gives length-scale dependent metrics. Both PSD and VBM allow the study of a surface's potential self-affinity, with self-affine surfaces following $\phi(q) \propto q^{-2-2H}$ or $h_{\text{rms}}(\ell) \propto \ell^H$ over the self-affine bandwidth.

We also compute the dissipated plastic energy $d_k^p$ during successive indentations, given by

$$d_k^p = \sum_{i=1}^{k} \int_{B_L} (\sigma_i : \Delta\epsilon_i^p - E_h e_i^p \Delta e_i^p) dV, \quad (14)$$

for the $J_2$ model ($B_L = A_L \times [0, +\infty[$), and the irrecoverable energy, which includes the dissipated plastic energy and the stored elastic energy due to residual stresses, and can be computed directly on the surface for both models:

$$\delta_k^p = \sum_{i=1}^{k} \int_{A_L} p_i \Delta h_i^p dA. \quad (15)$$

Note that the saturated model lacks a description of eigenstresses [30], thus $\delta_k^p \equiv d_k^p$ is the dissipated plastic energy for this model.

## 3 Results

We simulate the repeated indentation of 90 different surfaces with the same statistical properties and the following PSD parameters: $q_0L/2\pi = 2$, $q_1L/2\pi = 2$, $q_2L/2\pi = 128$, $H \in \{0.5, 0.8\}$, $h'_{\text{rms}} = 0.4$, where $L$ is the horizontal period of the surfaces. Since the computational cost of the saturated model is significantly lower than the $J_2$ model, we simulate it in 270 repeated indentation steps. Contact is made against a half-space with material properties $v = 0.3$, $\sigma_0 = 0.1E$,

$E_h = 0.1E$, with $E$ the Young modulus. All properties are nondimensionalized by $E$ and by the horizontal period $L$ and only those parameters control the outcome of the simulations.

Figure 2 shows a rendered view of the emerging surface roughness due to the accumulation of plastic residual deformation, for a normal force of $F_n/(Eh'_{\text{rms}}L^2) = 0.025$ (which corresponds to a contact area ratio of 7%), and a discretization size of 729×729×55 points for a domain size of $L \times L \times 0.075L$ (total of 175M degrees-of-freedom). The highlighted areas in fig. 2 show where plastic energy is dissipated on the surface for the current indentation step. We also show select one-dimensional cross-sectional profiles of these topography maps below them. After enough steps, no remaining “undeformed” area of the surface can be identified by eye, and the surface appears entirely rough from plasticity.

We first evaluate macroscopic scalar roughness parameters on the emerging roughness. Figure 3 shows the evolution, for the indented roughness, of the root-mean-square of heights $h_{k,\text{rms}}^p$ (panel a), the RMS of slopes $h_{k,\text{rms}}^{p'}$ (b), the dissipated energy $d_k^p$ and the irrecoverable energy $\delta_k^p$ normalized by the elastic conforming energy of the counterface $w^e = [h_{\text{rms}}^{(1/2)}]^2 E/(4(1-\nu^2))$ (c), and Nayak's parameter, $\alpha^p$ (d). The dynamics of the first three parameters is similar in the two plasticity models, $J_2$ (von Mises), and saturated. Both RMS metrics initially scale as a power-law with an exponent close to 1/2 before plateauing at steady-state values, with the $J_2$ model transitioning earlier to lower values than the saturated model. Interestingly, the saturated model, which does not have hardening, also reaches a plateau. The overall scaling of RMS metrics is independent of $p_m$ (not shown), and $p_m$ only controls the contact area at each indentation: A lower contact ratio delays the appearance of a steady-state roughness. While energy dissipation in the saturated models appears to be constant per indentation step, the $J_2$ model shows lower energy dissipation rates after the first ~20 contact cycles. The irrecoverable energy shows a trend similar to the dissipated energy. It has a larger value, which indicates a build-up of residual stresses with successive indentations. Residual stresses require work for their formation, and while the stored elastic energy is in principle reversible, it is not recoverable in practice. The ratio $\delta_k^p/d_k^p$ increases from 1.22 to 1.34 with $k$, reflecting an increase in the proportion of elastic stored energy in $\delta_k^p$, consistent
<figure><img src="image_2.png" /><figcaption>Fig. 2:</figcaption></figure>

Evolution of surface roughness through repeated elastic-plastic indentation. Highlighted areas on the three-dimensional views show where plastic energy is dissipated at the specific step that is shown. The panels show the topography after (a) $k = 1$, (b) 10 and (c) 90 repeated contacts. The bottom row shows the evolution of one-dimensional cross-sectional profiles of the above surfaces. After sufficient number of iterations, no undeformed (flat) areas can be identified in the surface.

with hardening in the $J_2$ model.

As a geometric measure, Nayak's parameter evolves differently in the $J_2$ and the saturated plasticity model. In the saturated plasticity model, residual displacements are non-zero only in areas where the contact pressure $p = p_m$. This criterion creates a slope discontinuity at the edge of zones where $p = p_m$ which dominates the high-frequency part of the PSD. They are reflected in large values of $h''_{rms}$, contributing to the sustained increase in $\alpha_k^p$ in the saturated model.

The influence of the Hurst exponent $H$ disappears when we nondimensionalize lengths, slopes and energies by RMS heights, slopes and conforming elastic energy, respectively, of the counterface. However, normalization does not collapse the steady-state value of RMS heights and Nayak's parameter, indicating that details of the counterface self-affine structure non-trivially affect the long-wavelength part of the emerging roughness spectrum.

Having established the evolution of macroscopic parameters with contact cycle, we now detail the evolution of the distribution of heights. Figure 4 shows, for the $J_2$ (a), saturated (b), and rigid-plastic (c) models, the height probability density function of the emerging roughness through the indentation steps, for the counterface with $H = 0.8$.

While the distributions are always non-Gaussian, there is a clear evolution in the shape of the distribution with $k$: For low indentation count, the distribution obtained from the $J_2$ model is dominated by a well-

defined peak, similar to a rectified distribution, with an exponential tail. This peak occurs for a height which corresponds to the residual displacement in the zones that have not been in contact yet. The peak has a finite width due to the non-local influence of subsurface residual strains. As more surfaces come in contact with increasing $k$, the “undeformed” area of the surface shrinks: the peak height decreases and its breadth increases due to a larger variation of residual displacements, caused by ever-increasing, in number and size, plastic zones. The peak eventually disappears, and the height distribution settles on a more regular shape.

Unlike the $J_2$ model, the saturated model gives a true rectified distribution, where the peak in fig. 4a becomes a Dirac distribution in fig. 4b, whose magnitude is given by the remaining (truly) undeformed area. The tail of the distribution at large (positive) heights is similar for both models, although the saturated model has not reached a steady-state RMS of heights at the last iteration shown in fig. 4b. Similarly to the saturated, the rigid-plastic model (fig. 4c) also shows a rectified distribution, except that the residual height range is larger by a factor of approximately 2. This is due to the initial elastic deformation needed to overcome the plastic yield stress, absent in the rigid-plastic approach.

Figure 4 shows fits of the simulation data with two extreme value distributions: the Weibull<fnref n="1" /> distribution, which Silva Sabino et al. [49] used in an elastic contact

<footnote n="1">¹cumulative distribution cdf(x) = 1 - exp((-x/λ)⁵), λ = 0.3, κ = 1.9</footnote>
<figure><img src="image_3.png" /></figure>

**Fig. 3:**

Evolution of the normalized root-mean-square (RMS) of (a) heights, (b) slopes, (c) dissipated and irrecoverable energy, and (d) Nayak's parameter with indentation step. While, for the $J_2$ model, the RMS of heights and slopes show similar evolution, with a steady-state after ~ 40 contact cycles, Nayak's parameter converges earlier and the dissipated energy does not reach a steady-state value within 90 contacts. The scaling of the saturated results with $k$ is similar to the $J_2$ model, with a steady state reached after a larger number of iterations. Nayak's parameter (d) differs qualitatively between calculations with $J_2$ and saturated plasticity models due to slope discontinuities in the residual displacement.
<figure><img src="image_4.png" /><figcaption>Fig. 4:</figcaption></figure>

Evolution of the probability density function (PDF) of residual roughness heights for (a) the $J_2$, (b) the saturated plasticity model and (c) the rigid-plastic model. The peak corresponds to the residual deformation in untouched areas. As more areas come in contact, the peak broadens and eventually disappears. The height distribution then tends towards a non-Gaussian, asymmetric form, which resembles common extreme value distributions.

analysis of non-Gaussian surfaces, and the Gumbel distribution<fnref n="2" />. We use extreme value distributions because, owing to the relatively low contact area and to long-range elastic interactions—which makes contact near large asperities unlikely, only the largest asperities of the counterface come into contact, leading to a biased sampling of the counterface roughness by the deformable elastic-plastic manifold, which retains deformation from the roughness it “saw” in contact. Neither of the two distribution functions appears to be a perfect fit, but the Gumbel distribution captures the largest values better than Weibull. This clearly shows that the overall distribution has non-Gaussian tails. We note that the saturated plasticity model also has non-Gaussian tails, as it samples the peaks of the counterface. In contrast, the rigid-plastic (bearing area) model simply reflects the Gaussian distribution of the counterface.

From the prior results is it clear that the height distribution is not symmetric, which we now characterize from its skewness. Figure 5 shows, for the $J_2$ model only, that the skewness of the height distribution, $s_k^p$, decreases with $k$ as a power-law that is independent of the counterface’s Hurst exponent. Unlike the RMS of heights or slopes, there is no observable transition from a power-law to a constant regime. We obtain a power-law exponent of $z = -0.61(1)$ from a least-squares fit

<figure><img src="image_5.png" /><figcaption>Fig. 5:</figcaption></figure>

Evolution of the skewness of the height distribution with indentation step $k$ for the $J_2$ plasticity model. The skewness follows a power-law $k^z$, whose dynamics exponent $z$ is independent of the counterface’s Hurst exponent. This exponent describes the dynamics of the roughness evolution.

<footnote n="2">²cumulative distribution cdf($x$) = $\exp(-\exp(-x/\beta))$, $\beta$ = 1/9</footnote>
to the logarithm of the data. Although $-1 < z < 0$ indicates a slow evolution towards a symmetric height distribution, the skewness may converge to a finite value if the surface hardens enough to prevent further plastic dissipation, which would indicate elastic shake-down. Nonetheless, as we show next, the exponent $z$ is representative of the evolution dynamics of the roughness at all scales.

Figure 6 (a) shows, for each indentation step $k$ and for $H = 0.8$, the power-spectral density (PSD) of the emerging roughness scaled by $k^z$: individual PSDs then collapse onto a singular master curve above a wavenumber of $\sim 10$. The scaled PSD shows a constant regime up to that wavenumber, where it transitions to a regime with a Hurst exponent of $H^P \approx 0.8$, which corresponds to the counterface. Beyond the counterface's short wavelength cutoff (at wavenumber 128), the PSD scales as $q^{-8}$, which corresponds to $H^P = 3$. This scaling regime is due to the mechanical response at small scales: fig. 6b shows the PSDs of the last indentation for counterfaces with $H = 0.5$ and $H = 0.8$, and both have identical scaling beyond the short wavelength cutoff. The self-affine regime does not exist in the simulations we conducted in the case of a counterface with $H = 0.5$. The transition from a constant (white noise) to power-law (self-affine or elastic response) occurs on a range of scales too large to be captured within the counterface's spectral bandwidth ($q_2/q_1 = 64$). We discuss below whether both resulting surfaces are in fact self-affine in the regime where the counterface is self-affine.

Note that only “vertical” scaling of the PSD was required to obtain a collapse, indicating that the final surface’s scaling behavior is already encoded in the first indentations. This is corroborated by the behavior of Nayak’s parameter in fig. 3 which is independent of the PSD’s absolute magnitude, and shows little variation with $k$.

In fig. 7, we show a multi-scale VBM analysis— illustrated by the schematic inset in (a)—of two roughness parameters: the scale-dependent RMS of heights (a), $h_{k,\text{rms}}^P(\ell)$, and skewness (b), $s_k^P(\ell)$, where $\ell$ is the window size. An infinite self-affine surface shows a power-law behavior of $h_{\text{rms}}(\ell) \propto \ell^H$. Due to the restricted bandwidth of the spectrum, the counterface follows this scaling approximately (but not perfectly). In contrast to this, the emerging roughness is clearly not self-affine.

The skewnews, shown in fig. 7b, is normalized by $k^{-z}$. The collapse of skewness for all scales of curves

with different $k$ corroborates the observation that $z$ describes the roughness evolution dynamics at all scales. As magnification increases, the skewness decreases, i.e. the height distributions become more symmetric at smaller scales.

This behavior is directly shown in fig. 8, which depicts the height distribution at magnifications $L/\ell = 1$ and 128 for both values of the Hurst exponent at the last step $k = 90$. The plot shows that the Hurst exponent has no influence on the shape of the height distribution. At large magnification, corresponding to a small wavelength cutoff of the counterface, the height distribution is symmetric, although non-Gaussian, with approximately exponential tails. This is roughly consistent with the Gumbel tails observed at the length scale $L$ of the overall simulations.

# 4 Discussion

Our simulation show that some roughness properties, such as the RMS height or slope, reach a steady-state long before the surface actually “shakes downs”, in the sense that it stops deforming plastically and reacts only elastically and reversibly. Indeed, our simulations never reach a true shakedown, but our analysis of the dissipated energy indicates that the rate of plastic deformation decreases with contact cycle, but only for the $J_2$ model that explicitly considers hardening. Shakedown is never reached in the saturated plasticity model without hardening, even though the roughness may appear to have reached a state-state.

In fig. 3a, the RMS of heights for both models initially scales as $\sqrt{k}$. Despite memory (from plasticity), and long-range correlations (from elasticity), the RMS of heights behaves similarly to a simple ballistic deposition model at short times [50], indicating that early iterations are essentially uncorrelated. This can also be observed in the dissipated energy, as early indentations are unaffected by hardening, while later indentations occur on already hardened areas. However, since a perfect-plasticity hypothesis also yields a steady-state, we conclude that hardening only plays a minor role in the number of iterations required to observe a steady-state.

Another interesting aspect of our simulations is the emergence of non-Gaussian height distributions. Non-Gaussianity is manifested in skewness as well as exponential tails. The skewness has a well-defined monotonous dynamics described by a simple power-
<figure><img src="image_6.png" /><figcaption>**Fig. 6:**</figcaption></figure>

(a) Evolution of the power-spectral density (PSD) with indentation step for $H = 0.8$, and (b) final power-spectral density for both $H = 0.5$ and $H = 0.8$. Scaling the PSD with $k^z$ collapses all curves to a single master curve. The master curve is constant at low wavenumbers (large wavelengths) and has two distinct power-law regimes at intermediate and small wavelengths. The intermediate scaling regime follows self-affine scaling of the counterface only for $H = 0.8$. The topographies obtained after indentation with $H = 0.5$ and $H = 0.8$ show universal $q^{-8}$ scaling beyond the counterface's short wavelength cutoff.

<figure><img src="image_7.png" /><figcaption>**Fig. 7:**</figcaption></figure>

Variable bandwidth (VBM) analysis of (a) RMS heights and (b) skewness, for the simulations with a counterface with $H = 0.8$. The inset in the left plot illustrates the successive surface subdivision which the VBM employs.
<figure><img src="image_8.png" /><figcaption>**Fig. 8:**</figcaption></figure>

Height distributions after 90 repeated contacts for magnifications $L/\ell = 1$ and $L/\ell = 128$. The plot shows the results obtained for counterfaces with $H = 0.5$ and $H = 0.8$. Normalization by $h_{\text{rms}}^{\text{p}}(\ell)$ shows that the Hurst exponent does not affect the distribution shape. At small scales, the height distribution appears symmetric but non-Gaussian with exponential tails.

law $s \propto k^z$ with $-1 < z < 0$. This type of dynamics implies that the emerging rough geometry slowly becomes symmetric.

More insights are obtained by looking at the skewness as a function of scales. Figure 7b clearly shows that the small-scale roughness is more symmetric than the larger scales, and it appears that a factor of $k^{-z}$ collapses the dynamic evolution of the skewness even at small scales. This seems to indicate that during repeated indentation, the small-scale roughness diffuses towards larger scales, or in other words that the small-scales are representative of height-distribution at larger scales obtained at later indentation cycles. The height distributions hence become symmetric in the asymptotic steady-state.

However, the small-scale roughness does not follow a Gaussian distribution but has exponential tails (see fig. 8). The above argument seems to indicate, that even in the steady-state where height distributions are symmetric they are still not Gaussian. It is noteworthy that height distributions with exponential tails are typically assumed in independent-asperity models of rough contacts to make them analytically tractable [10, 51].

The height distribution has also been shown to have implications in the way the true contact area depends on the applied normal force [49], which could have consequences in contact sealing [52, 53] and wear particle formation [54–56], although some of these applications would also see the development of anisotropic roughness, due to sliding for example [13].

# 5 Summary & Conclusions

In summary, we have shown how surface roughness of an initially flat elastic-plastic substrate evolves during repeated contact with a rigid self-affine counterface. The resulting rough topography is neither self-affine nor Gaussian, even if the rigid body has those properties. Common scalar roughness measures, such as the RMS height, converge quickly to a steady-state value. Yet, subsurface plastic deformation still occurs long after the topography appears to have shaken-down geometrically. The topography continues to evolve, but this can only be seen at intermediate scales obtained from a scale-dependent geometric analyses, such as PSD.

The question of elastic shakedown has preoccupied tribologists for more than half a century [10, 57, 1]. Shakedown occurs during complex phenomena such as friction [58], fretting [59] and rolling [60]. Our simplified simulations are directly representative of processes such as sheet metal forming, where a die is brought into repeated contact with a counterface [61]. However, they certainly also hold insights for roughness evolution during sliding but clearly cannot explain the development of anisotropic topographies [13]. Our results contribute to the recent discussion on the emergence of surface roughness on solid but deformable bodies [14, 19, 16].

**Acknowledgments** The authors acknowledge the financial support of the Deutsche Forschungsgemeinschaft (DFG Grant No. 461911253, AWEARNESS). Computational resources from NEMO, at the University of Freiburg (DFG Grant No. INST 39/963-1 FUGG) are also acknowledged. The library TAMAAS [34] was used for contact simulation, roughness post-processing was done with TAMAAS and SURFACETOPOGRAPHY [62]. BLENDER and MATPLOTLIB [63] were used for visualization.
**Author Contributions** L.F. and L.P. designed the study, discussed results and edited the final manuscript. L.F. carried out calculations, analyzed the data and wrote the first manuscript draft.

**Data Availability** Simulation and figure codes are available on Zenodo [64].

**Competing Interests** The authors declare no competing interests.

## References

[1] Vakis, A. I. et al. “Modeling and Simulation in Tribology across Scales: An Overview”. In: *Tribology International* 125 (2018), pp. 169–199.
ISSN: 0301-679X. DOI: 10.1016/j.triboint.2018.02.005.

[2] Renard, F., Candela, T., and Bouchaud, E. “Constant Dimensionality of Fault Roughness from the Scale of Micro-Fractures to the Scale of Continents”. In: *Geophysical Research Letters* 40.1 (2013), pp. 83–87. ISSN: 1944-8007. DOI: 10.1029/2012GL054143.

[3] Persson, B. N. J., Albohr, O., Tartaglino, U., Volokitin, A. I., and Tosatti, E. “On the Nature of Surface Roughness with Application to Contact Mechanics, Sealing, Rubber Friction and Adhesion”. In: *Journal of Physics: Condensed Matter* 17.1 (2005), R1. ISSN: 0953-8984. DOI: 10.1088/0953-8984/17/1/R01.

[4] Binder, L. “Der Widerstand von Kontakten”. In: *Elektrotechnik und Maschinenbau* 30 (1912), pp. 781–788.

[5] Holm, R. “Eine Bestimmung der wirklichen Berührungsfläche eines Bürstenkontakes”. In: *Wissenschaftliche Veröffentlichungen aus den Siemens-Werken* 17.4 (1938), pp. 405–409.

[6] Holm, R. “Über die auf die wirkliche Berührungsfläche bezogene Reibkraft”. In: *Wissenschaftliche Veröffentlichungen aus den Siemens-Werken* 17.4 (1938), pp. 400–404.

[7] Holm, R. “Beitrag zur Kenntnis der Reibung”. In: *Wissenschaftliche Veröffentlichungen aus den Siemens-Werken* 20.1 (1941), pp. 68–84.

[8] Holm, R. *Electric Contacts: Theory and Applications*. 4th ed. Berlin; New York: Springer, 2000. 482 pp. ISBN: 978-3-540-03875-7.

[9] Johnson, K. L. “A Review of the Theory of Rolling Contact Stresses”. In: *Wear* 9.1 (1966), pp. 4–19. ISSN: 0043-1648. DOI: 10.1016/0043-1648(66)90010-X.

[10] Greenwood, J. A. and Williamson, J. B. P. “Contact of Nominally Flat Surfaces”. In: *Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences* 295.1442 (1966), pp. 300–319. ISSN: 1364-5021, 1471-2946. DOI: 10.1098/rspa.1966.0242.

[11] Mandelbrot, B. B., Passoja, D. E., and Paullay, A. J. “Fractal Character of Fracture Surfaces of Metals”. In: *Nature* 308.5961 (1984), pp. 721–722. ISSN: 0028-0836. DOI: 10.1038/308721a0.

[12] Ponson, L. “Statistical Aspects in Crack Growth Phenomena: How the Fluctuations Reveal the Failure Mechanisms”. In: *International Journal of Fracture* 201.1 (2016), pp. 11–27. ISSN: 1573-2673. DOI: 10.1007/s10704-016-0117-7.

[13] Candela, T. and Brodsky, E. E. “The Minimum Scale of Grooving on Faults”. In: *Geology* 44.8 (2016), pp. 603–606. ISSN: 0091-7613, 1943-2682. DOI: 10.1130/G37934.1.

[14] Milanese, E., Brink, T., Aghababaei, R., and Molinari, J.-F. “Emergence of Self-Affine Surfaces during Adhesive Wear”. In: *Nature Communications* 10.1 (2019), p. 1116. ISSN: 2041-1723. DOI: 10.1038/s41467-019-09127-8.

[15] Pham-Ba, S. and Molinari, J.-F. “Creation and Evolution of Roughness on Silica under Unlubricated Wear”. In: *Wear* 472–473 (2021), p. 203648. ISSN: 0043-1648. DOI: 10.1016/j.wear.2021.203648.

[16] Aghababaei, R., Brodsky, E. E., Molinari, J.-F., and Chandrasekar, S. “How Roughness Emerges on Natural and Engineered Surfaces”. In: *MRS Bulletin* 47.12 (2022), pp. 1229–1236. ISSN: 1938-1425. DOI: 10.1557/s43577-022-00469-1.
[17] Garcia-Suarez, J., Brink, T., and Molinari, J.-F. *Roughness Evolution Induced by Third-Body Wear*. 2023. DOI: 10.48550/arXiv.2306.08993. arXiv: 2306.08993 [cond-mat]. (Visited on 07/21/2023). preprint.

[18] Thomson, P. F. and Nayak, P. U. "The Effect of Plastic Deformation on the Roughening of Free Surfaces of Sheet Metal". In: *International Journal of Machine Tool Design and Research* 20.1 (1980), pp. 73–86. ISSN: 0020-7357. DOI: 10.1016/0020-7357(80)90020-7.

[19] Hinkle, A. R., Nöhring, W. G., Leute, R., Junge, T., and Pastewka, L. "The Emergence of Small-Scale Self-Affine Surface Roughness from Deformation". In: *Science Advances* 6.7 (2020), eaax0847. DOI: 10.1126/sciadv.aax0847.

[20] Nöhring, W. G., Hinkle, A. R., and Pastewka, L. "Nonequilibrium Plastic Roughening of Metallic Glasses Yields Self-Affine Topographies with Strain-Rate and Temperature-Dependent Scaling Exponents". In: *Physical Review Materials* 6.7 (2022), p. 075603. DOI: 10.1103/PhysRevMaterials.6.075603.

[21] Almqvist, A., Sahlin, F., Larsson, R., and Glavatskih, S. "On the Dry Elasto-Plastic Contact of Nominally Flat Surfaces". In: *Tribology International* 40.4 (2007), pp. 574–579. ISSN: 0301679X. DOI: 10.1016/j.triboint.2005.11.008.

[22] Tiwari, A., Wang, A., Müser, M. H., and Persson, B. N. J. "Contact Mechanics for Solids with Randomly Rough Surfaces and Plasticity". In: *Lubricants* 7.10 (10 2019), p. 90. DOI: 10.3390/lubricants7100090.

[23] Bonnet, M. and Mukherjee, S. "Implicit BEM Formulations for Usual and Sensitivity Problems in Elasto-Plasticity Using the Consistent Tangent Operator Concept". In: *International Journal of Solids and Structures* 33.30 (1996), pp. 4461–4480. ISSN: 0020-7683. DOI: 10.1016/0020-7683(95)00279-0.

[24] Jacq, C., Nélias, D., Lormand, G., and Girodin, D. "Development of a Three-Dimensional Semi-Analytical Elastic-Plastic Contact Code". In: *Journal of Tribology* 124.4 (2002), p. 653. ISSN: 07424787. DOI: 10.1115/1.1467920.

[25] Frérot, L., Bonnet, M., Molinari, J.-F., and Anciaux, G. "A Fourier-accelerated Volume Integral Method for Elastoplastic Contact". In: *Computer Methods in Applied Mechanics and Engineering* 351 (2019), pp. 951–976. ISSN: 0045-7825. DOI: 10.1016/j.cma.2019.04.006.

[26] Stanley, H. M. and Kato, T. "An FFT-Based Method for Rough Surface Contact". In: *Journal of Tribology* 119.3 (1997), pp. 481–485. ISSN: 0742-4787. DOI: 10.1115/1.2833523.

[27] Campañá, C., Müser, M. H., and Robbins, M. O. "Elastic Contact between Self-Affine Surfaces: Comparison of Numerical Stress and Contact Correlation Functions with Analytic Predictions". In: *Journal of Physics: Condensed Matter* 20.35 (2008), p. 354013. ISSN: 0953-8984. DOI: 10.1088/0953-8984/20/35/354013.

[28] Yastrebov, V. A., Anciaux, G., and Molinari, J.-F. "Contact between Representative Rough Surfaces". In: *Physical Review E* 86.3 (2012). ISSN: 1539-3755, 1550-2376.

[29] Weber, B., Suhina, T., Junge, T., Pastewka, L., Brouwer, A. M., and Bonn, D. "Molecular Probes Reveal Deviations from Amontons' Law in Multi-Asperity Frictional Contacts". In: *Nature Communications* 9.1 (2018), p. 888. ISSN: 2041-1723. DOI: 10.1038/s41467-018-02981-y.

[30] Frérot, L., Anciaux, G., and Molinari, J.-F. "Crack Nucleation in the Adhesive Wear of an Elastic-Plastic Half-Space". In: *Journal of the Mechanics and Physics of Solids* 145 (2020), p. 104100. ISSN: 0022-5096. DOI: 10.1016/j.jmps.2020.104100. arXiv: 1910.05163.

[31] Simo, J. C. and Hughes, T. J. R. *Computational Inelasticity*. Interdisciplinary Applied Mathematics v. 7. New York: Springer, 1998. 392 pp. ISBN: 978-0-387-97520-7.

[32] Anderson, D. G. "Iterative Procedures for Nonlinear Integral Equations". In: *Journal of the ACM* 12.4 (1965), pp. 547–560. ISSN: 0004-5411, 1557-735X. DOI: 10.1145/321296.321305.

[33] Eyert, V. "A Comparative Study on Methods for Convergence Acceleration of Iterative Vector Sequences". In: *Journal of Computational*
*Physics* 124.2 (1996), pp. 271–285. ISSN: 0021-9991. DOI: 10.1006/jcph.1996.0059.

[34] Frérot, L., Anciaux, G., Rey, V., Pham-Ba, S., and Molinari, J.-F. “Tamaas: A Library for Elastic-Plastic Contact of Periodic Rough Surfaces”. In: *Journal of Open Source Software* 5.51 (2020), p. 2121. ISSN: 2475-9066. DOI: 10.21105/joss.02121.

[35] Tabor, D. *The Hardness of Metals*. Monographs on the Physics and Chemistry of Materials. Oxford: Clarendon Press, 1951.

[36] Gao, Y. F., Bower, A. F., Kim, K. -., Lev, L., and Cheng, Y. T. “The Behavior of an Elastic–Perfectly Plastic Sinusoidal Surface under Contact Loading”. In: *Wear* 261.2 (2006), pp. 145–154. ISSN: 0043-1648. DOI: 10.1016/j.wear.2005.09.016.

[37] Ghaednia, H., Wang, X., Saha, S., Xu, Y., Sharma, A., and Jackson, R. L. “A Review of Elastic–Plastic Contact Mechanics”. In: *Applied Mechanics Reviews* 69.6 (2017), pp. 060804-060804–30. ISSN: 0003-6900. DOI: 10.1115/1.4038187.

[38] Pei, L., Hyun, S., Molinari, J.-F., and Robbins, M. O. “Finite Element Modeling of Elastoplastic Contact between Rough Surfaces”. In: *Journal of the Mechanics and Physics of Solids* 53.11 (2005), pp. 2385–2409. ISSN: 0022-5096. DOI: 10.1016/j.jmps.2005.06.008.

[39] Wu, J.-J. “Simulation of Rough Surfaces with FFT”. In: *Tribology International* 33.1 (2000), pp. 47–58. ISSN: 0301-679X. DOI: 10.1016/S0301-679X(00)00016-5.

[40] Ramisetti, S. B., Campañá, C., Anciaux, G., Molinari, J.-F., Müser, M. H., and Robbins, M. O. “The Autocorrelation Function for Island Areas on Self-Affine Surfaces”. In: *Journal of Physics: Condensed Matter* 23.21 (2011), p. 215004. ISSN: 0953-8984. DOI: 10.1088/0953-8984/23/21/215004.

[41] Jacobs, T. D. B., Junge, T., and Pastewka, L. “Quantitative Characterization of Surface Topography Using Spectral Analysis”. In: *Surface Topography: Metrology and Properties* 5.1 (2017), p. 013001. ISSN: 2051-672X. DOI: 10.1088/2051-672X/aa51f8.

[42] Bui, H. D. “Some Remarks about the Formula-
tion of Three-Dimensional Thermoelastoplastic
Problems by Integral Equations”. In: *International Journal of Solids and Structures* 14.11
(1978), pp. 935–939. ISSN: 0020-7683. DOI: 10.1016/0020-7683(78)90069-0.

[43] Mindlin, R. D. “Force at a Point in the Interior of a Semi-Infinite Solid”. In: *Journal of Applied Physics* 7.5 (1936), pp. 195–202. ISSN: 0021-8979, 1089-7550. DOI: 10.1063/1.1745385.

[44] Frérot, L. *The Mindlin Fundamental Solution - A Fourier Approach*. Zenodo, 2018. DOI: 10.5281/zenodo.1492149.

[45] Nayak, P. R. “Random Process Model of Rough Surfaces”. In: *Journal of Lubrication Technology* 93.3 (1971), pp. 398–407. ISSN: 0742-4787. DOI: 10.1115/1.3451608.

[46] Cannon, M. J., Percival, D. B., Caccia, D. C., Raymond, G. M., and Bassingthwaighte, J. B. “Evaluating Scaled Windowed Variance Methods for Estimating the Hurst Coefficient of Time Series”. In: *Physica A: Statistical Mechanics and its Applications* 241.3 (1997), pp. 606–626. ISSN: 0378-4371. DOI: 10.1016/S0378-4371(97)00252-5.

[47] Peng, C.-K., Buldyrev, S. V., Havlin, S., Simmons, M., Stanley, H. E., and Goldberger, A. L. “Mosaic Organization of DNA Nucleotides”. In: *Physical Review E* 49.2 (1994), pp. 1685–1689. DOI: 10.1103/PhysRevE.49.1685.

[48] Peng, C.-K., Havlin, S., Stanley, H. E., and Goldberger, A. L. “Quantification of Scaling Exponents and Crossover Phenomena in Nonstationary Heartbeat Time Series”. In: *Chaos: An International Journal of Nonlinear Science* 5.1 (1995), pp. 82–87. ISSN: 1054-1500. DOI: 10.1063/1.166141.

[49] Silva Sabino, T., Couto Carneiro, A. M., Pinto Carvalho, R., and Andrade Pires, F. M. “Evolution of the Real Contact Area of Self-Affine Non-Gaussian Surfaces”. In: *International Journal of Solids and Structures* 268 (2023), p. 112173. ISSN: 0020-7683. DOI: 10.1016/j.ijsolstr.2023.112173.
[50] Barabási, A. L. and Stanley, H. E. *Fractal Concepts in Surface Growth*. 1st ed. Cambridge University Press, 1995. ISBN: 978-0-521-48308-7 978-0-521-48318-6 978-0-511-59979-8. DOI: 10.1017/CB09780511599798.

[51] Persson, B. N. J. *Sliding Friction*. Red. by Von Klitzing, K. and Wiesendanger, R. NanoScience and Technology. Berlin, Heidelberg: Springer Berlin Heidelberg, 2000. ISBN: 978-3-642-08652-6 978-3-662-04283-0. DOI: 10.1007/978-3-662-04283-0.

[52] Dapp, W. B., Lücke, A., Persson, B. N. J., and Müser, M. H. "Self-Affine Elastic Contacts: Percolation and Leakage". In: *Physical Review Letters* 108.24 (2012), p. 244301. DOI: 10.1103/PhysRevLett.108.244301.

[53] Shvarts, A. G. and Yastrebov, V. A. "Trapped Fluid in Contact Interface". In: *Journal of the Mechanics and Physics of Solids* 119 (2018), pp. 140-162. ISSN: 0022-5096. DOI: 10.1016/j.jmps.2018.06.016.

[54] Frérot, L., Aghababaei, R., and Molinari, J.-F. "A Mechanistic Understanding of the Wear Coefficient: From Single to Multiple Asperities Contact". In: *Journal of the Mechanics and Physics of Solids* 114 (2018), pp. 172-184. ISSN: 0022-5096. DOI: 10.1016/j.jmps.2018.02.015.

[55] Popov, V. L. and Pohrt, R. "Adhesive Wear and Particle Emission: Numerical Approach Based on Asperity-Free Formulation of Rabinowicz Criterion". In: *Friction* (2018). ISSN: 2223-7704. DOI: 10.1007/s40544-018-0236-4.

[56] Brink, T., Frérot, L., and Molinari, J.-F. "A Parameter-Free Mechanistic Model of the Adhesive Wear Process of Rough Surfaces in Sliding Contact". In: *Journal of the Mechanics and Physics of Solids* 147 (2021), p. 104238. ISSN: 0022-5096. DOI: 10.1016/j.jmps.2020.104238. arXiv: 2004.00559.

[57] Kapoor, A., Williams, J. A., and Johnson, K. L. "The Steady State Sliding of Rough Surfaces". In: *Wear* 175.1 (1994), pp. 81-92. ISSN: 0043-1648. DOI: 10.1016/0043-1648(94)90171-6.

[58] Flicek, R. C., Hills, D. A., Barber, J. R., and Dini, D. "Determination of the Shakedown Limit for Large, Discrete Frictional Systems". In: *European Journal of Mechanics - A/Solids* 49 (2015), pp. 242-250. ISSN: 0997-7538. DOI: 10.1016/j.euromechsol.2014.08.001.

[59] Fouvry, S., Kapsa, Ph., and Vincent, L. "An Elastic-Plastic Shakedown Analysis of Fretting Wear". In: *Wear* 247.1 (2001), pp. 41-54. ISSN: 0043-1648. DOI: 10.1016/S0043-1648(00)00508-1.

[60] Berthe, L., Sainsot, P., Lubrecht, A. A., and Baietto, M. C. "Plastic Deformation of Rough Rolling Contact: An Experimental and Numerical Investigation". In: *Wear* 312.1-2 (2014), pp. 51-57. ISSN: 0043-1648. DOI: 10.1016/j.wear.2014.01.017.

[61] Lo, S.-W. and Horng, T.-C. "Surface Roughening and Contact Behavior in Forming of Aluminum Sheet". In: *Journal of Tribology* 121.2 (1999), pp. 224-233. ISSN: 0742-4787. DOI: 10.1115/1.2833925.

[62] Röttger, M. C., Sanner, A., Thimons, L. A., Junge, T., Gujrati, A., Monti, J. M., Nöhring, W. G., Jacobs, T. D. B., and Pastewka, L. "ContactEngineering—Create, Analyze and Publish Digital Surface Twins from Topography Measurements across Many Scales". In: *Surface Topography: Metrology and Properties* 10.3 (2022), p. 035032. ISSN: 2051-672X. DOI: 10.1088/2051-672X/ac860a.

[63] Hunter, J. D. "Matplotlib: A 2D Graphics Environment". In: *Computing in Science & Engineering* 9.3 (2007), pp. 90-95. ISSN: 1558-366X. DOI: 10.1109/MCSE.2007.55.

[64] Frérot, L. and Pastewka, L. *Supplementary Codes and Data to "Elastic Shakedown and Roughness Evolution in Repeated Elastic-Plastic Contact"*. Zenodo, 2023. DOI: 10.5281/zenodo.8280362.
<figure><img src="image_9.png" /><figcaption>**Fig. 9: Graphical Abstract**</figcaption></figure>