On the Role of Roughness in the Indentation of Viscoelastic Solids

A numerical boundary element methodology is employed to understand how fractality intervenes when a 1D rigid rough profile indents a linear viscoelastic half-plane. The focus is, in particular, on the viscoelastic dissipation and how this is influenced by the profile statistical parameters, namely, the mean square roughness hrms\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{\rm rms}$$\end{document} of the profile, the mean square slope m2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_{2}$$\end{document} and the Hurst exponent H. Our numerical investigation, properly supported by a dimensional analysis, reveals that, in the one-dimensional case under investigation, the leading role is played by hrms\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{\rm rms}$$\end{document} and, thus, mainly by the large scales of the rough spectrum. Clearly, on an experimental level, this implies that a simple measure of the roughness parameter hrms\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{\rm rms}$$\end{document} is sufficient to determine the viscoelastic dissipation.


Introduction
Surface topography governs a variety of phenomena on several orders of magnitude, from the macro to the nano-scales, and, thus, constitutes the focus of most of the research efforts currently carried out in Tribology. Indeed, it is nowadays fundamental to understand how roughness is related to aspects like friction [1,2], wear [3][4][5], hysteretic dissipation [6][7][8]48], fluid percolation [9][10][11] and lubrication [12,13]. A clear comprehension of these phenomena is needed not only due to reasons related to pure theoretical interest, but is also driven by a very ambitious goal: current research trends in surface engineering aim, in fact, at shaping the surface to control frictional, wear and adhesive properties [14,15]. A field, where these research efforts have reached significant advancements, includes surely the laser microtextured surfaces: under lubricated conditions, a pattern of micro-dimples has proved to be effective in reducing friction and increasing the loading capacity [16][17][18]. On the other hand, mushrooms at the micro-scales have been successfully employed to enhance adhesion [19,20]. These are just examples of the groundbreaking potential implied in understanding and controlling roughness, but yet explain all the attention paid to surface tribology. In this context, theoretical and numerical models play a major role as they constitute a virtual laboratory to obtain rapid information on the phenomena previously discussed.
Hence, it is straightforward to understand the particular efforts made, in the last decades, to develop multi-scale simulations approaches which, at the same time, account for macro-features and micro-and nano-roughness. This is really crucial when dealing with materials with a timedependent rheology: in this case, morphology and materials properties combine to determine the solution in terms of stresses, strains and, ultimately, friction. The problem was pioneered in the 1960s by Lee and Radok [21] and, later, by Hunter [22] to study the contact mechanics of smooth solids on viscoelastic layers. Only in the last two decades, Persson has proposed an analytical theory which takes into account, with a good level of approximation, the role of roughness in this kind of contacts [1,2]. Parallelly, from a numerical point of view, a variety of methodologies has been developed. These include Finite Element Methods (FEM) [23]: in spite of their generality, however, FE approaches imply a really high computational load that is likely to become unfeasible when considering the entire rough spectrum. On the other hand, there exist well-established atomistic and molecular techniques (MD) that provide an accurate description at the nano-scales, but are hardly scalable at the macro-dimensions [24,25]. A good compromise is represented by the Boundary Element Methods (BEM): these approaches, developed in the last decade by different research groups in the world [26][27][28][29][30][31], allow to account for a large number of rough scales, considering at the same time what happens at the mesoscale. They have proved to be fully effective when dealing with materials exhibiting a time-dependent rheology: indeed, several investigations have been conducted to understand how roughness is related to hysteretic dissipation in sliding contacts [7]. Sliding motion has been widely studied due to its applicative importance and, in particular, to the so-called viscoelastic friction: the latter is related to the periodic deformation occuring during the sliding mechanics and is shown to be dependent on the entire roughness spectrum. To this extent, numerical calculations have qualitatively corroborated what was predicted on analytical basis by Persson's theory [1] and have generalized the analysis to different geometries, including thin layers [31], to several kinetic conditions, including the reciprocation of rough solids [31], and also to biphasic interfaces, with a number of studies enlightening the so-called visco-elasto-hydrodynamic regime [32][33][34][35].
Surprisingly, no such attention has been paid to evaluating the indentation, in the normal direction, of viscoelastic solids. In fact, following the original study by Lee and Radok on the indentation of Hertzian profiles, few contributions have been brought to this field. Normal contacts have, however, a specific applicative prominence: this is, for example, the case of creep tests [36][37][38]. Furthermore, let us cite, as sketched in Fig. 1, also the case of soft grasping, where human fingers -or, equivalently, artificial prostheses -are in normal contact with surfaces, which are nominally flats, but actually are rough on several orders of magnitude. Similar issues are found also in different areas, like industrial engineering, where, for example, vibrations in viscoelastic media are related, inter alia, to normal contacts [39][40][41][42]. In spite of the theoretical and practical interest in the theme, how rough scales influence normal contact is, however, something that has only recently started to be systematically studied [43], and multiple points are still obscure and out of our comprehension. This study aims at filling this gap in the understanding of the mechanics of rough interfaces: specifically, we will investigate how fractal parameters, including the fractal dimension, the root-meansquare and the mean slope of the heights, influence the normal contact features and, ultimately, the hysteretic losses. The paper is structured as follows. Section 2 includes the numerical methodology employed to carry out our analysis, while Sect. 3 shows the main numerical results and discuss how fractality is related to viscoelasticity in normal contacts. Final remarks close the paper.

Methodologies
In this section, we describe the numerical methodology implemented to analyze the viscoelastic indentation. Let us start by recalling that a linear viscoelastic material presents timedependent strain and stress distributions, being respectively equal to (t) and (t) and related throughout the following integral relation [44,45]: where the symbol ' ⋅ ' refers to the time derivative and J(t) is the so-called creep function. It is well know that the latter may be related to the material viscoelastic properties and specifically to the rubbery and glassy elastic moduli, being respectively defined as E 0 and E ∞ , to the creep spectrum C( ) and to the relaxation spectrum distribution . In fact, we can write J(t) as: where the Heaviside step function H(t) is introduced to satisfy the causality principle as J(t < 0) = 0 . Interestingly, Eq.
(2) can be easily discretized by redefining the creep where C k and k are the creep coefficients and the relaxation times. Now, we can deal with the contact problem and, specifically, on a 1D geometrical domain, where a rigid punch indents a viscoelastic half-plane. Let us notice that focusing on one-dimensional problem will alleviate computational efforts without penalizing the generality of the investigation; furthermore, 1D profiles are a qualitatively valid approach to deal with anisotropic surfaces, which are a common result of various manufacturing processes [48]. Now, given such a domain, we can rely on the translational invariance for the geometrical domain. Furthermore, the system is linear [46]. All this allows us to write the following integral equation between the normal surface displacement u(x, t) and the normal interfacial stress derivative ̇(x, ): where x is the position coordinate, t is the time, G tot is the global Green's function depending on the space and the time domains. Under the assumption of a perfectly homogenous solid, the integral equation kernel in Eq. (3) can be factorized in two terms, thus writing: where J(t) is the aforementioned creep function, whereas G(x) is a spatial Green's function, related to the distance between the points involved in the contact area. Different approaches can be followed to solve Eq. (4) depending on the kinetic conditions set in the contact problem. In general, as detailed in Ref. [46], in order to be solved, Eq. (4) requires to discretize both the time and the space domains, but, when analyzing particular sliding conditions, as for example in the case of steady-state [26] or reciprocating contacts [31], Eq. (4) can be rewritten as a convolution integral where there exists just a Green's function spatial term parametrically dependent on the steady-state speed [26] or on the time step during the reciprocating cycle [31]. Computational complexity is dramatically reduced in this way. In this paper, however, since we are interested in the normal indentation, we are not following a similar strategy and we will solve directly Eq. (4). Interestingly, as in an indentation process, the contact region is well defined,the spatial domain remains contained: the computational load is, consequently, still low enough to allow the study of rough interfaces. As suggested in Ref. [49], to approach Eq. (4), it is necessary to further develop the Green's term G(x) . In particular, as the aim of our investigation is to study rough profiles, it is convenient to impose periodic conditions: the most effective way is, to this aim, to define a periodic G (x) . This can be done by starting to Fourier-transform Eq. (4) as: and, further, again in the case of a linear viscoelastic homogeneous solid, as [50]: where and E( ) are respectively the Poisson ratio and the viscoelastic modulus in the Fourier domain, being easily related to the Fourier transformed creep function J( ) as We note that such a formulation is fully general and can be applied to any linear viscoelastic material, with any viscoelastic modulus E( ) and Poisson ratio . Furthermore, let us notice that S(q) is a correcting factor which takes into account the boundary conditions imposed on the viscoelastic substrate [50]: clearly, S(q) is unitary in the case of a half-plane. Now, as we need to impose periodic boundary conditions to a system that will have a spatial periodicity , as anticipated, it is necessary to define and determine a spatial periodic function G (x) . The latter is the displacement that will be found by imposing to the system a distribution of uniform pressures (x) : each term of the distribution will act over a strip with length a and will be periodically applied at a wavelength . On this basis, by adopting the numerical procedure detailed in Refs. [47][48][49], it is possible to determine G (x) as the following Fourier series: This may be easily computed by means of an IFFT (Inverse Fast Fourier Transform) algorithm. Now, Eq. (4), where the spatial Green's function can be set equal to G (x) , may be solved numerically. More details can be found in Appendix 1.
Now, once the mathematical formulation has been developed, it is crucial to define the features of the rough profiles under investigation. The 1D rough profile is described by the following Fourier series: where q is the long-distance cut-off q = 2 ∕ with being the fundamental periodicity length. The profile is generated numerically in a fractal self-affine form. To this end, let us assume for the power spectral density C(q), that is the Fourier-transfrom of the autocorrelation function of the profile, the following form: where q r and q c are the roll-off and the cut-off wavevector, being equal respectively to q r = N 0 q and to q c = Nq r : N 0 and N are, thus, the roll-off number and the number of scales. The presence of a roll-off vector is taken to improve the ergodicity of the system. Finally, C 0 is a constant to be chosen accordingly with the desired properties of the profile and, specifically, as shown later, with its statistical moments. For time being, let us focus on the numerical generation of the profile: we have to determine the amplitudes h k and the phases k . For the phases k , it is possible to choose random values uniformly distributed in the interval [− , [ : as shown in detail in Ref. [50], this ensures the invariance of the profile statistical properties. To obtain the amplitude, let us write the power spectral density (PSD) of the profile in Eq. (8): . Incidentally, we observe that the profile generated accordingly to this profile is periodic, thus allowing us to avoid any finite-size effect in our analyses.
Now, given the specific purpose of this paper, that is, to understand, on statistical basis, the influence of the profile features on the indentation solution, it is crucial to observe that, given the self-affinity of the process, just three independent parameters are necessary to determine its statistics. Our choice is to employ the following ones: (a) the mean square roughness ⟨ h 2 ⟩ of the profile, which coincides with the zeroth order moment of the PSD, that is, , which is the second order moment of the profile PSD, i.e., ⟨ h �2 ⟩ = m 2 = ∫ dqq 2 C(q) ; (c) the Hurst exponent H, that is related to the profile fractal dimension D f as D f = 2 − H . Incidentally, the choice of these parameters is, to some extent, arbitrary, but it has been driven by their straightforward relation to the PSD of the profile. Ultimately, now, our aim is to understand how the contact properties change with these parameters.

Results and Discussion
We study the indentation process of a 1D profile into a linear viscoelastic solid by controlling the mean separation s = s∕h rms throughout the following time law: where s 0 and Δs are respectively the dimensionless mean value s 0 = s 0 ∕h rms and amplitude Δs = Δs∕h rms of the oscillating indentation, and Ω and t are the dimensionless oscillation frequency Ω = and time t = t∕ with being the material relaxation time. In fact, in order to illustrate the main features of the system, but without any loss of generality, we employ a one-relaxation-time material with the glassy and the rubbery moduli respectively equal to E ∞ = 10 7 Pa and E 0 = 10 6 Pa, the Poisson's ratio = 0.5 and the relaxation time = 0.1 s. Regarding the rough profiles, in order to reach a statistical relevance, results are obtained by averaging the numerical outcomes of 150 different realizations and setting the roll-off number N 0 equal to N 0 = 4 ; the number of scales, included in the analysis, is N = 128 and, to obtain numerical convergence, the smallest scale is sampled with 32 points. Finally, we observe that the results are reported when all the transient effects have expired.
We start our study with a profile with h rms = 4.5 m, m 2 = 4.5 × 10 −4 and H = 0.6 : the analysis is carried out for s 0 = 0.35 , Δs = 0.1 ⋅s 0 and Ω = 6.8 . Figure 2 reports the deformed profile in terms of dimensionless displacement ũ = u∕h rms as a function of the dimensionless abscissa x = x∕ , for different time steps t after enough cycles have been completed to consider the transient concluded. Let us focus, now, on the main features marking the viscoelastic indentation compared to what would happen for a non-timedependent rheology: indeed, both in the full-size graph and in the zoom frames, we observe that, when a region is locally indented, once the contact has been released, the viscoelastic substrate keeps on showing there the shape assumed during the indentation. This memory effect is typical for a viscoelastic contact [49] and, crucially, is preserved also at the rough scales during the indentation cycle: as we will later see, this is intrinsically related to the work hysteretically dissipated in the indentation cycle. To this extent, in Fig. 3, for a given time-step and, precisely, for t = 3.25 , we can observe the pressure distribution: interestingly, the methodology is fully accurate and allow to obtain, in each contact cluster, a fully convergent distribution.
In order to assess the evolution of the entire indentation cycle, we can look at the global quantities and, specifically, at how the dimensionless separation s , the dimensionless real contact area ã = a∕ and the dimensionless global load per unit of lengthP = P∕ E 0 h rms evolve over the time. In Fig. 4a, we observe s as a function of P : as the load increases, the mean separation decreases, but correspondingly, with a diminishing load, separation soars. In this loading/unloading cycle, it is well clear that some energy is hysteretically dissipated. This is evident also in Fig. 4b and c, where we plot the contact area ã versus the load P and the separation s : the paths followed during the loading and the unloading process are clearly different due to the energy dissipation. Interestingly, at the two extreme deadpoints of the cycle, the trends are non-symmetrical own to the well-known non-linear relation between area and mean separation [48].
From this results, it is clear that the key quantity to consider during viscoelastic indentation is the dissipated work per unit of length W . The latter can be defined as a dimensionless quantity and, specifically, as the integral W = ∮Pds . In Fig. 5, where we plot W as a function of the frequency Ω , we observe the bell-shaped curve that is expected in processes involving viscoelastic dissipation: at the very low and very high frequencies, where the material retrieves an elastic behavior, W vanishes whereas, in the middle of the frequency range, we find a maximum for the dissipation.
At this stage, we can raise a question that is fundamental for the purposes of this paper and, specifically, we may wonder how the fractal features of the rough profiles influence the dissipation. As discussed in the previous section, a rough profile can be univocally identified on statistically basis by the mean square roughness h rms , the mean slope m 2 and the Hurst exponent H. In Fig. 6, we plot the maximum dissipated power W max during the cycle as a function of the aforementioned parameters. The values are chosen over a wide range and, specifically, are equal to h rms = 1 , 4.5, 20, 100 m, m 2 = 10 −4 , 4.5 × 10 −4 , 20 × 10 −4 , 10 −2 and H = 0.6 , 0.7, 0.8, 0.9. Incidentally, let us observe that the frequency corresponding to W max may vary slightly, but, as we will see, this will be due to the statistical scattering that affects the different realizations. In fact, the whisker boxes in Fig. 6 reveal that, despite a relatively large dispersion, mainly due to the different maximum height values in each numerical realization of the profile, there is no statistically significant difference for W max . This is a crucial result as it implies that, for rough profiles, the curve in Fig. 5 does not depend on the statistical properties of the profile. To better understand this, let us carry out a simple dimensional analysis. The total load per unit of length P can be estimated as P = m = |E( )| m with m being the mean stress and m the corresponding deformation. It is well known that the stress distribution mainly depends on the large scales, while the small ones produce just a perturbation of such a distribution. This implies that the total load P is mainly related to the large scales contribution: consequently, we can estimate m as m ≈ h rms ∕ and, on this basis, the load can be quantified as P ≈ |E( )|h rms . Now, as the separation s is on its turn proportional to h rms , i.e., s ≈ h rms , the dissipated work per unit of length W can be written as: and, since the dimensionless work W is in practise equal to W = W∕h 2 rms , we understand why the dimensionless dissipated power W does not appear to be influenced by the statistical parameters of the profile. In fact, for 1D profiles, the root-mean-square of the heights seems to be the only parameter that counts for the dissipated work. This is certainly important as h rms is mainly determined by the large scales, thus meaning that the small scales, in the case of alternate indentation, do not play a significant role in the viscoelastic dissipation. Such a behavior marks a significant difference with sliding contacts, where the leading parameter for viscoelastic friction is m 2 and, therefore, all the scales contribute to the viscoelastic friction [7]. This has interesting practical implications. In fact, an open problem in Tribology is to determine the cut-off vector when carrying out calculations on viscoelastic friction in sliding conditions: possible approaches, which relate the cut off to the material degradation or to the presence of interfacial micrometric particles, have been proposed in literature [1,51], but a definitive theory is still missing. On the contrary, this paper shows that, when dealing with normal indentation, the process is governed by the large scales and, thus, by the roughness h rms : the latter can be easily determined by an experimental profilometer measure.

Conclusions
In this paper, we investigate the normal indentation of a rough profile on a viscoelastic half-plane with the specific aim to understand how the fractal parameters influence the process and, precisely, the viscoelastic dissipation. To this end, we employ a boundary element methodology being capable of discretizing, with high accuracy, both the time and the space domains: this allows us to obtain the contact solution, at different time-steps, for profiles being rough on a large multi-scale range.
In our analysis, we focus on an indentation occuring with a sinusoidal time law. Interestingly, unlike what happens in the linear elastic case, we observe that, after being indented, the material needs time to relax and recover the past deformation: such a time-dependent behavior is strictly related to the energy hysteresis occuring during the periodic indentation cycle. This is well clear when studying the separation as a function of the load (or, equivalently, of the area): the loading path is different from the unloading one, as expected in a hysteric process. We have, therefore, computed the dissipated energy as a function of the indentation frequency and, for a simple one-relaxation-time materials, we have found the classic viscoelastic bell-shaped trend: at very low and very high frequencies, the material behaves elastically, while for intermediate values, we have the proper viscoelastic trend. At this stage, we can focus on the maximum dissipation value and on the crucial aim of this paper, that is, the influence of fractality on dissipation. Numerical results show that the dissipated work W, defined per unit of length in the 1D contact, depends exclusively on h 2 rms and, thus, mainly, on the large scales, which determine the root mean square of the height distribution. To this extent, the phenomenon is clearly different from what happens in sliding contacts, where the leading role is taken by m 2 and, so, by the entire rough spectrum. Interestingly, this results implies that, once made non-dimensional with respect to h 2 rms , the bell-shaped dissipation curve in the indentation of a 1D profile can be considered universal as it is independent on the mean square slope and the Hurst exponent. Clearly, this has experimental consequences as it means that dissipation will depend just on this simple roughness parameter.

Appendix 1: Numerical Implementation of the Boundary Element Approach
In order to implement a numerical approach to solve Eq. (4) and, specifically, a Boundary Element approach, let us mesh the one-dimensional spatial domain with N equallyspaced elements, having a length a, whilst the time domain is discretized with a step Δt . Thus, we can discretize Eq. (4) and, at a certain time step Δt , write the displacement u of the centroid x i as: To efficiently solve Eq. (14), let us isolate two terms that are already known at the time step Δt at which the equation is formulated. Specifically, these are b(x i , Δt) and c(x i , Δt): x j , kΔt − x j , (k − 1)Δt   included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.