Fractal Model of Contact Thermal Stiffness

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.


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.

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]: 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: 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.
the maximum frequency index related to the distance cut-off length; overlapping uplifts on the structural surface; G is the fractal roughne mined by frequency; and D (1 < D < 2) is the fractal dimension of the p When M=1, m=1, the W-M function that becomes a single indepe The three-dimensional contact profiler is used to measure the m sheet (MCS) of a GMT (see [7]) with different roughnesses, as show yellow line is the center of the measurement position and is measure material is electrical pure iron, the elastic modulus is 195 MPa, Poisso roughness of the measurement sample is 0.356 μm, the minimum sa the rough surface measuring instrument is 0.25 μm, and the sampling results of D and G can be obtained using the structural function met mary data is shown in table 1.   Figure 2 illustrates the conceptual approach to the contact betwe a rough surface. According to Hertz' contact theory, the contact betw face asperities is equivalent to the contact between the rigid plane and To obtain the contact deformation force under a given size, the cont   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]:

Fractal Model of Contact Mechanics of Asperity
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: tween spherical asperities and a rigid plane must be first obtained. If the maximum tact area of the asperity is l a , the number of asperities larger than the contact area be expressed as N, and the area distribution density of the contact point can be expre as [27]:  The interference of an asperity under the force of the rigid plate δ is equal t peak-to-valley amplitude of the z(x) function (see Figure 2). δ is given as [26]: According to the relationship between the deformation δ of the asperity an curvature radius R, considering R≫δ the curvature radius of the asperity is [16]: To determine the contact parameters of asperity, the contact state of asperity mu determined. According to the finite element results in the literature [9], the deform The interference of an asperity under the force of the rigid plate δ is equal to the peak-to-valley amplitude of the z(x) function (see Figure 2). δ is given as [26]: According to the relationship between the deformation δ of the asperity and the curvature radius R, considering R δ the curvature radius of the asperity is [16]: 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 (δ < δ c ), elastic-plastic deformation 1 (δ c δ 6δ c ), elastic-plastic deformation 2 (6δ c δ 110δ c ), and complete plastic deformation (110δ c < δ). δ 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: where H is the hardness of the soft material, K is the hardness coefficient, K = 0.454 + 0.41ν, and ν is Poisson's ratio. When δ c = δ, the corresponding contact area is the critical contact area of the asperity a c . Then, substituting Equations (3) and (4) into Equation (7) gives: 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 (δ < δ 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: If the deformation of the asperity is in the range of 6 −(D−1) −1 a c < a < a c (δ c δ 6δ c ), the contact mechanics model of elastic-plastic deformation zone 1 is expressed in fractal form as: Moreover, if the deformation of the asperity is in the range of 110 , the contact mechanics model of elastic-plastic deformation zone 2 is expressed in fractal form as: When 110δ c < δ, the contact mode is complete plastic deformation, and the contact mechanics model is expressed in fractal form as:

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.

Thermal-Mechanical Coupling Contact Interface
We now present a modified asperity contact stiffness model that tact mechanics model of asperity in Section 2.2, only because the ag propagates between rough contact interfaces, as shown in figure 3. contact model in Section 2.2, which assumes that the asperity is sub compressive stress, this model considers contact thermal resistance a thermal stress on contact deformation when heat flow passes throug face. This model especially considers the influence of thermoelastic d practical interface stiffness model of ultrasonic vibration is establish contact interface's static stiffness with the interface's thermal stiffnes sure, namely the interface composite stiffness.

Thermal Contact Resistance of Rough Surface Interface
When the heat flow generated by a magnetostrictive material rough contact interface of the MCS and the horn, the contact thermal ated at the interface under the action of the interface contact pressur contact force between the asperity and the interface is given in Section face contact thermal resistance mainly includes interface material ma resistance, heat flow shrinkage thermal resistance, and interface gap The contact model of heat flow propagation between the contact inter follows: the temperature difference between contact interfaces is the s fer mode of fluid medium in the contact interface gap is heat conducti convection and thermal radiation; the height area distribution function

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 Sections 2.2 and 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]: 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 * 0.5 r ; A * r is the dimensionless real contact area; ψ(c/b) is the interface shrinkage, which can be approximately expressed as ψ(c/b) = (1 − c/b) 1.5 ; k s is the equivalent thermal conductivity, which can be calculated by 2/k 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: 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: where a e , a e−1 , a e−2 , and a p represent the contact area of the asperity in the elastic, elasticplastic 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:

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: where α is the linear expansion coefficient of the material, and E is the equivalent elastic modulus; ∆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: In Equation (27), the dp Te /dδ value in the elastic deformation regime can be determined by differentiating Equation (28), which produces: (28) 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: After dimensionless processing, (29) is expressed as: Machines 2022, 10, 464 9 of 15 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: where A a is the nominal contact area, A r is the real contact area, a c is the critical contact area, and the dimensionless processing can be obtained by:

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 (∆T * ). D is taken as 1.1-1.9, G is taken as 10 −10 , 10 −11 , and 10 −12 . 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 ∆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 −12 , while it becomes an inclined line with a small slope at G = 10 −10 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 K * with P * under different D and G when the dimensionless temperature ∆T * is 300. The dimensional normal contact stiffness K * of the joint surface, considering the thermal effect, increases with the increase in the dimensional normal total contact load 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, K * and 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 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 ∆T * = 300, G = 10 −12 . As shown in the figure, within a specific range of contact force (1.1 D 1.5), the curve shows a nonlinear solid relationship. When D ≥ 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. ≥ , the linear interface relationship is apparent, which is consistent with the corresponding relationship in Figure 5. The nonlinear interface relationship is again verified.

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: 1.
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. 2.
When 1.1 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 G = 10 −12 .

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

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
Φ 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 r real contact area of joint surface δ peak-to-valley amplitude of the z(x) function R curvature radius of the asperity δ c elastic critical deformation K the hardness coefficient ν 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 c thermal shrinkage resistance at contact point of a single asperity h c contact thermal conductance of asperity H c contact thermal conductivity of real area of joint surface p T normal load caused by thermal stress p Te normal force of a single asperity containing thermal stress in elastic deformation stage P * total contact force containing thermoelastic stress K * total contact stiffness of whole joint surface containing thermal stiffness A a nominal contact area