A Review on the Contact Mechanics Modeling of Rough Surfaces in the Elastic Regime: Fundamentals, Theories, and Numerical Implementations

Abstract

This chapter reviews advanced models for solving the normal contact problem of two elastic bodies with rough boundaries. Starting from the fundamental formulation of Greenwood and Williamson, an extension is proposed with details on the possible algorithmic implementation to consider the interactions between asperities. A second multi-scale-based approach, considering the self-affine nature of the rough surface, also known as Persson’s theory, is briefly discussed. As a third method, special attention is given to review the standard Boundary Element Method (BEM). Finally, all the mentioned methods are applied to a rough gold surface measured by Atomic Force Microscope (AFM) and the evolution of the real contact area with loading is analyzed. The aim of this contribution is to present the basic guidelines to tackle the problem of contacting rough surfaces, accounting for the real surface topography.

1. Introduction

Among the most complicated problems in mechanics, there is the contact modeling of rough surfaces where one seeks to predict the stresses, strains, and true contact area. These predictions play a very important role in studying many phenomena such as friction [1], wear [2], thermal [3] and electrical conductivity [4], sealing [5], squeal [6, 7], etc. Understanding the micromechanical characteristics of these phenomena leads to a robust design of mechanical systems by increasing their performances. For example, in micro-systems with brittle behavior such as (i) radio frequency micro-electro-mechanical-systems (RF MEMS) with silicon-to-Silicon (Si-to-Si) contacts [8], (ii) micro-turbines using Si or Polydimethylsiloxane (PDMS)-based micro-valve and [9] (iii) Si-to-Si wafer bonding [10, 11, 12], the performances and the reliability of the device rely on the quality of the contact [13, 14, 15, 16, 17]. An accurate prediction of the pressure and the true contact area leads to a good dimensioning of the adhesion forces and thus, to an increase in the reliability of components.

As part of the design process of these systems, the contact problem is addressed by assuming that the contact interfaces are perfectly flat. This assumption neglects the notion of roughness, which often leads to an overestimation of the contact surface and an underestimation of the contact pressure. In reality, every surface, even mirror polished, still exhibits roughness at the micro- and nano-scales. A rough surface can be seen as the superposition of several wavelengths, forming hills (asperities or local maxima) and valleys (local minima), and it is fully described by two statistical functions: (i) the height probability density function (HPDF) and (ii) the power spectral density function (PSDF). The former describes the randomness of the rough surface, while the latter quantifies the contribution of the different wavelengths to the surface topography. Both HPDF and PSDF are widely used to characterize engineering surfaces [18, 19, 20] as well as to generate numerically artificial surfaces [21, 22].

When two rough surfaces come in contact, the mechanical load is borne by the top of the highest asperities. Thus, the real contact area is only a small fraction of the nominal area. The first and most popular theory that elucidated this point was made in the 1960s by the pioneers’ Greenwood and Williamson (GW) [23]. In such theory, the rough surface is represented by a set of spherically shaped asperities, with the same radius of curvature and whose height varies randomly following a given probability density function. Then, the Hertz theory of one elastic sphere contacting a rigid flat plane is extended to the entire surface. In the same period, Bush et al. [24] proposed a refined model approximately equivalent to its counterpart (i.e. GW theory). The special point with respect to GW theory is that Bush et al. assumed that asperities have a paraboloidal shape. Then, they used the random process theory [25] to describe the statistics of the rough surface. Other models quite similar to those described above can be found in Refs. [26, 27, 28]. Using the above theories, it has been found that, at a low squeezing load, the predicted true contact area $A$ increases linearly with the applied load.

The aforementioned contact models assume that all the contacting asperities behave independently. Thus, the interaction between the different contact region are neglected (i.e. each asperity is locally deformed without inducing a displacement of the neighboring asperity). This last assumption is not met for medium to high loading cases (i.e. nearly full contact case) since the contacting asperities become tighter and hence, the contribution of the surface deflections produced by each of them, at the level of the neighboring micro-contacts, becomes significant [29, 30]. On the other hand, O’Callaghan et al. [31] and Hendriks et al. [32] demonstrated experimentally that asperities can interact strongly and even coalesce at low loads questioning the reliability of the classic asperity-based models. To consider the interaction between asperities, authors in the literature have proposed several models based on GW theory. For instance, Paggi and Ciavarella [33, 34] included interaction between asperities by means of semi-analytical modeling to simulate the effect of the bulk. In the same context, Afferrante et al. [35] have introduced a multi-asperity model that takes into account the coalescence between two contacting asperities. Despite all these advanced models, the semi-analytical approach remains very approximate compared to numerical methods such as the finite or the boundary element method [36].

The original GW theory has also been criticized since it deals with special cases where the height distribution follows a Gaussian or exponential distribution. Especially, it is well admitted by the scientific community that the roughness can appear at many scales exactly like the fractal patterns, but down to a wavelength that corresponds to, perhaps, some lattice constant. The use of GW theory with the “fractal” characteristic has been questioned by many authors, in particular Persson [37]. The latter comes out with an ingenious idea, completely different from the asperity-based model, to tackle the scaling effect on realistic rough surfaces [37, 38, 39]. In contrast to asperity-based models, Persson’s theory considers the stress probability distribution as a function of the surface resolution. The idea is to determine the evolution of the contact stress density function when the PSDF is extended to cover a wide frequency range. Initially, the problem is solved by assuming a full-contact assumption (i.e. smooth surfaces assumption) under initial uniform pressure $p_0$. Then, step by step, the roughness is added according to the self-affine character carried by the PSDF which leads to the variations in the contact pressure [40]. According to Persson, the contact stress distribution satisfies a diffusion-like equation. Persson’s approach was also criticized by the leaders of the asperity-based models. They argued that Persson’s proof is not rigorous while deriving the diffusion equation [40]. Hence, the predicted contact area, using Persson’s theory, is smaller compared to the available results in the literature.

In parallel with analytical models, several research groups have proposed numerical simulations reproducing the behavior of rough surfaces in contact. Among these methods, there is the finite element method (FEM) introduced for the first time by Hyun et al. [41]. This category of method is free of any assumption. The only disadvantage it presents is the convergence which is related to (i) the numerical treatement (the penalty method for instance) and (ii) the discretization error. Indeed, to have significant results it is necessary that the size of the discretization is much smaller than the shortest wavelength of the rough surface [5]. However, the finest mesh leads to a prohibitive computational cost. The latest point was clearly demonstrated by Yastrebov et al. [36]. It was shown that the numerical error of the discretization depends on two ratios; (i) $\mathrm{\Delta }x/{\lambda }_s$ and (ii) ${\lambda }_l/L$ where $\mathrm{\Delta }x$ is the distance between the mesh points, ${\lambda }_s$ and ${\lambda }_l$ refer to the shortest and the longest wavelength of the surface PSDF, respectively. $L$ is the length of the representative elementary volume (i.e. the simulation box). Authors in [36] argued that previous ratios should tend towards zero for an accurate prediction and a good representation of the mechanical behavior of each asperity. However, if these criteria are no longer met, in particular the second ratio criterion, it is necessary to conduct several stochastic simulations so that the average prediction is relevant. To overcome these drawbacks, authors developed more attractive numerical techniques allowing to discretize only the rough surface but not the bulk. These methods are based on the Green Fuctions (GF) (fundamental solution) and were able to refine the mesh in order to perform more accurate simulation and thus, predict a reliable contact area—pressure law [42, 43, 44]. The first development of efficient Boundary Element Method (BEM) solvers goes back to the 1990s when Brandt et al. [45] suggested the multi-level multi-summation (MLMS) scheme allowing to simulate realistic rough surface within a reasonable computational cost. In the same context, Venner et al. [46] combined the MLMS solver with the full multi-grid (FMG) iteration method to solve elasto-hydrodynamically lubricated contact problem. It should be noted that the coupled scheme MLMS-FMG solver failed to converge when solving the dry contact problem with roughness [47] since the corresponding convergence theorems do not handle problems with multiple inequality constraints. In order to enhance the efficiency of the numerical scheme and to optimize the CPU cost, Polonsky et al. [42] enhanced MLMS solver by using a iteration scheme based on the conjugate gradient (CG) method showing good convergence properties for a 3D rough surface with large numbers of nodes (approximately ${10}^6$ nodes). Nogi et al. [48] and Polonsky et al. [49] have considered the Fast Fourier Transform (FFT) for solving rough contact problems for both homogeneous and layered solids. They demonstrated that the FFT-based method can increase dramatically the computational speed. Recently, Bemporad et al. [44] introduced the Non-Negative Least Squares (NNLS) algorithm with a suitable warm start based on accelerated gradient projections (GPs). The results in [44] clearly prove that the NNLS-GPs solver is faster than the previously mentioned solvers. It should be noted that the standard BEM is based on the linear elasticity assumption for a homogeneous medium. Hence, it can not handle the contact problem of rough boundaries including geometrical, interfacial, and material nonlinearities such as large deformation, adhesion, or plasticity. Sometimes, BEM can be adjusted to tackle the above-mentioned nonlinearities [50, 51, 52, 53] but for convenience in this type of problem, authors preferred FEM instead of BEM [54, 55].

This review paper provides a comprehensive comparison of the above-mentioned models with particular attention to the standard Boundary Element Method (BEM), which is based on the Green functions. The concern of the present study is limited to solve the frictionless rough contact problem of micro-systems with purely brittle behavior (i.e. in the elastic regime). Thus, the chapter is organized as follows. First, we start to recall and describe 4 theories, namely (i) the original GW model in its continuous and discrete form, (ii) the modified GW to take into account interaction between the contacting asperities, (iii) Persson theory, and (iv) the BEM. Then, the methods will be applied to a 1 μm-thick gold rough surfaces. The first one is measured by the atomic force microscopy (AFM) technique while the others are numerically generated using a suitable algorithm taking as an input the statistical functions of the measured Au rough surface (i.e. HDF and PSDF). Finally, the chapter will be closed with an analysis and discussion of the predicted contact area and pressure.

2. Theoretical methods of contact mechanics

2.1 The original Greenwood and Williamson theory

2.1.1 Continuous form

Consider a rigid perfectly smooth plane moving towards a rough surface of an elastic body with a motion $\delta$ (see Figure 1). In the GW theory, the rough surface is idealized by assuming that all summits (i.e. local maxima) are spherical with the same radius of curvature $R$. These are refereed to as asperities. The height of each asperity is treated as a random variable which follows a given probability law denoted by $h\left(z\right)$. The idea of GW theory is to solve the frictionless contact between the rough surface and the rigid flat plan by using the Hertz theory. In fact, for a spherical elastic body in contact with a rigid plan, Hertz theory gives the exact solution under the framework of the small strains. Let us denote $E_1$ and $E_2$ the Young modulus of the rough surface and the rigid plane, respectively. ${\nu }_1$ and ${\nu }_2$ are the corresponding Poisson coefficient. The contact area $A$ and load $F$ can be expressed in terms of the rigid body motion $\delta$ as,

Figure 1.

The Greenwood and Williamson theory. The rough surface is described by several spherical asperities of radius $R$. The function $g_n$ is the normal gap between the rigid plane and the mean one for the rough surface.

$$A=\mathit{\pi R\delta },$$E1

$$F=\frac{4}{3}{E}^{\ast }{R}^{1/2}{\delta }^{3/2},$$E2

where $E^{\ast }$ refers to the composite Young modulus. It’s given by the following expression,

$$\frac{1}{E^{\ast }}=\frac{1-{\nu }_1^2}{E_1}+\frac{1-{\nu }_2^2}{E_2}.$$E3

For a rigid plan, $E^{\ast }$ denotes the plan stress modulus.

Following the configuration in Figure 1, there will be contact between the rigid plane and the rough surface if and only if, the height $z$ of any asperity is greater than the normal gap $g_n$. This event $\left\{z>g_n\right\}$ has it own probability that can be written as,

$$\mathcal{P}\left(\left\{z>g_n\right\}\right)={\int }_{g_n}^{+\infty }h\left(z\right)\mathit{dz}.$$E4

Now, if the rough surface has $N$ asperities, then the expected number of contacting asperities $n$ will be,

$$n=N{\int }_{g_n}^{+\infty }h\left(z\right)\mathit{dz}.$$E5

The total number of asperities $N$ can be expressed as: $A_{0}D$, with $A_0$ is the nominal area of contact and $D$ represents the density of summits.

If we define the displacement of the asperities as $\omega =z-g_n$, then the true contact area $A$ and the expected total load $F$ are given by,

$$A=\mathit{\pi NR}{\int }_{g_n}^{+\infty }\left(z-g_n\right)h\left(z\right)\mathit{dz},$$E6

and

$$F=\frac{4}{3}{\mathit{NE}}^{\ast }{R}^{1/2}{\int }_{g_n}^{+\infty }{\left(z-g_n\right)}^{3/2}h\left(z\right)\mathit{dz}.$$E7

Note that one needs to know $A_0$, $D$, $R$ and the probability density of asperities $h$ to accurately apply the original GW model. While $A_0$ is obtained from the sample size, the topography parameters (i.e. $D$, $R$ and $h$) can not be directly measured. In the same framework, Nayak used the random process theory to statistically evaluate the topography parameters by the means of three spectral moments $m_0$, $m_2$ and $m_4$, also known as the mean square height, slope and curvature of the surface [25]. In the case of homogeneous, Gaussian and isotropic rough surfaces, the expression of the spectral moments is given by the following equation:

$$m_i={\int }_{-\infty }^{+\infty }{q}^i\text{PSDF}\left(q\right)\mathit{dq}=\frac{1}{L}{\int }_0^L{\left(\frac{d^{i}z}{{\mathit{dx}}^i}\right)}^2\mathit{dx},$$E8

where $L$ is the surface length, $i$ is a positive integer referring to the order of the moment and $q$ is the wavenumber (i.e. spacial frequency). Since data from measurement tools are discrete, McCool [56] proposed to use the finite difference approximations to evaluate the spectral moments. While $m_0$ is the variance,

$$m_0=\frac{1}{N_p-1}\sum _{i=1}^{N_p}{\left(z_i-\overline{z}\right)}^2,$$E9

the second and fourth spectral moments can be calculated as,

$$m_2=\frac{1}{N_p-1}\sum _{i=1}^{N_p-1}{\left(\frac{z_{i+1}-z_i}{\mathrm{\Delta }}\right)}^2,$$E10

$$m_4=\frac{1}{N_p-2}\sum _{i=1}^{N_p-2}{\left(\frac{z_{i+2}-2z_{i+1}+z_i}{\mathrm{\Delta }}\right)}^2,$$E11

with $N_p$ the total sampling points and $\mathrm{\Delta }$ the sampling length (i.e. $N\mathrm{\Delta }=L$) and $\overline{z}$ the mean surface height. Finally, the parameters $D$, $R$ and $h$ can be obtained using the random process theory as described by Nayak in [25].

2.1.2 Discrete form

GW theory can be rewritten in a discrete form. To do so, let us assume now that each asperity is characterized by its own radius $R_i$ and height $z_i$. Then, the displacement of the $i$-th asperity can be expressed using Macaulay bracket as,

$${\omega }_i=<\delta -\left(z_M-z_i\right)>,$$E12

where $z_M$ is the maximum height of the rough surface. We have $z_M=max\left(z_i\right)$.

The true contact area $A$ as well as the expected total load $F$ are given by,

$$A=\pi \sum _{i=1}^N{R}_i{\omega }_i,$$E13

$$F=\frac{4}{3}{E}^{\ast }\sum _{i=1}^N{R}_i^{1/2}{\omega }_i^{3/2}.$$E14

Note that, in the discrete GW formulation, the gap function is given by $g_n=\left(z_M-\overline{z}\right)-\delta$ where $\overline{z}$ is the mean plane of the rough surface.

2.2 Greenwood and Williamson model with interaction

The original GW model has been criticized for various reasons. First, in terms of the shape of the asperities, several authors have argued that in reality the asperity has a shape that tends towards a sinusoidal curve. In fact, experimental investigations has clearly shown that a sinusoidal description is more realistic than a circular one [57]. Second, in terms of computational point of view, sinusoidal asperities have interesting mathematical properties such as continuity and differentiability, which is not the case for the circular ones where stress concentrations occur and therefore, become problematic to tackle the plasticity behavior [58]. Third, with respect to asperities interaction, the GW theory assumes that each asperity deforms in an independent manner leading to an underestimation of the predicted contact pressure.

The main idea of this contribution is take into account the interaction between the asperities. All these developments will be with an asperity of sinusoidal form (see Figure 2).

Figure 2.

The modified Greenwood and Williamson theory. The rough surface is described by several sinusoidal asperities of a height $h_i$ and and wavelength ${\lambda }_i$. The function $g_n$ is the normal gap between the rigid plane and the mean one for the rough surface.

2.2.1 Modified Hertz theory for a single sinusoidal asperity

First, let us recall some basic results within the Hertz theory framework of two contacting bodies: flat surface and a sinusoidal body. The Young and the Poisson ratio for the deformable sinusoidal body are denoted $E_1$ and ${\nu }_1$, respectively. Assuming that a far field displacement $\delta$ is applied to the rigid plan. When the contact is activated, the rigid body motion $\delta$ will be accommodated by the displacement of the asperity $\omega$ and the displacement of the substrate $u_z$. We write,

$$\delta =\omega +u_z.$$E15

In Eq. (15), the $u_z$ term can be evaluated analytically [59] or numerically by approaching the displacement of the elastic half-space using curve fitting approach [58]. In this chapter, Boussinesq-like form will be used to compute $u_z$. It is given by:

$$u_z=\frac{1-{\nu }_1^2}{\pi E_{1}a}F,$$E16

where $a$ is the radius of the contact spot and $F$ is the contact load. Note that Eq. (16) is quite similar to the integral formulation of the Boundary Element Method.

Following Johnson work [60], the contact load, $F$, as well as the compression of the asperity, $\omega$, can be given through,

$$F=\sum _{i=1}^n\frac{4A_i{E}^{\ast }{\mathit{ia}}^{2i+1}2\cdot 4\cdots 2i}{\left(2i+1\right)1\cdot 3\cdots \left(2i-1\right)},$$E17

$$\omega =\sum _{i=1}^n\frac{2\cdot 4\cdots 2i}{1\cdot 3\cdots \left(2i-1\right)}{A}_i{a}^{2i},$$E18

where $E^{\ast }$ refers to the plan stress Young modulus and $A_i$ represents the expansion coefficient of Taylor series at the vicinity of the contacting point. They depend closely on the height $h_i$ and the wavelength ${\lambda }_i$ of each asperity (see Figure 2).

If we confine ourselves to the second order, the contact load and the compression of the asperity will be written as follows,

$$F=\frac{8}{3}{A}_1{E}^{\ast }{a}^3+\frac{64}{15}{A}_2{E}^{\ast }{a}^5,$$E19

$$\omega =2A_1{a}^2+\frac{8}{3}{A}_2{a}^4.$$E20

Note that when only the first order of the Taylor expansion is considered, Eq. (19) is simplified to Hertz formula given by Eq. (2).

2.2.2 Asperity interaction model

To account for the interactions, the Boussinesq solution will be added to Eq. (15). We write,

$${\delta }_i={\omega }_i+u_{z,i}+u_{z,i}^j.$$E21

The additional term, $u_{z,i}^j$, reefers to the interaction between the $i$-th and $j$-th asperities. We have,

$$u_{z,i}^j=\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^n\frac{1-{\nu }_1^2}{\pi E_1{r}_{\mathit{ij}}}{F}_j,$$E22

where $r_{\mathit{ij}}$ refers to the euclidean norm between the $i$-th and $j$-th asperities. $u_{z,i}^j$ may read as the displacement of the $i$-th asperity caused by the $j$-th contacting asperity whose contact force is $F_j$. The latter can be computed using Eq. (19). In this case, the contact radius will stand for the asperity $j$.

Therefore, the resolution of the contact problem with interaction is to find a set of contact spots radii $a=\left(a_1,\cdots ,a_n\right)$, for each load increment ${\delta }_i$, that minimizes the objective function $\mathrm{\Theta }\left(a\right)$ such as:

$$\left\{\begin{array}{l}\mathrm{\Theta }\left(a\right)=\frac{1}{2}\sum _{i=1}^n{\left(r_i\right)}^2=\frac{1}{2}\parallel r\left(a\right){\parallel }_2^2\\ r_i={\delta }_i-2A_1{a}_i^2-\frac{8}{3}{A}_2{a}_i^4-\frac{1-{\nu }_1^2}{\pi E_1{a}_i}{F}_i-\sum _{\underset{j\ne i}{j=1}}^n\frac{1-{\nu }_1^2}{\pi E_1{r}_{\mathit{ij}}}{F}_j\\ F_i=\frac{8}{3}{A}_1{E}^{\ast }{a}_i^3+\frac{64}{15}{A}_2{E}^{\ast }{a}_i^5\\ F_j=\frac{8}{3}{A}_1{E}^{\ast }{a}_j^3+\frac{64}{15}{A}_2{E}^{\ast }{a}_j^5\end{array}\right.$$E23

Since the problem in Eq. (23) is non linear, Levenberg-Marquardt (LM) algorithm is used to find the contact radii. In fact, the resolution of the contact problem consists in minimizing the objective function $\mathrm{\Theta }\left(a\right)$ defined in Eq. (23) using an iterative scheme. The global procedure can be split in $7$ phases summarized in the following flowchart:

For each contact step do,

1. Contact detection: detect all the contacting asperities;

2. Initialization: initialize the contact radii vector $a_i$. If $i=0$, It can be initialized by 0. Otherwise, it can be warm-started by taking into account the loading history of the previous iteration $i-1$;

3. Jacobian and Hessian matrices construction: compute the Jacobian, $\mathbf{J}$, and Hessian, $\mathbf{H}$, matrices such as the residual function at iteration $i+1$ is defined as follows:

$$\mathrm{\Theta }\left(a_{i+1}\right)=\mathrm{\Theta }\left(a_i\right)+\nabla \mathrm{\Theta }{\left(a_i\right)}^T\cdot \mathrm{\Delta }a+\frac{1}{2}\mathrm{\Delta }{a}^T\cdot \mathbf{H}\cdot \mathrm{\Delta }a,$$E24

where $\mathrm{\Delta }a$ is the radii step.

Since the contact problem is formulated as a nonlinear least-squares (NLS) problem, the gradient of the objective function can be written as,

$$\nabla \mathrm{\Theta }\left(a\right)=\mathbf{J}{\left(a\right)}^{T}r\left(a\right),$$E25

where $J$ denotes the Jacobian matrix. It is given by,

$$\mathbf{J}\left(a\right)=\nabla r{\left(a\right)}^T=\left(\begin{array}{c}\nabla r_1^T\\ \nabla r_2^T\\ ⋮\\ \nabla r_n^T\end{array}\right).$$E26

The Hessian matrix is the square matrix of second-order partial derivatives of $r\left(a\right)$. We write,

$$\mathbf{H}=\left(\begin{array}{cccc}\frac{{\partial }^2\mathrm{\Theta }}{\partial a_1^2}& \frac{{\partial }^2\mathrm{\Theta }}{\partial a_1\partial a_2}& \cdots & \frac{{\partial }^2\mathrm{\Theta }}{\partial a_1\partial a_n}\\ \frac{{\partial }^2\mathrm{\Theta }}{\partial a_2\partial a_1}& \frac{{\partial }^2\mathrm{\Theta }}{\partial a_2^2}& \cdots & \frac{{\partial }^2\mathrm{\Theta }}{\partial a_2\partial a_n}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{{\partial }^2\mathrm{\Theta }}{\partial a_n\partial a_1}& \frac{{\partial }^2\mathrm{\Theta }}{\partial a_n\partial a_2}& \cdots & \frac{{\partial }^2\mathrm{\Theta }}{\partial a_n^2}\end{array}\right).$$E27

For this case, the Hessian matrix can be rewritten as follows,

$$\mathbf{H}\left(a\right)=\mathbf{J}{\left(a\right)}^{T}J\left(a\right)+\sum _{i=1}^n{r}_i{\nabla }^2{r}_i.$$E28

LM algorithm differs from Newton’s or Gauss-Newton’s (GN) method by using a damped approximation of the Hessian matrix. In fact, in LM method GN approximation is damped using a suitable controlling strategy to avoid the weakness of the GN method when the Jacobian matrix is rank-deficient. We write,

$$\mathbf{H}\approx {\mathbf{J}}^T\mathbf{J}+\lambda \mathbf{I},$$E29

where $\lambda$ is the damping factor and $I$ is the identity matrix.

1. 4. Compute the radii step: the goal of this step is to find the direction $\mathrm{\Delta }a$ such as

   $$\left({\mathbf{J}}^T\mathbf{J}+\lambda \mathbf{I}\right)\mathrm{\Delta }a=-2{\mathbf{J}}^T{r}_i.$$E30

   Note that Eq. (30) is the normal equations for the following linear least-squares problem,

   $$\underset{\mathrm{\Delta }a}{min}\frac{1}{2}{‖\left[\begin{array}{c}\mathbf{J}\\ \lambda \mathbf{I}\end{array}\right]\mathrm{\Delta }a+\left[\begin{array}{c}r\\ 0\end{array}\right]‖}_2^2.$$E31

   Note also that the damping factor $\lambda$ affect both the search direction and the step size $\mathrm{\Delta }a$. In this research work, the trust-region (TR) method is used, after each iteration, to control and update the damping factor $\lambda$.

2. 5. Update the solution: compute $a_{i+1}$ such as

   $$a_{i+1}=a_i+\mathrm{\Delta }a.$$E32

3. 6. Evaluate the objective function: compute the residuals $r$ at the new solution $a_{i+1}$. Compute the new objective function ${\mathrm{\Theta }}_{\mathit{new}}=\mathrm{\Theta }\left(a_{i+1}\right)$ at the new solution $a_{i+1}$

   If ${\mathrm{\Theta }}_{\mathit{new}}<\mathrm{\Theta }\left(a_i\right)$, accept the new solution $a_{i+1}$. Otherwise, reject the radii step $\mathrm{\Delta }a$, keep the old parameter (the guess, the residual and the objective function) and adjust the damping factor using TR algorithm.

4. 7. Convergence verification: Compute the error. If the procedure has converged, return the contact radii $a_{i+1}$ for the current contact step. If the procedure had not yet converged but the radii step $\mathrm{\Delta }a$ was accepted, compute the Jacobian matrix at the new solution $a_{i+1}$, then go to the step 4.

Overall, the resolution of the contact problem with interaction falls into two loops. The first one handles the contact steps defined for each far field displacement $\delta$. Thus, for each step, the NLS problem is defined. The second loop solves the NLS problem by solving a sequence of linear least-squares sub-problem.

2.3 Persson scaling theory

In 2001, Persson [39] came out with an ingenious idea, completely different compared the asperity-based models, to tackle the micro-mechanic contact problem of rough surfaces. The Persson theory starts by considering the contact of the two rough bodies to be initially perfect (i.e. smooth contacting surfaces see Figure 3). For the full contact set up, the contact pressure $p_0$ is assumed uniform. Following the self-affine behavior carried by the PSDF, the roughness is progressively added which lead to a variation of the contact pressure and the gap. It’s important to note that when the PSDF covers a wide wavenumber range from the roll-off distance to the cut-off limit, Persson argued that the evolution of the pressure is given by its density function $P\left(p,\zeta \right)$, where $\zeta$ is the magnification factor which is related to the roll-off wavelength (or wavenumber). We have $\zeta =k/k_0={\lambda }_0/\lambda$. When the roughness is added, the function $P\left(p,\zeta \right)$ verifies the following diffusion equation,

Figure 3.

Two elastic bodies in contact. Illustration of the scaling theory. Both bodies has roughness on many different length scales.

$$\left\{\begin{array}{l}\frac{\partial P\left(p,\zeta \right)}{\partial \zeta }=f\left(\zeta \right)\frac{{\partial }^{2}P\left(p,\zeta \right))}{\partial p^2}\\ P\left(p,1\right)=\delta \left(p-p_0\right)\\ P\left(-\infty ,\zeta \right)=0\\ P\left(+\infty ,\zeta \right)=0\end{array}\right.$$E33

where,

$$f\left(\zeta \right)=G^{\prime }\left(\zeta \right)p_0^2.$$E34

$G^{\prime }$ is a function of the magnification. We write,

$$G\left(\zeta \right)=\frac{\pi }{4}{\left(\frac{E^{\ast }}{p_0}\right)}^2{\int }_{k_0}^{\zeta k_0}{k}^3\text{PSDF}\left(k\right)\mathit{dk},$$E35

from which its follows,

$$f\left(\zeta \right)=\frac{\pi }{4}{E}^{{\ast }^2}{k}_0{k}^3\text{PSDF}\left(k\right).$$E36

$k_0$ is the wavenumber related to the roll-off distance and the quantity $\pi {\int }_{k_0}^{\zeta k_0}{k}^3\text{PSDF}\left(k\right)\mathit{dk}$ is the second profile PSDF moment, $m_2$, for an isotropic surface.

It is important to note that Eq. (33) is supposed to hold for the full contact condition (first configuration of the Figure 3). For the partial contact, the 3-$\mathit{rd}$ boundary condition $P\left(-\infty ,\zeta \right)=0$ should be replaced by,

$$P\left(0,\zeta \right)=0.$$E37

2.4 Boundary element method

The Boundary Element Method (BEM) belongs to the most efficient numerical methods for solving normal contact problems between two linear elastic bodies with rough boundaries. Unlike asperity-based models, the BEM is free of any kind of assumption. Its only restriction is in the discretization of the rough surface. Indeed, an adequate mesh is required to guarantee the convergence of the predicted results. However, with the development of computer tools, the restriction of discretization is easily overcome.

The relevance of BEM is based on the treatment of the problem of contact of rough surfaces without the need to discretize the bulk. The fundamental BEM formulation deals with the frictionless contact of homogeneous and linear elastic bodies with rough boundaries. Several developments have been proposed to extend the BEM to address the contact problem of heterogeneous bodies with rough boundaries [61] including severe nonlinearities such as adhesion [50, 51], friction [52] or plasticity [53]. The purpose of this section is to present the basic BEM formulation to solve frictionless contact problem of rough surfaces.

Let us assume a system defined in Figure 4 which is equivalent of two bodies with two rough boundaries. The Young modulus $E_i$ and the Poisson ratio ${\nu }_i$ refer to the i-th body with $i=1$ or 2. In the configuration presented in Figure 4, an elastic half plan comes in contact with a rigid rough surface when the far field displacement $\delta$ is imposed. The profile of the rough surface is characterized by the function $\zeta \left(x\right)$ measured with respect to the mean plane. The peak of the rough surface is the ${\zeta }_{\max }$ such as ${\zeta }_{\max }={max}_{x\in S}\left(\zeta \left(x\right)\right)$. When the asperity comes in contact with the elastic half-plane, a displacement $\overline{u}$ is defined. It characterizes the distance between the local peak of the contacting asperity and the position of the elastic half-plane of the undeformed configuration, which is also equal to ${\zeta }_{\max }-\zeta$ with $\zeta$, in this case, denotes the local peak of the contacting asperity. Otherwise, a generic displacement of the elastic half-plane is simply denoted by $u\left(x\right)$.

Figure 4.

Contact problem configuration: rigid rough surface and an elastic half-plane. The latter moves normally towards the rigid rough surface according to a far field displacement $\delta$.

In the fundamental BEM, the normal displacement field, $u\left(x\right)$, at the coordinate $x$ is related to the contact pressure $p\left(y\right)$ at the coordinate $y$ by the following integral form:

$$u\left(x\right)={\int }_S{G}_r\left(x,y\right)p\left(y\right)\mathit{dy},$$E38

where $G_r\left(x,y\right)$ is the Green-function which links the displacement $u\left(x\right)$ at $x$ to the contact pressure $p\left(y\right)$ acting at y. It is given by the Boussinesq solution as follows:

$$G_r\left(x,y\right)=\frac{1}{\pi }\left(\frac{1-{\nu }_1^2}{E_1}+\frac{1-{\nu }_2^2}{E_2}\right)\frac{1}{r},$$E39

where $r$ is the distance between the points $x$ and $y$. It refers to the standard Euclidean norm. We have, $r=\parallel x-y\parallel$.

Therefore, the frictionless contact problem of the rough surfaces is solved by the following optimization problem,

$$\begin{array}{ll}\text{For}& \mathrm{a}\text{given}\mathrm{far}\text{field\hspace{0.17em}displacement}\delta \\ \underset{x\in S}{\text{find}}& u\left(x\right)\text{and}p\left(x\right)\\ \text{such}\mathrm{as}& u\left(x\right)={\int }_S\frac{1}{\pi }\left(\frac{1-{\nu }_1^2}{E_1}+\frac{1-{\nu }_2^2}{E_2}\right)\frac{1}{\parallel x-y\parallel }p\left(y\right)\mathit{dy}\\ \text{with}& \left\{\begin{array}{l}u\left(x\right)-\overline{u}\left(x,\delta \right)\ge 0\\ p\left(x\right)\ge 0\\ \left(u\left(x\right)-\overline{u}\left(x,\delta \right)\right)p\left(x\right)=0\end{array}\right.\end{array}$$E40

The discretized form of the above formulation can be given if we consider the following set of the discretized parameters,

1. the barycentric coordinate: $x_{i,j}=\frac{1}{S_{\mathit{ij}}}{\int }_{x\in S_{\mathit{ij}}}x\mathit{dA}$;

2. the average height: ${\zeta }_{i,j}=\frac{1}{S_{\mathit{ij}}}{\int }_{x\in S_{\mathit{ij}}}\zeta \left(x\right)\mathit{dA}$;

3. the resultant contact load: $p_{i,j}={\int }_{x\in S_{\mathit{ij}}}p\left(x\right)\mathit{dA}$;

4. the displacement: $u_{i,j}=\frac{1}{S_{\mathit{ij}}}{\int }_{x\in S_{\mathit{ij}}}u\left(x\right)\mathit{dA}$.

Note that $S_{\mathit{ij}}$ refers to a cell of the surface $S$. The latter can be seen as a square grid spacing by $\mathrm{\Delta }$ and containing $N\times N$ points. Therefore, the cell $S_{\mathit{ij}}$ is a square surface of dimension ${\mathrm{\Delta }}^2$.

Hence, the discretized form of Eq. (38) can be written as follows,

$$u_{i,j}=\sum _{k=1}^N\sum _{l=1}^N{G}_{r_{i-k,j-l}}{p}_{k,l},\forall \left(i,j\right)\in I_N=\left\{1,\cdots ,N\right\}\times \left\{1,\cdots ,N\right\}$$E41

where $G_{r_{i-k,j-l}}$ represents the averaged Green function over the cell $S_{\mathit{ij}}$. We write,

$$G_{r_{i-k,j-l}}=\frac{1}{{\mathrm{\Delta }}^2}{\int }_{S_{\mathit{ij}}}{\int }_{S_{\mathit{kl}}}{G}_r\left(x,y\right)\mathit{\text{dydx}}.$$E42

and

$$p_{k,l}\ge 0,\forall \left(k,l\right)\in I_N=\left\{1,\cdots ,N\right\}\times \left\{1,\cdots ,N\right\}.$$E43

To solve numerically the frictionless contact problem, it is convenient to recall some definitions introduced for the first time by Paggi et al. [44, 52]. First, let us define the set ${\overline{I}}_c$ such as,

$${\overline{I}}_c=\left\{\left(i,j\right)\in I_N:{\zeta }_{i,j}<{\zeta }_{\max }-\delta \right\},$$E44

which can be read as the set of the elements $S_{\mathit{ij}}$ that are certainly not in contact (see the configuration (iii) in Figure 4). Therefore, the complementary set of ${\overline{I}}_c$ is simply $I_c=I_N\backslash {\overline{I}}_c$ which represents the initial trial elements $S_{\mathit{ij}}$ that are in contact. Secondly, we define the real contact set $I_c^{\ast }$ which includes the actual contact elements $S_{\mathit{ij}}$. It should be noted that, due the interaction between the asperities (see the configuration (i) in Figure 4), $I_c^{\ast }$ is a subset of $I_c$. We may have,

$$u_{i,j}>{\overline{u}}_{i,j},\text{where}{\overline{u}}_{i,j}=\delta -{\zeta }_{\max }+{\zeta }_{i,j}.$$E45

It should be noted that the condition,

$$\left(u_{i,j}-{\overline{u}}_{i,j}\right)p_{i,j}=0,\forall \left(i,j\right)\in I_c,$$E46

must hold in the last case implying non contact between the boundaries and thus, a zero contact pressure. In the opposite case, where the contact occurs, the contact pressure is strictly positive and hence, $u_{i,j}={\overline{u}}_{i,j}$.

If we denote $u_{i,j}-{\overline{u}}_{i,j}$ by ${\omega }_{i,j}$, the complementary condition in Eq. (45) will read,

$${\omega }_{i,j}{p}_{i,j}=0,\forall \left(i,j\right)\in I_c.$$E47

Thus, for all elements $S_{\mathit{ij}}$ such as $\left(i,j\right)\in I_c$, the discretized form in Eq. (45) can be reduced to the set $I_c$ as,

$${\omega }_{i,j}+{\overline{u}}_{i,j}=\sum _{\left(k,l\right)}{G}_{r_{i-k,j-l}}{p}_{k,l},\forall \left(i,j\right)\in I_c.$$E48

Using the matrix form, the frictionless contact problem can be expressed following the classic Linear Complementarity Problem (LCP) as follows,

$$\begin{array}{ll}\text{For}& \mathrm{a}\text{given}\mathrm{far}\text{field\hspace{0.17em}displacement}\delta \\ \text{find}& \omega \text{and}p\\ \text{such}\mathrm{as}& \omega =\mathit{Ap}-\overline{u}\\ \text{with}& \omega \ge 0,p\ge 0,{\omega }^{T}p=0.\end{array}$$E49

where $\omega$, $p$ and $\overline{u}$ are vectors in ${\mathbb{R}}^n$, $n$ is the number of elements of $I_c$. They represent, respectively, the unknown elastic correction, the unknown contact force and the compenetration which is related to the far field displacement $\delta$ and the distance ${\zeta }_{\max }-\zeta$. $A$ is matrix obtained through the Green-functions. Note that problem (48) is equivalent to the following convex quadratic program (QP),

$$\begin{array}{l}\underset{p}{min}\frac{1}{2}{p}^T\mathit{Ap}-\overline{u}p\\ \mathrm{s}.\mathrm{t}.p\ge 0,\end{array}$$E50

where its solution $p$ and its dual $\omega$ solve the LCP defined in Eq. (48). The opposite remains true.

Solving the discretized form of LCP is not trivial since neither the contact pressure nor the displacement are known in priory. In the literature, authors developed several numerical schemes such as (i) pivoting methods (known also as Lemke’s algorithm), (ii) Non-Negative Least Squares (NNLS) solver and (iii) the FFT-based algorithm to solve efficiently the contact problem. The pivoting-based method reaches the exact solution after a finite number of pivoting. However, when the matrix A is large, pivots number can be quite large and hence, we lost the efficiency [62]. Regarding NNLS, it takes advantage of the linear elasticity property to recast the LCP as a NNLS problem with a particular factorization. In other words, Eq. (49) will be rewritten as,

$$\begin{array}{r}\underset{p}{min}\frac{1}{2}\parallel \mathit{Cp}-C^{-T}\overline{u}{\parallel }_2^2\\ \mathrm{s}.\mathrm{t}.p\ge 0,\end{array}$$E51

where $C$ is the so-called Cholesky matrix derived through Cholesky factorization of the matrix A. We write $A=C^{T}C$. Note that, using the algorithm developed in [44], the QP in Eq. (50) is solved without an explicit computation of the matrix $C$ or its inverse $C^{-1}$. This last point makes the NNLS solver optimal and efficient to solve the fictionless contact problem of rough surfaces. In addition to the above methods, the FFT-based scheme can be an alternative solution to solve the LCP while reducing the CPU time [49]. In fact, the matrix $A$ in Eqs. (48) and (50) of size $N_x^2\times N_y^2$, with $N_i$ denotes the number of the contacting nodes following i-th direction. Unlike the sparsity property within FEM, the resulted $A$ matrix is symmetric completely full. Furthermore, performing operations with it to solve the contact problem, especially for rough surfaces with high resolution, will be heavy and will require a lot of memory and computational power. In the formulation of the FFT-based solver, the convolution theorem is applied to Eq. (38) in order to compute the displacement. First, the Discrete Fourier Transform (DFT) algorithm is used to expand the matrix $A$ as well as the pressure $p$ into Fourier’s space. Then, the inverse DFT is applied to solve Eq. (38) following an optimization scheme to compute the displacement $u$. In this chapter, the NNLS and FFT-based solver will be considered to reveal how the true contact area and the distribution of local pressures vary with the applied load, the measured and generated surface roughness.

3. Experimental rough surface

In the previous section, several theories have been reviewed with details and comments on the possible algorithmic implementation to tackle the frictionnless contact problem of two elastic bodies with rough boundaries. We have, however, not yet discussed the roughness from an experimental and numerical point of view. In this section, we thus present an AFM measurement of Au surface. Special attention will be given to the characterization of the Au measured surface following the PSDF in order to generate numerically rough surfaces similar to what has been observed experimentally. For this, we introduce first basic concepts to characterize the measured rough surface by the means of the PSDF. Then, a brief comment will be given on the comparison between the measured and the generated Au rough surface.

3.1 Theoretical background

Every surface exhibits roughness at the micro- or nano-scales. The surface variations can be fully characterized by their height probability density and power spectral density functions (i.e. HPDF and PSDF). The former holds the out-of-plane information while the latter describes the spatial arrangement. The HPDF can take various forms (i.e. exponential, Gaussian, Weibull…) depending on the fabrication and the post-processing techniques. Surfaces fabricated using a random deposition process, such as evaporation and magnetron sputtering in the microelectronics industry, generally lead to surface heights with near-Gaussian distribution.

The relevant statistical quantities that characterize the HPDF are: (i) the mean $\overline{z}$, (ii) the variance ${\sigma }^2$, (iii) the skew ($S_k$) and (iv) the kurtosis ($K_u$) of surface heights. They are defined by Eqs. (51), (52), (53) and (54), respectively.

$$\overline{z}=\frac{1}{N_p}\sum _{i=1}^{N_p}{z}_i$$E52

$${\sigma }^2=\frac{1}{N_p}\sum _{i=1}^{N_p}{\left(z_i-\overline{z}\right)}^2$$E53

$$S_k=\frac{1}{N_p{\sigma }^3}\sum _{i=1}^{N_p}{\left(z_i-\overline{z}\right)}^3$$E54

$$K_u=\frac{1}{N_p{\sigma }^4}\sum _{i=1}^{N_p}{\left(z_i-\overline{z}\right)}^4$$E55

On the other hand, the PSDF quantifies the contribution of the different wavelengths to the surface topography. It is directly calculated from the Fast Fourier Transform (FFT) of measurement data $z\left(x_k,y_l\right)$ using Eq. (55):

$$\text{PSDF}\left(k_{x_k},k_{y_l}\right)=\frac{{\mathrm{\Delta }}^2}{{\left(2\pi \right)}^2{N}_p}{\left|\mathrm{FFT}\left(z\left(x_k,y_l\right)\right)\right|}^2,$$E56

where $\mathrm{\Delta }$ refers to the sampling length, ($x_k$,$y_l$) are the discrete ($x$,$y$) coordinates and ($k_{x_k}$,$k_{y_l}$) represent the discrete spatial frequencies.

According to Persson’s work on the nature of rough surfaces [63], the PSDF follows a power law for $k\in \left[k_r,k_s\right]$ and reaches a plateau $C_0$ for $k\in \left[k_l,k_r\right]$, where $k_r$ is the roll off frequency, $k_s=2\pi /\mathrm{\Delta }$ is the highest measurable frequency and $k_l$ represents the lowest measurable frequency and is related to the scan length $L$ as $k_l=2\pi /L$. The theoretical form of the PSDF is given as follows,

$$\text{PSDF}\left(k\right)=C_0\times \left\{\begin{array}{ll}1& \text{if}{k}_l<k\le k_r\\ {\left(\frac{k}{k_r}\right)}^a& \text{if}{k}_r<k\le k_s\\ 0& \text{else}\end{array}\right.,$$E57

where the exponent $a$ can be expressed depending on the Hurst component as $a=-2\left(1+H\right)$.

3.2 Rough surface topography measurement

This section presents the characterization of a real rough surface through its HPDF and PSDF. The surface of a 1 μm-thick Au film deposited on a silicon wafer using radio frequency magnetron sputtering technique is considered for the study of a real surface topography. The AFM measurements are performed under Veeco Dimension 3100 instrument using a super sharp silicon (SSS) probe from NANOSENSORS, the scanned area is $1.5\times 1.5{\mathrm{\mu m}}^2$ with a resolution of $512$ in both directions.

Figure 5a presents the 3D view of the measured surface. The standard deviation of surface heights $\sigma$ is $2.7$ nm and the mean height $\overline{z}$ is 11 nm. The surface HPDF is presented in normalized form and compared to the normal distribution in Figure 6. It’s characterized by a skewness of $-0.08$ and a kurtosis of 2.7.

Figure 5.

3D view of (a) Au micro-surface measured using AFM and (b) numerically generated rough surface. The resolution of the real and artificial rough surface is $512$ in both directions.

Figure 6.

Normalized HDFs of AFM and artificial surfaces.

The 2D PSDF computed from the AFM data using Eq. (55) shows that the surface is highly isotropic (see Figure 7a). It can be seen from Figure 8 that the radially averaged 2D PSDF increases with the increase of $\lambda$ ($\lambda =2\pi /k$) up to a roll of wavelength ${\lambda }_r$ and then it reaches a plateau. The linear region in the PSDF plot presents two slopes, the first one is proportional to $k^{-6.5}$ (i.e. $\text{PSDF}\propto k^{-6.5}$) for $k\in \left[k_r,k_c\right]$ while the second one is proportional to $k^{-4.6}$ (i.e. $\text{PSDF}\propto k^{-4.6}$) for $k>k_c$. The second slope is probably attributed to the effect of the finite AFM tip radius that filters out wavelengths smaller than a critical wavelength ${\lambda }_c=2\pi /k_c$ [64]. Note that the PSDF becomes constant for high values of $k$ indicating a reliability cutoff, where data close this constant tail region line are affected by instrumental noise and should be discarded [64]. The magnitude of the PSDF plateau $C_0$, the roll-off wavelength ${\lambda }_r$ and the critical wavelength ${\lambda }_c$ are obtained by curve fitting of the different regions, their respective values are 3.8 × 10² nm⁴, 90 nm and 50 nm.

Figure 7.

2D PSDFs of (a) Au micro-surface measured using AFM and (b) numerically generated rough surface.

Figure 8.

Radially averaged PSDF of from AFM measurements on Au surface and the corresponding fit.

3.3 Rough surface generation

The artificial surface is aimed to mimic the behavior of Au surface. Thus, a Gaussian surface is generated using an implementation inspired by the Wu method [21] and based on the fitted PSDF of the measured Au surface (see Figure 8—red dashed line). The artificial surface is characterized by geometric and statistical quantities similar to those of Au surface, specifically (i) the standard deviation of surface heights $\sigma =2.7\mathrm{nm}$, (ii) the roll-off frequency $k_r=2\pi /{\lambda }_r$ with ${\lambda }_r=90\mathrm{nm}$, (iii) the critical frequency ${\lambda }_c=50\mathrm{nm}$ (iv) the scanning length $L=1.5\mathsf{\mu m}$, (v) the sampling length $\mathrm{\Delta }=2.93\mathrm{nm}$ and (vi) the exponent $a$ of the PSDF. The latter takes values of $-6.5$ for $k\in \left[k_r,k_c\right]$ and $-4.6$ for $k¿k_c$. The plots of the 3D view, HDF and PSDF of the generated surface are presented respectively in Figure 5b, 6 and 7b along with those from Au measured surface. Even if both surfaces are mathematically equivalent (i.e. similar HDF and PSDF), the small discrepancies on the HDF and PSDF lead to slightly different topographies, with Au surface having wider summits. However, as discussed in 4.3, the area-load mechanical behavior of both surface is similar with a maximum error of $9\%$.

4. Results and discussion

4.1 Convergence of the reference solution: BEM

The accuracy of the estimated contact area through the BEM is highly linked to the discretized rough surfaces and thus, strongly depends on the number of cells chosen to model the roughness. To highlight this problem, also known as the mesh-convergence in the literature, several computations were carried out on the $12$ generated rough surfaces using Wu’s procedure [21]. It should be noted that the artificial rough surfaces are generated using the fitted PSDF based on the AFM measurement of the real Au micro-surface (see Figure 8). The numerical discretization of the $12$ generated rough surfaces varies from $4\times 4$ to $8192\times 8192$ pixels. For each case, the normalized contact area $A/A_0\times 100$ is plotted as a function of the normalized contact pressure $P/E^{\ast }$, where $A$ refers to the real contact area, $A_0$ denotes the nominal area, $P$ represents the contact pressure and $E^{\ast }$ is the composite Young modulus. The results are briefly summarized in Figure 9 for both FFT and NNLS solvers.

Figure 9.

Effect of the mesh discretization on the true contact area—contact pressure law for a generated rough surface. (a) FFT solver. (b) NNLS solver.

One can observe from Figure 9 that the contact law, $A/A_0\times 100=f\left(P/E^{\ast }\right)$, is clearly mesh sensitive for both FFT and NNLS solvers. In fact, the predicted results for a coarse mesh, i.e. $4\times 4$, $8\times 8$, $16\times 16$ and $32\times 32$, are not consistent for the NNLS solver (see Figure 9b) and the resulting real contact area is overestimated. However, it can be seen clearly that, for the same solver, the convergence is reached with $64\times 64$ pixels. On the other hand, the FFT solver has an identical behavior with a lower convergence rate. Indeed, the convergence of the FFT solver starts from a $256\times 256$ pixel mesh which is finer than the converging mesh of the NNLS (i.e. $64\times 64$). As regard to the CPU time, no comparison can be made as the two solvers have been implemented on two different codes. However, the FFT-based solver is sufficiently fast and requires less memory. For instance, on a computing station equipped with $32\mathrm{GB}$ RAM and $12$ logical cores, the FFT-based solver manages to solve the frictionless contact problem of a generated rough surface of $512\times 512$ pixels, while the NNLS solver ran out of memory.

In the following, the predicted results using $512\times 512$ pixels by means of the FFT-based solver is considered as a reference solution to study the above theories.

4.2 Linearity between the contact area and pressure

Let us consider a generated rough surface in Figure 5b which is in contact with a flat plane. The latter moves normally towards the rough surface according to a far field displacement $\delta$. In this section, the predicted contact area using the above theories are compared and discussed in respect to the reference solution which is the BEM prediction on the finest discretization of the generated rough surface. The analysis mainly covers the evolution of the normalized true contact area as a function of the normalized contact pressure using five models,

1. The Boundary Element Method on the generated rough surface. The resolution is performed on the finest discredited rough surface which contains $512\times 512$ pixels using three codes: the first is the NNLS and it is denoted by BEM-NNLS. The second is based on the FFT. Thus, it is denoted by BEM-FFT. The third is Lars Pastewka’s code available to the public on http://contact.engineering/. The later is denoted by BEM-L. Patewska.

2. The Discrete Greenwood and Williamson model using spherical asperities without interaction. It is denoted by GW theory;

3. The modified Greenwood and Williamson model using sinusoidal asperities with interaction. It is denoted by GW-interaction;

4. The analytic solution of Persson Theory which is denoted by Persson theory;

5. The Bush-Gibson-Thomas asperity contact theory, known as BGT theory.

The predicted results are presented in Figure 10. Overall, the above theories predict the same behavior for a small pressure margin. The normalized contact area—normalized contact pressure relationship is almost linear whether the interaction is taken into account or not. For an accurate modeling, the contact area vs. contact pressure may be approximated by a power law function rather than a linear function. This first results is well known in the literature [65]. However, the coefficient of the proportionality (i.e. the slop of Figure 10) varies from one theory to another. For instance, the BGT’s slope is quite different from that predicted by the reference solution and the Persson theory. With regards to the asperity based models, the coefficient of proportionality of GW theory with or without interaction are quite similar. They vary from $1.36$ for the GW theory to 1.60 for GW-interaction model. In particular, the enhanced GW model gives similar results to those obtained by BEM using the FFT-based solver and Patewska application. Note, however, that these last coefficients are lower than the asymptotic BGT theory which predict a coefficient of 2.5. On the other hand, Persson theory prediction is lower than the asperity based models; for the same pressure, Persson’s theory predicts a smaller contact area. The same observations can be made for the BEM solution using the NNLS solver, which is not the case in the studies available in the literature. This non-conventional result can be explained by the fact that the BEM - NNLS is performed on 1 randomly generated rough surface. For this type of problem, it is more convenient to run Monte-Carlo like simulation on several realization using the same PSDF. The last point will be addressed in Section 4.3.

Figure 10.

The comparison between the BGT, GW model with and without interaction, Persson theory and BEM. Observe that the predicted laws area-contact pressure are power law function except BGT prediction which involves linear approximation. The inset show the evolution of the slop, for each prediction, as a function of the dimensionless contact pressure.

Given that all contact models presented above, the pioneering theory of GW [23] gives an accurate results for a fraction of contact area lower than 20%. The predicted true contact area are very similar to the other enhanced models. With regards to the literature, the obtained results can be analyzed using Nayak’s parameter $\alpha$ which is related to the $i$-th order moment $m_i$ of the PSDF. For our case, the measured and generated surfaces have Nayak’s parameter roughly equal to 2.2 which is close to the analysis conducted by Carbone et al. [65, 66]. For $\alpha =2$, they showed that GW theory gives reasonable predictions for a low squeezing pressure. Under the same conditions, they showed that the linearity between true contact area and contact pressure holds only for small contact area fractions. However, when the squeezing load is higher, the prediction of the above theories deviates from linearity which is in line with the findings in this contribution.

4.3 Comment on the predicted contact area of measured and generated rough surfaces

In the previous sections, all the above methods were successfully applied on numerically generated rough surfaces. Recall that these were generated using Wu’s algorithm [21]. It takes as input the fitted PSDF based on the measured one. Let us, however, consider the case of a measured Au rough surface. The objective of this section is to conduct a BEM simulation on the measured Au rough surface and compare the results with those obtained on randomly generated rough surfaces using Monte-Carlo simulation on 1000 realizations.

The results for the evolution of the real contact area in dimensionless form for Au measured rough surface are depicted in Figure 11. They are compared to BEM solution over 1000 equivalent rough surfaces in order to obtain statistical relevant results. The average prediction over 1000 realizations is also computed (see Figure 11 red dashed line). The obtained predictions of Au measured rough surface are a little bit higher than the average BEM predictions on randomly generated rough surfaces. The maximum error does not exceed 9% for the considered contact pressure margin. Note that this difference is observed despite the fact that the PSDFs of the measured and generated rough surfaces are equivalent (see Figures 6 and 8). The observed discrepancy can be attributed to (i) the measured HDF of Au surface which not a perfect Gaussian (i.e the measured kurtosis is 2.7 instead of 3), (ii) the fluctuations on the Au PSDF due to local irregularities and (iii) the AFM white noise that causes the Au PSDF to flatten at high spatial frequencies $k$ [64].

Figure 11.

True contact area vs. contact pressure in dimensionless form: comparison between BEM solution on 1000 realizations (e.g. 1000 randomly generated rough surfaces—gray dashed lines), average curve over the 1000 realization (red dashed line) and BEM solution on the measured $512\times 512$ Au rough surface using AFM technique (black dashed line).

5. Conclusion

The aim of this contribution is to present the basic guidelines to tackle the problem of contacting rough surfaces accounting for the real surface topography. For this purpose, the authors presented different theories of contact mechanics and introduced the procedure to characterize the surface topography from AFM measurements by the mean of the height distribution and power spectral density functions (HDF and PSDF).

First, a brief description has been devoted to synthesize the fundamentals of Greenwood’s theory with and without interaction, Persson’s theory and the Boundary Element Method. Next, a rough surface of length 1.5 μm, which is comparable to the contact electrode surfaces for ohmic contact micro-switches, was fabricated by sputtering 1 μm thick gold over silicon wafer and measured by AFM technique with 512 heights per side.

The measured HDF and PSDF was approximated and then, used to generate several artificial surfaces to study the mesh effect of the BEM, considered as the reference solution, by means of two solvers (i.e. FFT and NNLS). The convergence study showed that the predicted results using NNLS solver converge from an approximation of $64\times 64$ pixels, while $256\times 256$ pixels are needed for the FFT-based solver to converge. However, the latter is fast and requires less computational power compared to the NNLS solver.

The Greenwood and Williamson’s theory, with and without interaction, and Persson’s theory were applied to the generated rough surface and compared to the BEM. All these models predict approximately a power-law evolution of the contact surface as a function of dimensionless pressure. A notable difference was observed in the slope of each prediction, and the proportionality coefficients found are consistent with the literature. It is worth mentioning that the pioneering theory of GW gives accurate results for small contact area fraction and is in a good agreement with GW-interaction model and BEM-FFT predictions. However, GW results start to deviate when the pressure load increases since the interactions between asperities have a larger effect.

Finally, the Au rough surface was studied using the BEM and the results were successfully compared with intensive Monte-Carlo simulations on 1000 generated rough surface populations. The generation of artificial surfaces from the approximated HDF and PSDF succeed to catch the real surface topography and allows to estimate the mechanical behavior of contacting rough surfaces with a high level of confidence.

Acknowledgments

The first author would like to warmly thank Prof. Vladislav Yastrebov and Prof. Marco Paggi for their immense help and extremely relevant explanations during the workshops organized at Mines-Paristech France and IMT-Lucca Italy, respectively. This research work could not have been done without their help.

Conflict of interest

The authors declare no conflict of interest.