The role of the roughness spectral breadth in elastic contact of rough surfaces

TL;DR

This paper investigates how the spectral breadth of surface roughness influences the true contact area in elastic contact, revealing a logarithmic decrease with the Nayak parameter and proposing a phenomenological model linking contact area and pressure.

Contribution

It introduces a highly accurate computational method to analyze the effect of spectral breadth on contact area evolution and proposes a new phenomenological equation incorporating the Nayak parameter.

Findings

Contact area decreases logarithmically with the Nayak parameter at constant pressure.

The effect of spectral breadth is significant for surfaces with rich spectra.

A phenomenological model for contact area growth is proposed, depending on the Nayak parameter.

Abstract

We study frictionless and non-adhesive contact between elastic half-spaces with self-affine surfaces. Using a recently suggested corrective technique, we ensure an unprecedented accuracy in computation of the true contact area evolution under increasing pressure. This accuracy enables us to draw conclusions on the role of the surface's spectrum breadth (Nayak parameter) in the contact area evolution. We show that for a given normalized pressure, the contact area decreases logarithmically with the Nayak parameter. By linking the Nayak parameter with the Hurst exponent (or fractal dimension), we show the effect of the latter on the true contact area. This effect, undetectable for surfaces with poor spectral content, is quite strong for surfaces with rich spectra. Numerical results are compared with analytical models and other available numerical results. A phenomenological equation for the contact area growth is suggested with coefficients depending on the Nayak parameter. Using this equation, the pressure-dependent friction coefficient is deduced based on the adhesive theory of friction. Some observations on Persson's model of rough contact, whose prediction does not depend on Nayak parameter, are reported. Overall, the paper provides a unifying picture of rough elastic contact and clarifies discrepancies between preceding results.

Tables

Table 1. Table 1: Coefficients a , b , a ′ , b ′ 𝑎 𝑏 superscript 𝑎 ′ superscript 𝑏 ′ a,b,a^{\prime},b^{\prime} for the linear interpolation for the contact area slope d ​ A ′ / d ​ p ′ ≈ a − 2 ​ b ​ p ′ 𝑑 superscript 𝐴 ′ 𝑑 superscript 𝑝 ′ 𝑎 2 𝑏 superscript 𝑝 ′ dA^{\prime}/dp^{\prime}\approx a-2bp^{\prime} and for the contact area secant A ′ / p ′ ≈ a ′ − b ′ ​ p ′ superscript 𝐴 ′ superscript 𝑝 ′ superscript 𝑎 ′ superscript 𝑏 ′ superscript 𝑝 ′ A^{\prime}/p^{\prime}\approx a^{\prime}-b^{\prime}p^{\prime} found by the least square fit of numerical results in the interval p ′ ∈ [ 0.006 , 0.06 ] superscript 𝑝 ′ 0.006 0.06 p^{\prime}\in[0.006,0.06] identified for H = 0.4 𝐻 0.4 H=0.4 .

$L/{\lambda }_l$	$L/{\lambda }_s$	$\alpha \pm \mathrm{std}\left(\alpha \right)$	$b$	$a$	$b^{\prime }$	$a^{\prime }$
2	16	4.64 $\pm$ 0.24	3.725	2.355	4.909	2.426
2	32	8.05 $\pm$ 0.55	1.761	2.181	3.581	2.293
2	64	14.08 $\pm$ 1.12	1.589	2.154	2.47	2.208
2	128	24.75 $\pm$ 2.08	1.541	2.134	2.194	2.174
2	256	43.52 $\pm$ 3.81	1.4	2.115	1.744	2.136
2	512	76.24 $\pm$ 6.83	1.332	2.099	1.518	2.111
4	16	2.66 $\pm$ 0.09	2.652	2.32	3.732	2.38
4	32	4.44 $\pm$ 0.14	1.889	2.243	2.786	2.301
4	64	7.69 $\pm$ 0.22	1.802	2.199	2.317	2.231
4	128	13.52 $\pm$ 0.43	1.644	2.166	2.083	2.193
4	256	23.84 $\pm$ 0.79	1.523	2.133	1.897	2.156
4	512	41.86 $\pm$ 1.41	1.465	2.112	1.675	2.125
8	16	1.75 $\pm$ 0.01	3.152	2.444	3.294	2.446
8	32	2.59 $\pm$ 0.04	2.308	2.329	2.749	2.357
8	64	4.29 $\pm$ 0.08	1.89	2.242	2.648	2.29
8	128	7.46 $\pm$ 0.13	1.909	2.205	2.559	2.246
8	256	13.14 $\pm$ 0.22	1.797	2.168	2.148	2.19
8	512	23.13 $\pm$ 0.4	1.553	2.124	1.98	2.15
16	32	1.75 $\pm$ 0.01	2.905	2.436	2.962	2.442
16	64	2.58 $\pm$ 0.02	2.494	2.339	2.808	2.359
16	128	4.3 $\pm$ 0.05	2.314	2.267	2.997	2.308
16	256	7.49 $\pm$ 0.1	1.949	2.199	2.525	2.236
16	512	13.18 $\pm$ 0.2	1.63	2.142	2.27	2.181

Table 2. Table 2: Coefficients a , b , a ′ , b ′ 𝑎 𝑏 superscript 𝑎 ′ superscript 𝑏 ′ a,b,a^{\prime},b^{\prime} for the linear interpolation for the contact area slope d ​ A ′ / d ​ p ′ ≈ a − 2 ​ b ​ p ′ 𝑑 superscript 𝐴 ′ 𝑑 superscript 𝑝 ′ 𝑎 2 𝑏 superscript 𝑝 ′ dA^{\prime}/dp^{\prime}\approx a-2bp^{\prime} and for the contact area secant A ′ / p ′ ≈ a ′ − b ′ ​ p ′ superscript 𝐴 ′ superscript 𝑝 ′ superscript 𝑎 ′ superscript 𝑏 ′ superscript 𝑝 ′ A^{\prime}/p^{\prime}\approx a^{\prime}-b^{\prime}p^{\prime} found by the least square fit of numerical results in the interval p ′ ∈ [ 0.006 , 0.06 ] superscript 𝑝 ′ 0.006 0.06 p^{\prime}\in[0.006,0.06] identified for H = 0.8 𝐻 0.8 H=0.8 .

$L/{\lambda }_l$	$L/{\lambda }_s$	$\alpha \pm \mathrm{std}\left(\alpha \right)$	$b$	$a$	$b^{\prime }$	$a^{\prime }$
2	16	5.84 $\pm$ 0.35	3.244	2.331	4.63	2.421
2	32	13.23 $\pm$ 0.89	2.619	2.213	3.786	2.286
2	64	32.45 $\pm$ 2.72	2.201	2.144	3.104	2.199
2	128	84.83 $\pm$ 8.0	1.971	2.088	2.78	2.136
2	256	231.58 $\pm$ 24.32	1.605	2.032	2.205	2.067
2	512	651.48 $\pm$ 73.05	1.261	1.981	1.82	2.014
4	16	2.82 $\pm$ 0.11	2.896	2.318	4.236	2.391
4	32	5.51 $\pm$ 0.26	1.952	2.215	3.468	2.306
4	64	12.26 $\pm$ 0.55	1.802	2.152	3.126	2.231
4	128	30.05 $\pm$ 1.33	1.744	2.115	2.437	2.157
4	256	78.58 $\pm$ 3.82	1.474	2.057	2.081	2.094
4	512	214.15 $\pm$ 10.83	1.245	2.008	1.723	2.037
8	16	1.76 $\pm$ 0.01	3.325	2.454	3.383	2.452
8	32	2.75 $\pm$ 0.04	2.202	2.312	2.963	2.354
8	64	5.26 $\pm$ 0.13	1.955	2.237	2.561	2.276
8	128	11.69 $\pm$ 0.32	1.883	2.18	2.541	2.221
8	256	28.58 $\pm$ 0.77	1.772	2.127	2.181	2.152
8	512	74.54 $\pm$ 2.03	1.551	2.067	1.909	2.089
16	32	1.75 $\pm$ 0.01	2.938	2.441	3.012	2.447
16	64	2.73 $\pm$ 0.02	2.546	2.343	2.821	2.361
16	128	5.26 $\pm$ 0.07	2.437	2.271	3.063	2.307
16	256	11.71 $\pm$ 0.18	2.086	2.188	2.801	2.232

Table 3. Table 3: Sampled numerical results for the contact area fraction A ′ superscript 𝐴 ′ A^{\prime} growth with the normalized nominal pressure p ′ = p 0 / ( 2 ​ m 2 ​ E ∗ ) superscript 𝑝 ′ subscript 𝑝 0 2 subscript 𝑚 2 superscript 𝐸 p^{\prime}=p_{0}/(\sqrt{2m_{2}}E^{*}) for surfaces with H = 0.4 𝐻 0.4 H=0.4 and different cutoffs in the spectrum, mean results over 10 realizations is shown as well as the standard deviation denoted std ( A ′ ) superscript 𝐴 ′ (A^{\prime}) . The contact area is corrected using ( 1 ).

	$L/{\lambda }_s=4,L/{\lambda }_l=32$		$L/{\lambda }_s=8,L/{\lambda }_l=128$		$L/{\lambda }_s=16,L/{\lambda }_l=512$		$L/{\lambda }_s=4,L/{\lambda }_l=512$
	$\alpha =4.45\pm 0.15$		$\alpha =7.46\pm 0.14$		$\alpha =13.19\pm 0.20$		$\alpha =41.87\pm 1.41$
$p^{\prime }$	$A^{\prime }$	std$\left(A^{\prime }\right)$	$A^{\prime }$	std$\left(A^{\prime }\right)$	$A^{\prime }$	std$\left(A^{\prime }\right)$	$A^{\prime }$	std$\left(A^{\prime }\right)$
0.001	0.245	0.031	0.232	0.007	0.224	0.002	0.216	0.003
0.004	0.942	0.053	0.912	0.009	0.883	0.006	0.854	0.006
0.007	1.628	0.066	1.579	0.011	1.530	0.009	1.486	0.006
0.010	2.299	0.076	2.234	0.015	2.172	0.009	2.114	0.008
0.013	2.962	0.075	2.887	0.018	2.806	0.011	2.738	0.011
0.016	3.621	0.070	3.533	0.021	3.437	0.013	3.360	0.012
0.019	4.273	0.061	4.175	0.026	4.065	0.012	3.978	0.012
0.022	4.919	0.067	4.812	0.030	4.686	0.012	4.593	0.013
0.025	5.563	0.073	5.447	0.032	5.304	0.013	5.209	0.013
0.028	6.203	0.079	6.079	0.034	5.918	0.014	5.818	0.012
0.031	6.840	0.090	6.704	0.036	6.531	0.015	6.425	0.014
0.034	7.475	0.102	7.328	0.040	7.142	0.015	7.029	0.014
0.037	8.110	0.113	7.950	0.042	7.747	0.017	7.631	0.014
0.040	8.745	0.121	8.567	0.042	8.352	0.019	8.231	0.014
0.043	9.377	0.131	9.180	0.043	8.952	0.019	8.828	0.015
0.046	10.003	0.137	9.789	0.043	9.551	0.019	9.423	0.015
0.049	10.621	0.145	10.396	0.044	10.148	0.021	10.014	0.015
0.052	11.235	0.152	11.000	0.045	10.743	0.021	10.604	0.016
0.055	11.845	0.159	11.601	0.047	11.333	0.022	11.192	0.016
0.058	12.453	0.165	12.198	0.048	11.921	0.024	11.775	0.018

Table 4. Table 4: Sampled numerical results for the contact area fraction A ′ superscript 𝐴 ′ A^{\prime} growth with the normalized nominal pressure p ′ = p 0 / ( 2 ​ m 2 ​ E ∗ ) superscript 𝑝 ′ subscript 𝑝 0 2 subscript 𝑚 2 superscript 𝐸 p^{\prime}=p_{0}/(\sqrt{2m_{2}}E^{*}) for surfaces with H = 0.8 𝐻 0.8 H=0.8 and different cutoffs in the spectrum, mean results over 10 realizations is shown as well as the standard deviation denoted std ( A ′ ) superscript 𝐴 ′ (A^{\prime}) . The contact area is corrected using ( 1 ).

	$L/{\lambda }_s=4,L/{\lambda }_l=32$		$L/{\lambda }_s=8,L/{\lambda }_l=128$		$L/{\lambda }_s=16,L/{\lambda }_l=512$		$L/{\lambda }_s=4,L/{\lambda }_l=512$
	$\alpha =5.51\pm 0.26$		$\alpha =11.69\pm 0.32$		$\alpha =28.67\pm 0.55$		$\alpha =214.16\pm 10.83$
$p^{\prime }$	$A^{\prime }$	std$\left(A^{\prime }\right)$	$A^{\prime }$	std$\left(A^{\prime }\right)$	$A^{\prime }$	std$\left(A^{\prime }\right)$	$A^{\prime }$	std$\left(A^{\prime }\right)$
0.001	0.250	0.029	0.228	0.006	0.221	0.004	0.207	0.005
0.004	0.933	0.068	0.902	0.015	0.870	0.006	0.817	0.015
0.007	1.615	0.103	1.562	0.020	1.512	0.012	1.427	0.020
0.010	2.295	0.104	2.209	0.029	2.142	0.014	2.029	0.023
0.013	2.959	0.102	2.852	0.033	2.769	0.015	2.627	0.024
0.016	3.617	0.094	3.495	0.036	3.391	0.016	3.220	0.025
0.019	4.264	0.089	4.133	0.040	4.008	0.016	3.810	0.029
0.022	4.905	0.091	4.762	0.044	4.620	0.016	4.396	0.032
0.025	5.535	0.098	5.386	0.049	5.228	0.017	4.981	0.033
0.028	6.158	0.106	6.008	0.056	5.835	0.021	5.563	0.035
0.031	6.778	0.120	6.626	0.063	6.437	0.021	6.141	0.036
0.034	7.400	0.133	7.243	0.068	7.034	0.021	6.717	0.038
0.037	8.020	0.146	7.856	0.071	7.629	0.023	7.293	0.039
0.040	8.637	0.153	8.467	0.072	8.221	0.024	7.864	0.039
0.043	9.253	0.157	9.072	0.070	8.811	0.024	8.435	0.040
0.046	9.867	0.162	9.676	0.071	9.398	0.025	9.005	0.041
0.049	10.477	0.172	10.277	0.071	9.981	0.026	9.571	0.039
0.052	11.085	0.181	10.875	0.072	10.562	0.028	10.137	0.039
0.055	11.691	0.188	11.470	0.073	11.137	0.031	10.701	0.040
0.058	12.293	0.190	12.061	0.075	11.713	0.032	11.264	0.041

Equations

$$A_{\ast }=A_{\text{mbox}sim}-\frac{\pi -1+\ln 2}{24}S\Delta x.$$

$$⟨x⟩=\frac{1}{N^2}\sum _{i=0}^{N-1}\sum _{j=0}^{N-1}{x}_{ij}.$$

$${\hat{w}}_{mn}=\sum _{i=0}^{N-1}\sum _{j=0}^{N-1}w(x_i,y_i)\exp \left[-2\pi \mathrm{i}(m{x}_i+n{y}_j)/L\right],$$

$${\hat{f}}_{ij}=\hat{f}(K_i,K_j)=\{\begin{array}{l}{\left[\frac{K_i^2+K_j^2}{k_l^2}\right]}^{-(1+H)/2},\text{mbox}for1\le \frac{\sqrt{K_i^2+K_j^2}}{k_l}\le \zeta \\ 0,\text{mbox}elsewhere,\end{array}$$

$${\hat{z}}_{ij}=\hat{z}(k_x,k_y)=\mathrm{Re}({\hat{f}}_{ij})\left[\mathrm{Re}({\hat{w}}_{ij})+\mathrm{i}\mathrm{Im}({\hat{w}}_{ij})\right].$$

$$z(x_i,y_j)=\sum _{n=0}^{N-1}\sum _{m=0}^{N-1}{\hat{z}}_{nm}\exp [\mathrm{i}2\pi (n{x}_i+m{y}_j)/L].$$

$$R(\Delta x,\Delta y)=\frac{1}{N^2}\sum _{i=0}^{N-1}\sum _{j=0}^{N-1}z(x_i+\Delta x,y_j+\Delta y)z(x_i,y_j),$$

$$\Phi (k_x,k_y)=\hat{R}(k_x,k_y)=\mathrm{DFT}\left[z(x+\Delta x,y+\Delta y)\ast z(x,y)\right],$$

$$\Phi (k_x,k_y)=\hat{z}(k_x,k_y){\hat{z}}^{\ast }(k_x,k_y)={\hat{f}}^2(k_x,k_y)\hat{w}(k_x,k_y){\hat{w}}^{\ast }(k_x,k_y)={\hat{f}}^2(k_x,k_y){\hat{w}}^2(k_x,k_y).$$

$$⟨\Phi (k_x,k_y)⟩=⟨{\hat{w}}^2(k_x,k_y)⟩{\hat{f}}^2(k_x,k_y)=\{\begin{array}{l}{\Phi }_0{\left[\frac{\sqrt{K_x^2+K_y^2}}{k_l}\right]}^{-2(1+H)},\text{mbox}for1\le \frac{\sqrt{K_x^2+K_y^2}}{k_l}\le \zeta \\ 0,\text{mbox}elsewhere,\end{array}$$

$$⟨\Phi (K)⟩=\{\begin{array}{ll}{\Phi }_0(K/k_l{)}^{-2(1+H)},& \text{mbox}if1\le K/k_l\le \zeta \\ 0,& \text{mbox}otherwise.\end{array}$$

$$m_{pq}={\left[\frac{2\pi }{L}\right]}^{p+q}\sum _{i=0}^{N-1}\sum _{j=0}^{N-1}{i}^p{j}^q\Phi (2\pi i/L,2\pi j/L),$$

$$m_{0p}\approx m_{p0}\approx {\Phi }_0\underset{k_l}{\overset{k_s}{\int }}\underset{0}{\overset{2\pi }{\int }}{\left[k\cos (\phi )\right]}^p(k/k_l{)}^{-2(1+H)}kdkd\phi ={\Phi }_0{k}_l^{p+2}\frac{{\zeta }^{p-2H}-1}{p-2H}T(p),$$

$$T(p)=\underset{0}{\overset{2\pi }{\int }}{\cos }^p(\phi )d\phi =\{\begin{array}{ll}2\pi ,& \text{mbox}ifp=0;\\ \pi ,& \text{mbox}ifp=2;\\ 3\pi /4,& \text{mbox}ifp=4.\end{array}$$

$$m_2=\frac{m_{20}+m_{02}}{2},m_4=\frac{m_{40}+3m_{22}+m_{04}}{3}.$$

$$\alpha =m_0{m}_4/m_2^2$$

$$⟨(z-⟨z⟩{)}^2⟩=m_0,⟨(\nabla z-⟨\nabla z⟩{)}^2⟩=2m_2,⟨(\nabla \cdot \nabla z-⟨\nabla \cdot \nabla z⟩{)}^2⟩=m_4,$$

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

$$p^{\prime }=\frac{p_0}{E^{\ast }\sqrt{⟨|\nabla z^2|⟩}}$$

$$D=\frac{\sqrt{3}}{18\pi }\frac{m_4}{m_2},$$

$$S\sim L\sqrt{AD}.$$

$$S^{\prime }=S/(\sqrt{D}{L}^2)=\frac{S}{k_l{L}^2}\sqrt{\frac{24\pi }{\sqrt{3}}\cdot \frac{2-H}{1-H}\cdot \frac{{\zeta }^{2-2H}-1}{{\zeta }^{4-2H}-1}}.$$

$$A^{\prime }=\frac{\sqrt{\alpha }}{6}\underset{\hat{g}=0}{\overset{\infty }{\int }}\underset{\xi =\hat{s}}{\overset{\infty }{\int }}P(\xi ,\hat{g})\frac{(\xi -\hat{s})}{\hat{g}}d\xi d\hat{g},$$

$$p_0=E^{\ast }\sqrt{m_2/\pi }\frac{2{\alpha }^{3/4}}{9\sqrt{\pi }}\underset{\hat{g}=0}{\overset{\infty }{\int }}\underset{\xi =\hat{s}}{\overset{\infty }{\int }}P(\xi ,\hat{g})\frac{(\xi -\hat{s}{)}^{3/2}}{\sqrt{\hat{g}}}d\xi d\hat{g},$$

$$P(\xi ,\hat{g})=\frac{9}{2\sqrt{2\pi }}\sqrt{\frac{\alpha }{\alpha -1}}{\hat{g}}^3\exp \left[-\frac{\alpha }{2(\alpha -1)}{\xi }^2+\frac{3{\hat{g}}^2}{2}\right]\mathrm{erfc}\left(3\sqrt{\frac{(\alpha -1)}{2(2\alpha -3)}}\left\{\hat{g}+\frac{\sqrt{\alpha }}{3(\alpha -1)}\xi \right\}\right).$$

$$A_P^{\prime }/p^{\prime }=\mathrm{erf}(p^{\prime }\sqrt{2})/p^{\prime },d{A}_P^{\prime }/d{p}^{\prime }=\sqrt{8/\pi }\exp (-2p^{\prime 2})$$

$$A^{\prime }/p^{\prime }=a^{\prime }-b^{\prime }{p}^{\prime }.$$

$$\begin{array}{rl}& a^{\prime }(\alpha )=2.35-0.057\ln (\alpha -1.5),\\ & b^{\prime }(\alpha )=2.85-0.24\ln (\alpha -1.5),\end{array}$$

$$\begin{array}{rl}& a(\alpha )=2.31-0.057\ln (\alpha -1.5),\\ & b(\alpha )=2.17-0.2\ln (\alpha -1.5).\end{array}$$

$$A^{\prime }(p^{\prime },m_2,\alpha )=a^{\prime }(\alpha )p^{\prime }-b^{\prime }(\alpha )p^{\prime 2},\text{mbox}with{p}^{\prime }=\frac{p_0}{\sqrt{2m_2}{E}^{\ast }},$$

Full text

The role of the roughness spectral breadth

in elastic contact of 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, F 91003 Evry, France

bComputational Solid Mechanics Laboratory (LSMS, IIC-ENAC, IMX-STI),

Ecole Polytechnique Fédérale de Lausanne (EPFL), Bât. GC, Station 18, CH 1015 Lausanne, Switzerland*)

Abstract.

We study frictionless and non-adhesive contact between elastic half-spaces with self-affine surfaces. Using a recently suggested corrective technique, we ensure an unprecedented accuracy in computation of the true contact area evolution under increasing pressure. This accuracy enables us to draw conclusions on the role of the surface’s spectrum breadth (Nayak parameter) in the contact area evolution. We show that for a given normalized pressure, the contact area decreases logarithmically with the Nayak parameter. By linking the Nayak parameter with the Hurst exponent (or fractal dimension), we show the effect of the latter on the true contact area. This effect, undetectable for surfaces with poor spectral content, is quite strong for surfaces with rich spectra. Numerical results are compared with analytical models and other available numerical results. A phenomenological equation for the contact area growth is suggested with coefficients depending on the Nayak parameter. Using this equation, the pressure-dependent friction coefficient is deduced based on the adhesive theory of friction. Some observations on Persson’s model of rough contact, whose prediction does not depend on Nayak parameter, are reported. Overall, the paper provides a unifying picture of rough elastic contact and clarifies discrepancies between preceding results.

Keywords. roughness, contact area, Nayak parameter, spectrum breadth, pressure-dependent friction, Hurst exponent.

1 Introduction

Many engineering systems include components with contacts: rolling bearings, tire/road and wheel/rail, pieces assembled by bolts and rivets, gears, electric switchers, vehicle and aircraft brakes, NEMS and MEMS, etc. Macroscopic behavior of these components are often determined by complex electro-thermo-chemico-mechanical interactions at contact interfaces at small scales. For many materials and structures the mechanical behavior at such scales is microstructure-dependent and near the surface can differ significantly from in-bulk behavior due to surface tension, coatings/tribo-films and oxide layers. Moreover, in most metallic devices, near-surface layers are cold hardened and recrystallized giving considerably different plastic properties. In addition, in systems with high interface stresses and/or high temperatures due to friction or fracture dissipation, near-surface microstructure may change in operation. An interplay between complex thermo-mechanical behavior at small scale and surface roughness coupled to other relevant physics determine the macroscopic response of the system, its life cycle and failure modes. Thus, understanding mechanical behavior of rough surfaces in mechanical contact is a key point in understanding numerous tribological systems: macroscopic friction and wear laws but also mass and heat transport along and across contacts.

Regardless the fact that numerous analytical, semi-analytical and numerical models are available nowadays, a systematic and accurate analysis of the effect of roughness parameters on contact properties is still missing. In analytical and semi-analytical models strong assumptions, introduced to make models traceable, do not allow to make quantitative predictions beyond small intervals of model validity. Multiple models relying on roughness representation by a set of individual non-interacting asperities (or with first-order interaction included in model extensions) are available and termed here asperity-based models, pioneered by works of Archard [5] and Greenwood & Williamson [28], more recent models [11, 61, 43, 26, 15, 1] are based on random process theory of rough surfaces established in [38, 44]. The main drawback of these models is that elastic interaction (long-range) between asperities is missing. In reality, even if a truly rough surface could be represented by a set of asperities [29], every asperity coming in contact affects vertical positions of all the surface points as $~{}F/r$, where $F$ is the contact reaction force and $r$ is the distance from the center of asperity. Such an interaction is not trivial to handle in analytic models, yet some models including zero-order interaction were developed [18, 46]. Another drawback of asperity models is that they do not include possibility of merging contact area spots associated with different asperities, which becomes a bottle-neck limitation for surfaces with high asperity densities [27]. This drawback was also partly attenuated by including the possibility of merging and creating new macroscopic asperities [1]. However, to include this effect, the model should become deterministic loosing its stochastic nature, and thus it has to be integrated numerically, which requires to keep track of every individual asperity at given location and its interactions with all other asperities. Nevertheless, this model being sufficiently accurate is much faster than full numerical models based on boundary or finite element methods. Finally, the curvature of each individual asperity is assumed to be constant (parabolic surface), which cannot be ensured in smooth surfaces as the sign of curvature has to change on a path connecting neighboring asperity tips, like for example in a wavy surface.

Another class of models of rough contact termed as Persson’s model [51] is based on the evolution of contact pressure probability density (PDF) with increasing “magnification”, the latter controls the root mean squared surface gradient. This model, being exact at full contact (which would require infinite pressure for surfaces with continuous spectra, and thus, for the height distribution defined on an infinite support) was extended to partial contacts by a boundary condition for zero-pressure. However, the diffusion equation obtained by the author for the full contact and linking the PDF of pressure with the magnification, was not rigorously proved to hold at partial contacts [42, 20]. In addition, the author introduced a corrective function for the contact area [50], which should be a pressure or contact-area dependent function, to take into account differences in elastic energy for the case of partial contact compared to full contact. It is worth also mentioning fractal models of contact [40, 41] and discretely-continuous multi-level models with embedded long range elastic interaction between asperities [23, 21, 22], which both also rely on several approximations and strong simplifications of realistic roughness.

In this light, direct numerical simulations of rough contacts present a promising approach for solving non-linear contact problems between rough surfaces [31, 48, 65, 32, 12, 46, 71, 3, 67, 54, 47, 56, 53, 4, 69]. Different methods based on boundary or finite element methods are available. On the one hand they are free of the aforementioned assumptions and limitations, but on the other hand they are subject to numerical errors due to 1) possible inaccurate treatment of mechanical behavior or contact (a prominent example would be the penalty or barrier method for constraint satisfaction); 2) discretization error and mesh convergence rate. Even though the importance of accurate discretization in contact interface was identified some time ago [71, 55, 67, 47], a true mesh convergence study was missing with notable exceptions [13, 71, 57, 53]. In our previous work [69] we made an effort to estimate the discretization-induced numerical error in contact area computation: it was shown that it depends not only on the ratio $\Delta x/\lambda_{s}$, where $\Delta x$ is the distance between grid or mesh points and $\lambda_{s}$ is the shortest wavelength in the surface spectrum, but that the error also depends on the ratio $\lambda_{l}/L$, where $L$ is the length of the simulation box and $\lambda_{l}$ is the longest wavelength in the surface spectrum for spectra without plateau. The first ratio $\Delta x/\lambda_{s}$ should tend to zero to capture accurately the mechanical behavior of each asperity, thus at grid size the surface should remain smooth contrary to what was done in early numerical studies [31, 48, 46]. The second ratio $\lambda_{l}/L$ determines representativity of the surface [67, 69], i.e. the proximity of the surface height distribution to a Gaussian one, so the representative surface limit can be reached for $\lambda_{l}/L\to 0$, which is also a numerically unreachable limit. High values of this ratio (of the order of unity) thus require an extensive statistical study and averaging over many roughness realizations to extract accurate mean values, similar to what is done for homogenization of mechanical properties [36]. However, in many studies $\lambda_{l}/L=1$, i.e. the longest wavelength in the surface spectrum is equivalent to the simulation box (which is often assumed to be periodic), which implies a considerable dispersion in results and requires many simulations to obtain a statistically meaningful mean behavior. In [69] it was also shown that when the ratio $\lambda_{l}/L$ decreases the true contact area increases: this scenario appeared reasonable to the authors, as it was argued [67] that the periodicity coupled to long range elastic interactions and non-representativity of the system may affect significantly the results. What was inconsistent in [69] is that no convergence was observed even for quite representative surfaces with $\lambda_{l}/L=1/16$. In other words, the surface becomes more and more representative and isotropic but the true contact area continues to increase when $\lambda_{l}/L$ decreases further. This spurious behavior was accounted for discretization errors, and thus all numerical results contaminated by this error could not be considered ultimately correct.

In our recent study [70] we performed an accurately designed mesh convergence test and suggested a specific technique to correct the measure of the contact area. The technique is based on the error estimation suggested in [69] and some geometric arguments. The resulting equation for the true contact area includes the measured contact area in simulations $A_{\mbox{\tiny sim}}$ and the contact perimeter $S$, i.e. the number of switches between contact and non-contact points along vertical and horizontal lines of the grid, scaled with the grid size $\Delta x$. Then the true contact area can be estimated as

[TABLE]

This form provides a very accurate estimation of the contact area: even for poorly discretized surfaces $\Delta x/\lambda_{s}=1/2$ this equation provides as accurate estimation as a simulation carried out on a much finer mesh $\Delta x/\lambda_{s}=1/32$ for unchanged spectral content of the surface.

With Eq. (1) in hand, in this paper we carry out numerous numerical simulation and analyze the results to obtain contact area evolution under increasing pressure with unprecedented accuracy. Based on these results, several important conclusions are made concerning the role of the two roughness parameters: first, the spectral breadth (so-called Nayak parameter) and, second, the fractal dimension or equivalently the Hurst exponent. Both parameters are strongly linked, but surprisingly the role of the former one has not yet been thoroughly studied in the framework of full numerical models.

The paper is organized as follows: in Section 2 the numerical methods and roughness models are briefly outlined, in Section 3 numerical results are presented, namely the contact perimeter, contact area, inverse mean pressure and derivative of the contact area with respect to the nominal pressure. A simple phenomenological equation is suggested for the contact area evolution with universal coefficients depending on the Nayak parameter. Finally, the role of Nayak parameter in contact between rough surfaces is put in evidence. Implications of our results for Persson’s model of contact, which does not depend on Nayak parameter, is also discussed. The results are summarized and interpreted in Section 4. In Appendix some additional plots and tables are presented. All numerical data can be found in Supplemental material [66].

2 Methods

2.1 Synthetic rough surfaces

To generate self-affine rough surfaces, we use a filtering-in-Fourier-space technique [30]. Four parameters and a seed for a random number generator fully determine the roughness: $k_{l}=2\pi/\lambda_{l}$ and $k_{s}=2\pi/\lambda_{s}$ are wavenumbers determining the lower and higher cutoffs in the discrete spectrum. We deal with a periodic square surface of side $L$, thus the shortest and the longest cutoff wavelengths should be defined as $\lambda_{s}=L/n_{s}$,$\lambda_{l}=L/n_{l}$, where $n_{s},n_{l}\in\mathbb{N}$. The third parameter $\Phi_{0}$ is the amplitude scaling parameter and $H$ is the Hurst exponent. The method consists in the following. We generate a white noise $w_{ij}=w(x_{i},y_{j})$ on a grid with $N\times N$ points, such that $\langle w\rangle=0$ and $\langle w^{2}\rangle=\Phi_{0}$, where $\langle\bullet\rangle$ denotes the mean value over the domain:

[TABLE]

It is transformed in Fourier space (Discrete Fourier Transform):

[TABLE]

so that $\langle\hat{w}\hat{w}^{*}\rangle=\langle w^{2}\rangle=\Phi_{0}$, where $\hat{w}^{*}$ denotes complex conjugate of $\hat{w}$. In Fourier space we create a power-law decaying filter which retains wavenumbers only in the interval $k_{l}\leq k\leq k_{s}$ and scales their amplitude according to a decaying power-law:

[TABLE]

where $K_{i}=(s+1)\pi/L-2\pi si/L$, $K_{j}=(t+1)\pi/L-2\pi tj/L$, for $s,t\in\{-1,1\}$, $\zeta=k_{s}/k_{l}$ is referred as magnification in Persson’s model [51], and the parameter $H$ is the Hurst exponent determining the rate of decay, which in fractal description of 2D surfaces is linked to the fractal dimension $D_{f}$ as $H=3-D_{f}$. Note that imaginary part of the filter is everywhere zero $\mathrm{Im}{f}=0$. Next, we take a product of the white noise with the filter, which in Fourier space results in

[TABLE]

The field is transformed back to the real space to provide the targeted rough surface

[TABLE]

Note that the imaginary part of $z$ is zero. The discrete auto-correlation function is defined for the periodic $L\times L$ surface as

[TABLE]

and finally the power spectral density (PSD) of the roughness defined as the discrete Fourier transform (DFT) of the auto-correlation function, is given by:

[TABLE]

which according to the convolution theorem can be found by point-wise product of the Fourier transform of $z$ as

[TABLE]

Thus the PSD in average decays proportionally to $\hat{f}^{2}(k_{x},k_{y})$, because the random value $\hat{w}^{2}(k_{x},k_{y})$ in average does not depend on $k_{x},k_{y}$. It is also clear that for random realizations of the white noise the following statement is ensured:

[TABLE]

here $\langle\bullet\rangle$ denotes the average over multiple realizations of the white noise. It is the targeted form of the PSD, which follows the power-law-decay with the wavenumber with the exponent $-2(1+H)$. If the resulting surface is isotropic (spectral content is independent in statistical sense of the profile orientation on the surface), and it is supposed to be so if the filter is axi-symmetric and the $k_{l}$ parameter is big enough, then the PSD can be represented as a function of a single wavenumber $K=\sqrt{K_{x}^{2}+K_{y}^{2}}$:

[TABLE]

Spectral moments $m_{pq}$ for a surface with a discrete spectrum are computed as follows:

[TABLE]

under the assumption of isotropic surface and for a broad spectrum, spectral moments can be approximated by an integral [69]:

[TABLE]

where $T(p)$

[TABLE]

Isotropy of the surface, in the first approximation, can be verified by the following relationships between moments: $m_{2}=m_{20}=m_{02}$, $m_{4}=3m_{22}=m_{40}=m_{04}$. For a surface with a discrete spectrum, these equalities are satisfied approximately, thus to even out errors one can rather use the following expressions for the spectral moments:

[TABLE]

The Nayak parameter $\alpha\in[1.5,\infty]$ introduced in [44] determines the breadth of the surface spectrum and is given by the following expression:

[TABLE]

Below, we recall some basic definitions and relationships between spectral moments and the variances of surface heights, of surface gradient and of surface curvature:

[TABLE]

where mean values are computed over a surface period $L\times L$, i.e. $\langle\bullet\rangle=\frac{1}{L^{2}}\int\limits_{0}^{L}\int\limits_{0}^{L}\bullet\,dxdy$. It is important to note that in the case of a discrete spectrum, these relationships are not exact and are strongly dependent on the surface discretization and on the way how gradients and Laplacians are computed [24, 46]. Thus to obtain a discretization-independent measurement of Nayak parameter and of the root mean squared (rms) gradient, which is abbreviated hereinafter as $\sqrt{\langle|\nabla z|^{2}\rangle}$, assuming that $\langle\nabla z\rangle=\boldsymbol{0}$, we evaluate them through spectral moments (2) and use average value (4): $\sqrt{\langle|\nabla z|^{2}\rangle}=\sqrt{m_{02}+m_{20}}$.

The following combinations of parameters were considered in this study: $L/\lambda_{l}=Lk_{l}/2\pi=\{2,4,8,16\}$, $L/\lambda_{s}=Lk_{s}/2\pi=\{16,32,64,128,256,512\}$, $H=\{0.4,0.8\}$. For every combination of parameters (46 combinations excluding cases $k_{l}=k_{s}$) ten different roughness were generated. For all surfaces the rms surface gradient is kept constant $\sqrt{\langle|\nabla z|^{2}\rangle}=1$. All generated rough surfaces are periodic and are constructed on the grid $N\times N=2048\times 2048$. Some examples of generated rough surfaces, which are used in numerical simulations, are depicted in Fig. 1.

2.2 Comments on the roughness and contact models’ validity

Note that contrary to the zeroth moment (or the variance of height), spectral moments $m_{2}$ and $m_{4}$ diverge when $Lk_{s}\to\infty$. Thus their determination from experimental roughness measurements depend on the magnification and measurement sampling, i.e. on the grid step. Roughly speaking, the value of $m_{2}$ is measurement-dependent, and since $m_{2}$ is the main scaling parameter in rough contact models of elastic contact, the output of these models is also measurement-dependent. Since many surfaces demonstrate self-affine scaling and power-law decrease of the power spectral density down to interatomic distances [37], where continuum mechanics and especially contact interaction [9, 39, 8] are rigorously not applicable, then the study of elastic contact becomes questionable. Moreover, at atomic scale the definition of contact becomes ambiguous as solid-solid interaction is determined by interactions with infinite support: short-range interaction in absence of Coulomb electrostatic forces, and long-range interaction in their presence. Therefore to use all existing models of elastic rough contact one needs to assume that the surface is smooth below a certain wavelength $\lambda_{s}$, which certainly is quite a rough approximation. A different approach would be to consider that elastic interaction does not longer holds at shorter wavelengths and other scale-related interface features are dominant: atomic or capillary adhesion, surface energy, etc. Moreover, the cutoff problem is easy to address if elasto-plastic material behavior is considered for which too sharp asperities coming in contact are readily flattened out in fully plastic regime, for which the contact pressure saturates at material hardness [10]. However, not all materials exhibit similar-to-metals plastic behavior (e.g. pressure dependent plasticity in rocks or amorphous solids), moreover at smaller scales hardness becomes size-dependent [45, 60], which might bias the conclusions of models based on macroscopic hardness. In summary, for the time being the question of physically motivated high-frequency cutoff remains open and will not be addressed in this study.

2.3 Mechanical contact between rough surfaces

We use a spectral boundary element method proposed in [59], which is rather similar to [33], the method is based on the Kalker’s formulation of the minimization problem under constraints [35]. We solve the normal contact problem between a synthetic rough surface, assumed to be rigid and a linearly elastic half-space. The roughness of the rigid surface should be considered as an effective roughness of two contacting solids $z(x,y)=z_{1}(x,y)-z_{2}(x,y)$. The slopes are assumed to be very small: $|\nabla z|\ll 1$ so that tangential displacements are negligible and the classical approach [34], where the loads are assumed to be applied on a flat surface, is ensured. Only the normal, frictionless and non-adhesive contact is considered. The only material parameter used, is the effective elastic modulus given by:

[TABLE]

where $E_{i},\nu_{i}$ are Young’s modulus and Poisson’s ratio of material $i=1,2$, respectively.

The normalized nominal pressure applied at infinity is computed as

[TABLE]

and is increased from zero up to $p^{\prime}=2$ in 120 load steps; to ensure accurate results at small contacts, 60 load steps are used to reach $p^{\prime}=0.06$. We store spectral moments of the surface and for every load step we store the perimeter of contact clusters $S$, and the contact area $A_{\mbox{\tiny sim}}$. The perimeter is computed as a number of switches from contact to non-contact and vice versa along all vertical and horizontal lines on the grid multiplied by the discretization step $\Delta x$, the contact area is computed as the number of grid points in contact multiplied by the area corresponding to every point $\Delta x^{2}$. In the paper we present the true contact area fraction given by the ratio $A^{\prime}=A/A_{0}$, where $A_{0}=L^{2}=N^{2}\Delta x^{2}$.

3 Results and Discussions

3.1 Contact perimeter

Since the contact perimeter plays an important role in the accurate estimation of the true contact area [70], it is of interest to know how it changes under increasing load. Evolution of the contact perimeter (normalized and non-normalized data) with respect to the true contact area and the nominal pressure is presented in Figs. 2,3 for $H=0.4,0.8$, respectively. The normalization $S^{\prime}=S/(\sqrt{D}L^{2})$ is found from the following consideration: if the contact is formed by identical circles of radius $r$, then the total area is given by $A\sim\pi DL^{2}r^{2}$, where $L$ is the side length of the considered region and $D$ is the total asperity density given in Nayak’s paper [44]: which is given by

[TABLE]

which should be proportional to the density of contacting asperities $D_{c}$. At the same time the total perimeter of contact spots is given by $S\sim 2\pi DL^{2}r$. Expressing $r$ through $A$ and $D$, we obtain that

[TABLE]

It means that for equivalent contact area $A$, which depends rather weakly on parameters $k_{l}$ and $k_{s}$, the perimeter scales as $S\sim L^{2}\sqrt{D}$. Then expressing asperity density through $k_{l},k_{s}$ and $H$ using Eq. (3), we finally obtain the following normalization:

[TABLE]

As seen in Figs. 2,3, the perimeter is not symmetric with respect to $A^{\prime}=A/A_{0}=50$ %, demonstrating the asymmetry in evolution of contact and non-contact clusters (see also [68]). Interestingly, for the chosen normalization, the slope of the perimeter approaching full contact is almost the same for all roughnesses [Figs. 2(b),3(b)]. On the contrary, for all cutoffs the normalized perimeter evolves differently near zero contact. The maximum of the perimeter, which corresponds to the maximal absolute error222Note that the maximal relative error is always reached at minimal computed contact area, see Fig. [69, Fig. 12] in contact area computation, can be easily identified in Figs. 2,3. Note that the normalized perimeter decreases with Nayak parameter $\alpha$, whereas the area corresponding to the peak perimeter increases with increasing $\alpha$. For surfaces with high enough $\alpha$ the contact area corresponding to the maximum of the perimeter seems to converge toward $A^{\prime}=A/A_{0}\approx 57$ % [see Fig. 3(b)]. It could seem to be meaningful to link this contact area with the percolation area of contact clusters, which corresponds to the appearance of an infinite contact cluster. However, this value appears to be considerably higher than the percolation limits found i) for a bi-wavy surface $A^{\prime}_{p}\approx 40.2$ % [68] and ii) for random rough surfaces $A^{\prime}_{p}\approx 40$ % [2], $A_{p}\approx 42.5$ % [19], $A^{\prime}_{p}\approx 35.5$ % [57], a broader range of percolation limit was reported in [58].

Another remark on the perimeter deals with a dimensionless characteristic of the shape, compactness [68], which can be defined as a ratio of the perimeter to square root of the contact area: $C=S/\sqrt{A}$. In the simplest case of identical circular contact spots, following (5) the compactness is a constant proportional to the square root of asperity density $S/\sqrt{A}\sim L\sqrt{D}$, however, in a realistic rough contact it is not the case. Evolution of the compactness for different roughnesses is shown in Fig. 4 for both $H=0.4,0.8$. Note that if normalized perimeter is used to determine the compactness, as it is done in Fig. 4, $C^{\prime}=S/\sqrt{DA^{\prime}}L^{2}$, then the normalization $DA^{\prime}$ can be interpreted as the density of contacting asperities, whereas $D$ is the total density of asperities. Interestingly, the area value corresponding to the maximum of the compactness also increases with the Nayak parameter $\alpha$ and saturates at $A/A_{0}\approx 42$ %, which is intriguingly close to the above-mentioned percolation limits found by other authors, for instance $A_{p}\approx 42.5$ % from [19]. In this paper we limit ourselves to this short remark on the link found between the percolation limit and the contact area at which the maximal value of compactness is reached. A further study is needed to confirm this correlation and to find a rigorous explanation for it.

3.2 Contact area

Here, we present numerical results of the contact area evolution computed for surfaces with different spectra. Some numerical data are presented in Tables 3,4 in the Appendix, all data can be found in supplementary material [66], which we provide for reader’s convenience. All results were post-processed using the suggested area-correcting technique Eq. (1) (see [70]). In Figs. 5,6 the true contact area evolution with normalized nominal pressure is shown for $H=0.4,0.8$, respectively.

Looking on raw non-corrected results Figs. 5(a),6(a), it can be concluded that Nayak parameter $\alpha$ has a non-monotonous effect on the contact area, which contradicts the prediction of all asperity-based models [16] that higher $\alpha$ results in a smaller contact area for the same normalized pressure $p^{\prime}$. Compare, for example, in Fig. 5 results obtained for $\alpha=1.75$, $\alpha=13.52$ and $\alpha=41.86$, the corresponding contact areas found at $p^{\prime}\approx 0.06$ first decreases and next increases with increasing $\alpha$. It also could seem from the raw data that the lower cutoff $L/\lambda_{l}$ plays an important role in contact area determination. Compare, for example, in Fig. 6(a) surfaces $L/\lambda_{l}=16,L/\lambda_{s}=128$ and $L/\lambda_{l}=8,L/\lambda_{s}=64$ which have the same $\alpha=5.26$ and the same magnification $\zeta=k_{s}/k_{l}$, but different lower frequency cutoff $\lambda_{l}$.

After the correction (Figs. 5(b),6(b)], clearly, this spurious dependence of results on the lower cutoff in the surface spectrum $\lambda_{l}$ and on the Nayak parameter disappears. The former was already identified in [69] (see also Fig. 17 in Appendix for post-processed results of [69]). In corrected results, close values of Nayak parameter result in close values of the contact area: compare, for example, in Fig. 6(a,b) $\alpha=5.26$ for $L/\lambda_{l}=16,L/\lambda_{s}=128$ and $L/\lambda_{l}=8,L/\lambda_{s}=64$, or $\alpha=28.6,L/\lambda_{l}=8,L/\lambda_{s}=256$ and $\alpha=30.05,L/\lambda_{l}=4,L/\lambda_{s}=128$. Complementary figures 18,19 demonstrating the effect of the correction can be found in Appendix. Simulation results are also compared with analytical models: asymptotic solution of Nayak’s based asperity models $A^{\prime}=\sqrt{2\pi}p^{\prime}$ [11, 43, 26, 15, 61], original Persson’s model $A^{\prime}=\mathrm{erf}(p^{\prime}\sqrt{2})$ [51], and numerically integrated Greenwood’s simplified elliptic model [26] for relevant values of Nayak parameter. We recall that in the later model the contact area and the nominal pressure are defined as functions of dimensionless separation between the effective rough surface and a rigid flat $\hat{s}=s/\!\!\sqrt{m_{0}}$ as follows:

[TABLE]

[TABLE]

where $\hat{g}=g/\!\!\sqrt{m_{4}}$ is the normalized geometric mean curvature. The joint probability density $P(\xi,\hat{g})$ of asperities with the summit at normalized elevation $\xi=z/\!\!\sqrt{m_{0}}$ and the normalized mean geometric curvature $\hat{g}$ was obtained [26] to be

[TABLE]

From Figs. 5(b),6(b) it follows that for a given pressure, the true contact area for a surface with lower value of $\alpha$ is bigger than this for a surface with bigger $\alpha$, i.e. the contact area decreases with increasing value of Nayak parameter. This is also predicted by asperity-based models [26, 16]. Second, it is clear that the effect of $\alpha$ in asperity based models is exaggerated, and that accurate simulations predict a softer dependence on this parameter. On the other hand, Persson’s model, which is independent of $\alpha$, strongly underestimates the true contact area. As was already reported [67, 69] the dispersion of results (see error bars in Figs.5,6) is higher for smaller $k_{l}$, smaller $k_{s}/k_{l}$ and bigger $H$; for the latter, compare Figs. 18 and 19 in Appendix. In all following figures only corrected contact area is used.

The normalized contact area derivative with respect to the normalized nominal pressure $dA^{\prime}/dp^{\prime}$ is presented in Fig. 7 for $H=0.4$ and $H=0.8$. The normalized inverse mean contact pressure $A^{\prime}/p^{\prime}$ is presented in Fig. 8 for the same Hurst exponents. These results are compared with i) asymptotic prediction of Nayak’s based asperity model: $\lim_{A^{\prime}\to 0}(dA^{\prime}/dp^{\prime})=\lim_{A^{\prime}\to 0}(A^{\prime}/p^{\prime})=\sqrt{2\pi}$, ii) Persson’s model:

[TABLE]

and iii) with the simplified elliptic model [26], which was numerically integrated for corresponding Nayak parameters. As seen from the figures, the contact area slope in the interval $p^{\prime}\in[0,0.06]$ and $A^{\prime}\in[0,\approx 13\%]$ decreases by approximately $10-16$ %, i.e. the contact area can be nowhere approximated by an affine function of pressure. For a given pressure the slope value decreases with increasing $\alpha$. Simulation results show higher slopes and softer dependence on Nayak parameter than all asperity-based models, among which the simplified elliptic model has the softest dependence on $\alpha$ [16], that was the reason we picked it up for our comparison. Persson’s model predicts very low slope and its decrease with increasing pressure in the interval $p^{\prime}\in[0,0.06]$ is also considerably smaller than what is numerically obtained. Results of the inverse normalized mean pressure (the secant of the area-pressure curve) $A^{\prime}/p^{\prime}$, which are presented in Fig. 8, demonstrate qualitatively the same results as for the slope, the same conclusions hold.

It could be also remarked that qualitatively the simulation results are consistent with predictions of the simplified elliptic model [26] computed for small $\alpha$ as was also shown in [70]: the slope decreases almost linearly starting from very small pressures, i.e. roughly it behaves as $dA^{\prime}/dp^{\prime}\sim a(\alpha)-2b(\alpha)p^{\prime}$. Hence, by integrating it, the contact area evolution can be approximated by a second order polynomial without zero-order term: $A^{\prime}=ap^{\prime}-bp^{\prime 2}$. As seen in Fig. 8, the contact-area secant can be also approximated by an affine function of pressure with approximately the same constants $a^{\prime}$ and $b^{\prime}$:

[TABLE]

In Fig. 7,8 these linear fits are also plotted. If the contact area evolution would really be simply the second order polynomial, then $a^{\prime}=a$ and $b^{\prime}=b$. Comparison of these coefficients found by the least squares fit in the interval $p^{\prime}\in[0.006,0.06]$ does demonstrate that $a^{\prime}\approx a$, but shows a considerable difference in coefficients $b$ and $b^{\prime}$ (see Fig. 9 and Tables 1,2 in Appendix). All coefficients can be approximated by a monotonically decreasing logarithmic functions of Nayak parameter, which should be rather universal, thus we give here the identified numerical values (see fit curves in Fig. 9):

[TABLE]

[TABLE]

The coefficients for the contact-area secant (12) $a^{\prime},b^{\prime}$ were identified for $\alpha<600$, but could be extrapolated up to $\alpha\leq 1.5+\exp(2.85/0.24)\approx 143\,600$ where $b^{\prime}$ turns zero, the following increase in $\alpha$ would result in a negative $b^{\prime}$, which is probably nonphysical. It can be seen in Fig. 9 that for the linear coefficients $a,a^{\prime}$ no significant difference between surfaces with different Hurst exponent can be seen, whereas for the quadratic coefficients $b,b^{\prime}$ the fit parameters seemingly cluster by the Hurst exponent. An additional study including more surfaces with broader spectrum is needed to draw more accurate conclusions.

Finally, in practice, for a known nominal pressure in elastic contact, the real contact area fraction can be accurately estimated at least in the interval $A^{\prime}\in[1,15]$ % by the following phenomenological formula:

[TABLE]

which requires knowing the second spectral moment $m_{2}$, Nayak parameter $\alpha$ of the effective roughness and the effective elastic modulus of contacting materials.

3.3 Pressure-dependent coefficient of friction

Knowing the real contact-area evolution with pressure, one can readily improve an estimation of the static coefficient of friction based on the adhesive theory of friction [10], in which the true contact interface can bear the maximum shear traction $\sigma_{t}$ independently of pressure. Assuming that at mesoscale a part of the nominally flat frictional interface is subject to pressure $p_{0}$, the friction coefficient $\mu$ can be defined as the maximal value of the mesoscopic shear traction $t^{\max}$, which the interface can sustain before sliding starts, normalized by the contact pressure $p_{0}$:

[TABLE]

which can be represented as

[TABLE]

where $\mu_{0}=a^{\prime}(\alpha)\sigma_{t}/E^{*}\sqrt{2m_{2}}$ is pressure-independent coefficient of friction and $\mu_{0}b^{\prime}(\alpha)/a^{\prime}(\alpha)\sqrt{2m_{2}}$ corresponds to the non-linear friction part, which is proportional to $\mu_{0}$ with a factor depending only on roughness characteristics $m_{2}$ and $\alpha$. Alternatively, using (12) the pressure-dependent coefficient of friction can be estimated as:

[TABLE]

with $\beta=[2.85-0.24\ln(\alpha-1.5)]\;/[2.35-0.057\ln(\alpha-1.5)]$, which can be estimated as follows: $\beta\approx 0.84$ for $\alpha=10^{2}$, $\beta\approx 0.61$ for $\alpha=10^{3}$ and $\beta\approx 0.35$ for $\alpha=10^{4}$.

3.4 Role of Nayak parameter

Change of the inverse value of the normalized mean contact pressure $A^{\prime}/p^{\prime}$ (the secant of the normalized contact-area/pressure curve) with respect to Nayak parameter $\alpha$ was computed for different normalized nominal pressures $p^{\prime}=0.005,0.02,0.04,0.06$ is presented in Fig. 10(a,b,c,d), respectively. These nominal pressure values roughly correspond to contact areas $A^{\prime}\approx 1,5,9,13$ %, respectively, computed for a surface with $\alpha\approx 1.75$. The value of the secant decays with Nayak parameter following for each pressure a logarithmic curve independently of the Hurst exponent, which is well approximated by:

[TABLE]

with constants $d(p^{\prime}),c(p^{\prime})$ which can be expressed through the logarithmic approximations of $a^{\prime},b^{\prime}$ (12):

[TABLE]

Note that points corresponding to $L/\lambda_{l}=2,\;L/\lambda_{s}=16$, which are not representative and have a rather poor spectrum are located quite far from the master curve. The remaining points lie very close to (15) and we expect that the for higher Nayak parameters this logarithmic trend would be preserved.

The contact area secant $A^{\prime}/p^{\prime}$ evaluated at nominal pressure $p^{\prime}=0.02$ is compared with numerical results reported in [46] and with predictions of the simplified elliptic model evaluated for comparable pressures (Fig. 11). Asymptotic solution for asperity based model [11, 16] and Persson’s solution, both independent of Nayak parameter, are also plotted for reference purpose. The inverse normalized mean pressure changes very slightly with Nayak parameter $\alpha$, nevertheless this change takes place and is well defined. Contrary to what was argued in [11], realistic rough surfaces posses very broad spectra, hence $\alpha$ can be very high [16].

The discrepancy between our results and results obtained in [46] (also presented in Fig. 11) comes from the fact that in the latter the rms gradient was underestimated because it was computed directly on the roughly discretized geometry in the real space:

[TABLE]

contrary to the more accurate spectral based computation used here $\sqrt{\langle|\nabla z|^{2}\rangle}=\sqrt{m_{20}+m_{02}}$, which for a given spectrum does not depend on discretization as far as there are enough modes to represent the spectrum. The same remark concerns the second and the fourth spectral moments needed to determine the Nayak parameter. In addition, no correction of the contact area was used in [46], and since the authors considered roughnesses with shortest wavelength equivalent to a double spacing between grid points ($\lambda_{s}=2\Delta x$), the contact area was overestimated. In this sense, the results from [46] are rather similar to those reported in [31], where a too strong dependence of the secant $A^{\prime}/p^{\prime}$ on the Hurst exponent was found (see Fig. 14). In this light, our results obtained with an accurately evaluated rms gradient and Nayak parameter, and also using the corrective procedure for the contact area computation [70] approach better the true solution than the data points extracted from [46].

3.5 Effect of Nayak parameter vs Hurst exponent

In the considered interval of loads the inverse normalized mean pressure (contact area secant, see Fig. 11) lies in the interval between Persson’s solution and asperity based asymptotic solution as it was already shown by many authors for models of various accuracy and at different pressures [31, 13, 54, 53]. However, all these studies report the contact area secant as a function of the Hurst exponent. To the best of our knowledge, there is a unique study [46] in which this secant was plotted with respect to the Nayak parameter (see open points in Fig. 11). Hurst exponent of rough surfaces is an independent parameter, but when it changes, the Nayak parameter changes in response. And thus, in all reported studies the change of the secant as a function of the Hurst exponent was probably induced by the change in the Nayak parameter. The link between the Hurst exponent and the Nayak parameter for surfaces without plateau can be expressed as follows [69]:

[TABLE]

which is a monotonically increasing function of $H$ in the interval $[0,1]$, i.e. higher $H$ corresponds to higher $\alpha$. In Fig. 12(a) the link between the Hurst exponent and the Nayak parameter is shown. The theoretical prediction (17) is compared with Nayak parameters evaluated through spectral moments on numerous randomly generated rough surfaces with discretization $4096\times 4096$ and parameters $L/\lambda_{l}=6$, $L/\lambda_{s}=\{192,512,2048\}$ corresponding to $\zeta=\{32,85.33,341.33\}$ and for changing $H\in[0,1]$; for every value of $H$, 20 surfaces were generated. A good agreement between Eq. (17) and the direct evaluation of the Nayak parameter is observed. The second and the zero-th spectral moments $m_{2}$ and $m_{0}$ take the form [69]:

[TABLE]

where $\Phi_{0}=\Phi(k_{l})$ is the value of PSD calculated at $k_{l}$. These estimations are in good agreement with the spectral moments directly evaluated on similarly generated random rough surfaces: the ratio $m_{2}/(m_{0}k_{l}^{2})$ is plotted in Fig. 12(b).

Using Eq. (17), it can be shown that at extreme values $H=0$ and $H=1$, parameter $\alpha$ takes the following values:

[TABLE]

[TABLE]

Therefore, the change in $\alpha$ over the entire range of allowed Hurst exponents $0\leq H\leq 1$ can be expressed as

[TABLE]

For small values of $\zeta$ the change in $\alpha$ with $H$ is very small, thus the change in area $A^{\prime}$ with $H$ is also small. This can be seen in Fig. 13. That was the reason that led authors of [54, 67] to an inaccurate conclusion that contact area is independent (or very weakly dependent) of the Hurst exponent. See Fig. 14 for comparison of results from [54] with our phenomenological results based on numerical simulations. In reality, the contact area depends on the Hurst exponent through Nayak parameter, but this dependence is very sensitive to the magnification $\zeta$: for small $\zeta$ the role of $H$ on the contact area is hardly detectable, but for high $\zeta$ the dependence is significant.

In Fig. 13(a) we plot the inverse normalized mean pressure $A^{\prime}/p^{\prime}$ for different values of the Hurst exponent $H\in[0,1]$ and for different magnifications $\zeta=k_{s}/k_{l}=\{4,16,64,256,1024,2048,4096\}$ for two normalized nominal pressures $p^{\prime}=0.005,0.06$. For this, we used Eq. (15),(16) with $\alpha$ expressed through the Hurst exponent via Eq. (17):

[TABLE]

If $p^{\prime}$ is expressed through the rms gradient (square root of doubled $m_{2}$) and, next, through the Hurst exponent and the magnification $\zeta$ using Eq. (18):

[TABLE]

then the true contact area can be found by substituting $p^{\prime}$ in form (20) in (19):

[TABLE]

where the coefficients can be found using (16):

[TABLE]

[TABLE]

The true contact areas computed for different pressures $p_{0}/E^{*}=0.0002,0.001$ and fixed $\Phi_{0}k_{l}^{2}=10^{-6}$ are plotted in Fig. 13(b) as functions of the Hurst exponent for different magnifications $\zeta$.

When reformulated in terms of the Hurst exponent and magnification, our phenomenological results based on very accurate and statistically sound numerical simulations can be compared with similar results of other authors (Fig. 14). Our data is plotted as bounds on the secant evaluated at nominal pressures $p^{\prime}=0.005,0.06$ and computed for surfaces with different magnifications $\zeta$. In one of the pioneering works in numerical analysis of rough contact [31, Fig. 7] it was shown that the contact-area secant $A^{\prime}/p^{\prime}$ (and thus the contact area) evaluated at certain nominal pressure, decreases with increasing Hurst exponent. It was shown that for $H=0.3$ it approached closely the asymptotic limit of asperity models and for $H=0.9$ it was just above Persson’s asymptotic value. Other authors reproduced this study and plotted the secant as a function of the Hurst exponent [13, 52, 54, 53]. As already mentioned, in some studies [54, 67], however, it was argued that for properly discretized surfaces the secant does not depend (or depends very weakly) on the Hurst exponent. As seen from our results, indeed the dependence on the Hurst exponent is the weakest for $H$ approaching unity. In this light, the numerical results from [54] show a reasonable agreement with our predictions. Results from [13, 53] seem to overestimate the contact area as no correction technique was used, which is crucial for accurate numerical simulations, especially for the high magnifications $\zeta=1024$ used. Finally, the results from [31] suffer both from inaccurate estimation of rms gradient, too coarse mesh and the absence of area correction.

3.6 From infinitesimal to full contact

Contact area evolution up to full contact is depicted in Fig. 15 and compared with Persson’s prediction (10). Qualitatively Persson’s model describes properly the evolution of the true contact area but quantitatively, simulations predict higher values of area and thus full contact occurs earlier in simulations 333Since Persson’s model deals with a continuous spectrum and thus with a truly Gaussian surface with an infinite support, full contact happens only at infinite nominal pressure, but since error function tends very rapidly to unity, a certain threshold can be used to identify full contact, i.e. it could be $A^{\prime}_{\mbox{\tiny th}}=1-1/N^{2}$, which corresponds to the last element coming in contact, where $N$ is the number of grid points per side.. Near full contact, shown in the inset in Fig. 15, all results collapse in a single master curve, which is well described by the queue of an error function.

The results near full contact can be of interest in metal forming applications [7, 62, 6] but also in theoretical findings [25, 63, 17]. However, near full contact our results have to be handled with precaution, because there, the area correction technique [70] might work less properly. Due to insufficient discretization if a non-contact patch, surrounded by a contact zone, decreases in size under increasing load, it will soon disappear when reached the diameter of one-two elements (in raw simulations the contact area is overestimated). When such a patch disappears, the associated contact area cannot be corrected because there is no perimeter associated with this patch.

3.7 Implications for Persson’s model

Persson’s model, derived for full contact under infinite nominal pressure, deals with the probability density of contact pressures $P(p)$, which takes the forms of a diffusion equation, where $P$ acts as the concentration of diffusing quantity and the variance of the surface gradient $m_{2}$ acts as time. In a simplified form obtained in [42], it writes as:

[TABLE]

This interesting result, was extended to partial contacts and thus finite pressures, by imposing the probability density of zero pressure to be zero, i.e. $P(0,m_{2})=0$. This extension, however, is the weakest point of the model as this transition from full to partial contact is not justified [42]. Numerous comparisons of Persson’s model with numerical simulations [31, 48, 13, 53, 14, 54, 67, 69], especially see [20], revealed that the contact area predicted by Persson’s model is underestimated. Meanwhile, Persson argued that introducing the boundary condition is not sufficient for extending his model to partial contacts, but that it is also needed to take into account the interfacial elastic energy change due to detachment [50, 64]. It was argued that the “scaling factor” needed to adjust the energy is of “order of unity”. In Fig. 16 we plot the ratio of the numerically predicted contact area to the prediction of the original Persson’s model. Three surfaces are considered $L/\lambda_{l}=4$, $H=0.8$ and $L/\lambda_{s}=\{32,128,512\}$ corresponding to $\alpha\approx\{5.5,30,214\}$. This ratio evolves non-linearly with increasing pressure and at the same time it strongly depends on Nayak parameter.

An approximate form of this ratio takes the following form:

[TABLE]

and could be used as an alternative phenomenological correction [50, 49]. Note that as $p^{\prime}\to\infty$ then $f(p^{\prime})\to 1$, and as $p^{\prime}\to 0$ then $f(p^{\prime})\to\pi/2$, i.e. $f(p^{\prime})$ tends to the ratio between asperity models’ asymptote to the Persson’s asymptote. Eq. (22) includes two exponential decays: the first one, decaying with the square root of pressure, is relevant to the interval of small pressure, especially for surfaces with high Nayak parameters, the second decay partly balanced by a power-law with an exponent $3/2$ is responsible for the long range decay up to full contact. This non-trivial two-stage correction demonstrates that the “factor of order of unity” is not simply a factor but a complex function depending on Nayak parameter, whose introduction would be hard to justify. In Fig. 16 the following coefficients were used to fit the ratio of the simulated true contact area to the Persson’s prediction:

• •

$L/\lambda_{s}=32$, $\alpha\approx 5.5$: $c_{1}=0.905$, $c_{2}=15$, $c_{3}=5.7$, $c_{4}=6.5$;

• •

$L/\lambda_{s}=128$, $\alpha\approx 30$: $c_{1}=0.855$, $c_{2}=30$, $c_{3}=5.2$, $c_{4}=5.0$;

• •

$L/\lambda_{s}=512$, $\alpha\approx 214$: $c_{1}=0.815$, $c_{2}=50$, $c_{3}=5.2$, $c_{4}=5.4$.

all three roughnesses have $L/\lambda_{l}=4$ and $H=0.8$.

4 Conclusion

By carrying out statistically sound high accuracy numerical simulations of mechanical contact between rough surfaces with determined spectral content, we demonstrated how Nayak parameter $\alpha$ affects the true contact area evolution with increasing nominal pressure. We showed that in the interval $\alpha\in[1.75,\;\approx 600]$, an increase in $\alpha$ results in a decrease in the contact area fraction $A^{\prime}$. This trend was already predicted by advanced asperity models based on Nayak’s random process model, but due to inherent approximations of asperity models, this dependence was exaggerated. On the contrary, the original Persson’s model of rough contact does not take consideration of Nayak parameter. We deduced that the contact area for a given normalized pressure decreases logarithmically with Nayak parameter.

In several previous studies of the rough contact, attempts ware undertaken to deduce the role of the Hurst exponent (or fractal dimension) in the true contact area. In some studies it was argued that the contact area might be independent or too weakly dependent on the Hurst exponent. Here we show clearly the link between the Nayak parameter and the Hurst exponent, and also we demonstrate that for surfaces with high “magnification”, i.e. surfaces containing many modes $\zeta=k_{s}/k_{l}\gg 1$, Nayak parameter changes drastically with the Hurst exponent. Thus, the latter has a strong impact on the contact area for surfaces with a rich spectrum, but this impact, we believe, is determined mainly by the corresponding change in Nayak parameter, and not by the Hurst exponent itself. A study of the rough contact for higher Nayak parameters, which possibly can discriminate independent effects of Nayak parameter and Hurst exponent, is still missing.

As a by-product of this study, we suggested a phenomenological second order polynomial estimation of the contact area evolution, which is valid at least in the interval $A^{\prime}\in[1,15]$ %. The coefficients of this equation decay logarithmically with the Nayak parameter, the associated universal constants are also determined. Based on the contact area shape, we gave an explicit formula for a pressure-dependent coefficient of friction. Finally, we showed that the correction of Persson’s model for partial contact in a “phenomenological way” is not an easy task, because of the lack of dependence on Nayak parameter and because of a complex nonlinear evolution of the true contact area, which cannot be accurately described by the error function.

**Acknowledgment

**

This paper was partly presented at 24${}^{\mbox{\tiny th}}$ International Congress of Theoretical and Applied Mechanics (ICTAM) which was hold in Montréal, Canada, and VY is grateful to Association Française de Mécanique (AFM), Comité National Français de Mécanique (CNFM) and to CNRS for the travel support.

**Appendix

**

Appendix A Numerical data

In Tables 1,2, the coefficients for an affine fit of the inverse mean pressure $A^{\prime}/p^{\prime}$ evolution (up to $\approx 15$ % of contact area) are listed for $H=0.4$ and $H=0.8$, respectively. The coefficients are found by the mean least squared error method in the interval $p^{\prime}\in[0.006,0.06]$. These coefficients, depicted in Fig. 9 are used in Eq. (11) and (12) to define the contact area dependence on Nayak parameter.

In Tables 3,4 we list some numerical results obtained for different rough surfaces for $H=0.4$ and $H=0.8$, respectively. Only every third point is shown, and the data is limited to the interval $p^{\prime}<0.06$. The entire set of numerical data is provided in the supplementary material [66].

In Fig. 17, the results from [69] are post-treated using Eq. (1) to obtain accurate results on the contact area evolution, both raw and corrected results are shown. The latter demonstrate a reasonable dependence on the Nayak parameter, the spurious dependence on the cutoff wavenumbers, present in the raw data, is corrected. In Figs. 18,19, some additional results on the contact area evolution are presented for $H=0.4$ and $H=0.8$, respectively. Again, both the raw and corrected data are presented, demonstrating the correcting effect of Eq. (1) suggested in [70].