A Revised Continuous Observation Length Model of Rough Contact without Adhesion

: The real contact area of rough surfaces has significant importance in many engineering applications, such as tribology, wear, lubrication and seals. A continuous observation length-de-pendent mechanic model of rough contact without adhesion is proposed, which assumes that the rough surface is divided into ideal subplanes. However, the model ignores the elastoplastic deformation of asperities, and the standard deviation of ideal subplanes’ heights is assumed to vary linearly with continuous observation length, which is not precise for all the surface fractal dimensions. In this work, a revised continuous observation length model is proposed with elastic, elastoplastic and fully plastic stages. The expressions of force and real contact areas are derived. For surfaces with different fractal dimensions, the quadratic polynomial, cubic polynomial and power relation-ships between standard deviation and observation length are proposed, respectively . In addition, the influences of the dimensionless observation length, fractal dimension and equivalent elastic modulus on the real contact areas in different contact stages are also analyzed. It can be concluded that the quadrate real contact area decreases as the dimensionless observation length decreases, which can be applied to the percolation theory for leak seal problems. J.W.; software, J.W. and M.L.; validation, J.W. and G.X.; formal analysis, J.W.; investigation, J.W.; resources, L.Z.; data curation, J.W. and G.X.; writing—original draft preparation, J.W.; writing—review and editing, M.L.; visualization, L.Z.; supervision, L.Z.; project administration, L.Z. All authors have read and agreed to the published version of the manuscript.


Introduction
The dimensions of these deviations and the surface topography vary in a big range from macro (millimeter) to micro (nanometer) and even down to the level of interatomic distance. Consequently, the contact tends to be localized and incomplete when the two surfaces are engaged with each other at the beginning. The real contact area is of great importance in many engineering applications, such as structural adhesives, tribology, wear, lubrication, as well as seals.
Several theoretical methods were employed to deal with surface contact. An important pioneer in this field was Hertz, who developed a classic elastic contact model for both macroscopic and microscopic contact problems. A common critical criterion of plastic contacts was proposed by Tabor [1]. The pioneering work of the rough surface contact was contributed by Greenwood and Williamson, who proposed a single-length elastic contact model (GW model) for a rough surface based on Hertz's theory and statistical method [2]. Whitehouse and Archar assumed that the curvature radius of the asperity is dependent on its height and proposed the WA model based on the GW model. The CEB model was also based on the GW model with the assumption that the volume of plastic deformation is conserved [3]. Xu et al. developed statistical models based on the GW model and extended the range of application of the classic models to the case of near complete contact [4]. Majumdar and Bhushan developed a very popular analytic method (MB model) for the contact problem of fractal rough surfaces, in which the real contact area

The OLD Model
The contact area is different when two contact surfaces are studied at different observation lengths. The two surfaces appear to be in complete contact at a large observation length because the surface roughness is neglected due to its small length. The incomplete contact emerges in reality with a small observation length, even in extremely tight contact conditions.
As the observation length decreases, the incomplete contact becomes gradual. While the observation length tends to a very small value (atomic length), the contact area of the two surfaces tends to zero [25].
Recently, Liu et al. proposed an observation length-dependent (OLD) mechanic model to calculate the real contact area for a frictionless contact between an elastic-perfectly plastic rough surface and a rigid plane [29]. In the model, the surface profile was determined by the WM function, and the contact of two rough surfaces can be idealized as a rigid plane pressed on a deformable rough surface, as shown in Figure 1a. The sample surface with the sample length L0 was regarded as a discrete surface composed of ideal subplanes with the observation length λ, λ0 < λ < L0, as shown in Figure 1b. The ideal subplanes were regarded as asperities of the same curvature radius at a certain observation length when the contact force and the real contact area were calculated, as shown in Figure 1c. Therefore, the contact of a rigid plane and a deformable rough surface can be idealized as the contact between a rigid plane and deformable ideal subplanes. The height of the asperity is the average height of the asperity at the observation length division plane. Assuming that the heights of the ideal subplanes follow the Gaussian distribution in all observation lengths, the standard deviation of the ideal subplane height can be obtained. To obtain the relationship between the observation length and standard deviation, a linear relationship was assumed (as shown in Section 2.1.1). However, there are three problems in this model: (a) assuming that the standard deviation of the ideal subplane height varies linearly with successive observed lengths is not accurate for all the surface fractal dimensions; (b) no dimensionless parameters were used; (c) the model only investigated the elastic and plastic deformation of the asperities but ignored the elastic-plastic deformation. This leads to errors in solving the real contact area and contact force. To solve these problems, a revised continuous observation length model is presented.

The Standard Deviation
According to the multiscale characteristics of a rough surface, it is assumed that the rough surface is isotropic, and the W-M function can be used to define the surface profile [5,14]. The height of the surface profile z(x) is where D is the fractal dimension of the surface profile, and 1 < D < 2, G is the fractal roughness parameter of the surface; γ n determines the frequency spectrum of the surface roughness, and γ ≥ 1.0. For most typical surfaces, γ = 1.5. nl is the minimum frequency index, and nl = logγ(1/L0). As can be seen in Figure 1, the height of an arbitrary ideal subplane is where h(λ) is the height of the ideal subplane; ω is the contact deformation; and d(λ) is the distance between the reference plane and the rigid plane after contact. λ is the observation length.
The standard deviation of the ideal subplane is where x represents the coordinate of the ideal subplane; h(λ) is the average height of the ideal subplane.
Assuming that the height of the ideal subplane obeys the Gaussian distribution on all observation lengths, the probability density of ideal subplane height distribution at observation length λ is where σ is the standard deviation of the ideal subplane height distribution, h(λ) = 0. In the OLD model, Liu et al. assume a linear relationship applies, and (λ) can be expressed as 0 ( )=k y σ λ λ + , where k is the linear fitting slope; y0 is the standard deviation of the surface profile at the minimum observation length.
It can be seen that the standard deviation was assumed to vary linearly with the observed lengths, which is not accurate. Additionally, the observation length is not dimensionless to eliminate the influence of sample length.

The Revised Continuous Observation Length Model
According to the OLD model, it is assumed that the rough surfaces are composed of ideal subplanes with the same observation length and asperity curvature radius, and every subplane obtains a specific height with regard to statistics. The multilength contact mechanics model of the rough surface is established based on the above assumptions.

The Standard Deviation
Surface profiles with various D and G, listed in Table 1, are used to calculate the standard deviation (λ) of the ideal subplane height using Equations (1) and (2). According to Equation (1), D has a significant effect on the height of the rough surface profile. To simplify the problem, it is assumed that the observation length λ is continuous. The variation of the standard deviation σ with λ * at different D is shown in Figure 2, where λ * is the dimensionless observation length (the observation scale) and λ * = λ/L0. The observation length is nondimensionalized to eliminate the influence of sampling length.
The surface profiles with various D, listed in Table 1, are used to calculate the standard deviation σ of the ideal subplane height. λmin is the minimum observation length. According to the W-M function, D has a significant effect on the height of the rough surface profile [38].
It can be seen from Figure 2a that σ increases monotonically with the decrease in λ * . Moreover, the σ of the surface with a larger D is smaller when the λ * is the same. It is obvious that different D has a great impact on the order of magnitude of σ. In order to facilitate observation, the standard deviation is magnified by a certain multiple, as shown in Figure 2b. From Figure 2b, when D = 1.5, the relationship between λ * and σ is approximately linear. When D < 1.5, the rising rate of σ slows down at a small λ * . When D > 1.5, the rising rate speeds up. This is because the value of D has an influence on the variation relationship between σ and λ * . It can be seen from Equation (1) that z(x) is directly proportional to G (D−1) , and z(x) becomes smaller as D becomes bigger. When the λ * is small, the change of z(x) is significant due to the increase in D. Therefore, the standard deviation σ also changes greatly, resulting in different rising rates.
After analyzing the calculation results under different D, the fitting is divided into two parts: D ≤ 1.5 and D > 1.5. The best fitting methods under different D are also different.
When D ≤ 1.5, linear fitting (σ = y0 + B1λ * ), quadratic polynomial fitting (σ = y0 + B1λ * + B2λ *2 )and cubic polynomial fitting (σ = y0 + B1λ * + B2λ *2 + B3λ *3 ) are adopted, respectively. The fitting curves of λ * and σ are shown in Figure 3. It can be seen that the linear fitting curves have a large deviation from the data points compared to the polynomial fitting curves. The quadratic fitting curves and the cubic fitting curves are more consistent with the data points. Thus, the polynomial fitting is better than the linear fitting.  When D > 1.5, it can be seen from Figure 2b that the variation between λ * and σ is obvious. The other values of D, 1.7, 1.8 and 1.9, are taken for detailed analysis. At the small λ * , the rising rate of σ becomes much larger, so the cubic fitting relationship is no longer good. Finally, the ExpDec2 function (σ = y0 +B1e -x/A 1 +B1e -x/A 2) is used for fitting. The fitting curves of λ * and σ are shown in Figure 4.  The coefficient of determination R 2 is important to evaluate the fitting effect. The closer the value of R 2 is to 1, the better the fit of the model. The relationship between R 2 and D in different fitting methods is shown in Figure 5. It can be concluded that: (a) when D ≤ 1.5, the power fitting is not suitable; (b) when D ≤ 1.5, the R 2 of all the three fitting methods is above 0.99, which is well fitted. Among the three methods, there is little difference between the two polynomial fitting methods, both of which are better than linear fitting, and the difference gradually decreases with the increase in D; (c) when D = 1.6, the R 2 of the three fitting methods is basically the same. However, from the figure, the quadratic polynomial fitting is more in line with the requirements; (d) when D ≥ 1.7, the R 2 of the linear fitting and polynomial fitting decreases with the increase in D. When D = 1.9, the R 2 of the linear fitting, quadratic fitting and cubic fitting is 0.82, 0.94, 0.96, respectively. This indicates that the fitting is poor. In addition, the R 2 of the ExpDec2 function fitting is bigger than 0.99, which is better than the other fitting methods.   Figure 5. The relationship between R 2 and D.

Establishment of Contact Mechanics Model
The W-M function provides the asperity profile at any λ [38]: Therefore, the asperity curvature radius is The critical contact depth ωe of the plastic contact is derived by Jackson based on Von Mises' yield criterion, as shown below: where K is the hardness coefficient related to the Poisson ratio ν of the material by K = 0.454 + 0.41 ν. H is the hardness of the softer material, H = 2.8σY, and σY is the yield strength. ϕ is the material property, ϕ = H/E. E is the equivalent modulus of Hertzian elasticity and equals where E1, E2, ν1 and ν2 are the Young's modulus and Poisson's ratios of two contacting surfaces, respectively.

Contact Areas and Mechanics
At the macrolevel, the machined surface is a smooth surface, but it is a rough surface at the microlevel. Therefore, the contact between the two surfaces, which are assembled together, is actually the contact between asperities. According to the ideal subplane asperity deformation, there are three contact states: elastic, elastoplastic and plastic deformation. The elastic contact state indicates that the ideal subplanes are fully in elastic contact. The plastic contact state indicates that the plastic contact dominates, and the elastic contact is extremely small. The elastoplastic deformation stage is in between and is divided into two stages. In the OLD model, the elastoplastic deformation was ignored, which makes the results inaccurate. Kogut and Etsion used the finite element method to analyze the mechanical behavior of elastoplastic contact of a single spherical asperity [8].When 0 < ω < ωe, the asperity is at elastic deformation, and the ideal subplane height of elastic deformation should be d(λ) < h(x, λ) ≤ d(λ) + ωe. When ωe < ω < 6ωe, it is the first stage of elastoplastic deformation, and the height should be d(x, λ) + ωe < h(x, λ) ≤ d(x, λ) + 6ωe. When 6ωe < ω < 110ωe, it is the second stage of elastoplastic deformation, and the height should be d(x, λ) + 6ωe < h(x, λ) ≤ d(x, λ) + 110ωe. In addition, when ω > 110 ωe, the contact is completely plastic.

The Areas in Different Stages
The real area Ar of contact is The real plastic contact area Arp is The dimensionless areas in different stages can be, respectively, written as where A0 is the nominal area, A0 = L 2 0 . The total dimensionless contact area A * r is obtained by accumulating the different stages' contributions. Hence, from Equations (14)-(17)

The Relationship between F and ωe of Single Asperity
The relationship between the contact of a single asperity and the critical contact depth is [8] 1 3 2 2 re where ωe is the critical contact depth; Fre, Frep1, Frep2 and Frp are the contact forces of a single asperity in four stages; Fec is the critical contact force of the elastic contact. The critical contact force is

The Contact Force of the Surfaces in Different Stages
The surface contact forces consist of elastic, elastoplastic and plastic parts and can be described as where Pre, Prep and Prp are the elastic, elastoplastic and plastic contact forces with the observation length λ, which can be calculated as the sum of forces in different stages of asperities.
The mean contact stress σ0 is The dimensionless contact forces are The total dimensionless contact force is * * * * * c re rep1 rep2 rp P P P P P = + + + The surface contact force PN and the nominal area A0 are usually known for a certain contact system and will not change with the observation length λ. The separation d(λ) at any observation length can be solved by the numerical method. Consequently, the real contact area, elastic contact area and plastic contact area can be calculated with d(λ).

Numerical Results and Discussion
In order to study multilength contact characters, the following common parameters are selected: the surface fractal dimension D = 1.3, 1.4 and 1.5; the equivalent modulus of Hertzien E = 50 GP, 75 GPa and 100 GPa; the fractal roughness parameter G = 1 × 10 −12 m; the Poisson's ratio ν = 0.3 and the yield strength σY = 345 MPa. The mean contact stress σ0 is 10 MPa. By numerical analysis in Python, the changing rules between the contact area and λ * are obtained. The numerical analysis is based on the above equation through circular iterations and judgments.

The Influence of the Fractal Dimension
The fractal surface profiles obtained by Equation (5) with G = 1 × 10 −12 m and D = 1.3, 1.4, 1.5 are plotted in Figure 6. After abundant and repetitive simulations, we found out that a high value of D results in more surface details.  Figure 6 shows that A * r can be reduced by the decrease in λ * with the same D. Additionally, the observation scale decreases when A * r begins to reduce with the rising D. It is worth pointing out that A * r grows as D increases with the same λ * . This is because the surface with a smaller value of D reaches the incomplete contact.
The effects of λ * on the dimensionless contact areas are shown in Figure 7 with G = 1 × 10 −12 m and E = 100 GPa. When λ * decreases, the plastic contact dominates, and otherwise, the elastic contact dominates. The dimensionless elastic contact area A * re decreases as λ * reduces, and in the other three stages, the areas rise to the peak before decreasing to a small value. The maximum areas reduce with the increase in D. Moreover, the λ * where peaks occur reduce with the increase in D. In the elastic deformation stage, when λ * > 0.5, the proportion of A * re is 1, and there is only elastic deformation, as shown in Figure 8. The proportion gradually decreases when a) D = 1.3 and λ * < 0.4; b) D = 1.4 and λ * < 0.02; c) D = 1.5 and λ * < 0.0006. This is because the elastic deformation transforms into the first elastoplastic deformation with the decrease in λ * . For the elastoplastic deformation stage, the trend is similar to that of their dimensionless areas, and the value of the peaks increases as D increases. For the plastic deformation stage, the proportion of the plastic dimensionless area is 0 when a) D = 1.3 and λ * > 0.002; b) D = 1.4 and λ * > 9 × 10 −5 ; c) D = 1.5 and λ * > 7 × 10 −6 . Additionally, when λ * is small enough, there is only plastic deformation.
From the above discussion, it can be seen that there is only elastic deformation at a large λ * . However, with the decrease in λ * , the first elastoplastic and the second elastoplastic contact states become the main contact states successively. Therefore, it is more consistent with reality to consider elastic-plastic deformation.

The Influence of the Equivalent Modulus
In order to study the effects of the equivalent modulus of elasticity, it is assumed that D = 1.3, G = 1 × 10 −12 m, ν = 0.3. Figure 9 is the dimensionless real contact area versus λ * with various E:E = 50 GPa, E = 75 GPa, E = 100 GPa. As shown in Figure 9, the dimensionless real contact areas decrease with λ * . When λ * < 0.2, the areas with various E are nearly the same. Additionally, when λ * > 0.2, a high value of E produces a big slope. Figure 10 presents the real contact areas in different stages versus λ * with various E. The area A * re decreases with λ * , and a high value of E produces a big slope. As λ * reduces, the areas A * rep1 , A * rep2 and A * rp rise to the peak before decreasing to a small value. Additionally, the maximum λ * where the area begins to reduce is different with different E. For example, when the area A * rep2 begins to reduce, E = 50 GPa and λ * = 7 × 10 −5 , E = 75 GPa and λ * = 2 × 10 −4 , E = 100 GPa and λ * = 1 × 10 −4 .  Figure 11 shows the influence of the equivalent modulus of elasticity on multilength contact characteristics. The area A * re decreases with λ * , and a high value of E produces a big slope. Similarly, as λ * reduces, the proportion of A * rep rises to the peak before decreasing to a small value, and the maximum λ * , where the peaks exist, becomes bigger when the value of E becomes higher. As E increases, the peak becomes lower in the first stage of elastoplastic contact, but it remains approximately unchanged in the other stage. The proportion of plastic contact area reduces as λ * increases, and when λ * > 0.001, the proportion is greater than 0.999.

Model Comparison
The  Figure 12. In order to compare with the Persson model, the magnification is defined [26]. The magnification ζ is defined as ζ = λ0/λ where λ0 is the critical value of observation length. When ζ = λ0/λ, the two contact surfaces are in full contact, and A * r = 1. When ζ = λ0/λ, the two contact surfaces are in full contact, and A * r = 1. When λ0 < λ, the contact is incomplete, and A * r < 1. . This is because the present model replaces the contact areas of the asperities with the ideal subplane areas and obtains a relatively large result. The contact area of the present model is smaller compared to the Persson model at a smaller observation length (ζ > 10 1.5 ). This is because Persson assumed that when the local contact is in a plastic state, the plastic contact area will no longer change as the observation length decreases. However, in the present model, after the contact enters into a plastic state, the local surface is still in an incomplete contact. It was suggested that the minimum observation length of the incomplete contact can reach the atomic level.
To compare the present model more specifically with the Liu model, Figure 13 is obtained. It shows the error ratio between the Liu model and the present model. It can be seen that the error fluctuates within 0~44.8%. For example, when log ζ = 1.5, the error ratio is 44.6. It can be seen that whether to ignore the elastic-plastic deformation has a large impact on the contact area, especially when the observation length is small. Therefore, the present model is more accurate than the Persson model, and it is meaningful.

Conclusions
In this work, a revised continuous observation length model of rough contact without adhesion based on the statistics method is proposed. The W-M function is used to simulate the rough surface and establish the ideal subplanes. The relationship between the standard deviation and the observation length at different fractal dimensions is fitted in different methods. Moreover, the transition from elastic deformation to plastic deformation is also considered. The mathematical models of forces and areas in different states are proposed, and the change of the areas with λ * are discussed. Additionally, the present model is compared to the Persson model and the Liu model to verify this model. The detailed conclusions are as follows: • When the fractal dimension D is different, the fitting relationship between the standard deviation and observation length is also different. For D ≤ 1.5, the polynomial fitting is the best. For D = 1.6, the quadratic polynomial fitting is more in line. For D ≥ 1.7, the ExpDec2 function fitting is better than the other fitting methods. The larger the D, the larger the error of the linear fitting.

•
As λ * reduces, the areas A * r and A * re decrease, and in the other three stages, the areas rise to the peaks before decreasing to a small value. The values of the peaks and the values of λ * where peaks occur are related to D and E. As λ * reduces, the proportion of A * re decreases, the proportion of A * rep rises to the peak before similarly decreasing to a small value, and the proportion of A * rp increases. When λ * reduces to a certain value, the proportion of A * rp approaches 1, and the proportion of A * rep1 is larger than A * rep2 . • As D increases, the area A * r increases. The values of the peaks and λ * where peaks occur reduce with rising D. As E increases, the area A * r drops quicker when λ * decreases. Additionally, the values of the peaks and λ * where peaks occur decrease with the increase in E.

•
The dimensionless real contact area of the present model is in good agreement with the Persson model and the Liu model in terms of variation patterns. • The present model can accurately describe the contact characteristics of rough surfaces and the monotonicity and continuity in contact processes. The research results provide a theoretical basis for analyzing the contact characteristics on rough surfaces and designing a contact surface topography.

The Future Work
In the next stage, the contact force and the contact area of the sealing device of engineering equipment will be studied. To ensure the accuracy of the data, multiple specimens will be used for contact performance testing. Firstly, the surface morphology of the experimental sample of the sealing device is measured to obtain the 2D contour profile. The fractal parameters of the contact surface are calculated, such as the fractal dimension D and the fractal roughness G. Secondly, the experimental setup will be set up, and normal loading experiments will be carried out to obtain the results of contact force, contact depth and contact area. Then, the present theoretical model, material parameters and fractal parameters will be input into Python for programing the calculation. Numerical results are obtained for the contact force and contact area as the contact depth increases. Finally, the experimental results are compared with the theoretical numerical results to verify the rationality of the present model. After modifying the model according to the experimental results, the model can be directly used to calculate the contact force without any experiment, which will save time and costs considerably. Data Availability Statement: The study did not report any data.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
The notations that were frequently used throughout the article are listed.