A Fractal Contact Model for Rough Surfaces considering the Variation of Critical Asperity Levels

Abstract

A contact model for rough surfaces based on the fractal theory is proposed in the present work. Firstly, the deformation of the material is divided into four stages: elastic deformation, the first elastoplastic deformation, the second elastoplastic deformation, and full plastic deformation. And the variation of material hardness is considered when analyzing the contact characteristics of a single asperity within the first and second elastoplastic deformation stages. Secondly, the size distribution function of contact spots at different frequency levels is derived. And the expressions of asperity critical frequency levels are rederived. Lastly, the feasibility and credibility of the proposed model are verified by comparison with other contact models and experimental data. The results show that when the variation of the material hardness is considered, the contact area of a single asperity in the first elastoplastic deformation stage becomes larger, while the contact area of a single asperity in the second elastoplastic deformation stage becomes smaller. Moreover, the critical asperity frequency levels of the rough surface are not constant, but the variables are related to the total real contact area of the rough surface and decrease as the real contact area increases. The proposed model is a modification and improvement of the existing fractal contact models, which can lead to a more accurate relationship between the contact load and the total real contact area of the rough surface.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

The aero-engine external piping system is mainly used for the transmission of fuel, lubricating oil, hydraulic oil, and air and is an important part of the external accessory device [1–3]. Hundreds of pipelines are installed on an aero-engine. Due to their heat resistance and corrosion resistance, tube connectors in the form of metal-to-metal seals (i.e., no dedicated sealing component used) are often used for the connection between pipelines [4]. Tube connectors are the weakest link of the pipeline system’s sealing performance, and tube connector sealing failure has become one of the pipeline system’s main failure modes [4–6]. Once the metal-to-metal seals are out of work, leakage will be formed. Both stability and reliability of the aero-engine will be affected.

The metal-to-metal seal is performed by a direct-metal/metal-tight contact of rough surfaces as shown in Figure 1. Although it is simple in structure, the sealing behavior of a metal-to-metal seal is affected by a variety of factors [7–11], among which the surface topography, which usually has a microstructure given by machining processes, is thought to be one of the most important factors [12–14]. Clearly, engineering materials are known to have rough surfaces, and the full control of surface topography at all scales during manufacturing processes is still out of reach [9]. When two rough surfaces come into contact, the topography of the surface leads to imperfect contact and makes the real contact area only a fraction of the nominal area [4, 15]. Noncontact areas communicate with each other to form leakage channels. Thus, the research on leakage mechanism and sealing performance of metal-to-metal seals involves the resulting geometry of contact between the rough surfaces. Accurately characterizing the contact state and extracting the relationship between the contact load and the real contact area of the two sealing surfaces are necessary for the study of the leakage mechanism and sealing performance of metal-to-metal seals [6, 16–18].

[figure(s) omitted; refer to PDF]

The research on the contact problem of rough surfaces started as early as the Coulomb’s friction proposed by the French engineer Coulomb [19] in 1781. Hertz [20] gave an analytical solution to the contact problem of frictional elastomers in 1882, which opened up the study of modern contact mechanics. Subsequently, researchers have proposed a variety of contact models to describe the contact behavior of rough surfaces. These models mainly include numerical contact models, statistical contact models, and fractal contact models [15, 21–25]. The numerical contact models generally use digital technology (i.e., SEM, AFM, and STM) to obtain the specific parameters of the surface topography and use the finite element method to simulate the contact behavior of two rough surfaces. The advantage of numerical contact models is that they can obtain simulation results that are closer to the actual situation. But they might require, in general, a large and dense grid, and the computational efficiency might be unacceptable [18]. Statistical contact models use statistical parameters to characterize the rough surfaces. The GW contact model proposed by Greenwood [26] in 1966 is a typical representative of the statistical contact models, and subsequent scholars have improved the GW model from different aspects [27, 28]. The advantage of statistical contact models is that their expressions are simple and clear, which greatly simplifies the derivation of the contact equation between rough surfaces and is conducive to rapid contact analysis between rough surfaces. However, the statistical contact models all simplified the characterizations of the surface topography to various degrees, which consequently caused a large error between the calculated and the actual results. Engineering surfaces are found to have fractal characteristics [29, 30]. When fractal parameters instead of statistic parameters are adopted to characterize rough surfaces, there is no need to make too many assumptions about the surface topography, and the fractal characteristics of rough surfaces can be preserved [31–35]. As a result, although there are some skeptical opinions on the fractal approaches [36], contact models based on fractal theory have received more and more attention from researchers. Bhushan and Bhushan [37] proposed one of the first fractal contact models (MB model) based on the fractal theory and Weierstrass-Mandelbrot (WM) function in 1991. Komvopoulos and Komvopoulos [38, 39] believed that the truncated area of an asperity should not be equal to its real contact area. They presented an improved contact model (WK model) in which the deformation mode of an asperity is divided into three stages: elastic deformation, elastoplastic deformation, and full plastic deformation. Considering the elastoplastic deformation of asperities, Komvopoulos and Komvopoulos [40] established a 3D fractal contact model (YK model) for rough surfaces, which is suitable for both isotropic and anisotropic surfaces. The MB model predicts that the asperity first deforms plastically and then elastically during the loading process. This conclusion is contrary to classical contact mechanics and contradicts people’s intuitive feelings. However, it has never been challenged until Etsion and Etsion [41] published their work (ME model) in 2007. In the ME model, the concept of critical asperity frequency level was proposed for the first time. Subsequently, the ME model was developed into a complete contact model by Huang and Huang (MH model) [42]. Yuan et al. [43] proposed a revised contact model (YC model) based on the MB model and the ME model. Based on the YC model, Yuan et al. have successively proposed a loading-unloading contact model for rough surfaces [35, 44], as well as a normal contact stiffness model for joint surfaces [45].

The YC model has solved several deficiencies of the MB model and played a significant role in promoting the research of contact modeling based on the fractal theory. However, the elastic critical frequency level of the YC model is derived on the condition that the asperity height is not greater than itself critical interference. The consequence of this condition in the YC model is that the variation of the critical frequency levels with the total real contact area of the rough surface is not taken into account. In addition, most fractal contact models do not consider the variation of material hardness with deformation. The variation of material hardness has a vital influence on the mechanical properties of the asperities on the rough surface [46], and it hence affects the contact characteristics of the entire rough surface.

Therefore, the main purpose of this paper is to propose a fractal contact model considering the variation of the critical asperity levels as well as the variation of the material hardness to extract the relationship between the contact load and the total real contact area of the rough surface. The present work will lay the foundation for the subsequent analysis of the sealing performance of static metal seals.

2. Theory and Method

2.1. The WM Function

An engineering surface profile often appears random, multiscale, and disordered. The mathematical properties of such a profile are as follows: it is continuous, nondifferentiable, and statistically self-affine [37]. It is found that the WM function satisfies all these properties and is often used to create 2D rough surface profiles. The WM function is given as [37]$$\begin{array}{}\text{(1)}& zx=G^{D-1}{\displaystyle \sum }_{n=n_{\min }}^{\infty }\frac{\cos 2\pi {\gamma }^{n}x}{{\gamma }^{2-Dn}},1<D<2,\gamma >1,\end{array}$$where $zx$ represents the height of the profile along with the $x$ direction. $D$ and $G$ are fractal parameters, representing fractal dimension and characteristic length scale (roughness parameter), respectively. ${\gamma }^n$ determines the frequency spectrum of the surface roughness, and $n$ indicates the frequency level of the asperities. Generally, in order to meet the requirements for high spectral density and phase randomization, $\gamma$ is taken as 1.5. $n_{\min }$ indicates the low cut-off frequency of the profile.

2.2. The Contact of a Single Asperity

The idea of contact modeling based on the fractal theory is to obtain the “contact area-contact load” relationship of the entire rough surface through integration based on the “contact area-contact load” relationship of a single asperity. Therefore, the first step is to obtain the “contact area-contact load” relationship of a single asperity.

2.2.1. The Existing Elastic Microcontacts

According to the WM function, the profile of the fractal asperity at a certain level $n$ before deformation is [37]$$\begin{array}{}\text{(2)}& z_{n}x=G^{D-1}{l}_n^{2-D}\cos \frac{\pi x}{l_n},-\frac{l_n}{2}<x<\frac{l_n}{2},\end{array}$$where $l_n$ represents the length scale of the asperity of level $n$; that is, the base diameter of the asperity of level $n$, $l_n=1/{\gamma }^n$.

The asperity height is [37]$$\begin{array}{}\text{(3)}& {\delta }_n=z_{n}0=G^{D-1}{l}_n^{2-D}.\end{array}$$

The contact between two rough surfaces can be simplified to an equivalent rough surface in contact with a rigid flat surface [42]. Figure 2 shows the relationship between the geometric parameters when a fractal asperity is in contact with a rigid smooth surface. According to Figure 2, the interference ${\omega }_n$ of the asperity with a length scale $l_n$ is$$\begin{array}{}\text{(4)}& {\omega }_n={\delta }_n-d=G^{D-1}{l}_n^{2-D}-d,\end{array}$$where $d$ is the separation distance between the rigid flat surface and the equivalent rough surface. If the surface height $zx,y$ follows the Gaussian distribution, separation distance $d$ can be obtained by the relationship [37]$$\begin{array}{}\text{(5)}& \frac{A_r}{A_a}=\frac{1}{\sqrt{2\pi }}{\int }_{d/\sigma }^{\infty }\exp -x^2/2dx=\frac{1}{2}\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\frac{d}{\sqrt{2}\sigma },\end{array}$$where $A_r$, $A_a$are the total real contact area and the nominal contact area, respectively. erfc($x$) represents the complementary error function.

[figure(s) omitted; refer to PDF]

According to the Hertz theory [37], the contact area and contact load of the asperity within the elastic deformation stage can be obtained as$$\begin{array}{}\text{(6)}& a_n=\pi R_n{\omega }_n,\text{(7)}& f_{ne}=\frac{4}{3}E{R}_n^{1/\mathit{2}}{\omega }_n^{3/\mathit{2}},\end{array}$$where $R_n$ is the radius of curvature at the asperity summit, $R_n=\mathit{1}/{d^2{z}_n/d{x}^2}_{x=0}=l_n^D/{\pi }^2{G}^{D-1}$. $E$ is the equivalent elastic modulus, $\mathit{1}/E=1-{\nu }_1^2/E_1+1-{\nu }_2^2/E_2$, and $E_1$, $E_2$, ${\nu }_1$, and ${\nu }_2$ are elastic modulus and Poisson’s ratios of the two rough surfaces. According to equations (6) and (7), the relationship between contact load and real contact area of a single asperity within the elastic deformation stage is$$\begin{array}{}\text{(8)}& f_{ne}=\frac{4E{\pi }^{1/2}{G}^{D-1}}{3l_n^D}{a}_n^{3/2}.\end{array}$$

2.2.2. The Existing Condition for Initial Yield

The internal stress of the asperity increases with the increase of the contact load or contact interference. An initial yield point will eventually be generated inside the asperity due to excessive stress. The interference corresponding to the initial yield point is termed as the critical interference ${\omega }_{nec}$, and is given as [47]$$\begin{array}{}\text{(9)}& {\omega }_{\mathrm{n}\mathrm{e}\mathrm{c}}={\frac{\pi KH}{2E}}^2{R}_n.\end{array}$$where $K$ is the hardness factor of the softer material and is given by $K=0.454+0.41\nu$. $H$ denotes the hardness of the softer material and is given by $H=2.8Y$, and $Y$ is the yield strength of the softer material. The elastic critical contact area of the asperity $a_{\mathrm{n}\mathrm{e}\mathrm{c}}$ can be obtained as [43]$$\begin{array}{}\text{(10)}& a_{\mathrm{n}\mathrm{e}\mathrm{c}}=\pi R_n{\omega }_{\mathrm{n}\mathrm{e}\mathrm{c}}=\frac{1}{\pi }{\frac{KH}{2E}}^2{\frac{l_n^D}{G^{D-1}}}^2.\end{array}$$

According to equations (7), (9), and (10), we can get the expression of the elastic critical contact load of the asperity $f_{\mathrm{n}\mathrm{e}\mathrm{c}}$, which is$$\begin{array}{}\text{(11)}& f_{\mathrm{n}\mathrm{e}\mathrm{c}}=\frac{2}{3}KH{a}_{\mathrm{n}\mathrm{e}\mathrm{c}}.\end{array}$$

2.2.3. The Revised Elastoplastic Microcontacts

As the load or interference of the asperity increases, the plastic part inside the asperity will gradually expand to the contact surface, and an annular plastic part will be formed on the contact surface, while the rest of the contact surface surrounded by the plastic region remains elastic deformation. This transition stage, with the transition interference ratio of $1<{\omega }_n/{\omega }_{\mathrm{n}\mathrm{e}\mathrm{c}}\le 6$($a_{nec}<a_n\le a_{\mathrm{n}\mathrm{e}\mathrm{p}\mathrm{c}}$), is known as the first elastoplastic deformation stage. The first elastoplastic critical contact area $a_{\mathrm{n}\mathrm{e}\mathrm{p}\mathrm{c}}$ and the first elastoplastic contact load $f_{\mathrm{n}\mathrm{e}\mathrm{p}1}$ of the asperity are derived by Yuan et al. [43] based on the ME model, which are$$\begin{array}{}\text{(12)}& a_{\mathrm{n}\mathrm{e}\mathrm{p}\mathrm{c}}=7.1197a_{\mathrm{n}\mathrm{e}\mathrm{c}},\text{(13)}& f_{\mathrm{n}\mathrm{e}\mathrm{p}1}=\frac{2}{3}KH\times 1.1282a_{\mathrm{n}\mathrm{e}\mathrm{c}}^{-0.2544}{a}_n^{1.2544}.\end{array}$$

In the first elastoplastic deformation stage, the hardness of the material will change with the deformation rather than remaining a constant value [25]. According to equations (12)∼(13), we assume that the material hardness within the first elastoplastic deformation stage $H_{G1}{a}_n$ satisfies the following relationship:$$\begin{array}{}\text{(14)}& H_{G1}{a}_n=c_{11}Y{\frac{a_n}{a_{nec}}}^{c_{12}},a_{nec}<a_n\le a_{nepc},\end{array}$$where $c_{11}$, $c_{12}$are the parameters need to be solved.

Equation (14) should satisfy two boundary conditions:$$\begin{array}{}\text{(15)}& H_{G1}{a}_{nec}=p_{em}{a}_{nec},\text{(16)}& H_{G1}{a}_{nepc}=p_{epm1}{a}_{nepc},\end{array}$$where $p_{em}{a}_n$ is the average contact pressure of the asperity in the elastic deformation stage, which is given by $p_{em}{a}_n=f_{ne}/a_n$. Hence, $p_{em}{a}_{nec}$ can be obtained as $p_{em}{a}_{nec}=2/3KH$ asperity in the first elastoplastic deformation stage and is given by $p_{epm1}{a}_n=f_{nep1}/a_n$. Substituting equation (14) and equation $p_{em}{a}_{nec}=KH$ into equation (15) yields $c_{11}Y=2/3KH$. Thus, the parameter $c_{11}$ can be obtained as$$\begin{array}{}\text{(17)}& c_{11}=\frac{28}{15}K.\end{array}$$

Substituting equations (12)–(14) and (17) into equation (16) yields$$\begin{array}{}\text{(18)}& \frac{2}{3}KH\times {7.1197}^{c_{12}}=\frac{2}{3}KH\times 1.1282a_{nec}^{-0.2544}{7.1197a_{nec}}^{0.2544}.\end{array}$$

And parameter $c_{12}$ can be obtained as$$\begin{array}{}\text{(19)}& c_{12}=\frac{\ln 1.1282\times {7.1197}^{0.2544}}{\ln 7.1197}.\end{array}$$

Therefore, the contact load of a single asperity in the first elastoplastic deformation stage is revised as$$\begin{array}{}\text{(20)}& f_{nep1}^{\prime }=H_{G1}{a}_{nep1}·a_{nep1}=\frac{28}{15}KY{a}_{nec}^{-c_{12}}{a}_n^{c_{12}+1}.\end{array}$$

If the contact load or interference further increases, the plastic part of the asperity gradually expands to envelop the shrinking elastic core. According to Etsion and Etsion [48], the transition interference ratio of this stage is $6<{\omega }_n/{\omega }_{nec}\le 110$($a_{nepc}<a_n\le a_{npc}$). Some researchers named this stage the second elastoplastic deformation stage. The second elastoplastic critical contact area $a_{npc}$ and the second elastoplastic contact load $f_{nep2}$ of the asperity are also given by Yuan et al. [43]:$$\begin{array}{}\text{(21)}& a_{npc}=205.3827a_{nec},\text{(22)}& f_{nep2}=\frac{2}{3}KH\times 1.4988a_{nec}^{-0.1021}{a}_n^{1.1021}.\end{array}$$

The hardness of the material will also change with the deformation in the second elastoplastic deformation stage [25]. Similarly, we assume that the material hardness within the second elastoplastic deformation stage $H_{G2}{a}_n$ satisfies the following relationship:$$\begin{array}{}\text{(23)}& H_{G2}{a}_n=c_{21}Y{\frac{a_n}{a_{nec}}}^{c_{22}},a_{nepc}<a_n\le a_{npc},\end{array}$$where $c_{21}$ and $c_{22}$ are the two parameters need to be solved.

And equation (23) should also satisfy two boundary conditions:$$\begin{array}{}\text{(24)}& H_{G2}{a}_{nepc}=p_{em1}{a}_{nepc},\text{(25)}& H_{G2}{a}_{npc}=p_{epm2}{a}_{npc},\end{array}$$where $p_{epm2}{a}_n$ is the average contact pressure of the asperity in the second elastoplastic deformation stage and is given by $p_{epm2}{a}_n=f_{nep2}/a_n$.

Substituting equations (12), (13), and (23) into equation (24) yields$$\begin{array}{}\text{(26)}& c_{21}Y{7.1197}^{c_{22}}=\frac{2}{3}KH\times 1.1282\times {7.1197}^{0.2544}.\end{array}$$

Substituting equations (21), (22), and (23) into equation (25) yields$$\begin{array}{}\text{(27)}& c_{21}Y{205.3827}^{c_{22}}=\frac{2}{3}KH\times 1.4988\times {205.3827}^{0.1021}.\end{array}$$

According to equations (26) and (27), parameter $c_{22}$ can be expressed as$$\begin{array}{}\text{(28)}& c_{22}=\frac{\ln 1.4988+0.1021\ln 205.3827-\ln 1.1282-0.2544\ln 7.1197}{\ln 205.3827-\ln 7.1197}.\end{array}$$

And parameter $c_{21}$ can be obtained as$$\begin{array}{}\text{(29)}& c_{21}=K\times \frac{28}{15}\times 1.1282\times {7.1197}^{0.2544-c_{22}}.\end{array}$$

Therefore, the contact load of a single asperity in the second elastoplastic deformation stage is revised as$$\begin{array}{}\text{(30)}& f_{nep2}^{\prime }=H_{G2}{a}_{nep2}·a_{nep2}=K\times \frac{28}{15}\times 1.1282\times {7.1197}^{0.2544-c_{22}}Y{a}_{nec}^{-c_{22}}{a}_n^{1+c_{22}}.\end{array}$$

2.2.4. The Existing Plastic Microcontacts

When ${\omega }_n/{\omega }_{nec}>110$($a_n>a_{npc}$), the asperity undergoes full plastic deformation. The contact area and contact load of the asperity in this deformation stage are [43]$$\begin{array}{c}\text{(31)}& a_n=2\pi R_n{\omega }_n,f_{np}=H{a}_n.\end{array}$$

For the sake of clarity, Figure 3 shows the deformation law of a single asperity with the fractal asperity frequency level $n$, as well as the relationship between each critical contact area.

[figure(s) omitted; refer to PDF]

2.3. The Revised Size Distribution Function of Contact Spots

In the MB model, the size distribution function of contact spots is defined as [37]$$\begin{array}{}\text{(32)}& na=\frac{D}{2}{a}_l^{\frac{D}{2}}{a}^{-\frac{D+2}{2}},0<a<a_l,1<D<2,\end{array}$$where $a_l$ denotes the largest contact area of asperity. And the total real contact area is [37]$$\begin{array}{}\text{(33)}& A_r={\int }_0^{a_l}naada=\frac{D}{2-D}{a}_l.\end{array}$$

According to equation (2), the period of asperity at frequency level $n$ is $T_n=2\pi /2\pi {\gamma }^n$, and the period of asperity at frequency level $n+1$ is $T_{n+1}=2\pi /2\pi {\gamma }^{n+1}$. Therefore, $T_{n+1}/T_n=1/\gamma$. Hence, we can obtain the relationship between the size distribution function of the asperity at frequency level $n$ and the size distribution function of the asperity at frequency level $n+1$, which is $n_{n+1}a=\gamma n_{n}a$. Letting the size distribution function of the asperity at the initial frequency level $n_{\min }$ be $n_{n_{\min }}a=Mna$, then the contact spot size distribution function at each frequency level can be obtained as$$\begin{array}{}\text{(34)}& n_{n}a=M{\gamma }^{n-n_{\min }}na,n_{\min }\le n\le n_{\max },\end{array}$$where $M$ is a parameter that needs to be solved. According to equations (33) and (34), the total real contact area of asperities at each frequency level is$$\begin{array}{}\text{(35)}& A_{nr}={\int }_0^{a_{nl}}{n}_{n}aada=\frac{MD}{2-D}{\gamma }^{n-n_{\min }}{a}_{nl},\end{array}$$where $a_{nl}$ is the largest contact area of contact spots at frequency level $n$. According to the geometry relationship shown in Figure 2, we assume that the largest contact spot of frequency level $n_{\min }$ is the largest contact spot of the entire rough surface, i.e., $a_{n_{\min }l}=a_l$. And we assume that the relationship between the largest contact spots of two adjacent frequency levels is $a_{n+1l}/a_{nl}=1/{\gamma }^2$. Then, the largest contact area of the asperity at any frequency level would be $a_{nl}=1/{\gamma }^{2n-n_{\min }}{a}_l$. Thus equation (35) can be reexpressed as$$\begin{array}{}\text{(36)}& A_{nr}=\frac{M}{{\gamma }^{n-n_{\min }}}\frac{D}{2-D}{a}_l=\frac{M}{{\gamma }^{n-n_{\min }}}{A}_r.\end{array}$$

And the relationship between $a_{nl}$ and $A_r$ can be obtained as$$\begin{array}{}\text{(37)}& a_{nl}=\frac{2-D}{D}\frac{1}{{{\gamma }^{n-n_{\min }}}^2}{A}_r.\end{array}$$

The parameter $M$ can be determined according to the relationship $A_{n_{\min }r}+A_{n_{\min +1}r}+\cdots +A_{n_{\max }r}=A_r$. Substituting equation (36) into it yields $M=1-1/\gamma /1-{1/\gamma }^{n_{\max }-n_{\min }+1}$.

2.4. The Revised Critical Asperity Frequency Levels

For a certain frequency level $n$, the relationship between the interference ${\omega }_n$ and the asperity height ${\delta }_n$ is ${\omega }_n\le {\delta }_n$ [46]. If the height of an asperity is less than its elastic critical interference, inelastic deformation occurs. Thus, the condition of elastic deformation for an asperity given by the YC model is ${\delta }_n\le {\omega }_{nec}$ [43]. However, according to Etsion and Etsion [41], the condition for elastic deformation of an asperity should be ${\omega }_n\le {\omega }_{nec}$, which indicates that the asperity occurs elastic deformation when the interference is not greater than its critical interference. Etsion and Etison [41] further derived the inequality used to solve the elastic critical length scale $l_{ec}$(or the elastic critical asperity frequency level $n_{ec}$), which is $1-\cos \pi l_{nt}/2l_n\le {KH/2E\cdot {l_n/G}^{D-1}}^2$. Unfortunately, we cannot get the analytical solution of $l_{ec}$ according to this inequality, nor can we get an analytical solution of $n_{ec}$ according to the relationship $l_n=1/{\gamma }^n$.

Now, we reanalyze the conditions for judging the elastic deformation of the asperity. At a mean surface separation distance $d$, if the largest interference of the asperity at frequency level $n$ is not greater than its critical interference, all of the contact asperities at frequency level $n$ deform elastically. Thus, the condition we give to judge the elastic deformation of the asperity is$$\begin{array}{}\text{(38)}& {\omega }_{nl}\le {\omega }_{nec},\end{array}$$where ${\omega }_{nl}$ is the largest interference of the contact spot at frequency level $n$. According to equation (6), equation (39) can be re-written as$$\begin{array}{}\text{(39)}& a_{nl}\le a_{nec}.\end{array}$$

Substituting equations (10) and (37) into the above inequality, the elastic critical frequency level can be obtained from the following equation:$$\begin{array}{}\text{(40)}& \frac{2-D}{D}\frac{1}{{{\gamma }^{n_{ec}-n_{\min }}}^2}{A}_r=\frac{1}{\pi }{\frac{KH}{2E}}^2\frac{1}{{\gamma }^{2n_{ec}D}}\frac{1}{G^{2D-1}}.\end{array}$$

According to the above inequality, the elastic critical frequency level $n_{ec}$ is not only related to the material parameters $H,E$ and the topography parameters $D,G$ but also related to the total real contact area of the two rough surfaces $A_r$. Similarly, the first elastoplastic critical frequency level $n_{epc}$ and the second elastoplastic critical frequency level $n_{pc}$ can be obtained from the following two equations:$$\begin{array}{c}\text{(41)}& \frac{2-D}{D}\frac{1}{{{\gamma }^{n_{epc}-n_{\min }}}^2}{A}_r=7.1197\frac{1}{\pi }{\frac{KH}{2E}}^2\frac{1}{{\gamma }^{2n_{epc}D}}\frac{1}{G^{2D-1}},\frac{2-D}{D}\frac{1}{{{\gamma }^{n_{pc}-n_{\min }}}^2}{A}_r=205.3827\frac{1}{\pi }{\frac{KH}{2E}}^2\frac{1}{{\gamma }^{2n_{pc}D}}\frac{1}{G^{2D-1}}.\end{array}$$

According to equations (41)∼(43), the equations used to calculate the critical frequency levels $n_{ec},n_{epc},n_{pc}$ are not only related to the material parameters $H,E$ and topography parameters $D,G$, but also related to the total real contact area of the rough surface $A_r$. As the given parameter $A_r$ changes, so do the calculation results of the critical frequency levels $n_{ec},n_{epc},n_{pc}$. That is, the critical frequency levels $n_{ec},n_{epc},n_{pc}$ variate with the change of the $A_r$ value. However, in the equations for calculating the critical frequency levels given by the YC model, the critical frequency levels are only related to the material parameters and topography parameters, and the influence of the total real contact area $A_r$ on the critical frequency levels is not considered.

2.5. The Revised Real Contact Area and Contact Load of the Rough Surface

When the asperity frequency level ranges from $n_{\min }$ to $n_{ec}$, i.e., $n_{\min }\le n\le n_{ec}$, only elastic deformation takes place. And the real contact area $A_{r1}$ and contact load $F_{r1}$ are obtained as$$\begin{array}{c}\text{(42)}& A_{r1}=A_{re}={\displaystyle \sum }_{n=n_{\min }}^{n_{ec}}{\int }_0^{a_{nl}}M{\gamma }^{n-n_{\min }}naada,F_{r1}=F_{re}={\displaystyle \sum }_{n=n_{\min }}^{n_{ec}}{\int }_0^{a_{nl}}{f}_{ne}M{\gamma }^{n-n_{\min }}nada.\end{array}$$

When the asperity frequency level belongs to $n_{ec}<n\le n_{epc}$, elastic or the first elastoplastic deformation can take place. The real contact area $A_{r2}$ is obtained as$$\begin{array}{}\text{(43)}& A_{r2}=A_{re}+A_{rep1}={\displaystyle \sum }_{n=n_{ec}+1}^{n_{epc}}{\int }_0^{a_{nec}}M{\gamma }^{n-n_{\min }}naada+{\displaystyle \sum }_{n=n_{ec}+1}^{n_{epc}}{\int }_{a_{nec}}^{a_{nl}}M{\gamma }^{n-n_{\min }}naada.\end{array}$$

And the contact load $F_{r2}$ is obtained as$$\begin{array}{}\text{(44)}& F_{r2}=F_{re}+F_{rep1}={\displaystyle \sum }_{n=n_{ec}+1}^{n_{epc}}{\int }_0^{a_{nec}}{f}_{ne}M{\gamma }^{n-n_{\min }}nada+{\displaystyle \sum }_{n=n_{ec}+1}^{n_{epc}}{\int }_{a_{nec}}^{a_{nl}}{f}_{nep1}^{\prime }M{\gamma }^{n-n_{\min }}nada.\end{array}$$

When the asperity frequency level is $n_{epc}<n\le n_{pc}$, elastic deformation, the first elastoplastic deformation, or the second elastoplastic deformation will occur. The real contact area $A_{r3}$ is evaluated as$$\begin{array}{}\text{(45)}& \begin{array}{c}{A}_{r3}=A_{re}+A_{rep1}+A_{rep2}={\displaystyle \sum }_{n=n_{epc}+1}^{n_{pc}}{\int }_0^{a_{nec}}M{\gamma }^{n-n_{\min }}naada+{\displaystyle \sum }_{n=n_{epc}+1}^{n_{pc}}{\int }_{a_{nec}}^{a_{nepc}}M{\gamma }^{n-n_{\min }}naada\\ +{\displaystyle \sum }_{n=n_{epc}+1}^{n_{pc}}{\int }_{a_{nepc}}^{a_{nl}}M{\gamma }^{n-n_{\min }}naada,\end{array}.\end{array}$$

And the contact load $F_{r3}$ is evaluated as$$\begin{array}{}\text{(46)}& \begin{array}{c}{F}_{r3}=F_{re}+F_{rep1}+F_{rep2}={\displaystyle \sum }_{n=n_{epc}+1}^{n_{pc}}{\int }_0^{a_{nec}}{f}_{ne}M{\gamma }^{n-n_{\min }}nada+{\displaystyle \sum }_{n=n_{epc}+1}^{n_{pc}}{\int }_{a_{nec}}^{a_{nepc}}{f}_{nep1}^{\prime }M{\gamma }^{n-n_{\min }}nada\\ +{\displaystyle \sum }_{n=n_{epc}+1}^{n_{pc}}{\int }_{a_{nepc}}^{a_{nl}}{f}_{nep2}^{\prime }M{\gamma }^{n-n_{\min }}nada,\end{array}.\end{array}$$

When the asperity frequency level belongs to $n_{pc}<n\le n_{\max }$, all of the four types of deformation can take place. The real contact area $A_{r4}$ is calculated as$$\begin{array}{}\text{(47)}& \begin{array}{c}{A}_{r4}=A_{re}+A_{rep1}+A_{rep2}+A_{rp}={\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\int }_0^{a_{nec}}M{\gamma }^{n-n_{\min }}naada+{\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\int }_{a_{nec}}^{a_{nepc}}M{\gamma }^{n-n_{\min }}naada,\\ +{\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\int }_{a_{nepc}}^{a_{npc}}M{\gamma }^{n-n_{\min }}naada+{\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\int }_{a_{npc}}^{a_{nl}}M{\gamma }^{n-n_{\min }}naada\end{array}.\end{array}$$

And the contact load $F_{r4}$ is estimated as$$\begin{array}{c}\text{(48)}& F_{r4}=F_{re}+F_{rep1}+F_{rep2}+F_{rp}={\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\int }_0^{a_{nec}}{f}_{ne}M{\gamma }^{n-n_{\min }}nada+{\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\int }_{a_{nec}}^{a_{nepc}}{f}_{nep1}^{\prime }M{\gamma }^{n-n_{\min }}nada+{\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\int }_{a_{nepc}}^{a_{npc}}{f}_{nep2}^{\prime }M{\gamma }^{n-n_{\min }}nada+{\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\int }_{a_{npc}}^{a_{nl}}{f}_{np}M{\gamma }^{n-n_{\min }}nada.\end{array}$$

The detailed calculation results of the real contact area and contact load for the rough surface are provided in Appendix A. The total real contact area and contact load of all asperity levels are estimated by$$\begin{array}{c}\text{(49)}& A_r=A_{r1}+A_{r2}+A_{r3}+A_{r4}.F_r=F_{r1}+F_{r2}+F_{r3}+F_{r4}.\end{array}$$

The nondimensional forms of $A_r$ and $F_r$ are defined as [43]$$\begin{array}{c}\text{(50)}& A_r^{\ast }=\frac{A_r}{A_a},F_r^{\ast }=\frac{F_r}{A_{a}E},\end{array}$$where the nominal contact area $A_a$ is given by $A_a=L^2$, $L=1/{\gamma }^{n_{\min }}$. For the sake of clarity, the flowchart of the solution procedure for extracting the relationship between the total real contact area and contact load is shown in Figure 4.

[figure(s) omitted; refer to PDF]

3. Results and Discussion

In order to obtain the contact parameters, the equivalent rough surface parameters [40, 46, 49] used in the present work are shown in Table 1. Let the topography parameters be $G=2.5\times {10}^{-9}$ m, $D=1.5$. Let the sampling length and the resolution be $L=0.5\times {10}^{-4}$ m ($L=1/{\gamma }^{n_{\min }}$), $L_s=1.5\times {10}^{-9}$m ($L_s=1/{\gamma }^{n_{\max }}$), respectively. The nominal contact area, the range of asperity frequency levels, and the critical asperity frequency levels, hence, can be calculated as $A_a=L^2=0.25\times {10}^{-8}{\mathrm{m}}^2$, $n_{\min }=24$, $n_{\max }=50$, $n_{ec}=32$, $n_{epc}=37$, $n_{pc}=45$. It should be noted that parameters $D$ and $G$ refer to the fractal parameters of the equivalent rough surface, which can be calculated from the topography parameters of the two contact surfaces. The calculation procedure is given in Appendix B.

Table 1

Equivalent rough surface parameters.

Parameters	Values
Equivalent elastic modulus $E$	$7.2\times {10}^{10}$ N/㎡
Poisson’s ratio $\nu$	0.17
Initial hardness $H$	$5.5\times {10}^9$ N/㎡
Profile scale parameter $G$	$1\times {10}^{-9}\sim 1\times {10}^{-6}$ m
Fractal dimension $D$	$1<D<2$

Figure 5 shows the comparison of the proposed model and the MB model on the critical contact areas and the largest contact area at different asperity levels. The value of nondimensional total real contact area $A_r^{\ast }$($A_r^{\ast }=A_r/A_a$) is 0.4. It is can be seen from Figure 5 that the value of critical contact areas ($a_{nec}$, $a_{nepc}$, $a_{npc}$) and the largest contact area $a_{nl}$ of the present work decrease with the increase of frequency level $n$. And when $n$ is taken as a specific value, the relationship between $a_{nec}$, $a_{nepc}$ and $a_{npc}$ of the present work is $a_{nec}<a_{nepc}<a_{nnc}$. Moreover, in the present work, the largest contact area is smaller than the elastic critical contact area when the value of $n$ belongs to $0,32$. The largest contact area is larger than the elastic critical contact area and smaller than the first elastoplastic critical contact area if the asperity level $n$ belongs to $33,37$. The largest contact area is larger than the first elastoplastic critical contact area and smaller than the second elastoplastic critical contact area if the asperity level $n$ belongs to $38,45$. The largest contact area is larger than the second elastoplastic critical contact area if the asperity level $n$ belongs to $46,50$. However, the largest contact area $a_{nl}$ and the elastic critical contact area $a_{nec}$ of the MB model do not change with the increase of $n$. And the largest contact area of the MB model is always larger than its elastic critical contact area. In addition, Figure 5 also shows the relationship between the asperity frequency level and the critical contact areas in the YC model. And the relationships between $n$ and critical contact areas in the YC model are the same as those in the proposed model. The analytic expression of $a_{nl}$ is not given by the YC model, so we cannot give a comparison of $a_{nl}$ between the proposed model and the YC model in Figure 5.

[figure(s) omitted; refer to PDF]

The relationships between the fractal dimension $D$ and the critical asperity frequency levels of the proposed model are shown in Figure 6. Figure 6 also shows the changing trend of critical asperity frequency levels with fractal dimension $D$ in the YC model. The characteristic length scale is $G=2.5\times {10}^{-9}$m, and the nondimensional total real contact area $A_r^{\ast }$ is 0.4. It should be noted that, to facilitate comparison, the critical frequency level is set to be 0 if it is not within the range of $n_{\min }$ and $n_{\max }$. Figure 6 indicates that the elastic critical level, elastoplastic critical level, and plastic critical level all increase with the increase of fractal dimension in both the proposed model and the YC model. Moreover, the relationship between $n_{ec}$, $n_{epc}$ and $n_{pc}$ in the proposed model is the same as that in the YC model, which is $n_{ec}<n_{epc}<n_{pc}$. In addition, when $D$ is taken as a specific value, the critical levels calculated by the proposed model are larger than the critical levels of the same type calculated by the YC model. The reason for this difference is that the method of solving the elastic critical level in the proposed model is different from that in the YC model. The conditions for the proposed model and the YC model to solve the elastic critical level are ${\omega }_{nl}\le {\omega }_{nec}$ and ${\delta }_n\le {\omega }_{nec}$, respectively. And the asperity interference will not be greater than its height, i.e., ${\omega }_{nl}\le {\delta }_n$, so the value of $n_{ec}$ calculated by the YC model is smaller than that calculated by the proposed model. It is also due to the same reason that the value of $n_{epc}$ and $n_{pc}$ calculated by the YC model are smaller than those calculated by the proposed model.

[figure(s) omitted; refer to PDF]

The relationships between nondimensional total real contact area $A_r^{\ast }$ and the critical asperity frequency levels of the proposed model and the YC model are presented in Figure 7. Figure 7 indicates that the critical frequency levels of the proposed model decrease with the increase of $A_r^{\ast }$ value. However, as the $A_r^{\ast }$ value increases, the critical frequency levels of the YC model remain unchanged. According to equations (41)∼(43), the asperity critical frequency levels of the proposed model are related to the total real contact area of the rough surface, and the values of $n$ and the total real contact area $A_r$ are negatively correlated. Therefore, the critical frequency levels in the proposed model decrease with the increase of $A_r^{\ast }$. The influence of the total real contact area on the asperity critical frequency levels is not considered in the YC model, and the change of $A_r^{\ast }$ will consequently not affect the calculation results of the critical frequency levels.

[figure(s) omitted; refer to PDF]

Figure 8 shows the “contact load–contact area” relationships of a single asperity both with and without considering the variation of material hardness within the first elastoplastic deformation stage. The frequency level $n$ is 33, and the value of $A_r^{\ast }$ is 0.4. According to Figure 8, whether the variation of the material hardness is considered or not, the “contact load–contact area” relationship of a single asperity shows the same trend; that is, the contact area increases with the increase of contact load. In addition, the contact area with the consideration of material hardness variation is larger than that without the consideration of material hardness variation with the same contact load. And the difference between them gradually shrinks as the contact load increases.

[figure(s) omitted; refer to PDF]

According to the previous analysis, the material hardness within the first elastoplastic deformation stage can be expressed as a function of the contact area of a single asperity. Figure 9(a) presents the relationships between the material hardness and the contact area of a single asperity for different values of $n$. It is can be seen from Figure 9(a) that the material hardness in the first elastoplastic deformation stage increases with the increase of contact area. And when the contact area is taken as a fixed value, the larger the value of $n$, the larger the material hardness. Figure 9(b) shows that when the contact area is taken as a fixed value, the smaller the value of $D$, the larger the material hardness. Figure 9(c) demonstrates that when the contact area is taken as a fixed value, the larger the value of $G$, the larger the material hardness.

[figure(s) omitted; refer to PDF]

The “contact load–contact area” relationships of a single asperity with and without the consideration of material hardness variation within the second elastoplastic deformation stage are shown in Figure 10. The asperity frequency level is 38, and the value of $A_r^{\ast }$ is 0.4. In the second elastoplastic deformation stage, regardless of whether the variation of material hardness is considered or not, the contact area increases with the increase of contact load. And the real contact area without considering the variation of material hardness is larger than the contact area considering the variation of material hardness for the same contact load. Moreover, as the contact load increases, the difference between the contact areas tends to increase.

[figure(s) omitted; refer to PDF]

Figure 11 shows the relationship between the material hardness and the contact area of a single asperity for different values of asperity level $n$, fractal dimension $D$, and roughness parameter $G$ in the second elastoplastic deformation stage. Similar to the first elastoplastic deformation stage, the material hardness increases as the contact area of a single asperity increases. Moreover, the influences of different values of asperity level, fractal dimension, and roughness parameter on the relationship between the contact area and material hardness are the same as those in the first elastoplastic deformation stage.

[figure(s) omitted; refer to PDF]

Figure 12(a) indicates the comparison of the proposed model and the MB model on the relationship between nondimensional total contact load $F_r^{\ast }$ and nondimensional total real contact area $A_r^{\ast }$ for different roughness parameters. The comparison of the proposed model and the MB model on the relationship between nondimensional total real contact area $A_r^{\ast }$ and the fraction of elastic contact area for different roughness parameters is shown in Figure 12(b). In order to obtain reasonable comparison results, some surface topography parameters and material parameters are given, such as $D=1.5$, $Y/E=0.01$. It is can be seen from Figure 12(a) that, in the MB model, as $G/\sqrt{A_a}$ increases, the $F_r^{\ast }$ value required to produce a specific $A_r^{\ast }$ value increases. The same changing trend is also found in the presented model. In addition, the $F_r^{\ast }$ value of the proposed model required to produce a particular $A_r^{\ast }$ value is less than that of the MB model. Figure 12(b) indicates that, in the MB model, the fraction of elastic contact area increases as $A_r^{\ast }$ increases. This is the contradiction between the MB model and classical contact mechanics. Classical contact mechanics believes that as the real contact area increases, the contact load increases, which means that more asperities will undergo inelastic deformation. As a result, the fraction of elastic contact area should decrease. In addition, it is can be obtained from Figure 12(b) that the fraction of elastic contact area in the MB model decreases with the increase of $G/\sqrt{A_a}$. In the proposed model, the fraction of elastic contact area is approximately equal to 1, and it hardly changes with the change of $A_r^{\ast }$ value. It does not change with the change of $G/\sqrt{A_a}$ value either.

[figure(s) omitted; refer to PDF]

Figure 13(a) shows the comparison of the present model and the MB model on the relationship between nondimensional total contact load $F_r^{\ast }$ and nondimensional total real contact area $A_r^{\ast }$ for different values of $D$. Figure 13(a) indicates that, in the MB model, the $F_r^{\ast }$ value required for a particular $A_r^{\ast }$ value decreases with the increase of $D$ when $D$ belongs to $1.1,1.5$. When $D$ is between 1.5 and 1.9, as the value of $D$ increases, the $F_r^{\ast }$ value required to produce a specific $A_r^{\ast }$ value increases. However, in the proposed model, the $F_r^{\ast }$ value needed for a specific $A_r^{\ast }$ value decreases as the value of $D$ increases within the global value range of $D$. Furthermore, it is also can be seen from Figure 13(a) that the $F_r^{\ast }$ value corresponding to the same $A_r^{\ast }$ value in the proposed model is less than that in the MB model when the fractal dimension takes the same value. Figure 13(b) indicates that, in the MB model, the fraction of elastic contact area is 0 when $D=1.1$, which means that no asperity deforms elastically. And when $D=1.5$ and $D=1.9$, the fraction of elastic contact area increases with the increase of $A_r^{\ast }$ value. As mentioned above, this phenomenon is contrary to classical contact mechanics, while in the presented model, when $D=1.1$, the fraction of elastic contact area decreases as $A_r^{\ast }$ value increases. This trend of change is consistent with classical contact mechanics. When $D$ takes the value of 1.5 and 1.9, the fraction of elastic contact area in the presented model hardly changes with the increase of $A_r^{\ast }$ value and is approximately equal to 1. In addition, the fraction of elastic contact area in the proposed model is greater than that in the MB model when the fractal dimension takes the same value. It should be noted that the YC model does not give the calculation formula of the largest contact area of the asperity $a_{nl}$ at frequency level $n$, which makes it impossible for us to program the YC model. Thus, we do not give a comparison between the proposed model and the YC model regarding the $A_r/A_a\sim F_r/A_{a}E$ relationship and the $A_{re}/A_r\sim A_r/A_a$ relationship. Nevertheless, the idea of contact modeling on rough surfaces in reference [43] still has a vital reference significance.

[figure(s) omitted; refer to PDF]

Figure 14 presents the comparison results of the proposed model and other contact models with experimental data on the relationship between the nondimensional total contact load and the nondimensional total real contact area. In order to ensure the reliability of the comparison results, some parameters used in the present work are the same as those adopted by the MB model and the YC model, i.e., $D=1.49$, $G/\sqrt{A_a}=1\times {10}^{-10}$, $Y/E=0.05$. It should be noted that the critical contact areas of the proposed model are scale-dependent. Therefore, the values of the critical contact areas are related to the sampling length and the resolution of the measuring instrument. Different sampling lengths and instrument resolutions will result in different critical contact areas, which will lead to the change in the relationship between nondimensional total contact load $F_r^{\ast }$ and nondimensional total real contact area $A_r^{\ast }$. Figure 14(a) shows the comparison results for $n_{\min }=67$, $n_{\max }=90$. At lower loads, the MB model is in good agreement with the experimental data. Under the same load within the lower loads, the results of the proposed model and the YC model are larger than the experimental data, and the proposed model is closer to experimental data than the YC model. At higher loads, the MB model deviates seriously from the experimental results. The explanation given by Majumdar and Bushan is that the interaction of asperities is not considered. However, both the proposed model and the YC model are close to the experimental values at higher loads and are better than the results of the GW model. Figure 14(b) indicates that the contact areas of the proposed model and the YC model are larger than the experimental value for the same load at lower loads, and the YC model is closer to the experimental results. At higher loads, the results of the proposed model and the YC model are less than the experimental values, and the proposed model is in better agreement with the experimental results. In addition, it can be seen from Figures 14(a) and 14(b) that, in the entire loading range, the agreement between the proposed model and the experimental data is better than that of other models.

[figure(s) omitted; refer to PDF]

4. Conclusions

An improved fractal contact model considering the variation of the critical asperity levels as well as the variation of the material hardness is proposed in the present work. The main conclusions are as follows:

(1) The real contact area of a single asperity obtained by considering the variation of material hardness is greater than that without considering the variation of material hardness within the first elastoplastic deformation stage, while in the second elastoplastic deformation stage, the real contact area of a single asperity considering the variation of material hardness is less than that without considering the variation of material hardness.

(2) The size distribution functions of the contact spots at different frequency levels are derived. The expressions of asperity critical frequency levels are rederived. The results show that the critical asperity levels are not constant values, but variable values related to the total real contact area of the rough surface and decrease with the increase of the total real contact area.

(3) The proposed model is a modification and improvement of the existing fractal contact models, which can lead to a more accurate relationship between the contact load and the total real contact area of the rough surface.

The proposed model is helpful for the analysis of the sealing performance of static metal seals under different contact pressures in our subsequent studies.

Acknowledgments

This work was supported by the Fundamental Research Funds for the Central Universities of the Civil Aviation University of China (Grant number: 3122019093) and the Scientific Research Project of Tianjin Municipal Education Commission (Grant number: 2020KJ017).

Glossary

Subscripts

$a$:Apparent or nominal

$c$:Critical

$e$:Elastic

$l$:Largest

$m$:Mean or average

$p$:Plastic

$r$:Real

$s$:Smallest.

Appendix

A. Detailed Calculation Results of Real Contact Area and Contact Load for the Rough Surface

A.1. When the Asperity Frequency Level Belongs to $n_{\min }\le n\le n_{ec}$

The real elastic contact area $A_{re}$ and elastic contact load $F_{re}$ are evaluated as$$\begin{array}{}\text{(A.1)}& A_{re}={\displaystyle \sum }_{n=n_{\min }}^{n_{ec}}{\int }_0^{a_{nl}}M{\gamma }^{n-n_{\min }}naada=\frac{MD}{2-D}{\displaystyle \sum }_{n=n_{\min }}^{n_{ec}}{\gamma }^{n-n_{\min }}{a}_{nl},\text{(A.2)}& F_{re}={\displaystyle \sum }_{n=n_{\min }}^{n_{ec}}{\int }_0^{a_{nl}}{f}_{ne}M{\gamma }^{n-n_{\min }}nada=\frac{MD}{3-D}{\displaystyle \sum }_{n=n_{\min }}^{n_{ec}}\frac{4E{\pi }^{1/2}{G}^{D-1}}{3{1/{\gamma }^n}^D}{\gamma }^{n-n_{\min }}{a}_{nl}^{\mathit{3}/\mathit{2}}.\end{array}$$

A.2. When the Asperity Frequency Frequency Level Belongs to $n_{ec}<n\le n_{epc}$

The real elastic contact area $A_{re}$ and the real first elastoplastic contact area $A_{rep1}$ are evaluated as$$\begin{array}{}\text{(A.3)}& A_{re}={\displaystyle \sum }_{n=n_{ec}+1}^{n_{epc}}{\int }_0^{a_{nec}}M{\gamma }^{n-n_{\min }}naada=\frac{MD}{2-D}{\displaystyle \sum }_{n=n_{ec}+1}^{n_{epc}}{\gamma }^{n-n_{\min }}{a}_{nec}^{2-D/2}{a}_{nl}^{D/2},\text{(A.4)}& A_{rep1}={\displaystyle \sum }_{n=n_{ec}+1}^{n_{epc}}{\int }_{a_{nec}}^{a_{nl}}M{\gamma }^{n-n_{\min }}naada=\frac{MD}{2-D}{\displaystyle \sum }_{n=n_{ec}+1}^{n_{epc}}{\gamma }^{n-n_{\min }}{a}_{nl}^{D/2}{a}_{nl}^{2-D/2}-a_{nec}^{2-D/2}.\end{array}$$

And the elastic contact load $F_{re}$ and the first elastoplastic contact load $F_{rep1}$ are evaluated as$$\begin{array}{c}\text{(A.5)}& F_{re}={\displaystyle \sum }_{n=n_{ec}+1}^{n_{epc}}{\int }_0^{a_{nec}}{f}_{ne}M{\gamma }^{n-n_{\min }}nada=\frac{MD}{3-D}{\displaystyle \sum }_{n=n_{ec}+1}^{n_{epc}}\frac{4E{\pi }^{1/2}{G}^{D-1}}{3{1/{\gamma }^n}^D}{\gamma }^{n-n_{\min }}{a}_{nec}^{3-D/2}{a}_{nl}^{D/2},\text{(A.6)}& F_{rep1}={\displaystyle \sum }_{n=n_{ec}+1}^{n_{epc}}{\int }_{a_{nec}}^{a_{nl}}{f}_{nep1}^{\prime }M{\gamma }^{n-n_{\min }}nada=\frac{1.8667KYMD}{2c_{12}-D+2}{\displaystyle \sum }_{n=n_{ec}+1}^{n_{epc}}{\gamma }^{n-n_{\min }}{a}_{nec}^{-c_{12}}{a}_{nl}^{c_{12}+1}-a_{nec}^{2-D/2}{a}_{nl}^{D/2}.\end{array}$$

A.3. When the Asperity Frequency Level Is $n_{epc}<n\le n_{pc}$

The real elastic contact area $A_{re}$, the real first elastoplastic contact area $A_{rep1}$, and the real second elastoplastic contact area $A_{rep2}$ are evaluated as$$\begin{array}{}\text{(A.7)}& A_{re}={\displaystyle \sum }_{n=n_{epc}+1}^{n_{pc}}{\int }_0^{a_{nec}}M{\gamma }^{n-n_{\min }}naada=\frac{MD}{2-D}{\displaystyle \sum }_{n=n_{epc}+1}^{n_{pc}}{\gamma }^{n-n_{\min }}{a}_{nec}^{2-D/2}{a}_{nl}^{D/2},\text{(A.8)}& A_{rep1}={\displaystyle \sum }_{n=n_{epc}+1}^{n_{pc}}{\int }_{a_{nec}}^{a_{nepc}}M{\gamma }^{n-n_{\min }}naada=\frac{MD}{2-D}{\displaystyle \sum }_{n=n_{epc}+1}^{n_{pc}}{\gamma }^{n-n_{\min }}{a}_{nl}^{D/2}{a}_{nepc}^{2-D/2}-a_{nec}^{2-D/2},\text{(A.9)}& A_{rep2}={\displaystyle \sum }_{n=n_{epc}+1}^{n_{pc}}{\int }_{a_{nepc}}^{a_{nl}}M{\gamma }^{n-n_{\min }}naada=\frac{MD}{2-D}{\displaystyle \sum }_{n=n_{epc}+1}^{n_{pc}}{\gamma }^{n-n_{\min }}{a}_{nl}^{D/2}{a}_{nl}^{2-D/2}-a_{nepc}^{2-D/2}.\end{array}$$

And the elastic contact load $F_{re}$, the first elastoplastic contact load $F_{rep1}$, and the second elastoplastic contact load $F_{rep2}$ are evaluated as$$\begin{array}{c}\text{(A.10)}& F_{re}={\displaystyle \sum }_{n=n_{epc}+1}^{n_{pc}}{\int }_0^{a_{nec}}{f}_{ne}M{\gamma }^{n-n_{\min }}nada=\frac{MD}{3-D}{\displaystyle \sum }_{n=n_{epc}+1}^{n_{pc}}\frac{4E{\pi }^{1/2}{G}^{D-1}}{3{1/{\gamma }^n}^D}{\gamma }^{n-n_{\min }}{a}_{nec}^{3-D/2}{a}_{nl}^{D/2},\text{(A.11)}& \begin{array}{c}{F}_{rep1}={\displaystyle \sum }_{n=n_{epc}+1}^{n_{pc}}{\int }_{a_{nec}}^{a_{nepc}}{f}_{nep1}^{\prime }M{\gamma }^{n-n_{\min }}nada\\ =\frac{1.8667KYMD}{2+2c_{12}-D}{\displaystyle \sum }_{n=n_{epc}+1}^{n_{pc}}{\gamma }^{n-n_{\min }}{a}_{nec}^{-c_{12}}{a}_{nepc}^{1+c_{12}-\frac{D}{2}}-a_{nec}^{1+c_{12}-\frac{D}{2}}{a}_{n\mathrm{l}}^{D/2},\end{array}\text{(A.12)}& F_{rep2}={\displaystyle \sum }_{n=n_{epc}+1}^{n_{pc}}{\int }_{a_{nepc}}^{a_{nl}}{f}_{nep2}^{\prime }M{\gamma }^{n-n_{\min }}nada=\frac{28}{15}\times 1.1282\times {7.1197}^{0.2544-c_{22}}\frac{KYMD}{2c_{22}-D+2}{\displaystyle \sum }_{n=n_{epc}+1}^{n_{pc}}{\gamma }^{n-n_{\min }}{a}_{nec}^{-c_{22}}{a}_{nl}^{1+c_{22}-\frac{D}{2}}-a_{nepc}^{1+c_{22}-\frac{D}{2}}{a}_{nl}^{D/2}.\end{array}$$

A.4. When the Asperity Frequency Level Belongs to $n_{pc}<n\le n_{\max }$

The real elastic contact area $A_{re}$, the real first elastoplastic contact area $A_{rep1}$, the real second elastoplastic contact area $A_{rep2}$, and the real full plastic contact area $A_{rp}$ are evaluated as$$\begin{array}{}\text{(A.13)}& A_{re}={\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\int }_0^{a_{nec}}M{\gamma }^{n-n_{\min }}naada=\frac{MD}{2-D}{\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\gamma }^{n-n_{\min }}{a}_{nec}^{2-D/2}{a}_{nl}^{D/2},\text{(A.14)}& A_{rep1}={\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\int }_{a_{nec}}^{a_{nepc}}M{\gamma }^{n-n_{\min }}naada=\frac{MD}{2-D}{\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\gamma }^{n-n_{\min }}{a}_{nepc}^{2-D/2}-a_{nec}^{2-D/2}{a}_{nl}^{D/2},\text{(A.15)}& A_{rep2}={\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\int }_{a_{nepc}}^{a_{npc}}M{\gamma }^{n-n_{\min }}naada=\frac{MD}{2-D}{\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\gamma }^{n-n_{\min }}{a}_{npc}^{2-D/2}-a_{nepc}^{2-D/2}{a}_{nl}^{D/2},\text{(A.16)}& A_{rp}={\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\int }_{a_{npc}}^{a_{nl}}M{\gamma }^{n-n_{\min }}naada=\frac{MD}{2-D}{\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\gamma }^{n-n_{\min }}{a}_{nl}^{2-D/2}-a_{npc}^{2-D/2}{a}_{nl}^{D/2}.\end{array}$$

And the elastic contact load $F_{re}$, the first elastoplastic contact load $F_{rep1}$, the second elastoplastic contact load $F_{rep2}$, and the full plastic contact load $F_{rp}$ are evaluated as$$\begin{array}{c}\text{(A.17)}& F_{re}={\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\int }_0^{a_{nec}}{f}_{ne}M{\gamma }^{n-n_{\min }}nada=\frac{MD}{3-D}{\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}\frac{4E{\pi }^{1/2}{G}^{D-1}}{3{1/{\gamma }^n}^D}{\gamma }^{n-n_{\min }}{a}_{nec}^{3-D/2}{a}_{nl}^{D/2},\text{(A.18)}& F_{rep1}={\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\int }_{a_{nec}}^{a_{nepc}}{f}_{nep1}^{\prime }M{\gamma }^{n-n_{\min }}nada=\frac{1.8667KYMD}{2+2c_{12}-D}{\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\gamma }^{n-n_{\min }}{a}_{nec}^{-c_{12}}{a}_{nepc}^{c_{12}-\frac{D}{2}+1}-a_{nec}^{c_{12}-\frac{D}{2}+1}{a}_{nl}^{D/2},\\ \text{(A.19)}& \begin{array}{c}{F}_{rep2}={\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\int }_{a_{nepc}}^{a_{npc}}{f}_{nep2}^{\prime }M{\gamma }^{n-n_{\min }}nada\\ =\frac{28}{15}\times 1.1282\times {7.1197}^{0.2544-c_{22}}\frac{KYMD}{2c_{22}-D+2}{\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\gamma }^{n-n_{\min }}{a}_{nec}^{-c_{22}}{a}_{npc}^{c_{22}-\frac{D}{2}+1}-a_{nepc}^{c_{22}-\frac{D}{2}+1}{a}_{nl}^{D/2},\end{array}\text{(A.20)}& F_{rp}={\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\int }_{a_{npc}}^{a_{nl}}{f}_{np}M{\gamma }^{n-n_{\min }}nada=\frac{MHD}{2-D}{\displaystyle \sum }_{n=n_{pc}+1}^{n_{\max }}{\gamma }^{n-n_{\min }}{a}_{nl}^{2-D/2}-a_{npc}^{2-D/2}{a}_{nl}^{D/2}.\end{array}$$

B. Calculation Procedure for the Fractal Parameters $D$ and $G$ of Equivalent Rough Surface

At present, there are many methods for calculating fractal parameters of fractal rough surfaces, and the structure function method is more commonly used [4]. For a fractal rough surface profile, the structure function is$$\begin{array}{}\text{(B.1)}& S\tau ={zx+\tau -zx}^2=\underset{-\infty }{\overset{+\infty }{{\displaystyle \int }}}P\omega e^{j\omega \tau }-1d\omega ,\end{array}$$where $\tau$ is an arbitrary increment of $x$. $$ refers to spatial average. $P\omega$ is the power spectrum of WM function. The discrete power spectrum of WM function can be approximated by a continuous spectrum, which is given by [50]$$\begin{array}{}\text{(B.2)}& P\omega =\frac{G^{2D-1}}{2\ln \gamma }\frac{1}{{\omega }^{5-2D}}.\end{array}$$

Substituting equation (B.2) into equation (B.1) yields$$\begin{array}{}\text{(B.3)}& S\tau ={zx+\tau -zx}^2=\frac{\Gamma 2D-3\sin 2D-3\pi /2}{4-2D\ln \gamma }{G}^{2D-1}{\tau }^{4-2D}=C{\tau }^{4-2D},\end{array}$$where $C=\mathrm{\Gamma }2D-3\sin 2D-3\pi /24-2D\ln \gamma G^{2D-1}$, and $\mathrm{\Gamma }x$ is the second type Euler Gamma function. Using the log algorithm on both sides of equation (B.3) yields$$\begin{array}{}\text{(B.4)}& \mathrm{l}\mathrm{g}S\tau =\mathrm{l}\mathrm{g}C+4-2D\mathrm{l}\mathrm{g}\tau .\end{array}$$

Therefore, there is a linear relationship between $\mathrm{l}\mathrm{g}S\tau$ and $\mathrm{l}\mathrm{g}\tau$. The fractal dimension $D$ can be obtained from the conversion relation between $D$ and the slop of the line $k$$$\begin{array}{}\text{(B.5)}& D=\frac{4-k}{2}.\end{array}$$

The relationship between parameter $C$ and the intercept of the line on the ordinate $B$ is $C={10}^B$. The characteristic length scale $G$ thus can be obtained as$$\begin{array}{}\text{(B.6)}& G={4-2D\cdot {10}^B\cdot \ln \gamma \Gamma 2D-3\cdot \sin \frac{2D-3\pi }{2}}^{\mathit{1}/\mathit{2}D-\mathit{1}}.\end{array}$$

When rough surface 1 comes into contact with rough surface 2, the structure function of the equivalent rough surface profile is given by [45].$$\begin{array}{}\text{(B.7)}& S_e\tau =S_1\tau +S_2\tau ,\end{array}$$where $S_1\tau$ and $S_2\tau$ are the structure functions of the profile of the rough surface 1 and the rough surface 2, respectively. Thus, the fractal parameters $D$ and $G$ of the equivalent rough surface can be estimated according to equations (B.4)∼(B.6).

References

[1] A. Adib, C. Baptista, M. Barboza, C. Marques, "Aircraft engine bleed system tubes: material and failure mode analysis," Engineering Failure Analysis, vol. 14 no. 8, pp. 1605-1617, DOI: 10.1016/j.engfailanal.2006.11.053, 2007.

[2] X. M. Guo, C. L. Xiao, H. Ge, H. Ma, H. Li, W. Liu, Z. H. Liu, "Dynamic modeling and experimental study of a complex fluid-conveying pipeline system with series and parallel structures," Applied Mathematical Modelling, vol. 109, pp. 186-208, DOI: 10.1016/j.apm.2022.04.003, 2022.

[3] X. D. Liu, W. Gao, Z. H. Gao, "Optimization of hoop layouts for reducing vibration amplitude of pipeline system using the semi-analytical model and genetic algorithm," IEEE Access, vol. 8, pp. 224394-224408, DOI: 10.1109/access.2020.3044087, 2020.

[4] Z. Chen, Y. Liu, P. Zhou, "A novel method to identify the scaling region of rough surface profile," Fractals Complex Geometry Patterns and Scaling in Nature and Society, vol. 27 no. 02,DOI: 10.1142/s0218348x19500117, 2019.

[5] Y. Chai, M. Chai, "Sealing failure and fretting fatigue behavior of fittings induced by pipeline vibration," International Journal of Fatigue, vol. 136,DOI: 10.1016/j.ijfatigue.2020.105602, 2020.

[6] Y. Yan, J. Zhai, P. Han, Q. Han, "A multi-scale finite element contact model for seal and assembly of twin ferrule pipeline fittings," Tribology International, vol. 125, pp. 100-109, DOI: 10.1016/j.triboint.2018.04.028, 2018.

[7] T. Z. Wen, F. Guo, Y. J. Huang, S. X. Jia, X. H. Jia, "Analysis of static sealing rules of foamed silicone rubber based on a porous media model," Cellular Polymers, vol. 39 no. 3, pp. 101-116, DOI: 10.1177/0262489319890076, 2020.

[8] G. R. Murtagian, V. Fanelli, J. A. Villasante, D. H. Ernst, H. A. Ernst, "Sealability of stationary metal-to-metal seals," Journal of Tribology, vol. 126 no. 3, pp. 591-596, DOI: 10.1115/1.1715103, 2004.

[9] Y. Ledoux, D. Lasseux, H. Favreliere, S. Grandjean, J. Grandjean, "On the dependence of static flat seal efficiency to surface defects," International Journal of Pressure Vessels and Piping, vol. 88 no. 11-12, pp. 518-529, DOI: 10.1016/j.ijpvp.2011.06.002, 2011.

[10] C. Lasseux, D. Lasseux, "Experimental leak-rate measurement through a static metal seal," Journal of Fluids Engineering, vol. 129 no. 6, pp. 799-805, DOI: 10.1115/1.2734250, 2007.

[11] Y. P. Shao, Y. X. Yin, S. Du, L. F. Xi, "A surface connectivity-based approach for leakage channel prediction in static sealing interface," Journal of Tribology, vol. 141 no. 6,DOI: 10.1115/1.4043123, 2019.

[12] M. Beghini, L. Bertini, C. Santus, A. Mariotti, G. Mariotti, "Partially open crack model for leakage pressure analysis of bolted metal-to-metal flange," Engineering Fracture Mechanics, vol. 144, pp. 16-31, DOI: 10.1016/j.engfracmech.2015.06.005, 2015.

[13] S. K. S. Kambhammettu, A. P. Chebolu, L. R. Chebolu, "A compressible porous media model to estimate fluid leak through a metal-elastomer interface," Transport in Porous Media, vol. 136 no. 1, pp. 191-215, DOI: 10.1007/s11242-020-01507-9, 2021.

[14] F. Pérez-Ràfols, R. Larsson, E. J. van Riet, A. Almqvist, "On the loading and unloading of metal-to-metal seals: a two-scale stochastic approach," Proceedings of the Institution of Mechanical Engineers-Part J: Journal of Engineering Tribology, vol. 232 no. 12, pp. 1525-1537, DOI: 10.1177/1350650118755620, 2018.

[15] M. H. Müser, W. B. Dapp, R. Bugnicourt, P. Sainsot, N. Lesaffre, T. A. Lubrecht, B. N. J. Persson, K. Harris, A. Bennett, K. Schulze, S. Rohde, P. Ifju, W. G. Sawyer, T. Angelini, H. Ashtari Esfahani, M. Kadkhodaei, S. Akbarzadeh, J.-J. Wu, G. Vorlaufer, A. Vernes, S. Solhjoo, A. I. Vakis, R. L. Jackson, Y. Xu, J. Streator, A. Rostami, D. Dini, S. Medina, G. Carbone, F. Bottiglione, L. Afferrante, J. Monti, L. Pastewka, M. O. Greenwood, J. A. Greenwood, "Meeting the contact-mechanics challenge," Tribology Letters, vol. 65 no. 4,DOI: 10.1007/s11249-017-0900-2, 2017.

[16] Z. Xianjun, Z. Haodong, W. Weipeng, W. Yanze, "Study on the compression-resilience and sealing performance of new metal-to-metal contact gasket," Journal of Pressure Vessel Technology, vol. 140 no. 2,DOI: 10.1115/1.4039122, 2018.

[17] Q. Zhang, X. Chen, Y. Zhang, X. Zhang, "An experimental study of the leakage mechanism in static seals," Applied Sciences, vol. 8 no. 8,DOI: 10.3390/app8081404, 2018.

[18] F. Pérez-Ràfols, R. Almqvist, A. Almqvist, "Modelling of leakage on metal-to-metal seals," Tribology International, vol. 94, pp. 421-427, DOI: 10.1016/j.triboint.2015.10.003, 2016.

[19] V. L. Popov, "Contact mechanics and friction," Physical Principles and Applications,DOI: 10.1007/978-3-662-53081-8, 2017.

[20] L. Violano, G. Violano, "On the effective surface energy in viscoelastic Hertzian contacts," Journal of the Mechanics and Physics of Solids, vol. 158,DOI: 10.1016/j.jmps.2021.104669, 2022.

[21] B. Bhushan, Contact Mechanics of Rough Surfaces in Tribology: Single Asperity Contact, 1998.

[22] B. Bhushan, "Contact mechanics of rough surfaces in tribology: multiple asperity contact," Tribology Letters, vol. 4 no. 1,DOI: 10.1023/A:1019186601445, 1998.

[23] G. Liu, Q. Wang, C. Lin, "A survey of current models for simulating the contact between rough surfaces," Tribology Transactions, vol. 42 no. 3, pp. 581-591, DOI: 10.1080/10402009908982257, 1999.

[24] J. R. Barber, M. Ciavarella, "Contact mechanics," International Journal of Solids and Structures, vol. 37 no. 1-2, pp. 29-43, DOI: 10.1016/s0020-7683(99)00075-x, 2000.

[25] H. Ghaednia, X. Wang, S. Saha, Y. Xu, A. Jackson, R. L. Jackson, "A review of elastic-plastic contact mechanics," Applied Mechanics Reviews, vol. 69 no. 6,DOI: 10.1115/1.4038187, 2017.

[26] J. A. Greenwood, "Contact of nominally fiat surfaces," Proc.r.soc.a, vol. 295, pp. 300-319, 1966.

[27] G. Zhao, S. X Li, Z. l Xiong, W. d Gao, Q. k. Han, "A statistical model of elastic-plastic contact between rough surfaces," Transactions of the Canadian Society for Mechanical Engineering, vol. 43 no. 1, pp. 38-46, DOI: 10.1139/tcsme-2017-0052, 2019.

[28] R. L. Streator, J. L. Streator, "A multi-scale model for contact between rough surfaces," Wear, vol. 261 no. 11-12, pp. 1337-1347, DOI: 10.1016/j.wear.2006.03.015, 2006.

[29] A. Majumdar, B. Bhushan, "Role of Fractal Geometry in Roughness Characterization and Contact Mechanics of Surfaces," Journal of Tribology-transactions of The Asme, vol. 112, pp. 205-216, 1990.

[30] Z. Chen, Y. Zhou, P. Zhou, "A comparative study of fractal dimension calculation methods for rough surface profiles," Chaos, Solitons & Fractals, vol. 112, pp. 24-30, DOI: 10.1016/j.chaos.2018.04.027, 2018.

[31] H. Zhu, S. Ge, X. Huang, D. Zhang, J. Liu, "Experimental study on the characterization of worn surface topography with characteristic roughness parameter," Wear, vol. 255 no. 1-6, pp. 309-314, DOI: 10.1016/s0043-1648(03)00215-1, 2003.

[32] J. Song, W. Wang, Y. Lang, Y. Luo, G. Luo, "Double rough surface contact model and finite element simulation based on fractal theory," Journal of Physics: Conference Series, vol. 1877 no. 1,DOI: 10.1088/1742-6596/1877/1/012016, 1877.

[33] Z. Zhao, H. Han, P. Wang, H. Ma, S. Yang, Y. Yang, "An improved model for meshing characteristics analysis of spur gears considering fractal surface contact and friction," Mechanism and Machine Theory, vol. 158,DOI: 10.1016/j.mechmachtheory.2020.104219, 2021.

[34] G. Lan, W. Sun, X. Zhang, Y. Chen, W. Li, X. Li, "A three-dimensional fractal model of the normal contact characteristics of two contacting rough surfaces," AIP Advances, vol. 11 no. 5,DOI: 10.1063/5.0045151, 2021.

[35] Y. Yuan, K. Zhao, K. Zhao, "A fractal model of contact between rough surfaces for a complete loading-unloading process," Proceedings of the Institution of Mechanical Engineers-Part C: Journal of Mechanical Engineering Science, vol. 234 no. 14, pp. 2923-2935, DOI: 10.1177/0954406220910440, 2020.

[36] D. J. Whitehouse, "Fractal or fiction," Wear, vol. 249 no. 5-6, pp. 345-353, DOI: 10.1016/s0043-1648(01)00535-x, 2001.

[37] A. Bhushan, B. Bhushan, "Fractal model of elastic-plastic contact between rough surfaces," Journal of Tribology, vol. 113 no. 1,DOI: 10.1115/1.2920588, 1991.

[38] S. Komvopoulos, K. A. Komvopoulos, "A fractal theory of the interfacial temperature distribution in the slow sliding regime: Part I-elastic contact and heat transfer analysis," Journal of Tribology, vol. 116 no. 4, pp. 812-822, DOI: 10.1115/1.2927338, 1994.

[39] S. Komvopoulos, K. Komvopoulos, "A fractal theory of the interfacial temperature distribution in the slow sliding regime: Part II-multiple domains, elastoplastic contacts and applications," Journal of Tribology, vol. 116 no. 4, pp. 824-832, DOI: 10.1115/1.2927341, 1994.

[40] W. Komvopoulos, K. Komvopoulos, "Contact analysis of elastic-plastic fractal surfaces," Journal of Applied Physics, vol. 84 no. 7, pp. 3617-3624, DOI: 10.1063/1.368536, 1998.

[41] Y. Etsion, I. Etsion, "Resolving the contradiction of asperities plastic to elastic mode transition in current contact models of fractal rough surfaces," Wear, vol. 262 no. 5-6, pp. 624-629, DOI: 10.1016/j.wear.2006.07.007, 2007.

[42] X. Huang, X. Huang, "A complete contact model of a fractal rough surface," Wear, vol. 309 no. 1-2, pp. 146-151, DOI: 10.1016/j.wear.2013.10.014, 2014.

[43] Y. Yuan, Y. Cheng, K. Gan, L. Gan, "A revised Majumdar and Bushan model of elastoplastic contact between rough surfaces," Applied Surface Science, vol. 425, pp. 1138-1157, DOI: 10.1016/j.apsusc.2017.06.294, 2017.

[44] Y. Yuan, J. Zhang, L. Zhang, "Loading-unloading contact model between three-dimensional fractal rough surfaces," AIP Advances, vol. 8 no. 7,DOI: 10.1063/1.5027437, 2018.

[45] K. Xu, Y. Zhang, L. Zhang, "A normal contact stiffness model of joint surface based on fractal theory," Computer Modeling in Engineering and Sciences, vol. 119 no. 3, pp. 459-479, DOI: 10.32604/cmes.2019.04677, 2019.

[46] Y. Wang, X. Zhang, S. Chen, Y. Chen, "Fractal loading model of the joint interface considering strain hardening of materials," Advances in Materials Science and Engineering, vol. 2019,DOI: 10.1155/2019/2108162, 2019.

[47] W. R. Chang, I. Bogy, D. B. Bogy, "An elastic-plastic model for the contact of rough surfaces," Journal of Tribology, vol. 109 no. 2, pp. 257-263, DOI: 10.1115/1.3261348, 1987.

[48] L. Etsion, I. Etsion, "Elastic-plastic contact analysis of a sphere and a rigid flat," Journal of Applied Mechanics, vol. 69 no. 5, pp. 657-662, DOI: 10.1115/1.1490373, 2002.

[49] J. L. Liou, C. M. Lin, J.-F. Lin, "A microcontact model developed for sphere- and cylinder-based fractal bodies in contact with a rigid flat surface," Wear, vol. 268 no. 3-4, pp. 431-442, DOI: 10.1016/j.wear.2009.08.033, 2010.

[50] A. Tien, C. Tien, "Fractal characterization and simulation of rough surfaces," Wear, vol. 136 no. 2, pp. 313-327, DOI: 10.1016/0043-1648(90)90154-3, 1990.