A computational method for equivalent area distribution characterized by the normal contact behavior of two rough surfaces

Abstract

The contact state of surfaces significantly influences the dynamic properties of assemblies, where the contact area is the key parameter to describe the contact state of assemblies. To address the computational challenges associated with contact area distribution, an improved nominal contact area distribution model is proposed for accurately characterizing the contact state of surfaces. Based on the distribution pattern of contact area, a cumulative quantity function for nominal contact area in logarithmic form is established. And its superiority is proved by comparing with two traditional functions. Based on this, a fractal dimension conversion formula was derived to equate the contact of two rough surfaces to the contact of single rough surface, thereby the optimized expression for the nominal contact area distribution function is obtained. The proposed model is verified through actual contact analysis of the anchorage system (anchor plate and anchor base plate). The results show that the proposed cumulative quantity function for nominal contact area and equivalent fractal dimension function exhibit excellent fitting accuracy and feasibility. In the contact between two rough surfaces, the surface with the smaller fractal dimension dominates the contact behavior. The effectiveness and applicability of the proposed model in complex contact problems are verified through case studies. This research provides new theoretical tools for multiscale mechanical analysis of complex contact surfaces, which has important application value in these fields such as precision manufacturing and tribological optimization.

Graphical abstract

Keywords

fractal theory / rough surface / pressure grade / equivalent fractal dimension / normal contact area distribution function

1 Introduction

The contact mechanics analysis of rough surfaces is key in the fields of friction, wear and lubrication, electrical and thermal resistance [1,2]. Early studies were mainly based on statistical methods, such as Greenwood and Williamson (G-W model) assuming that the surface asperity height obeys a Gaussian distribution and using Hertz contact theory to describe the single-point contact behavior [3]. However, the parameters of the statistical model (roughness, slope, curvature) are dependent on the resolution of the measuring instrument and the sample size, which cannot characterize the multiscale fractal properties of the surface. In 1991, Majumdar and Bhushan (M−B model) [4] first introduced the fractal theory and proposed a contact model based on fractal parameters by using the Weierstrass−Mandelbrot (W−M) function to describe the surface morphology. Nevertheless, it assumes that the contact points obey the power-law distribution rule and does not consider the multi-stage characteristics of elastic-plastic deformation.

With further research, Yan and Komvopoulos [5] extended the W−M function from two- to three-dimensional to derive the solution of elastic−plastic contact force and contact area, which revealed the influence of surface morphology parameters on the contact characteristics. To further solve the problem that the size distribution of contact points in the M−B model does not match with the reality, Miao and Huang [6] proposed a complete fractal contact model. According to the modified asperity model, the critical area of a single asperity was scale dependent and that the asperity's plastic to elastic mode transition agreed with classical contact mechanics. Liou et al. [7,8] modified the elastic-plastic transition condition by combining the non-Gaussian probability density function and the classical contact mechanics theory for spherical and cylindrical fractal surfaces, and found that the fractal dimension D and the scale coefficient G significantly affect the contact stiffness.

Even now, the simplified method for inter-contact (which equates the contact between two rough surfaces to the contact between a rough surface and a rigid plane) proposed by Majumdar and Bhushan [4] continues to be utilized in the framework of the fractal contact model. The contact area density function used in these fractal contact models is constructed based on the scale independence of fractal geometry, commonly expressed as n(a) = (D–1)a_L^(D–1)/2a^(1–D)/2/2. This area distribution function is widely employed in various contact models, such as contact models considering scale effects [9–11], strain hardening [12], or interactions [13,14]; contact models under high-speed rotation [15] or cyclic vibration [16]; and contact models involving loading and unloading cycles [17–20], among others. With further research, the domain expansion factor ψ serves as a key parameter for modifying the contact area distribution range in the fractal contact model, which makes the area distribution function better suit the actual contact characteristics of rough surfaces. For specific geometric surfaces such as curved surfaces [21], spherical surfaces [22,23], and cylindrical surfaces [24], it is now common practice to incorporate ψ into the contact model to calculate the contact area and force at each deformation stage. Additionally, the area distribution function involving ψ has been applied in contact models considering scale effects [25,26] or strain hardening [27], contact models under loading and unloading cycles [28], and contact models for bolted connections [29–33]. In summary, rough surface contact models based on fractal theory have been gradually developed. By introducing elastic-plastic deformation, self-affine properties, material anisotropy and other factors, and combining numerical simulation and test validation to continuously optimize the distribution function of the contact points, a more accurate and widely applicable framework of fractal contact mechanics is finally formed.

The current studies of surface contact are based on simplifying the contact of two rough surfaces to the contact of a single rough surface and an ideal rigid plane. However, the actual contact situation is that the morphologies of both rough surfaces have fractal characteristics, and the distribution characteristics of normal contact area should be affected by the fractal dimensions of these two rough surfaces, so the characteristic influence of the other rough surface cannot be neglected. In this paper, first, the three-dimensional (3D) surface function (Y−K function) and two traditional area distribution functions are introduced (Section 2). Second, a logarithmic cumulative area function containing parameter ξ is proposed to characterize the normal contact area distribution, and it is compared with the two traditional area distribution functions. Based on the equivalence relations, a fractal dimension conversion formula between the contact of single rough surface and that of dual rough surfaces is acquired, thereby deriving an optimized expression for the fractal surface area distribution function (Section 3). Finally, the effectiveness and feasibility of this computational method are validated through contact analysis of the measured morphologies of anchor plates (APs) and anchor bearing plates (ABPs) (Section 4).

2 Surface morphology description and simplified contact method

2.1 Description of rough surface morphology

Previous scholars have found that the rough surface profile curves are statistically similar, i.e., irregular fractal characteristics, by studying the measured profile of the rough surface. Ausloos and Berman extended the two-dimensional W−M function in the polar coordinate system to a three-dimensional W−M function [34], on which Yan and Komvopoulos proposed a formula (Y−K function) [5] for the three-dimensional surface morphology in the Cartesian coordinate system by introducing the cutoff length, which is denoted by

(1)\begin{array}{rl}& \mathit{\text{z}}\left(x,y\right)=L{\left(\frac{G}{L}\right)}^{D-2}{\left(\frac{\ln \gamma }{M}\right)}^{\frac{1}{2}}\sum _{m=1}^M\sum _{n=0}^{n_{max}}{\gamma }^{\left(D-3\right)n}\\ & \times \left\{\cos {\varphi }_{m,n}-\cos \left[\frac{2\text{π}{\gamma }^n\sqrt{x^2+y^2}}{L}\cos \left(\arctan \left(\frac{y}{x}\right)-\frac{\text{π}m}{M}\right)+{\varphi }_{\mathrm{m},n}\right]\right\},\end{array}

where z(x,y) is the surface morphology height, x and y are the corresponding coordinates, D (2 < D < 3) is the fractal dimension, G is the fractal scale coefficient, L is the sampling length, L_s is the cutoff length, φ_m,n is the random phase, φ_m,n∈[0, 2π], M is the number of superposed rumbles on surface morphology, γ_n is the spatial frequency, which is generally taken to be 1.5 for the surface satisfying the normal distribution, and n is the frequency index of the asperity, n_max = n_min + int(lg (L/L_s)/lg γ).

2.2 Traditional area distribution function

The rough surface contact model proposed by Majumdar and Bhushan [4] is to simplify the contact of two rough surfaces to the contact of a rough surface and a rigid plane (Fig. 1). And it ignores the variation of material hardness with surface depth and the interaction generated by the deformation of neighboring asperities. This M−B model is based on the W−M fractal function and the global island cross-sectional area distribution rule, and it is assumed that the curvature radius of the asperity is related to the surface profile, which makes the model scale-independent.

The cumulative quantity function N(a) of the normal contact area proposed in the M−B model is shown in Eq. (2), and its differentiation leads to the normal contact area distribution function n(a) [35] shown in Eq. (3):

(2)N\left(A\ge a\right)={\left(a_{\text{L}}/a\right)}^{\left(D-1\right)/2},

(3)n\left(a\right)=\left|\frac{\text{d}N\left(a\right)}{\text{d}a}\right|=\frac{D}{2}{a}_{\text{L}}^{\frac{D}{2}}{a}^{-\frac{D+2}{2}},

where a and A are the normal contact areas, a_L is the maximum normal contact area, D is fractal dimension (1 < D < 2).

Komvopoulos proposed a reasonable approximate extension of the nonzero domain of the normal contact area distribution function [36], which in turn led to a more accurate modified function for the normal contact area distribution on fractal surfaces. The normal contact area distribution function expressed by Eq. (4) is still only related to the fractal dimension D and the maximum normal contact area a_L.

(4)n\left(a\right)=\frac{D}{2}{\psi }^{\frac{2-D}{2}}{a}_{\text{L}}^{\frac{D}{2}}{a}^{-\frac{D+2}{2}},

where ψ is the regional expansion factor.

3 Contact analysis of single/dual rough surfaces

In this section, the cut-off area is regarded as the contact area, that is the normal contact area is regarded as the actual contact area, which is consistent with the assumption of the current fractal contact model. First, the calculation method of the equivalent area distribution is introduced; then, the area distributions of single and dual rough surfaces are analyzed, respectively; finally, the optimized equivalent area distribution function is established through the equivalent fractal dimension conversion matrix.

3.1 The calculation method of the equivalent area distribution

When two rough surfaces are in contact, it is not only related to their own morphology parameters, but also affected by the pressure grade. As shown in Fig. 2(a), the ratio of the contact distance Δh to the height h_max of the maximum asperity on one surface is used to express the contact grade of surface-to-plane or surface-to-surface, i.e., λ = Δh/h_max.

1) First, the contact analysis of single rough surface is performed. The area a_i of each contact region with different fractal dimensions D_1,s-p at different pressure grades λ is counted to calculate a_L/a_i and N (inverse rank sum of a_i).

2) The parameter ξ_s-p that characterizes the relationship of N–a_L/a_i is obtained by fitting this data set, which in turn establishes the relationship between the parameter ξ_s-p and the fractal dimensions D_1,s-p and the pressure grade λ.

3) Then, the contact analysis of dual rough surfaces is performed. The area a_i of each contact region with different fractal dimensions D_1,s-s and D_2,s-s at different pressure grades λ is counted to obtain a_L/a_i and N.

4) The parameter ξ_s-s that characterizes the relationship of N–a_L/a_i is obtained by fitting this data set, which in turn establishes the relationship between the parameter ξ_s-s and the two fractal dimensions D_1,s-s, D_2,s-s, and the pressure grade λ.

5) Finally, the equivalent fractal dimension is obtained by the equivalent conversion of D_1,s-p, D_1,s-s, and D_2,s-s according to ξ_s-p = ξ_s-s, which in turn achieves the improved area distribution function n(a) that satisfies the equivalent contact calculations.

3.2 Results analysis on the contact between single/dual rough surfaces

3.2.1 Contact analysis of single rough surface

When surface 1 is in contact with a rigid plane (Fig. 3), the area a_i of each contact region can be obtained by calculating the connectivity domain, and the number N of corresponding areas not less than a_i can be achieved by solving the inverse rank sum of each area ai, and the data set of N–a_L/a_i can be established from this.

Figure 4 gives the cumulative quantitative relationship of contact area for Y−K surfaces with fractal dimension D_1,s-p = 2.3 at different pressure grades (λ = 30%, 50%, and 70%). Table 1 shows the fitting results, deviations and root mean square error (RMSE) of these two traditional area distribution functions in the effective range. The fractal dimension D involved in the calculations takes the value range of (2,3). When the value range of a_L/a is large, the obtained fitted fractal dimension D deviates from (2,3), that is, ineffective fitting. Therefore, the value range of a_L/a needs to be reduced to satisfy the effective fitting of these two formulas. Although the fitted fractal dimension D of these two formulas is able to approach the target fractal dimension D to some extent better, it is still limited by the smaller value range of a_L/a, which ignores the main role of the larger asperity in the contact process.

According to the distribution characteristics of N–a_L/a_i and the requirement that it must pass through (1,1), a logarithmic-form cumulative contact area function containing parameter ξ is proposed as shown in Eq. (5). The fitting of N–a_L/a_i using this function is shown in Fig. 4, and this function has good fitting effect for all value ranges of a_L/a_i. The RMSEs of the fitting results at different pressure grades (λ = 30%, 50%, and 70%) are 0.716, 1.881, and 1.073, respectively, which have higher confidence compared to the two traditional formulas (Table 1).

(5)N=\xi \cdot \ln \left(\frac{a_L}{a}\right)+1,

where ξ is the parameter to be determined.

Let G = 1 × 10⁻¹⁰ m, n_max = 64, M = 50, L = 10⁻³ m, L_s = 10⁻⁶ m to solve the cumulative quantitative relationship of contact area of Y−K surfaces with different fractal dimensions D_1,s-p and pressure grades λ. And Eq. (5) is fitted to derive the corresponding parameter ξ_s-p, which in turn achieves the relationship between the parameter ξ_s-p and the fractal dimension D_1,s-p and the pressure grade λ as shown in Fig. 5.

It can be seen that the parameter ξ_s-p exhibits the relationship that increases first and decreases later with the pressure grade λ. The peak point of the parameter ξ_s-p corresponding to the pressure grade λ is mainly concentrated around 60%, with a sensitivity range of 50%−80% (Fig. 5(a)). This is approximately consistent with a normal distribution (Fig. 5(b)), specifically the relationship between the parameter ξ_s-p and the pressure grade λ satisfying Eq. (6) for the surface with the same fractal dimension D_1,s-p. Figure 6 shows the relationship between the normalized parameter ξ_s-p/ξ_s-p,max and the pressure grade λ, which shows a stronger regularity.

(6)\xi /{\xi }_{\text{max}}=N_{\lambda }\left(\mu ,\sigma \right)\cdot \sqrt{2\text{π}}\sigma ={N^{\mathrm{\prime }}}_{\lambda }\left(\mu ,\sigma \right)={\text{e}}^{-\frac{{\left(\lambda -\mu \right)}^2}{2{\sigma }^2}},

where ξ_max is the maximum value of parameter; N_λ(μ,σ) is the normal distribution function for the pressure grade λ, where μ and σ are the mean and standard deviation, respectively; and N_{\lambda }^{\prime }(μ,σ) is the normalized normal distribution function.

In addition, the parameter ξ_s-p at the same pressure grade λ shows an approximately exponential growth relationship with the fractal dimension D_1,s-p (Fig. 5(d)), and the closer the pressure grade λ is to 60%, the larger the parameter ξ_s-p is for the surface with the same fractal dimension D_1,s-p. This is consistent with the actual situation, the larger the fractal dimension, the finer the asperity on the surface, the greater the contact number. Figures 5(c) and 5(e) show the partial fitting results of N–a_L/a_i in Figs. 5(b) and 5(d), respectively. The proposed function has a high agreement in the larger value range of a_L/a, and still has the good compatibility in the smaller value range of a_L/a, which can greatly reflect the main role of the larger asperity in the contact process.

In summary, it can be obtained that the parameter ξ_s-p in the contact analysis of single rough surface shows exponential growth with the fractal dimension D_1,s-p, and normal distribution characteristics with the pressure grade λ. As these results, Eq. (7) is established to express the relationship between the parameter ξ_s-p and the fractal dimension D_1,s-p and the pressure grade λ. Figure 7 shows the fitting results of Eq. (7) for the data in Fig. 5(a), which demonstrates high agreement (Goodness of fit R² = 0.9687), thereby validating the accuracy of the proposed function.

(7)\xi ={\text{e}}^{aD+b}\cdot N{{}^{\prime }}_{\lambda }\left(\mu ,\sigma \right)={\text{e}}^{aD+b}\cdot {\text{e}}^{-\frac{{\left(\lambda -\mu \right)}^2}{2{\sigma }^2}},

where ξ is the parameter, D and λ are the fractal dimension and pressure grade, respectively, N_{\lambda }^{\prime }(\mu ,\sigma ) is the normalized normal distribution function, a, b, μ, and σ are the parameters to be determined.

3.2.2 Contact analysis of dual rough surfaces

When surface 1 is in contact with surface 2 (Fig. 8), the area a_i of each contact region can be obtained by calculating the connectivity domain, and the number N of corresponding areas not less than a_i can be achieved by solving the inverse rank sum of each area ai, and the data set of N–a_L/a_i can be established from this.

The contact of dual rough surfaces is more stochastic compared to the contact of single rough surface, and the former has one more fractal dimension for the characterization of surface 2 compared to the latter. To study the relationship between these two more accurately, the same surface 1 is still used and other parameters are kept the same as before, i.e., G = 1 × 10⁻¹⁰ m, n_max = 64, M = 50, L = 10⁻³ m, and L_s = 10⁻⁶ m. On this basis, the cumulative quantitative relationship of the contact area is solved, and the parameter ξ_s-s is obtained by fitting using Eq. (5), and then the relationship between the parameter ξ_s-s and the different fractal dimensions D_1,s-s, D_2,s-s and the pressure grade λ is obtained, as shown in Fig. 9.

Figures 9(a)–9(c) show the relationship between the parameter ξ_s-s and the fractal dimension D_1,s-s and the pressure grade λ when the fractal dimension D_2,s-s = 2.5. The pressure grade λ corresponding to the maximum of parameter ξ_s-s is mainly concentrated around 60%, with a sensitivity range of 40%–70% (Fig. 9(a)), which is similar to the conclusions of the previous section. The fluctuation of the parameter ξ_s-s decreases as the fractal dimension D_1,s-s increases (Fig. 9(b)). This is because when D_1,s-s ≥ D_2,s-s, surface 1 tends to be more plane compared to surface 2, that is, the contact between surfaces 1 and 2 tends to be more similar to the contact between surface 2 and plane. Whereas, when D_1,s-s ≤ D_2,s-s, surface 2 tends to be more close to the plane compared to surface 1, which is in agreement with the contact between surface 1 and plane in the previous section, and the parameter ξ_s-s shows the relationship with the pressure grade λ that increases first and then decreases, which is similar to the normal distribution.

Figures 9(d)–9(f) show the relationship between the parameters ξ_s-s and the fractal dimensions D_1,s-s, D_2,s-s when the pressure grade λ = 60%. For similar values of fractal dimensions D_1,s-s and D_2,s-s, the parameter ξ_s-s increases as the fractal dimension D_1,s-s or D_2,s-s increases (Fig. 9(d). This is due to the fact that the finer the asperity on the surface with larger fractal dimension, the larger the contact region number. Furthermore, Fig. 9(c) and 9(f) show the partial fitting results of N–a_L/a_i in Fig. 9(b) and 9(c), respectively. It can be seen that the proposed function exhibits high consistency for all ranges of a_L/a.

3.2.3 Dimensional transformation relation and area distribution function

The dual rough surfaces have two different fractal dimensions, whereas the original cumulative quantity function N(a) [35,36] has only one fractal dimension. Therefore, the key issue is how to utilize the fractal dimensions of dual rough surfaces to describe the area distribution, that is how to equate the fractal dimensions of dual rough surfaces to the fractal dimension of single rough surface.

To make the contact more general, the contact behavior of 1000 sets of single rough surface (surface 1 and plane) and dual rough surfaces (surface 1 and 2) was investigated. It is important to emphasize that the surface 1 involved in both types of contact for each set is the same. The parameters ξ resulting from these two types of contacts (single and dual rough surfaces) are fitted using Eq. (7) and the results obtained are shown in Eqs. (8) and (9), respectively. These two have high agreement, which shows that the proposed calculation method has good applicability. In addition, it can be found that the mean of these two are similar (μ_s-p = 0.65, μ_s-s = 0.68) and the standard deviation of these two is nearly twice as large (σ_s-p = 0.22, σ_s-s = 0.39). This is due to the fact that the contact of dual rough surfaces is more random than the contact of single rough surface, which leads to an increase in the standard deviation.

(8){\xi }_{\text{s-p}}=e^{4.32D_{1,\text{s-p}}-9.60}\cdot N{{}^{\prime }}_{\lambda }\left(0.68,0.22\right),\left(R^2=0.9055\right),

(9)\begin{array}{rl}{\xi }_{\text{s-s}}=& e^{-6.26D_{1,\text{s-s}}{}^2+14.27D_{1,\text{s-s}}{D}_{2,\text{s-s}}-6.26D_{2,\text{s-s}}{}^2-2.77D_{1,\text{s-s}}-2.77D_{2,\text{s-s}}+3.89}\\ & \cdot N{{}^{\prime }}_{\lambda }\left(0.65,0.39\right),\left(R^2=0.8801\right).\end{array}

Figures 10(a) and 10(b) show the relationship between the parameter ξ and the fractal dimension D_i for the contact of single and dual rough surfaces at the pressure grade λ = 65%, respectively. According to ξ_s-p = ξ_s-s, the fractal dimensions D_1,s-s, D_2,s-s used for the contact between dual rough surfaces are equivalent to the fractal dimensions D_1,s-p used for the contact between surface and plane, whose expression is shown in Eq. (10), and the relationship is shown in Fig. 10(c). When the fractal dimensions D_1,s-s and D_2,s-s take similar values, the equivalent fractal dimension D_eff increases with the increase of the fractal dimension D_1,s-s or D_2,s-s; while when the fractal dimensions D_1,s-s and D_2,s-s take slightly different values, the equivalent fractal dimension D_eff decreases.

(10)\begin{array}{rl}{D}_{\text{eff}}=& D_{1,\text{s-p}}=\left[\begin{array}{ccc}1& D_{\text{1,s-s}}& D_{\text{1,s-s}}{}^2\end{array}\right]\mathit{C}\left[\begin{array}{c}1\\ D_{\text{2,s-s}}\\ D_{\text{2,s-s}}{}^2\end{array}\right]\\ =& \left[\begin{array}{ccc}1& D_{\text{1,s-s}}& D_{\text{1,s-s}}{}^2\end{array}\right]\left[\begin{array}{ccc}3.1236& -0.6415& -1.4501\\ -0.6415& 3.3046& 0\\ -1.4501& 0& 0\end{array}\right]\\ & \cdot \left[\begin{array}{c}1\\ D_{\text{2,s-s}}\\ D_{\text{2,s-s}}{}^2\end{array}\right],\end{array}

where D_eff is the equivalent fractal dimension, C is the third-order transformation matrix.

Since the fractal dimensions D_1,s-p, D_1,s-s, D_2,s-s∈(2, 3), two ineffective ranges (Fig. 10(c)) arise from this Eq. (10), which is due to the large difference between the two fractal dimensions D_1,s-s and D_2,s-s. When the difference between the two fractal dimensions is large, the surface with larger fractal dimension tends to be more plane compared to the surface with smaller fractal dimension, at which time the surface with smaller fractal dimension dominates the contact. Equation (11) is the improved equivalent fractal dimension D_eff. Based on this, Fig. 10(d) shows the comparison between the actual ξ generated by the contact of dual rough surfaces and the equivalent ξ generated by the contact of a rough surface with fractal dimension D_eff and a plane. The actual ξ and equivalent ξ are essentially comparable, with a RMSE of 1.0789, which indicates that the proposed method is effective in terms of equivalent contact.

(11)\begin{array}{rl}& D_{\text{eff}}=D_{1,\text{s-p}}\\ =& \{\begin{array}{l}\begin{array}{rl}& \left[\begin{array}{ccc}1& D_{\text{1,s-s}}& D_{\text{1,s-s}}{}^2\end{array}\right]\\ & \cdot \left[\begin{array}{ccc}3.1236& -0.6415& -1.4501\\ -0.6415& 3.3046& 0\\ -1.4501& 0& 0\end{array}\right]\left[\begin{array}{c}1\\ D_{\text{2,s-s}}\\ D_{\text{2,s-s}}{}^2\end{array}\right],\\ & \left|D_{1,\text{s-s}}-D_{2,\text{s-s}}\right|⩽0.2\end{array}\\ min\left(D_{1,\text{s-s}},D_{2,\text{s-s}}\right),\left|D_{1,\text{s-s}}-D_{2,\text{s-s}}\right|>\mathrm{0.2.}\end{array}\end{array}

The cumulative quantity function N(a) (Eq. (12)) of the contact area for the contact between surface and plane can be obtained by bringing Eq. (8) into Eq. (5). Then the distribution function n(a) of contact area as shown in Eq. (13) can be obtained by differentiation of Eq. (12). The area distribution function n(a) (Eq. (14)) for the contact between two rough surfaces equivalent to the contact between surface and plane can be obtained by bringing Eq. (11) into Eq. (13).

(12)N\left(a\right)=e^{4.32D_{1,\text{s-p}}-9.60}\cdot N{{}^{\prime }}_{\lambda }\cdot \ln \left(a_{\text{L}}/a\right)+1,

(13)n\left(a\right)=\left|\text{d}N\left(a\right)/\text{d}a\right|={\text{e}}^{4.32D_{1,\text{s-p}}-9.60}\cdot N{{}^{\prime }}_{\lambda }/a,

(14)n\left(a\right)=\{\begin{array}{l}{\text{e}}^{4.32\times \left[1D_{1,\text{s-s}}{D}_{1,\text{s-s}}^2\right]\mathit{C}{\left[1D_{2,\text{s-s}}{D}_{2,\text{s-s}}^2\right]}^{\text{T}}-9.60}\\ \cdot N{{}^{\prime }}_{\lambda }/a,\left|D_{1,\text{s-s}}-D_{2,\text{s-s}}\right|⩽0.2\\ {\text{e}}^{4.32\times min\left(D_{1,\text{s-s}},D_{2,\text{s-s}}\right)-9.60}\cdot N{{}^{\prime }}_{\lambda }/a,\left|D_{1,\text{s-s}}-D_{2,\text{s-s}}\right|>\mathrm{0.2.}\end{array}

4 Numerical contact analysis of anchorage systems

To verify the applicability of the above calculation method, the contact between the AP and the ABP will be analyzed in the following.

4.1 Surface morphologies of anchorage system

The anchor bearing system consists of an AP and an ABP, of which the material of the AP is 45# steel (specification: M15-7BC) and the material of the ABP is HT-200 gray cast iron (specification: YM15-7D). The surface morphologies of the anchorage system (AP and ABP) were measured using the auto-zoom 3D surface measuring instrument (Alicona IFMG5). The realistic degree of surface morphology is affected by the sampling distance, and the smaller the sampling distance is, the closer it is to the real surface. Referring to the similar method in Ref. [37]. and considering the sampling parameters of the measuring instrument, it was determined that the sampling distance was defined as 1.763212 μm, and the test region was 5.924 mm × 5.924 mm. These new morphologies shown in Fig. 11 were obtained after preprocessing the initial morphology (plane rotation, noise reduction by wavelet method and smoothing by Savitzky-Golay filter, etc.).

4.2 Calculation of the fractal dimension

The Variogram method (Eq. (15)) [38,39] for calculating the fractal dimension is an analytical method based on the variability of spatial data, where the fractal dimension is estimated by the power-law property of the variogram function. The fractal dimension is calculated as D = 3–H for surfaces.

(15)\text{log}\left\{E\left[{\left(B\left(x_2,y_2\right)-B\left(x_1,y_1\right)\right)}^2\right]\right\}=2H\log \left(k\right)+\log \left(C\right),

where E(∙) is the cumulative contribution function, B(∙) is the fractal Brown function; and the parameters H and C are usually estimated by the fitted slope and intercept of the double logarithmic curve of E[(B(x₂,y₂)–B(x₁,y₁))²] on k (the Euclidean distance k = [(x₂–x₁)² + (y₂–y₁)²]^1/2) is estimated.

Figures 12(a) and 12(b) show the double logarithmic relationships (lgE–lgk) for the morphology of AP and ABP after the variational function method, respectively. It can be seen that the double logarithmic relationships of the ABP morphology for the x- and y-directions are extremely close to each other, while the double logarithmic relationships of the AP morphology for the x- and y-directions are somewhat different. This is due to the specific profile shape (ripple shape, Fig. 11(a)) of the AP. The fractal dimensions of the AP and the ABP are calculated to be 2.428 and 2.307, respectively.

4.3 Analysis of contact result

Based on the conclusions of the previous section, when two rough surfaces are in contact, if the fractal dimension of one surface is lower (that is, the surface is rougher), the distribution of contact area is primarily determined by the surface with the lower fractal dimension. When a relatively smooth surface and a rough surface are in contact, the distribution of actual contact area is primarily influenced by the rough surface, which is in accordance with reality. Therefore, when the ABP and the AP are in contact, the ABP with a lower fractal dimension plays the dominant role.

Figures 13(a) and 13(b) show the contact between the ABP and the plane, and between the ABP and the AP, respectively. Both of them exhibit similar distribution characteristics. The parameter ξ_s-p generated by the contact between ABP and plane are numerically slightly lower than the parameter ξ_s-s generated by the contact between AP and ABP. This is because the contact between two surfaces has more contact possibilities compared to the contact between a surface and a plane. Bringing D_AP = 2.428 and D_ABP = 2.307 into Eq. (11) yields D_eff = 2.330. Figure 14 shows the comparison relationship between the parameters ξ before and after equivalence transformation and the actual ξ. The RMSEs for the non-equivalent ξ_s-p and equivalent ξ_s-p are 1.2386 and 0.3175, respectively, which shows that the proposed method is reliable for equivalent contact.

5 Conclusions

An improved area distribution function for characterizing the contact between two rough surfaces was established based on the proposed equivalent fractal dimension conversion relationship. This method was validated through contact analysis of the anchorage system. The main conclusions are as follows.

1) The proposed cumulative number function N(a) of the contact area in logarithmic form has a high agreement in the whole range. It overcomes the failure problem of the traditional function for larger area ratios a_L/a, and is especially effective in characterizing the dominant role of large asperity.

2) In the contact between surface and plane, the parameter ξ shows an exponential relationship with the fractal dimension D, and obeys a normal distribution rule with the pressure grade λ. In the contact between two rough surfaces, the distribution characteristics of the parameter ξ are influenced by the difference in the fractal dimensions of these two surfaces, which dominates the contact behavior by the surface with the smaller fractal dimension.

3) The equivalent fractal dimension model is realized for the contact between single/dual rough surfaces by means of a third-order transformation matrix. On this basis, the area distribution function n(a) considering the fractal dimension D_i and the pressure grade λ is established.

4) The contact analysis of the anchorage system shows that the distribution rule of the parameter ξ is highly consistent with the conclusion of the numerical prediction, which verifies the applicability of the equivalent fractal dimension model in engineering contact analysis.

This study provides a more universal fractal model framework for complex surface contact behavior. It has important theoretical value for engineering fields such as solid mechanics, tribology, and seal design. In the future, the construction of the dynamic contact model will be further improved by combining multi-scale conditions.

References

[1] Cheng S , Meng X , Li R , Liu R , Zhang R , Sun K , Ye W , Zhao F . Rough surface damping contact model and its space mechanism application. International Journal of Mechanical Sciences, 2022, 214: 106899

[2] Sun Y , Xiao H , Xu J . Investigation into the interfacial stiffness ratio of stationary contacts between rough surfaces using an equivalent thin layer. International Journal of Mechanical Sciences, 2019, 163: 105147

[3] Chen J , Liu D , Wang C , Zhang W , Zhu L . A fractal contact model of rough surfaces considering detailed multi-scale effects. Tribology International, 2022, 176: 107920

[4] Majumdar A , Bhushan B . Role of fractal geometry in roughness characterization and contact mechanics of surfaces. Journal of Tribology, 1990, 112(2): 205–216

[5] Yan W , Komvopoulos K . Contact analysis of elastic–plastic fractal surfaces. Journal of Applied Physics, 1998, 84(7): 3617–3624

[6] Miao X , Huang X . A complete contact model of a fractal rough surface. Wear, 2014, 309(1-2): 146–151

[7] Liou J , Tsai C , Lin J . A microcontact model developed for sphere- and cylinder-based fractal bodies in contact with a rigid flat surface. Wear, 2010, 268(3-4): 431–442

[8] Liou J L , Lin J F . A modified fractal microcontact model developed for asperity heights with variable morphology parameters. Wear, 2010, 268(1): 133–1449

[9] Xu K , Yuan Y , Chen J . The effects of size distribution functions on contact between fractal rough surfaces. AIP Advances, 2018, 8: 075317

[10] Yuan Y , Cheng Y , Liu K , Gan L . A revised Majumdar and Bushan model of elastoplastic contact between rough surfaces. Applied Surface Science, 2017, 425: 1138–1157

[11] Zhang X , Shen H , Zhang F , Zhang G , Xia H , Liu J , Ao X , Luo J , Zhu X , Zhang J . et al. A multi-scale contact model for rough surfaces under loading and unloading conditions. Tribology Letters, 2025, 73(2): 70

[12] Wang Y H , Zhang X L , Wen S H , Chen Y . Fractal loading model of the joint interface considering strain hardening of materials. Advances in Materials Science and Engineering, 2019, 2019(1442): 1–14

[13] Tikhomirov V P , Izmerov M A . Fractal model of contact interaction of nominally flat surfaces. AIP Conference Proceedings, 2021, 2340: 060004

[14] Wang R , Zhu L , Zhu C . Research on fractal model of normal contact stiffness for mechanical joint considering asperity interaction. International Journal of Mechanical Sciences, 2017, 134: 357–369

[15] Zhao Y S , Song X L , Cai L G , Liu Z , Cheng Q . Surface fractal topography-based contact stiffness determination of spindle–toolholder joint. Proceedings of the Institution of Mechanical Engineers. Part C: Journal of Mechanical Engineering Science, 2016, 230(4): 602–610

[16] Pan W , Li X , Wang L , Guo N , Yang Z . Influence of contact stiffness of joint surfaces on oscillation system based on the fractal theory. Archive of Applied Mechanics, 2017, 88(4): 525–541

[17] Pan W , Sun Y , Li X , Song H , Guo J . Contact mechanics modeling of fractal surface with complex multi-stage actual loading deformation. Applied Mathematical Modelling, 2024, 128: 58–81

[18] Yuan Y , Xu K , Zhao K . The loading–unloading model of contact between fractal rough surfaces. International Journal of Precision Engineering and Manufacturing, 2020, 21(6): 1047–1063

[19] Yuan Y , Chen J , Zhang L . Loading–unloading contact model between three-dimensional fractal rough surfaces. AIP Advances, 2018, 8(7): 075017

[20] Yuan Y , Xu K , Zhao K . A fractal model of contact between rough surfaces for a complete loading–unloading process. Proceedings of the Institution of Mechanical Engineers. Part C: Journal of Mechanical Engineering Science, 2020, 234(14): 2923–2935

[21] Yu X , Sun Y , Zhao D , Wu S . A revised contact stiffness model of rough curved surfaces based on the length scale. Tribology International, 2021, 164: 107206

[22] Liu Y , Wang Y , Chen X , Yu H . A spherical conformal contact model considering frictional and microscopic factors based on fractal theory. Chaos, Solitons, and Fractals, 2018, 111: 96–107

[23] Shen J , Xu S , Liu W , Yang J . Fractal model of normal contact stiffness between two spheres of joint interfaces with simulation. Mechanics, 2017, 23(5): 703–713

[24] Li X , Yue B , Wang D , Liang Y , Sun D . Dynamic characteristics of cylinders’ joint surfaces considering friction and elastic–plastic deformation based on fractal theory. Australian Journal of Mechanical Engineering, 2017, 15(1): 11–18

[25] Yuan Y , Gan L , Liu K , Yang X . Elastoplastic contact mechanics model of rough surface based on fractal theory. Chinese Journal of Mechanical Engineering, 2016, 30(1): 207–215

[26] Sun J , Ji Z , Zhang Y , Yu Q , Ma C . A contact mechanics model for rough surfaces based on a new fractal characterization method. International Journal of Applied Mechanics, 2018, 10(06): 1850069

[27] Tian H , Chen B , He K , Dong Y , Zhong X , Wang X , Xi N . A new fractal model of elastic, elastoplastic and plastic normal contact stiffness for slow sliding interface considering dynamic friction and strain hardening. Journal of Shanghai Jiaotong University (Science), 2017, 22(5): 589–601

[28] Lan G , Sun W , Zhang X , Chen Y , Tan W , Li X . A three-dimensional fractal model of the normal contact characteristics of two contacting rough surfaces. AIP Advances, 2021, 11: 055023

[29] Zhou Y , Xiao Y , He Y , Zhang Z . A detailed finite element analysis of composite bolted joint dynamics with multiscale modeling of contacts between rough surfaces. Composite Structures, 2020, 236: 111874

[30] Zhao Y , Yang C , Cai L , Shi W , Liu Z . Surface contact stress-based nonlinear virtual material method for dynamic analysis of bolted joint of machine tool. Precision Engineering, 2016, 43: 230–240

[31] Zhang Z , Xiao Y , Xie Y , Su Z . Effects of contact between rough surfaces on the dynamic responses of bolted composite joints: Multiscale modeling and numerical simulation. Composite Structures, 2019, 211: 13–23

[32] Zhou H , Long X , Meng G , Liu X . A stiffness model for bolted joints considering asperity interactions of rough surface contact. Journal of Tribology, 2022, 144(1): 011501

[33] Zhao Y , Xu J , Cai L , Shi W , Liu Z . Stiffness and damping model of bolted joint based on the modified three-dimensional fractal topography. Proceedings of the Institution of Mechanical Engineers. Part C: Journal of Mechanical Engineering Science, 2017, 231(2): 279–293

[34] Ausloos M , Berman D H . A multivariate Weierstrass-Mandelbrot function. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 1819, 1985(400): 331–350

[35] Mandelbrot B. The Fractal Geometry of Nature. New York, NY: W.H. Freeman and Company, 1977, 116–121

[36] Li X , Wang W , Zhao M . Fractal prediction model for tangential contact damping of jointsurface considering friction factors and its simulation. Journal of Mechanical Engineering, 2012, 48(23): 46–50

[37] Xie Y , Chen S , Wan X , Tse P . A preliminary numerical study on the interactions between nonlinear ultrasonic guided waves and a single crack in bone materials with motivation to the evaluation of micro cracks in long bones. IEEE Access, 2020, 8: 169169–169182

[38] Xu T , Moore L D , Gallant J C . Fractals, fractal dimensions and landscapes-a review. Geomorphology, 1993, 8(4): 245–262

[39] Kulatilake P H S W , Du S G , Ankah M L Y , Yong R , Sunkpal D T , Zhao X , Liu G , Wu R . Non-stationarity, heterogeneity, scale effects, and anisotropy investigations on natural rock joint roughness using the variogram method. Bulletin of Engineering Geology and the Environment, 2021, 80(8): 6121–6143

RIGHTS & PERMISSIONS

Higher Education Press