Fractal Model of Contact Thermal Stiffness

Abstract

The continuity, self-similarity, and self-affinity of a microscopic contact surface can be described by the Weierstrass–Mandelbrot (W–M) function in fractal theory. To address the problems that the existing normal contact load fractal model does not take into account the effect of thermal stress and is not applicable to the temperature variation in the joint surface of the giant magnetostrictive ultrasonic vibration systems, a fractal model of thermal–elastic–plastic contact normal load fractal is established based on fractal theory. The model is an extension of the traditional model in terms of basic theory and application scope, and it takes into account the effects of temperature difference, linear expansion coefficient, fractal dimension, and other parameters. Finally, the effect of the temperature difference at the joint surface on the normal load of the thermoelastic contact is revealed through numerical simulations. The results show that the nonlinearity of the contact stiffness of the thermoelastic joint surface is mainly related to the surface roughness and the fractal dimension, while the effect of the temperature change on the joint surface properties within a certain range is linear.

1. Introduction

The giant magnetostrictive transducer has the advantages of higher energy density, lower sound speed, no overheating failure, fast response speed, and a strong loading force. These characteristics determine the research direction of developing an ultrasonic machining system with high power, large-amplitude, and wide frequency band [1,2]. Rotational ultrasonic vibration processing can be carried out by giant magnetostrictive transducers made from materials with coupled magnetic, electric, elastic, mechanical, and thermal multiphysics fields [3,4,5,6]. When the GMT works in a high-frequency magnetic field, the eddy current loss in the material is tremendous. The existence of an eddy current weakens the magnetic field provided by the drive coil and changes the uniform distribution of the magnetic field. In addition, frequency-related hysteresis and eddy current reduce the energy conversion efficiency of the giant magnetostrictive materials and change the temperature distribution of the transducer. Temperature rise also causes changes in the contact stiffness of each transducer interface, changes the resonant frequency, and decreases the amplitude stability [2,7]. Zhou et al. [2] developed an equivalent dynamics model of a rotating ultrasonic transducer to investigate the effect of the joint surface stiffness between the horn and the preload force mechanism on the amplitude by a nonlinear least-square fitting method. The results showed that the joint stiffness has an important influence on the dynamic characteristics of the ultrasonic vibration system, but the influence of the thermal effect of the internal giant magnetostrictive material on the interface contact stiffness was not considered. Therefore, determining the interface contact stiffness of the ultrasonic vibration system under the thermal effect and predicting the dynamic characteristics of the ultrasonic system under the change in interface stiffness are urgent problems that need to be solved.

The fractal characteristics of rough surfaces have objective uniqueness or scale independence, providing all surface morphology information within all scale ranges on fractal surfaces [8]. Based on the Weierstrass–Mandelbrot function and statistics, the fractal model and statistical model of rough surface contact were established by early scholars. The direct relationship between the actual contact area, contact load force, contact stiffness, and surface deformation in the elastic, plastic, and elastic–plastic stages of rough surface asperity have been studied. However, these contributions are based on the static characteristics of pure mechanical joints [8,9,10,11,12,13,14,15,16,17,18,19,20]. Based on the fractal theory of rough surface and measurement experiments, Kogut and Komvopoulosa [21] yield insight into the influence of rough surface morphology, current, and contact force on electrical contact thermal resistance. The results showed that the actual contact area of the interface plays a decisive role in the influence of electrical contact thermal resistance. Wang and Komvopoulos [22] established the rough surface fractal theory of interface temperature distribution. They studied the relationship between the maximum temperature rise in the fractal surface domain and the real contact area, thermal–mechanical properties, and fractal parameters. Shi et al. [23] numerically simulated the maximum contact temperature at the interface and the temperature distribution in the contact area of the plate at different sliding speeds. Due to the friction heat source, they studied the thermal–mechanical contact characteristics between the elastic–plastic asperity and the rigid plate. Song et al. [24] studied the thermo–mechanical contact responses of a rigid sphere sliding on an elastic–plastic sphere with a radius more significant than that of the rigid sphere using the finite element method. They obtained the transient values of the interface contact force, the maximum temperature rise, and the actual contact area. Horovistiz et al. [25] used the fractal method in tribological phenomena to fractalize the elevation map of sliding surfaces. They proposed a robust method to determine the threshold between the microscopic scale and the macroscopic scale of fractal behavior; the results showed that the fractal dimension data depend on the tribological conditions and the position along the sliding surface, which have research significance for the accurate determination of material surface fractal parameters.

However, these studies were based on traditional pure mechanical models and did not consider the effect of thermal stress. Therefore, the traditional mechanical interface model cannot truly reflect the actual situation of thermal–mechanical products, such as GMTs working in a long-term changing thermal environment. A fractal model of thermal–elastic–plastic contact normal stiffness, taking into account the thermal effect of rough surface, is proposed in this article by using the superposition method based on the normal contact stiffness model of an anisotropic material. The influence of each main parameter on the normal stiffness of thermal–elastic–plastic contact of the joint surface is analyzed. Fractal theory is used to study the contact behavior of rough surfaces. Then, the theoretical model of the thermal–elastic–plastic contact problem of the rough surface is established. This thermo–mechanical contact characteristic of the joint surface can be utilized to improve the structural design of thermo–mechanical products such as GMTs.

2. Materials and Methods

2.1. Fractal Characterization of Surface Morphology

The rough surface has the mathematical characteristics of statistical self-radioactivity and the multiscale fractal. The Weierstrass–Mandelbrot (W–M) function with double independent variables can fully characterize the mathematical characteristics of a rough surface, and its expression is [26]:

(1)$\begin{array}{c}z(x,y)=L{(\mathrm{G}/L)}^{\mathrm{D}-1}{(\mathrm{In}(\gamma )/\mathrm{M})}^{0.5}{\displaystyle \sum _{\mathrm{m}=1}^{\mathrm{M}}{\displaystyle \sum _{n=0}^{n_{\max }}{\gamma }^{(\mathrm{D}-2)n}}}\times \hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\left\{\cos {\Phi }_{\mathrm{m},n}-\cos \left[2\pi {\gamma }^n{\left(x^2+y^2\right)}^{0.5}/L\times \cos (\arctan (y/x)-\pi \mathrm{m}/\mathrm{M})+{\Phi }_{\mathrm{m},n}\right]\right\}\\ \times n_{\max }=\mathrm{int}\left[\mathrm{In}\left(L/L_S\right)/\mathrm{In}\gamma \right]\hfill \end{array}$

where z(x) is the contour height of the asperity; x, y is the coordinate displacement distance; L is the sample length; L_s is the minimum cut-off length; Φ_1,n is the random phase; γ is greater than 1 scale parameter, generally γ = 1.5; n is the spatial frequency index; n_max is the maximum frequency index related to the distance cut-off length; M is the number of overlapping uplifts on the structural surface; G is the fractal roughness coefficient determined by frequency; and D (1 < D < 2) is the fractal dimension of the profile.

When M = 1, m = 1, the W–M function that becomes a single independent variable is:

(2)$\begin{array}{c}z(x)=L{(\mathrm{G}/L)}^{\mathrm{D}-1}{(In(\gamma ))}^{0.5}{\displaystyle \sum _{n=0}^{n_{\max }}{\gamma }^{(\mathrm{D}-2)n}}\left[\cos {\Phi }_{1n}-\cos \left(\frac{2\pi {\gamma }^n}{L}-{\Phi }_{1n}\right)\right]\\ (1<\mathrm{D}<2,\gamma >1)\end{array}$

The three-dimensional contact profiler is used to measure the magnetic conductive sheet (MCS) of a GMT (see [7]) with different roughnesses, as shown in Figure 1. The yellow line is the center of the measurement position and is measured at the center. The material is electrical pure iron, the elastic modulus is 195 MPa, Poisson’s ratio is 0.29, the roughness of the measurement sample is 0.356 μm, the minimum sampling distance of the rough surface measuring instrument is 0.25 μm, and the sampling length is 4 mm. The results of D and G can be obtained using the structural function method [27]. The summary data is shown in Table 1.

2.2. Fractal Model of Contact Mechanics of Asperity

Figure 2 illustrates the conceptual approach to the contact between a rigid plane and a rough surface. According to Hertz’ contact theory, the contact between two rough -surface asperities is equivalent to the contact between the rigid plane and spherical asperities. To obtain the contact deformation force under a given size, the contact deformation between spherical asperities and a rigid plane must be first obtained. If the maximum contact area of the asperity is $a_l$, the number of asperities larger than the contact area a can be expressed as N, and the area distribution density of the contact point can be expressed as [27]:

(3)$n(a)=\left|\frac{dN}{da}\right|=\frac{\mathrm{D}}{2}\frac{a_l^{0.5\mathrm{D}}}{a^{(0.5\mathrm{D}+1)}}$

According to the area distribution density function of the contact point in Equation (3), the real contact area of the joint surface is derived by:

(4)$A_{\mathrm{r}}={\displaystyle {\int }_0^{a_l}n}(a)ada=\frac{\mathrm{D}}{2-\mathrm{D}}{a}_l$

The interference of an asperity under the force of the rigid plate $\delta$ is equal to the peak-to-valley amplitude of the z(x) function (see Figure 2). $\delta$ is given as [26]:

(5)$\delta =2^{3-\mathrm{D}}{\mathrm{G}}^{\mathrm{D}-1}{(\ln \gamma )}^{0.5}{\pi }^{0.5\mathrm{D}-1}{a}^{0.5\mathrm{D}-1}$

According to the relationship between the deformation $\delta$ of the asperity and the curvature radius R, considering R≫δ the curvature radius of the asperity is [16]:

(6)$R=\frac{a^{0.5\mathrm{D}}}{2^{4-\mathrm{D}}{\pi }^{0.5\mathrm{D}}{\mathrm{G}}^{\mathrm{D}-1}{(\ln \gamma )}^{0.5}}$

To determine the contact parameters of asperity, the contact state of asperity must be determined. According to the finite element results in the literature [9], the deformation of asperity is divided into elastic deformation ($\delta <{\delta }_c$), elastic–plastic deformation 1 (${\delta }_c⩽\delta ⩽6{\delta }_c$), elastic–plastic deformation 2 ($6{\delta }_c⩽\delta ⩽110{\delta }_c$), and complete plastic deformation ($110{\delta }_c<\delta$). ${\delta }_c$ is the elastic critical deformation. When the deformation of the asperity exceeds the elastic critical deformation, it begins to yield and enters the elastic–plastic deformation region, which can be expressed as:

(7)${\delta }_c={\left(\frac{\pi KH}{2E}\right)}^{2}R$

where H is the hardness of the soft material, K is the hardness coefficient, K = 0.454 + 0.41$\nu$, and $\nu$ is Poisson’s ratio.

When ${\delta }_c=\delta$, the corresponding contact area is the critical contact area of the asperity $a_c$. Then, substituting Equations (3) and (4) into Equation (7) gives:

(8)$a_c={\left[2^{9-2\mathrm{D}}{\mathrm{G}}^{2\mathrm{D}-2}{\pi }^{\mathrm{D}-3}\mathrm{In}\gamma {\left(\frac{E}{KH}\right)}^2\right]}^{1/\mathrm{D}-1}$

The real contact areas arising in the elastic, elasto–plastic (1), elasto–plastic (2), and fully plastic deformation regimes can be determined if the deformation of asperity for these four deformation regimes are available. Then, the fractal theory of contact area can be obtained in different stages.

If the deformation of the asperity is in the range of $a_c<a<a_l$($\delta <{\delta }_c$), the elastic contact load $P_e$ and contact stiffness $K_e$, produced by a sphere with a radius of curvature of R in contact with a rigid plane with an elastic interference, are given based on fractal theory as follows:

(9)$a_e=\frac{a^{\mathrm{D}-1}}{2}$

(10)$p_e=\frac{4}{3}E{\left(\frac{a^{0.5\mathrm{D}}}{2^{4-\mathrm{D}}{\pi }^{0.5\mathrm{D}}{\mathrm{G}}^{\mathrm{D}-1}{(\ln \gamma )}^{0.5}}\right)}^{0.5}\times {\left(2^{3-\mathrm{D}}{\mathrm{G}}^{\mathrm{D}-1}{(\ln \gamma )}^{0.5}{\pi }^{0.5\mathrm{D}-1}{a}^{0.5\mathrm{D}-1}\right)}^{1.5}$

(11)$k_e=\frac{d{p}_e}{d\delta }=\frac{4E(3-\mathrm{D})}{3\sqrt{2\pi }(2-\mathrm{D})}{a}^{0.5}$

If the deformation of the asperity is in the range of $6^{-{(\mathrm{D}-1)}^{-1}}{a}_c<a<a_c$ (${\delta }_c⩽\delta ⩽6{\delta }_c$), the contact mechanics model of elastic–plastic deformation zone 1 is expressed in fractal form as:

(12)$p_{1ep}=\frac{1.03}{3}\times 2^{3.825-0.85\mathrm{D}}{(KH)}^{0.15}{\pi }^{0.425\mathrm{D}-1.275}{\mathrm{G}}^{0.85\mathrm{D}-0.85}{(\mathrm{In}\gamma )}^{0.45}{E}^{0.85}{a}^{1.425-0.425\mathrm{D}}$

(13)$k_{1ep}=\frac{d_{1ep}}{d\delta }=\frac{1.03\times 2^{1.825+0.15\mathrm{D}}{(KH)}^{0.15}{E}^{0.85}{a}^{0.075\mathrm{D}+0.425}}{3{\pi }^{0.075\mathrm{D}+0.275}{\mathrm{G}}^{0.15\mathrm{D}-0.15}{(\mathrm{In}\gamma )}^{0.45}(1.425-0.425\mathrm{D}(2-\mathrm{D}))}$

Moreover, if the deformation of the asperity is in the range of ${110}^{-{(\mathrm{D}-1)}^{-1}}{a}_c<a<6^{-{(\mathrm{D}-1)}^{-1}}{a}_c$ ($6{\delta }_c\le \delta ⩽110{\delta }_c$), the contact mechanics model of elastic–plastic deformation zone 2 is expressed in fractal form as:

(14)$p_{2ep}=\frac{1.4}{3}\times 2^{2.367-0.526\mathrm{D}}{(KH)}^{0.474}{\pi }^{0.263\mathrm{D}-0.789}{\mathrm{G}}^{0.526\mathrm{D}-0.526}{(\mathrm{In}\gamma )}^{0.263}{E}^{0.526}{a}^{1.263-0.263\mathrm{D}}$

(15)$k_{2ep}=\frac{d_{2ep}}{d\delta }=\frac{1.4\times 2^{0.367+0.474\mathrm{D}}{(KH)}^{0.474}{E}^{0.526}{a}^{0.237\mathrm{D}+0.263}}{3{\pi }^{0.237\mathrm{D}-0.211}{\mathrm{G}}^{0.474(\mathrm{D}-1)}{(\mathrm{In}\gamma )}^{0.237}(1.263-0.263\mathrm{D})(2-\mathrm{D})}$

When $110{\delta }_c<\delta$, the contact mode is complete plastic deformation, and the contact mechanics model is expressed in fractal form as:

(16)$a_p=2\pi R\delta =a$

(17)$p_p=2\pi R\delta H=H{a}_p$

(18)$k_p=2\pi RH$

2.3. Thermal–Mechanical Coupling Contact Interface

We now present a modified asperity contact stiffness model that is similar to the contact mechanics model of asperity in Section 2.2, only because the aggregated heat flow propagates between rough contact interfaces, as shown in Figure 3. Unlike the interface contact model in Section 2.2, which assumes that the asperity is subjected to pure static compressive stress, this model considers contact thermal resistance and the influence of thermal stress on contact deformation when heat flow passes through the contact interface. This model especially considers the influence of thermoelastic deformation. A more practical interface stiffness model of ultrasonic vibration is established by coupling the contact interface’s static stiffness with the interface’s thermal stiffness under static pressure, namely the interface composite stiffness.

2.3.1. Thermal Contact Resistance of Rough Surface Interface

When the heat flow generated by a magnetostrictive material passes through the rough contact interface of the MCS and the horn, the contact thermal resistance is generated at the interface under the action of the interface contact pressure. Additionally, the contact force between the asperity and the interface is given in Section 2.2 and Section 2.3. Interface contact thermal resistance mainly includes interface material matrix contact thermal resistance, heat flow shrinkage thermal resistance, and interface gap thermal resistance. The contact model of heat flow propagation between the contact interfaces is assumed as follows: the temperature difference between contact interfaces is the same; the heat transfer mode of fluid medium in the contact interface gap is heat conduction, without thermal convection and thermal radiation; the height area distribution function of contact interface asperity obeys the fractal model. Based on the GMT of the ultrasonic machining system, the contact thermal resistance of the material matrix and the contact thermal resistance of the interface gap are negligible compared with the heat flux shrinkage thermal resistance. Therefore, we only considered the heat flux shrinkage thermal resistance.

Microscopically, the rough interface contact between the MCS and the horn manifests as mutual contact between the asperities, which is further manifested as peak-to-peak contact between the asperities. We did not consider other contact types, such as side contact or mixed contact. When the heat flow propagates at two rough contact interfaces, it produces thermal shrinkage, namely shrinkage thermal resistance. The degree of heat flow contraction is inversely proportional to the number of contacted asperities. According to the classical truncated cone contact model, the thermal shrinkage resistance at the contact point of a single asperity can be expressed as [28]:

(19)$r_{\mathrm{c}}=\frac{\psi (c/b)}{2k_{\mathrm{s}}c}$

where b and c are the radius of the heat flow channel and the contact point, respectively, which can be solved by the equation $c/b=A_{\mathrm{r}}^{*0.5}$; $A_{\mathrm{r}}^{*}$ is the dimensionless real contact area; $\psi (c/b)$ is the interface shrinkage, which can be approximately expressed as $\psi (c/b)={(1-c/b)}^{1.5}$; $k_{\mathrm{s}}$ is the equivalent thermal conductivity, which can be calculated by $2/k_{\mathrm{s}}=1/k_1+1/k_2$, where $k_1$ and $k_2$ are the thermal conductivity of the two contact materials, respectively. Substituting contact radius R of asperity into Equation (19), contact thermal conductance of asperity becomes:

(20)$h_{\mathrm{c}}=1/r_{\mathrm{c}}=\frac{2k_{\mathrm{s}}R}{\left(1-{\left(A_{\mathrm{r}}^{*}0.5\right)}^{1.5}\right)}$

Then, the contact thermal conductivity of asperity in the plastic stage, elastic–plastic stage 1 and 2, and elastic deformation stage can be further expressed as:

(21)$h_{\mathrm{ce}}=\frac{2k_{\mathrm{s}}{a}_{\mathrm{e}}^{*0.5}}{\sqrt{2\pi }{\left(\left(1-A_{\mathrm{r}}^{*0.5}\right)\right)}^{1.5}}$

(22)$h_{\mathrm{cep}-1}=\frac{2k_{\mathrm{s}}{a}_{\mathrm{ep}-1}^{*0.5}}{\sqrt{2\pi }{\left(\left(1-A_{\mathrm{r}}^{*0.5}\right)\right)}^{1.5}}$

(23)$h_{\mathrm{cep}-2}=\frac{2k_{\mathrm{s}}{a}_{\mathrm{ep}-2}^{*0.5}}{\sqrt{2\pi }{\left(\left(1-A_{\mathrm{r}}^{*0.5}\right)\right)}^{1.5}}$

(24)$h_{\mathrm{cp}}=\frac{2k_{\mathrm{s}}{a}_{\mathrm{p}}^{*0.5}}{\sqrt{2\pi }{\left(\left(1-A_{\mathrm{r}}^{*0.5}\right)\right)}^{1.5}}$

where $a_{\mathrm{e}}^{}$, $a_{\mathrm{e}-1}^{}$, $a_{\mathrm{e}-2}^{}$, and $a_p^{}$ represent the contact area of the asperity in the elastic, elastic–plastic 1 and 2, and plastic stages, respectively, which can be obtained based on the fractal model of contact mechanics in each stage.

According to the statistical principle, the contact thermal conductivity of the real area of the joint surface can be obtained based on the idea of a single microconvex body to the entire interface:

(25)$\begin{array}{l}{H}_{\mathrm{c}}={\displaystyle {\int }_{a_c}^{a_l}{h}_{\mathrm{ce}}}n(a)da+{\displaystyle {\int }_{6^{-{(D-1)}^{-1}}{a}_c}^{a_c}{h}_{\mathrm{cep}-1}}n(a)da\\ \hspace{1em} +{\int }_{{110}^{-{(\mathrm{D}-1)}^{-1}}{a}_c}^{6^{-{(\mathrm{D}-1)}^{-1}}{a}_c}{h}_{\mathrm{cep}-1}n(a)da+{\displaystyle {\int }_0^{{110}^{-{(\mathrm{D}-1)}^{-1}}{a}_c}{h}_{\mathrm{cp}}}n(a)da\end{array}$

2.3.2. Composite Stiffness

A single asperity’s steady thermoelastic deformation under thermal effect is analyzed based on thermoelastic strain. Although the thermal strain of a contact interface can also be accurately expressed from the perspective of the one-dimensional height of the asperity, it requires more assumptions and more complex mathematical calculations. At the same time, in order to ensure the consistency and universality of the calculation of the normal contact force in Section 2.2, the calculation of the thermal–elastic force is simplified, that is, the contact thermal–elastic deformation force between the asperities can be expressed as:

(26)$p_{\mathrm{T}}=\alpha E\Delta \mathrm{T}{a}_{\mathrm{e}}$

where α is the linear expansion coefficient of the material, and E is the equivalent elastic modulus; $\Delta T$ is the temperature difference with the environment; and $a_c$ is the actual deformation area of elastic deformation stage. Combined with (10), the normal force of a single asperity containing thermal stress in the elastic deformation stage can be obtained as follows:

(27)$\begin{array}{l}{p}_{T\mathrm{e}}=\frac{4}{3}E{\left(\frac{a^{0.5\mathrm{D}}}{2^{4-\mathrm{D}}{\pi }^{0.5\mathrm{D}}{\mathrm{G}}^{\mathrm{D}-1}{(\ln \gamma )}^{0.5}}\right)}^{0.5}\\ \hspace{1em} \times {\left(2^{3-\mathrm{D}}{\mathrm{G}}^{\mathrm{D}-1}{(\ln \gamma )}^{0.5}{\pi }^{0.5\mathrm{D}-1}{a}^{0.5\mathrm{D}-1}\right)}^{1.5}+\alpha E\Delta \mathrm{T}{a}_e\end{array}$

In Equation (27), the $d{p}_{\mathrm{Te}}/d\delta$ value in the elastic deformation regime can be determined by differentiating Equation (28), which produces:

(28)$\frac{d{p}_{\mathrm{Te}}}{d\delta }={\mathrm{k}}_{\mathrm{Te}}=\frac{4E(3-\mathrm{D})}{3\sqrt{2\pi }(2-\mathrm{D})}{a}^{0.5}+\frac{{\pi }^{1-0.5\mathrm{D}}\alpha E\Delta \mathrm{T}{a}^{\mathrm{D}/2}}{2^{3-\mathrm{D}}(1-0.5\mathrm{D}){(\mathrm{In}\gamma )}^{1/2}{\mathrm{G}}^{\mathrm{D}-1}}$

where the former is the static elastic stiffness of the mechanical load, and the latter is the elastic deformation stiffness of the thermal response, namely, the thermal stiffness. According to statistics from local to global, the thermal–elastic–plastic normal stiffness of the joint surface is composite stiffness. Regarding coupled thermal stiffness, we take into account only the elastic stage of deformation without considering the thermal effects of other stages of deformation. We also ignore the impact of coupled thermal stiffness on plastic deformation. This is very important. Although it is rough, it dramatically reduces the complexity and improves computational efficiency. The thermal stiffness of the joint surface can be expressed as:

(29)${\mathrm{K}}_{\mathrm{T}}={\displaystyle {\int }_{a_c}^{a_l}\frac{{\pi }^{1-0.5\mathrm{D}}\alpha E\Delta \mathrm{T}{a}^{\mathrm{D}/2}}{2^{3-\mathrm{D}}(1-0.5\mathrm{D}){(\mathrm{In}\gamma )}^{1/2}{\mathrm{G}}^{\mathrm{D}-1}}}n(a)da$

After dimensionless processing, (29) is expressed as:

(30)$\begin{array}{l}{\mathrm{K}}_{\mathrm{T}}^{*}={\displaystyle {\int }_{a_c^{*}}^{a_l^{*}}\frac{{\pi }^{1-0.5\mathrm{D}}\alpha E\Delta \mathrm{T}{a}^{\mathrm{D}/2}}{2^{3-\mathrm{D}}(1-0.5\mathrm{D}){(\mathrm{In}\gamma )}^{1/2}{\mathrm{G}}^{*\mathrm{D}-1}}}n\left(a^{*}\right)d{a}^{*}\\ \hspace{1em} =\frac{{\pi }^{1-0.5\mathrm{D}}\alpha \Delta \mathrm{T}}{2^{3-\mathrm{D}}{(\mathrm{In}\gamma )}^{1/2}{\mathrm{G}}^{*}{}^{\mathrm{D}-1}}{g}_4\mathrm{In}\left(\frac{A_{\mathrm{r}}^{*}}{a_c^{*}}\frac{2-\mathrm{D}}{\mathrm{D}}\right)\end{array}$

By applying Equations (3), (10) and (18), the deformation force and contact stiffness function from a single asperity to the whole interface based on statistical principles can be deduced. After dimensionless processing, the total contact force containing thermoelastic stress can be calculated as follows when the fractal dimension D is not equal to 1.5:

(31)$\begin{array}{l}{\mathrm{P}}^{*}={\displaystyle {\int }_{a_c^{*}}^{a_l^{*}}{p}_{\mathrm{e}}^{*}}n\left(a^{*}\right)d{a}^{*}+f_{\mathrm{T}}^{*}+{\displaystyle {\int }_{6^{-(D-1)-1}{a}_c^{*}}^{a_c^{*}}{p}_{1\mathrm{ep}}^{*}}n\left(a^{*}\right)d{a}^{*}\\ \hspace{1em} +{\displaystyle {\int }_{{110}^{-(\mathrm{D}-1)-1}{a}_c^{*}}^{6^{-(\mathrm{D}-1)-1}{a}_c^{*}}{p}_{2\mathrm{ep}}^{*}}n\left(a^{*}\right)d{a}^{*}+{\displaystyle {\int }_0^{{110}^{-(\mathrm{D}-1)-1}{a}_c^{*}}{p}_{\mathrm{p}}^{*}}n\left(a^{*}\right)d{a}^{*}\\ \hspace{1em} =\frac{2^{4.5-\mathrm{D}}{\mathrm{G}}^{*\mathrm{D}-1}{(\mathrm{In}\gamma )}^{0.5}\mathrm{D}{g}_1{A}_{\mathrm{r}}^{*0.5\mathrm{D}}}{3{\pi }^{1.5-0.5\mathrm{D}}(3-2\mathrm{D})}\cdot \left(g_2{A}_{\mathrm{r}}^{*1.5-\mathrm{D}}-a_c^{*1.5-\mathrm{D}}\right)\\ \hspace{1em} +\frac{1.03\times 2^{2.825-0.85\mathrm{D}}{(2.8K\phi )}^{0.15}{\pi }^{0.425\mathrm{D}-1.275}{\mathrm{G}}^{*0.85\mathrm{D}-0.85}{(\mathrm{In}\gamma )}^{0.425}}{3(1.425-0.925\mathrm{D})}\\ \hspace{1em} \times \mathrm{D}{g}_1{A}_{\mathrm{r}}^{*}{}^{0.5\mathrm{D}}{a}_c^{*1.425-0.925\mathrm{D}}\left[1-6^{-(1.425-0.925\mathrm{D})/\mathrm{D}-1}\right]\\ \hspace{1em} +1.4\times 2^{1.367-0.526\mathrm{D}}{(2.8K\phi )}^{0.474}{\pi }^{0.263\mathrm{D}-0.789}{\mathrm{G}}^{*0.526(\mathrm{D}-1)}{(\mathrm{In}\gamma )}^{0.263}\\ \hspace{1em} \times \left[6^{-(1.263-0.763\mathrm{D})/(\mathrm{D}-1)}-{110}^{-(1.263-0.763\mathrm{D})/(\mathrm{D}-1)}\right]\\ \hspace{1em} \times \frac{ \mathrm{D}{g}_1{A}_{\mathrm{r}}^{*0.5\mathrm{D}}{a}_c^{*1.263-0.763\mathrm{D}} }{ 3(1.263-0.763\mathrm{D})}\\ \hspace{1em} +\frac{ 2.8\phi \mathrm{D}{g}_1{A}_{\mathrm{r}}^{{*}^{{}_{0.5\mathrm{D}}}}{a}_c^{*1-0.5\mathrm{D}}}{ 2-\mathrm{D}}\cdot {110}^{-(2-\mathrm{D})/(2\mathrm{D}-2)}\\ \hspace{1em} +\alpha \Delta \mathrm{T}\left[A_{\mathrm{r}}^{*}-A_{\mathrm{r}}^{*0.5\mathrm{D}}{\left(\frac{2-\mathrm{D}}{\mathrm{D}}\right)}^{0.5\mathrm{D}-1}{a}_c^{*(1-0.5\mathrm{D})}\right]\end{array}$

At a given mean surface separation distance, when the fractal dimension D is equal to 1.5, the total contact force containing thermoelastic stress can be calculated as follows:

(32)$\begin{array}{l}{\mathrm{P}}^{*}={\displaystyle {\int }_{a_c^{*}}^{a_l^{*}}{p}_{\mathrm{e}}^{*}}n\left(a^{*}\right)d{a}^{*}+f_{\mathrm{T}}^{*}+{\displaystyle {\int }_{6^{-(\mathrm{D}-1)-1}{a}_c^{*}}^{a_c^{*}}{p}_{1\mathrm{ep}}^{*}}n\left(a^{*}\right)d{a}^{*}\\ \hspace{1em} +{\displaystyle {\int }_{{110}^{-(\mathrm{D}-1)-1}{a}_c^{*}}^{6^{-(\mathrm{D}-1)-1}{a}_c^{*}}{p}_{2\mathrm{ep}}^{*}}n\left(a^{*}\right)d{a}^{*}+{\displaystyle {\int }_0^{{110}^{-(\mathrm{D}-1)-1}{a}_c^{*}}{p}_{\mathrm{p}}^{*}}n\left(a^{*}\right)d{a}^{*}\\ \hspace{1em} =\frac{2^{4.5-\mathrm{D}}{\mathrm{G}}^{*\mathrm{D}-1}{(\mathrm{In}\gamma )}^{0.5}\mathrm{D}{g}_1{A}_{\mathrm{r}}^{*0.5\mathrm{D}}}{3{\pi }^{1.5-0.5\mathrm{D}}(3-2\mathrm{D})}\cdot \left(g_2{A}_{\mathrm{r}}^{**1.5-\mathrm{D}}-a_c^{*1.5-\mathrm{D}}\right)\\ \hspace{1em} +\frac{1.03\times 2^{2.825-0.85\mathrm{D}}{(2.8K\phi )}^{0.15}{\pi }^{0.425\mathrm{D}-1.275}{\mathrm{G}}^{*0.85\mathrm{D}-0.85}{(\mathrm{In}\gamma )}^{0.425}}{3(1.425-0.925\mathrm{D})}\\ \hspace{1em} \times \mathrm{D}{g}_1{A}_{\mathrm{r}}^{*0.5\mathrm{D}}{a}_c^{*1.425-0.925\mathrm{D}}\left[1-6^{-(1.425-0.925\mathrm{D})/\mathrm{D}-1}\right]\\ \hspace{1em} +1.4\times 2^{1.367-0.526\mathrm{D}}{(2.8K\phi )}^{0.474}{\pi }^{0.263\mathrm{D}-0.789}{\mathrm{G}}^{*0.526(\mathrm{D}-1)}{(\mathrm{In}\gamma )}^{0.263}\\ \hspace{1em} \times \left[6^{-(1.263-0.763\mathrm{D})/(\mathrm{D}-1)}-{110}^{-(1.263-0.763\mathrm{D})/(\mathrm{D}-1)}\right]\\ \hspace{1em} \times \frac{ \mathrm{D}{g}_1{A}_{\mathrm{r}}^{*0.5\mathrm{D}}{a}_c^{*1.263-0.763\mathrm{D}} }{ 3(1.263-0.763\mathrm{D})}\\ \hspace{1em} +\frac{ 2.8\phi \mathrm{D}{g}_1{A}_{\mathrm{r}}^{*0.5\mathrm{D}}{a}_c^{*1-0.5\mathrm{D}} }{ 2-\mathrm{D} }\cdot {110}^{-(2-\mathrm{D})/(2\mathrm{D}-2)}\\ \hspace{1em} +\alpha \Delta \mathrm{T}\left[A_{\mathrm{r}}^{*}-A_{\mathrm{r}}^{*0.5\mathrm{D}}{\left(\frac{ 2-\mathrm{D} }{ \mathrm{D}}\right)}^{0.5\mathrm{D}-1}{a}_c^{*(1-0.5\mathrm{D})}\right]\end{array}$

Then, the total contact stiffness of the whole joint surface containing thermal stiffness can be unified as:

(33)$\begin{array}{l}{\mathrm{K}}^{*}={\displaystyle {\int }_{a_c^{*}}^{a_l^{*}}{\mathrm{k}}_{\mathrm{Te}}^{*}}n\left(a^{*}\right)d{a}^{*}+{\displaystyle {\int }_{6^{-{(\mathrm{D}-1)}^{-1}}{a}_c^{*}}^{a_c^{*}}{\mathrm{k}}_{1\mathrm{ep}}^{*}}n\left(a^{*}\right)d{a}^{*}+{\displaystyle {\int }_{{110}^{-{(\mathrm{D}-1)}^{-1}}{a}_c^{*}}^{6^{-{(\mathrm{D}-1)}^{-1}}{a}_c^{*}}{\mathrm{k}}_{2\mathrm{ep}}^{*}}n\left(a^{*}\right)d{a}^{*}\\ \hspace{1em} =\frac{4(3-\mathrm{D})\mathrm{D}{g}_1{A}_{\mathrm{r}}^{*\mathrm{D}/2}}{3\sqrt{2\pi }(2-\mathrm{D})(1-\mathrm{D})}\cdot \left[g_3{A}_{\mathrm{r}}^{*(1-\mathrm{D})/2}-a_c^{*(1-\mathrm{D})/2}\right]\\ \hspace{1em} +\frac{{\pi }^{1-0.5\mathrm{D}}\alpha \Delta \mathrm{T}}{2^{3-\mathrm{D}}{(\mathrm{In}\gamma )}^{1/2}{\mathrm{G}}^{*-1}}{g}_4{\left(A_{\mathrm{r}}^{*}\right)}^{\mathrm{D}/2}\mathrm{In}\left(\frac{A_{\mathrm{r}}^{*}}{a_c^{*}}\frac{2-\mathrm{D}}{\mathrm{D}}\right)\\ \hspace{1em} +\frac{1.03\times 2^{0.825+0.15\mathrm{D}}{(2.8K\phi )}^{0.15}\mathrm{D}{g}_1{A}_{\mathrm{r}}^{*0.5\mathrm{D}}{a}_c^{*(1-\mathrm{D})0.425}\left(1-6^{0.425}\right)}{3{\pi }^{0.075\mathrm{D}+0.275}{\mathrm{G}}^{*0.15(\mathrm{D}-1)}{(\mathrm{In}\gamma )}^{0.075}(1.425-0.425\mathrm{D})(2-\mathrm{D})\times 0.424(1-\mathrm{D})}\\ \hspace{1em} +1.4\times 2^{0.474\mathrm{D}-1.633}{(2.8K\phi )}^{0.474}{E}^{0.526}\mathrm{D}{g}_1{A}_{\mathrm{r}}^{*0.5\mathrm{D}}{a}_c^{*(1-\mathrm{D})0.263}\\ \hspace{1em} \times \frac{6^{0.263}-{110}^{0.263}}{3{\pi }^{0.237\mathrm{D}-0.211}{\mathrm{G}}^{*0.474(\mathrm{D}-1)}{(\mathrm{In}\gamma )}^{0.237}(1.263-0.263\mathrm{D})(2-\mathrm{D})\times 0.263(1-\mathrm{D})}\end{array}$

where $A_{\mathrm{a}}$ is the nominal contact area, $A_{\mathrm{r}}$ is the real contact area, $a_c$ is the critical contact area, and the dimensionless processing can be obtained by:

(34)$$\begin{gathered} g_1={\left(\frac{2-\mathrm{D}}{\mathrm{D}}\right)}^{\mathrm{D}/2};g_2={\left(\frac{2-\mathrm{D}}{\mathrm{D}}\right)}^{(3-2\mathrm{D})/2};g_3={\left(\frac{2-\mathrm{D}}{\mathrm{D}}\right)}^{(1-\mathrm{D})/2}; \\ {g}_4={\left(\frac{\mathrm{D}}{2-\mathrm{D}}\right)}^{(2-\mathrm{D})/2};{\mathrm{P}}^{*}=\frac{\mathrm{P}}{E{A}_a};{\mathrm{G}}^{*}=\frac{\mathrm{G}}{\sqrt{A_{\mathrm{a}}}};A_{\mathrm{r}}^{*}=\frac{A_{\mathrm{r}}}{A_{\mathrm{a}}};a_{\mathrm{c}}^{*}=\frac{a_{\mathrm{c}}}{A_{\mathrm{a}}};\phi =\frac{{\sigma }_y}{E} \end{gathered}$$

3. Results and Discussion

Numerical Simulation

Regarding the effects of temperature on the elastic modulus E and linear expansion coefficient, it is necessary to take corresponding values according to the variable physical parameters of related materials. The results, K^*, as shown in Figure 4 are expressed as a function of the dimensionless temperature change ($\Delta {\mathrm{T}}^{*}$). D is taken as 1.1–1.9, G is taken as 10⁻¹⁰, 10⁻¹¹, and 10⁻¹². As shown in the Figure 4, when the modified fractal model considering temperature deformation is considered, the dimensionless contact stiffness K^* increases with the increase in $\Delta {\mathrm{T}}^{*}$ and presents a robust linear relationship. For the same change temperature, K^* decreases with the increase in G within a specific range, and this change becomes particularly significant at G = 10⁻¹², while it becomes an inclined line with a small slope at G = 10⁻¹⁰ $\mathrm{D}>1.5$. The increase in G indicates the increase in roughness of a rough surface, and the elastic contact percentage of the thermal effect between rough contact surfaces decreases, indicating that the minor roughness of the joint surface has a significant influence on the thermal effect of the interface.

Figure 5 shows the variation in $\mathrm{K}*$ with $\mathrm{P}*$ under different D and G when the dimensionless temperature $\Delta {\mathrm{T}}^{*}$ is 300. The dimensional normal contact stiffness $\mathrm{K}*$ of the joint surface, considering the thermal effect, increases with the increase in the dimensional normal total contact load $\mathrm{P}*$; that is, the normal deformation of the joint surface, the normal contact load, and the normal contact stiffness are proportional. In the range of fractal dimension D = 1.1~1.5, in the small range of the contact force, $\mathrm{K}*$ and $\mathrm{P}*$ show a significant nonlinear relationship, as shown in Figure 5a–d, that is, strong interface nonlinearity. The contact force and stiffness in this range have an essential influence on the stability of the ultrasonic vibration system. In the design, the influence of the nonlinearity of this process on the stability of ultrasonic vibration should be avoided. When $\mathrm{D}⩾1.7$, the linear relationship of interface characteristic parameters is increasingly prominent, as shown in Figure 5e,f. This range is beneficial to the stability of the ultrasonic vibration system due to the linear interface characteristics. In short, under the thermal effect of a specific temperature, the actual elastic contact area and normal contact stiffness of the rough contact interface increase with the increase in the contact load.

Figure 6 is the influence curve of fractal dimension D on normal contact stiffness K* when $\Delta {\mathrm{T}}^{*}$ = 300, $\mathrm{G}={10}^{-12}$. As shown in the figure, within a specific range of contact force ($1.1⩽\mathrm{D}⩽1.5$), the curve shows a nonlinear solid relationship. When $\mathrm{D}\ge 1.6$, the linear interface relationship is apparent, which is consistent with the corresponding relationship in Figure 5. The nonlinear interface relationship is again verified.

4. Conclusions

In this study, the composite contact stiffness model of an ultrasonic vibration system under a preload due to the thermal effect was established. The following conclusions were obtained:

• The change in interface stiffness under the thermal effect increases with temperature, and the linear relationship between the temperature and interface stiffness is apparent. Therefore, a reasonable cooling method must be adopted in the GMT system in the operation process so that the system can work in a reasonable temperature range.

• When $1.1⩽\mathrm{D}⩽1.5$, the interface nonlinearity is pronounced. For the same change temperature, K* decreases with the increase in G within a specific range, and this change becomes particularly significant at $\mathrm{G}={10}^{-12}$.

• The proposed interface thermal stiffness contact model can provide theoretical guidance for the study and optimal design of nonlinear dynamic characteristics of ultrasonic vibration systems (e.g., rotating giant magnetostrictive transducers) or other thermal products from the point of view of the contact interface.

Author Contributions

Conceptualization, Y.C.; methodology, P.L.; software, J.S; validation, J.S.; formal analysis, Y.C.; investigation, J.S., M.S. and L.S.; resources, P.L.; data curation, Y.C.; writing—original draft preparation, Y.C.; writing—review and editing, P.L.; visualization, J.S.; supervision, L.S.; project administration, P.L.; funding acquisition, P.L. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data are contained within the article as figures and tables.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

${\Phi }_{1,n}$	random phase
n	spatial frequency index
M	number of overlapping uplifts on the structural surface
G	fractal roughness coefficient determined by frequency
L	sample length
$n(a)$	area distribution density of the contact point
N	number of asperities larger than contact area a
$A_{\mathrm{r}}$	real contact area of joint surface
$\delta$	peak-to-valley amplitude of the z(x) function
R	curvature radius of the asperity
${\delta }_c$	elastic critical deformation
K	the hardness coefficient
$\nu$	Poisson’s ratio
$a_c$	critical contact area of asperity
$P_e$	elastic contact load
$p_{1ep}$	contact mechanics of elastic–plastic deformation zone 1
$p_{2ep}$	contact mechanics of elastic–plastic deformation zone 2
$p_p$	contact mechanics of complete plastic deformation
$r_{\mathrm{c}}$	thermal shrinkage resistance at contact point of a single asperity
$h_{\mathrm{c}}$	contact thermal conductance of asperity
$H_{\mathrm{c}}$	contact thermal conductivity of real area of joint surface
$p_{\mathrm{T}}$	normal load caused by thermal stress
$p_{T\mathrm{e}}$	normal force of a single asperity containing thermal stress in elastic deformation stage
${\mathrm{P}}^{*}$	total contact force containing thermoelastic stress
${\mathrm{K}}^{*}$	total contact stiffness of whole joint surface containing thermal stiffness
$A_{\mathrm{a}}$	nominal contact area

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Table

Figure 1. Micro-surface morphology characteristics of a magnetic conductive sheet.

Figure 2. Schematic diagrams of the geometry of contact; (a) the rigid plane and single asperity (b) the rigid plane and actual rough surface.

Figure 3. Ultrasound system–interface stiffness model under thermal effect.

Figure 4. Influence of length scale parameter G* on [Forumla omitted. See PDF.]-K* curve; (a) D = 1.1 (b) D = 1.3 (c) D = 1.5 (d) D = 1.7 (e) D = 1.9.

Figure 5. Influence of length scale parameter G* on [Forumla omitted. See PDF.]-K* curve with [Forumla omitted. See PDF.] = 300; (a) D = 1.1 (b) D = 1.2 (c) D = 1.3 (d) D = 1.5 (e) D = 1.7 (f) D = 1.9.

Figure 6. Influence of fractal dimensions D on [Forumla omitted. See PDF.]-K* curve with [Forumla omitted. See PDF.] = 300; (a) D = 1.1, 1.2 (b) D = 1.2, 1.3 (c) D = 1.3, 1.4 (d) D = 1.4, 1.5 (e) D = 1.5, 1.6 (f) 1.6, 1.7 (g) 1.7, 1.8, 1.9.

Table 1

Fractal parameter values under the same measured surface.

Experimental Times	1	2	3	4	5	6	7	8
D	1.4007	1.4061	1.4021	1.4613	1.4026	1.4525	1.4533	1.4060
G (10−10 mm)	6.5340	8.6558	9.2133	24.005	25.682	8.9967	41.003	31.017

References

1. Ma, K.; Wang, J.; Zhang, J.; Feng, P.; Yu, D.; Ahmad, S. A highly temperature-stable and complete-resonance rotary giant magnetostrictive ultrasonic system. Int. J. Mech. Sci.; 2022; 214, 106927. [DOI: https://dx.doi.org/10.1016/j.ijmecsci.2021.106927]

2. Zhou, H.; Zhang, J.; Feng, P.; Yu, D.; Wu, Z. An amplitude prediction model for a giant magnetostrictive ultrasonic transducer. Ultrasonics; 2020; 108, 106017. [DOI: https://dx.doi.org/10.1016/j.ultras.2019.106017] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32690341]

3. Zhan, Y.-S.; Lin, C.-H. A Constitutive Model of Coupled Magneto-thermo-mechanical Hysteresis Behavior for Giant Magnetostrictive Materials. Mech. Mater.; 2020; 148, 103477. [DOI: https://dx.doi.org/10.1016/j.mechmat.2020.103477]

4. Wang, T.-Z.; Zhou, Y.-H. Nonlinear dynamic model with multi-fields coupling effects for giant magnetostrictive actuators. Int. J. Solids Struct.; 2013; 50, pp. 2970-2979. [DOI: https://dx.doi.org/10.1016/j.ijsolstr.2013.05.012]

5. Xiao, Y.; Gou, X.F.; Zhang, D.G. A one-dimension nonlinear hysteretic constitutive model with elasto-thermo-magnetic coupling for giant magnetostrictive materials. J. Magn. Magn. Mater.; 2017; 441, pp. 642-649. [DOI: https://dx.doi.org/10.1016/j.jmmm.2017.06.050]

6. Zhan, Y.-S.; Lin, C.-H. Micromechanics-based constitutive modeling of magnetostrictive 1–3 and 0–3 composites. Compos. Struct.; 2021; 260, [DOI: https://dx.doi.org/10.1016/j.compstruct.2020.113264]

7. Li, P.; Liu, Q.; Zhou, X.; Xu, G.; Li, W.; Wang, Q.; Yang, M. Effect of Terfenol-D rod structure on vibration performance of giant magnetostrictive ultrasonic transducer. J. Vib. Control; 2020; 27, pp. 573-581. [DOI: https://dx.doi.org/10.1177/1077546320932029]

8. Chen, Y.; Zhang, X.; Wen, S.; Lan, G. Fractal Model for Normal Contact Damping of Joint Surface Considering Elastoplastic Phase. J. Mech. Eng.; 2019; 55, pp. 58-68.

9. Kogut, L.; Etsion, I. Elastic-Plastic Contact Analysis of a Sphere and a Rigid Flat. J. Appl. Mech.; 2002; 69, pp. 657-662. [DOI: https://dx.doi.org/10.1115/1.1490373]

10. Kogut, L.; Etsion, I. A Finite Element Based Elastic-Plastic Model for the Contact of Rough Surfaces. Tribol. Trans.; 2003; 46, pp. 383-390. [DOI: https://dx.doi.org/10.1080/10402000308982641]

11. Brake, M.R. An analytical elastic-perfectly plastic contact model. Int. J. Solids Struct.; 2012; 49, pp. 3129-3141. [DOI: https://dx.doi.org/10.1016/j.ijsolstr.2012.06.013]

12. Xu, C.; Wang, D. An Improved Analytical Model for Normal Elastic-Plastic Contact of Rough Surfaces. J. Xi’an Jiaotong Univ.; 2014; 48, pp. 115-121.

13. Jiang, S.; Zheng, Y.; Zhu, H. A Contact Stiffness Model of Machined Plane Joint Based on Fractal Theory. J. Tribol.-Trans. Asme; 2010; 132, 011401. [DOI: https://dx.doi.org/10.1115/1.4000305]

14. Raffa, M.L.; Lebon, F.; Vairo, G. Normal and tangential stiffnesses of rough surfaces in contact via an imperfect interface model. Int. J. Solids Struct.; 2016; 87, pp. 245-253. [DOI: https://dx.doi.org/10.1016/j.ijsolstr.2016.01.025]

15. Liu, W.; Yang, J.; Xi, N.; Shen, J.; Li, L. A study of normal dynamic parameter models of joint interfaces based on fractal theory. J. Adv. Mech. Des. Syst. Manuf.; 2015; 9, JAMDSM0070. [DOI: https://dx.doi.org/10.1299/jamdsm.2015jamdsm0070]

16. Liou, J.L.; Tsai, C.M.; Lin, J.-F. A microcontact model developed for sphere- and cylinder-based fractal bodies in contact with a rigid flat surface. Wear; 2010; 268, pp. 431-442. [DOI: https://dx.doi.org/10.1016/j.wear.2009.08.033]

17. Liou, J.L.; Lin, J.F. A modified fractal microcontact model developed for asperity heights with variable morphology parameters. Wear; 2010; 268, pp. 133-144. [DOI: https://dx.doi.org/10.1016/j.wear.2009.07.003]

18. Wang, Z.; Liu, X. Model for Elastic-Plastic Contact between Rough Surfaces. J. Tribol.; 2018; 140, 051402. [DOI: https://dx.doi.org/10.1115/1.4040385]

19. He, L.G.; Zuo, Z.X.; Xiang, J.H. Normal Contact Stiffness Fractal Model Considering Asperity Elastic-Plastic Transitional Deformation Mechanism of Joints. J. Shanghai Jiaotong Univ.; 2015; 49, pp. 116-121.

20. Li, L.; Liang, X.; Xing, Y.Z.; Yan, D.; Wang, G. Measurement of Real Contact Area for Rough Metal Surfaces and the Distinction of Contribution From Elasticity and Plasticity. J. Tribol.; 2021; 143, 071501. [DOI: https://dx.doi.org/10.1115/1.4048728]

21. Kogut, L.; Komvopoulos, K. Electrical contact resistance theory for conductive rough surfaces separated by a thin insulating film. J. Appl. Phys.; 2004; 95, pp. 576-585. [DOI: https://dx.doi.org/10.1063/1.1629392]

22. Wang, S.; Komvopoulos, K. Closure to “Discussion of ‘A Fractal Theory of the Interfacial Temperature Distribution in the Slow Sliding Regime: Part I—Elastic Contact and Heat Transfer Analysis’” (1994, ASME J. Tribol., 116, p. 822). J. Tribol.; 1994; 116, pp. 822-823. [DOI: https://dx.doi.org/10.1115/1.2927340]

23. Shi, X.; Wu, A.; Jin, C.; Qu, S. Thermomechanical modeling and transient analysis of sliding contacts between an elastic–plastic asperity and a rigid isothermal flat. Tribol. Int.; 2015; 81, pp. 53-60. [DOI: https://dx.doi.org/10.1016/j.triboint.2014.08.004]

24. Song, W.; Li, L.; Ovcharenko, A.; Talke, F.E. Thermo-mechanical contact between a rigid sphere and an elastic–plastic sphere. Tribol. Int.; 2016; 95, pp. 132-138. [DOI: https://dx.doi.org/10.1016/j.triboint.2015.11.019]

25. Horovistiz, A.; Laranjeira, S.; Davim, J.P. 3-D reconstruction by extended depth-of-field in tribological analysis: Fractal approach of sliding surface in Polyamide66 with glass fiber reinforcement. Polym. Test.; 2019; 73, pp. 178-185. [DOI: https://dx.doi.org/10.1016/j.polymertesting.2018.11.017]

26. Yan, W.; Komvopoulos, K. Contact analysis of elastic-plastic fractal surfaces. J. Appl. Phys.; 1998; 84, pp. 3617-3624. [DOI: https://dx.doi.org/10.1063/1.368536]

27. Majumdar, A.; Bhushan, B. Fractal Model of Elastic-Plastic Contact Between Rough Surfaces. J. Tribol.; 1991; 113, pp. 1-11. [DOI: https://dx.doi.org/10.1115/1.2920588]

28. Cooper, M.G.; Mikic, B.B.; Yovanovich, M.M.; Transfer, M. Thermal Contact Conductance. Int. J. Heat Mass Transf.; 1969; 12, pp. 279-300. [DOI: https://dx.doi.org/10.1016/0017-9310(69)90011-8]