Elastic shakedown and roughness evolution in repeated elastic-plastic contact

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.


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 , 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][5][6][7][8].Assuming a surface hardness  m , this means that for a normal force  n the true area of contact is given by  =  n ∕ m , because the mean pressure within contacts equals  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][14][15][16][17], but plastic deformation can also contribute to the emergence of roughness [18][19][20].In this work, we focus on the roughness created through quasistatic repeated indentation (i.e.without sliding), at a small contact area ratio, of an initially flat elasticplastic half-space with a number of rough, rigid, periodic, isotropic counter-faces, each different.Our work stands apart from prior theoretical work on elasticplastic 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 steadystate 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.

Elastoplastic Contact
Figure 1 illustrates the repeated indentation procedure described in the introduction.We simulate the contact between an elastic-plastic half-space  and a periodic rigid rough surface using a volume-integral equation approach [23][24][25].In these repeated contact simulations, without horizontal movement, we denote ℎ  the rigid surface of the -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 ( 2 ) yield criterion,  y : We have defined  the Cauchy stress tensor,  the total strain,  the small-strain elasticity tensor,   the plastic strain and  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  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 materialunspecific, 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, , we describe plastic strain rates as increments [31].The equivalent cumulated plastic strain is defined as: The plasticity conditions can be written as, in terms of  p and : where  h ( p ) =  0 +  h  p is the linear hardening function with initial yield stress  0 and hardening modulus  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: The usual return-mapping algorithm [31] is used to compute Δ 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  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  m is a debated question: While  m ≈ 3 y is a common choice stemming from spherical indentation [35], finite element simulations of sinusoidal [36,37] and rough surfaces [38]  with the  2 model requires a surface-specific adjustment of   [30].Here we nonetheless set  m = 3 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,  n , divided by  m .Note that these two models corresponds roughly to historic plastic contact models by Holm [8].

Roughness Statistics
Gaussian, self-affine rough surfaces with a period , are generated using a random-phase algorithm [39,40,3,41] for each indentation step.They follow the isotropic power-spectral density (PSD) where  is the Fourier transform,   are angular wavenumbers for the cutoffs in the surface spectrum,  is the Hurst exponent, and  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  is noted  p  and is computed using a Green's function [42], which accounts for "spring-back" effects, where  is the Mindlin tensor [43,44,25].We note ℎ p  the trace on , the boundary of , of the component of  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 ℎ the following statistical properties [41], where   = [0, ] 2 , ∇  is the fractional Riesz derivative, computed in the Fourier domain, which links ℎ ()   rms to Nayak's definition of PSD moments [45].The parameters ℎ rms , ℎ ′ rms and ℎ ′′ rms are the root-mean-square (RMS) of heights, RMS of slopes and RMS of curvatures, respectively.The parameter  is the skewness of the surface height distribution.
Rather than reporting the curvature directly, we report it in form of Nayak's parameter [45], 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 () = | [ℎ]| 2 .We also apply a variable bandwidth method (VBM) in form of detrended fluctuation analysis [46][47][48]41].Specifically, the surface with side length  is divided into square patches of side length  (see inset in fig.7 for an illustration).Each patch is then independently detrended by subtracting the plane that minimizes ℎ rms in that patch.Evaluating the above roughness parameters on the detrended patches for decreasing  gives length-scale dependent metrics.Both PSD and VBM allow the study of a surface's potential self-affinity, with self-affine surfaces following () ∝  −2−2 or ℎ rms () ∝   over the self-affine bandwidth.
We also compute the dissipated plastic energy  p  during successive indentations, given by for the  2 model (  =   ×[0, +∞[), 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: Note that the saturated model lacks a description of eigenstresses [30], thus  p  ≡  p  is the dissipated plastic energy for this model.

Results
We simulate the repeated indentation of 90 different surfaces with the same statistical properties and the following PSD parameters: , where  is the horizontal period of the surfaces.Since the computational cost of the saturated model is significantly lower than the  2 model, we simulate it in 270 repeated indentation steps.Contact is made against a half-space with material properties  = 0.3,  0 = 0.1,  ℎ = 0.1, with  the Young modulus.All properties are nondimensionalized by  and by the horizontal period  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  n ∕(ℎ ′ rms  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 ××0.075(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 ℎ The dynamics of the first three parameters is similar in the two plasticity models,  2 (von Mises), and saturated.Both RMS metrics initially scale as a powerlaw with an exponent close to 1∕2 before plateauing at steady-state values, with the  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  m (not shown), and  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  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  with hardening in the  2 model.As a geometric measure, Nayak's parameter evolves differently in the  2 and the saturated plasticity model.In the saturated plasticity model, residual displacements are non-zero only in areas where the contact pressure  =  m .This criterion creates a slope discontinuity at the edge of zones where  =  m which dominates the high-frequency part of the PSD.They are reflected in large values of ℎ ′′ rms , contributing to the sustained increase in  p  in the saturated model.The influence of the Hurst exponent  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 nontrivially 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  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  = 0.8.
While the distributions are always non-Gaussian, there is a clear evolution in the shape of the distribution with : For low indentation count, the distribution obtained from the  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 , 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  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 Weibull1 distribution, which Silva Sabino et al. [49]   analysis of non-Gaussian surfaces, and the Gumbel distribution 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  2 model only, that the skewness of the height distribution,  p  , decreases with  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 powerlaw exponent of  = −0.61(1)from a least-squares fit Evolution of the skewness of the height distribution with indentation step  for the  2 plasticity model.The skewness follows a power-law   , whose dynamics exponent  is independent of the counterface's Hurst exponent.This exponent describes the dynamics of the roughness evolution.
to the logarithm of the data.Although −1 <  < 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 shakedown.Nonetheless, as we show next, the exponent  is representative of the evolution dynamics of the roughness at all scales.Figure 6 (a) shows, for each indentation step  and for  = 0.8, the power-spectral density (PSD) of the emerging roughness scaled by   : individual PSDs then collapse onto a singular master curve above a wavenumber of ∼ 10.The scaled PSD shows a constant regime up to that wavenumber, where it transitions to a regime with a Hurst exponent of  p ≈ 0.8, which corresponds to the counterface.Beyond the counterface's short wavelength cutoff (at wavenumber 128), the PSD scales as  −8 , which corresponds to  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  = 0.5 and  = 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  = 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 ( 2 ∕ 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 .
In fig. 7, we show a multi-scale VBM analysisillustrated by the schematic inset in (a)-of two roughness parameters: the scale-dependent RMS of heights (a), ℎ p ,rms (), and skewness (b),  p  (), where  is the window size.An infinite self-affine surface shows a power-law behavior of ℎ rms () ∝   .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  − .The collapse of skewness for all scales of curves with different  corroborates the observation that  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 ∕ = 1 and 128 for both values of the Hurst exponent at the last step  = 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  of the overall simulations.

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  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 √ .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 steadystate.
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-    (𝓁) 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  ∝   with −1 <  < 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  − 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][55][56], although some of these applications would also see the development of anisotropic roughness, due to sliding for example [13].

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 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, AWEAR-NESS).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 MAT-PLOTLIB [63] were used for visualization.

Self-affine rough surface Repeated indentations
Elastoplastic substrate Steady-state roughness appears before "elastic shakedown"

Fig. 2 :
Fig. 2: Evolution of surface roughness through repeated elastic-plastic indentation.Highlighted areas on the threedimensional views show where plastic energy is dissipated at the specific step that is shown.The panels show the topography after (a)  = 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.

Fig. 4 :
Fig. 4:Evolution of the probability density function (PDF) of residual roughness heights for (a) the  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.

2
Fig. 5:Evolution of the skewness of the height distribution with indentation step  for the  2 plasticity model.The skewness follows a power-law   , whose dynamics exponent  is independent of the counterface's Hurst exponent.This exponent describes the dynamics of the roughness evolution.

Fig. 6 :
Fig. 6: (a) Evolution of the power-spectral density (PSD) with indentation step for  = 0.8, and (b) final powerspectral density for both  = 0.5 and  = 0.8.Scaling the PSD with   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  = 0.8.The topographies obtained after indentation with  = 0.5 and  = 0.8 show universal  −8 scaling beyond the counterface's short wavelength cutoff.

Fig. 7 :
Fig. 7:Variable bandwidth (VBM) analysis of (a) RMS heights and (b) skewness, for the simulations with a counterface with  = 0.8.The inset in the left plot illustrates the successive surface subdivision which the VBM employs.

Fig. 9 :
Fig. 9: Graphical Abstract 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  is similar to the  2 model, with a steady state reached after a larger number of iterations.Nayak's parameter (d) differs qualitatively between calculations with  2 and saturated plasticity models due to slope discontinuities in the residual displacement.
k / d 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  2 model, the RMS of heights and slopes show similar evolution, with a steady-state after ∼ 40 contact cycles, Height distributions after 90 repeated contacts for magnifications ∕ = 1 and ∕ = 128.The plot shows the results obtained for counterfaces with  = 0.5 and  = 0.8.Normalization by ℎ p rms