Contact analysis of measured surfaces based on modified interacting and coalescing Hertzian asperities (ICHA) model with oblique contact asperities

The contact mechanism for joint surfaces is important for predicting the loading process and dynamic properties of precision machinery products. A multiasperity model considering the lateral interaction and coalescence of contact regions with oblique asperity contact, the interacting and coalescing Hertzian asperities (ICHA) model, is introduced in this work. Contact angle α and radius of curvature R ( r ) are vital parameters that are directly calculated based on grid points measured from the milling and grinding surfaces. Numerical simulations of α and R ( r ) were analysed. The results show that the proposed model is more effective at improving the contact area versus force relationship than various asperity models for low and middle contact forces. α and R ( r ) follow a Gaussian distribution. Because of the existence of coalescence, they decrease with increasing contact force. Furthermore, the α and R ( r ) in the startup phase are crucial data that can be used to judge the effects of oblique asperity contact. They are greatly influenced by fractal dimension D and fractal roughness G. From the comparisons with the surfaces generated by Weierstrass-Mandelbrot function, it is suggested that to obtain an accurate contact prediction in the loading process, oblique contact is as important as lateral interaction and coalescence.


Introduction
The existence of joint surfaces between assemblies has enormous implications on the dynamic properties of precision machinery products. 1,2 The mechanism of contact action involves friction, wear, lubrication, adhesion, thermal conductivities and other contact statuses. [3][4][5][6] Many analytical methods have been proposed and improved to model machined surfaces. However, it is difficult to establish a complete, accurate and universal contact model owing to complex statuses in multiphysics. Most theoretical models are derived with some assumptions under given constraints and particular conditions. The classical Hertz theory 7 assumes that the deformation of the contact body is small and the relationship between the stress and strain is linear. The Greenwood-Williamson (GW) model 8 is an asperity-based contact theory that treats the shape of the asperity as spherical or elliptical. The radius of curvature, size distribution and height distribution are defined and measurable. Fractal theory 9 uses the self-similarity and self-affine characteristics of the machined surface and uses the size distribution of contact spots to follow the cumulative size distribution of islands. Wang and Komvopoulos 10 improved the microcontact distribution function with a domain extension factor based on the continuous approximation of the Dirac delta function. Persson's theory, 11 which deduces the surface separation and the contact pressure distribution from full contact to partial contact through boundary conditions, is practical and has grown popular in recent years. Based on these basic contact theories, in this study, many extended and modified versions are proposed, while some arguments are also raised that the conditions needed to obtain the exact relationship between the real contact area and load are harsh and restrictive. The interacting and coalescing Hertzian asperities (ICHA) approach has been shown to be extremely accurate at solving complex contact problems.
Yastrebov et al. 12 introduced a correction function and smoothing technique to deal with the contact boundary of a discrete mesh. The real contact area of a random rough surface can be calculated using this method. Putignano et al. 13 found that the contact area is only related to the root mean square slope of the asperity height distribution. Numerical simulations also show that the relationship between the contact area and the load is nonlinear. Paggi and Ciavarella 14 compared the accuracy of the coefficient of proportionality k between the real contact area and load using various statistical theories, such as the Bush, Gibson and Thomas model and the GW model, using a range of generated fractal surfaces. The efficiency and accuracy can be significantly improved by considering the interaction effects and simplifying the mesh discretisation method. Most related studies on contact mechanics focused on load, true contact area, surface separation, gap distribution, interfacial stress, multi-scale characteristics and related parameters of generated height profiles by Fourier transform. 15 Moreover, few studies have introduced a detailed method that can measure and analyse real machined surfaces.
However, the micro-scale and nano-scale accuracy of precision equipment such as optical, capacitive, tomography and tactile probing systems have been greatly improved and are widely used. In addition, the reliability of data obtained from machined surfaces has been improved. Thus, it is reasonable and feasible to analyse real machined joint surfaces by measuring the surface point cloud. Bora et al. 16 constructed a finite element numerical simulation model using surface data obtained via atomic force microscopy and estimated the relationship between the force and displacement of normal degrees of freedom using the Boussinesq solver. Campan˜a´and Mu¨ser 17 measured the surface morphology of refined steel via atomic force microscopy and calculated the area occupied by each atom using the Voronoi structure. The results showed that the contact pressure of randomly rough surfaces follows a Gaussian distribution, while the experimentally determined topographies obey an exponential distribution. The true contact area of the measured surface is 20% larger than that of the random surface. Putignano et al. 18 regarded the effect between tire and pavement as being a grinding action that has fractal self-similar characteristics. The contour of the concrete pavement was measured by contact-less optical methods and compared with the surface generated by the power spectrum method. The variation of the limit contact area corresponding to the percolation threshold was analysed.
Some previous studies also found pronounced differences between real machined surfaces and generated random rough ones. One important reason is that the two rough surface contacts are usually assumed to be an elastic half-space or a rigid plane in contact with an isotropic rough surface (Figure 1(a)). However, the joint surfaces between real parts in precision machinery products usually have the same finish roughness. The oblique contact between asperities is unavoidable when the scales of the upper and lower rough surfaces are not significantly different (Figure 1(b)). The effect of oblique contact is simplified by various degrees or neglected except in cases involving friction, sliding and adhesion. However, it is clear from the visualisation and principle 19,20 that shoulder-shoulder contact between the asperities of rough surfaces is not only persistent but also exists on scales that differ from the amplified area in shown in Figure 1. Thus, a multi-scale analysis is necessary.
Sepehri and Farhang 19 presented the normal and tangential load of single asperity based on the Hertz theory, in which the asperity height sum is assumed to be a Gaussian distribution. The minimum error of the total normal and tangential loads can reach 0.2%. However, the assumed Gaussian distribution was not verified, and the multi-step hybrid interactive and optimisation procedure was not clear. Jackson et al. 21 described the average tangential and normal forces from the asperity interactions using semi-analytical and finite element simulations. The properties of typical metals were the study subjects to help verify their effectiveness. Mazeran and Beyaoui 22 depicted the sliding of nanoscopic elastic contact under a scanning force microscope Si 3 N 4 probe using an extended Hertz theory (Derjaguin, Muller and Toporov theory). The experimental results were in relatively good agreement with the Mindlin and Savkoor theories. Zhao et al. 23 extended the elastoplastic model considering the surface hardness proposed by Jackson and Green using the oblique contact condition. A white-light interferometer was used to measure the surface profile, and the asperity shape was fitted to an elliptical paraboloid. Because of the quadratic approximation of the radius of curvature and the offset of the intersection midpoint from the peak, the results still had a deviation error.
In this study, multiasperity contact analysis considering lateral interaction and the coalescence of contact regions is extended with oblique asperity contact. A suitable method was applied to the data clouds of two machined surfaces obtained by a three-coordinate measuring machine. As a result, a new method to calculate the contact angle and radius of curvature of the oblique contact point is discussed in detail. Contact stiffness measurement is used to indirectly verify the accuracy and efficiency of the proposed model and simulation method. Experimental data from milling and grinding surfaces were simulated and compared with generated topographies for qualitative analysis. The results show that the oblique contact can be considered to be the effect of interaction and coalescence, especially for low and medium loads on the machined surface. The fractal parameters play important roles in the contact angle and radius of curvature.

Contact model
In this work, frictionless non-adhesive contact was studied without considering sliding, wear and other mechanisms that occur in the tangential direction. Although oblique contact always has a tangential force, the normal force is the key behaviour. The material of the machined surface is homogenous. Asperity-based theory and Persson's theory are the most commonly used surface contact theories. 24 Persson's theory 11 is practical for adhesion, lubrication because the multiscale and long-range elastic coupling properties enable the approximate process to obtain highly accurate simulation results. However, this theory extends the differential diffusion equation under full contact to partial contact without involving asperity topography. Asperity-based theory is a multi-asperity theory and has more advantages for situations of partial contact and low load. Various approaches have been introduced to consider the defects that may lead to significant errors. ICHAs were applied and extended to a newly proposed approach with oblique contact. Figure 2 presents the geometric diagram of a single oblique asperity contact pair. In the contact model, the shape of the asperity is approximately elliptic paraboloid. The function of the contact force, f a , equivalent radius of curvature, R(r) and asperity interference, w, of the contact point can be given by the Hertz theory: where E 0 is the equivalent elastic modulus. According to the geometric relationships and derivative results of Sepehri and Farhang, 19 the normal components can be obtained as w = w n cos a ð2Þ where a is the contact angle, w n is the normal interference of the asperities and f n is the normal component of the contact force. The method 19 proposed to calculate R(r) uses the fourth-order approximation as R(r) = R 1 + r 2 R 2 s À Á 1:5 , and a can be given as cos a = 1 + r 2 R 2 s À Á À0:5 , in which r 1 = R 1 r=R s , Þ =R 1 R 2 , and r 1 ,R 1 and r 2 ,R 2 are the tangential deviations and radii of two summits, respectively. R is the equivalent asperity radius for the contact of two hemispheres obtained by the arithmetic mean of the two principal radii. However, the elliptical contact patches are closer to the truth than the circular ones. Greenwood 20 demonstrated that the geometric mean value can obtain the plausible equivalent asperity radius of an elliptical solution with good approximation, which can be written as: To consider the lateral interaction and coalescence of contact regions, the Interacting and Coalescing Hertzian Asperities (ICHA) multiasperity contact model proposed by Afferrante et al. 25 was applied to calculate the normal displacement, w ni , of a generic point in a non-contact region of asperity ASP i : where a i and a j are the ith and jth contact radii of asperities ASP i and ASP j , R i is the geometric mean of the principal radii of curvature, d ij is the tangential distance between asperities ASP i and ASP j , and n a is the total quantity of asperities. Meanwhile, two contact asperities for two surfaces are defined as a new asperity with equivalent contact radii a eqv , geometrical volume centroid, x eqv , y eqv À Á and radius of curvature, R eqv , which can be given as: Radius of curvature based on grid points One practical problem is that it is difficult to fit all asperities in real rough surfaces using elliptic paraboloids. Even if it works, the errors are great because the profile of the surface is composed of multi-scale asperities. Moreover, the valleys and peaks are determined by the machine tools and scratches. Surface defects are inevitable. Thus, the calculation method for the radius of curvature based on the profile of a real machined surface can reduce such errors and become more reliable. 26 The data of the profile are usually made up of discrete mesh points that coincide at regular intervals with differing normal interfacial separations. Figure 3(a) shows a three-dimensional (3-D) image of regular non-obtuse triangulated asperity ASP i meshes, where v i is the vertex and its 1-ring neighbours v j + k : 0 ł k ł 7 È É . The unit length normal vector N vi of v i (Figure 3(b)) can be obtained as: 27 where F i is the surface list of the triangle faces formed by v i and its 1-neighbourhood points, and f k 1 ł k ł 8 ð Þ is a 3-D triangle face belonging to F i . N f k is the unit length normal vector for each face. The area f k j j is treated as the weighting coefficient introduced by Taubin. kk means norm operation. Vector angle, centroid and triangle geometrical features also affect the normal vector, and many studies have proposed improved algorithms. Because the precision profiler can achieve uniform space sampling and the machined surface has regular textures and patterns, the accuracy of the surface area is sufficient to be the weighting coefficient.
There are several classical methods of discrete curvature estimation based on grid surfaces, such as tensorbased techniques, derivative methods with local quadric approximation and mathematical formulae based on operators. 28 Meyer et al. 29 introduced a flexible tool that approximates arbitrary triangle meshes by averaging Voronoi cells. This method is applied to calculate the equivalent radius of curvature R(r) of the contact point in equation (1) directly without the tangential deviations and fitting process of the asperity shape. However, R(r) is the radius of the mean curvature rather than the radius of the normal curvature, principal curvature or Gaussian curvature. The grey parts in Figure 3(c) depict top view of the Voronoi region for a vertex with v i and the neighbours v j . O j is the circumcentre of face f k . The non-obtuse Voronoi area A Vori, j in triangle Dv i v j v j + 1 and the total Voronoi area A Vori of v i can be given as: where n is the amount of face f j and b i, j , g i, j , u i, j and l i, j are the angles between vectors, as shown in Figure  3(b). In addition, A Vori, j = 0: is an obtuse angle. The radius of the mean curvature R i of the asperity ASP i and equivalent radius of curvature R(r) can be obtained by equation (4) as: Where R 1i and R 2i are the radii of the mean curvature of the two asperities at the contact point, as shown in Figure 2(b).

Contact area and contact angle
Although the contact region of the oblique contact has a slope, the contact area from the projection view is treated as the true contact area. Figure 4(a) describes a 3-D mesh of arbitrary asperity ASP i , in which the grey zone is the contact region. The contact area extracted from the coarse grid can be computed simply by counting the number of grid points N g . Each point contributes an area, Dx 2 , where Dx is the grid side length. The nominal contact area is a i = N g Dx 2 . From the top view in Figure 4(b), a i there are significant errors because the interface contact points are in partial contact rather than full contact. Yastrebov et al. 30 introduced a useful corrective factor to compensate for this deviation by: where S d = MDx is the perimeter of the contact/noncontact boundaries and M is the number of switches from contact to non-contact and vice versa along grid points both in the vertical and horizontal directions, as shown by the red line in Figure 4(b). As mentioned, cos a = 1 + r 2 R 2 s À Á À0:5 can be used to calculate the contact angle. However, for real profiles measured by the machined surface, it is difficult to obtain r and other tangential deviations for the same reason that it is difficult to obtain the radius of curvature R i . Figure 5 is the left view of Figure 4(a), and it presents the deformation of a single asperity with oblique contact. Grid points B x B , y B ð Þ and D x D , y D ð Þ are the crossover points of the contact boundaries and the curve between summit A and contact point C corresponding to the position in Figure 4(a). Thus, the contact angle a can be calculated between the angle of BD and substrate as follows:

Measurement and simulation
Surface measurement and data processing Precision instruments, such as atomic force microscopy (AFM), scanning force microscopy (SFM) and surface force apparatus (SFA) optical, capacitive, tomography and tactile probing systems, are becoming more and more accurate at surface measurement. 31 A threecoordinate measuring machine (ACCURA II AKTIV, Carl Zeiss AG, GER) was used to measure the surface topography, as shown in Figure 6. Four samples formed by milling and grinding forming methods were selected as the upper and lower joint surfaces. The sample material was C45E4 (ISO) steel. To take the effect of the fixture, machine tool and processing strain into consideration, four test areas (1-4) of each sample were collected to compare and optimise the measured data, as shown in Figure 6. The holes were utilised as the benchmark to calibrate the position of the test areas.
The sample size was 20 3 20 mm, and the test area was 5 3 5 mm. The sample interval Dx was 50 mm, and the measurement resolution was 2 mm. It is necessary to process the measured data with random errors, systematic errors and noise. Amongst these, systematic errors can be ignored because the relative positions of the data points remain unaffected. Furthermore, the flatness and parallelism of the upper and lower surfaces were analysed to correct the deflection. For random errors and noise, general data processing methodology is not appropriate. When these errors are reduced, the multi-scale features are also removed. The wavelet-based method is feasible for machined surface denoising. 32 By decomposing and reconstructing the measured data, the noise can be reduced by setting low-pass and high-pass frequency thresholds at different levels. However, because the measurement interval is still large relative to the size of asperity, it is generally difficult to continuously take a plurality of points on an asperity slope. The surface obtained is zigzag, and the asperity is tapered. To compute the parameter of oblique contact for a smooth asperity surface, data interpolation is required. 33 Detailed illustrations of wavelet denoising and data interpolation are not outlined in this work as they can be implemented easily .  Because machined surfaces usually have multi-scale and self-similar profiles 34,35 fractal dimension D and fractal roughness G were used to quantify the surface features. The fractal parameters of each sample are evaluated using the structure function (SF) method, which can obtain a good result for two-dimensional fractal dimensions D f D = D f + 1 À Á and G compared to methods of the same type. 36

Surface profiles generated by WM function
Regenerated surfaces are used to evaluate the relationship between the fractal parameters and the characteristics of the oblique contact. Multi-scale and fractal profiles can be synthesised using many techniques, such as the Fourier transform filtering algorithm with prescribed Hurst exponent or power spectrum density, the Weierstrass-Mandelbrot (WM) function method and the midpoint displacement method. 37,38 The surface height function z(x, y) of the 3-D fractal surfaces generated by the WM function method is composed of simultaneous sinusoidal equations. A modified WM function can be given as where x and y are the planar Cartesian coordinates, L is the sample length, g is a scale parameter that is, usually set at 1.5, M is the number of superposed ridges, f is a random phase and n (n max is the upper limit) is a frequency index. Figures 7(b) and 8(b) show the surface profiles generated by the WM function with the same fractal parameters as those of the milling and grinding  surfaces. The asperities of real machined surfaces are smoother, while the shapes of asperity are similar to those of the cone for the generated surfaces. Furthermore, it has been pointed out that the WM function has some limitations. 35,39 For instance, the approximate continuous power spectrum contains limited bandwidth, which means that the function can generate limited characteristics of specific scales. Some studies also found that the regenerated profiles based on the self-affine WM function are not consistent with the machined real surface profiles. 35 However, this function can reflect the variation trend between the fractal parameters and contact mechanics to some extent. Thus, the WM function is sufficient for reliable qualitative analyses.
Procedure for numerical simulation Figure 9 shows the contact procedure of two surfaces. The x and y coordinates of the data points in lower surface S l and upper surface S u are initially fixed with z coordinates z u .z l (z u 2 S u , z l 2 S l ). The contact iterative operation is the downward displacement of z u with a fixed z l . When z u .z l , the asperity pair begins to contact. Before the iteration, some data need to be calculated, including the coordinates and radii of curvature R(r) of the initial contact points of all asperity contact pairs, PA i (lASP i , uASP ui ) : lASP i 2 S l , uASP i 2 S u , 0\i ł n a f g , using equations (13) and (14). After that, the surface S u starts to descend step by step as z u(next) = z u À Dz according to the height gap, Dz which can be arranged by the values of z u . The overall procedures of each step to simulate contact asperity pair PA i of the upper and lower surface profiles are detailed as follows.
(1) a i and a of PA i can be calculated using equations (15) and (16) if w n .0 z u \z l ð Þ. Otherwise, consider the next asperity pair, PA i + 1 .
(2) Determine contact areas a i for all asperity pairs PA j using equation (15). if a j .0, asperity pair PA j contributes to w n . (3) (a) If j = i, it means that PA i is the first contact pair on the surface without lateral interaction, and the contribution to w n is a 2 i R i . (b) If j 6 ¼ i and there is no overlap between PA i and PA j , PA j contributes to w n is given by the term in the brackets in equation (5). (c) If j 6 ¼ i and there is overlapping between PA i and PA j , it means that lASP i and lASP j coalesce in contact spots. Thus, they are replaced by a new asperity with a eqv , x eqv , y eqv À Á , and R eqv by equations (6)- (9). Then, return to Procedure 1. Through these procedures, the total normal force F n and total real contact area A n for each step can be obtained by F n = P n a i = 1 f n and A n = P n a i = 1 a i . Moreover, the distribution of contact angle n a ð Þ, distribution of contact radius n R ð Þ, total real contact area fraction A n =A a and distribution of contact clusters can be determined, where A a = 5 3 5 mm is the apparent area of the contact region. The dimensionless total normal contact force F Ã n can be given by Results and discussion

Contact area
To demonstrate the accuracy of the contact area evaluation by the proposed model, which is the modified version of ICHA (m-ICHA), the total real contact area fraction A n =A a and dimensionless total normal contact force F Ã n relative to different contact models are plotted in Figure (10) by averaging the four pairs of different test areas in milling and grinding surfaces. The model (IHA), which considers elastic coupling, was introduced by Ciavarella et al. 40 Some studies have already found that the ICHA model is the most accurate model. From results relative to GW, IHA, ICHA and m-ICHA models, it also reveals that the ICHA model can obtain good results and it is critically essential to lateral interaction and coalescence of contact regions. Otherwise, the area versus load relation would be predicted with poor accuracy, as shown for the GW and IHA models. Furthermore, the results of the m-ICHA model are not in perfect agreement with those of the ICHA model, especially when F Ã n \0:3. The results of the milling surface are clearer. Evidently, the effect of oblique asperity contact leads to bias. Equations (1)-(3) yield: f n = 4( cos a) 0:5 E 0 R(r) 0:5 w 1: 5 3, ð19Þ which shows that the contact angle a and radius of curvature R(r) are vital influence coefficients. Generally, a ł 1, and the R(r) of the slope is larger than that of the summits of the same asperity pair. Under the comprehensive effect of both coefficients, the bias becomes smaller and almost equals zero with increasing contact force. This exhibits that the contact model for oblique contact can lead to more accurate results. Figure 11 shows the joint distribution of contact angle a and real contact area fraction A n =A a for milling and grinding surfaces, and the insets the amplified regions of A n =A a & ½0, 0:3. Because of the existence of coalescence in equations (6)-(8), contact spots can merge forming contact patches. The contact angle a becomes smaller (the amount of contact pairs with big a decreases relatively) with the increase in contact area. For arbitrary A n =A a , a presents a Gaussian-like distribution, and its mean value and variance decrease rapidly until A n =A a ' 0:6. Moreover, the a of the grinding surface is almost equal to zero, corresponding to the overlapping of the ICHA and m-ICHA models in Figure 10(b). This suggests that the oblique contact asperity pair almost vanishes when the contact area is large enough. However, coalescence plays a very important role in the contact mechanism. In particular, for low and middle contact forces, a large number of micro-asperities that exist in a macro-asperity region gradually become part of the macro-asperity. Finally, the dominant macro-asperities are almost turned into the base of the surface. Seen in another light, this explains that the oblique contact takes effect mainly in the low and middle contact forces.

Radius of curvature
The results are presented for the joint distribution of the radius of curvature R(r) and real contact area fraction A n =A a , as shown in Figure 12. Equation (9) indicates that R(r) also increase with an increase in the contact area. The distribution of R(r) is also Gaussianlike for an arbitrary A n =A a . Its mean value undergoes a slight change in the logarithmic coordinate, while the variance of the trend becomes much smaller in processes. The distribution width of the milling and grinding surfaces is almost the same, while the distribution of the grinding surface is more concentrated. This is because the fractal dimension D on the grinding surface is larger, which means that the ratio of micro-asperities is larger. Although the fractal roughness G of the grinding surface is larger, its effect is much slighter. Some approaches 35 also observed that fractal roughness G depends on fractal dimension D, so it is meaningless to discuss G without D.
The comprehensive analysis of Figures 10 to 12 indicates that lateral interaction and the coalescence of contact regions can lead to definitive improvement of the contact model. The oblique asperity contact can also further enhance the prediction accuracy, especially for low and middle contact forces. Under practical circumstances, it is difficult to reach a value of one for real contact area fraction A n =A a for the machined joint surfaces. Hence, taking full consideration of the oblique Figure 9. The contact procedure of two surfaces. asperity contact with contact angle a and radius of curvature R(r) can optimise the loading process.

Initial contact angle and radius of curvature
From the above conclusions, it is clear that contact angle a and radius of curvature R(r) in the startup phase seem to be the crucial points that can judge the effect of oblique asperity contact. Since the A n =A a ,a and R(r) of the grinding surface tend to stabilise more quickly, as shown in Figures 10 to 12, it can be deduced that the deviation ratio of the mean value has a direct relationship with the variability rate of the contact area. Fractal dimension D and fractal roughness G represent the multi-scale characteristics of machined rough surfaces. Hence, the profiles generated by the WM function were simulated to explore the comparison with different fractal parameters.
The value range of the fractal dimension D is defined as 2\D\3. However, previous studies have shown that the D of real machined surfaces is only possible in a particular range. For instance, Persson 41 found out that D is usually \2:3 for polished or sandblasted surfaces by the power spectral density (PSD) method. Zhang et al. 36 demonstrated that such results occur when nominally self-similar real rough surfaces are assumed to be self-affined. The PSD method may not provide an adequate approach for real surfaces, while the SF method can be regarded as a reliable one of various methods. In the present work, D : 2:4, 2:5, 2:6 f gis set as the simulation variable according to the fractal parameters of the measured milling and grinding surfaces. Figure 13 presents the probability density as a function of contact angle a for G = 3 3 10 25 mm ( Figure  13(a)) and G = 6.5 3 10 24 mm (Figure 13(b)). The results were obtained by averaging on eight different measured and numerically generated surfaces. The proposed comparison makes clear that fractal parameters have a tremendous influence on the initial a. The  probability density distribution of generated surfaces is well described by Gaussian distribution. Measured surfaces do not strictly obey Gaussian distribution and are not in good agreement with the generated surfaces for the same fractal parameters, especially for a.p=4. Close observation reveals that the mean values of measured surfaces are smaller. Furthermore, the distribution for smaller D and bigger G is much closer to zero with a smaller variance.
In addition, the contact angle a of the oblique contact for elliptic parabolic asperities should be 0 ł a\p=4, 42 but it was a.p=4 in the numerical simulations, as presented in Figure 13. This contradiction is caused by the theoretical error that the rough surfaces are fitted by discrete grid points at regular intervals Dx. The grids in Figure 14 depict different length scale asperity pairs obtained from the same surface. The small length-scale asperity corresponds to Dx, which is determined by the resolution of the instrument and the sampling interval of the simulation. In this case, the transition curve is replaced with a straight line, and the shape of the elliptic parabolic asperity is turned into a cone. Hence, somewhere in the surface occurs as a contact pair with a.p=4. Irrespective of the small Dx, the asperities shown in Figure 14(b) are always existential owing to the multiscale characteristics of the machined surface. Moreover, almost all initial contact angles a of the black line in Figure 14(a) are .p=4. This indicates that most asperity scales on the generated surface are almost approximate with the interval in the startup phase when D = 2:6 and G = 3 3 10 25 mm. Thus, it can be concluded that there are some limitations in analysing the proposed contact model based on measured grid points: the accuracy of the simulation is constrained by the resolution of the instrument at very low contact forces. This is one of the main criticisms of asperity-based contact theory.   Figure 15 shows the probability density for the initial radius of curvature R(r). Similarly, the distribution of probability density versus the R(r) relationship follows the Gaussian distribution, and a larger D results in a smaller R(r). Compared with contact angle a in Figure  13, the differences amongst R(r) of various fractal parameters are smaller. Actually, a and R(r) are not completely independent. When a converges towards zero, the contact point is near the summit and R(r) is noticeably close to the radius of curvature for the summit, which is the smallest in the asperity. However, this dependency is strong only for a small a. As shown in Figure 16, the joint distribution of R(r) and a corresponds to their distribution in Figures 13 and 15, respectively. And R(r) and a appear to be mutually independent.

Contact stiffness experimental verification
Direct measurement of contact area, contact angle and radius of curvature is difficult. Therefore, contact force and contact stiffness of the samples are taken as the indirect variates to verity the theory.
The data of the contact stiffness experiment in Zhang et al. 42  Based on the procedure for m-ICHA numerical simulation proposed in the paper, we can obtain accurate A i =A n as listed in Tables 1 and 2. Figures 17 and 18 show the dimensionless normal stiffness and the errors of experiment, simulations of the modified WK models with different size distributions (Gaussian distribution, uniform distribution and the simulated distribution).
Combined with the discussion in Zhang et al. 42 the dimensionless normal pressure and the dimensionless  normal stiffness of milling surface and grinding surface are calculated based on the simulated distribution of Real contact area fraction A i =A n . It can be seen from the figures that in the entire process, the stiffness is almost in line with the actual experiment results and the variance is small. The maximum errors of Gaussian distribution and uniform distribution reach more than 10% while the errors of the simulated distribution are less than 5%. Especially during the light load process between the milling surfaces, the accuracy is improved by more than two times. Moreover, the both errors for milling and grinding surfaces are less than 5% which shows that the simulation can reflect the actual contact principles of machined rough surfaces.
Thus, it can be summarized that the modified WK model considering the oblique contact can achieve accurate prediction of the dimensionless normal stiffness when the distributions between contact force, contact angle a and real contact area fraction A i =A n are estimated precisely. And it also verified the correctness of the proposed model in the paper.

Conclusion
We analysed the contact characteristics of real machined rough surfaces by extending a multiasperity model considering lateral interaction and the coalescence of contact regions with oblique asperity contact. Numerical simulations were applied with directly calculated contact angles a and radii of curvature R(r) based on grid points. The prediction abilities of the proposed method were assessed. Compared with several contact models, the contact area evaluation of the suggested model can lead to more accurate results, especially at low loads. The contact angle a and radius of curvature R(r) are vital contact parameters that are affected by coalescence in the loading process. Both follow a Gaussian distribution and decrease with an increase in the contact force. Furthermore, the rate of the reduction is influenced by a and R(r) in the initial state. The effect of oblique contact is less with a smaller fractal dimension D and larger fractal roughness G, as seen from the comparisons with the measured surfaces and generated surfaces by the WM function. However, the joint surfaces of precision machinery parts usually have relatively large D and small G values. Based on the contact stiffness experiment, the accuracy and efficiency of the proposed model and simulation procedure are verified. Thus, to accurately predict the loading process, it is essential to take oblique contact, lateral interaction and coalescence into consideration.