Open access peer-reviewed chapter

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

Written By

Farouk Maaboudallah, Mohamed Najah and Noureddine Atalla

Reviewed: 22 December 2021 Published: 22 February 2022

DOI: 10.5772/intechopen.102358

From the Edited Volume

Tribology of Machine Elements - Fundamentals and Applications

Edited by Giuseppe Pintaude, Tiago Cousseau and Anna Rudawska

Chapter metrics overview

306 Chapter Downloads

View Full Metrics


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.


  • contact mechanics
  • roughness
  • elastic contact
  • true contact are
  • atomic force microscopy

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 p0. 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) Δx/λs and (ii) λl/L where Δx is the distance between the mesh points, λs and λ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 106 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 δ (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 hz. 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 E1 and E2 the Young modulus of the rough surface and the rigid plane, respectively. ν1 and ν2 are the corresponding Poisson coefficient. The contact area A and load F can be expressed in terms of the rigid body motion δ as,

Figure 1.

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


where E refers to the composite Young modulus. It’s given by the following expression,


For a rigid plan, E 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 gn. This event z>gn has it own probability that can be written as,


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


The total number of asperities N can be expressed as: A0D, with A0 is the nominal area of contact and D represents the density of summits.

If we define the displacement of the asperities as ω=zgn, then the true contact area A and the expected total load F are given by,




Note that one needs to know A0, D, R and the probability density of asperities h to accurately apply the original GW model. While A0 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 m0, m2 and m4, 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:


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 m0 is the variance,


the second and fourth spectral moments can be calculated as,


with Np the total sampling points and Δ the sampling length (i.e. NΔ=L) and 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 Ri and height zi. Then, the displacement of the i-th asperity can be expressed using Macaulay bracket as,


where zM is the maximum height of the rough surface. We have zM=maxzi.

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


Note that, in the discrete GW formulation, the gap function is given by gn=zMz¯δ where 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 hi and and wavelength λi. The function gn 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 E1 and ν1, respectively. Assuming that a far field displacement δ is applied to the rigid plan. When the contact is activated, the rigid body motion δ will be accommodated by the displacement of the asperity ω and the displacement of the substrate uz. We write,


In Eq. (15), the uz 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 uz. It is given by:


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, ω, can be given through,


where E refers to the plan stress Young modulus and Ai represents the expansion coefficient of Taylor series at the vicinity of the contacting point. They depend closely on the height hi and the wavelength λ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,


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,


The additional term, uz,ij, reefers to the interaction between the i-th and j-th asperities. We have,


where rij refers to the euclidean norm between the i-th and j-th asperities. uz,ij may read as the displacement of the i-th asperity caused by the j-th contacting asperity whose contact force is Fj. 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=a1an, for each load increment δi, that minimizes the objective function Θa such as:


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 Θa 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 ai. 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 i1;

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


where Δ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,


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


The Hessian matrix is the square matrix of second-order partial derivatives of ra. We write,


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


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,


where λ 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 Δa such as


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


    Note also that the damping factor λ affect both the search direction and the step size Δa. In this research work, the trust-region (TR) method is used, after each iteration, to control and update the damping factor λ.

  2. 5. Update the solution: compute ai+1 such as


  3. 6. Evaluate the objective function: compute the residuals r at the new solution ai+1. Compute the new objective function Θnew=Θai+1 at the new solution ai+1

    If Θnew<Θai, accept the new solution ai+1. Otherwise, reject the radii step Δ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 ai+1 for the current contact step. If the procedure had not yet converged but the radii step Δa was accepted, compute the Jacobian matrix at the new solution ai+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 δ. 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 p0 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 Ppζ, where ζ is the magnification factor which is related to the roll-off wavelength (or wavenumber). We have ζ=k/k0=λ0/λ. When the roughness is added, the function Ppζ 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.




G is a function of the magnification. We write,


from which its follows,


k0 is the wavenumber related to the roll-off distance and the quantity πk0ζk0k3PSDFkdk is the second profile PSDF moment, m2, 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-rd boundary condition Pζ=0 should be replaced by,


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 Ei and the Poisson ratio ν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 δ is imposed. The profile of the rough surface is characterized by the function ζx measured with respect to the mean plane. The peak of the rough surface is the ζmax such as ζmax=maxxSζx. When the asperity comes in contact with the elastic half-plane, a displacement 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 ζmaxζ with ζ, in this case, denotes the local peak of the contacting asperity. Otherwise, a generic displacement of the elastic half-plane is simply denoted by ux.

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 δ.

In the fundamental BEM, the normal displacement field, ux, at the coordinate x is related to the contact pressure py at the coordinate y by the following integral form:


where Grxy is the Green-function which links the displacement ux at x to the contact pressure py acting at y. It is given by the Boussinesq solution as follows:


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

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

Foragivenfarfield displacementδfindxSuxandpxsuchasux=S1π1ν12E1+1ν22E21xypydywithuxu¯xδ0px0uxu¯xδpx=0E40

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

  1. the barycentric coordinate: xi,j=1SijxSijxdA;

  2. the average height: ζi,j=1SijxSijζxdA;

  3. the resultant contact load: pi,j=xSijpxdA;

  4. the displacement: ui,j=1SijxSijuxdA.

Note that Sij refers to a cell of the surface S. The latter can be seen as a square grid spacing by Δ and containing N×N points. Therefore, the cell Sij is a square surface of dimension Δ2.

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


where Grik,jl represents the averaged Green function over the cell Sij. We write,




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 I¯c such as,


which can be read as the set of the elements Sij that are certainly not in contact (see the configuration (iii) in Figure 4). Therefore, the complementary set of I¯c is simply Ic=IN\I¯c which represents the initial trial elements Sij that are in contact. Secondly, we define the real contact set Ic which includes the actual contact elements Sij. It should be noted that, due the interaction between the asperities (see the configuration (i) in Figure 4), Ic is a subset of Ic. We may have,


It should be noted that the condition,


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, ui,j=u¯i,j.

If we denote ui,ju¯i,j by ωi,j, the complementary condition in Eq. (45) will read,


Thus, for all elements Sij such as ijIc, the discretized form in Eq. (45) can be reduced to the set Ic as,


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

Foragivenfarfield displacementδfindωandpsuchasω=Apu¯withω0,p0,ωTp=0.E49

where ω, p and u¯ are vectors in Rn, n is the number of elements of Ic. They represent, respectively, the unknown elastic correction, the unknown contact force and the compenetration which is related to the far field displacement δ and the distance ζmaxζ. A is matrix obtained through the Green-functions. Note that problem (48) is equivalent to the following convex quadratic program (QP),


where its solution p and its dual ω 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,


where C is the so-called Cholesky matrix derived through Cholesky factorization of the matrix A. We write A=CTC. 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 C1. 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 Nx2×Ny2, with Ni 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 z¯, (ii) the variance σ2, (iii) the skew (Sk) and (iv) the kurtosis (Ku) of surface heights. They are defined by Eqs. (51), (52), (53) and (54), respectively.


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 zxkyl using Eq. (55):


where Δ refers to the sampling length, (xk,yl) are the discrete (x,y) coordinates and (kxk,kyl) represent the discrete spatial frequencies.

According to Persson’s work on the nature of rough surfaces [63], the PSDF follows a power law for kkrks and reaches a plateau C0 for kklkr, where kr is the roll off frequency, ks=2π/Δ is the highest measurable frequency and kl represents the lowest measurable frequency and is related to the scan length L as kl=2π/L. The theoretical form of the PSDF is given as follows,


where the exponent a can be expressed depending on the Hurst component as a=21+H.

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×1.5μm2 with a resolution of 512 in both directions.

Figure 5a presents the 3D view of the measured surface. The standard deviation of surface heights σ is 2.7 nm and the mean height 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 λ (λ=2π/k) up to a roll of wavelength λr and then it reaches a plateau. The linear region in the PSDF plot presents two slopes, the first one is proportional to k6.5 (i.e. PSDFk6.5) for kkrkc while the second one is proportional to k4.6 (i.e. PSDFk4.6) for k>kc. The second slope is probably attributed to the effect of the finite AFM tip radius that filters out wavelengths smaller than a critical wavelength λc=2π/kc [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 C0, the roll-off wavelength λr and the critical wavelength λc are obtained by curve fitting of the different regions, their respective values are 3.8 × 102 nm4, 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 σ=2.7nm, (ii) the roll-off frequency kr=2π/λr with λr=90nm, (iii) the critical frequency λc=50nm (iv) the scanning length L=1.5μm, (v) the sampling length Δ=2.93nm and (vi) the exponent a of the PSDF. The latter takes values of 6.5 for kkrkc and 4.6 for k¿kc. 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×4 to 8192×8192 pixels. For each case, the normalized contact area A/A0×100 is plotted as a function of the normalized contact pressure P/E, where A refers to the real contact area, A0 denotes the nominal area, P represents the contact pressure and E 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/A0×100=fP/E, is clearly mesh sensitive for both FFT and NNLS solvers. In fact, the predicted results for a coarse mesh, i.e. 4×4, 8×8, 16×16 and 32×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×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×256 pixel mesh which is finer than the converging mesh of the NNLS (i.e. 64×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 32GB RAM and 12 logical cores, the FFT-based solver manages to solve the frictionless contact problem of a generated rough surface of 512×512 pixels, while the NNLS solver ran out of memory.

In the following, the predicted results using 512×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 δ. 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×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 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 α which is related to the i-th order moment mi 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 α=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×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×64 pixels, while 256×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.



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.


  1. 1. Panagouli OK, Mastrodimou K. Dependence of friction coefficient on the resolution of asperities in metallic rough surfaces under cyclic loading. International Journal of Solids and Structures. 2017;108:85-97
  2. 2. Wang H, Zhou C, Wang H, Hu B, Liu Z. A novel contact model for rough surfaces using piecewise linear interpolation and its application in gear wear. Wear. 15 July 2021;476:203685
  3. 3. Panagouli OK, Margaronis K, Tsotoulidou V. A multiscale model for thermal contact conductance of rough surfaces under low applied pressure. International Journal of Solids and Structures. 2020;200-201:106-118
  4. 4. Zhai C, Hanaor D, Proust G, Brassart L, Gan Y. Interfacial electro-mechanical behaviour at rough surfaces. Extreme Mechanics Letters. 2016;9:422-429
  5. 5. Yastrebov VA, Durand J, Proudhon H, Cailletaud G. Rough surface contact analysis by means of the finite element method and of a new reduced model. Comptes Rendus Mcanique. 2011;339:473-490
  6. 6. Magnier V, Brunel JF, Dufrnoy P. Impact of contact stiffness heterogeneities on friction-induced vibration. International Journal of Solids and Structures. 2014;51:1662-1669
  7. 7. Maaboudallah F, Atalla N. An efficient numerical strategy to predict the dynamic instabilities of a rubbing system: Application to an automobile disc brake system. Computational Mechanics. 31 March 2021;67:1465-1483
  8. 8. Rogozhin A, Miakonkikh A, Tatarintsev A, Lebedev K, Kalnov V, Rudenko K. Silicon Ohmic Lateral-Contact MEMS Switch for RF Applications. Proceedings Vol. 10224. Zvenigorod, Russian Federation: International Conference on Micro- and Nano-Electronics. 2016. pp. 1-7. DOI: 10.1117/12.2267093
  9. 9. Qian J-y, Hou C-w, Li X-j, Jin Z-j. Actuation mechanism of microvalves: A review. Micromachines (Basel) – MDPI. 2020;172(11):2. Available from:
  10. 10. Schmidt MA. “SiliconWafer Bonding for Micro mechanical Devices”. Hilton Head Island, South Carolina: Solid-State Sensors, Actuators, and Microsystems Workshop; 12-16 June 1994. pp. 127–131. DOI: 10.31438/trf.hh1994.30
  11. 11. Manneborg A, Nese M, Bhlckers P. “Silicon-to-silicon anodic bonding with a borosilicate glass layer.” Journal of Micromechanics and Microengineering. 1991;1:139
  12. 12. Gui C, Albers H, Gardeniers JGE, Elwenspoek M, and Lambeck PV. “Fusion bonding of rough surfaces with polishing technique for silicon micromachining”. Microsystem Technologies. 1997;3:122–128
  13. 13. Ma Q, Tran Q, Chou T-KA, Heck J, Bar H, Kant R, et al. Metal contact reliability of RF MEMS switches. In: Reliability, Packaging, Testing, and Characterization of MEMS/MOEMS VI. Proceedings. Vol. 6463. San Jose, California, United States: MOEMS-MEMS 2007 Micro and Nanofabrication. 2007. p. 646305. DOI: 10.1117/12.702177
  14. 14. Zhao Y-P, Wang L, Yu T. Mechanics of adhesion in MEMS—a review. Journal of Adhesion Science and Technology. 2003;17(4):519-546. Cited by 366
  15. 15. Nguyen N, Schubert S, Richter S, Dtzel W. Hybrid-assembled micro dosing system using silicon-based micropump/valve and mass flow sensor. Sensors and Actuators A: Physical. 1998;69(1):85-91
  16. 16. Hayamizu S, Higashino K, Fujii Y, Sando Y, Yamamoto K. Development of a bi-directional valve-less silicon micro pump controlled by driving waveform. Sensors and Actuators A: Physical. 2003;103(1):83-87
  17. 17. Sensortechnics uk. Silicon micro valve technology. Vacuum. 1994;45(12):I
  18. 18. Gong Y, Misture ST, Gao P, Mellott NP. Surface roughness measurements using power spectrum density analysis with enhanced spatial correlation length. Journal of Physical Chemistry C. 2016;120(39):22358-22364
  19. 19. Panda S, Panzade A, Sarangi M, Roy Chowdhury SK. Spectral approach on multiscale roughness characterization of nominally rough surfaces. Journal of Tribology. 2017;139(3):1-10
  20. 20. Chen JQ, Huang QS, Qi RZ, Feng YF, Feng JT, Zhang Z, et al. Effects of sputtering power and annealing temperature on surface roughness of gold films for high-reflectivity synchrotron radiation mirrors. Nuclear Science and Techniques. 2019;30(7):1-6
  21. 21. Wu JJ. Simulation of rough surfaces with FFT. Tribology International. 2000;33(1):47-58
  22. 22. Pérez-Ràfols F, Almqvist A. Generating randomly rough surfaces with given height probability distribution and power spectrum. Tribology International. 2019;131(October 2018):591-604
  23. 23. Greenwood JA, Williamson JBP. Contact of nominally flat surfaces. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences. 1966;295(1442):300-319
  24. 24. Bush AW, Gibson RD, Thomas TR. The elastic contact of a rough surface. Wear. 1975;35:87-111
  25. 25. Nayak PR. Random process model of rough surfaces. Journal of Lubrication Technology. 1971;93:398-407
  26. 26. Miki BB. Thermal contact conductance; theoretical considerations. International Journal of Heat and Mass Transfer. 1974;17:205-214
  27. 27. Cooper MG, Mikic BB, Yovanovich MM. Thermal contact conductance. International Journal of Heat and Mass Transfer. 1969;12:279-300
  28. 28. Sridhar MR, Yovanovich MM. Review of elastic and plastic contact conductance models—Comparison with experiment. Journal of Thermophysics and Heat Transfer. 1994;8(4):633-640
  29. 29. Berthe D, Vergne P. An elastic approach to rough contact with asperity interactions. Wear. 1987;117:211-222
  30. 30. Yongsheng L, Guiping Y, Yan H, Linqing Z. The characteristics of elastically contacting ideal rough surfaces. Journal of Tribology. 1996;118:90-97
  31. 31. O’Callaghan PW, Probert SD. Real area of contact between a rough surface and a softer optically flat surface. Journal of Mechanical Engineering Science. 1970;12:259-267
  32. 32. Hendriks CP, Visscher M. Accurate real area of contact measurements on polyurethane. Journal of Tribology. 1995;117:607-611
  33. 33. Paggi M, Ciavarella M. The coefficient of proportionality k between real contact area and load, with new asperity models. Wear. March 2010;268(7):1020-1029
  34. 34. Ciavarella M, Greenwood JA, Paggi M. Inclusion of interaction in the Greenwood and Williamson contact theory;265(5):729-734
  35. 35. Afferrante L, Carbone G, Demelio G. Interacting and coalescing hertzian asperities: A new multiasperity contact model. 278-279:28-33
  36. 36. Yastrebov VA, Anciaux G, Molinari J-F. From infinitesimal to full contact between rough surfaces: Evolution of the contact area. International Journal of Solids and Structures. 1 January 2015;52:83-102
  37. 37. Persson BNJ. Relation between interfacial separation and load: A general theory of contact mechanics. Physical Review Letters. 2007;99(12):1-4
  38. 38. Persson BNJ, Albohr O, Tartaglino U, Volokitin AI, Tosatti E. On the nature of surface roughness with application to contact mechanics, sealing, rubber friction and adhesion. Journal of Physics: Condensed Matter. 2004;17:R1-R62
  39. 39. Persson BNJ. Contact mechanics for randomly rough surfaces. Surface Science Reports. 2006;61:201-227
  40. 40. Manners W, Greenwood JA. Some observations on Persson’s diffusion theory of elastic contact. Wear. 2006;261:600-610
  41. 41. Hyun S, Pei L, Molinari J-F, Robbins MO. Finite-element analysis of contact between elastic self-affine surfaces. Physical Review E. 2004;70:026117
  42. 42. Polonsky IA, Keer LM. A numerical method for solving rough contact problems based on the multi-level multi-summation and conjugate gradient techniques. Wear. July 1999;231(2):206-219
  43. 43. Stanley HM, Kato T. An FFT-based method for rough surface contact. Journal of Tribology. 1997;119:481-485
  44. 44. Bemporad A, Paggi M. Optimization algorithms for the solution of the frictionless normal contact between rough surfaces. International Journal of Solids and Structures. 2015;69-70:94-105
  45. 45. Brandt A, Lubrecht AA. Multilevel matrix multiplication and fast solution of integral equations. Journal of Computational Physics. 1990;90:348-370
  46. 46. Venner CH, Lubrecht AA. Numerical analysis of the influence of waviness on the film thickness of a circular EHL contact. Journal of Tribology. 1996;118:153-161
  47. 47. Lubrecht AA, Ioannides E. A fast solution of the dry contact problem and the associated sub-surface stress field, using multilevel techniques. Journal of Tribology. 1991;113:128-133
  48. 48. Nogi T, Kato T. Influence of a hard surface layer on the limit of elastic contact part I: Analysis using a real surface model. Journal of Tribology. 1997;119:493-500
  49. 49. Polonsky IA, Keer LM. A fast and accurate method for numerical analysis of elastic layered contacts. Journal of Tribology. 1999;122:30-35
  50. 50. Popov VL, Pohrt R, Li Q. Strength of adhesive contacts: Influence of contact geometry and material gradients. Friction. 2017;5:308-325
  51. 51. Li Q, Argatov I, Popov VL. Onset of detachment in adhesive contact of an elastic half-space and flat-ended punches with non-circular shape: Analytic estimates and comparison with numeric analysis. Journal of Physics D: Applied Physics. 2018;51:145601
  52. 52. Paggi M, Bemporad A, Reinoso J. Computational methods for contact problems with roughness. In: Paggi M, Hills D, editors. Modeling and Simulation of Tribological Problems in Technology CISM International Centre for Mechanical Sciences. Cham: Springer International Publishing; 2020. pp. 131-178
  53. 53. Frrot L, Bonnet M, Molinari J-F, Anciaux G. A Fourier-accelerated volume integral method for elastoplastic contact. Computer Methods in Applied Mechanics and Engineering. 2019;351:951-976
  54. 54. Seitz A, Wall WA, Popp A. Nitsches method for finite deformation thermomechanical contact problems. Computational Mechanics. 2019;63:1091-1110
  55. 55. Seitz A, Wall WA, Popp A. A computational approach for thermo-elasto-plastic frictional contact based on a monolithic formulation using non-smooth nonlinear complementarity functions. Advanced Modeling and Simulation in Engineering Sciences. 2018;5:5
  56. 56. McCool JI. Finite difference spectral moment estimation for profiles the effect of sample spacing and quantization error. Precision Engineering. 1982;4(4):181-184
  57. 57. Poon CY, Bhushan B. Comparison of surface roughness measurements by stylus profiler, AFM and non-contact optical profiler. Wear. 1995;190:76-88
  58. 58. Zhang S, Song H, Sandfeld S, Liu X, Wei YG. Discrete Greenwood Williamson modeling of rough surface contact accounting for three-dimensional sinusoidal asperities and asperity interaction. Journal of Tribology. 2019;141(12):121401 (11 pages). DOI: 10.1115/1.4044635
  59. 59. Ciavarella M, Delfine V, Demelio G. A “re-vitalized” Greenwood and Williamson model of elastic contact between fractal surfaces. Journal of the Mechanics and Physics of Solids. 2006;54(12):2569-2591
  60. 60. Johnson K. Contact Mechanics. Cambridge: Cambridge University Press; June 2021, 1985. ISBN: 9781139171731. DOI: 10.1017/CBO9781139171731
  61. 61. Leroux J, Fulleringer B, Nlias D. Contact analysis in presence of spherical inhomogeneities within a half-space. International Journal of Solids and Structures. 2010;47:3034-3049
  62. 62. Cottle RW, Pang J-S, Stone RE. The Linear Complementarity Problem. Society for Industrial and Applied Mathematics. Classics in Applied Mathematics 2009. ISBN:978-0-89871-686-3. eISBN:978-0-89871-900-0. DOI: 10.1137/1.9780898719000
  63. 63. Persson BN, Albohr O, Tartaglino U, Volokitin AI, Tosatti E. On the nature of surface roughness with application to contact mechanics, sealing, rubber friction and adhesion. Journal of Physics: Condensed Matter. 2005;17(1):R1. Available from:
  64. 64. Jacobs TD, Junge T, Pastewka L. Quantitative characterization of surface topography using spectral analysis. Surface Topography: Metrology and Properties. 2017;5(5):013001. Available from:
  65. 65. Carbone G, Bottiglione F. Asperity contact theories: Do they predict linearity between contact area and load? Journal of the Mechanics and Physics of Solids. 2008;56:2555-2572
  66. 66. Carbone G, Bottiglione F. Contact mechanics of rough surfaces: A comparison between theories. Meccanica. 2011;46:557-565

Written By

Farouk Maaboudallah, Mohamed Najah and Noureddine Atalla

Reviewed: 22 December 2021 Published: 22 February 2022