Modeling and simulation of dynamical contact of asperities between flat rough surfaces

Dynamic contact analysis on the level of surface asperities is a powerful means to study the wear behavior of surface and the physical contact properties of interface. This article established an analytical model of asperity contacts between nominal rough planes during sliding process based on the statistical contact analysis of rough surfaces. The analytical expression describes the mathematical relationship between the total number of asperity contacts and its main influencing factors including the nominal contact area, sliding distance, the spacing between two contact planes, and surface topographical parameters. A simulation model of dynamic contact of asperities on rough planes was established to rapidly calculate the number of asperity contacts during the sliding process. The comparison between the calculation results and simulation results of the number of asperity contacts obtained through multifactor and multilevel analysis verified the rationality of the analytical model.


Introduction
Since machining surfaces are composed of randomly distributed asperities with different heights, the contact of rough planes in movement will lead to the elastic or elastic-plastic dynamic contact of the asperities. Studies of rough surface contact were mainly focused on the features of the pure contact (static contact), including the contact area, deformation, and stress distribution of single asperity or multiple asperities under a certain external force. Relevant study methods included several theoretical modeling methods, 1 the numerical analysis methods, [2][3][4] and experimental study methods. 5,6 Current studies of surface dynamic contact were mainly focused on the stress-strain of asperities during the sliding contact process 7,8 and unstable friction (dynamic friction) caused by the changes in the micro-contact spacing during the surface sliding process. [9][10][11] The study of dynamic contact parameter of asperities on rough planes during sliding process, i.e. total number of asperity contact, was not reported. From the microscopic perspective, fatigue wear of sliding friction surfaces is the accumulated fatigue damage caused by repeated elastic contacts of asperities. Therefore, total number of asperity contact during sliding process is an important factor in microscopic fatigue wear of the surfaces. Moreover, because the current and heat of two contact surfaces is delivered 1 School of Control Engineering, Northeastern University at Qinhuangdao, Qinhuangdao, China through asperity contact, the dynamic contact parameter of asperities during sliding process can be used to analyze the dynamic contact resistance and thermal contact resistance of the two surfaces during the sliding movement. The study of dynamic contact parameter of rough surfaces can provide new research idea for modeling of the friction and wear behavior of contact surfaces and the physical contact properties of interface.
In order to establish the dynamic contact model of rough planes, it is necessary to obtain the rational model descriptions of rough surfaces. In the past decades, the modeling of the contact between rough surfaces was extensively studied. Due to the randomness and anisotropy of the surface topography, the theoretical modeling is very difficult. The classical statistical model for the analysis of the elastic and elastic-plastic hybrid contact between a rough surface and a smooth surface was proposed by Greenwood and Williamson (GW), 1 it is assumed that the rough surface is composed of hemispherical asperities with the same curvature radius, and that the height of hemispherical asperities follows the Gaussian distribution. In addition, it is also assumed that the deformations of different asperities are independent from each other during the asperity contact. The model was amended by Whitehouse and Archard 12 and Onions and Archard. 13 With the development of surface topography measurement techniques and numerical analysis techniques, based on GW model, Bush et al., 14 Gibson,15 and McCool 16 established a new elastic asperity contact model, in which the asperities were treated as elliptical paraboloids with randomly oriented elliptical contact areas. Kotwal and Bhushan 17 established a statistical non-Gaussian surface contact model to enhance the application of the statistical contact model. Majumdar and Bhushan 18,19 established a fractal contact model of rough surfaces, in which fractal geometrical parameters were used to describe the surface topography, and the uncertainty of topographical parameters was avoided. The development of the computer technology allows the digital treatment and computer simulation of the surface contour [20][21][22] to obtain the digital surface topography, in which the shape, size, and height distribution of the asperities can be arbitrarily designated to more accurately indicate the actual surface topography. Digital topography of different contact surfaces can also be used to analyze the surface asperity contacts and obtain the contact contours. Although statistical contact models, compared to other contact models, have some limitations and a certain deviation from real rough surfaces, a number of simple topographical parameters and simple mathematical expressions are adopted to predict the real contact area, contact load, and the parameters of other surface contact properties, and such prediction results have been proved by many theoretical and practical research. 23 Compared to other contact models, statistical contact models have theoretical and practical significances. Therefore, statistical contact models are still widely used in the surface micro-contact analysis of engineering machinery. In this article, based on the statistical contact theory, the asperity contact model of nominal rough planes was established. The rough surface asperity was approximated to the spherical body with the same curvature radius. Thus, the contact between two planes could be considered as the contact among a series of spherical bodies with different heights.
In this article, the dynamic contact model was established according to the following two steps. First, based on the statistical analysis theory, a surface asperity model was established to analyze the static elastic contact parameters of two rough planes, including the contact area and the number of the contacted asperities. Second, based on the fact that the transient contact parameter of rough planes was the same during sliding process, the dynamic contact model of orderly asperities was established to obtain the analytical expression of the number of asperity contacts through statistical integration.
To validate the dynamic contact model, a simulation model of the dynamic contact of asperities between rough planes was established in this article. In the simulation model, based on the established mathematical model of rough planes, the relative sliding process of surfaces was simulated, and the contact number of asperities during the sliding process was cumulatively calculated. The comparison results between the theoretical value and the simulated value of the number of asperity contacts verified the rationality of the theoretical analysis model.

Statistical contact modeling of nominally flat rough surfaces
The topography of nominally flat rough surfaces, as shown in Figure 1(a), is composed of a series of randomly distributed asperities. The surface topography can be characterized with a group of three-dimensional (3D) statistical parameters obtained through topography measurement. These 3D statistical parameters include the arithmetical mean m of asperity height, the standard deviation s of asperity height, the distribution function c(z) of asperity height (asperity height z is the distance from the top of an asperity to the reference plane), the mean radius R of asperities, the surface density D of asperities, and so on. As the elastic contact between rough surfaces occurs at the tips of the asperities on reference planes and the tip edge of the asperities is circular arc, the rough surfaces can be represented by an asperity model which is defined by a series of elastic spheroids, as shown in Figure 1(b). The radii of the spheroids are equal to the mean radius R of asperities on the rough surface, and the height of each spheroid is the same as the height of the corresponding asperity. The asperity model in Figure 1(b) will be used to analyze the elastic contact of nominally flat rough surfaces.
When two flat rough surfaces are in the elastic contact, the contacts occur on part of the asperities on the surfaces, as shown in model I of Figure 2. The rough surfaces are constructed by the asperity model of Figure 1(b). The contact spots (the top view of model I) present the contact state of asperities, including the number of contact asperities and the contact area of each contact asperity. According to the statistical contact model proposed by Greenwood and Tripp, 24 the contact of two rough surfaces can be represented as the contact between a rigid plane and an equivalent rough surface, as shown in model II in Figure 2. The contact spots of model II are consistent with those of model I. The equivalent rough surface (the lower surface in model II) is defined with the composite standard deviation s# of asperity height and the composite mean radius R# of asperities of the two rough surfaces in model I. s# and R# are expressed by equations (1) and (2) as follows where s 1 and s 2 are the standard deviations of asperity height distribution of the two original rough surfaces where R 1 and R 2 are the mean radii of asperities of the two original rough surfaces. Let n be the number of asperities and c(z) be the distribution function of asperity height; then the number of contact asperities m and the real contact area A of model II can be obtained by equations (3) and (4) as described by GW 1 where h is the distance between the reference planes.
The upper rigid plane of model II can be further equivalently represented as a rigid flat surface with uniform and orderly asperities, as shown as the upper surface of model III, which is based on the following conditions: (a) the height of the asperities is uniform, and their top is leveled with the original rigid plane in model II and (b) the radius and surface density of the asperities are consistent with the asperities of the lower surface of model II. The two upper surfaces of model II and model III have the same contact effect.
The asperities on the lower surface of model II can also be orderly distributed in lines without changing the topography parameters of the surface, which brings about the lower surface of model III. This equivalent transformation guarantees the contact parameters of model II and model III are in consistency.
As shown in Figure 2, the number of contact asperities (m) and the contact area of each contact asperity (also the total contact area A) of model III and model II are the same, although the contact spots of model III are orderly distributed, and those of model II is stochastic. The three models in Figure 2 have the equivalent contact since they present identical contact states of asperities. The orderliness of model III makes it efficient to analyze the dynamic contact of two rough surfaces.
Dynamic contact modeling of nominally flat rough surfaces for calculation of number of asperity contacts during sliding process As shown in Figure 3, when a nominally flat rough surface slides on another nominally flat rough surface, the asperities on the interface contact each other randomly due to the stochastic surface topography. The contact state of asperities (shown as contact spots) on the contact interface at every moment during sliding can be treated as the static elastic contact of nominally flat rough surfaces represented as the models in Figure 2. Therefore, the contact parameters (m, A) at every moment during sliding can be calculated by equations (3) and (4). If the two nominally flat rough surfaces are machined uniformly and the parameters of topography are consistent throughout the surfaces, the contact parameters (m, A) calculated by equations (3) and (4) for every instant remain constant during sliding. That is to say, the contact of asperities for every instant during sliding process can be evaluated by the static elastic contact model III with the identical surface topography parameters. Thus, when a sliding rough surface with the nominal area S slides on a fixed rough surface for the distance l, as shown in Figure 3, the dynamic contacts of asperities during sliding process can be described by a Model III when its upper rigid surface (with nominal area S) slides on the lower surface for the distance l, as shown in Figure 4. Figure 4 shows the arrangement of the asperities on the upper rigid surface of model III with the nominal area S (l 3 w). The asperities are packed closely in the first column from the left margin. The other asperities are equidistantly arranged in columns along the length direction. If n is the number of asperities on the upper surface, D is the surface density of asperities of the upper surface, R# is the mean radius of the asperities, i is the number of asperities in one column, and j is the number of columns, then Substituting equations (5) and (6) into equation (7), the number of columns j will be expressed as   When the upper rigid surface slides on the lower surface for the distance l, as shown in the top view in Figure 4, each asperity column will also slide over the length of l. The trajectory of each asperity column will form a rigid plane that contacts with the lower surface, and the contact state is equivalent to that of model II in Figure 2. The number of the contacted asperities is m in model II, so the number of contacts between asperities during the stroke of one column is also m. Since there are totally j columns on the upper surface, the total number N 1 of contacts between asperities in the sliding distance l will be Substituting equations (3) and (8) into equation (9) gives Substituting equation (5) into equation (10) gives Equation (11) presents the number N 1 of dynamic contacts between asperities of the two flat rough surfaces when sliding for one cycle length l. For a given sliding stroke L, the number N of dynamic contacts between asperities through the sliding process will be After substituting equation (11) into equation (12), it turns out The distribution function of asperity height z, c(z) in equation (13), usually follows the normal distribution with the arithmetic mean of m and standard deviation of s# (s# is expressed by equation (1)). It can be seen in Figure 5 that the height of asperities is more likely in the exponential distribution in the range of h ; +N; h is the spacing distance between the two reference planes on the contact rough surfaces; therefore, c(z) in equation (13) can be expressed approximately by a exponential function as And the relative error of the approximation is not more than 9.3% in the range of h;+N. Then, the integration expressions in equation (13) can be calculated analytically as Substituting equation (15) into equation (13) leads to Equation (16) presents the analytical expression of the number N of elastic dynamic contacts between asperities of two flat rough surfaces during sliding process. As shown in equation (16), the number of elastic contacts of asperities depends on the surface topography parameters R, D, and s; nominal contact area S; sliding distance L; and spacing h between two rough surfaces. Since the elastic contact states of asperities are the same during the forward and backward movements of the sliding surface when it slides on the fixed surface, equation (16) is applicable to both one-way and reciprocating sliding process.

Mathematical descriptions of nominal rough planes
Two rough surfaces and the 3D coordinate system are shown in Figure 6(a). For the upper surface, the length is l and the width is w. The contact area between the upper and lower surfaces is S = lw. The spacing between the surfaces is h. The sliding distance of the upper surface is L. As described as Figure 1, the asperities of rough surfaces in Figure 6(a) are composed of a series of spherical segments with the same radius and different heights. For the upper and lower surfaces, the radii of asperity are R 2 and R 1 , respectively; the distribution density of asperity is D 2 and D 1 , respectively; the standard deviation of the asperity height distribution is s 2 and s 1 , respectively. As shown in the top view of asperities (Figure 6(b)), the asperities are the circles with different radii, and the asperity position can be determined by the coordinate values of its vertex in the x-axis and y-axis position of the rough surface. The size of each asperity can be determined by the apex height z. The greater the height z, the greater the radius of spherical segment c, as shown in Figure 6(c). Therefore, all the asperities of the upper and lower surfaces can be uniquely determined by their 3D coordinate values of the vertices. The asperities on the upper and lower surfaces are indicated by the arrays of (X, Y, Z) and (x, y, z), respectively. For the surface asperities are randomly distributed, the position coordinates (x and y) of asperities are the random variables in the length direction and the width direction. The height z is also a random variable with normal distribution. As shown in Figure  5, the normal distribution is determined by the topographical parameters, m (the arithmetic mean of asperity height) and s (the standard deviation of asperity height). The position coordinates and height coordinates of the asperities are generated with the random number generation function and the normal distribution function in MATLAB program.
Determination method of the dynamic contact of surface asperities during sliding process When the sliding surface (the upper surface in Figure 6(a)) slides along the fixed surface (the lower surface in Figure  6(a)) for the length L, each asperity on the sliding surface also slides for the length L, and the sliding trajectory of each asperity forms the contact domain of the asperity, as shown in Figure 7. The asperity only contacts with the lower asperities which are located in the contact domain of the asperity. The range of the contact domain of the asperity is shown in Figure 8. When the position coordinates of the asperities on the lower surfaces meet equations (17) and (18), it means that the asperities are in the contact domain  where X and Y are the position coordinates of an asperity on the upper surface, c is the bottom radius of the asperity, and x and y are the position coordinates of an asperity on the lower surface. Let Z be the height of the asperity on the upper surface and z be the height of the asperity on the lower surface. h is the spacing between the two planes, as shown in Figure 6(a). If the upper asperity contacts with the lower asperity, it means that In conclusion, only when the asperities on the lower surface satisfy equations (17)- (19), they are contacted with the asperities on the upper surface during sliding process.

Simulation process of the dynamic contact of surface asperities
The process of the simulation calculation was mainly composed of two loops, as shown in Figure 9. One loop is used to define the contact domain of the asperities on the upper sliding surface. Another loop is used to determine whether the asperities on the lower fixing surface are in the contact domain and whether they are contacted with the upper asperities. Ultimately, the dynamic contact between asperities during sliding process is cumulatively calculated to obtain the dynamic contact parameter, i.e. total number of contacts between asperities.

Simulation experiments
The topography of the machined rectangular sample was measured with Taylor Hobson surface topography instrument (Figure 10(a)) to obtain various statistical parameters of surface topography, including height parameter, functional parameters, spatial parameters, hybrid parameters, feature parameters, and so on. With the statistical parameters, automatic fitting with 3D microscopic surface topography is available, as shown in Figure 10(b). For the testing rectangular surface, the 3D topographical height parameters (the average height m and the standard deviation s) and feature parameters (the average radius of the asperities R and the distribution density of asperities D) are respectively 0.6107 mm, 0.7607 mm, 85.4 mm, and 111.6690 mm 22 .
Based on these topography parameters, the following simulation experiments a-f with single variable were designed to analyze the influence of each factor (as shown in equation (16)) on the number of asperity contacts between two rough planes, wherein the topographical parameters (R, D, and s) of the upper and lower surfaces are the same. Since the rough surface was randomly generated with the statistical topography parameters in the simulation model, even the same topography parameters will give different simulation   results within a certain range for every calculation. Therefore, the simulation calculation was performed for 50 times with each set of parameters, and the average value was selected as the final simulation result, that is, the total number of contacts.  Table 4.  Table 5.
Case f, s is a variable. The invariants are S = 100 3 10 4 mm 2 , L = 500 mm, h = 1.1 mm, R = 85.4 mm, and D = 111.6690 mm 22 . The values of variable s and the simulation results are shown in Table 6.

Results and discussion
The theoretical values of number of asperity contacts for each simulation case in section ''Simulation experiments'' were calculated by analytical expression equation (16). With the comparison between the theoretical and simulation results, the rationality of the analytical model of asperity contacts can be validated.

Comparison between the theoretical values and the simulation values of case a
As shown in Figure 11, the curve of theoretical values is very close to that of simulation values with the relative deviation of 3.6%-14.6%, and they both indicate that N is in the linear increasing relationship with S. The consistent results indicate that equation (16) can reasonably reflect the relationship between N and S.  Comparison between the theoretical values and the simulation values of case b As shown in Figure 12, under different sliding distances, the theoretical values are very close to the simulation values, and the relative deviation is 6.3%-11.7%. The theoretical and simulation curves show that N is in the linear increasing relationship with L. The consistent results indicate that equation (16) can reasonably reflect the relationship between N and L.

Comparison between the theoretical values and the simulation values of case c
As shown in Figure 13, the theoretical curve and the simulation curve show some exponential relationship between N and h. When the spacing h increases more than a certain value (the value in Figure 13 is 1.1, which is approximately equal to the standard deviation s#), the theoretical values are very close to the simulation values, and the minimum relative deviation is 2.9%. When the spacing is less than this value, the deviation between the theoretical values and the simulation values is increased with the decreased spacing, and the maximum relative deviation is 82.7%. The large deviation is caused by the fact that the asperity height during the simulation calculation follows the normal distribution strictly. Contrastively, in the theoretical model, in order to obtain analytical expression equation (16), the exponential function equation (14) was used to express the distribution of asperity height during the integration. As shown in Figure 5, the normal distribution curve tends to be the exponential distribution when the spacing is larger than a certain value h; when the spacing is less than a certain value h, the normal distribution curve deviates from the exponential distribution. Thus, the deviation between the theoretical value and the simulation value is large in the small value range of h. Actually, when the two surfaces are in elastic contact, the spacing h between the two surfaces is approximately equal to the standard deviation s#, 15 and it cannot obtain those small values as shown in Figure 13. It is indicated that the theoretical model can be used to calculate the dynamic contacts of two rough planes under the elastic contact.

Comparison between the theoretical values and the simulation values of case d
As shown in Figure 14, under the condition of different asperity radii, the theoretical values of the number of asperity contacts are very close to the simulation values, and the relative deviation is 5.3%-24.7%. Equation (16) indicates that the number of asperity contacts N is in the linear increasing relationship with the composite   radius R#. R# can be calculated with Equation (2) and it is in the linear relationship with R. So N is in the linear increasing relationship with R in Equation (16). The simulation results shown in Figure 14 also indicate that N is in the linear increasing relationship with R. The consistent results indicate that equation (16) can reasonably reflect the relationship between N and R.
Comparison between the theoretical values and the simulation values of case e As shown in Figure 15, under the conditions of different asperity distribution densities D, the theoretical values are close to the simulation values. When the asperity distribution density is small, the theoretical values are very close to the simulation values, and the minimum relative deviation is 17.2%. When the asperity distribution density is large, the deviation between the theoretical values and the simulation values is increased, and the maximum relative deviation is 20.0%. It may be interpreted as follows: when the density is large, the asperity number is also large in the given surface area, and it is prone to produce contact interference (one asperity contacts with multiple asperities at the same time, and this contact state is considered as one contact) during the simulation of dynamic contact. Thus, the accumulated contact number is reduced. On the contrary, in the theoretical model, it is assumed that the contact occurs between two asperities, and the contacts among several asperities at the same time are considered as multiple contacts. Thus, the theoretical value is large. The simulation value curve and the theoretical value curve show the same exponential relationship, indicating that equation (16) can reasonably reflect the relationship between N and D.

Comparison between the theoretical values and the simulation values of case f
According to equation (16), the number of asperity contacts is in the exponential relationship with composite standard deviation s# of asperity height distribution. s# and s show the linear relationship (see equation (1)). Thus, the theoretical values of N are also in the exponential relationship with s, as shown in Figure 16. When s is greater than a certain value (about 0.5), the theoretical value curve is very close to the simulation value curve, indicating that equation (16) can reasonably reflect the relationship between N and s. When s is \0.5, the theoretical values decrease rapidly, and the deviation between the theoretical value and the simulation value is large.
The deviation can be interpreted as follows: According to equation (13), during the theoretical calculation, the theoretical value is cumulatively integrated with the asperity height distribution function c(z) within the range (h ; +N). As shown in Figure 17, with the decreasing s, the normal distribution peak   rapidly increases and the integration of c(z) within the range (h ; +N) rapidly decreases. Thus, the theoretical values calculated using equation (13) rapidly decrease. For simulation calculation, as discussed in section ''Mathematical descriptions of nominal rough planes,'' the asperity height follows the normal distribution. As shown in Figure 17, when s is decreased, the peak of the height distribution curve becomes high, indicating that the height of more asperities is concentrated to m and the probability that 2Z . h (see equation (19), which is a condition of asperity contacts) becomes small. It means that the probability of asperity contact becomes small. Thus, in the small value range of s, the simulation values also decrease with the theoretical values. But as can be seen in Figure 17, the difference in the concentration degree of asperity height distribution c(z) decreases with the decrease in s, which means the influence of the change in c(z) on the simulation results is decreased with the decrease in s. Thus, the simulation results obtained under the low s will not show large variation as the theoretical results.
For the commonly machined rough surfaces, the standard deviation s of the height distribution of asperities is not small. According to the measurement results of the metal surface machined with precise grinder, the asperity height is 0.4314 mm, and the standard deviation is 0.5106 mm. The surface machined with other conventional machining methods shows the higher standard deviation. Moreover, the spacing h is decreased with the standard deviation s of asperity height when the two surfaces are in elastic contact. According to Figure 17, when s and h are decreased simultaneously, the integration of c(z) in the range (h ; +N) in equation (13) will not change greatly. Therefore, the theoretical results will not show the large deviation as presented in Figure 16 for actual rough surfaces with elastic contact.

Conclusion
1. Based on the statistical contact theory, a dynamic contact model of rough planes under elastic contact state was established to calculate the total number of asperity contacts for engineering rough planes during sliding process. 2. A computer simulation model of the sliding contact between the rough planes was established to rapidly calculate the number of dynamic contacts during the sliding process. 3. Through the comparison between the theoretical values and simulation values of the dynamic contacts under different parameter conditions, it was verified that the modeling of the dynamic contact of asperities was reasonable, and the theoretical model correctly described the mathematical relationship between the number of asperity contacts and the major influencing factors. The dynamic contact model can provide a reference for the analysis of dynamic contact performance and microscopic fatigue damage of mechanical surfaces.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work received the financial support from the National Natural Science Foundation of China under grant no. 51275083 and the Fundamental Research Funds for the Central Universities under grant N152303011.