Fractal variation of three-dimensional surface topography during sliding wear under mixed elastohydrodynamic lubrication

The surface morphology variation of the friction pair during the wear process directly affects the local lubrication effect and contact strength, and this variation is difficult to quantify in mixed lubrication. In this study, a numerical simulation and characterization model method of sliding wear surface morphology evolution based on deterministic elastohydrodynamic lubrication (EHL) model and surface three-dimensional morphology fractal method is established. The effectiveness of the lubrication model is verified by the film thickness measurement, and the actual surface wear coefficient is obtained through the wear experiment. On this basis, the variation of surface morphology and fractal parameters in the sliding wear process is investigated, and the effect of contact load and sliding speed on surface lubrication characteristics and fractal parameters is analyzed. The results show: wear will cause the higher wave peaks on the ground surface to be gradually worn away and the fine structure of the surface is reduced, and the fractal dimension of the morphology is therefore reduced. In addition, the contact area ratio and contact load ratio, as well as the coefficient of friction and surface evolution under different loads and sliding speeds are consistent.


Introduction
Generally, excessive sliding wear is one of the most common surface failure mechanisms in machine components [1,2]. In fact, wear can also affect many operating characteristics of mechanical systems, such as vibration and noise levels, load-bearing capacity and friction [3][4][5]. Surface morphology is an important feature of wear parts, and its evolution during wear will directly affect the contact performance between the surfaces of mechanical parts, and have an important impact on the stability and reliability of the mechanical system [6,7]. A comprehensive grasp of the changes in the surface micro-geometric morphology during the wear process will provide a basis for the fault diagnosis, state recognition and wear resistance evaluation of mechanical equipment.
Shyu et al [16] investigated the changes in the surface characterization parameters of engine piston cylinders during wear. Mezghani et al [17] analyzed the changes in the characteristic parameters such as the root-meansquares of plateau(S pq ) and valleys height (S vq ) of the topography height of the cylindrical surface during the wear process. Cabanettes et al [18] investigated the changes in the topography of the camshaft during the running-in period, and compared the difference between the 3D standard roughness parameters (amplitude, spacing, mixing, volume and functional parameters) of the unworn and worn surfaces. In general, the study on the evolution of the wear surface morphology under dry contact provides a basis for the identification of wear conditions and revealing the wear mechanism.
However, only a few mechanical parts in the engineering work under dry wear, and most of the mechanical parts work under lubrication. Although the wear of properly lubricated contact points is not as severe as dry contact, there is still wear at the contact Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
points. Zuo et al [19] conducted a staged wear test of GCr15 steel and H70 brass on a disc tribometer, and discussed the changes of the three-dimensional morphology during the wear process using the multifractal method. Nyman et al [20] used surface topography parameters arithmetic mean peak curvature (S sc ) and root mean square slope (S dq ) to predict the topography changes of friction materials in wet clutches. Reizer et al [21] discussed the arithmetic mean deviation S a , root mean square deviation S q , maximum valley depth S v and maximum peak height S p of the surface morphology of cast iron and 42CrMo4 steel during wear under lubricated conditions. Zhang et al [22] established a worn surface morphology prediction model, and analyzed the variation of the surface morphology parameters S dq and S dr before and after wear. Zhang et al [23] predicted the evolution of the probability density function of the wear profile and surface height during the running-in process under the condition of line contact mixed lubrication. Pei et al [24]analyzed the root mean square (RMS) change of roughness in the process of point contact mixed lubrication wear through experiments and simulations. In all the above applications, statistical parameters such as surface height, slope and standard deviation of curvature are used to characterize surface roughness. However, these parameters strongly depend on the resolution and sampling length of the roughness measurement [25,26]. Due to the lack of intrinsic scale-independent parameters to characterize the surface evolution process, this has caused difficulties in further research on the lubricating contact and wear performance problems associated with machining rough surfaces. A comprehensive understanding of the law of surface wear under lubrication conditions will help improve the wear and lubrication performance of rough surfaces and optimize the mechanical structure and operating parameters, which is essential for gears, bearings, clutches and other mechanical parts that are prone to wear [27][28][29][30][31]. Fractal theory has good applicability in describing natural phenomena with the characteristics of scale laws and can effectively solve the above problems. Recently, many fractal models have been developed to characterize the surface under dry wear or the lubrication characteristics of a constant surface [32][33][34]. However, the evolution of the fractal dimension of the surface morphology under mixed lubrication wear has not been fully discussed.
The main purpose of this paper is to apply fractal theory to characterize the evolution of surface morphology in mixed lubrication sliding wear. For this reason, this article is based on the surface threedimensional morphology fractal method and deterministic mixed elastohydrodynamic lubrication (EHL). The study established a numerical simulation model of sliding wear. The validity of the lubrication model is verified through experiments, and the wear coefficient of the actual surface is obtained through wear experiments. On this basis, the changes of surface fractal parameters in the process of sliding wear and the influence of speed and contact force on surface fractal parameters are studied. This research is expected to provide theoretical support for reducing wear and improving lubrication performance.

Mixed EHL model
In mixed lubrication, the Hertzian contact area is divided into the EHL area where the two surfaces are separated by the lubricating film, and the rough contact area where the rough peaks directly contact. In the EHL regions, the hydrodynamic pressure is governed by the Reynolds equation [35].
1 3 3 In the contact area, as the film thickness becomes zero, the pressure flow of the Reynolds equation disappears. Therefore, the Reynolds equation can be simplified to the following form [36]: The local lubricant film thickness is calculated as follows: where h 0 (t) is the center film thickness of the rigid body. δ 1 and δ 2 is the roughness amplitudes of surface1 and surface2, respectively. υ is the elastic deformation of the surfaces. Boussinesq integral formula for solving elastic deformation terms, which is calculated as follows: Generally, discrete convolution and fast Fourier transform (DC-FFT) algorithms can be used to speed up the elastic deformation calculation process [37].
The density of lubricating oils is adopt the Dowson-Higgison dense-pressure relationship. The load applied in the EHL area is balanced by the pressure integral in the entire area. The load balance equation is shown as follows: , 7 In this study, the temperature rise analytical formula proposed by Zhu et al [38]. was used to solve the temperature rise characteristics of rough surfaces.
where q is the heat generated by either interface friction and lubricant shearing, Since heat generation and surface temperature are interdependent, surface temperatures T 1 and T 2 require an iterative process to solve.
In the present model, it is assumed that material removal occurs at the local spots where the contact pressure exceeds the yielding limit of the material [39], and wear rate follows the Archard's law shown below: where D is wear rate, p is contact pressure, u is sliding speed, H is hardness, and k is wear coefficient, which is measured by the wear experiment. The purpose of this study is to explore the changes in the morphology during the wear process. The timevarying wear coefficient caused by different wear mechanisms is ignored. The Archard model is used to approximate the wear evolution process. Its wear volume is: where V is the total volume; F is the contact load; S is the total slip of wear. In the ball-disk experiment, the total sliding distance is related to time, disk rotation speed and radius.
In this study, the wear calculation process of the rough surface is illustrated in figure 1. First, the film pressure and film thickness with the smooth surface are calculated according to the initial contact parameters such as load, radius of curvature, and velocity which are used as a basis to calculate the film thickness and pressure of the rough surface. When solving the stable solution of the smooth surface film thickness and film pressure based on the contact parameters, we set δ 1 and δ 2 in the film thickness equation to 0. When solving the stable solution of the film thickness and film pressure of the rough surface, we assign the surface height information δ 1 and δ 2 to the film thickness equation. After the iterative process, the stable solution of the film thickness and film pressure of the rough surface is obtained. Then the film thickness and pressure of the rough surface are used as initial input to start the calculation of the wear results. The direct contact area is determined based on the film thickness distribution. Film thickness of zero is considered as direct contact, and pressure (film pressure and rough contact pressure) is obtained at the same time. Then according to the Archard method, the wear height of each rough contact area is calculated. In the wear calculation, the amount of wear is calculated and the surface is updated in each time step, and the new surface is used as the starting point for the next calculation. The result is output after a certain amount of calculation is completed.

Three-dimensional fractal dimension
The microscopic morphology of the engineering surface is often irregular, with disorder, randomness, and multi-scale. In other words, the surface behavior is fractal. Fractal dimension can be used to measure the characteristics of irregular surface topography. The traditional characterization of the roughness surface has always been expressed by the fractal dimension D of the profile line of the roughness surface. The characterization of the three-dimensional roughness surface is generally set as D s =D+1. In fact, the morphology of the rough surface is very complex, manifested in the spatial distribution of anisotropy, variability and local characteristics. The two-dimensional fractal dimension on a specific line of the roughness morphology cannot be used to characterize the three-dimensional roughness. At present, the methods that can calculate the three-dimensional fractal dimension include Hausdorff, Correlation method, similarity method, box counting method, etc Among them, the box counting method has a clear physical meaning, has the characteristics of less calculation, high precision, and wide application range [40,41]. The box counting method is used in this study to characterize and describe the three-dimensional shape.
Based on the definition of fractal dimension, the fractal dimension D of set A can be defined as: where N is the measurement result when r measures A. Different r measures A will have different N. From this, the fractal dimension D can be calculated. The schematic diagram of the box dimension method to calculate the fractal dimension is shown in figure 2. For an image with a size of M×M, it is assumed that the entire image can be covered with a grid of size e×e×e', and e' = e×G/M, where G is the number of gray levels of the entire image. ue and be are the maximum and minimum gray levels of the image in the area e×e, respectively. The number of grids containing at least one gray level between e and be is denoted as N, and the value of N varies with the value of e. Therefore, the relationship between e and N can be used to calculate the fractal dimension of image gray data. Equation (13) is derived by simple mathematics: where D is the dimension of the image.

Film thickness measurement
The elastohydrolubricating oil film measuring instrument LFT 3.0 is used to measure the film thickness at different speed. The schematic diagram of the test device is shown in figure 3. The ring made of optical glass has a diameter of 150 mm and a thickness of 10 mm, and the inner cylindrical surface is covered with a semi-reflective chromium film. A high-quality steel ball with a diameter of 25.4 mm is in contact with the inner cylindrical surface of the ring. Lubricant is supplied via an oil-air lubrication system through a grease nipple near the entrance area.

Wear coefficient measurement
Wear coefficient is of great significance for assessing the degree of wear of mechanical parts and predicting the service life, operational safety and reliability of mechanical equipment. Since the wear coefficient under lubrication conditions is affected by factors such as lubricating oil viscosity, speed, pressure, etc, it is not representative, and the wear coefficient under dry friction conditions is used. In this study, a rotating module was used to conduct a wear test under dry friction conditions to measure the wear coefficient. The wear coefficient measurement experiment was carried out on a commercial friction tester (MFT3000, Rtec Instruments, USA) as shown in figure 4. In the experiment, the three-dimensional rough surface was    In the test, the disc material is GH15 high temperature alloy, and the ball is made of Cr12MoV steel and carburized to enhance the hardness. The chemical composition of the material is shown in table 1. The surface of the steel ball is a polished surface, while the test surface of the disc is a ground surface. In the test, the normal load is 60 N and the disc speed is 600 r min −1 .

Results and discussion
4.1. Numerical simulation and experimental verification Figure 5 shows the comparison between the lubrication model and the experimental center film thickness at different entrainment speeds. In the experiment, the Hertzian contact stress is 0.59 GPa, and the roll-slip rate is 0.4. It can be found from figure 5 that when the entrainment speed is low, the central contact area in the interference image has almost the same color, which means that the film in the central contact area is approximately flat. When the entrainment speed gradually increases, the interference image shows the typical horseshoe-shaped EHL film characteristics. And the innermost color has gradually moved from both sides to the exit area. In order to obtain effective data on film thickness and coefficient of friction, the measurement result of film thickness and coefficient of friction are obtained by averaging the three measurement results. The maximum relative error of the three sets of center film thickness measurement results is 2.5%, and the maximum relative error of the friction coefficient measurement results is 2.1%. The error bar in figure 5 represents the standard deviation of the three measurement results. The three sets of measurement results are very similar, which means that the measurement results are stable and reliable. It can be seen from the figure 5 that in the range of 0.01∼5 m s −1 , the center film thickness increases with the increase of entrainment speed, and the coefficient of friction gradually decreases with the increase of the speed. The maximum relative error between the center film thickness measured by the experiment and the film thickness calculated by the simulation is 8.4%, and the maximum relative error between the coefficient of friction measured by the experiment and the coefficient of friction calculated by the simulation is 6.8%. In general, the calculated results of the lubrication model are in strong agreement with the experimental results. The similarity of  the two curves shows the effectiveness of the model built in this article. In order to obtain the wear coefficient, the wear volume needs to be determined. Because the accurate wear volume is difficult to calculate, this study uses the average cross-sectional area, which is measured by a three-dimensional surface profilometer, and the product of the wear track circumference to obtain the amount of wear. Figure 6(a) is the wear surface morphology at the time of 0 min, 1 min, 2 min and 3 min respectively, and figure 6(b) is the cross-sectional view of the wear depth at different moments. After preliminary calculation and correction using simulation results, the estimated value of the wear coefficient is determined to be 0.018.

The fractal variation of the three-dimensional surface topography in the process of wear
In this section, the lubrication wear model is used to analyze the changes of the grinding surface lubrication and fractal parameters during the wear process. Taking the typical grinding texture as the initial surface morphology, it has clear directionality and uniform height distribution. In this section, the changes in surface contact parameters under the conditions of  Figure 7 shows the evolution of the film thickness in the wear process.
It can be found from figure 7 that the change in the contact area due to wear has a significant impact on the contact stress. In the initial operation stage from no wear to 1000 cycles, the larger rough peaks are removed and the smaller high-frequency rough peaks become more common, resulting in an increase in the frequency of the corresponding stress peaks. At the same time, the rough peaks gradually become uniform, resulting in a uniform overall appearance of contact pressure and film thickness distribution. When the rough surface is worn in 1500-3000 wear cycles, the deepening of the wear leads to an increase in the contact area and a redistribution of the oil film isolation zone and the contact area, which further increases the overall contact stress. Due to the close relationship between the subsurface stress and the surface contact stress, wear will further affect the actual wear process due to the change of the contact shape and the corresponding redistribution of the entire stress field [39]. Figure 8 shows the evolution of surface accumulated wear and fractal dimension in the wear process.
It can be seen from figures 8(a) and (b) that there is a short running-in period. During the first 1500 cycles, the amount of wear increased steadily, and after 1500 cycles, the amount of wear increased rapidly. One possible reason for this is the change of lubrication state, that is, after 1500 cycles of wear, the surface is severely damaged and the lubrication state shifts from mixed lubrication to boundary lubrication. Another reason is that the deepening of the degree of wear changes the geometric shape of the contact pair, and thus the change in the amount of wear.
The fractal dimension can reflect the complexity of the three-dimensional morphology of the worn surface, the degree of irregularity, and the degree of selfsimilarity between the whole and the part [42,43]. The larger the fractal dimension, the more complex the three-dimensional surface morphology, the more fine structures, and the stronger the overall and local selfsimilarity. It can be seen from figure 8(b) that the fractal dimension presents a trend opposite to the amount of wear with the number of cycles, that is, it decreases gradually and then rapidly decreases. In the initial wear stage, the surface morphology is ground texture, and the fine structure is more, the fractal dimension is larger. At this time, the higher wave crests on the surface are gradually worn away and become a wider bearing plane, and as the wear time increases, the area of the bearing plane gradually increases, which leads to a reduction in the fine structure of the surface and a reduction in the fractal dimension of the topography. After a period of time, the surface is severely worn, the depth and width of the wave troughs increase, the detailed surface structure is destroyed, and the overall and local self-similarity weakens, and the fractal dimension decreases faster. Figure 9 shows the simulation results of the evolution of the contact area load ratio and contact load area ratio under different loads. The typical ground surface at 0 min in figure 6(a) was used for simulation. It can be seen from the figure that the initial surface morphology has a clear directionality and a uniform height distribution. The entrainment speed is 0.1 m s −1 , the sliding to rolling ratio is 0.4, and the wear coefficient is set to 0.018.

Effect of contact load on surface morphology and lubrication characteristics
It can be observed from figure 9 that the contact load ratio and the contact area ratio are not uniformly increased. In the initial stage of wear, the contact area increases slowly, and then increases sharply. When the wear reaches a certain level, the increase in contact load ratio and contact area ratio becomes slower, which is similar to the experimental phenomenon observed by Shirong et al [44]. This may be due to the presence of obvious high roughness peaks in the initial stage, resulting in slow increase in contact load ratio and contact area ratio. When these rough peaks are worn away, more small rough peaks are exposed, and the wear rate increases at this time. Finally, when the wear reaches a certain level, because there are more contacts, assuming the same wear rate, the increase in wear will slow down again. In addition, it can be seen from the figure 9 that when the wear reaches a certain level, the increase in load leads to a significant increase in the contact load ratio and the contact area ratio. This phenomenon may be due to the thicker film and more oil in the contact area under light load conditions, resulting in slow wear and slower increase in contact load ratio and contact area ratio. However, under heavy load conditions, due to the thin film thickness, the surface is severely damaged under the same cycle, and the lubrication state shifts from mixed lubrication to boundary lubrication, resulting in a significant increase in contact load ratio and contact area ratio [24]. Figure 10 shows the evolution of the surface fractal dimension under different loads.
It can be seen from the figure 10 that the higher the load, the more severe the decrease in the fractal dimension. The effect of acceleration at the beginning is not obvious. This may be because their surfaces are the same at first, but as the wear progresses, the greater the load, the greater the wear and the greater the difference between the rough surfaces. After 1500 cycles, under high load conditions, more rough peaks were worn, reflecting that the detailed structure and complexity of the worn surface were damaged more severely, and the effect of the load seemed very significant. Therefore, reducing the load is of great significance for reducing the degree of surface wear and increasing the working life of the equipment. Figure 11 shows the evolution of surface coefficient of friction COF and surface max flash temperature rise(MFTR) under different contact loads.
It can be seen from figure 11(a) that in the initial 1000 cycles, the contact load may have a small effect on the friction coefficient, because the area ratio and the contact load ratio increase slowly at this time. But with the passage of time, the load effect gradually became significant. After 1000 cycles, the friction coefficient of the case when the load was 1.82 GPa gradually increased to 0.12 or more, and then the rate of change became stable. The friction coefficient with a load of 0.59 GPa shows an approximately linear increase. The friction coefficient under mixed lubrication conditions is determined by the friction coefficient of the oil film and the friction coefficient of the contact. The change in the surface morphology has a significant effect on the friction coefficient. Figure 11(b) shows the evolution of the maximum flash temperature rise.  It can be seen that the MFTR of the high load condition is significantly higher than that of the light load condition. It should be pointed out that there is no obvious correlation between the MFTR and the wear process due to the introduction of wear, since the surface temperature rise is determined by the combination of lubricating oil, surface morphology, and matrix temperature. It is difficult to distinguish the relationship between the MFTR and the wear period. Figure 12 shows the simulation results of the evolution of the contact area load ratio and contact load area ratio under different velocity. The contact load is 0.59 GPa, the sliding to rolling ratio is 0.4, and the wear coefficient is set to 0.018.

Effect of sliding speed on surface morphology and lubrication characteristics
It can be seen from figure 12 that the sliding speed have a significant effect on the wear performance of the rough surface. When the sliding speed is 0.01 m s −1 , there is almost no oil film between the contact surfaces, and the two contact surfaces are mainly in contact. When the speed increases to 1 m s −1 , most areas of the two surfaces are separated by the oil film, and the surface wear becomes very slow. Figure 13 shows the evolution of the surface fractal dimension under different velocity.
As can be seen from figure 13, according to Archard's law, wear has a strong relationship with sliding speed. When the sliding speed is 0.01 m s −1 , because the sliding speed is small, the amount of surface wear is small, and the change in surface morphology is very small. When the speed increases to 1 m s −1 , the two contact surfaces are separated by an oil film, and the contact area is small, which leads to a reduction in wear and a small change in the fractal dimension of the surface. Only when the sliding velocity is large enough and the contact area is not separated by the oil film at 0.1 m s −1 , the large wear depth results in a relatively large change in the surface morphology. Therefore, under the current working conditions, the speed of 0.1 m s −1 should be avoided, and a higher speed  should be selected to increase the film thickness to improve the wear. Figure 14 shows the evolution of surface friction coefficient and MFTR under different sliding speed.
It can be seen from figure 14(a) that the friction coefficient decreases significantly with the increase of the sliding speed. And when the sliding speed is 0.1 m s −1 and 1m s −1 , the friction coefficient changes more obviously, which corresponds to the evolution of the surface fractal dimension with the progress of wear at different speed in figure 13. When the sliding speed is 0.1 m s −1 and 1m s −1 , the surface fractal dimension decreases more drastically, and the surface morphology changes more obviously, which leads to more drastic changes in the friction coefficient. As can be seen in figure 14(b), the highest flash temperature shows the opposite change law. When the sliding speed is 0.01 m s −1 , the contact area is all rough peak contact at this time, and the relative friction causes the MFTR to be very significant. When the speed is 1m/s, the surface is separated by the oil film, the rough peaks are less in contact, and the variation of MFTR is not obvious.

Conclusions
Based on the surface three-dimensional morphology fractal method and the deterministic mixed EHL wear model, this paper characterizes the evolution process of the wear surface morphology during the mixed lubrication process. First, the effectiveness of the lubrication model is verified through the optical interference oil film measurement experiment, and the wear coefficient of the actual surface is obtained through wear experiments. On this basis, the evolution  process of surface lubrication characteristics and fractal parameters in the process of sliding wear is investigated. Finally, the effect of contact load and sliding speed on the surface fractal parameters is discussed. The following conclusions can be drawn.
(1) In the initial wear stage, the surface morphology is ground texture, and the fine structure is more, the fractal dimension is larger. When the higher wave peaks on the surface are gradually worn away, it becomes a wider bearing plane, which causes the surface roughness to increase, the fine structure to decrease, and the fractal dimension of the morphology to decrease.
(2) Due to the dual effects of changes in lubrication status and wear, the increase in contact load results in a significant increase in contact load ratio and contact area ratio. And the higher the load, the more the fractal dimension decreases, that is, more rough peaks are worn out, and the detailed structure and complexity of the surface are destroyed more seriously.
(3) Both higher sliding speed and lower sliding speed will reduce the wear. Only when the two contact surfaces are not completely isolated by oil and there is a certain sliding speed, will there be obvious wear and the surface fractal dimension will be significantly reduced.