On the accurate computation of the true contact-area in mechanical contact of random rough surfaces

TL;DR

This paper presents a correction method for accurately computing the true contact area in rough surface contact problems, even with coarse meshes, by using geometrical insights and interface length measurements.

Contribution

A novel correction function is introduced that improves contact area calculations on coarse meshes, validated across various surface types and spectral contents.

Findings

Accurate contact area estimation with coarse meshes achieved.

Correction method validated for different fractal dimensions.

Topology smoothing enhances contact cluster analysis.

Abstract

We introduce a corrective function to compensate errors in contact area computations coming from mesh discretization. The correction is based on geometrical arguments and requires only one additional quantity to be computed: the length of contact/non-contact interfaces. The new technique enables us to evaluate accurately the true contact area using a coarse mesh for which the shortest wavelength in the surface spectrum reaches the grid size. The validity of the approach is demonstrated for surfaces with different fractal dimensions and different spectral content using a properly designed mesh convergence test. In addition, we use a topology preserving smoothing technique to adjust the morphology of contact clusters obtained with a coarse grid.

Equations

$${\Phi }_p(k)=\{\begin{array}{ll}{\Phi }_0,& \text{mbox}ifk\le k_l\text{mbox}(plateau)\\ {\Phi }_0(k/k_l{)}^{-2(H+1)},& \text{mbox}if{k}_l<k\le k_s\\ 0,& \text{mbox}else\end{array},{\Phi }_p(k)=\{\begin{array}{ll}{\Phi }_0(k/k_l{)}^{-2(H+1)},& \text{mbox}if{k}_l<k\le k_s\\ 0,& \text{mbox}else\end{array},$$

$${\Phi }_p(k)=\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}{\Phi }_s(k_x,k_y)d\varphi =\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}{\Phi }_s(k\cos (\varphi ),k\sin (\varphi ))d\varphi$$

$${\Phi }_p(k)=\frac{1}{n(k)}\sum _{\forall k_x,k_y:k_x^2+k_y^2=k^2}{\Phi }_s(k_x,k_y),$$

$$\tilde{k}=kL/2\pi =L/\lambda ,$$

$$p(x,y)=p_0\cos (k_{x}x)\cos (k_{y}y)\iff u^z(x,y)=\frac{2p_0}{E^{\ast }\sqrt{k_x^2+k_y^2}}\cos (k_{x}x)\cos (k_{y}y),$$

$$E^{\ast }=\frac{E_1{E}_2}{(1-{\nu }_1^2)E_2+(1-{\nu }_2^2)E_1},$$

$$[u^z]={\mathcal{F}}^{-1}\left\{[w]:\mathcal{F}([p])\right\},$$

$$[u^z]=u_{ij}^z{\mathbf{e}}_i\otimes {\mathbf{e}}_j,[p]=p_{ij}{\mathbf{e}}_i\otimes {\mathbf{e}}_j$$

$$[w]=w_{kl}{\delta }_{ki}{\delta }_{lj}{\hat{\mathbf{e}}}_k\otimes {\hat{\mathbf{e}}}_l\otimes {\hat{\mathbf{e}}}_i\otimes {\hat{\mathbf{e}}}_j,$$

$$w_{kl}=\frac{L}{\pi E^{\ast }\sqrt{k^2+l^2}}$$

$$\underset{[p]}{\min }\mathbb{F},\mathbb{F}=\frac{1}{2}[p]:[u]+[p]:[g],\text{mbox}underconstraints\frac{1}{N^2}\sum p_{ij}=p_0,p_{ij}\ge 0$$

$$\frac{[p]:[g]}{N^2{p}_0\sqrt{⟨z^2⟩}}\le \epsilon ,$$

$$p\ge 0,g\ge 0,pg=0.$$

$$\frac{A_{\ast }}{A_0}\approx \frac{N_a^i+(1-\beta )N_a^b}{N^2}=\frac{N_a}{N^2}-\beta \frac{N_a^b}{N^2},$$

$$S=\frac{\pi }{4}{S}_d.$$

$$\frac{A_{\ast }}{A_0}\approx \frac{A_d}{A_0}-\beta \frac{\pi }{4}\frac{S_d\Delta x}{A_0},$$

$$\frac{A_d}{A_0}\underset{\Delta x\to 0}{\to }\frac{A_{\ast }}{A_0}.$$

$$a_{\text{mbox}sm}(x_i,\varphi )=\{\begin{array}{ll}\frac{1}{2}(\Delta x-x_i{)}^2\tan (\varphi ),& \text{mbox}if0\le \varphi <\mathrm{atan}(\Delta x/(\Delta x-x_i))\\ \Delta x(\Delta x-x_i)-\frac{\Delta x^2}{2\tan (\varphi )},& \text{mbox}if\mathrm{atan}(\Delta x/(\Delta x-x_i))\le \varphi <\mathrm{atan}(\Delta x/(\Delta x-2x_i))\text{mbox}andx<\Delta x/2\\ \Delta x{x}_i+\frac{\Delta x^2}{2\tan (\varphi )},& \text{mbox}if\mathrm{atan}(\Delta x/(\Delta x-2x_i))\le \varphi <\pi /2\text{mbox}andx<\Delta x/2\\ \Delta x(\Delta x-x_i)-\frac{\Delta x^2}{2\tan (\varphi )},& \text{mbox}if\mathrm{atan}(\Delta x/(\Delta x-x_i))\le \varphi <\pi /2\text{mbox}andx\ge \Delta x/2\\ -x_i^2\tan (\varphi ),& \text{mbox}if\varphi >\pi -\mathrm{atan}(\Delta x/x_i)\\ \Delta x{x}_i+\frac{\Delta x^2}{2\tan (\varphi )},& \text{mbox}if\pi -\mathrm{atan}(\Delta x/(2x_i-\Delta x))\le \varphi <\pi -\mathrm{atan}(\Delta x/x_i)\text{mbox}andx\ge \Delta x/2\\ \Delta x(\Delta x-x_i)-\frac{\Delta x^2}{2\tan (\varphi )},& \text{mbox}if\pi /2\le \varphi <\pi -\mathrm{atan}(\Delta x/(2x_i-\Delta x))\text{mbox}andx\ge \Delta x/2\\ \Delta x{x}_i+\frac{\Delta x^2}{2\tan (\varphi )},& \text{mbox}if\pi /2\le \varphi <\pi -\mathrm{atan}(\Delta x/x_i)\text{mbox}andx<\Delta x/2\end{array}$$

$$⟨a_{\text{mbox}sm}⟩=\underset{0}{\overset{\Delta x}{\int }}\underset{0}{\overset{\pi }{\int }}{a}_{\text{mbox}sm}(x_i,\varphi )P_x(x_i)P_{\varphi }(\varphi )d\varphi d{x}_i=\frac{1}{\pi \Delta x}\underset{0}{\overset{\Delta x}{\int }}\underset{0}{\overset{\pi }{\int }}{a}_{\text{mbox}sm}(x_i,\varphi )d\varphi d{x}_i$$

$$\frac{⟨a_{\text{mbox}sm}⟩}{\Delta x^2}=\frac{\pi -1+\ln 2}{6\pi }\approx 0.150387619$$

$$\boxed{A_{\ast }\approx A_d-\frac{\pi -1+\ln 2}{24}{S}_d\Delta x}$$

$${\frac{A}{A}}_0=\kappa \frac{p_0}{E^{\ast }\sqrt{⟨|\nabla z{|}^2⟩}},$$

$$\frac{A}{A_0}=\mathrm{erf}\left(\frac{\sqrt{2}{p}_0}{E^{\ast }\sqrt{⟨|\nabla z{|}^2⟩}}\right).$$

$$\sqrt{⟨|\nabla z{|}^2⟩}=\sqrt{2m_2},$$

$${\kappa }_i=\kappa (p_0^i)=E^{\ast }\sqrt{⟨|\nabla z{|}^2⟩}{\frac{\partial A/A_0}{\partial p_0}|}_{p_0^i}\approx \frac{E^{\ast }\sqrt{⟨|\nabla z{|}^2⟩}}{A_0}\frac{A^{i+1/2}-A^{i-1/2}}{p_0^{i+1/2}-p_0^{i-1/2}}.$$

$$\underset{p_0\to 0}{\lim }\frac{E^{\ast }\sqrt{⟨|\nabla z{|}^2⟩}}{A_0}\frac{\partial A}{\partial p_0}={\kappa }_0.$$

$${\kappa }^c=\frac{E^{\ast }\sqrt{⟨|\nabla z{|}^2⟩}}{A_0}\frac{\partial A_c}{\partial p_0}=\frac{E^{\ast }\sqrt{⟨|\nabla z{|}^2⟩}}{A_0}\left(\frac{\partial A}{\partial p_0}-\frac{(\pi -1+\ln 2)\Delta x}{24}\frac{\partial S}{\partial p_0}\right)$$

$${\bar{\kappa }}_i=\frac{1}{2}{\kappa }_i+\frac{1}{4}\left({\kappa }_{i-1}+{\kappa }_{i+1}\right).$$

$${\kappa }_P(p_0)=\frac{\partial \mathrm{erf}(\sqrt{2}{p}^{\prime })}{\partial p^{\prime }}=\sqrt{\frac{8}{\pi }}\exp \left[-\frac{p_0^2}{E^{\ast 2}⟨|\nabla z{|}^2⟩}\right],$$

$${\frac{A}{A}}_0\approx {\kappa }_0{p}^{\prime }-c{p}^{\prime 2}.$$

Full text

On the accurate computation of the true contact-area in mechanical contact of random rough surfaces

Vladislav A. Yastrebova111Corresponding author $\langle$vladislav.yastrebov@mines-paristech.fr$\rangle$, Guillaume Anciauxb, Jean-François Molinarib

(*a**MINES ParisTech, PSL Research University, Centre des Matériaux

CNRS UMR 7633, BP 87, 91003, Evry, France

bCivil Engineering Department, Materials Science Department,

Ecole Polytechnique Fédérale de Lausanne (EPFL), Station 18, 1015, Lausanne, Switzerland*)

Abstract.

We introduce a corrective function to compensate errors in contact area computations coming from mesh discretization. The correction is based on geometrical arguments and requires only one additional quantity to be computed: the length of contact/non-contact interfaces. The new technique enables us to evaluate accurately the true contact area using a coarse mesh for which the shortest wavelength in the surface spectrum reaches the grid size. The validity of the approach is demonstrated for surfaces with different fractal dimensions and different spectral content using a properly designed mesh convergence test. In addition, we use a topology preserving smoothing technique to adjust the morphology of contact clusters obtained with a coarse grid.

Keywords: roughness, contact, simulations, boundary element method

1 Introduction

The roughness of natural and industrial surfaces determines properties of the mechanical contact between solids: interface stiffness, true contact area and the morphology of the free interface volume. Thus the roughness governs many interface phenomena such as contact electrical resistance [67, 76], thermal contact resistance [41, 4], friction [18, 6, 61], adhesion [23, 45], wear [8, 2], as well as fluid transport at contact interfaces [42, 51, 62, 15, 16, 47]. For most of these phenomena it is critical to accurately estimate the true contact area for given thermo-electro-mechanical loads and given roughness of contacting surfaces. It is now well known that a simple load bearing area model, relying on a geometrical overlap of two rough surfaces considerably overestimates the true contact area [46, 60], and for equal contact areas, the former results in a much higher transmissivity in transport problems [15]. Existing analytical models, asperity-based [27, 9, 66, 38, 26, 1], as well as Persson’s model [50, 52] with its adjusted version [49], rely on a few approximations and thus cannot provide very accurate results in terms of true contact area over a wide interval of loading conditions (for a detailed discussion and comparison see [37, 12, 44, 73, 75]).

For these reasons, a numerical analysis, free of restrictive assumptions, is now widely used for the study of rough contact: the finite element method [46, 77], a wide class of continuum boundary element methods [65, 54, 55, 36], discrete methods based on molecular dynamics [10, 11, 71] or basic molecular dynamics [3, 64, 45]. Continuum models are particularly attractive since they permit to cover a large spectrum of length scales. However, they are subject to discretization and convergence errors. The former is related to the finite size of the used grid/mesh, whereas the latter is related to the strong non-linearity and discontinuous nature [35, 21] of contact problems requiring iterative procedures (Newton-Raphson method, iterative solvers) to achieve convergence, or in the case of explicit techniques for both finite element [29] and the Green’s function molecular dynamics [10, 11] obtaining the results requires damping of elastic vibrations. Various continuous numerical methods in contact mechanics exist [70], among them penalty and barrier methods, Lagrange multipliers, augmented Lagrangian and other techniques, which convert the constrained optimization problem to an unconstrained one (or at least partially unconstrained). Some of these methods allow accurate satisfaction of contact constraints for a given discretization whereas others (like penalty or barrier type method) only approach the exact solution with an accuracy that depends on the choice of parameters. In addition, different contact discretization techniques, which integrate contact tractions in the weak form in the finite-element method, provide varying accuracy and convergence rates. These depend on the interpolation order of elements, mesh densities on contacting surfaces and mesh curvature. For details, see [70, 58, 22, 69, 56, 72, 19] and references therein.

Assuming that the convergence is ensured and that the numerical method is accurate enough, the discretization remains the sole limiting point in achieving accurate results in contact problems. From experience, it is known that for non-conformal but simple geometries (1D or 2D wave, circle or sphere against a flat surface), a rather dense mesh is required at the contact interface to track the contact area evolution, see e.g. [39, 40, 31, 20, 77, 63]. We also refer to a study of a bi-wavy surface in contact with a rigid flat [74], in which a very dense grid (4096 points per side per wavelength) was used, which allowed the authors to reveal peculiar mean contact pressure behaviour near the percolation point, which was missed in previous studies. Extrapolating to the case of “rough” surfaces, we need to ensure accurate discretization of all harmonics present in the surface spectrum, i.e. the shortest wavelength $\lambda_{s}$ should be sufficiently discretized to ensure accurate estimation of the contact area. This statement implies that the surface “roughness” should be smooth enough (or should be smoothly interpolated from experimental data  [30, 77]) and the ratio $\lambda_{s}/\Delta x$, where $\Delta x$ is the grid step, should be kept rather big [73, 57]. Importantly, as was demonstrated in [75], the discretization error is affected not only by the shortest wavelength $\lambda_{s}$ but also by the longest wavelength $\lambda_{l}$ in the spectrum222For surfaces with a plateau in the surface spectrum, $\lambda_{l}$ corresponds to the shortest wavelength present in the plateau as in examples shown in Fig. 1(c-e)., since it determines the number of individual macro-contact clusters (see [48, 75]). Thus, since the error in the true contact area is proportional to the length of the boundary between contact and non-contact zones, a shorter longer-wavelength, results in more contact clusters and leads to a higher discretization error. Hence, for any rough surface with a sufficiently large discrete spectrum, the discretization requirement may rapidly become a bottle neck in terms of computational resources. Two alternative solutions can be used: sacrifice the accuracy by using a coarse mesh and/or use new more efficient numerical methods [54, 5] and adapted hardware combining CPU and GPU and/or parallel algorithms. A notable example, of combining a Green’s function molecular dynamics with computations on GPU enables researchers to make rough contact simulations on a grid with more than 17 billion(!) grid points [57].

The true contact area in numerical simulations of contact can be computed as the total area of surface-faces being in contact plus the areas associated with nodes on the contact-non contact border. In spectral methods [65, 36] using FFT techniques and in discrete techniques [10, 11] requiring regular discretizations (equally spaced grid points), the true contact area fraction can be simply computed as a ratio between all points in contact (points with zero gap and non-zero pressure) to the total number of points. However, this estimation of the contact area slightly overestimates the true solution. The coarser the discretization, the higher the error: the convergence rate with element size is linear [75].

Here, we suggest an alternative approach that allows an estimation of the true contact area with high accuracy on a reasonably coarse mesh, i.e. enabling full representation of the surface spectrum. The approach is based on a corrective function which uses the length of the contact/non-contact boundary (or simply the perimeter of contact clusters). The method was already introduced in [75], but the previous study lacked a mesh convergence analysis and the corrective factor was not evaluated. Here, we correct these shortcomings and demonstrate the accuracy of the suggested techniques on numerous cases.

The paper is organized as follows: in Section 2, all numerical methods and models are presented: a method used to generate model roughness in Subsection 2.1, a spectral method used to solve contact mechanical problems in Subsection 2.2. The area-correction method is described in detail in Subsections 2.3 and 2.4. In Section 3, we present mesh convergence tests and demonstrate the performance of the approach. A morphology smoothing of contact clusters is briefly outlined in Subsection 3.3. A short conclusion is drawn in Section 4.

2 Methods

2.1 Rough surface generation

We use an FFT filtering technique [28] to synthesize an isotropic 2D random rough surface with prescribed Hurst exponent $H$ (or, equivalently, fractal dimension $D=3-H$) and cut-offs in the surface spectrum as was done in [73, 75]. Note that we use a periodic surface resulting in a discrete spectrum. Surfaces with and without plateau in the power spectral density (PSD) are considered (see [32] for a detailed discussion of PSD measurements and interpretation). One can compare surfaces with plateau depicted in Fig. 1(c),(e) with those without it, see Fig. 1(d),(f). The surface PSD is characterized by the following parameters: a scaling factor $\Phi_{0}$, Hurst exponent $H$, absence or presence of a plateau and also two wavenumbers $k_{l}$ and $k_{s}$: $k_{l}$ the lowest wavenumber determining the start of a power-law decaying PSD and $k_{s}$ the highest wavenumber determining the shortest wavelength present in the spectrum. The mean profile PSD333Lower index “p” is used for profile PSD and lower index “s” is used for surface spectrum. for a surface with and without plateau can be formalized in the following forms, respectively:

[TABLE]

where $k=\sqrt{k_{x}^{2}+k_{y}^{2}}$, with $k_{x},k_{y}$ are wavenumbers in $x$ and $y$ directions ($k_{x,y}=2\pi/\lambda_{x,y}$), respectively, and

[TABLE]

is the mean profile PSD for a given absolute value of the wavevector $k$ in all directions determined by $\phi$. The surface is considered to be isotropic, so that statistical characteristics of every profile should be independent of the direction $\phi$. Since we deal with a discrete spectrum the integral in the previous expression should be rewritten as:

[TABLE]

where $n(k)$ is the number of elements satisfying $k_{x}^{2}+k_{y}^{2}=k^{2}$ for all $k_{x},k_{y}\in\mathbb{Z}$. Because of this discreteness, the accuracy of definition (1) cannot be satisfied for a single surface. Thus, multiple realizations of rough surfaces for a given set of parameters $k_{l},k_{s},H$ is needed to capture the statistical nature of the roughness. It is convenient to introduce dimensionless wavenumbers $\tilde{k}$, which will be used throughout the paper:

[TABLE]

where $\lambda$ is the wavelength and $L$ is the length of a side of simulation box. It means that $\tilde{k}$ is the number of waves per length, and because of the periodicity this quantity takes only integer values. Some examples of generated rough surfaces with and without plateau are depicted in Fig. 1.

2.2 Spectral boundary element method

The method that we use here is based on a variational formulation of contact problems given by Kalker [34] and on a simple relationship between a wavy pressure $p(x,y)$ profile and the resulting wavy vertical displacement profile $u^{z}(x,y)$ due to [68, 33, 65], which can be summarized as:

[TABLE]

where $x,y$ are in-plane coordinates, $E^{*}$ is the effective elastic modulus of a pair of contacting elastic materials, given by:

[TABLE]

where $E_{i},\nu_{i}$ are Young’s moduli and Poisson’s ratios of materials $i=1,2$, respectively. Relationship (2) is easy to generalize for an arbitrary superposition of modes using a discrete Fourier transformation $\mathcal{F}$ [65]:

[TABLE]

where by colon we denote a double scalar product, and

[TABLE]

are $N\times N$ matrices of vertical displacements and pressures, respectively, defined on a discrete grid $\{x_{i},y_{j}\}$, for $i,j\in[0,N-1]$, with $N$ be the number of grid points per period $L$ and $\boldsymbol{e}_{i}$ is an orthonormal vector basis of dimension $N$, i.e. $\boldsymbol{e}_{i}\cdot\boldsymbol{e}_{j}=\delta_{ij}$, where $\delta_{ij}$ is the Kroneker symbol. By $[w]$ we denote an $N\times N\times N\times N$ matrix of the form:

[TABLE]

where $\hat{\boldsymbol{e}}_{i}$ is another orthonormal basis in Fourier space such that $\mathcal{F}(p_{ij}\boldsymbol{e}_{i}\otimes\boldsymbol{e}_{j})=\hat{p}_{ij}\hat{\boldsymbol{e}}_{i}\otimes\hat{\boldsymbol{e}}_{j}$ and

[TABLE]

A frictionless non-adhesive normal contact between two linearly elastic solids and roughness $z_{1}(x,y)$ and $z_{2}(x,y)$ can be mapped on a contact problem between a rigid rough surface with an effective roughness $z=z_{2}-z_{1}$ and an elastic half space with an effective elastic modulus $E^{*}$ introduced above. To use linear elastic approximation and the boundary element method, we assume that the roughness slope remains everywhere vanishingly small $|\nabla z|\ll 1$. The problem can be formulated as a constrained optimization problem:

[TABLE]

where $p_{0}$ is the applied pressure and $[g]=g_{ij}\boldsymbol{e}_{i}\otimes\boldsymbol{e}_{j}$ is the matrix of gaps $g_{ij}=z_{ij}-u_{ij}+z_{0}$, where $z_{0}$ is an offset value. This optimization problem can be solved by different optimization techniques [7, 5]. Here we use the methodology suggested in [54]. The solution of the optimization problem is a matrix of pressures $[p*]$, which takes positive values in contact zones, where the gap vanishes $g=0$, and zero values $p=0$ where gap is positive. At solution, the orthogonality of the gap and pressure is ensured with a given precision $\varepsilon$:

[TABLE]

where $\sqrt{\langle z^{2}\rangle}$ is the root mean squared roughness. Thus, in average we approximately satisfy Hertz-Signorini-Moreau contact conditions [70] at every grid point:

[TABLE]

2.3 Computation of the contact area

To compute the discrete contact area $A_{d}$ we simply count the number of grid points $N_{a}$ with positive pressure $p_{ij}>0$, each point contributing an area $\Delta x^{2}$, where $\Delta x=L/N$, $L$ is the side length and $N$ is the number of grid points per side. The ratio $A_{d}/A_{0}=N_{a}/N^{2}$ is then the approximate contact area fraction, since $A_{0}=N^{2}\Delta x^{2}$ is the nominal contact area. The grid points in contact are split into two types: internal points $N_{a}^{i}$, which do not have non-contact grid points in their nearest neighbourhood and boundary contact points $N_{a}^{b}$, which have non-contact points in their nearest neighbourhood; $N_{a}=N_{a}^{i}+N_{a}^{b}$. We suppose that the true (obtained for $N\to\infty$) boundary between contact and non-contact zones is located somewhere in between of these discrete measurements. Thus the true contact area $A_{*}$ can be estimated as

[TABLE]

where $0<\beta<1$. The number of boundary contact points $N_{a}^{b}$ correlates directly with the perimeter of contact/non-contact boundaries $S_{d}=M\Delta x$, where $M$ is the number of switches from contact to non-contact and vice versa along grid points both in vertical and in horizontal directions (Fig. 2). Since for a finite $k_{s}$ and in the limit of zero grid spacing $N\to\infty$ the contact clusters are smooth figures, the discrete perimeter $S_{d}$ should be converted into a smooth perimeter $S$, which for an arbitrary isotropic shape can be obtained as

[TABLE]

The factor $\pi/4$ appears in the conversion from a discrete perimeter measurement to a continuous one because a straight line inclined by an arbitrary angle to a regular Cartesian grid has in average a true length which is $4/\pi$ times shorter than the length computed on a Cartesian grid with vertical and horizontal lines444In [45] for the same correction the authors used a slightly different factor $4\mathrm{arcsinh}(1)/\pi\approx 1.1222$, whereas our factor $4/\pi\approx 1.2732$.. It is similar to the conversion of $L_{1}$ (or Manhattan) metrics to $L_{2}$ (or Euclidean) metrics. To give an example, a circle with radius $R$ has a discrete perimeter $S_{d}\approx 8R$, which is $4/\pi$ times bigger than the true perimeter $S=2\pi R$. Finally, the true contact area fraction can be estimated as

[TABLE]

where the factor $\beta$ will be determined in the following section. An alternative formulation of this correction can be found in [75, Sect.9]. Since $S_{d}$ is almost independent of $\Delta x$, and thus $S_{d}\Delta x\to 0$ as $\Delta x\to 0$, from Eq. (5) it follows that

[TABLE]

Note that the value of the correction $\beta\pi S_{d}\Delta x/4A_{0}$, which is nothing but the absolute error, is proportional to discretization step $\Delta x$.

2.4 Corrective factor

Let us take a square with side $\Delta x$ and cut it into two parts (left and right) by a randomly oriented straight line. Let us make the following assumptions without loss of generality (see Fig. 3(a)):

Let us consider a square occupying $x\in[0,\Delta x]$, $y\in[0,\Delta x]$. 2. 2.

Let one of the intersections $\{x_{i},y_{i}\}$ of the cut line with the sides of the square lie on the side with $y=0$, so that $x_{i}\in[0,\Delta x]$, $y_{i}=0$ with a uniform probability distribution $P_{x}(x_{i})=1/\Delta x$. 3. 3.

The angle of inclination555A seemingly simpler approach would be to introduce a cut position $x$ on one side of the square and a cut position $y$ on another side, but in this case the probability density of the second cut position depends on the first cut and thus cannot be considered uniform which would complicate the following analysis. of the line with respect to axis $0X$ is denoted $\phi\in[0,\pi]$ with a uniform probability distribution $P_{\phi}(\phi)=1/\pi$.

The area of the smaller part $a_{\mbox{\tiny sm}}$ can be easily computed for any given $x_{i}$ and $\phi$, the remaining (bigger) part has area $\Delta x^{2}-a_{\mbox{\tiny sm}}$:

[TABLE]

The mean area of the smaller part can be computed as follows

[TABLE]

The integrand of (6) is plotted in Fig. 3(b) and can be evaluated analytically in closed form (see Appendix A) as:

[TABLE]

We assume that in our numerical code we overestimate the contact area at border grids by this smaller portion of a square, so replacing factor $\beta$ in (5) by expression in (7) we obtain the ultimate equation for the true contact area estimation:

[TABLE]

3 Mesh convergence tests

3.1 Contact area

To validate our area correction technique, we present a mesh convergence test. We carried out several simulations on surfaces for two Hurst exponents $H=0.4,\;0.8$ and the following combinations of cut-off wavenumbers666Here, by wavenumber we mean a normalized dimensionless wavenumber $\tilde{k}=L/\lambda$. $\{\tilde{k}_{l},\tilde{k}_{s}\}=\{1,32\},\{1,64\}$ with discretizations $N\times N$ with $N=128,256,512,1024,2048,4096$, and also for $\{\tilde{k}_{l},\tilde{k}_{s}\}=\{4,128\},\{4,256\}$ with discretizations $N=256,512,1024,2048,4096$. Coarse roughness discretizations were obtained from a fine discretization $N=8192$ by sampling every $2,4,8,16$ or $32$ points. In all considered cases, the formula (8) for corrected true contact area gives accurate results independent of discretization. This suggests that simulations of non-adhesive contact can be done on a rough surface with a very coarse discretization $N=2\tilde{k}_{s}=2L/\lambda_{s}$ but with almost the same accuracy(!) as for $N=32L/\lambda_{s}$. This represents a gain factor of $16^{2}=256$ in terms of discretization points. It is worth mentioning that the accuracy can be ensured solely for the estimation of the global true contact area but not for other important quantities such as local pressure or gap distribution, which would still require fine discretizations or different corrective techniques. This correction will allow to perform without loss of accuracy rough contact simulations with very broad surface spectra (high Nayak’s parameter [43]) encountered in real engineering surfaces. To demonstrate the performance of our technique, we show the contact area evolution in Fig. 4 ((a) raw data, (b) corrected data using Eq. (8)) for $H=0.8$, $\tilde{k}_{l}=1$, $\tilde{k}_{s}=32$ and different discretizations; results for $H=0.8$, $\tilde{k}_{l}=4$, $\tilde{k}_{s}=128$ are depicted in Fig. 5.

According to analytical models based on asperity considerations [27, 9, 38, 66, 26, 12] or based on the evolution of the probability density of contact pressures under increasing magnification, an idea developed by Persson [50, 49, 37, 17], the contact area fraction $A/A_{0}$ is proportional to the applied pressure $p_{0}$ at infinitesimal pressure $p_{0}$. More precisely, the relation stands as:

[TABLE]

where $\kappa$ is the coefficient of proportionality and $\sqrt{\langle|\nabla z|^{2}\rangle}$ is the root mean squared roughness gradient. The coefficient $\kappa$ was found to be $\sqrt{2\pi}$ in all asperity based models (which use Nayak’s random process model for roughness description [9, 38, 66, 26, 12]), and was deduced to be $\sqrt{8/\pi}$ in Persson’s model. In Figs. 4 and 5, to provide bounds, we plot also the asymptotic solution (9) of asperity based models and Persson’s solution for the contact area evolution given by:

[TABLE]

Note that in Fig. 5 the results are obtained for a surface with a plateau in the PSD similar to the one shown in Fig. 1(c).

It is important to remark that to carry out this convergence study, the measure of the rms surface gradient should be done using the following equivalence:

[TABLE]

where $m_{2}$ is the second spectral moment [43, 26]. Since the spectrum and thus the power spectral density of every surface with a cut-off wavenumber $\tilde{k}_{s}<N/2$ is independent of the discretization $N$, by measuring the rms surface gradient as the average second spectral moment $2m_{2}=m_{02}+m_{20}$ we can find the exact value of $\sqrt{\langle|\nabla z|^{2}\rangle}$ in the continuous limit. Therefore we can ensure the proper normalization of the pressure. If the rms surface gradient is measured on the real space geometry using forward finite difference (FFD) or central finite difference (CFD) schemes, it becomes a discretization dependent parameter [25, 44] and thus the convergence study on different grids becomes ambiguous as both the normalized pressure and the contact area change simultaneously.

3.2 Contact area slope

To demonstrate the accuracy of the contact area evaluation using the suggested technique, we plot the derivative of the contact area with respect to the applied normalized pressure $p_{0}/E^{*}\sqrt{\langle|\nabla z|\rangle}$. The slope denoted $\kappa(p_{0})$ is evaluated as a finite difference

[TABLE]

Defined in this way, the slope should tend to the area-pressure proportionality coefficient as pressure tends to zero:

[TABLE]

This proportionality coefficient $\kappa_{0}$, predicted by asperity based model [9, 66, 38, 26, 12] to be $\sqrt{2\pi}$ and by Persson’s model [50] to be $\sqrt{8/\pi}$, was a topic of numerous numerical studies [29, 30, 44, 11, 59, 73, 53, 57, 75]. To evaluate the slope correctly we use the corrected contact area and thus compute it as follows (the upper index $c$ indicates the corrected value):

[TABLE]

To slightly reduce oscillations we smooth the slope using:

[TABLE]

In Figs. 6 and 7 the slopes are computed for $\tilde{k}_{l}=4$, $\tilde{k}_{s}=128$ for $H=0.8$ and $H=0.4$, respectively and for different grids with $N=256,512,1024,2048,4096$ using raw data with Eq. (11) and corrected data with Eq. (12), both results were slightly smoothed using Eq. (13). The results are compared with the slope of Persson’s model [50] given by

[TABLE]

where $p^{\prime}=p_{0}/E^{*}\sqrt{\langle|\nabla z|^{2}\rangle}$. Note that the effect of the area correction is considerable even for the finest discretization $N=4096$ especially in the interval of small pressures.

The corrected data of the slopes can be accurately fitted by a linear function suggesting that the contact area evolves approximately as a polynomial of the second degree:

[TABLE]

For the presented results, $\kappa_{0}\approx 2.05$, $c\approx 2.35$ for $\tilde{k}_{l}=4$, $\tilde{k}_{s}=128$, $H=0.8$ and $\kappa_{0}\approx 2.15$, $c\approx 2.80$ for $\tilde{k}_{l}=4$, $\tilde{k}_{s}=128$, $H=0.4$. These values of $\kappa_{0}$ are comparable with those found in the literature [29, 30, 44, 11, 59, 73, 53, 57, 75]. However, a more rigorous statistical analysis would be needed to investigate the separate roles of the fractal dimension (or Hurst exponent) and Nayak’s parameter.

3.3 Morphology of contact clusters

The integral measure of the contact area is an important but not the only characteristic that determines physical properties of contact interfaces. The morphology of contact clusters and their distribution also play an important role in many interfacial phenomena [24, 16]. In Fig. 8 we show maps of contact clusters computed at different discretizations for applied pressure $p_{0}/E^{*}\sqrt{\langle|\nabla z|\rangle}\approx 0.12$ giving $A/A_{0}\approx 22.5$ % for $H=0.8$, $\tilde{k}_{l}=4$, $\tilde{k}_{s}=128$. Three discretization-dependent morphological effects can be distinguished: 1) loss of small contact spots at coarser grids (marked with A in Fig. 8) (see discussion in [17, 5]), 2) loss of small non-contact spots at coarser grids (marked with B in Fig. 8), 3) ambiguity in connecting or non-connecting between contact clusters connected by a single grid point (marked with C in Fig. 8) or a single next-nearest connection (marked with D in Fig. 8). The error 3) is probably the most critical for mass transport analysis in near-percolation regime. However, at reasonable discretizations $N=4L/\lambda_{s}$ ($\tilde{k}_{s}=N/4$) in most cases this problem does not persist. Overall, even the coarsest discretization $N=2L/\lambda_{s}$ ($\tilde{k}_{s}=N/2$) shows a reasonable estimation of contact clusters with most morphological features present like on much finer discretizations.

To adjust better the morphology of contact clusters obtained with a coarse grid, we use a smoothing algorithm that preserves topology [14] (see Figs. 9, 10) of the contact cluster map, which is an important feature for rough contact problems. The general method is based on the simultaneous smoothing by metric disks of increasing radius and homotopic cutting/filling, which allow to preserve the topology. For the case of zoomed pixelized images (macropixels of minimal size contain several pixels, i.e. $\times 2$ zoom results in $2\times 2$ pixels) the method preserves centres of macropixels and smooths the image; the chessboard metric $d(x,y)=\max\{|x_{1}-y_{1}|,|x_{2}-y_{2}|\}$ is used to construct the topological skeleton of the original image and the Euclidean metric is used for the homotopic filter. All details of the method can be found in [14]. Moreover, a source code to deal with raster images is provided by the authors of the algorithm on their web-page [13].

Figs. 9 and 10 present the results of the smoothing procedure, in which every contact point (black pixel) in the raw data was replaced by $4\times 4$ pixels before running the smoothing procedure. The figures correspond to $H=0.8$, $\tilde{k}_{l}=1$, $\tilde{k}_{s}=32$ and $H=0.8$, $\tilde{k}_{l}=4$, $\tilde{k}_{s}=128$, respectively. Even though such a topology-preserving smoothing cannot resolve the aforementioned problems 1)-3) in cluster computations, it nevertheless delivers improved contact clusters very similar to those obtained with a much finer discretization. Contact clusters connected only through a next-nearest neighbour connection, can be considered by the algorithm to be topologically disjoint (connection=0) or joint (connection=1) (compare (b) and (c) in Fig. 10), which will result in different topological connections critical for the study of percolation. According to our analysis (see example in Fig. 11), it is roughly equivalently probable to find in finely discretized simulations next-nearest neighbours of the coarse mesh to be connected or disconnected. Thus, in percolation analysis assuming connection=0 and =1 could provide us with upper $A_{\mbox{\tiny sm}}^{u}$ and lower $A_{\mbox{\tiny sm}}^{l}$ bounds of the percolation limit. However, an accurate comparison between these bounds and the true contact area revealed that this topological smoothing gives a contact area estimation close to the reference coarse simulation without correction, i.e. strongly overestimates the true contact area given by (8). It follows that morphological smoothing, even if it provides more realistic micro-contact morphologies, cannot be used for quantitative estimation of the true contact area.

4 Conclusion

In this paper we suggested a correction for the contact area computation in numerical analysis of non-adhesive frictionless contact between elastic rough surfaces. Simple geometrical considerations and the evaluation of the perimeter of contact clusters inspired from our previous work [75] enables us to obtain an accurate estimation of the contact area. The technique is validated using a properly designed mesh convergence test on a particular spectral boundary element method [65], which uses a regular grid combined with conjugate gradient method from [54]. We believe that this approach will be helpful for all spectral based techniques. We demonstrate on several examples that this correction, which vanishes as the grid spacing tends to zero, ensures an accurate contact area estimation even for a coarse grid with grid spacing $\Delta x=\lambda_{s}/2$, where $\lambda_{s}$ is the shortest wavelength present in the roughness spectrum. Moreover, the correction enables us to compute accurately the contact area slope. These results suggest that when accurately computed, the derivative of the contact area with respect to the contact pressure can be accurately approximated by a linear decreasing function at least up to $50$ % of contact area, thus the contact area evolution can be properly described by a second order polynomial without a zero-degree coefficient.

In addition to the integral measure of the contact area, we used a topology preserving smoothing technique developed in [14], which improves the appearance of contact clusters morphology. However, this smoothing does not permit to obtain an accurate quantitative measure of the true contact area.

Overall, using very moderate grid sizes the area correction technique allows to access accurate integral measure of the true contact area in the mechanics of rough contact. Thus it opens new perspectives in numerical studies of contact properties of rough surfaces.

Acknowledgement

Enriching discussions with Valentine Rey and Lucas Frérot are greatly appreciated.

Appendix A Integration of the corrective factor

To evaluate the size of the smallest part of a square cut by a randomly oriented line we split the integral (16) into four parts:

[TABLE]

where integrals can be evaluated separately yielding the following expressions:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Which finally gives the following expression:

[TABLE]

We used the following relationships:

[TABLE]