Abstract
The contact stiffness of the machined surface has an important effect on the performance of the complex mechanical product. A modified fractal model based on oblique asperity contact is proposed in this research. First, the contact radius and the critical contact area are analyzed based on oblique contact condition. The normal contact stiffness and elastic–plastic force are calculated. The ratio of the actual contact area and a new parameter related to the current contact angle are introduced. Second, numerical simulations indicate the difference. The results show that the stiffness of the oblique contact is smaller, and with the increment of the fractal dimension, the extent of the stiffness reduction is larger. In contrast, the uniform distribution has the lower proportion of the elastic force in the total normal contact force. Finally, experiments including specimen surface observation and load-deformation measurement are utilized to obtain contact stiffness of the machined surface. To some extent, the modified fractal stiffness model is more reasonable and accurate from the result.
Introduction
The contact stiffness of joint surface in a complex mechanical product is a vital component of overall stiffness.1,2 At the same time, it shows that contact stiffness varies with contact load, contact state, material, processing method, surface roughness, and surface structure. These factors result in complicated mechanism and many nonlinear variables when studying the contact stiffness.3,4
Establishing the appropriate contact stiffness model of the machined interface from microscopic scale is a main theoretical method. Hertz’s theory 5 focuses on the small elastic deformation of isotropic material by the normal load. When the deformation is small and the curved surface is smooth from the macro perspective, Hertz’s theory is convenient and practical.
Based on statistical information and parameters of surface morphology, the statistical contact model was proposed and continuously improved. The Greenwood–Williamson model (GW) 6 is a typical statistical model. The asperity summits are supposed to have the same curvature radius. And the asperities are separated from each other without interaction. Their contour heights obey the Gaussian distribution. A discrete version of the GW model is put forward by using an approximate solution based on the geometric-mean summit curvature rather than the statistics of the asperities distribution. 7 Bush, Gibson, and Thomas (BGT) 8 modeled the asperities as paraboloids with two different radii of curvature. Greenwood and Tripp (GT) 9 pointed out that the shape and the relative position of the asperities have no effect on the model. Whitehouse and Archard (WA) 10 studied the correlation and the joint probability density between the heights and the curvatures of the summits. Chang, Etsion and Bogy (CEB) 11 proposed the concept of critical contact deformation and described the statistical model from GW model.
Persson’s theory12,13 generates another important contact models which do not rely on the notion of asperities. These models use a diffusion equation for the probability density of the contact pressure of full contact and extend to finite pressures and partial contacts by imposing a boundary condition. Yastrebov et al. 14 studied the evolution of the true contact area on self-affine rough surfaces and compared the numerical results both with statistical models at light pressures and with Persson’s contact model for the entire range of pressures. Putignano et al. 15 investigated the influence of surface parameters as fractal dimensions, mean square slope, and mean square roughness on the relation between the contact area, the load, and the average separation. Based on Persson’s theory, numerical results of the load, contact pressure distribution, and contact area are compared by many researchers.16–18
When establishing the GW model, Greenwood argues that the average curvature radii of asperity summits are related to the instrument resolution . Stanley 19 found that the surface topography has a non-stable random characteristic and pointed out that the height distribution is related to the sample length. It shows that statistical parameters can only reflect surface roughness information related to the measuring instrument and sampling interval. Furthermore, the approximations in the GW model arise some limitations such as the interaction between asperities and the coalescence between adjacent growing contact zones, and to considering the elastic coupling between contact regions, Ciavarella et al. 20 depicted the interacting Hertzian asperities (IHA) by using the Johnson formulas. Aferrante et al. 21 improved the GW model by considering interacting and coalescing Hertzian asperities (ICHA). Many other works also discussed and extrapolated further and definitive improvements of GW model from different points of view.16,22 On the other hand, fractal characterization of surface roughness is independent of the instrument resolution. 19
Majumdar and Bhushan23,41 (MB) first describe the rough surface profile curve with Weierstrass–Mandelbrot function (WM function) and propose fractal contact model. After that, the fractal contact model is developed continuously. Wang and Komvopoulos24,25 (WK) propose truncated area and the domain extension factor. The formula of the curvature radius in the MB model is improved and the contact model is proposed. Yan and Komvopoulos 26 (YK) use a two-variable WM function to establish a three-dimensional fractal model. Although WM function and its modified versions have been widely used to generate the fractal surfaces in many simulations, they have some limitations. The scale-invariant characteristic of WM function only covers a finite range of length scales.27,28 Its power spectral density has limited bandwidth, and frequencies outside the range determined by the lower and upper wavelengths do not contribute to the generated rough surface. Furthermore, the self-similarity behavior is approximate, and the self-affine property is asymptotical. 29 These limitations imply that WM function does not correctly reproduce rough surface topography, and the features are close resemblance to that of the real surface. However, it is suitable to constructing surfaces with the correct fractal parameters when the instrument resolution and the length of the profile correspond to the values in WM function.
According to these classical and basic theories, scholars have modified, extended, and applied the fractal models in many fields. Chen et al. 30 discussed the effect of the friction coefficient and roughness amplitude. They found that the relationship between load and normal contact stiffness is direct or inverse ratio decided by the value of fractal dimension. Miao and Huang 31 extended the Morag–Etsion (ME) asperity contact model into a complete contact model of a fractal surface and discussed the impact of the fractal parameters. Wan et al. 32 proposed a bi-fractal method and surface impact coefficient to solve the weak points of GW and MB models. The contact body is approximated to a number of asperities stacked on a base, similar to a series structure of two springs. Gong et al. 33 developed a coupling fractal model considering adhesive wear with three-body abrasive wear, and the fractal parameters are time-varying.
Most of these papers follow the hypothesis in the classic references that the two rough interfaces are simplified to a plane and a rough interface. Nevertheless, oblique contact is more similar to the actual condition. For the oblique asperity contact interaction, contact force has an oblique angle at the contact point of asperities.
Jackson et al. 34 researched the elasto-plastic hemispherical force between spheres based on Kogut–Etsion (KE) model and characterized the component of friction. Sepehri and Farhang 35 extended the force between paraboloidal asperities defined in GT model and proposed oblique asperity contact model. Zhao et al. 36 reported a stiffness model in view of shoulder–shoulder contact with Jackson–Green (JG) model. Furthermore, they depicted the influence of the oblique asperity stiffness by material properties and surface state. 37 The sliding interaction between asperities was also introduced considering substrate deformation and power-law hardening materials. 38
According to the experiments and conclusions of these publications, a modified WK fractal model is proposed based on shoulder-shoulder asperity contact. The numerical simulations are adopted to observe the difference between the modified model and the original model. Gaussian distribution and uniform distribution are used as the assumption of the oblique contact distribution. Furthermore, experiments including specimen surface observation and load-deformation measurement are made.
Modified fractal model
Normal stiffness of single asperity pair
In this article, basic model assumptions and approximations should be made:
The asperity of the surface is considered as elliptic paraboloid because the elliptical contact patches are closer to the truth than the circular one. However, the WK model is derived from the contact of two hemispheres. To solve this contradictory, many pioneering researchers have tried diverse methods and obtained satisfying results.7,39 Greenwood
7
introduced a classical method called “mildly elliptical contacts” and applied to the full elliptical solution of BGT model. By replacing the arithmetic mean of two principal radii with the geometric mean, the equivalent radius R is more plausible and BGT model is reproduced more conveniently. This means that we can acquire a good approximation if equation
Oblique asperity contact is usually considered and calculated with the effect of adhesion, sliding, and the friction.19,34,37 Our study merely focused on the normal elasto-plastic behavior under frictionless non-adhesive condition and there is no relative increment of tangential displacement. So the interfacial shear strength, tangential friction force, and elastic energy are ignored. The deformation of substrate which the asperities attach to is also ignored.
The long-range elastic coupling between the contact regions is a very strong assumption for traditional asperity models, for example, GW model or BGT model. Some theories and simulations have shown that this approximation almost always leads to qualitatively wrong results.16,20,39 However, fractal theory has the characteristic of multi-scale like Persson’s theory. The size distribution density of fractal models holds (approximately) for all magnifications of both the micro asperities and the macro ones. So the lateral interaction between the contact asperities is not considered.
The rough surface is isotropic and the material is homogeneity.
According to the assumptions of WK model, the asperities are treated as hemisphere and the mating interfaces can be described as spherical asperities condensed by a rigid plane (Figure 1). R is the equivalent asperity radius.

Single asperity pair and its simplification.
The WK fractal model for the normal pressure p and the normal stiffness
where
The superscript of
where
Figure 2 shows the simplified contact between surfaces. The asperities in Figure 2(a) are compressed by rigid plane as described in WK model. However, the asperities of two real surfaces are deformed more complicated as shown in Figure 2(b). And the amplified area implies that the oblique asperity contact should be considered in different scales. If the asperity has nonzero contact angle
The subscript a means the component along the normal direction of the contact point, while the subscripts n and t mean the normal and tangential components of the mean plane. Some expressions are obtained by geometric relationship 35
Because two asperities are not simplified in Aferrante et al.
35
and
Therefore, the normal load
The normal stiffness

Simplified contact between (a) rough surfaces and (b) two rough surfaces with oblique contact asperities.

Oblique asperity contact.
Multi-asperity normal fractal model
The size distribution density expression
where
The contact area ratio is
When
When
where
Similarly, the dimensionless total normal plastic pressure
where
Due to the difference in curvature radius and offset of the contact point, the contact angle is not the same. The contact angle set can be defined as
As
The sum of the asperities of all contact angles is the total load of the joint surface
The dimensionless form is
Similarly, when the contact angle is
where
Using equation
Accordingly, the total normal stiffness can be written as
The dimensionless form is
Numerical simulations
Simulation preparation
According to equations (22) to (24), (28), and (34), once dimensionless total normal pressure
Domain extension factor
In this work, C45E4 (ISO) steel is chosen as the simulation and experiment sample material. Table 2 shows some relevant material properties. Simulations are given for 9 mm × 9 mm isotropic round surface.
Material properties.
The oblique contact distributions are assumed as Gaussian distribution and uniform distribution. The values of parameter
Comparison of dimensionless normal contact stiffness
Figure 4 presents the dimensionless normal contact stiffness and dimensionless normal contact pressure. The stiffness of the oblique contact is smaller. Furthermore, with the increase of the fractal dimension, the extent of the stiffness reduction is larger. Meanwhile, the stiffness of the proposed model with uniform distribution (m-WK-uniform) is smaller compared with the assumed Gaussian distribution (m-WK-Gaussian). Because of the existence of shoulder–shoulder contact among asperity contact pairs, the normal load on the joint surface is converted into tangential and normal component, and bigger contact angles have bigger tangential pressure. The tangential components of all asperities offset each other through the common base plane, which leads to the reduction of actual load and bearing capacity of the joint surface. On the other hand, the increase of D indicates the increment of surface dimension, which means the irregularity of the height scale increases. This causes higher probability of shoulder–shoulder contact, so the stiffness reduction is more significant.

Comparison of normal contact stiffness of different models.
Comparison of the elastic pressure ratio
The purpose to analyze the elastic pressure ratio is to study the component of the joint surface stiffness when loading. As illustrated in Figure 5, the ratio of the elastic pressure increases quickly and stabilizes at a certain value with the increase of dimensionless normal contact stiffness. We can find out several results from the simulation results. First, the elastic pressure ratio decreases with the increment of fractal dimension. Because of the increment of the surface irregularity, the load on asperities becomes more uneven and some asperities have to bear concentrated force, which lead to the plastic deformation. Second, the modified models have smaller elastic pressure ratio and that with uniform distribution has the smallest. The reason is that when

Comparison of elastic pressure ratio of different models.
Experiment validation
Samples preparation
The micro-scale surface of manufacturing process and finishing technologies has self-affine and self-similar profiles. The cylindrical samples of C45E4 (ISO) steel with milling and grinding surfaces are made (Figure 6). When two samples are compressed as shown in Figure 10, the total displacement of the samples are caused by the cylindrical base and the joint surface. In order to catch the deformation of the contact surface, the single equivalent sample with double length is prepared to subtract the base displacement. Their sizes are Φ18 mm × 40 mm and Φ18 mm × 40 mm, respectively.

Samples of C45E4 steel.
Fractal parameters computation
Fractal dimension and fractal roughness are core parameters for fractal theory. There are many methods to measure these parameters such as the box counting method and the structure function method.
40
Here, the structure function method is the chosen one. The point data of the contact surfaces were measured by three-coordinate measuring machine (ACCURA II AKTIV; Carl Zeiss AG, Oberkochen, Germany) illustrated in Figure 7 and the sensor precision is 3 µm. Each contact surface is divided into four test areas and tested severally according to the sequence in the figure. In this way, fractal parameters can be obtained by the average of the values of all areas to reduce the errors. The mating test areas were marked as couples, for example, Mu1 and Mv1. The point spacing

Measuring field for point cloud data of the contact surface.

Three-dimensional plot of two surfaces: (a) Mu1 and (b) Gu1.
The function of the structure function method is
where
where

The results of structure function method for calculating fractal parameters: (a) Mu1 and (b) Gu1.
Contact stiffness measurement
The measurements of deformation and the force of the samples are tested by microcomputer control electron universal testing machine (WDW-100; Ji’nan Meitesi Co., Jinan, China). The linear degree and sensitivity of the deformation sensor are 0.12 (%FS) and 2% mV/V. The single equivalent sample testing experiment equipment is shown in Figure 10, while the short samples are condensed together. The force sensor is above the clamp. The normal contact stiffness in unit area
where

The single equivalent sample testing experiment equipment.

Normal deformations and pressures of different testing samples.

Normal deformations and pressures of the joint surfaces: (a) milling surface and (b) grinding surface.
Results and discussion
The simulations of WK fractal model and the modified models with different size distributions are compared with experiment results (Figure 13).

Comparison of dimensionless normal stiffness by experiment results and fractal models: (a) milling surface and (b) grinding surface.
The dimensionless total normal contact stiffness is calculated by

Errors of dimensionless normal stiffness between experiment results and fractal models: (a) milling surface and(b) grinding surface.
Combined with the analysis of the simulations before, the reason is that there are many plastic deformations occur when loading. The modified models can describe the contact state more reasonable for considering the shoulder–shoulder contact. The stiffness of experiment is lower than the modified models at the beginning, then higher after certain point. It indicates that the actual contact angle and size distribution of the shoulder–shoulder contact are more complicated than the assumed distributions. Actually, from the feature of the m-WK-Uniform model and m-WK-Gaussian in the simulations and experiments, it indicates that the assumed uniform distribution is closer to the real surface topography in the low pressure. With the increasing of the pressure, the asperities with small radii become part of the asperities with big radii gradually. Thus, the shoulder–shoulder contact becomes less in large pressure, and the WK model can predict the stiffness accurately relatively.
In the entire process, the contact angle is increasing and the amount of oblique asperities is reducing. Consequently, the stiffness of the modified model deviates from the experiment when the pressure is high enough. However, the modified model is well approximated to the experiment relatively. The key factor to predict accurate stiffness in tribological applications is based on the precision of fractal dimension, actual contact angle, and size distribution of the contact surfaces.
Conclusion
A modified fractal model was raised based on WK model. In this article, oblique contact was utilized and analyzed. The ratio of the actual contact area and a new parameter related to the current contact angle were introduced when calculating normal contact stiffness. Then, numerical simulations showed that the stiffness of the oblique contact was smaller, and with the increment of the fractal dimension D, the extent of the stiffness reduction was larger. Comparing with the Gaussian distribution, the uniform distribution had greater errors, because the normal load on the joint surface is converted into tangential and normal component, and bigger contact angles have bigger tangential pressure. And the increase of D means the irregularity of the height scale increases. This causes higher probability of shoulder–shoulder contact, so the stiffness reduction is more significant. Furthermore, the uniform distribution had lower proportion of the elastic force in the total normal contact force because the coefficient of elastic pressure
Footnotes
Handling Editor: Michal Kuciej
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 was supported by National Natural Science Foundation of China (Grant No. 51775552) and Key Research Project of the Education Department of Hunan Province (Grant No. 15A107). The authors gratefully acknowledge these supports.
