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

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 rst 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 rst and second elastoplastic deformation stages. Secondly, the size distribution function of contact spots at dierent frequency levels is derived. And the expressions of asperity critical frequency levels are rederived. Lastly, the feasibility and credibility of the proposed model are veried by comparison with other contact models and experimental data. e results show that when the variation of the material hardness is considered, the contact area of a single asperity in the rst 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. e proposed model is a modication 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.


Introduction
e 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][2][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][5][6]. Once the metalto-metal seals are out of work, leakage will be formed. Both stability and reliability of the aero-engine will be a ected. e 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 a ected by a variety of factors [7][8][9][10][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][13][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. us, 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-tometal seals [6,[16][17][18].
e 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. ese models mainly include numerical contact models, statistical contact models, and fractal contact models [15,[21][22][23][24][25]. e 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. e 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. e 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]. e 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][32][33][34][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. ey 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. e MB model predicts that the asperity first deforms plastically and then elastically during the loading process. is 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 loadingunloading contact model for rough surfaces [35,44], as well as a normal contact stiffness model for joint surfaces [45].
e 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. e 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. e 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. erefore, 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. e present work will lay the foundation for the subsequent analysis of the sealing performance of static metal seals.

Theory and Method
2.1. e WM Function. An engineering surface profile often appears random, multiscale, and disordered. e 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. e WM function is given as [37] z where z(x) 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. c 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, c is taken as 1.5. n min indicates the low cut-off frequency of the profile.

e Contact of a Single Asperity.
e 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. erefore, the first step is to obtain the "contact area-contact load" relationship of a single asperity.

e Existing Elastic Microcontacts.
According to the WM function, the profile of the fractal asperity at a certain level n before deformation is [37] 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/c n . e asperity height is [37] e 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 ω n of the asperity with a length scale l n is where d is the separation distance between the rigid flat surface and the equivalent rough surface. If the surface height z(x, y) follows the Gaussian distribution, separation distance d can be obtained by the relationship [37] where A r , A a are the total real contact area and the nominal contact area, respectively. erfc(x) represents the complementary error function. According to the Hertz theory [37], the contact area and contact load of the asperity within the elastic deformation stage can be obtained as a n � πR n ω n , where R n is the radius of curvature at the asperity summit, R n � |1/|d 2 z n /dx 2 | x�0 | � l D n /π 2 G D− 1 . E is the equivalent elastic modulus, 1/E � 1 − ] 2 1 /E 1 + 1 − ] 2 2 /E 2 , and E 1 , E 2 , ] 1 , and ] 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  Figure 2: e contact between a fractal asperity at level n and a rigid smooth plane, where l n , l nr , and l nt are the base diameter, the real contact diameter, and the truncation diameter, respectively [43].

e Existing Condition for Initial
Yield. e 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. e interference corresponding to the initial yield point is termed as the critical interference ω nec , and is given as [47] where K is the hardness factor of the softer material and is given by 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. e elastic critical contact area of the asperity a nec can be obtained as [43] a nec � πR n ω nec � 1 π According to equations (7), (9), and (10), we can get the expression of the elastic critical contact load of the asperity f nec , which is

e 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. is transition stage, with the transition interference ratio of 1 < ω n /ω nec ≤ 6(a nec < a n ≤ a nepc ), is known as the first elastoplastic deformation stage. e first elastoplastic critical contact area a nepc and the first elastoplastic contact load f nep1 of the asperity are derived by Yuan et al. [43] based on the ME model, which are 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: H G1 a n � c 11 Y a n a nec c 12 , a nec < a n ≤ a nepc , (14) where c 11 , c 12 are the parameters need to be solved. Equation (14) should satisfy two boundary conditions: 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/3)KH 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/3)KH. us, the parameter c 11 can be obtained as Substituting equations (12)- (14) and (17) into equation (16) yields And parameter c 12 can be obtained as erefore, the contact load of a single asperity in the first elastoplastic deformation stage is revised as 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 < ω n /ω nec ≤ 110(a nepc < a n ≤ a npc ). Some researchers named this stage the second elastoplastic deformation stage. e 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]: e 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: where c 21 and c 22 are the two parameters need to be solved.
And equation (23) should also satisfy two boundary conditions: 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 .

e Existing Plastic Microcontacts.
When ω n /ω nec > 110(a n > a npc ), the asperity undergoes full plastic deformation. e contact area and contact load of the asperity in this deformation stage are [43] a n � 2πR n ω n , 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.

e Revised Size Distribution Function of Contact Spots.
In the MB model, the size distribution function of contact spots is defined as [37] where a l denotes the largest contact area of asperity. And the total real contact area is [37] A r � According to equation (2), the period of asperity at frequency level n is T n � 2π/2πc n , and the period of asperity at frequency level n + 1 is T n+1 � 2π/2πc n+1 . erefore, T n+1 /T n � 1/c. 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 f ne (a n ) The first elastoplastic deformation: 0 The second elastoplastic deformation: f nep1 '(a n ) a nec a nepc a npc a n f nep2 '(a n ) f np (a n ) Plastic deformation: Elastic deformation: Figure 3: e diagram of the deformation law of a single asperity and the relationship between each critical contact area.
Advances in Materials Science and Engineering asperity at frequency level n + 1, which is n n+1 (a) � cn n (a). Letting the size distribution function of the asperity at the initial frequency level n min be n n min (a) � Mn(a), then the contact spot size distribution function at each frequency level can be obtained as 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 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+1)l /a nl � 1/c 2 . en, the largest contact area of the asperity at any frequency level would be a nl � (1/c 2(n− n min ) )a l . us equation (35) can be reexpressed as And the relationship between a nl and A r can be obtained as e parameter M can be determined according to the relationship A n min r + A n min +1 r + · · · + A n max r � A r . Substituting equation (36) into it yields M � (1 − (1/c))/(1 − (1/c) n max − n min +1 ).

e Revised Critical Asperity Frequency Levels.
For a certain frequency level n, the relationship between the interference ω n and the asperity height δ n is ω n ≤ δ n [46]. If the height of an asperity is less than its elastic critical interference, inelastic deformation occurs. us, the condition of elastic deformation for an asperity given by the YC model is δ n ≤ ω nec [43]. However, according to Etsion and Etsion [41], the condition for elastic deformation of an asperity should be ω n ≤ ω 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 . 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/c 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. us, the condition we give to judge the elastic deformation of the asperity is where ω nl is the largest interference of the contact spot at frequency level n. According to equation (6), equation (39) can be re-written as Substituting equations (10) and (37) into the above inequality, the elastic critical frequency level can be obtained from the following equation: 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: 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 . at 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.

e 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 ≤ n ≤ n ec , only elastic deformation takes place. And the real contact area A r1 and contact load F r1 are obtained as 6 Advances in Materials Science and Engineering When the asperity frequency level belongs to n ec < n ≤ n epc , elastic or the first elastoplastic deformation can take place. e real contact area A r2 is obtained as And the contact load F r2 is obtained as When the asperity frequency level is n epc < n ≤ n pc , elastic deformation, the first elastoplastic deformation, or the second elastoplastic deformation will occur. e real contact area A r3 is evaluated as If n min ≤ n ≥ n ec If n ec < n ≤ n epc If n epc < n ≤ n pc If n pc < n ≤ n max Critical contact areas of the asperities a nec , a nepc , a npc (Equation (10), (12), (21)) The largest contact area of the asperities a nl (Equation (38)) Critical asperity frequency levels n ec, n epc, n pc (Equation (41) And the contact load F r3 is evaluated as When the asperity frequency level belongs to n pc < n ≤ n max , all of the four types of deformation can take place. e real contact area A r4 is calculated as And the contact load F r4 is estimated as       e detailed calculation results of the real contact area and contact load for the rough surface are provided in Appendix A. e total real contact area and contact load of all asperity levels are estimated by e nondimensional forms of A r and F r are defined as [43] where the nominal contact area A a is given by A a � L 2 , L � 1/c 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.

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 × 10 − 9 m, D � 1.5. Let the sampling length and the resolution be L � 0.5 × 10 − 4 m (L � 1/c n min ), L s � 1.5 × 10 − 9 m (L s � 1/c n max ), respectively. e 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 × 10 − 8 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. e calculation procedure is given in Appendix B. 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. e value of nondimensional total real contact area A * r (A * r � 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]. e 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]. e 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]. e 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. e 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. e 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. e characteristic length scale is G � 2.5 × 10 − 9 m, and the nondimensional total real contact area A * r 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. e 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. e conditions for the proposed model and the YC model to solve the elastic critical level are ω nl ≤ ω nec and δ n ≤ ω nec , respectively. And the asperity interference will not be greater than its height, i.e., ω nl ≤ δ 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. e relationships between nondimensional total real contact area A * r 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 value. However, as the A * r 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. erefore, the critical frequency levels in the proposed model decrease with the increase of A * r . e 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 will consequently not affect the calculation results of the critical frequency levels. 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. e frequency level n is 33, and the value of A * r 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.
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.
e "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. e asperity frequency level is 38, and the value of A * r 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 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 12(a) indicates the comparison of the proposed model and the MB model on the relationship between nondimensional total contact load F * r and nondimensional total real contact area A * r for different roughness parameters. e comparison of the proposed model and the MB model on the relationship between nondimensional total real contact area A * r 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/ �� � A a increases, the F * r value required to produce a specific A * r value increases. e same changing trend is also found in the presented model. In addition, the F * r value of the proposed model required to produce a particular A * r 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 increases. is 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/ �� � 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 value. It does not change with the change of G/ �� � A a value either. Figure 13(a) shows the comparison of the present model and the MB model on the relationship between nondimensional total contact load F * r and nondimensional total real contact area A * r for different values of D. Figure 13(a) indicates that, in the MB model, the F * r value required for a particular A * r 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 value required to produce a specific A * r value increases. However, in the proposed model, the F * r value needed for a specific A * r 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 value corresponding to the same A * r 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

12
Advances in Materials Science and Engineering fraction of elastic contact area increases with the increase of A * r 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 value increases. is 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 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. us, we do not give a comparison between the proposed model and the YC model regarding the A r /A a ∼ F r /(A a E) relationship and the A re /A r ∼ A r /A a relationship. Nevertheless, the idea of contact modeling on rough surfaces in reference [43] still has a vital reference significance. 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/ �� � A a � 1 × 10 − 10 , Y/E � 0.05. It should be   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. e 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.

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. e main conclusions are as follows: (1) e 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) e size distribution functions of the contact spots at different frequency levels are derived. e expressions of asperity critical frequency levels are rederived. e 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) e 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. e proposed model is helpful for the analysis of the sealing performance of static metal seals under different contact pressures in our subsequent studies.

Nomenclature a l :
Area of the largest contact spot a n : Contact area of a single asperity at frequency level n a nec : Elastic critical contact area of the asperity at frequency level n a nepc : First elastoplastic critical contact area of the asperity at frequency level n a nl : Area of the largest contact spot at frequency level n a npc : Second elastoplastic critical contact area of the asperity at frequency level n A a : Nominal contact area of the rough surface A nr : Real contact area of the asperities at frequency level n A r : Real contact area of the rough surface  A.4. When the Asperity Frequency Level Belongs to n pc < n ≤ n max . e 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 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

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 where τ is an arbitrary increment of x. 〈 〉 refers to spatial average. P(ω) is the power spectrum of WM function. e discrete power spectrum of WM function can be approximated by a continuous spectrum, which is given by [50] P(ω) � .

(B.6)
When rough surface 1 comes into contact with rough surface 2, the structure function of the equivalent rough surface profile is given by [45]. where S 1 (τ) and S 2 (τ) are the structure functions of the profile of the rough surface 1 and the rough surface 2, respectively. us, the fractal parameters D and G of the equivalent rough surface can be estimated according to equations (B.4)∼(B.6). Data Availability e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors have no conflicts of interest to disclose.