Mechanics Analysis of Rough Surface Based on Shoulder-Shoulder Contact

: Proper methods and models for mechanical analysis of rough surface can improve the theory of surface contact. When the topography parameters of two rough surfaces are similar, the contact should be considered shoulder-shoulder rather than top-top. Based on shoulder-shoulder contact and fractal characteristics, the geometric model for asperity and contact mechanics model for rough surfaces are established, and the deformation of asperity, the real contact area and contact load of sealing surface are discussed. The effects of contact pressure p and topography parameters (fractal dimension D and fractal roughness G ) on the variation of porosity and contact area ratio A r / A 0 are achieved. Results show that with the increase of p , larger D and smaller G corresponds to larger initial porosity but faster and larger decrease of porosity; with the increment of D , porosity increases first and then decreases, and smaller G corresponds to larger porosity reduction; as G becomes bigger, porosity increases, and larger D corresponds to larger porosity difference and change. With the addition of p , A r / A 0 increases, and the variation of A r / A 0 is closer to linearity and less at smaller D and larger G ; with the increase of D , A r / A 0 increases gradually, and the growth rate is bigger at smaller G and bigger p ; as G becomes bigger, A r / A 0 declines, and it declines more gently at smaller D and p . The influence of D on A r / A 0 is greater than that of G . The results can provide the theoretical basis for the design of sealing surfaces and the research of sealing or lubrication technologies of rough surfaces.


Introduction
At a microscopic level, the rough surface consists of numerous asperities of different sizes which may be caused by form errors, waviness and roughness. Although the waviness affects the contact of surfaces [1], the influence of surface roughness is more significant for rough machined fractal sealing surfaces. When two contacting surfaces are under load, these asperities will deform, affecting the pressure bearing and sealing performance of the surfaces. Mechanical analysis of rough surface is the main content of surface contact theory, which involves real contact area, contact pressure and the relationship between them. A variety of contact mechanics models were established by predecessors, promoting the continuous development of the surface contact theory.
Hertz studied the contact mechanics of two spherical elastomers without friction and laid the theoretical foundation of modern contact mechanics [2]. Greenwood and Williamson combined the simulated surface with Hertz contact theory and established the GW model [3] based on statistical characterization parameters. Based on GW model, Greenwood [4] and Bush [5] further modified and improved the height distribution, shape and curvature radius of asperities. Assumed that all the asperities are purely plastic deformed, Pullen [6] established the corresponding plastic contact model of rough surfaces. On the principle of constant volume of asperities during plastic deformation, Chang [7] proposed the CEB model to analyze the elastic and plastic contact characteristics of rough surfaces. However, the CEB model did not consider the elastic-plastic deformation, so the results are discontinuous. Zhao et al. re-analyzed the deformation process of asperities under load [8,9] and put forward a contact model including the elastic, elastic-plastic and plastic deformation of asperities. Ciavarella [10], Vakis [11], Tian [12] and Song [13] took the asperity base deformation caused by the asperity deformation into consideration and established various statistical contact models for the interaction of asperities. In addition, many researchers established the loading-unloading contact model [14,15], rough surface contact stiffness model [16,17] and viscoelastic contact model [18,19] to meet various working requirements. Based on the GW statistical model and probability density function of height distribution for asperity, Xiao [20] established the total stiffness model of rough surface contact by considering the contact state of elastic, elastic-plastic and purely plastic deformation. Zhu proposed an unlubricated contact model for random rough surfaces [21] and analyzed the contact between equivalent rough surfaces and rigid planes based on statistical theory. Contact mechanical models of rough surface with statistical characterization were established under the premise that "The contact area is less than the cross-sectional area of the profile opening area", however, this assumption is not consistent with the principle of volume conservation during deformation. Moreover, the statistical parameters in these models are greatly influenced by the sample length and resolution of the instrument, which is highly scale-dependent and not universal.
Majumdar and Bhushan [22] established the rough surface contact model (M-B model) based on fractal theory for the first time. They deduced the elastic and plastic contact deformation of the single asperity, but they did not consider the elastic-plastic contact problem. Wang and Komvoulos modified the M-B model by considering the effect of the expansion factor of the size distribution domain of the micro contact on the deformation of asperities [23]. Yan and Komvoulos initially put forward a three-dimensional rough surface contact model with fractal characteristics [24] and obtained the total contact load and contact area of the rough surface by derivation. Based on the modified M-B model, Zhu [25], Miao [26] and Tian [27] derived the relationship between the contact area of a single asperity and its normal deformation. Liu [28] thought that the contact area of a single asperity increases with the increment of its normal deformation, and when the contact area is less than the critical contact area the deformation is plastic, otherwise it is elastic, which is consistent with the deformation mechanism of the M-B model. Morag [29] and Liou [30] proposed using the opening size of asperity's base to express the contour curve. They defined the index n of asperities at different levels to describe the deformation of asperity, and pointed out that under certain contact loads, the high-level asperity will transform to a low level. Huang [31] combined fractal theory with a finite element model to study the contact characteristics of three-dimensional rough surface under thermal-mechanical coupling. By deducing the W-M function in the spherical coordinate system and conducting numerical simulation, Liu [32] obtained the profile of rough spherical surface, established the fractal rough spherical contact model and analyzed the influence of fractal parameters on the relationship between contact load and contact area on rough spherical surface. According to the contact mechanism of spherical surface, Yuan [33] established the contact model of rough spherical surface, derived the analytical expression of real contact area and contact load, and obtained the contact pressure distribution of rough spherical contact area at different deformation stages. In the study of contact mechanics of rough surfaces based on fractal theory, fractal parameters are used to describe rough surfaces, it can overcome the influence of measuring instrument resolution and sampling length. However, the conclusion that "asperity appears plastic, elastic-plastic and elastic deformation in turn" is against reality.
Zhang [34] fully considered the influence of surface roughness and analyzed the contact characteristics of rough surface when establishing the contact model. Wang [35] established a loading-unloading contact analytical model between two cylindrical rough surfaces based on fractal theory and investigated the effect of parameters in the model on the loading-unloading contact performance. In addition, some scholars use the shoulder-shoulder contact model [36][37][38][39] rather than top-top contact model when studying the contact force, contact stiffness and friction heat of two surfaces. However, these models do not have explicit applications for shoulder-shoulder contact, and do not specifically study the fractal machined surface that conforms to the W-M function. Therefore, the research on the contact mechanics for rough surfaces, especially the shoulder-shoulder contact, needs in-depth investigation.
In the present paper, the contact state and porosity of rough sealing surfaces are analyzed; the geometric model for asperity and contact mechanics model for rough surfaces are established; and the deformation of asperity, the real contact area and contact load of rough sealing surfaces are discussed. Based on these, the effects of contact pressure and surface topography parameters (fractal dimension D and fractal roughness G) on the porosity and contact area of sealing surface are studied.

Contact State of Rough Surface
For the research of contacting rough surfaces, the two surfaces are mostly simplified. The smoother surface is usually simplified as a smooth rigid plane, and the contact state of the two surfaces is shown in Figure 1. However, in reality the sealing interface is not completely smooth and there are many asperities on both surfaces, as shown in Figure 2.  The red areas in Figures 1 and 2 represent the total volume of pores and skeletons between the two surfaces, and the skeleton volume can be obtained by function integral.
Porosity φ is the proportion of pore volume in total volume, subscript r-f is the contact between rough and flat surfaces, and r-r is the contact between two rough surfaces. The porosity in these two cases is as follows: In contact analysis, the object of study is usually a single asperity or a pair of asperities. The surface contact shown in Figures 1 and 2 can be represented by the largest asperity, as shown in Figures 3 and 4.  According to the W-M function of fractal surface and the two different contact states, the initial contact state of two rough surfaces can be simulated by MATLAB (as shown in Figure 5), and then the initial porosity can be obtained to verify which contact state is more practical.

The Initial Porosity of Rough Sealing Surface
Zuo [40] used fractal parameters to describe rough surfaces and found they can comprehensively characterize the overall, local and internal structure of the surface. For the fractal machined rough surface, W-M function can be used to describe the contour curve of any asperity before deformation where D is the fractal dimension of the surface, G is the fractal roughness, l is the base diameter of the asperity which can be given by l = 1/γ n , γ n is the spatial frequency of the asperity. When n = nmin, the base diameter of the largest asperity is l = 1/γ nmin , when n = nmin+1, the base diameter of the corresponding asperity is l= 1/γ nmin+1 , etc.
For the asperity satisfying Equation (1), its maximum height h = G D-1 l 2-D corresponds to the maximum surface roughness, and the corresponding l is the base diameter of the maximum asperity.
The number of asperities on the sealing surface is related to the maximum contact area [23], which satisfies where n(a) is the distribution density function for the contact area, subscript 1 and 2 represents two rough surfaces, respectively, n(a) takes the smaller value of n(a)1 and n(a)2, a is the contact area of asperity, aL is the maximum area of micro contact, ψ is the correction factor of Ar/aL, which is the ratio of the real contact area to the maximum area of micro contact. This can be obtained from the transcendental equation [41]: For the sealing surface with known surface topography parameters, the nominal contact area A0 is equal to the sum of the cross-sectional area of the largest asperity and all other smaller asperities on the surface: where aL1,2 is the largest area of micro contact, aL1,2 = πl1,2 2 /4, m 2 ; A0 takes the larger value of A01 and A02.
So far, the numerical expression can be used to calculate the porosity for sealing surface under certain topography parameters. Equation (5) is the initial porosity when a plane contacts a rough surface, and Equation (6) is the initial porosity when two rough surfaces contact.
In Equation (5), h is the height of the largest asperity on the rough surface; in Equation (6), h' is the maximum distance between asperities on two surfaces, h1 and h2 is the maximum height of asperity on two rough surfaces, respectively. Once the rough surface is determined, the specific value of fractal dimension D and fractal roughness G, h', h1 and h2 can be obtained through assigning or measurement. According to Equation (1), the corresponding l, l1 and l2 can be obtained, so aL, aL1 and aL2 can also be achieved, and then the initial porosity can be calculated.
Sealing surfaces with specific topography parameters (D: 1.3, 1.35, 1.4, 1.45, 1.5, 1.55, 1.6; G:25 × 10 −11 m) were selected to achieve the simulation value for initial porosity based on W-M function and MATLAB programming. Then these values were compared with the numerical calculation value obtained by Equations (5) and (6). The results are shown in Figure 6.  When the fractal parameters of the two rough surfaces are the same, the calculated value based on the rough-rough contact is closer to the simulated value; the maximum relative error is less than 5%, while the relative error between the calculated and simulated value based on the rough-flat contact is 10.5~15.5%. Therefore, when the topography parameters of two surfaces are the same or similar, the contact should be considered shoulder-shoulder rather than top-top. Shoulder-shoulder contact should be used for the modeling and mechanical analysis of the rough surface.

Porosity of Rough Sealing Surface after Loading
Under the action of contact pressure, the asperities on the rough surface will deform (i.e., produce the compression w'), which will change the porosity of the sealing interface. Taking the calculation unit of sealing interface (the total volume is V0 = L × L × h', L is the length and width of the calculation unit, and h' is the height of the calculation unit, namely the height of red area in Figure 2) as the research object, after loading the height of the sealing interface will change from h' to h' − w'. Then, the volume of the calculation unit changes from According to the principle that the skeleton volume remains unchanged before and after deformation, the skeleton volume is still Vs = (1 − φ0) × V0 and the pore volume is Vp' = V0' − Vs, then the porosity of sealing surfaces becomes: where φ0 is the initial porosity of the sealing interface.

Model Simplification of the Shoulder-Shoulder Contact
When two rough surfaces contact, the actual contact situation is shown in Figure 7a, which can be simplified as the contact of the two largest asperities, as shown in Figure 7b. For the single asperity, its contour curve satisfies Equation (1), and the curve equations of two asperities in Figure 7b can be expressed as: H is the distance between the two surfaces and d is the distance between the rotation axes of the two asperities.
The curvature radius of the two asperities can be obtained by: represents the x coordinate value at the contact point. The equivalent curvature radius at the contact point is R = 1/(1/R1 + 1/R2). The stress analysis at the contact area is shown in Figure 8. When Fn = Fa × cosθ, w = w′ × cosθ, the contact area perpendicular to the direction of Fn is equal to the contact area at inclined direction (1-2 direction) multiplied by cosθ.

Area Solution for the Shoulder-Shoulder Contact
The coordinates of point 1 and point 2 in Figure 8 can be obtained by setting up the two curve equations of asperities simultaneously (x1 < x2 is the default condition), According to Figure 9, coordinate transformation is carried out, and finally the inclined plane contact part, as shown in Figure 10, is obtained and its area can be calculated.
Therefore, the intersection area S equals to where Ld is the linear distance between point 1 and point 2,

Analysis of Deformation State
With the increment of applied load, the deformation of asperity will increase, resulting in elastic, elastic-plastic and purely plastic deformation.
When analyzing the contact mechanics of asperities on rough surface, the assumptions are as follows: ① Deformations only occur on asperities on the rough surface, while the macro base does not deform. The volume of asperities remains unchanged before and after deformation.
② In the process of surface contact, there is no interaction between asperities, and neither surface material strengthens.
(1) Elastic deformation According to Hertz theory, the contact force in the direction of Fa during elastic deformation is When the average contact pressure on the contact surface is less than the yield limit σy of the material, the asperities are in an elastic state. Under this condition, the elastic contact pressure can be expressed by pn. When pn equals to σy, the inception of yield occurs and the critical area aec can be given by 2 y 1.5 Purely plastic deformation The asperity yields fully when the average pressure is equal to 3σy [42], in this condition σ After simplification, the critical plastic deformation area apc can be expressed as (3) Elastic-plastic deformation When aec ≤ a ≤ apc, the asperity has elastic-plastic deformation. In this contact state, the relationship between the contact area and contact pressure becomes extremely complex. Considering that variations of the contact area and contact pressure should be continuous and smooth at the critical point of initial yield and be fully plastic, a template function f(a) [8] can be constructed. The template function is only used to distribute the weight of elastic area and plastic area in the elastic-plastic interval and does not involve the derivation conditions and process of the model, so it does not affect the compatibility of the model.
Then the average contact pressure in the elastic-plastic stage is given by (14)

Real Contact Area and Contact Load of Rough Surface
For two rough surfaces in contact when the number of asperities on the interface satisfies Equation (2), its nominal contact area A0 can be expressed by Equation (4).
When aL < aec, all the contact asperities are in elastic deformation, therefore, the real contact area Ar is the elastic contact area Are and the contact force Fn is the elastic contact force Fne, which can be given by  , aL is the contact area of the largest asperity; Ar is the larger value between Ar1 and Ar2; Fn takes the larger value of Fn1 and Fn2.
When aec ≤ aL ≤ apc, the contact asperities are in elastic-plastic deformation. In this scenario, the total contact area Ar is the sum of elastic contact area Are and elastic-plastic contact area Arep, and the total force Fn is the sum of elastic contact force Fne and elastic-plastic contact force Fnep When aL > apc, the contact asperities are in plastic deformation. In this condition, the total contact area Ar is the sum of elastic contact area Are, elastic-plastic contact area Arep and plastic contact area Arp, and the total force Fn is the sum of elastic contact force Fne, elastic-plastic contact force Fnep and plastic contact force Fnp.   (20) Ar is the larger value between Ar1 and Ar2; Fn takes the larger value of Fn1 and Fn2. Although the expressions of the real contact area Ar are the same in the three states, since the largest asperity has different deformation properties, aL in each formula is different, so the actual contact area is not the same. Therefore, when calculating the real contact area, it is necessary to determine the deformation properties of the largest asperity.

Model Validation
To verify the accuracy of the established model, results of Ar* − p* in Kucharski experiment [43], Y-K modified model [24] and Sahoo finite element analysis [44] were extracted, and the same material parameters and topography parameters were used in the present established model for calculation. These four sets of data were compared. The relationship between Ar* and p* in different methods are shown in Figure 11. It can be seen that the data of this model is in good agreement with Kucharski experimental data and Sahoo finite element analysis data.

Solution of Porosity and Real Contact Area
The calculation steps of porosity and real contact area are as follows: (1) For the rough sealing surface, determine the comprehensive elastic modulus E and the yield strength σy of the softer surface.
(2) By assigning or measuring, the detailed value of fractal dimension D and fractal roughness G, the maximum distance h' between the rough surfaces and the maximum height of the asperity on each rough surface h1 and h2 can be achieved. The corresponding l1 and l2 can be obtained according to Equation (1), A0 can be achieved through Equation (4), and the initial porosity can be gained by Equation (5). Under the action of contact pressure, there exists compression w', and the porosity can be calculated by Equation (10).
(3) Calculate the critical elastic deformation aec and critical plastic deformation apc of the asperity through Equations (11) and (12), obtain the contact area aL of the largest asperity on the sealing surface by Equations (16), (18) and (20) according to the given specific pressure pc, then compare aL with aec and apc to determine the deformation properties of the asperity.
(4) Obtain the real contact area Ar according to Equations (15), (17) and (19), and then the specific value of Ar/A0 is achieved.
The The porosity and contact area ratio Ar/A0 of the above sealing surfaces were calculated, and the influence of contact pressure (the unit is MPa) and surface topography parameters D and G (the unit is 10 −11 m) were analyzed.

Porosity
(1) Influence of contact pressure p on porosity The change of porosity with contact pressure p is shown in Figure 12. Since asperities will be squeezed and filled into the pores of the sealing interface under the action of contact pressure, porosity decreases with the increase of p. Larger D and smaller G corresponds to larger initial porosity, but the corresponding decrease of porosity is faster and larger in amplitude with the increase of p. When the rough surface has large fractal dimension and small fractal roughness, the surface is smooth and the roughness is small, so the ratio of the initial pore volume to the total volume of the sealing surface is high. However, with the increase of the contact pressure and the deformation of the asperities, the decrease rate of the pore volume is greater than that of the total volume of the sealing surface, so the porosity continues to decline. (2) Influence of topography parameters on porosity Seen in Figure 13(a1-a3), with the increase of fractal dimension D, the overall variation trend of porosity increased first and then reduced (when p = 0.1 MPa, all the porosity increases; the larger the p is, the earlier and faster the porosity decreases), and the smaller the G is, the larger the porosity decreases. At the beginning, with the increase of D, many small pores within the rough sealing surfaces are shown, so the porosity increases. With the further increase of D and p, it is easier for the asperities to deform and fill into the pores, so the porosity decreases. Seen in Figure 13 (b1-b3), with the increment of G, all the porosity increases, larger p corresponds to lower porosity and larger D corresponds to larger difference and amplitude of porosity change. Larger G means rougher surface and coarser asperity, and it is much more difficult for the asperity to deform and fill into the pores, so the porosity is bigger.

Contact Area Ratio
(1) Influence of contact pressure on Ar/A0 As shown in Figure 14, with the addition of contact pressure p, Ar/A0 increases, and the variation of Ar/A0 is closer to linearity with smaller changing amplitude at smaller fractal dimension D and larger fractal roughness G. With the increase of p, the deformation of the asperity increases, resulting in the augment of the real contact area. Smaller D and larger G means the asperity has smaller curvature radius, so the value and increment amplitude of real contact area after compression is smaller, therefore, the value and variation range of Ar/A0 is smaller. Contact area ratio Ar/A0 Figure 14. Change of Ar/A0 with contact pressure p.
(2) Influence of topography parameters on Ar/A0 As can be seen from Figure 15a1-a3, with the increase of D, Ar/A0 increases gradually, and the growth rate is bigger at smaller G and bigger p. Seen in Figure  15b1-b3, as G becomes bigger, Ar/A0 declines, and it declines more gently at smaller D and smaller p. The influence of surface topography parameters on Ar/A0 is mainly reflected by the curvature radius of surface asperity and the real contact area after deformation. From Figure 15 it is found that the influence of fractal dimension D on Ar/A0 is greater than that of fractal roughness G.
(a1)  Through the research, the contact condition and porosity of rough sealing surface can be better understood. In practice, by combining these findings and the specific working requirements (better lubrication or lower leakage), we can choose a rougher or smoother surface so as to determine the corresponding surface processing method. On the basis of the given surfaces, by measuring the relevant parameters and calculating according to the established model, the contact pressure under which the sealing surface has no percolation or has appropriate lubrication can be calculated.

Conclusions
(1) The contact state and porosity of rough sealing surfaces were analyzed. Through numerical calculation and simulation, it was found that when the topography parameters of two rough surfaces are the same or similar, the contact should be considered shoulder-shoulder rather than top-top.
(2) The geometric model for the shoulder-shoulder contact and contact mechanical model for rough surfaces were established, and the deformation of asperity, the real contact area and contact load of rough sealing surface were discussed. By comparison with other methods, the present contact mechanics model was verified.
(3) The change rules of porosity φ with contact pressure p and topography parameters were obtained: Porosity decreases with the increase of p; larger fractal dimension D and smaller fractal roughness G correspond to larger initial porosity but faster and larger decrease of porosity; with the increment of D, porosity increases first and then decreases on the whole, and smaller fractal roughness G corresponds to larger porosity reduction; as the fractal roughness G becomes bigger, all the porosity increases, and larger D corresponds to a larger difference and change of porosity .
(4) The change rules of contact area ratio Ar/A0 with contact pressure p and topography parameters were achieved: With the addition of contact pressure p, Ar/A0 increase, and the variation of Ar/A0 is closer to linearity with smaller changing amplitude at smaller fractal dimension D and larger fractal roughness G; with the increase of D, Ar/A0 increases gradually, and the growth rate is bigger at smaller G and bigger p; as G becomes bigger, Ar/A0 declines, and it declines more gently at smaller D and p; the influence of fractal dimension D on Ar/A0 is greater than that of fractal roughness G.
(5) The findings of this paper could be further verified by comparing the porosity obtained by the established model with the results obtained by 3D reconstruction and finite element analysis, and the leakage or lubrication of sealing surface also could be studied, all which need further in-depth research.