Three-dimensional fractal model of normal contact damping of dry-friction rough surface

Abstract

According to the influence of the normal contact damping of joint surfaces on the dynamic characteristics of high-precision machinery (machine tool, robot, etc.), in this article, a three-dimensional fractal model of normal contact damping of dry-friction rough joint surfaces based on Hertz theory and fractal theory is established. The three-dimensional surface topography is constructed, according to the modified double variable Weierstrass–Mandelbrot function. The fractal model of strain energy $E_{\mathrm{e}}$ , dissipated energy $E_{\mathrm{p}}$ , damping loss factor $\mathrm{\Delta }$ , and normal contact damping $C_{\mathrm{n}}$ are deduced in detail, and the influence of fractal parameters and dynamic friction coefficient $\mu$ on them are simulated. The simulation curves show that strain energy $E_{\mathrm{e}}$ , dissipated energy $E_{\mathrm{p}}$ , damping loss factor $\mathrm{\Delta }$ , and normal damping $C_{\mathrm{n}}$ increase with the increase in the fractal roughness G; the influence of fractal dimension D on $E_{\mathrm{e}}$ is more changeable, first $E_{\mathrm{e}}$ decreases with the increase in D and then increases. $E_{\mathrm{p}}$ , $\mathrm{\Delta }$ , and $C_{\mathrm{n}}$ increase with the increase in D and $\mu$ , respectively; the effect of $\mu$ on $E_{\mathrm{e}}$ is not obvious, so simple change in $\mu$ has no significant change in $E_{\mathrm{e}}$ ; a comparative analysis of the theoretical calculation of normal damping and experimental results show that their general trend is consistent, and they increase with the increase in total normal contact load P, the relative error is 5%–25%. The theoretical model can provide reference for the design of normal contact damping of the joint surfaces.

Introduction

Mechanical joint surfaces is widely existed in all kinds of mechanical systems, such as the joint surfaces of spindle and guide rail of high-end CNC machine tool, joint surface of robot’s moving pairs and rotating pairs, bolted joint surfaces, and so on. The identification and design of the contact stiffness and damping parameters of the mechanical joint surfaces have important influence on the pre-judgment and evaluation of the static and dynamic characteristics of the whole machine. According to the research, ¹ >90% of the total damping of the machine tool are contributed by joint surfaces, that is, the contact damping on the joint surfaces is much higher than the structure damping of material. Therefore, it is necessary to carefully make the modeling of contact damping on the joint surfaces, the accurate damping model not only can be used to analyze the contact mechanical behavior of joint surfaces but also provide a basis for the construction of the dynamics model of the whole machine’s, so as to better analyze the nonlinearity of system dynamics which is caused by the nonlinearity of the joint surfaces.^2,3

From the macroscopic view, the surfaces of mechanical components are smooth, while from the microscopic view, there are a large number of rough asperities on the surfaces. The analysis of contact behavior on joint surfaces is the study of contact problem of asperities on joint surfaces in essence. The contact between two elastic solids was first studied by Hertz in 1882, the equation of elastic contact mechanics ⁴ was given, which is a milestone in the study of solid contact. However, the elastic–plastic and plastic behavior of the contact of two solids and the friction factor are not considered in the Hertz theory. Thus, Zhao et al. ⁵ considered the whole process of elastic to complete plastic flow when asperities are loaded and established a micro contact model. But the elastic–plastic zone is not further divided. Kogut and Etsion^6–8 in detail researched the contact process of sphere and rigid plate under external loading by finite element method; the dimensionless deformation $\omega$ was given; the elastic, elastic–plastic, and plastic deformation zones of the sphere were carefully divided; and the model of loading and static friction was studied. Shi and Polycarpou ⁹ gave the equation of micro-scale contact stiffness and damping and compared them with the Hertz solution. Majumdar and Bhushan ¹⁰ found that there was statistical self-similarity on each scale through the observation of the surface topography of components, thereupon the classic fractal contact model (M-B model) was given first based on fractal theory. ¹¹ Tian et al. ¹² and Li et al.^13,14 deduced the contact damping and contact stiffness model based on the M-B model. Since the contact damping model is based on two-dimensional fractal, it is not reasonable^15,16 in theory to describe the three-dimensional surface topography based on two-dimensional fractal curve, and the elastic–plastic transition mechanism of asperities is not considered.

Therefore, this article is based on the literatures above, by considering the influence of dynamic friction coefficient, elastic–plastic deformation mechanism, and three-dimensional fractal dimension, the normal contact damping fractal model of dry-friction joint surfaces based on Hertz theory and fractal theory is derived in detail. The influence of the fractal parameters of joint surfaces on the damping model is simulated, which provides a basis for the design and modeling of contact damping of joint surfaces.

Basis of contact theory

Elastic deformation and plastic deformation of asperity

According to the study of three-dimensional rough topography of engineering surfaces, Yan and Komvopoulos modified the Weierstrass–Mandelbrot (W-M) function. The W-M function is a fractal function, the feature of the function is that any point is amplified, and the obtained regional images and overall image are similar. Moreover, in any case, the images of the function are not more smooth, and there are also no monotonic interval. In other words, the function is continuous and nonderivable everywhere. The modified W-M function’s expression ¹⁷ is

$$z\left(x,y\right)=L{\left(\frac{G}{L}\right)}^{(D-2)}{\left(\frac{\ln \gamma }{M}\right)}^{1/2}{\displaystyle \sum _{m=1}^M{\displaystyle \sum _{n=0}^{n_{\max }}{\gamma }^{(D-3)n}}}\times \left\{\cos {\phi }_{m,n}-\cos \left[\frac{2\pi {\gamma }^n{\left(x^2+y^2\right)}^{1/2}}{L}\times \cos \left({\tan }^{-1}\left(\frac{y}{x}\right)-\frac{\pi m}{M}\right)+{\phi }_{m,n}\right]\right\}$$

(1)

where L is the sampling length of surface topography; D is the three-dimensional fractal dimension ( $2<D<3$ ); G is fractal roughness; $\gamma$ ( $\gamma >1$ ) is a parameter that determines the density of frequency; M is the number of surface peak ridge superposition; n denotes the frequency index, $n_{max}=\mathrm{Int}[\log (L/L_{\mathrm{s}})/\log \gamma ]$ ; $L_{\mathrm{s}}$ is the minimum truncation length; and ${\phi }_{m,n}$ is the random phase. The contact surface topography has fractal characteristics, which mainly depends on two important fractal parameters, fractal dimension D and fractal roughness G. That the larger D means more high-frequency components of three-dimensional surface topography, so the surface is much finer in the macro view; that the larger G means amplitude of three-dimensional surface topography is larger, so the surface is more rough. These given parameters $D=2.75,G=5.42\times {10}^{-5}\mathrm{m}$ , $L=256,L_{\mathrm{s}}=1\mathrm{\mu }\mathrm{m},\gamma =1.5,\mathrm{and}M=5$ can be simulated to obtain the surface topography of Figure 1.

Figure 1.

Three-dimensional surface topography.

In order to determine the deformation of asperity, equation (1) is simplified as equation (2) by Yan et al.

$$\begin{array}{l}z\left(x\right)=G^{D-2}{\left(\ln \gamma \right)}^{0.5}{\left(2r\prime \right)}^{3-D}\\ \left[cos{\phi }_{1,n_0}-\cos \left(\frac{\pi x}{r\prime }-{\phi }_{1,n_0}\right)\right]\end{array}$$

(2)

Equation (2) is a cosine function, the height difference between the peak and valley of the function is defined as deformation amount $\omega$ of asperity after loading

$$\omega =2G^{D-2}{\left(\ln \gamma \right)}^{0.5}{(2r\prime )}^{3-D}$$

(3)

where $r\prime$ is the radius of the equivalent asperity’s cross section, which is cut by the equivalent rigid plane, in Figure 2. According to Hertz theory, the contact of the two hemispherical asperity after loading can be equivalent to the contact of a rigid plane and an equivalent asperity, as shown in Figure 2.

Figure 2.

Sketch of contact equivalent of asperity.

In Figure 2, the right-angled triangle ode exists in the following relationship by the Pythagorean theorem

$$R^2=r{\prime }^2+{(R-\omega )}^2$$

(4)

Equation (4) can be transformed into

$$R=\frac{\omega }{2}+\frac{{r\prime }^2}{2\omega }$$

(5)

As the R is much larger than the loaded deformation $\omega$ , so there is $R>>\omega /2$ , the approximate expression of equation (5) is

$$R\approx \frac{{r\prime }^2}{2\omega }$$

(6)

According to the classical Hertz elastic contact theory, the normal load of a single asperity in the elastic deformation zone is

$$f_{\mathrm{e}}\left(\omega \right)=\frac{4}{3}E{R}^{0.5}{\omega }^{1.5}$$

(7)

When the asperity is Hertz elastic contact, the radius ¹⁷ r of the actual contact area is

$$r={\left(\frac{3f_{e}R}{4E}\right)}^{1/3}$$

(8)

where E is the equivalent elastic modulus of asperity, it can be expressed as

$$E={(\frac{1-{\nu }_1^2}{E_1}+\frac{1-{\nu }_2^2}{E_2})}^{-1}$$

(9)

where E₁ and E₂ are the elastic modulus of two contact asperities, respectively, and ${\nu }_1$ and ${\nu }_2$ represent their Poisson’s ratio, respectively.

Substituting equation (7) into equation (8) can obtain

$$r={\left(R\omega \right)}^{0.5}$$

(10)

Substituting equation (6) into equation (10) can obtain

$$r=\frac{1}{\sqrt{2}}r\prime$$

(11)

According to equation (11), in Figure 2, the actual contact area of the equivalent asperity and the equivalent rigid plane is

$$a=\mathrm{\pi }{r}^2=\frac{1}{2}\mathrm{\pi }r{\prime }^2=\frac{1}{2}\mathrm{\pi }[R^2-{(R-\omega )}^2]\approx \mathrm{\pi }R\omega$$

(12)

The elastic critical deformation amount of a single asperity is ¹³

$${\omega }_{\mathrm{c}}={\left(\frac{33\mathrm{\pi }{k}_{\mu }\varphi }{40}\right)}^{2}R$$

(13)

where $\varphi ={\sigma }_{\mathrm{y}}/E$ is the material’s strain parameter; ${\sigma }_{\mathrm{y}}$ is the yield strength of the softer material; and $k_{\mu }$ is the piecewise factor of friction, and its expression is

$$k_{\mu }=\{\begin{array}{c}1-0.228\mu 0\le \mu \le 0.3\hfill \\ \multicolumn{1}{c}{0.932e^{-1.58(\mu -0.3)}0.3<\mu \le 0.94}\end{array}$$

(14)

Substituting equation (12) into equation (3) and eliminating $r\prime$ can obtain

$$\omega =2^{5.5-1.5D}{G}^{D-2}{\left(\ln \gamma \right)}^{0.5}{\mathrm{\pi }}^{0.5D-1.5}{a}^{1.5-0.5D}$$

(15)

Substituting equation (15) into equation (12) can obtain the curvature radius, which is expressed by the area a

$$R=2^{1.5D-5.5}{\mathrm{\pi }}^{0.5-0.5D}{G}^{2-D}{a}^{0.5D-0.5}{\left(\ln \gamma \right)}^{-0.5}$$

(16)

According to equations (13), (15), and (16), the elastic critical deformation area of a single asperity is

$$\begin{array}{l}{a}_{\mathrm{c}}=2^{(3D-11)/(2-D)}{\left(\frac{33k_{\mu }\varphi }{40}\right)}^{2/(2-D)}\\ {\pi }^{(4-D)/(2-D)}{\left(\ln \gamma \right)}^{1/(D-2)}{G}^2\end{array}$$

(17)

The normal load of a single asperity at fully plastic stage is

$$f_{\mathrm{p}}=\lambda {\sigma }_{\mathrm{y}}a$$

(18)

where $\lambda =H/{\sigma }_{\mathrm{y}}$ is the defined parameter and H is the hardness of softer material.

In order to properly associate the elastic deformation of the asperity with the plastic deformation of the asperity, the elastic contact area and the plastic contact area of the asperity can be expressed by a in the classical M-B model. ¹⁰ Furthermore, the value of the contact area a is the same in the M-B model, for the sake of simplicity, this article refers to the representation in the M-B model. So, we can substitute equation (12) into equation (18), and the relationship between the plastic load and deformation is

$$f_{\mathrm{p}}\left(\omega \right)=\lambda {\sigma }_{\mathrm{y}}\mathrm{\pi }R\omega$$

(19)

Elastic–plastic deformation of asperity

According to previous studies,^6,15,16 the difference between the two elastoplastic zones is deformation value, and according to the size of the deformation, the elastic–plastic zone of asperity can be further divided into elastic–plastic I zone and elastic–plastic II zone, and the critical deformation area of the two deformation zones are

$$a_{\mathrm{ep}1}={110}^{1/(2-D)}{a}_{\mathrm{c}}$$

(20)

$$a_{\mathrm{ep}2}=6^{1/(2-D)}{a}_{\mathrm{c}}$$

(21)

Similarly, the relationship between the contact load and deformation on the two elastic–plastic zones is as follows:

1. When $1\le \omega /{\omega }_{\mathrm{c}}\le 6$ , the relationship between the normal load f_ep1 of elastic–plastic I zone and deformation is

   $$\frac{f_{\mathrm{ep}1}}{f_{\mathrm{c}}}=1.03{\left(\frac{\omega }{{\omega }_{\mathrm{c}}}\right)}^{1.425}$$

   (22)

2. When $6\le \omega /{\omega }_{\mathrm{c}}\le 110$ , the relationship between the normal load f_ep2 of elastic–plastic II zone and deformation is

   $$\frac{f_{\mathrm{ep}2}}{f_{\mathrm{c}}}=1.40{\left(\frac{\omega }{{\omega }_{\mathrm{c}}}\right)}^{1.263}$$

   (23)

where $f_{\mathrm{c}}$ is the normal load, when $\omega ={\omega }_{\mathrm{c}}$ , that is, $f_{\mathrm{c}}=(4/3)E{R}^{0.5}{{\omega }_{\mathrm{c}}}^{1.5}$ .

Substituting equation (13) into equations (22) and (23), respectively, which can be simplified as

$$\begin{array}{l}{f}_{\mathrm{ep}1}\left(\omega \right)=\frac{1.03}{3}E{\left(\frac{33k_{\mu }\varphi }{40}\right)}^{0.15}{\left(\ln \gamma \right)}^{0.425}{2}^{6.75-1.275D}\\ \times G^{0.85D-1.7}{\pi }^{0.15}{R}^{1.85-0.425D}{\omega }^{1.85-0.425D}\end{array}$$

(24)

$$\begin{array}{l}{f}_{\mathrm{ep}2}\left(\omega \right)=\frac{1.4}{3}E{\left(\frac{33k_{\mu }\varphi }{40}\right)}^{0.474}{\left(\ln \gamma \right)}^{0.263}{2}^{4.893-0.789D}\\ \times G^{0.526D-1.052}{\pi }^{0.474}{R}^{1.526-0.263D}{\omega }^{1.526-0.263D}\end{array}$$

(25)

Damping loss factor of joint surfaces

When the external force acts on the asperity of rough metal surface, the elastoplastic transition deformation includes the elastic and plastic deformation, and they will appear at the same time, regardless of stage. So that when the strain energy and dissipated energy of a single asperity are established, a part of energy brought by the elastoplastic deformation needs to be taken into account.

The strain energy is obtained from elastic and elastic–plastic deformation stage of a single asperity, and it can be deduced by integrating of each stage

$$\begin{array}{c}\multicolumn{1}{c}{U_{\mathrm{e}}=\underset{0}{\overset{\omega }{\int }}[f_{\mathrm{ep}2}\left(\omega \right)+f_{\mathrm{ep}1}\left(\omega \right)+f_{\mathrm{e}}\left(\omega \right)]\mathrm{d}\omega }\\ =\frac{1.4}{3(2.526-0.263D)}E{\left(\frac{33k_{\mu }\varphi }{40}\right)}^{0.474}{\left(\ln \gamma \right)}^{0.263}{2}^{4.893-0.789D}\hfill \\ \times G^{0.526D-1.052}{\mathrm{\pi }}^{0.474}{R}^{1.526-0.263D}{\omega }^{2.526-0.263D}\hfill \\ +\frac{1.03}{3(2.85-0.425D)}E{\left(\frac{33k_{\mu }\varphi }{40}\right)}^{0.15}{\left(\ln \gamma \right)}^{0.425}{2}^{6.75-1.275D}\hfill \\ \times G^{0.85D-1.7}{\mathrm{\pi }}^{0.15}{R}^{1.85-0.425D}{\omega }^{2.85-0.425D}+\frac{8}{15}E{R}^{0.5}{\omega }^{2.5}\hfill \end{array}$$

(26)

The dissipated energy is obtained from elastic–plastic and plastic deformation stage of a single asperity, and it can be deduced by integrating of each stage

$$\begin{array}{c}\multicolumn{1}{c}{U_{\mathrm{p}}=\underset{0}{\overset{\omega }{\int }}[f_{\mathrm{p}}\left(\omega \right)+f_{\mathrm{ep}2}\left(\omega \right)+f_{\mathrm{ep}1}\left(\omega \right)]\mathrm{d}\omega }\\ =\frac{1}{2}\lambda {\sigma }_{\mathrm{y}}\mathrm{\pi }R{\omega }^2+\frac{1.4}{3(2.526-0.263D)}E{\left(\frac{33k_{\mu }\varphi }{40}\right)}^{0.474}\hfill \\ \times {\left(\ln \gamma \right)}^{0.263}{2}^{4.893-0.789D}{G}^{0.526D-1.052}{\mathrm{\pi }}^{0.474}{R}^{1.526-0.263D}{\omega }^{2.526-0.263D}\hfill \\ +\frac{1.03}{3(2.85-0.425D)}E{\left(\frac{33k_{\mu }\varphi }{40}\right)}^{0.15}{\left(\ln \gamma \right)}^{0.425}{2}^{6.75-1.275D}{G}^{0.85D-1.7}\hfill \\ \times {\mathrm{\pi }}^{0.15}{R}^{1.85-0.425D}{\omega }^{2.85-0.425D}\hfill \end{array}$$

(27)

The area distribution function $n(a)$ of asperity on the joint surfaces is ¹⁷

$$n\left(a\right)=\frac{D-1}{2}{a}_{\mathrm{l}}^{(D-1)/2}{a}^{-(D+1)/2},0<a\le a_{\mathrm{l}},2<D<3$$

(28)

where $a_{\mathrm{l}}$ is the maximum contact area of asperity.

According to the area distribution function $n(a)$ , the strain energy of the whole joint surfaces in a contact loading process is deduced by the integration of the strain energy which is obtained from elastic and elastic–plastic deformation stages of all asperities

$$\begin{array}{c}\multicolumn{1}{c}{E_{\mathrm{e}}=\underset{a_{\mathrm{ep}2}}{\overset{a_{\mathrm{l}}}{\int }}{U}_{\mathrm{e}}n\left(a\right)\mathrm{d}a}\\ =\frac{1.4E(D-1)}{3(2.526-0.263D)(3.526-1.263D)}{\left(\frac{33k_{\mu }\varphi }{40}\right)}^{0.474}{\left(\ln \gamma \right)}^{0.763}\hfill \\ \times 2^{9.393-2.289D}{G}^{1.526D-3.052}{\mathrm{\pi }}^{-2.522+0.763D}{a}_{\mathrm{l}}^{0.5D-0.5}(a_{\mathrm{l}}^{3.526-1.263D}-a_{\mathrm{ep}2}^{3.526-1.263D})\hfill \\ +\frac{1.03E(D-1)}{3(2.85-0.425D)(3.85-1.425D)}{\left(\frac{33k_{\mu }\varphi }{40}\right)}^{0.15}\hfill \\ \times {\left(\ln \gamma \right)}^{0.925}{2}^{11.25-2.775D}{G}^{1.85D-3.7}{\mathrm{\pi }}^{-3.2+0.925D}{a}_{\mathrm{l}}^{0.5D-0.5}(a_{\mathrm{l}}^{3.85-1.425D}-a_{\mathrm{ep}2}^{3.85-1.425D})\hfill & +\frac{D-1}{3(2.5-D)}E{\mathrm{\pi }}^{-2.+0.5D}{2}^{6.5-1.5D}\hfill \\ \times {\left(\ln \gamma \right)}^{0.5}{G}^{D-1}{a}_{\mathrm{l}}^{0.5D-0.5}(a_{\mathrm{l}}^{2.5-D}-a_{\mathrm{ep}2}^{2.5-D})\hfill \end{array}$$

(29)

Similarly, the dissipated energy of the whole joint surfaces in a contact loading process is deduced by the integration of the dissipated energy which is obtained from elastic–plastic and plastic deformation stages of all asperities

$$\begin{array}{c}\multicolumn{1}{c}{E_{\mathrm{p}}=\underset{0}{\overset{a_{\mathrm{c}}}{\int }}{U}_{\mathrm{p}}n\left(a\right)\mathrm{d}a}\\ =\frac{D-1}{3-D}\lambda {\sigma }_{\mathrm{y}}{2}^{3.5-1.5D}{\mathrm{\pi }}^{0.5D-1.5}{G}^{D-2}{\left(\ln \gamma \right)}^{0.5}{a}_{\mathrm{l}}^{0.5D-0.5}{a}_{\mathrm{c}}^{3-D}\hfill \\ +\frac{1.4E(D-1)}{3(2.526-0.263D)(3.526-1.263D)}{\left(\frac{33k_{\mu }\varphi }{40}\right)}^{0.474}{\left(\ln \gamma \right)}^{0.763}\hfill \\ \times 2^{9.393-2.289D}{G}^{1.526D-3.052}{\mathrm{\pi }}^{-2.522+0.763D}{a}_{\mathrm{l}}^{0.5D-0.5}{a}_{\mathrm{c}}^{3.526-1.263D}\hfill \\ +\frac{1.03E(D-1)}{3(2.85-0.425D)(3.85-1.425D)}{\left(\frac{33k_{\mu }\varphi }{40}\right)}^{0.15}\hfill \\ \times {\left(\ln \gamma \right)}^{0.925}{2}^{11.25-2.775D}{G}^{1.85D-3.7}{\mathrm{\pi }}^{-3.2+0.925D}{a}_{\mathrm{l}}^{0.5D-0.5}{a}_{\mathrm{c}}^{3.85-1.425D}\hfill \end{array}$$

(30)

The dimensionless forms of equations (29) and (30) are, respectively

$$E_{\mathrm{e}}^{*}=\frac{E_{\mathrm{e}}}{{\mathrm{EA}}_{\mathrm{a}}^{1.5}}$$

(31)

$$E_{\mathrm{p}}^{*}=\frac{E_{\mathrm{p}}}{{\mathrm{EA}}_{\mathrm{a}}^{1.5}}$$

(32)

The damping loss factor $\mathrm{\Delta }$ is a parameter, which is used to describe the hysteresis characteristic of material, and it can be expressed by the ratio of dissipated energy of the whole joint surfaces and strain energy after a single load

$$\mathrm{\Delta }=\frac{E_{\mathrm{p}}}{E_{\mathrm{e}}}$$

(33)

Total normal contact load

According to the continuous integration of the area distribution function on each deformation stage of asperity, the actual contact area can be calculated

$$A_{\mathrm{r}}=\underset{0}{\overset{a_{\mathrm{l}}}{\int }}n\left(a\right)a\mathrm{d}a==\frac{D-1}{3-D}{a}_{\mathrm{l}}$$

(34)

According to equation (7), equations (15)–(18), (20), (21), (24), and (25), multiplying the load by area distribution function on each stage and on each stage integrating on deformation area can obtain the total normal contact load P:

1. When D ≠ 2.5, there is

   $$\begin{array}{c}\multicolumn{1}{c}{P=\underset{0}{\overset{a_{\mathrm{ep}1}}{\int }}{f}_{\mathrm{p}}n\left(a\right)\mathrm{d}a+\underset{a_{\mathrm{ep}1}}{\overset{a_{\mathrm{ep}2}}{\int }}{f}_{\mathrm{ep}2}n\left(a\right)\mathrm{d}a+\underset{a_{\mathrm{ep}2}}{\overset{a_{\mathrm{c}}}{\int }}{f}_{\mathrm{ep}1}n\left(a\right)\mathrm{d}a+\underset{a_{\mathrm{c}}}{\overset{a_{\mathrm{l}}}{\int }}{f}_{\mathrm{e}}n\left(a\right)\mathrm{d}a}\\ =\frac{D-1}{3-D}\lambda {\sigma }_{\mathrm{y}}{110}^{(3-D)/(4-2D)}{a}_{\mathrm{l}}^{0.5D-0.5}{a}_{\mathrm{c}}^{1.5-0.5D}\hfill \\ +\frac{1.40}{3}E{\left(\frac{33k_{\mu }\varphi }{40}\right)}^{0.474}\frac{D-1}{2.026-0.763D}{2}^{3.893-0.789D}{\mathrm{\pi }}^{0.263D-1.052}{\left(\ln \gamma \right)}^{0.263}\hfill \\ \times G^{0.526D-1.052}{a}_{\mathrm{l}}^{0.5D-0.5}{a}_{\mathrm{c}}^{2.026-0.763D}[6^{(2.026-0.763D)/(2-D)}-{110}^{(2.026-0.763D)/(2-D)}]\hfill \\ +\frac{1.03}{3}E{\left(\frac{33k_{\mu }\varphi }{40}\right)}^{0.15}\frac{D-1}{2.35-0.925D}{2}^{5.75-1.275D}{\mathrm{\pi }}^{0.425D-1.7}\hfill \\ \times {\left(\ln \gamma \right)}^{0.425}{G}^{0.85D-1.7}{a}_{\mathrm{l}}^{0.5D-0.5}{a}_{\mathrm{c}}^{2.35-0.925D}[1-6^{(2.35-0.925D)/(2-D)}]\hfill \\ +\frac{(D-1)}{3(2.5-D)}{2}^{6.5-1.5D}E{\mathrm{\pi }}^{0.5D-2}{\left(\ln \gamma \right)}^{0.5}{G}^{D-2}{a}_{\mathrm{l}}^{0.5D-0.5}(a_{\mathrm{l}}^{2.5-D}-a_{\mathrm{c}}^{2.5-D})\hfill \end{array}$$

   (35)

2. When D = 2.5, there is

   $$\begin{array}{c}\multicolumn{1}{c}{P=\underset{0}{\overset{a_{\mathrm{ep}1}}{\int }}{f}_{\mathrm{p}}n\left(a\right)\mathrm{d}a+\underset{a_{\mathrm{ep}1}}{\overset{a_{\mathrm{ep}2}}{\int }}{f}_{\mathrm{ep}2}n\left(a\right)\mathrm{d}a+\underset{a_{\mathrm{ep}2}}{\overset{a_{\mathrm{c}}}{\int }}{f}_{\mathrm{ep}1}n\left(a\right)\mathrm{d}a+\underset{a_{\mathrm{c}}}{\overset{a_{\mathrm{l}}}{\int }}{f}_{\mathrm{e}}n\left(a\right)\mathrm{d}a}\\ =3\lambda {\sigma }_{\mathrm{y}}{110}^{-0.5}{a}_{\mathrm{l}}^{0.75}{a}_{\mathrm{c}}^{0.25}+0.7E{\left(\frac{33k_{\mu }\varphi }{40}\right)}^{0.474}{2}^{1.9205}{\mathrm{\pi }}^{-0.3945}{\left(\ln \gamma \right)}^{0.263}\hfill \\ \times G^{0.263}{a}_{\mathrm{l}}^{0.75}{a}_{\mathrm{c}}^{0.1185}(6^{-0.237}-{110}^{-0.237})\hfill \\ +0.515E{\left(\frac{33k_{\mu }\varphi }{40}\right)}^{0.15}{2}^{2.5625}{\mathrm{\pi }}^{-0.6375}{\left(\ln \gamma \right)}^{0.425}{G}^{0.425}{a}_{\mathrm{l}}^{0.75}{a}_{\mathrm{c}}^{0.0375}(1-6^{-0.075})\hfill \\ +0.125E{\mathrm{\pi }}^{-0.75}{2}^{4.75}{\left(\ln \gamma \right)}^{0.5}{G}^{0.5}{a}_{\mathrm{l}}^{0.75}(\ln a_{\mathrm{l}}-\ln a_{\mathrm{c}})\hfill \end{array}$$

   (36)

The dimensionless forms of equations (35) and (36) are

$$P^{*}=\frac{P}{E{A}_{\mathrm{a}}}$$

(37)

Fractal prediction model of normal contact damping

According to the vibration theory, the relationship between the damping loss factor $\mathrm{\Delta }$ and damping ratio $\xi$ is

$$\mathrm{\Delta }=2\xi =2\frac{C_{\mathrm{n}}}{C_{\mathrm{c}}}$$

(38)

where $C_{\mathrm{n}}$ is the normal damping of joint surfaces and $C_{\mathrm{c}}$ is the critical damping coefficient. $C_{\mathrm{c}}$ is expressed as

$$C_{\mathrm{c}}=2{\left(m{K}_{\mathrm{n}}\right)}^{0.5}$$

(39)

where m is the mass of the object and K_n is the normal stiffness of joint surfaces. ¹³

So, the normal damping fractal model of joint surfaces is

$$C_{\mathrm{n}}=\mathrm{\Delta }{\left(m{K}_{\mathrm{n}}\right)}^{0.5}$$

(40)

The dimensionless form of equation (40) is

$$C_{\mathrm{n}}^{*}=\mathrm{\Delta }{K}_{\mathrm{n}}^{*0.5}$$

(41)

where $C_{\mathrm{n}}^{*}=C_{\mathrm{n}}/(A_{\mathrm{a}}^{0.25}(\mathrm{mE}{)}^{0.5}),K_{\mathrm{n}}^{*}=K_{\mathrm{n}}/({\mathrm{EA}}_{\mathrm{a}}^{0.5})$ .

Numerical simulation analysis and experimental verification

The effect of fractal parameters and dynamic friction coefficient on strain energy

As is shown in Figure 3, the dimensionless strain energy $E_{\mathrm{e}}^{*}$ increases with the increase in dimensionless actual contact area $A_{\mathrm{r}}^{*}$ . In Figure 3(a) and (b), the dimensionless strain energy decreases with the increase in fractal dimension D and then increases with the increase in D. In Figure 3(c), the dimensionless strain energy decreases with the decrease in dimensionless fractal roughness G*. In Figure 3(d), the significant variation of the dynamic friction coefficient $\mu$ has a little influence on the change of dimensionless strain energy. On the whole, the influence of D and G* on $E_{\mathrm{e}}^{*}$ is great, where the influence of D on $E_{\mathrm{e}}^{*}$ is not always monotone increasing or decreasing in the range of the whole dimensionless actual contact area, and it is increasing at first and then decreasing; $\mu$ have little effect on $E_{\mathrm{e}}^{*}$ , that is, simple changing of $\mu$ of joint surfaces is invalid with significant changes in $E_{\mathrm{e}}^{*}$ .

Figure 3.

Effect of D, G*, and µ on $E_{\mathrm{e}}^{*}$ : (a) D = 2.15–2.45, (b) D = 2.55–2.85, (c) G* = 10⁻⁶−10⁻¹⁰, and (d) µ = 0.1–0.9.

The effect of fractal parameters and dynamic friction coefficient on dissipated energy

As is shown in Figure 4(a) and (b), the dimensionless dissipated energy $E_{\mathrm{p}}^{*}$ increases with the increase in fractal dimension D and dimensionless fractal roughness G*, but the rate of change of $E_{\mathrm{p}}^{*}$ following with the D and G* is greater than the rate of change of $E_{\mathrm{e}}^{*}$ in Figure 3. In Figure 4(c), $E_{\mathrm{p}}^{*}$ is compared to $E_{\mathrm{e}}^{*}$ , and $E_{\mathrm{p}}^{*}$ increases with the increase in $\mu$ . On the whole, the three parameters have significant effect on $E_{\mathrm{p}}^{*}$ .

Figure 4.

Effect of D, G*, and µ on $E_{\mathrm{p}}^{*}$ : (a) D = 2.05–2.17, (b) G* = 10⁻⁶−10⁻¹⁰, and (c) µ = 0.1–0.9.

The effect of fractal parameters and dynamic friction coefficient on damping loss factor

As is shown in Figure 5, the damping loss factor $\mathrm{\Delta }$ decreases with the increase in the dimensionless actual contact area $A_{\mathrm{r}}^{*}$ . In a small range of $A_{\mathrm{r}}^{*}$ (0–0.01), $\mathrm{\Delta }$ decreases linearly and then changes gently, that is, $A_{\mathrm{r}}^{*}$ has significant effect on $\mathrm{\Delta }$ in a very small range. In addition, $\mathrm{\Delta }$ increases with the increase of two fractal parameters D and G* and dynamic friction coefficient $\mu$ . Therefore, if we need to increase the damping loss factor of joint surfaces, we could increase fractal dimension, fractal roughness, and dynamic friction coefficient.

Figure 5.

Effect of D, G*, and µ on Δ: (a) D = 2.510–2.516, (b) G* = 1.15 × 10⁻¹⁰−1.95 × 10⁻¹⁰, and (c) µ = 0.1–0.9.

The effect of fractal parameters and dynamic friction coefficient on the normal contact damping

As is shown in Figure 6, the dimensionless normal contact damping of the joint surfaces increases with the increase in the dimensionless total normal load. The reason is that the total normal contact load increases and then the normal contact stiffness increases, which leads to the increase in the critical damping coefficient and finally represents the increase in the normal contact damping. In Figure 6(a)–(c), the dimensionless normal damping increases with the increase in fractal dimension, fractal roughness, and dynamic friction coefficient is also shown. Therefore, when the total normal load of the joint surfaces is constant, the change in normal contact damping can be achieved by changing fractal parameters and dynamic friction coefficient.

Figure 6.

Effect of D, G*, and µ on $C_{\mathrm{n}}^{*}$ : (a) D = 2.510–2.516, (b) G* = 1.342 × 10⁻¹⁰−6.710 × 10⁻¹⁰, and (c) µ = 0.1–0.3.

Comparison between theoretical and experimental value

In order to verify the reasonableness of the three-dimensional fractal model of normal contact damping of the dry-friction joint surfaces in this article, the comparison analysis of theoretical value and experimental value in Zhang et al. ¹⁸ is carried out. The engineering parameters needed for theoretical calculation are ${\sigma }_{\mathrm{y}}=235\mathrm{MPa}$ , $E=2.07\times {10}^{11}\mathrm{Pa}$ , $H=1.96\times {10}^9\mathrm{Pa}$ , $\nu =0.29$ , m = 2800 kg, D = 2.4241, and G = 5.1372 × 10⁻⁵ m. When the applied load is P = 2.5 kN, substitute the engineering parameters above into equations (29), (30), and (35) in this article and normal stiffness formula in the literature, ¹³ the damping loss factors, $\mathrm{\Delta }=5.4562\times {10}^{-4}$ and $K_{\mathrm{n}}=4\times {10}^8\mathrm{N}/\mathrm{m}$ , can be obtained. According to equation (40) in this article, when the applied load is P = 2.5 kN, we can find that the normal contact damping is $C_{\mathrm{n}}=18.2589\mathrm{kN}\mathrm{s}{\mathrm{m}}^{-1}$ . In the same way, we can find out the normal contact damping under different applied load P.

The experimental study on contact damping identification of machine tool ground foot was given in Zhang et al., ¹⁸ and the two tangential (x-direction and y-direction) and a normal (z-direction) contact damping test values of joint surfaces were obtained, as is shown in Figure 7. The overall trend of theoretical and experimental normal contact damping curves is consistent, and they both increase with the increase in the total normal contact load P, and the relative error between theoretical value and experimental value is −5% to 25%, <30%. Therefore, the normal contact damping fractal model can be used to predict the normal damping value of joint surfaces under different external loads to a certain extent.

Figure 7.

Comparison between theoretical and experimental damping: (a) comparison between theoretical and experimental values and (b) the relative error between theoretical and experimental normal contact damping.

Conclusion

1. The strain energy $E_{\mathrm{e}}$ and dissipated energy $E_{\mathrm{p}}$ increase with increasing fractal roughness G; $E_{\mathrm{p}}$ also increases with the increase in fractal dimension D, but the effect of D on $E_{\mathrm{e}}$ is more changeable; first $E_{\mathrm{e}}$ decreases with the increase in D and then increases; $\mu$ has certain effect on $E_{\mathrm{p}}$ , $E_{\mathrm{p}}$ increases with the increase in $\mu$ , but the effect of $\mu$ on $E_{\mathrm{e}}$ is not very obvious, so simple changing of $\mu$ is invalid for significant change in $E_{\mathrm{e}}$ .

2. The damping loss factor $\mathrm{\Delta }$ increases with the increase in two fractal parameters D and G and dynamic friction coefficient $\mu$ ; the normal contact damping $C_{\mathrm{n}}$ increases with the increase in $\mu$ , D, G, and the total normal contact load P.

3. The comparison between theoretical calculation of normal damping and experimental results show that their general trend is consistent, and they increase with the increase in total load P, and the relative error is −5% to 25%.