Normal contact stiffness of rough surfaces considering oblique asperity contact

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.

Keywords

Fractal model oblique asperity contact normal stiffness elastic–plastic force contact angle

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⁵ 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)⁶ 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.⁷ Bush, Gibson, and Thomas (BGT)⁸ modeled the asperities as paraboloids with two different radii of curvature. Greenwood and Tripp (GT)⁹ pointed out that the shape and the relative position of the asperities have no effect on the model. Whitehouse and Archard (WA)¹⁰ studied the correlation and the joint probability density between the heights and the curvatures of the summits. Chang, Etsion and Bogy (CEB)¹¹ proposed the concept of critical contact deformation and described the statistical model from GW model.

Persson’s theory^12,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.¹⁴ 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.¹⁵ 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¹⁹ 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.²⁰ depicted the interacting Hertzian asperities (IHA) by using the Johnson formulas. Aferrante et al.²¹ 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.¹⁹

Majumdar and Bhushan^23,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 Komvopoulos^24,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²⁶ (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.²⁹ 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.³⁰ 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³¹ 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.³² 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.³³ 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.³⁴ researched the elasto-plastic hemispherical force between spheres based on Kogut–Etsion (KE) model and characterized the component of friction. Sepehri and Farhang³⁵ extended the force between paraboloidal asperities defined in GT model and proposed oblique asperity contact model. Zhao et al.³⁶ 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.³⁷ The sliding interaction between asperities was also introduced considering substrate deformation and power-law hardening materials.³⁸

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⁷ 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 $1 / R = \sqrt{(1 / R_{upper}) (1 / R_{lower})}$ is used to calculate the equivalent radius rather than equation $1 / R = 1 / R_{upper} + 1 / R_{lower}$ .

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. $r'$ and r are the radii of the truncated area and the real micro-contact area, respectively. $δ$ is the asperity interference.

Figure 1.

Single asperity pair and its simplification.

The WK fractal model for the normal pressure p and the normal stiffness $k_{n}$ of a single pair are

p = {\begin{matrix} p_{e} = \frac{4}{3 \sqrt{2 π}} E' G^{D - 1} {(a')}^{1.5 - 0.5 D}, & a' > {a'}_{c} \\ p_{p} = k σ_{y} a', & a' < {a'}_{c} \end{matrix}

(1)

k_{n} = 2 E' \sqrt{\frac{a'}{2 π}}

(2)

where $E'$ is the equivalent elastic modulus, defined as $1 / E' = (1 - ν_{1}^{2}) / E_{1} + (1 - ν_{2}^{2}) / E_{2}$ , $a'$ is the truncated area, and $E_{1}$ , $ν_{1}$ and $E_{2}$ , $ν_{2}$ are the elastic moduli and Poisson ratios of two contact surfaces, respectively. D is the fractal dimension $(1 < D < 2)$ and G is the fractal roughness parameter of the surface profile. k is a coefficient and σ_y is the yield strength of the softer material. H=kσ_y where H is the hardness. $a_{c}$ is the critical area differentiating elastic and plastic regions. If it satisfies the relation $a' > a'_{c}$ , then the asperity behaves elasticity and $p_{e}$ is the elastic force. If $a' > a'_{c}$ , the spot is in plastic deformation and $p_{p}$ is the plastic force.

The superscript of $a'$ means the component of the truncated contact spot. The asperity radius R and the interference $δ$ are defined as²⁴

R = \frac{{(a')}^{0.5 D}}{2 π G^{D - 1}}

(3)

δ = G^{D - 1} {(a')}^{1 - 0.5 D}

(4)

$δ_{c}$ is the critical interference which can be obtained as⁹

δ_{c} = {(\frac{π k ϕ}{2})}^{2} R

(5)

where $ϕ = σ_{y} / E'$ . When $δ_{c} = δ$ , $a_{c}$ can be obtained by substituting equations (3) and (4) into equation (5)

G^{D - 1} {(a')}_{c}^{1 - 0.5 D} = {(\frac{π k ϕ}{2})}^{2} \frac{{(a')}_{c}^{0.5 D}}{2 π G^{D - 1}}

(6)

a'_{c} = G^{2} {(\frac{k ϕ}{2})}^{2 / (1 - D)} {(\frac{π}{2})}^{1 / (1 - D)}

(7)

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 $α$ as presented in Figure 3, the asperity contact force will yield two components. l is the length scale of an asperity and

l = 2 r = a^{0.5} = {(0.5 a')}^{0.5}

(8)

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³⁵

\cos α = {(1 + \frac{r_{a}^{2}}{R_{s}^{2}})}^{- \frac{1}{2}}

(9)

p_{n} = p_{a} \cos α

(10)

R (r) = R {(1 + \frac{r_{a}^{2}}{R_{s}^{2}})}^{\frac{3}{2}}

(11)

Because two asperities are not simplified in Aferrante et al.³⁵ and $R_{s}$ is the summation of radii of two asperities, r is the tangential deviation of two summits. In order to apply these expressions to the oblique contact as shown in Figure 3, R is equal to $R_{s}$ . $r_{a}$ is the tangential offset of the asperity summit and the contact point. Substituting equations (3), (8), and (11) into equations (5) to (7) results in

l'_{a} = l' {(1 + \frac{r_{a}^{2}}{R^{2}})}^{\frac{3}{2 D}}

(12)

\begin{matrix} δ_{a} = G^{D - 1} l'_{a}^{2 - D} = G^{D - 1} l'^{2 - D} {(1 + \frac{r_{a}^{2}}{R^{2}})}^{\frac{3 (2 - D)}{2 D}} \\ = δ {(1 + \frac{r_{a}^{2}}{R^{2}})}^{\frac{3 (2 - D)}{2 D}} \end{matrix}

(13)

δ_{ac} = {(\frac{π k ϕ}{2})}^{2} R {(1 + \frac{r_{a}^{2}}{R^{2}})}^{\frac{3}{2}}

(14)

Therefore, the normal load $p_{a}$ and critical area $a'_{ac}$ are

p_{a} = {\begin{matrix} p_{ae} = \frac{4}{3} E' R^{1 / 2} δ_{a}^{3 / 2} = \frac{4}{3} E' R^{1 / 2} δ^{3 / 2} {(1 + \frac{r_{a}^{2}}{R^{2}})}^{\frac{3 (3 - D)}{2 D}}, & a' > {a'}_{ac} \\ p_{ap} = k σ_{y} {a'}_{a} = k σ_{y} a' {(1 + \frac{r_{a}^{2}}{R^{2}})}^{\frac{3}{D}}, & a' < {a'}_{ac} \end{matrix}

(15)

a'_{ac} = G^{2} {(\frac{k ϕ}{2})}^{2 / (1 - D)} {(\frac{π}{2})}^{1 / (1 - D)} {(1 + \frac{r_{a}^{2}}{R^{2}})}^{\frac{3}{2 (1 - D)}}

(16)

The normal stiffness $k_{an}$ is obtained by

k_{an} = \frac{d p_{n}}{d δ_{n}} = \frac{d p_{n}}{d p_{a}} \frac{d p_{a}}{d δ_{a}} \frac{d δ_{a}}{d δ_{n}} = \frac{d p_{a}}{d δ_{a}} \cos^{2} α = k_{n} \cos^{2} α

(17)

Figure 2.

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

Figure 3.

Oblique asperity contact.

Multi-asperity normal fractal model

The size distribution density expression $n (a')$ from the WK model is given as²⁴

n (a') = 0.5 D ψ^{1 - 0.5 D} {(a')}_{l}^{0.5 D} {(a')}^{- 1 - 0.5 D}, 0 < a' \leq a'_{l}

(18)

where $ψ$ is the domain extension factor $(ψ > 1)$ and $a'_{l}$ is the truncated area of the largest micro-contact. The values of $ψ$ are related to the fractal dimension D $(ψ^{(2 - D) / 2} - (1 - ψ^{- D / 2})^{- (2 - D) / 2} = (2 - D) / D)$ . The real contact area $A_{r}$ is

\begin{matrix} A_{r} = \int_{0}^{{a'}_{l}} n (a') ada' = \int_{0}^{{a'}_{l}} n (a') 0.5 a' da' \\ = \frac{D ψ^{1 - 0.5 D}}{2 (2 - D)} a'_{l} = \frac{D ψ^{1 - 0.5 D}}{2 - D} a_{l} \end{matrix}

(19)

The contact area ratio is $A_{r}^{*} = A_{r} / A_{a}$ , where $A_{a}$ is an apparent area of the contact region. And then

\frac{a_{l}}{A_{a}} = \frac{a_{l}}{A_{r}} A_{r}^{*} = \frac{2 - D}{D ψ^{1 - 0.5 D}} A_{r}^{*}

(20)

When $a'_{l} > a'_{c}$ , the total normal elastic pressure $P_{e}$ can be calculated with equations (8), (15), and (18)

\begin{matrix} P_{e} = \int_{{a'}_{c}}^{{a'}_{l}} p_{e} (a') n (a') da' \\ = \frac{4}{3 \sqrt{2 π}} E' G^{D - 1} \int_{{a'}_{c}}^{{a'}_{l}} n (a') {(a')}^{1.5 - 0.5 D} da' \\ = {\begin{matrix} \frac{\sqrt{2} D ψ^{1 - 0.5 D}}{3 \sqrt{π} (1.5 - D)} E' G^{D - 1} {(a')}_{l}^{0.5 D} ({(a')}_{l}^{1.5 - D} - {(a')}_{c}^{1.5 - D}), D \neq 1.5 \\ \frac{{2.0007}^{0.25}}{\sqrt{2 π}} E' \sqrt{G} {(a')}_{l}^{0.75} \ln \frac{{a'}_{l}}{{a'}_{c}}, D = 1.5 \end{matrix} \\ = {\begin{matrix} \frac{2^{2 - 0.5 D} D ψ^{1 - 0.5 D}}{3 \sqrt{π} (1.5 - D)} E' G^{D - 1} a_{l}^{0.5 D} (a_{l}^{1.5 - D} - a_{c}^{1.5 - D}), D \neq 1.5 \\ {4.0014}^{0.25} E' \sqrt{\frac{G}{π}} a_{l}^{0.75} \ln \frac{a_{l}}{a_{c}}, D = 1.5 \end{matrix} \end{matrix}

(21)

When $A_{r}^{*} > A_{rc}^{*}$ , the dimensionless total normal elastic pressure $P_{e}^{*}$ can be written as

P_{e}^{*} = {\begin{matrix} \frac{2^{2 - 0.5 D} D ψ^{1 - 0.5 D}}{3 \sqrt{π} (1.5 - D)} {G^{*}}^{D - 1} {(\frac{2 - D}{D} A_{r}^{*})}^{0.5 D} [{(\frac{2 - D}{D ψ^{1 - 0.5 D}} A_{r}^{*})}^{1.5 - D} - a_{c}^{* 1.5 - D}], D \neq 1.5 \\ 2^{0.25} {2.0007}^{0.0625} \sqrt{\frac{G^{*}}{π}} (\frac{A_{r}^{*}}{3}) \ln \frac{A_{r}^{*}}{3 \times {2.0007}^{0.25}}, D = 1.5 \end{matrix}

(22)

where $P_{e}^{*} = P_{e} / E' A_{a}$ , $a_{c}^{*} = a_{c} / A_{a}$ , $A_{rc}^{*} = A_{rc} / A_{a}$ , $G^{*} = G / \sqrt{A_{a}}$ .

Similarly, the dimensionless total normal plastic pressure $P_{p}^{*}$ and dimensionless total normal stiffness $K_{n}^{*}$ , respectively, are

P_{p}^{*} = {\begin{matrix} 2 k ϕ {(A_{r}^{*})}^{0.5 D} {(\frac{D ψ^{1 - 0.5 D}}{2 - D} a_{c}^{*})}^{1 - 0.5 D}, D \neq 1.5 \\ 2 k ϕ {(A_{r}^{*})}^{0.75} {(3 \times {2.0007}^{0.25} a_{c}^{*})}^{0.25}, D = 1.5 \end{matrix}

(23)

\begin{matrix} K_{n}^{*} = \frac{2 D^{1 - 0.5 D} {(2 - D)}^{0.5 D}}{(1 - D) \sqrt{π}} ψ^{\frac{D^{2} - 4 D + 4}{4}} {(A_{r}^{*})}^{0.5 D} \\ [{(\frac{2 - D}{D})}^{0.5 - 0.5 D} ψ^{\frac{- D^{2} + 3 D - 2}{4}} {(A_{r}^{*})}^{\frac{1 - D}{2}} - {(a_{c}^{*})}^{\frac{1 - D}{2}}] \end{matrix}

(24)

where $K_{n}^{*} = K_{n} / (E' \sqrt{A_{a}})$ .

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 ${α_{1}, α_{2}, \dots, α_{i}, \dots, α_{i max}} (1 \leq i \leq i max)$ . If the contact point is the summit, $α_{i} = 0$ . While the contact point is in the limit position, the asperity does not bear the normal pressure and $r_{a} = R$ . Then, $\cos α_{i} = \sqrt{2} / 2$ and $α_{i} = π / 4$ . Hence, the contact angle range can be obtained as $0 \leq α_{i} < (π / 4)$ . Assume that the sum of the asperities contact area is $A_{i}$ when the contact angle is $α_{i}$ . Then, $\sum_{i = 1}^{i max} A_{i} = A_{r}$ . According to the definition of $n (a')$ , it can be seen that the area is distributed in all directions. So $\cos α_{i}$ is a constant for each contact angle and $1 + (r_{a}^{2} / R^{2})$ can be written as $c_{i}$ $(1 \leq c_{i} < 2)$ . When $a'_{l} > a'_{c}$ and the contact angle is $α_{i}$ , the total pressure $P_{ai}$ can be given as

\begin{matrix} P_{ai} = \frac{A_{i}}{A_{r}} (\int_{{a'}_{c}}^{{a'}_{l}} p_{aie} (a') n (a') da' + \int_{0}^{{a'}_{c}} p_{aip} (a') n (a') a' da') \\ = \frac{A_{i}}{A_{r}} (c_{i}^{\frac{3 (3 - D)}{2 D}} \int_{{a'}_{c}}^{{a'}_{l}} p_{e} (a') n (a') da' + c_{i}^{\frac{3}{D}} \int_{0}^{{a'}_{c}} p_{p} (a') n (a') a' da') \\ = \frac{A_{i}}{A_{r}} (c_{i}^{\frac{3 (3 - D)}{2 D}} P_{e} + c_{i}^{\frac{3}{D}} P_{p}) \end{matrix}

(25)

As $p_{n} = p_{a} \cos α$ , the total normal pressure $P_{ni}$ can be obtained as

P_{ni} = P_{ai} \cos α_{i} = P_{ai} c_{i}^{- \frac{1}{2}} = \frac{A_{i}}{A_{r}} (c_{i}^{\frac{9 - 4 D}{2 D}} P_{e} + c_{i}^{\frac{6 - D}{2 D}} P_{p})

(26)

The sum of the asperities of all contact angles is the total load of the joint surface

P_{na} = \sum_{i = 1}^{i max} P_{ni} = \sum_{i = 1}^{i max} \frac{A_{i}}{A_{r}} (c_{i}^{\frac{9 - 4 D}{2 D}} P_{e} + c_{i}^{\frac{6 - D}{2 D}} P_{p})

(27)

The dimensionless form is

P_{na}^{*} = \sum_{i = 1}^{i max} \frac{A_{i}}{A_{r}} (c_{i}^{\frac{9 - 4 D}{2 D}} P_{e}^{*} + c_{i}^{\frac{6 - D}{2 D}} P_{p}^{*})

(28)

Similarly, when the contact angle is $α_{i}$ , the total normal stiffness $K_{nai}$ is

K_{nai} = \frac{A_{i}}{A_{r}} \int_{{a'}_{c}}^{{a'}_{l}} k_{na} n (a') da' = \frac{A_{i}}{A_{r}} \cos^{2} α_{i} \int_{{a'}_{c}}^{{a'}_{l}} k_{a} n (a') da'

(29)

where $k_{a}$ is the single asperity contact pair stiffness along the normal direction of the contact point. According to the assumption, single oblique micro-contact pressure is the same as that in WK model. Then

p_{ae} = \frac{4}{3} E' R {(r)}^{1 / 2} δ_{a}^{3 / 2} = \frac{4}{3} E' R^{1 / 2} δ^{3 / 2} {(1 + \frac{r_{a}^{2}}{R^{2}})}^{\frac{3 (3 - D)}{2 D}}

(30)

\frac{\partial p_{ae}}{\partial δ_{a}} = 2 E' R^{1 / 2} δ^{1 / 2} {(1 + \frac{r_{a}^{2}}{R^{2}})}^{\frac{3 (3 - D)}{2 D}}

(31)

Using equation $r = R^{1 / 2} δ^{1 / 2}$ ,³⁵ $k_{a}$ is

\begin{matrix} k_{a} = \frac{\partial p_{ae}}{\partial δ_{a}} = 2 E' r_{a} = 2 E' r {(1 + \frac{r^{2}}{R^{2}})}^{\frac{3 (3 - D)}{2 D}} \\ = k_{n} {(1 + \frac{r_{a}^{2}}{R^{2}})}^{\frac{3 (3 - D)}{2 D}} \end{matrix}

(32)

Accordingly, the total normal stiffness can be written as

\begin{matrix} K_{na} = \sum_{i = 1}^{i max} \frac{A_{i}}{A_{r}} K_{nai} = \sum_{i = 1}^{i max} \frac{A_{i}}{A_{r}} c_{i}^{- 1} K_{ai} \\ = \sum_{i = 1}^{i max} \frac{A_{i}}{A_{r}} c_{i}^{- 1} \int_{{a'}_{c}}^{{a'}_{l}} k_{n} c_{i}^{\frac{3 (3 - D)}{2 D}} n (a') da' \\ = \sum_{i = 1}^{i max} \frac{A_{i}}{A_{r}} c_{i}^{\frac{9 - 5 D}{2 D}} \int_{{a'}_{c}}^{{a'}_{l}} k_{n} n (a') da' = \sum_{i = 1}^{i max} \frac{A_{i}}{A_{r}} c_{i}^{\frac{9 - 5 D}{2 D}} K_{n} \end{matrix}

(33)

The dimensionless form is

K_{na}^{*} = \sum_{i = 1}^{i max} \frac{A_{i}}{A_{r}} c_{i}^{\frac{9 - 5 D}{2 D}} K_{n}^{*}

(34)

Numerical simulations

Simulation preparation

According to equations (22) to (24), (28), and (34), once dimensionless total normal pressure $P_{n}^{*}$ is provided, dimensionless total normal stiffness $K_{na}^{*}$ can be obtained. From the modified model, it is obvious that the components $c_{i}^{\frac{9 - 5 D}{2 D}}$ and $c_{i}^{\frac{9 - 4 D}{2 D}}$ are the main difference compared with WK model. Thus, fractal dimension D and the shoulder–shoulder contact distribution are the simulation variables. The values of domain extension factor $ψ$ are obtained from literature and displayed in Table 1. G is given as $1.0 \times 10^{- 10} m$ .

Table 1.

Domain extension factor $ψ$ with the fractal dimension.

D	1.2	1.4	1.6	1.8
$ψ$	2.3014	2.0849	1.9283	1.8103

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.

Table 2.

Material properties.

	Elastic modulus,E(GPa)	Poissonratio, $ν$	Hardness,H (HV)
Values	180	0.3	255

The oblique contact distributions are assumed as Gaussian distribution and uniform distribution. The values of parameter $c_{i}$ are given as ${1, 1.1, 1.2, \dots, 1.8, 1.9}$ . For uniform distribution, $A_{i} / A_{r} = {0.1, 0.1, 0.1, \dots, 0.1, 0.1}$ . And for Gaussian distribution, when $c_{i} = 1.1$ is set as the average point and the standard deviation is given as 0.4, the contact area ratio can be evaluated to give

\begin{matrix} \frac{A_{i}}{A_{r}} = {0.2, 0.18, 0.16, 0.14, 0.1, 0.08, 0.06, \\ 0.035, 0.025, 0.02} \end{matrix}

(35)

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.

Figure 4.

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 $1 < D < 2$ , the coefficient of elastic pressure $c_{i}^{\frac{9 - 4 D}{2 D}}$ is smaller than that of plastic pressure $c_{i}^{\frac{6 - D}{2 D}}$ by equation (28). Third, the elastic pressure ratio rises to a constant rapidly. Once the asperities contact, the amount of the contact asperities is few and the contact area is much smaller than critical area $a'_{ac}$ . Consequently, plastic deformation occurs until the contact area is bigger than $a'_{ac}$ .

Figure 5.

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.

Figure 6.

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.⁴⁰ 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, M_u1 and M_v1. The point spacing $Δ x$ is 50 µm and each test area size is 3 mm × 3 mm. The topography height set is defined as $z (x_{i}) (i = 0, 1, \dots, N - 1)$ , and Figure 8 shows three-dimensional plot of M_u1 and G_u1 after filtering the rough points.

Figure 7.

Measuring field for point cloud data of the contact surface.

Figure 8.

Three-dimensional plot of two surfaces: (a) M_u1 and (b) G_u1.

The function of the structure function method is

S (τ) = \frac{1}{N - n} \sum_{i = 0}^{N - n} {[z (x_{i} + n) - z (x_{i})]}^{2}

(36)

where $N - n$ is the quantity of data pairs separated by an interval $τ$ and $τ = n \times Δ x$ . Then, the fractal dimension D and the fractal roughness parameter G can be extracted from the log–log plot. Thus, fractal dimension and fractal roughness are given by

D = \frac{4 - k_{τ}}{2}, \log G = \frac{h_{τ} - \log C}{2 (D - 1)}

(37)

where $k_{τ}$ and $h_{τ}$ are the slope and intercept of the fitting line between $\log (S (τ))$ and $\log (τ)$ . When $n = 500$ , the fitting line of structure function method of M_u1 for calculating fractal parameters is shown in Figure 9(a), and $k_{τ} = 1.2377$ and $h_{τ} = - 6.3714$ . Then, the fractal parameters can be obtained by equation (37). Each area is tested in this way and the fractal parameters are calculated. After averaging the fractal parameters of the test areas in the same machined surface, the mean values for the milling surface are found to be D = 1.4 and G = 7.9 × 10⁻¹⁰ mm, and the mean values for the grinding surface are found to be D = 1.9 and G = 4.8 × 10⁻⁵ mm.

Figure 9.

The results of structure function method for calculating fractal parameters: (a) M_u1 and (b) G_u1.

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 $K_{\exp}$ of the contact surface is calculated as

K_{ep} = - \frac{d P_{ep}}{d δ_{ep}} = - \frac{d (F_{epn} / A_{ep})}{d (δ_{sn} - δ_{\ln})}

(38)

where $P_{ep}$ is the total normal pressure of samples, and $δ_{ep}$ and $F_{epn}$ are the normal deformation and force, respectively. $δ_{sn}$ and $δ_{\ln}$ are the total deformation of the short and long samples, respectively. $A_{ep}$ is the apparent area of the contact region as $A_{ep} = π \times 9^{2} m m^{2}$ . Figure 11 is the measured data, and Figure 12 shows the deformation deviation of the short and long samples.

Figure 10.

The single equivalent sample testing experiment equipment.

Figure 11.

Normal deformations and pressures of different testing samples.

Figure 12.

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).

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 $K_{ep}^{*} = K_{ep} / (E_{ep} \sqrt{A_{ep}})$ , and the dimensionless total normal pressure is calculated by $P_{ep}^{*} = P_{ep} / (E_{ep} A_{ep})$ , where $E_{ep}$ is the elastic modulus of C45E4. Because the samples have to be preloaded to fix and joint with the test rig, the experiment data are not recorded in the initial stage. From Figure 13, we can see that the WK model derives a higher stiffness than the experiment over the entire range. In contrast, the modified WK models are closer to the experiment. Moreover, the errors of dimensionless normal stiffness between experiment results and fractal models can be calculated as $E r r o r (%) = | (K_{model}^{*} - K_{e p}^{*}) / K_{e p}^{*} |$ . Obviously, the errors of the WK model are large (Figure 14), while those of the modified models are smaller especially in the low pressure. When the pressure is high enough, the WK model has more accurate stiffness.

Figure 14.

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 $c_{i}^{\frac{9 - 4 D}{2 D}}$ is smaller than that of plastic pressure $c_{i}^{\frac{6 - D}{2 D}}$ . Finally, the modified fractal stiffness model was validated by experiments, including specimen surface observation and load-deformation measurement. The result revealed that stiffness of the proposed model was closer to the experiment, which indicated that the new model is appropriate and reasonable. On the other hand, it also indicates the limitation that the actual contact angle and size distribution of the shoulder–shoulder contact are more complicated than the assumed distributions, which need further research to find out the discipline and predict the normal contact stiffness more precisely.

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.

ORCID iD

Kai Zhang

References

Kono

Inagaki

Matsubara

et al . Stiffness model of machine tool supports using contact stiffness. Precis Eng 2013; 37: 650–657.

Soom

Serpe

CI.

Normal stiffness and damping at lightly loaded rough planar contacts. Tribol Int 2016; 100: 171–177.

Zhai

Huang

et al . Investigation of the contact performance of machined surface morphology. Tribol Int 2017; 107: 125–134.

Lundberg

Nordborg

Arteaga

IL.

The influence of surface roughness on the contact stiffness and the contact filter effect in nonlinear wheel-track interaction. J Sound Vib 2016; 366: 429–446.

Johnson

KL.

Contact mechanics. Cambridge: Cambridge University Press, 1987.

Greenwood

Williamson

Contact of nominally flat surfaces. Proc R Soc Lond A 1966; 295: 300–319.

Greenwood

JA.

A simplified elliptic model of rough surface contact. Wear 2006; 261: 191–200.

Bush

Gibson

Thomas

TR.

The elastic contact of a rough surface. Wear 1975; 35: 87–111.

Greenwood

Tripp

JH.

The contact of two nominally flat rough surfaces. Proc Inst Mech Eng 1970; 185: 625–633.

10.

Whitehouse

Archard

JF.

The properties of random surface of significance in their contact. Proc R Soc Lond A 1970; 316: 97–121.

11.

Chang

Etsion

Bogy

DB.

An elastic-plastic model for the contact of rough surfaces. J Tribol 1987; 109: 257–263.

12.

Persson

BNJ

. Elastoplastic contact between randomly rough surfaces. Phys Rev Lett 2001; a87: 116101.

13.

Persson

BNJ

. Theory of rubber friction and contact mechanics. J Chem Phys 2001; b115: 3840–3861.

14.

Yastrebov

Anciaux

Molinari

JF.

From infinitesimal to full contact between rough surfaces: evolution of the contact area. Int J Solids Struct 2015; 52: 83–102.

15.

Putignano

Afferrante

Carbone

et al . The influence of the statistical properties of self-affine surfaces in elastic contacts: a numerical investigation. J Mech Phys Solids 2012; 60: 973–982.

16.

Müser

Dapp

Bugnicourt

et al . Meeting the contact-mechanics challenge. Trib Lett 2017; 65: 118.

17.

Aferrante

Bottiglione

Putignano

et al . Carbone elastic contact mechanics of randomly rough surfaces: an assessment of advanced asperity models and Persson’s theory. Trib Lett 2018; 66: 75.

18.

Campañá

Müser

. Contact mechanics of real vs. randomly rough surfaces: a Green’s function molecular dynamics study. Europhys Lett 2007; 77: 38005.

19.

Stanley

Etison

Bogy

DB.

Adhesion of contacting rough surfaces in the presence of sub-boundary. J Tribol 1990; 112: 98–104.

20.

Ciavarella

Delfine

Demelio

A “re-vitalized” Greenwood and Williamson model of elastic contact between fractal surfaces. J Mech Phys Solids 2006; 54: 2569–2591.

21.

Aferrante

Carbone

Demelio

Interacting and coalescing Hertzian asperities: a new multiasperity contact model. Wear 2012; 278–279: 28–33.

22.

Carbone

A slightly corrected Greenwood and Williamson model predicts asymptotic linearity between contact area and load. J Mech Phys Solids 2009; 57: 1093–1102.

23.

Majumdar

Bushan

Role of fractal geometry in roughness characterization and contact mechanics of surfaces. J Tribol 1990; 112: 205–215.

24.

Wang

Komvopoulos

A fractal theory of the interfacial temperature distribution in the slow sliding regime: part I—elastic contact and heat transfer analysis. J Tribol 1994; 116: 812–823.

25.

Wang

Komvopoulos

A fractal theory of the interfacial temperature distribution in the slow sliding regime: part II—multiple domains, elastoplastic contacts and applications. J Tribol 1994; 116: 824–832.

26.

Yan

Komvopoulos

Contact analysis of elastic-plastic fractal surfaces. J Appl Phys 1998; 84: 3617–3624.

27.

Jana

Mitra

Sahoo

Dynamic contact interactions of fractal surfaces. Appl Surf Sci 2017; 392: 872–882.

28.

Hanaor

DAH

Gan

Einav

Contact mechanics of fractal surfaces by spline assisted discretization. Int J Solids Struct 2015; 59: 121–131.

29.

Xiao

Shao

Brennan

MJ.

On the contact stiffness and nonlinear vibration of an elastic body with a rough surface in contact with a rigid flat surface. Eur J Mech Solid 2015; 49: 321–328.

30.

Chen

Liu

et al . Research on fractal model of normal contact stiffness between two spheroidal joint surfaces considering friction factor. Tribol Int 2016; 97: 253–264.

31.

Miao

Huang

A complete contact model of a fractal rough surface. Wear 2014; 309: 146–151.

32.

Wan

Gong

et al . An improved algorithm for the normal contact stiffness and damping of a mechanical joint surface. Proc IMechE, Part B: J Engineering Manufacture 2014; 228: 751–765.

33.

Gong

Chen

et al . Coupling fractal model for adhesive and three-body abrasive wear of AISI 1045 carbon steel spool valves. Wear 2019; 418–419: 75–85.

34.

Jackson

Duvvuru

Meghani

et al . An analysis of elasto-plastic sliding spherical asperity interaction. Wear 2007; 262: 210–219.

35.

Sepehri

Farhang

On elastic interaction of nominally flat rough surfaces. J Tribol 2008; 130: 1–5.

36.

Zhao

Zhang

Man

et al . A modified normal contact stiffness model considering effect of surface topography. Proc IMechE, Part J: J Engineering Tribology 2015; 229: 677–688.

37.

Zhao

Zhang

Wang

et al . Loading-unloading normal stiffness model for power-law hardening surfaces considering actual surface topography. Tribol Int 2015; 90: 332–342.

38.

Zhao

Zhang

Keer

LM.

Semi-analytical and numerical analysis of sliding asperity interaction for power-law hardening materials. Wear 2016; 364–365: 184–192.

39.

Persson

BNJ

. Contact mechanics for randomly rough surfaces. Surf Sci Rep 2006; 61: 201.

40.

Zhang

Jackson

RL.

An analysis of generated fractal and measured rough surfaces in regards to their multi-scale structure and fractal dimension. Tribol Int 2017; 105: 94–101.

41.

Majumdar

Bhushan

Fractal model of elastic plastic contact between rough surfaces. J Tribol 1991; 113: 1–11.