Recent Advancements in the Tribological Modelling of Rough Interfaces

Abstract

This paper analyses some effective strategies proposed in the last few years to tackle contact mechanics problems involving rough interfaces. In particular, we present Boundary Element Methods capable of solving the contact with great accuracy and, at the same time, with a marked computational efficiency. Particular attention is paid to non-linearly elastic constitutive relations and, specifically, to a linearly viscoelastic rheology. Possible implications deal with all the tribological mechanical systems, where contact interactions are present, including, e.g., seals, bearings and dampers.

1. Introduction

In the last few decades, a variety of analytical, numerical and experimental approaches have been implemented to understand what happens when two solids, with real rough interfaces, come into contact. The reasons for such a marked attention in the rough contact mechanics are different and both theoretical and applicative. In fact, there is certainly a genuine academic interest in the problem due to its specific intricacy: the presence of the roughness, whose spectrum covers several orders of magnitude and goes down to the atomistic scales, introduces a huge number of space and time scales. On the other hand, rough contact mechanics has a practical importance in the optimized design of engineering systems and components: classical industrial applications include seals and dampers ([1,2]), but a constantly increasing interest is rising in the frontiers’ fields, such as bio-adhesive ([3,4,5,6]), cellular scaffolds ([7,8]) and even touch-screen devices ([9]).

The first theoretical answer to the rough contact problem was given by Greenwood and Williamson in 1966 in a pioneering paper where the first of the so-called multiasperity theories has beeen introduced. Following such a contribution, a variety of models have been proposed (see, e.g., [10,11,12,13]): basically, these models reduce the surface roughness to a discrete distribution of asperities behaving as independent Hertzian punches, thus neglecting the reciprocal interaction between the contact clusters. Due to this assumption, multiasperity theories cannot match experimental results in terms of applied load and contact area ([14]). In the last twenty years, Persson has developed a totally different approach, where the contact pressure probability distribution is demonstrated to be determined by a diffusive process, being dependent on the magnification at which the contact interface is observed. This theory can be considered exact in full contact conditions, but provides still qualitatively accurate information for partial contacts. On the other hand, to provide quantitativelly reliable predictions as needed in applications, a number of numerical methods and, in particular, several Boundary Element approaches, implemented either in the real ([15,16,17,18,19,20]) and in the Fourier space ([21,22,23,24]), have been developed. Nowadays, these techniques are extremely accurate, but at the same time, in most cases were affected by a significant limitation: they were developed for linear elastic contact mechanics.

However, it should be pointed out that, in many cases, and, crucially, when soft materials are involved in the problem, a marked non-elastic time-dependent mechanical behaviour is evident. Specifically, in a number of systems of applicative interest, including civil engineering ([2,25]) and biomechanics ([7,8,26]), soft matter can be described, with good approximation, by a linearly viscoelastic rheology. Recently, due to the theoretical and practical prominence of viscoelastic contact mechanics, a large amount of research activities was dedicated to the topic ([27,28,29,30,31,32,33,34,35]): different cases, including constant sliding velocity, reciprocating motion and lubricated contacts, have been investigated by some of the authors of this paper and other research groups in the world ([27,28,29,30,31,32,33,34,35,36,37,38,39,40,41]). There are still multiple issues to point out, but these studies have already properly demonstrated how viscoelasticity alters the contact solution in terms of contact area, load, separation, stiffness and, ultimately, friction compared to the elastic case. Results not only differ from the orders of magnitude from purely elastic conditions, but also have a qualitative impact on the global solution. To this extent, in this paper, we will review, inter alia, a point that, given its theoretical and practical importance, has to be properly accounted for: this is the anisotropy induced by the material viscoelastic rheology on the contact solution. As shown in Refs. [33,42], when a rigid isotropic rough punch slides over a viscoelastic layer, the contact solution, in terms of both contact clusters and displacement, is strongly anisotropic. In fact, in each contact cluster, we have a different behaviour between the leading edge and the trailing edge, where the material is not yet relaxed. Clearly when the contacting surfaces are already anisotropic, as it often happens in real interfaces due to the manufacturing treatment, the problem further complicates. All this has significant practical consequences: for example, this is crucial for sealing systems as it intervenes on the percolation phenomenon. In fact, most of the theories in the field ([43,44,45,46]) assume that the contact patches distribution is perfectly isotropic, but clearly this assumption fails when dealing with viscoelastic interfaces, thus leading to a potential underestimation of the percolating flow in rubber-based viscoelastic seals.

Another issue to properly account for when dealing with rough contacts is the case of coated bodies. Soft coatings offer the chance to tailor the resulting interface behavior in terms of adhesive toughness, local contact stiffness, frictional behavior, etc. Similarly, biological systems have also often evolved in order to exploit specific features of multi-layer tissues such as, for instance, human skin. For these reasons, besides the aforementioned classical investigations focusing on both adhesive [47,48,49,50,51,52,53] and adhesiveless [17,28,54,55,56,57,58,59,60,61] contacts involving half-spaces, detailed studies have been led specifically focusing on contacts of thin layers [35,62,63,64,65,66].

As a matter of fact, material dissimilarity between the contacting bodies (i.e., material coupling) is the only source of interaction between the displacement fields in the directions normal and tangential to the surface for half-space contacts [67,68]. This has been clearly pointed out in a series of studies dealing with both homogeneous [69,70,71] and graded [72,73] elastic materials. However, very little has been conducted with respect to the case of thin bodies. In Refs. [74,75,76], the aforementioned coupling description has been proven to hold true only for the case of semi-infinite contacts, whereas contacts involving deformable layers of finite thickness behave differently. Indeed, in this case, a thickness-related source of normal-tangential coupling exists, which has been defined as geometric coupling [77,78] (notably, the latter term is negligible for very thick bodies). In Refs [77,78], it has been demonstrated that the contact area is strongly affected by the geometric coupling in sliding frictional contacts, with predicted values up to 10% larger than the expected values for the uncoupled case. Electrical conductivity prediction [79], and wear estimation [80] are side problems which could result to be worsened by neglecting geometric coupling effects. Ref. [78] also shows that geometric coupling affects the overall frictional beahvior of the interface by inducing an asymmetric normal pressure distribution also in purely eleastic contacts. Similarly, the finite thickness of the deformable substrate is also expected to impact the adhesive performance of the contact interface. It is the case, for instance, of thin coatings in orthopedic implants [81,82], medical adhesive bands [83], and pressure sensitive adhesives [84,85,86], where the layer is comparable in thickness to the surface features size. It has been demonstrated that, in this case, the specific boundary conditions applied to the thin deformable layer may play a key role in determining the adhesive strength and toughness at the interface, both in the case of adhesive contact mechanics against wavy counterparts [64], and peeling detachment from flat substrates [87,88,89,90].

In this paper, we review the main approaches developed to tackle all the aforementioned issues and, in particular, to understand what happens in terms of pressure distribution, displacement and friction when rough interfaces are in contact. The paper is structured as follows. Section 2 include the methodology developed for elastic interfaces, while the following one is dedicated to the viscoelastic one. Results and final remarks complete the manuscript.

2. BEM Formulation for Elastic Sliding Contacts

The system under investigation is shown in Figure 1, where a randomly rough rigid surface is in steady-state sliding contact with a deformable solid backed onto a rigid substrate. In the same figure, $r\left(\mathbf{x}\right)$ represents the surface roughness (with $\mathbf{x}$ being the in-plane position vector) with periodicity $\lambda$, and h is, in general, the deformable thickness. Notably, for $h\to \infty$, the behavior of the contacting solids asymptotically approaches the half-space. Morevoer, as shown in Figure 1, we assume the rigid rough indenter to penetrate the solid surface by a quantity ${\delta }_z$, whereas with ${\overline{u}}_z$ and $\Delta$ we indicate the mean normal displacement of the solid surface and the mean penetration, respectively, so that:

$${\delta }_z=\Delta +{\overline{u}}_z.$$

(1)

In addition, $\mathrm{\Lambda }$ and $\lambda$ indicate the roughness peak and fundamental wave-length, respectively.

In order to generalize our contact model, we assume the presence of Amonton/Coulomb friction at the contact interface, with friction coefficient ${\mu }_c$, all over the contact domain $\mathrm{\Omega }$. Indeed, the stress distribution acting along the x direction is given by

$${\tau }_x\left(\mathbf{x}\right)={\mu }_{c}p\left(\mathbf{x}\right);\mathbf{x}\in \mathrm{\Omega },$$

(2)

where $p\left(\mathbf{x}\right)$ is the normal pressure distribution.

We assume that the relative sliding speed only occurs in the x direction. And that, ${\mu }_c$ does not depend on the relative sliding speed. Notably, the frictionless contact mechanics behavior (i.e., purely normal indentation) is easily reseambled for ${\mu }_c=0$.

Following the contact formulation given in Refs. [17,28,63,64], we formulate the contact problem in terms of interfacial surface displacement vector $\mathbf{v}=\left(v_x,v_y,v_z\right)$ and stress vector $\sigma =\left({\mu }_{c}p,0,-p\right)$. In the reference system joint to the moving indenter, we have

$$\mathbf{v}\left(\mathbf{x}\right)=\mathbf{u}\left(\mathbf{x}\right)-\overline{\mathbf{u}}={\int }_{\mathrm{\Omega }}{d}^{2}s\mathrm{\Theta }\left(\mathbf{x}-\mathbf{s},h\right)\sigma \left(s\right);\mathbf{x}\in \mathrm{\Omega },$$

(3)

where we applied the the coordinate transformation $x-Vt\to x$. In Equation (3), $\mathbf{u}$ and $\overline{\mathbf{u}}$ represent the total and mean surface displacement vectors, respectively, which can be calculated following the procedure indicated in Refs. [17,28,63,64]. Notably,

$$\overline{\mathbf{u}}=\frac{1}{{\lambda }^2}{\int }_{\mathrm{\Omega }}\mathbf{u}\left(\mathbf{x}\right)d^{2}x$$

(4)

with $\lambda$ being the size of the periodic calculation domain. Moreover, the term $\mathrm{\Theta }\left(x,h\right)$ in Equation (3) represents the Green’s tensor, which depends on the specific thickness h of the elastic body. Moreover, we define the mean contact pressure as

$$p_m=\frac{1}{{\lambda }^2}{\int }_{\mathrm{\Omega }}p\left(x\right)d^{2}x$$

(5)

Notably, by means of Equation (3), the contact problem can be reduced to a Fredholm equation of the first kind.Presently, let us focus our attention on the specific case of an half-space and a frictionless contact. To numerically solve the contact problem, the penetration depth $\Delta$ is controlled, while the computational domain D is discretized with small squares of non-uniform size. The unknown stress in each single square is assumed to be uniformly distributed on it. Thus, we can discretize Equation (3) as the following linear system:

$$v_i=L_{ij}{\sigma }_j$$

(6)

where ${\sigma }_i$ is the normal stress uniformely acting on the square, $v_i$ is the normal displacement at the centre of each square, and $L_{ij}$ is the elastic response matrix, that is, the matrix provided by the discretized version of Equation (3). In the case of, $L_{ij}$ can be computed by employing the Love solution (see [91]), which furnishes the elastic displacement due to a uniform pressure on a rectangular area, and adding up the contribution of each elementary cell D to account for the periodicity of the problem. Thus $L_{ij}=l_{ij}-l_m$, where:

$$\begin{array}{cc}\hfill l_{ij}& =\frac{1-{\nu }^2}{\pi E}\sum _{k=-\infty }^{+\infty }\sum _{h=-\infty }^{+\infty }\left\{\left({\xi }_{ij}+d_j\right)ln\left(\frac{\left({\eta }_{ij}+d_j\right)+{\left[{\left({\xi }_{ij}+d_j\right)}^2+{\left({\eta }_{ij}+d_j\right)}^2\right]}^{1/2}}{\left({\eta }_{ij}-d_j\right)+{\left[{\left({\xi }_{ij}+d_j\right)}^2+{\left({\eta }_{ij}-d_j\right)}^2\right]}^{1/2}}\right)\right.\hfill \\ \hfill & +\left({\eta }_{ij}+d_j\right)ln\left(\frac{\left({\xi }_{ij}+d_j\right)+{\left[{\left({\eta }_{ij}+d_j\right)}^2+{\left({\xi }_{ij}+d_j\right)}^2\right]}^{1/2}}{\left({\xi }_{ij}-d_j\right)+{\left[{\left({\xi }_{ij}-d_j\right)}^2+{\left({\eta }_{ij}+d_j\right)}^2\right]}^{1/2}}\right)\hfill \\ \hfill & +\left({\xi }_{ij}-d_j\right)ln\left(\frac{\left({\eta }_{ij}-d_j\right)+{\left[{\left({\xi }_{ij}-d_j\right)}^2+{\left({\eta }_{ij}-d_j\right)}^2\right]}^{1/2}}{\left({\eta }_{ij}+d_j\right)+{\left[{\left({\xi }_{ij}-d_j\right)}^2+{\left({\eta }_{ij}+d_j\right)}^2\right]}^{1/2}}\right)\hfill \\ \hfill & \left.+\left({\eta }_{ij}-d_j\right)ln\left(\frac{\left({\xi }_{ij}-d_j\right)+{\left[{\left({\xi }_{ij}-d_j\right)}^2+{\left({\eta }_{ij}-d_j\right)}^2\right]}^{1/2}}{\left({\xi }_{ij}+d_j\right)+{\left[{\left({\xi }_{ij}-d_j\right)}^2+{\left({\eta }_{ij}-d_j\right)}^2\right]}^{1/2}}\right)\right\}\hfill \end{array}$$

(7)

and

$$\begin{array}{cc}\hfill l_m& =\frac{1-{\nu }^2}{\pi E}{\left(\frac{d_j}{\lambda }\right)}^2\sum _{k=-\infty }^{+\infty }\sum _{h=-\infty }^{+\infty }\left\{\lambda (h+1)ln\left(\frac{k+1+{\left[{(h+1)}^2+{(k+1)}^2\right]}^{1/2}}{k-1+{\left[{(h+1)}^2+{(k-1)}^2\right]}^{1/2}}\right)\right.\hfill \\ \hfill & +\lambda (k+1)ln\left(\frac{h+1+{\left[{(k+1)}^2+{(h+1)}^2\right]}^{1/2}}{h-1+{\left[{(k+1)}^2+{(h-1)}^2\right]}^{1/2}}\right)\hfill \\ \hfill & +\lambda (h-1)ln\left(\frac{k-1+{\left[{(h-1)}^2+{(k-1)}^2\right]}^{1/2}}{k+1+{\left[{(h-1)}^2+{(k+1)}^2\right]}^{1/2}}\right)\hfill \\ \hfill & \left.+\lambda (k-1)ln\left(\frac{h-1+{\left[{(k-1)}^2+{(h-1)}^2\right]}^{1/2}}{h+1+{\left[{(k-1)}^2+{(h+1)}^2\right]}^{1/2}}\right)\right\}\hfill \end{array}$$

(8)

where $d_j$ is the size of the elementary boundary element and ${\xi }_{ij}=\left|x_j-x_i\right|+\lambda h$ and ${\eta }_{ij}=\left|y_j-y_i\right|+\lambda k$.

Clearly, Equation (6) will be exploited to compute the interfacial stresses once the displacements $v_i$ are known. As the problem under investigation belongs to the class of mixed boundary problems, we need to determine the real contact area. This can be performed iteratevely by implementing the following procedure: (1) fix the displacement ${\Delta }_i$, (2) evaluate the so-called bearing area as the intersection between the deformed elastic layer, calculated with respect to the elastic solution determined previously for the penetration ${\Delta }_{i-1}<{\Delta }_i$, and the rigid rough punch, (3) compute the displacements in the contact clusters as $v_i=h_i-h_{max}+\Delta$, where $h_i=h\left({\mathbf{x}}_i\right)$, $h_{max}$ the maximum height of the rough profile, (4) solve Equation (6) to assess the stress distribution ${\sigma }_j$ in the contact areas, (5) determine the displacements $v_i=L_{ij}{\sigma }_j$ out of the contact areas, (6) update the contact area at each iterative step by deleting the elements with negative pressure and adding those where there exists compenetration. It should be noted that we invert the matrix $L_{ij}$ only for those points belonging to the contact area: this leads to a strong reduction of the computational efforts. The numerical inversion of the Equation (6) is made by means of an iterative method based on a Gauss-Seidel scheme.

It should be noted that the scheme previously described can be employed only in the adhesiveless case as negative values for the pressure distribution are discarded during the iterative case. In the adhesive case, based on the energy balance defined in Ref. [63], under isothermal and frictionless (i.e., ${\mu }_c=0$) conditions, for any given value of the contact penetration $\Delta$, the contact domain $\mathrm{\Omega }$ can be calculated by requiring that, at equilibrium,

$${\left(\frac{\partial \mathcal{F}}{\partial \mathrm{\Omega }}\right)}_{\Delta }=0$$

(9)

where $\mathcal{F}=\mathcal{E}+\mathcal{A}$ is the total free energy, with

$$\mathcal{E}=\frac{1}{2}{\int }_{\mathrm{\Omega }}p\left(\mathbf{x}\right)v_z\left(\mathbf{x}\right)d^{2}x$$

(10)

being the interfacial elastic energy stored into the deformable body, and

$$\mathcal{A}=-\Delta \gamma A$$

(11)

being the adhesion energy with $\Delta \gamma$ being the work of adhesion, also referred to as the Duprè energy of adhesion.

Notably, Equation (9) can be numerically calculated across infinitesimal variations of the discretized contact domain $\mathrm{\Omega }$ in the direction normal to the local boundary $\partial \mathrm{\Omega }$.

We also observe that the present formulation can be extended to the case of frictional contacts by following the procedure defined in Ref. [49], and the fundamental solution derived in Ref. [92]. To this extent, let us focus on the case of a rigid 1D rough profile in contact with a thin elastic layer of thickness h.

Morevoer, the term $\mathrm{\Theta }\left(x\right)$ in Equation (3) takes the form

$$\mathrm{\Theta }\left(x\right)=\frac{1}{E}\left(\begin{array}{cc}{G}_{xx}& G_{xz}\\ G_{xz}& G_{zz}\end{array}\right),$$

(12)

where E is the Young’s modulus fo the elastic material, and according to Ref. [77] we have

$$\begin{array}{cc}\hfill G_{xx}\left(x\right)& =-\frac{2\left(1-{\nu }^2\right)}{\pi }\left[log\left|2sin\left(\frac{q_{0}x}{2}\right)\right|+\sum _{m=1}^{\infty }B\left(m{q}_{0}h\right)\frac{cos\left(m{q}_{0}x\right)}{m}\right],\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill G_{xz}\left(x\right)& =-G_{zx}\left(x\right)=\frac{1+\nu }{\pi }\left[\frac{1-2\nu }{2}\left[\mathrm{sgn}\left(x\right)\pi -q_{0}x\right]-\sum _{m=1}^{\infty }C\left(m{q}_{0}h\right)\frac{sin\left(m{q}_{0}x\right)}{m}\right],\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill G_{zz}\left(x\right)& =-\frac{2\left(1-{\nu }^2\right)}{\pi }\left[log\left|2sin\left(\frac{q_{0}x}{2}\right)\right|+\sum _{m=1}^{\infty }A\left(m{q}_{0}h\right)\frac{cos\left(m{q}_{0}x\right)}{m}\right],\hfill \end{array}$$

(15)

with $q_0=2\pi /\lambda$, and

$$\begin{array}{cc}\hfill A\left(m{q}_{0}h\right)& =1+\frac{2m{q}_{0}h-\left(3-4\nu \right)sinh\left(2m{q}_{0}h\right)}{5+2{\left(m{q}_{0}h\right)}^2-4\nu \left(3-2\nu \right)+\left(3-4\nu \right)cosh\left(2m{q}_{0}h\right)},\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill B\left(m{q}_{0}h\right)& =1-\frac{2m{q}_{0}h+\left(3-4\nu \right)sinh\left(2m{q}_{0}h\right)}{5+2{\left(m{q}_{0}h\right)}^2-4\nu \left(3-2\nu \right)+\left(3-4\nu \right)cosh\left(2m{q}_{0}h\right)},\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill C\left(m{q}_{0}h\right)& =\frac{4\left(1-\nu \right)\left[2+{\left(m{q}_{0}h\right)}^2-6\nu +4{\nu }^2\right]}{5+2{\left(m{q}_{0}h\right)}^2-4\nu \left(3-2\nu \right)+\left(3-4\nu \right)cosh\left(2m{q}_{0}h\right)}.\hfill \end{array}$$

(18)

3. BEM Formulation for Viscoelastic Sliding Contacts

In multiple applications, it is necessary to account for soft contacts and, specifically, for a linear viscoelastic rheology for the solids into contact. To this end, let us briefly recall the main features of linear viscoelasticity [93,94] throughout the following integral equation, which correlates two time-dependent quantities, that is, the strain $\epsilon \left(t\right)$ and the stress $\sigma \left(t\right)$:

$$\epsilon \left(t\right)={\int }_{-\infty }^{t}d\tau \mathcal{J}\left(t-\tau \right)\stackrel{·}{\sigma }\left(\tau \right),$$

(19)

where $\mathcal{J}\left(t\right)$ is the creep function and the symbol ‘·’ refers to the time derivative. Now, if we define the real quantities $E_0$ and $E_{\infty }$ respectively as the rubbery and glassy elastic moduli of the viscoelastic material, $\mathcal{C}\left(\tau \right)$ as a the positive defined as creep spectrum [93], and $\tau$ is the relaxation time distribution, the creep function $\mathcal{J}\left(t\right)$ can be written as:

$$\mathcal{J}\left(t\right)=\mathcal{H}\left(t\right)\left[\frac{1}{E_0}-{\int }_0^{+\infty }d\tau \mathcal{C}\left(\tau \right)exp\left(-t/\tau \right)\right]=\mathcal{H}\left(t\right)\left[\frac{1}{E_{\infty }}+{\int }_0^{+\infty }d\tau \mathcal{C}\left(\tau \right)\left(1-exp\left(-t/\tau \right)\right)\right]$$

(20)

where $\mathcal{H}\left(t\right)$ is the Heaviside step function introduced so that $\mathcal{J}\left(t\right)$ can satisfy the principle of causality and, thus, $\mathcal{J}\left(t<0\right)=0$.

Now, Equation (20), let us understand what happens when viscoelastic bodies are into contact. Due to the linearity of the system, as the geometrical domain is non-finite and, thus, translational invariant [91], let us focus on the following integral equation formulated to correlate the normal surface displacement $u\left(x,t\right)$ and the normal interfacial stress derivative $\dot{\sigma }\left(x^{\prime },\tau \right)$:

$$u\left(x,t\right)={\int }_{-\infty }^{t}d\tau \int d^2{x}^{\prime }{G}_{tot}\left(\mathbf{x}-{\mathbf{x}}^{\prime },t-\tau \right)\dot{\sigma }\left({\mathbf{x}}^{\prime },\tau \right),$$

(21)

where $\mathbf{x}$ is the position vector, t is the time, $G_{tot}$ is a global Green’s function. If we assume that we are dealing with a perfectly homogenous solid, we can factorize the integral equation kernel in Equation (21) in two terms, thus writing:

$$u\left(x,t\right)={\int }_{-\infty }^{t}d\tau \int d^2{x}^{\prime }\mathcal{J}\left(t-\tau \right)G\left(\mathbf{x}-{\mathbf{x}}^{\prime }\right)\dot{\sigma }\left({\mathbf{x}}^{\prime },\tau \right),$$

(22)

where $G\left(\mathbf{x}\right)$ is a spatial Green’s function, which is defined in detail in Ref. [95]. It is crucial to observe that, in Ref. [28], under steady-state assumptions, that means the punch is sliding at constant v over the viscoelastic substrate, Equation (22) has been further developed. In fact, we can introduce a speed parametrically dependent Green’s function $G\left(x,v\right)$, thus strongly simplifying Equation (22) as:

$$u\left(x\right)=\int d^2{x}^{\prime }\mathcal{G}\left(x-x^{\prime },v\right)\sigma \left(x^{\prime }\right)$$

(23)

Similarly, in Ref. [36], reciprocating conditions are exploited: the rough rigid punch is assumed in sinusoidal motion over the viscoelastic substrate. Indeed, in this case, we can introduce a Green’s function, being parametrically dependent on the time t, and we can write the following Equation:

$$u\left(x,t\right)=\int d^2{x}^{\prime }\mathcal{G}\left(x-x^{\prime },t\right)\sigma \left(x^{\prime },t\right)$$

(24)

When we reduce Equation (22) to Equation (23) or to Equation (24), respectively for the steady-state or the reciprocating conditions, a dramatic reduction in the computational complexity is obtained. As a result, multi-scale problems, where the roughness spectrum covers several orders of magnitude, can be investigated. On the other side, as these simplifications affect the generality of the kinematic conditions which can be investigated, in Ref. [95] Equation (22) has been directly tackled. Although the computational complexity has allowed the solution just for the 1D case, it has been possible to explore different conditions, including normal indentation and transient contacts [95].

Another aspect to consider when approaching viscoelastic contact problems deals with the case of thin viscoelastic layers. Let us focus here on the case of a 1D rough profile in sliding contact with a thin viscoelastic layer under the assumptions of frictional interactions at the interface (i.e., ${\mu }_c>0$). In this case, again other steady conditions, we can define a speed parametrically dependent tensor ${\mathrm{\Theta }}^V$ as conducted previously in Equation (3) for the purely elastic case. In particular, ${\mathrm{\Theta }}^V\left(x\right)$, with x being the position coordinate, can be defined as follows:

$${\mathrm{\Theta }}^V\left(x\right)=J\left(0^{+}\right)\left[\begin{array}{cc}{G}_{xx}\left(x\right)& G_{xz}\left(x\right)\\ G_{xz}\left(x\right)& G_{zz}\left(x\right)\end{array}\right]+{\int }_{0^{+}}^{+\infty }\left[\begin{array}{cc}{G}_{xx}\left(x+Vt\right)& G_{xz}\left(x+Vt\right)\\ G_{xz}\left(x+Vt\right)& G_{zz}\left(x+Vt\right)\end{array}\right]\dot{J}\left(t\right)dt,$$

(25)

which, as expected, parametrically depends on the sliding velocity V. Specifically, in Equation (25), the viscoelastic creep function $J\left(t\right)$ is given by (20), and the G terms are given by Equations (13)–(15).

Since the interface is adhesiveless, the solution strategy adopted to calculate the unknow contact area domain is the same as previously introduced.

Notably, a very recent study has also demonstrated that adhesive viscoelastic contacts in sliding conditions can be addressed by relying on an energy balance approach.

4. Results and Discussion

4.1. Two-Dimentional Viscoelastic Sliding Rough Contacts: The Role of Anisotropy

An aspect with really important implications from both a theoretical and an applicative point of view is the anisotropy related to the contact area and on the deformation field when a rigid surface is in sliding contact with a viscoelastic half-space [42]. A field, where this can become crucial, is the fluid leakage: percolation is really influenced by the contact anisotropic solution. To quantify all these effects, we can employ self-affine fractal rigid surfaces generated with spectral components in the interval $q_r<q<q_c$, where $q_r=2\pi /\lambda$, with $\lambda$ being the side of the square punch equal to $\lambda =0.01$$\mathrm{m}$, $q_r=N{q}_0$ and N number of scales (or wavelengths) [17]. Computations are carried out with $N=64$. Moreover, regarding the material properties, in the following developments, we employ a linear viscoelastic material with a single relaxation time $\tau =0.1$ s, with a high frequency modulus $E_{\infty }$ equal to $E_{\infty }={10}^8\mathrm{Pa}$, the ratio $E_{\infty }/E_0=11$ and the Poisson ratio equal to $\nu =0.5$.

Let us start showing in Figure 2 this anisotropic effect underlined for the first time in Ref. [36]. Basically, when we focus our attention on the sliding contact mechanics of a perfectly isotropic rough surface over a linearly viscoelastic solid, the contact area solution results are anisotropic: in detail, each contact cluster, as zoomed in the inset, tends to shrink at the trailing edge, where the material is still relaxing.Consequently, each contact area is stretched perpendicularly to the speed.

We can quantify the anisotropy degree by introducing, in a certain range of wave vectors ${\zeta }_1{q}_0<\left|\mathbf{q}\right|<{\zeta }_2{q}_0$, the symmetric anisotropy tensor for the deformed surface $\mathcal{M}\left({\zeta }_1,{\zeta }_2\right)$:

$$\mathcal{M}\left({\zeta }_1,{\zeta }_2\right)={\int }_{{\zeta }_1{q}_0<\left|\mathbf{q}\right|<{\zeta }_2{q}_0}{d}^{2}q\mathbf{q}\otimes \mathbf{q}{C}_d\left(\mathbf{q},{\zeta }_1,{\zeta }_2\right)$$

(26)

where $C_d\left(\mathbf{q};{\zeta }_1,{\zeta }_2\right)={\left(2\pi \right)}^{-2}\int d^{2}x〈u\left(\mathbf{0};{\zeta }_1,{\zeta }_2\right)u\left(\mathbf{x};{\zeta }_1,{\zeta }_2\right)〉exp\left(-i\mathbf{q}·\mathbf{x}\right)$ is the power spectral density of the deformed surface $u\left(\mathbf{x};{\zeta }_1,{\zeta }_2\right)$ that is filtered. Crucially, the band-pass filter in the interval $\left[{\zeta }_1,{\zeta }_2\right]$ is defined to highlight the frequencies where the viscoelastic effects are higher.

Now, let us notice that the quantity ${\mathcal{M}}_{ij}={\int }_{{\zeta }_1{q}_0<\left|\mathbf{q}\right|<{\zeta }_2{q}_0}{d}^{2}q{q}_i{q}_j{C}_d\left(q\right)$, with i and $j=1,2$, is the second order moments of the filtered surface power spectral density of the, i.e., ${\mathcal{M}}_{11}={\mu }_{20}=〈u_x^2〉$, ${\mathcal{M}}_{22}={\mu }_{02}=〈u_y^2〉$, ${\mathcal{M}}_{12}={\mu }_{11}=〈u_x{u}_y〉$, where $u_x=\partial u/\partial x$, $u_y=\partial u/\partial y$ (see also Ref. [96]). Thus, once we have defined the symmetric tensor $\mathcal{M}$, the quadratic form $\mathcal{Q}\left(\mathbf{x}\right)={\mathcal{M}}_{ij}{x}_i{x}_j$ can be introduced. In a polar reference system with $x=rcos\theta$, and $y=rsin\theta$, one may obtain:

$$\mathcal{Q}\left(\mathbf{x}\right)=r^2{\left|\nabla u·\mathbf{e}\left(\theta \right)\right|}^2=r^2{\mu }_2\left(\theta \right)$$

where $\mathbf{e}\left(\theta \right)$ is the unit vector $\left(cos\theta ,sin\theta \right)$ and

$${\mu }_2\left(\theta \right)={\mu }_{20}{cos}^2\left(\theta \right)+2{\mu }_{11}sin\left(\theta \right)cos\left(\theta \right)+{\mu }_{02}{sin}^2\left(\theta \right)$$

(27)

is, thus, the average square slope of the profile defined by carrying out a cut of the deformed surface $u\left(\mathbf{x};{\zeta }_1,{\zeta }_2\right)$ along the direction $\theta$ [33]. Thus, when plotting the quantity ${\mu }_2\left(\theta \right)$ in a polar diagram, if the system is isotropic, ${\mu }_2\left(\theta \right)$ has to be circular; otherwise, there exists a different elliptical shape. Furthermore, we can quantify the degree of anisotropy by looking at the ratio ${\gamma }_d={\mu }_{2min}/{\mu }_{2max}$ between the minimum ${\mu }_{2min}$ and the maximum ${\mu }_{2max}$ eigenvalues of the tensor $\mathcal{M}$; furthermore, the principal direction of anisotropy can be introduced by detecting the value of the angle ${\theta }_d$ maximizing ${\mu }_2\left(\theta \right)$, i.e., ${\mu }_2\left({\theta }_d\right)=$${\mu }_{2max}$.

Similarly, we can introduce the roughness anisotropy tensor $\mathbf{M}$ for the rigid surface as $\mathbf{M}\left({\zeta }_1,{\zeta }_2\right)={\int }_{{\zeta }_1{q}_0<\left|\mathbf{q}\right|<{\zeta }_2{q}_0}{d}^{2}q\mathbf{q}\otimes \mathbf{q}C\left(\mathbf{q},{\zeta }_1,{\zeta }_2\right)$: its components are, then, $M_{11}=m_{20}=〈h_x^2〉$, $M_{22}=m_{20}=〈h_y^2〉$, $M_{12}=m_{11}=〈h_x{h}_y〉$ with $h_x=\partial h/\partial x$ and $h_y=\partial h/\partial y$. As conducted before, it is possible to associate to the tensor $\mathbf{M}$ the quadratic form Q and the parameters ${\gamma }_s$ and ${\theta }_s$ for the quantity $m_2\left(\theta \right)$. The latter is the average square slope for the profile cut on the rough surface $h\left(\mathbf{x};{\zeta }_1,{\zeta }_2\right)$ along the direction $\theta$ [33].

If we focus on the isotropic surface sliding over a viscoelastic half-space as pointed out in Figure 2, it is possible to plot $m_2\left(\theta \right)$ and ${\mu }_2\left(\theta \right)$. In Figure 3, we observe that $m_2\left(\theta \right)$ is perfectly circular (with ${\gamma }_s=1$), while crucially ${\mu }_2\left(\theta \right)$ is elliptical with ${\gamma }_d=0.37$ and ${\theta }_d$ is approximately $\pi /2$. We quantify in this way what was clear in Figure 2: as the contact clusters are perpendicular to the velocity direction, the maximum anisotropy angle must be close to $\pi /2$. In Ref. [28], it was demonstrated that the contact area shrinkage and the anisotropic shape for the spectral moment ${\mu }_2\left(\theta \right)$ are correlated: the shrinkage of each contact cluster at the trailing edge applies a high-pass filter to the frequencies which corresponds to the scales along the velocity direction: thus, as observed in Figure 3, the spectral moment of the deformed profile reduces in such a direction. More details on the generalization of the anisotropy induced by viscoelasticity can be found in Ref. [42].

4.2. Adhesion in Elastic Contacts of Thin Layers

Another paradigmatic condition to consider when dealing with rough contact is the presence of thin layers. In this case, in order to simplify the calculations, without any loss of generality of the method employed, adhesive elastic conditions are investigated in the case of a 1D rigid wavy profile (with single wavelength $\lambda$ and amplitude $\mathrm{\Lambda }$) in adhesive contact against a linear elastic layer of thickness h. Moreover, since the elastic layer presents a finite thickness, specific boundary conditions can be imposed to the layer surface opposed to the contact interface. In this regard, we consider two different boundary conditions as reported in the insets of Figure 4: the confined case, where a rigid constraint is applied to the layer boundary (see Figure 4a); and the remote pressure case, which consider a free layer boundary with uniform pressure applied (see Figure 4b).

Figure 4 shows the dimensionless mean pressure ${\tilde{p}}_m$ as a function of the $\tilde{\Delta }$.

Referring to the confined configuration, Figure 4a shows that the contact behaves stiffer as the dimensionless thickness $\tilde{h}$ is reduced because of the presence of the upper rigid constraint which hampers the elastic deformation of the layer. A different scenario holds in the case of remote pressure configuration, which, at a relatively small value of $\tilde{h}$, behaves like an Euler-Bernoulli beam. Indeed, as shown in Figure 4b, the contact stiffness significantly reduces with $\tilde{h}$ reducing. This leads to a peculiar behavior as, in presence of external loads, a threshold thickness $h_{\mathrm{th}}$ exists below which partial contact cannot occur (i.e., jump into full-contact occurs). Notably, in the confined configuration, complete contact can occur only for $h>$$\mathrm{\Lambda }$.

Figure 4 also allows to appreciate the different behavior of the two configurations at pull-off, under load controlled conditions. Indeed, in the case of remote pressure configuration, reducing the layer dimensionless thickness $\tilde{h}$ leads to lower pull-off pressure and smaller (more negative) pull-off penetration. As shown in Ref. [64], this entails larger adhesive toughness, thus suggesting safety applications, where large amount of energy needs to be absorbed. On the contrary, the confined case shows increasing pull-off pressure with $\tilde{h}$ decreasing, which is peculiarly suited for structural applications (e.g., adhesives).

4.3. Frictional Elastic/Viscoelastic Sliding Rough Contacts of Thin Layers

In the case of sliding contacts with frictional interactions at the sliding interface, the elastic field in the deformable body depends on the distribution of the normal and tangential tractions at the interface. In the usual assumption of half-space contacts (i.e., $h\gg \lambda$) with rigid against incompressible (i.e., $\nu =0.5$) materials, the presence of in-plane stress (e.g., frictional) does not play any role, and the frictionless normal contact solutions still holds true (uncoupled case). However, Equation (3) shows that in the case of non-vanishing out-of-diagonal terms in the Green’s tensor $\mathrm{\Theta }$, the elastic fields caused by normal and in-plane stress interact with each other, and the linear superposition of the contact solutions is no longer possible (coupled case). In this case, the contact behavior depends on the specific distribution of both normal and tangential stress. This is the case, for instance, of frictional sliding contacts involving thin layers and/or compressible (i.e., $\nu <0.5$) materials.

The frictional contact of thin elastic layers can be investigated by exploiting the formalism given by Equations (3) and (12)–(18). The same approach can be adopted for the viscoelastic case, providing that the Green’s tensor component are calculated by means of Equation (25). Specifically, in what follows we focus on the case of a 1D rigid rough profile in contact with a thin layer of thickness h backed onto a rigid substrate (i.e., confined configuration). For most of the calculation, we assume incompressible material for the layer, so that the only source of normal-tangential coupling arises from the finite thickness of the deformable layer. In agreement with Refs. [77,78], we qualitatively refer to this as to geometric coupling. The rigid profile presents a self-affine roughness spanning over 100 scales with a periodicity wavelength $\lambda$ and Hurst exponent $H=0.8$. We also assume $r_{rms}=10$ $\mathsf{\mu }$m, and $\overline{g}=0.13$, as the profile root mean square height and slope, respectively.

Due to normal-tangential coupling, different interfacial displacements are expected for frictional and frictionless contact conditions, given the same contact configuration. In turn, this may give rise to a different value of the contact area a in the presence of interfacial friction compared to the value $a_0$ of the frictionless case. Indeed, Figure 5a shows the contact area ratio $a/a_0$ as a function of the dimensionless mean contact pressure ${\tilde{p}}_m$, in the presence of geometric coupling. We observe that, at low values of ${\tilde{p}}_m$, the effect of coupling on the contact area ratio is poor; whereas, for ${\tilde{p}}_m>2$, the contact area ratio significantly increases by increasing the value of ${\tilde{p}}_m$. Finally, at very high contact pressure, a saturation of the value of a occurs, as the full contact conditions is approached both for frictional and frictionless contacts. As expected, increasing the interfacial friction coefficient, leads to enhanced coupling effects. Indeed, in highly frictional contacts (e.g., rubber contacts with ${\mu }_c\approx 1$) the predicted contact area increase may raise up to 10% compared to the expected value in frictionless case. Moreover, a closer look at Equations (14) and (17) shows that the geometric coupling term is fast decaying with $\tilde{h}=q_{0}h$ increasing. This is confirmed by the data shown in Figure 5b, where we observe that the actual contact area increase predicted for $\tilde{h}=1.5$ is of only 2 %, compared to the uncoupled conditions. A further increase in the layer thickness leads to vanising coupling effects, as the confinement offered by the rigid substrate is very poor.

The frictional behavior of the interface is also affected by coupling. Indeed, from Equation (14) we observe that, in the presence of geometric coupling (i.e., for thin layers), the normal displacements on the layer surface under a normal point force are asymmetric even in the case of a purely elastic material. Hence, in rough contacts’ conditions, an asymmetric pressure distribution is expected on each contacting asperity, which eventually entails an additional friction force opposing the relative sliding between the indenter and the deformable body. In the case of viscoelastic materials, as clearly shown in Refs. [62,66], due to the delayed material relaxation, an asymmetric contact pressure distribution is expected even in the case of vanishing coupling (i.e., for $h\gg \lambda$). Therefore, in contacts involving thin viscoelastic layers, the pressure asymmetry can be ascribed to the combined actions of both coupling and viscoelasticity [78]. The friction term resulting from the degree of contact pressure asymmetry is usually referred to as interlocking friction. The overall friction coefficient $\mu$ experienced by the contacting bodies can be calculated as

$$\mu ={\mu }_c+{\mu }_a$$

where ${\mu }_c$ is the Coulomb friction coefficient at the sliding interface, and

$${\mu }_a=\frac{1}{\lambda p_m}{\int }_{\mathrm{\Omega }}p\left(x\right)u_z^{\prime }\left(x\right)dx$$

is the interlocking friction coefficient due to either coupling and/or viscoelasticity, with $u_z^{\prime }$ being the spatial first derivative of $u_z$.

Figure 6a shows the normalized friction coefficient ${\mu }_a/{\mu }_c\overline{g}$ induced by the asymmetry of the contact pressure distribution in purely elastic contacts as a function of the dimensionless mean contact presssure ${\tilde{p}}_m$, for different values of $\nu$. For $\nu =0.5$, due to geometric coupling, the degree of asymmetry of the contact pressure distribution is the highest possbile; moreover, since the pressure eccentricity is shifted in the direction of sliding, the resulting normalized friction coefficient ${\mu }_a/{\mu }_c\overline{g}>0$. Consequently, under these conditions, regardless of the specific value of ${\tilde{p}}_m$, the overall contact friction is higher than for uncoupled contacts (i.e., for incompressible half-space). For $\nu <0.5$, also material coupling occurs between normal and tangential displacements fields. The contact pressure eccentricity depends on the value of the contact mean pressure, therefore the frictional behavior of the contact may result in being increased or decreased with respect to the uncoupled corresponding case.

Figure 6b shows the normalized friction coefficient ${\mu }_a/\overline{g}$ as a function of the dimensionless sliding velocity $\zeta =V{t}_c{q}_0$, for the frictional (${\mu }_c=0.8$) and frictionless case (${\mu }_c=0$). In this case, the layer is assumed lineraly viscoelastic (with single relaxation time $t_c$). Firstly, we observe that at very high and very low values of the dimensionless sliding velocity, the normalized friction coefficient ${\mu }_a/\overline{g}$ presents its minima. Indeed, under these conditions, the viscoelastic hysteresis vanishes and the material response is barely elastic. We also observe that the presence of geometric coupling leads to higher values of ${\mu }_a/\overline{g}$ compared to uncoupled conditions. Moreover, since the coupling terms in Equations (14) and (17) do not explicitly depend on $\zeta$, non-vanishing values of ${\mu }_a/\overline{g}$ are reported even for $\zeta \to 0$ and $\zeta \to \infty$ when coupling occurs.

5. Conclusions

In this paper, we review a variety of methodologies dealing with rough contact mechanics. In particular, we focus on Boundary Element methods: these have really impacted the tribology community in the last fifteen years as they have provided a solution for both purely elastic and viscoelastic normal and sliding contacts. In particular, for viscoelastic solids, we find important practical implications in all the systems, where there occurs relative motion between the contacting viscoelastic bodies: these include countless possible applications, such as, for example, vibration isolators, dynamic seals, pick and place devices. From a numerical point of view, given the translation invariance and the linearity of these systems, these problems have been tackled by means of a convolution integral, possibly accounting for the time and the space domains or relying on particular kinetic conditions, such as steady-state or reciprocating motions. This has required one to develop the ad hoc defined Green’s functions. In a wider sense, and specifically, when thin layers are considered, an entire Green’s tensor has to be considered to account for the coupling between normal and tangential actions: this can lead to an increase in the friction force.

This paper demonstrates the necessity of developing proper numerical strategies, such as the BEM introduced in this paper, to have accurate interfacial information, while preserving computational efficiency.