**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 *h z*ð Þ. 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

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

*A Review on the Contact Mechanics Modeling of Rough Surfaces in the Elastic Regime… DOI: http://dx.doi.org/10.5772/intechopen.102358*

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*<sup>1</sup> and *E*<sup>2</sup> the Young modulus of the rough surface and the rigid plane, respectively. *ν*<sup>1</sup> and *ν*<sup>2</sup> are the corresponding Poisson coefficient. The contact area *A* and load *F* can be expressed in terms of the rigid body motion *δ* as,

$$A = \pi R \delta,\tag{1}$$

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

where *E*<sup>∗</sup> refers to the composite Young modulus. It's given by the following expression,

$$\frac{1}{E^\*} = \frac{1-\nu\_1^2}{E\_1} + \frac{1-\nu\_2^2}{E\_2}.\tag{3}$$

For a rigid plan, *E*<sup>∗</sup> 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,

$$\mathcal{P}(\{z > \mathbf{g}\_n\}) = \int\_{\mathbf{g}\_n}^{+\infty} h(z) \, dz. \tag{4}$$

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

$$m = N \int\_{\mathcal{g}\_u}^{+\infty} h(z) \, dz. \tag{5}$$

The total number of asperities *N* can be expressed as: *A*0*D*, with *A*<sup>0</sup> is the nominal area of contact and *D* represents the density of summits.

If we define the displacement of the asperities as *ω* ¼ *z* � *gn*, then the true contact area *A* and the expected total load *F* are given by,

$$A = \pi \text{NR} \int\_{\mathcal{g}\_n}^{+\infty} (z - \mathbf{g}\_n) h(z) \, dz,\tag{6}$$

and

$$F = \frac{4}{3} N E^\* \, R^{1/2} \Big|\_{\mathbf{g}\_n}^{+\infty} (\mathbf{z} - \mathbf{g}\_n)^{3/2} h(\mathbf{z}) \, d\mathbf{z}.\tag{7}$$

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*<sup>0</sup> 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*<sup>2</sup> 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:

*Tribology of Machine Elements - Fundamentals and Applications*

$$m\_i = \int\_{-\infty}^{+\infty} q^i \text{PSD} \mathbf{F}(q) dq = \frac{1}{L} \int\_0^L \left(\frac{d^i x}{dx^i}\right)^2 dx,\tag{8}$$

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*<sup>0</sup> is the variance,

$$m\_0 = \frac{1}{N\_p - 1} \sum\_{i=1}^{N\_p} \left( z\_i - \overline{z} \right)^2,\tag{9}$$

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}{\Delta} \right)^2,\tag{10}$$

$$m\_4 = \frac{1}{N\_p - 2} \sum\_{i=1}^{N\_p - 2} \left( \frac{z\_{i+2} - 2z\_{i+1} + z\_i}{\Delta} \right)^2,\tag{11}$$

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,

$$
\alpha\_i = <\delta - (z\_M - z\_i) > ,\tag{12}
$$

where *zM* is the maximum height of the rough surface. We have *zM* ¼ max ð Þ *zi* . The true contact area *A* as well as the expected total load *F* are given by,

$$A = \pi \sum\_{i=1}^{N} R\_i o\_i,\tag{13}$$

$$F = \frac{4}{3} E^\* \sum\_{i=1}^{N} R\_i^{1/2} \alpha\_i^{3/2}. \tag{14}$$

Note that, in the discrete GW formulation, the gap function is given by *gn* ¼ ð Þ� *zM* � *z δ* 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 *A Review on the Contact Mechanics Modeling of Rough Surfaces in the Elastic Regime… DOI: http://dx.doi.org/10.5772/intechopen.102358*

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**).

#### *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*<sup>1</sup> 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,

$$
\delta = \alpha + \mathfrak{u}\_x. \tag{15}
$$

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:

$$
\mu\_x = \frac{1 - \nu\_1^2}{\pi E\_1 a} F,\tag{16}
$$

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,

$$F = \sum\_{i=1}^{n} \frac{4A\_i E^\* \operatorname{i} a^{2i+1} \mathbf{2} \cdot 4 \cdots \mathbf{2i}}{(2i+1)\mathbf{1} \cdot \mathbf{3} \cdots (2i-1)},\tag{17}$$

$$\rho = \sum\_{i=1}^{n} \frac{2 \cdot 4 \cdots 2i}{1 \cdot 3 \cdots (2i - 1)} A\_i a^{2i},\tag{18}$$

#### **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.*

where *E*<sup>∗</sup> 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 *λ<sup>i</sup>* 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^\* \, a^3 + \frac{64}{15} A\_2 E^\* \, a^5,\tag{19}$$

$$
\omega = 2A\_1 a^2 + \frac{8}{3} A\_2 a^4. \tag{20}
$$

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 = a\_i + u\_{x,i} + u\_{x,i}^j. \tag{21}
$$

The additional term, *u <sup>j</sup> z*,*i* , reefers to the interaction between the *i*-th and *j*-th asperities. We have,

$$\mu\_{x,i}^{j} = \sum\_{j=1}^{n} \frac{\mathbf{1} - \nu\_1^2}{\pi E\_1 r\_{ij}} F\_j,\tag{22}$$

where *rij* refers to the euclidean norm between the *i*-th and *j*-th asperities. *u <sup>j</sup> z*,*i* may read as the displacement of the *i*-th asperity caused by the *j*-th contacting asperity whose contact force is *F <sup>j</sup>*. 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* ¼ ð Þ *a*1, ⋯, *an* , for each load increment *δi*, that minimizes the objective function Θð Þ *a* such as:

$$\begin{cases} \Theta(a) = \frac{1}{2} \sum\_{i=1}^{n} \left( r\_i \right)^2 = \frac{1}{2} \| r(a) \|\_{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\_{j=1 \atop j \neq i}^{n} \frac{1 - \nu\_1^2}{\pi E\_1 r\_{ij}} F\_j \\\\ F\_i = \frac{8}{3} A\_1 E^\* a\_i^3 + \frac{64}{15} A\_2 E^\* a\_i^5 \\\ F\_j = \frac{8}{3} A\_1 E^\* a\_j^3 + \frac{64}{15} A\_2 E^\* a\_j^5 \end{cases} \tag{23}$$

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:

*A Review on the Contact Mechanics Modeling of Rough Surfaces in the Elastic Regime… DOI: http://dx.doi.org/10.5772/intechopen.102358*

For each contact step do,


$$
\Theta(a\_{i+1}) = \Theta(a\_i) + \nabla \Theta(a\_i)^T \cdot \Delta a + \frac{1}{2} \Delta a^T \cdot \mathbf{H} \cdot \Delta a,\tag{24}
$$

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,

$$\nabla\Theta(a) = \mathbf{J}(a)^T r(a),\tag{25}$$

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

$$\mathbf{J}(\boldsymbol{a}) = \nabla r(\boldsymbol{a})^T = \begin{pmatrix} \nabla r\_1^T \\ \nabla r\_2^T \\ \vdots \\ \nabla r\_n^T \end{pmatrix}. \tag{26}$$

The Hessian matrix is the square matrix of second-order partial derivatives of *r a*ð Þ. We write,

$$\mathbf{H} = \begin{pmatrix} \frac{\partial^2 \Theta}{\partial a\_1^2} & \frac{\partial^2 \Theta}{\partial a\_1 \partial a\_2} & \cdots & \frac{\partial^2 \Theta}{\partial a\_1 \partial a\_n} \\ \frac{\partial^2 \Theta}{\partial a\_2 \partial a\_1} & \frac{\partial^2 \Theta}{\partial a\_2^2} & \cdots & \frac{\partial^2 \Theta}{\partial a\_2 \partial a\_n} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^2 \Theta}{\partial a\_n \partial a\_1} & \frac{\partial^2 \Theta}{\partial a\_n \partial a\_2} & \cdots & \frac{\partial^2 \Theta}{\partial a\_n^2} \end{pmatrix}. \tag{27}$$

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

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

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},\tag{29}$$

where *λ* is the damping factor and *I* is the identity matrix.

4. **Compute the radii step:** the goal of this step is to find the direction Δ*a* such as

$$(\mathbf{J}^T \mathbf{J} + \lambda \mathbf{I}) \Delta \mathbf{a} = -\mathbf{2} \mathbf{J}^T r\_i. \tag{30}$$

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

$$\min\_{\Delta a} \frac{1}{2} \left\| \begin{bmatrix} \mathbf{J} \\ \lambda \mathbf{I} \end{bmatrix} \Delta a + \begin{bmatrix} r \\ \mathbf{0} \end{bmatrix} \right\|\_{2}^{2}. \tag{31}$$

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 *λ*.

5. **Update the solution:** compute *ai*þ<sup>1</sup> such as

$$
\mathfrak{a}\_{i+1} = \mathfrak{a}\_i + \Delta \mathfrak{a}.\tag{32}
$$

6. **Evaluate the objective function:** compute the residuals *r* at the new solution *ai*þ1. Compute the new objective function Θ*new* ¼ Θð Þ *ai*þ<sup>1</sup> at the new solution *ai*þ<sup>1</sup>

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.

7. **Convergence verification:** Compute the error. If the procedure has converged, return the contact radii *ai*þ<sup>1</sup> 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 *p*<sup>0</sup> is assumed uniform. Following the self-affine behavior carried by the PSDF, the roughness is progressively added which

**Figure 3.**

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

*A Review on the Contact Mechanics Modeling of Rough Surfaces in the Elastic Regime… DOI: http://dx.doi.org/10.5772/intechopen.102358*

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 p*ð Þ , *ζ* , where *ζ* is the magnification factor which is related to the roll-off wavelength (or wavenumber). We have *ζ* ¼ *k=k*<sup>0</sup> ¼ *λ*0*=λ*. When the roughness is added, the function *P p*ð Þ , *ζ* verifies the following diffusion equation,

$$\begin{cases} \frac{\partial P(p,\zeta)}{\partial \zeta} = f(\zeta) \frac{\partial^2 P(p,\zeta)}{\partial p^2} \\ P(p,1) = \delta(p - p\_0) \\ P(-\infty, \zeta) = \mathbf{0} \\ P(+\infty, \zeta) = \mathbf{0} \end{cases} \tag{33}$$

where,

$$f(\zeta) = G'(\zeta)p\_0^2. \tag{34}$$

*G*<sup>0</sup> is a function of the magnification. We write,

$$\mathcal{G}(\zeta) = \frac{\pi}{4} \left(\frac{E^\*}{p\_0}\right)^2 \int\_{k\_0}^{\zeta k\_0} k^3 \text{PSD} \,\mathrm{F}(k) \,d\boldsymbol{k},\tag{35}$$

from which its follows,

$$f(\zeta) = \frac{\pi}{4} E^{\*^2} k\_0 k^3 \text{PSDF}(k). \tag{36}$$

*k*<sup>0</sup> is the wavenumber related to the roll-off distance and the quantity *π* Ð *<sup>ζ</sup>k*<sup>0</sup> *<sup>k</sup>*<sup>0</sup> *<sup>k</sup>*<sup>3</sup> PSDFð Þ*<sup>k</sup> dk* is the second profile PSDF moment, *<sup>m</sup>*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-*rd* boundary condition *P*ð Þ¼ �∞, *ζ* 0 should be replaced by,

$$P(\mathbf{0}, \zeta) = \mathbf{0}.\tag{37}$$

#### **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 *ν<sup>i</sup>* 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 ¼ max *<sup>x</sup>*∈*S*ð Þ *ζ*ð Þ *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 *u x*ð Þ.

In the fundamental BEM, the normal displacement field, *u x*ð Þ, at the coordinate *x* is related to the contact pressure *p y*ð Þ at the coordinate *y* by the following integral form:

$$u(\mathbf{x}) = \int\_{\mathcal{S}} G\_r(\mathbf{x}, \mathbf{y}) p(\mathbf{y}) d\mathbf{y},\tag{38}$$

where *Gr*ð Þ *x*, *y* is the Green-function which links the displacement *u x*ð Þ at *x* to the contact pressure *p y*ð Þ acting at y. It is given by the Boussinesq solution as follows:

$$G\_r(\mathbf{x}, \boldsymbol{y}) = \frac{1}{\pi} \left( \frac{\mathbf{1} - \nu\_1^2}{E\_1} + \frac{\mathbf{1} - \nu\_2^2}{E\_2} \right) \frac{\mathbf{1}}{r},\tag{39}$$

where *r* is the distance between the points *x* and *y*. It refers to the standard Euclidean norm. We have, *r* ¼ ∥*x* � *y*∥.

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

$$\begin{array}{ll}\textbf{For} & \textbf{a given}\,\textbf{far}\,\textbf{field}\,\textbf{displacement}\,\delta\\ \displaystyle\underset{\textbf{x}\in S}{\text{find}} & u(\textbf{x})\,\textbf{and}\,p(\textbf{x})\\ \textbf{such as}\,\textbf{a} & u(\textbf{x}) = \int\_{S} \frac{1}{\pi} \left(\frac{1-\nu\_{1}^{2}}{E\_{1}} + \frac{1-\nu\_{2}^{2}}{E\_{2}}\right) \frac{1}{\|\textbf{x}-\textbf{y}\|} p(\textbf{y}) \,d\textbf{y} \\ \textbf{with} & \begin{cases} u(\textbf{x}) - \overline{u}(\textbf{x},\delta) \ge 0\\ p(\textbf{x}) \ge 0\\ (u(\textbf{x}) - \overline{u}(\textbf{x},\delta)) p(\textbf{x}) = \textbf{0} \end{cases} \end{array} \tag{40}$$

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

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*,*<sup>j</sup>* <sup>¼</sup> <sup>1</sup> *Sij* Ð *x*∈*Sij x dA*;

2. the average height: *<sup>ζ</sup>i*,*<sup>j</sup>* <sup>¼</sup> <sup>1</sup> *Sij* Ð *x*∈*Sij ζ*ð Þ *x dA*;

3. the resultant contact load: *pi*,*<sup>j</sup>* <sup>¼</sup> <sup>Ð</sup> *x*∈*Sij p x*ð Þ *dA*;

4. the displacement: *ui*,*<sup>j</sup>* <sup>¼</sup> <sup>1</sup> *Sij* Ð *x*∈*Sij u x*ð Þ *dA*.

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 Δ<sup>2</sup> .

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

$$\mu\_{i,j} = \sum\_{k=1}^{N} \sum\_{l=1}^{N} \mathbf{G}\_{r\_{i-kj-l}} p\_{k,l}, \ \forall (i,j) \in I\_N = \{1, \dots, N\} \times \{1, \dots, N\} \tag{41}$$

where *Gri*�*k*,*j*�*<sup>l</sup>* represents the averaged Green function over the cell *Sij*. We write,

$$G\_{r\_{i-kj-l}} = \frac{1}{\Delta^2} \int\_{S\_{\vec{y}}} \int\_{S\_{kl}} G\_r(\mathbf{x}, y) \, dy d\mathbf{x}.\tag{42}$$

and

$$p\_{k,l} \ge 0, \ \forall (k,l) \in I\_N = \{1, \cdots, N\} \times \{1, \cdots, N\}. \tag{43}$$

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 *Ic* such as,

$$\overline{I}\_{\varepsilon} = \left\{ (i, j) \in I\_N : \zeta\_{i\circ} < \zeta\_{\max} - \delta \right\},\tag{44}$$

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 *Ic* is simply *Ic* ¼ *IN*n*Ic* which represents the initial trial elements *Sij* that are in contact. Secondly, we define the real contact set *I* <sup>∗</sup> *<sup>c</sup>* 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**), *I* <sup>∗</sup> *<sup>c</sup>* is a subset of *Ic*. We may have,

$$
\overline{u}\_{i\circ} > \overline{u}\_{i\circ}, \quad \text{where} \quad \overline{u}\_{i\circ} = \delta - \zeta\_{\text{max}} + \zeta\_{i\circ}.\tag{45}
$$

It should be noted that the condition,

$$(u\_{i,j} - \overline{u}\_{i,j})p\_{i,j} = \mathbf{0}, \ \forall (i,j) \in I\_c,\tag{46}$$

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*,*<sup>j</sup>* ¼ *ui*,*<sup>j</sup>*.

If we denote *ui*,*<sup>j</sup>* � *ui*,*<sup>j</sup>* by *ωi*,*j*, the complementary condition in Eq. (45) will read,

$$
\alpha\_{i,j} p\_{i,j} = 0, \ \forall (i,j) \in I\_c. \tag{47}
$$

Thus, for all elements *Sij* such as ð Þ *i*, *j* ∈*Ic*, the discretized form in Eq. (45) can be reduced to the set *Ic* as,

$$
\rho\_{i,j} + \overline{u}\_{i,j} = \sum\_{(k,l)} G\_{r\_{i-kj-l}} p\_{k,l}, \ \forall (i,j) \in I\_c. \tag{48}
$$

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} & \text{a given} \text{farfeld displacement} \,\delta\\ \text{find} & \omega \,\text{and} \, p\\ \text{such as} & \omega = Ap - \overline{u} \\ \text{with} & \omega \ge 0, p \ge 0, \,\omega^T p = 0. \end{array} \tag{49}$$

where *ω*, *p* and *u* are vectors in *<sup>n</sup>*, *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),

$$\begin{aligned} \min\_{p} & \frac{1}{2} p^T A p - \overline{u} p \\ & \text{s.t.} \, p \ge 0, \end{aligned} \tag{50}$$

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,

$$\begin{aligned} \min\_{\overline{p}} & \frac{1}{2} \|Cp - C^{-T}\overline{u}\|\_{2}^{2} \\ & \text{s.t.} \, p \ge 0, \end{aligned} \tag{51}$$

where *C* is the so-called Cholesky matrix derived through Cholesky factorization of the matrix *<sup>A</sup>*. We write *<sup>A</sup>* <sup>¼</sup> *<sup>C</sup>TC*. 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*�<sup>1</sup> . 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*<sup>2</sup> *<sup>x</sup>* � *<sup>N</sup>*<sup>2</sup> *y* ,

*A Review on the Contact Mechanics Modeling of Rough Surfaces in the Elastic Regime… DOI: http://dx.doi.org/10.5772/intechopen.102358*

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

$$\overline{z} = \frac{1}{N\_p} \sum\_{i=1}^{N\_p} z\_i \tag{52}$$

$$\sigma^2 = \frac{1}{N\_p} \sum\_{i=1}^{N\_p} \left(\mathbf{z}\_i - \overline{\mathbf{z}}\right)^2 \tag{53}$$

$$S\_k = \frac{1}{N\_p \sigma^3} \sum\_{i=1}^{N\_p} \left( z\_i - \overline{z} \right)^3 \tag{54}$$

*Tribology of Machine Elements - Fundamentals and Applications*

$$K\_u = \frac{1}{N\_p \sigma^4} \sum\_{i=1}^{N\_p} (z\_i - \overline{z})^4 \tag{55}$$

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 xk*, *yl* � � using Eq. (55):

$$\text{PSDF}\left(k\_{\mathbf{x}\_k}, k\_{\mathbf{y}\_l}\right) = \frac{\Delta^2}{\left(2\pi\right)^2 N\_p} \left|\text{FFT}\left(z\left(\mathbf{x}\_k, \mathbf{y}\_l\right)\right)\right|^2,\tag{56}$$

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 *k*∈ *kr*, *ks* ½ � and reaches a plateau *C*<sup>0</sup> for *k*∈ *kl*, *kr* ½ �, 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,

$$\text{PSDF}(k) = \mathbf{C}\_0 \times \begin{cases} 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{cases} \tag{57}$$

where the exponent *a* can be expressed depending on the Hurst component as *a* ¼ �2 1ð Þ þ *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*:*<sup>5</sup> � <sup>1</sup>*:*<sup>5</sup> <sup>μ</sup>m<sup>2</sup> 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.

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 *λ<sup>r</sup>* and then it reaches a plateau. The linear region in the PSDF plot presents two slopes, the first one is proportional to *k*�6*:*<sup>5</sup> (i.e. PSDF∝*k*�6*:*<sup>5</sup> ) for *<sup>k</sup>*<sup>∈</sup> *kr*, *kc* ½ � while the second one is proportional to *<sup>k</sup>*�4*:*<sup>6</sup> (i.e. PSDF∝*k*�4*:*<sup>6</sup> ) 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 *λ<sup>c</sup>* ¼ 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 *C*0, the roll-off

*A Review on the Contact Mechanics Modeling of Rough Surfaces in the Elastic Regime… DOI: http://dx.doi.org/10.5772/intechopen.102358*

#### **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.*

wavelength *λ<sup>r</sup>* and the critical wavelength *λ<sup>c</sup>* are obtained by curve fitting of the different regions, their respective values are 3.8 � <sup>10</sup><sup>2</sup> nm<sup>4</sup> , 90 nm and 50 nm.

#### **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*:*7 nm, (ii) the roll-off frequency *kr* ¼ 2*π=λ<sup>r</sup>* with *λ<sup>r</sup>* ¼ 90 nm, (iii) the critical frequency *λ<sup>c</sup>* ¼ 50 nm (iv) the scanning length *L* ¼ 1*:*5 μm, (v) the sampling length Δ ¼ 2*:*93 nm and (vi) the exponent *a* of the PSDF. The latter takes values of �6*:*5 for *k*∈ *kr*, *kc* ½ � 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 meshconvergence in the literature, several computations were carried out on the 12

*A Review on the Contact Mechanics Modeling of Rough Surfaces in the Elastic Regime… DOI: http://dx.doi.org/10.5772/intechopen.102358*

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

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=A*<sup>0</sup> � 100 is plotted as a function of the normalized contact pressure *P=E*<sup>∗</sup> , where *A* refers to the real contact area, *A*<sup>0</sup> denotes the nominal area, *P* represents the contact pressure and *E*<sup>∗</sup> is the composite Young modulus. The results are briefly summarized in **Figure 9** for both FFT and NNLS solvers.

One can observe from **Figure 9** that the contact law, *<sup>A</sup>=A*<sup>0</sup> � <sup>100</sup> <sup>¼</sup> *f P=E*<sup>∗</sup> ð Þ, 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 32 GB 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

*A Review on the Contact Mechanics Modeling of Rough Surfaces in the Elastic Regime… DOI: http://dx.doi.org/10.5772/intechopen.102358*

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

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 http://contact.engineering/. The later is denoted by BEM-L. Patewska.


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

#### **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.*

*A Review on the Contact Mechanics Modeling of Rough Surfaces in the Elastic Regime… DOI: http://dx.doi.org/10.5772/intechopen.102358*

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

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 FFTbased 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 *A Review on the Contact Mechanics Modeling of Rough Surfaces in the Elastic Regime… DOI: http://dx.doi.org/10.5772/intechopen.102358*

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.
