Non-Hertzian Elastohydrodynamic Contact Stress Calculation of High-Speed Ball Screws

: In order to eliminate the calculation error of the Hertzian elastohydrodynamic contact stress due to the asymmetry of the contact region of the helix raceway, a non-Hertzian elastohydrodynamic contact stress calculation method based on the minimum excess principle was proposed. Firstly, the normal contact stresses of the screw raceway and the nut raceway were calculated by the Hertzian contact theory and the minimum excess principle, respectively. Subsequently, the Hertzian solution and the non-Hertzian solution of the elastohydrodynamic contact stress could be determined by the Reynolds equation under different helix angles and screw speeds. Finally, the friction torque test of the double-nut ball screws was designed and implemented on a self-designed bed for validation of the proposed method. The comparison showed that the experimental friction torque was the good agreement with the simulated friction torque, which veriﬁed the effectiveness and correctness of the non-Hertzian elastohydrodynamic contact stress calculation method. Under the large helix angle, the calculation accuracy of asperity contact stress for the non-Hertzian solution was more accurate than that of the Hertzian solution at the contact region of ball screws. Therefore, the non-Hertzian elastohydrodynamic contact stress considering the asymmetry of the raceway contact region could more accurately analyze the wear depth of the high-speed ball screws.


Introduction
With the improvement of machining efficiency and machining accuracy of the NC machine tools, high-speed ball screws (BSs) are widely used in the machining field [1,2]. However, there is a calculation error problem of the Hertzian elastohydrodynamic contact stress due to the asymmetry of the contact region of the helix raceway. Thus, it is essential to investigate the non-Hertzian elastohydrodynamic contact stress of the BSs in these properties, such as kinematics, force balance, and elastohydrodynamic contact mechanics.
In terms of the force and deformation characteristics of the double-nut ball screws (DNBSs), many researchers studied the force balance equation and the torque balance equation of the working ball. Liu et al. [3] proposed a static load distribution model of the preloaded BSs considering geometric errors, which takes into account the influence of non-loaded balls and the interaction between the elastic deformation of a screw/nut and the Hertz contact forces of screw-ball/nut-ball contact areas. Lin et al. [4] established a low order static load distribution model for BSs that incorporates lateral deformation and geometric error effects to optimize the design parameters of BSs. Liu et al. [5] presented a novel method to analyze the static load distribution of BSs, which takes into account the nut position variation. The effects of axial load and original contact angle on load distribution are obtained. Huang et al. [6] used the machine learning method to diagnose the preload state of the BSs by monitoring the mechanical signals, such as the driving motion nut are in a fully preloaded state due to the contact deformation of the left nut and right nut, respectively ( Figure 1a). The static normal contact force at points A and B is expressed as  (2) where QA and QB denote the normal contact force of the ball at the screw contact point A and the nut contact point B, respectively. FSA and FSB indicate the friction of the ball at the screw contact point A and the nut contact point B, respectively. FIH is the inertia force. MSA and MSB are the friction torque of the ball at the screw contact point A and the nut contact point B, respectively. MIH is the inertia torque. The parameters of DNBSs in the operating state, such as βA, βB, QA, QB, FSA, and FSB, are obtained by Equation (2) using the Newton iteration method. The specific solution process of the parameters of DNBSs can be observed in the reference [20]. The structural parameters of DNBSs are shown in Table 1.  The DNBSs is operating at a constant speed; the working balls are in the balanced state of the force and torque, which comprise normal contact force, friction, inertia force, friction torque, and inertia torque. The balance equations of the force and torque of the ball are determined as where Q A and Q B denote the normal contact force of the ball at the screw contact point A and the nut contact point B, respectively. F SA and F SB indicate the friction of the ball at the screw contact point A and the nut contact point B, respectively. F IH is the inertia force. M SA and M SB are the friction torque of the ball at the screw contact point A and the nut contact point B, respectively. M IH is the inertia torque. The parameters of DNBSs in the operating state, such as β A , β B , Q A , Q B , F SA , and F SB , are obtained by Equation (2) using the Newton iteration method. The specific solution process of the parameters of DNBSs can be observed in the reference [20]. The structural parameters of DNBSs are shown in Table 1.

Non-Hertzian Normal Contact Stress Calculation
Due to the asymmetry of the contact region of the helical raceway of the high-speed BSs, there is a big error in the calculation of the contact stress distribution and contact region in the screw raceway contact point A and the nut raceway contact point B by using the traditional Hertzian contact theory. The non-Hertzian contact region between the ball and raceway is inclined compared with the Hertzian contact region, which is shown in Figure 1b. Based on the minimum excess principle, a method was proposed to analyze and calculate the normal contact stress between the ball and the raceway of the high-speed DNBSs.
In terms of the numerical solution process of the contact stress distribution between the ball and raceway, the contact region and deformation compatibility equation should be discretized. The contact gap and contact stress of the two contact surfaces should satisfy the following equation [21].
where c is the contact surface gap matrix. c i,j is the contact surface gap distribution in the mesh region, which is the element of the contact surface gap matrix c. d 0 denote the approaching distance of the contact surface. c g is the contact surface geometry matrix. c gi,j indicates the contact surface geometry in the mesh region, which is the element of the contact surface geometry matrix c g . c r is the contact surface roughness matrix. c ri,j is the contact surface roughness distribution, which is the element of the contact surface roughness matrix c r . V is the surface elastic deformation matrix. V i,j is the surface elastic deformation distribution in the mesh region, which is the element of the surface elastic deformation matrix V. I and J indicate the mesh grid number in the x and y directions, respectively. The discretized contact deformation equation is expressed as where p i,j is the element of the contact stress matrix P. K i − k,j − l is the element of deformation influence coefficient matrix K.
The boundary condition of the contact equation is expressed as In the contact region, the boundary condition of the contact equation is c i,j = 0, p i,j ≥ 0; Out of the contact region, the boundary condition of the contact equation is c i,j > 0, p i,j = 0. According to the above two conditions, c i,j , and p i,j cannot be both zero, thus The solution of the contact stress and the contact geometry is the key to the contact equation. The contact equation between the ball and raceway can be expressed by a matrix equation: where ⊗ indicates cross product. According to the variational principle, the contact problem described by Equations (7)-(9) can be transformed to the conditional extremum problem of quadratic function, which is derived by: The boundary condition of Equation (10) is given by: The fast convergence of the contact Equation (11) in the cyclic iteration is achieved by using the conjugate gradient method. The flowchart is shown in Figure 2. The detailed derived process and iteration solution process of the non-Hertzian contact stress is listed as follows [22]: The initial normal contact stress P is given, p i,j ≥ 0. The normal contact stress distribution p i,j is also satisfied in the mesh region.

2.
The discretized normal contact deformation V is determined by Equation (4).

3.
The contact gap c between the two contact surfaces is calculated by Equations (3)-(6), and the contact gap is satisfied with the boundary condition.

4.
According to the variational principle, the contact in Equation (7) can be transformed by the conditional extremum of the quadratic function, the normal contact stress p i , j is modified by Equation (10) satisfying the boundary condition in Equation (11). 5.
The current relative error is determined by: 6.
If ε ≥ ε 0 , the iteration solution process of the non-Hertzian contact stress continues from steps (2)-(5), otherwise the iteration solution process ends.

() 2
T T W = + P c P P KP (10 The boundary condition of Equation (10) is given by: The fast convergence of the contact Equation (11) in the cyclic iteration is achieve by using the conjugate gradient method. The flowchart is shown in Figure 2. The detaile derived process and iteration solution process of the non-Hertzian contact stress is liste as follows [22]: 1. The initial normal contact stress P is given, pi,j ≥ 0. The normal contact stress distribu tion pi,j is also satisfied in the mesh region. 2. The discretized normal contact deformation V is determined by Equation (4). 3. The contact gap c between the two contact surfaces is calculated by Equations (3)-(6 and the contact gap is satisfied with the boundary condition. 4. According to the variational principle, the contact in Equation (7) can be transforme by the conditional extremum of the quadratic function, the normal contact stress p is modified by Equation (10) satisfying the boundary condition in Equation (11). 5. The current relative error is determined by:  The classic solution between the ball and smooth plane can be determined by Hertz ian contact theory [23]; it is expressed as The comparison between the Hertzian solution and non-Hertzian solution of the con tact stress under the helix angle is zero is shown in Figure 3. As shown in Figure 3, th non-Hertzian solution is in good agreement with the Hertzian solution in the contac The classic solution between the ball and smooth plane can be determined by Hertzian contact theory [23]; it is expressed as The comparison between the Hertzian solution and non-Hertzian solution of the contact stress under the helix angle is zero is shown in Figure 3. As shown in Figure 3, the non-Hertzian solution is in good agreement with the Hertzian solution in the contact stress calculation results, which is verified the correctness of the non-Hertzian contact stress calculation method.
Appl. Sci. 2021, 11,12081 stress calculation results, which is verified the correctness of the non-Hertzian stress calculation method.

Non-Hertzian Elastohydrodynamic Contact Stress Calculation
The viscosity flow behavior of the lubricant film between the ball and racew be described by the Reynolds equation shown in Figure 4b. The initial film contac distribution in the Reynolds equation is the non-Hertzian contact stress distri which can further observe the variation of the elastohydrodynamic lubrication cha istics.
The assumptions of the Reynolds equation in the analysis are as follows: 1. The lubricant and the solid are isothermal, and the contact stress and the film ness do not vary with time. 2. The curvature radius of the solid surface is much larger than the lubricant film ness. 3. No relative sliding between the lubricant film and solid surface at the commo face. 4. Compared with the film shear stress, the inertial force and other bulk forces lubricant film are negligible. 5. Due to the lubricant film thickness being very thin, it can be assumed that t contact stress remains constant along the film thickness direction, namely no s film effect. 6. The inlet of the contact region is flooded fully by the lubricant film. 7. The lubricant is the Newtonian fluid and obeys Newton's law of viscosity.
The coupling relationship of the elastohydrodynamic contact stress distribu film thickness distribution h, film viscosity η, film density ρ, and film entrainment v us could be accurately determined by the Reynolds equation [24], which is written where p denotes the non-Hertzian elastohydrodynamic contact stress in the con gion. u1 and u2 indicate the linear velocity of the ball and the raceway in the tan direction (x-direction) of the helix raceway track, respectively. The entrainment v of the lubricant film is us = (u1 + u2)/2. The output parameters of point A and point B lated by Equation (2), such as QA, βA, QB, and βB, are substituted into the Reynolds Eq

Non-Hertzian Elastohydrodynamic Contact Stress Calculation
The viscosity flow behavior of the lubricant film between the ball and raceway can be described by the Reynolds equation shown in Figure 4b. The initial film contact stress distribution in the Reynolds equation is the non-Hertzian contact stress distribution, which can further observe the variation of the elastohydrodynamic lubrication characteristics.  In terms of the elastohydrodynamic lubrication theory [25], the film thickness distribution is expressed as The lubricant film viscosity proposed by Roelands [26] is written as The lubricant film density from Dowson and Higginson [27] is shown in Equation (17).
where h0 indicates the central film thickness in the contact region. η0 and ρ0 denote the The assumptions of the Reynolds equation in the analysis are as follows: 1.
The lubricant and the solid are isothermal, and the contact stress and the film thickness do not vary with time.

2.
The curvature radius of the solid surface is much larger than the lubricant film thickness.

3.
No relative sliding between the lubricant film and solid surface at the common interface.

4.
Compared with the film shear stress, the inertial force and other bulk forces of the lubricant film are negligible.

5.
Due to the lubricant film thickness being very thin, it can be assumed that the film contact stress remains constant along the film thickness direction, namely no squeeze film effect. 6.
The inlet of the contact region is flooded fully by the lubricant film.

7.
The lubricant is the Newtonian fluid and obeys Newton's law of viscosity.
The coupling relationship of the elastohydrodynamic contact stress distribution p, film thickness distribution h, film viscosity η, film density ρ, and film entrainment velocity u s could be accurately determined by the Reynolds equation [24], which is written as where p denotes the non-Hertzian elastohydrodynamic contact stress in the contact region. u 1 and u 2 indicate the linear velocity of the ball and the raceway in the tangential direction (x-direction) of the helix raceway track, respectively. The entrainment velocity of the lubricant film is u s = (u 1 + u 2 )/2. The output parameters of point A and point B calculated by Equation (2), such as Q A , β A , Q B , and β B , are substituted into the Reynolds Equation (14), the elastohydrodynamic characteristics of point A and point B can be determined, respectively. In terms of the elastohydrodynamic lubrication theory [25], the film thickness distribution is expressed as The lubricant film viscosity proposed by Roelands [26] is written as The lubricant film density from Dowson and Higginson [27] is shown in Equation (17).
where h 0 indicates the central film thickness in the contact region. η 0 and ρ 0 denote the initial value of the film viscosity and density, respectively. ROU(x,y) indicates the roughness of the contact surface calculated by the fractal theory [28] shown in Figure 4a. The major parameters of the rough fractal rough surface are the fractal dimension D and scale parameter G, D = 2.47, G = 0.41. The rough fractal surface is the contact region of the raceway. The hardness range of the ball is 62-64 HRC higher than that of the screw raceway 58-62 HRC; The processing technology of the working ball is mature, and the surface of the working ball is smooth. According to the above two reasons, the ball contact surface is regarded as a smooth surface. The surface roughness of the working ball is much less than that of the screw raceway. Therefore, the surface roughness of the ball has little effect on the calculation efficiency of this model. Due to the film contact stress being largely influenced by the surface roughness, the contact surface roughness should be considered in the elastohydrodynamic contact system of BSs. In terms of the mature surface manufacturing technology of the ball, the ball contact surface can be seen as the smooth surface, and the raceway contact surface is the rough fractal surface (Figure 4a).
The boundary conditions of the Reynolds equation are written as The load balance condition of the Reynolds equation is expressed as Substituting Equations (15)- (19) into Equation (14), the film contact stress distribution and the elastohydrodynamic contact stress distribution can be determined accurately by using the finite difference method [29,30], please see the Appendix A.

Experiments and Verification
To validate the correctness of elastohydrodynamic contact characteristics of the DNBSs by using the non-Hertzian contact stress calculation method. The friction torque of the DNBSs is measured by the self-designed bed, and the experimental friction torques are in good agreement with the simulated friction torques, which verifies the correctness of the non-Hertzian elastohydrodynamic contact stress calculation method.

Experiments Procedures
The experimental setup of the DNBSs is shown in Figure 5. The DNBSs was only applied with the preload in the experiment. The tension-compression sensor was mounted on the nut by the dowel bar. The friction torque of the BSs was measured under the experimental condition with constant preload F pre = 1500 N, constant lubricating oil viscosity η 0 = 0.05 Pa·s, and the different screw speeds. The simulated friction torque of the BSs can be expressed as where F SBt and F SBb represent friction components at the directions of t and b in the Frenet-Serret coordinate system Otnb at point B, respectively [31]. n bb denotes the number of working balls. where xo and yo are the entrance coordinates of the contact region of BSs, respectively, and xe and ye are the exit coordinates of the contact region of BSs, respectively. The values of these coordinate parameters have been given: xo = yo = −277 μm, xe = ye = 277 μm.
The load balance condition of the Reynolds equation is expressed as ( ) , d d p x y x y Q  =  (19) Substituting Equations (15)- (19) into Equation (14), the film contact stress distribution and the elastohydrodynamic contact stress distribution can be determined accurately by using the finite difference method [29,30], please see the Appendix A.

Experiments and Verification
To validate the correctness of elastohydrodynamic contact characteristics of the DNBSs by using the non-Hertzian contact stress calculation method. The friction torque of the DNBSs is measured by the self-designed bed, and the experimental friction torques are in good agreement with the simulated friction torques, which verifies the correctness of the non-Hertzian elastohydrodynamic contact stress calculation method.

Experiments Procedures
The experimental setup of the DNBSs is shown in Figure 5. The DNBSs was only applied with the preload in the experiment. The tension-compression sensor was mounted on the nut by the dowel bar. The friction torque of the BSs was measured under the experimental condition with constant preload Fpre = 1500 N, constant lubricating oil viscosity η0 = 0.05 Pa·s, and the different screw speeds. The simulated friction torque of the BSs can be expressed as where FSBt and FSBb represent friction components at the directions of t and b in the Frenet-Serret coordinate system Otnb at point B, respectively [31]. nbb denotes the number of working balls.

Experimental Results
The friction coefficient [32] of BSs can be written as

Experimental Results
The friction coefficient [32] of BSs can be written as where Q is the total normal contact force, including the film bearing force Q f , and the asperity bearing force Q a . The ratio of the asperity bearing stress is R S , R S = Q a /Q. µ a is the dry friction coefficient and µ a = 0.1. τ f is the film shear stress proposed by the Eyring model [33]. A h is the film hydrodynamic contact area, η is the film viscosity, ∆u is the relative sliding velocity, and ∆u = u 1 − u 2 . u 1 and u 2 are the same as that of Equation (14). h 0 is the film central thickness.
The results of the friction torque test are shown in Figure 6. The friction coefficient curve of the DNBSs is indicated by the blue lines in Figure 6b, which is in good agreement with the schematic Stribeck curve [34] shown in the green box (Figure 6a). The friction coefficient is composed of the asperity friction coefficient and the hydrodynamic lubrication friction coefficient. When the screw speed is lower than 502.65 mm/s, the BSs is in the mixed lubrication state. When the screw speed is higher than or equal to 502.65 mm/s, the BSs is in the hydrodynamic lubrication state. The comparison of the simulated and experimental friction torque is shown in Figure 6c. The comparisons showed that the experimental friction torque was consistent with the simulated friction torque, which verified the correctness of the non-Hertzian elastohydrodynamic contact stress calculation method.

Comparison of the Normal Contact Stress between Hertzian Solution and Non-Hertzian Solution
In order to determine the contact stress distribution of DNBSs accurately, the contact

Comparison of the Normal Contact Stress between Hertzian Solution and Non-Hertzian Solution
In order to determine the contact stress distribution of DNBSs accurately, the contact stress distributions between the ball and raceway are calculated by the minimum excess principle and Hertzian contact theory, respectively. The calculation results by the two methods are analyzed comparatively.
The contact stress distribution of DNBSs is calculated under the preload F pre is 1500 N, and the helix angle α is 21.7 • . The contact stress comparison between the Hertzian solution and non-Hertzian solution of BSs is shown in Figure 7. As shown in Figure 7, the difference contour maps between the Hertzian solution and non-Hertzian solution at contact point A and contact point B indicate the great difference between the two calculation methods. The difference of the contact stress in the edge of the contact region is the largest, and the one in the center of the contact region is the smallest. Thus, there is a big error of contact stress distribution calculated by the Hertzian contact theory, and the contact stress distribution calculated by the non-Hertzian solution is more accurate.

Analysis of the Non-Hertzian Contact Stress with Different Helix Angles
In order to further analyze the effect of contact stress error calculated by the Hertzian contact theory with the variation of the helix angle, the contact stress difference contour maps at the screw raceway contact point A under the helix angle α are 5.68°, 13.97°, and 21.71° are determined, respectively. The elliptical contact region of non-Hertzian solution PA at the contact point A generates a deflection angle θs shown in Figure 8a. The contact stress distribution difference PA-PAH between the non-Hertzian solution PA and the Hertzian solution PAH in the edge of the contact region is largest, and the one in the center of the contact region is smallest shown in Figure 8b.

Analysis of the Non-Hertzian Contact Stress with Different Helix Angles
In order to further analyze the effect of contact stress error calculated by the Hertzian contact theory with the variation of the helix angle, the contact stress difference contour maps at the screw raceway contact point A under the helix angle α are 5.68 • , 13.97 • , and 21.71 • are determined, respectively. The elliptical contact region of non-Hertzian solution P A at the contact point A generates a deflection angle θ s shown in Figure 8a. The contact stress distribution difference P A -P AH between the non-Hertzian solution P A and the Hertzian solution P AH in the edge of the contact region is largest, and the one in the center of the contact region is smallest shown in Figure 8b.

Analysis of the Elastohydrodynamic Contact Stress with Different Helix Angles
The elastohydrodynamic characteristics of point A under the preload are 1500 N and the entrainment velocity is 1005.3 mm/s, which are shown in Figure 10. When the lubricant film is entrained into the contact region, the secondary peak of the film contact stress in the exit coordinates of the contact region is generated, as shown in Figure 10a. Figure 10b shows the film thickness distribution of point A.

Analysis of the Elastohydrodynamic Contact Stress with Different Helix Angles
The elastohydrodynamic characteristics of point A under the preload are 1500 N and the entrainment velocity is 1005.3 mm/s, which are shown in Figure 10. When the lubricant film is entrained into the contact region, the secondary peak of the film contact stress in the exit coordinates of the contact region is generated, as shown in Figure 10a. Figure 10b shows the film thickness distribution of point A.

Effect of the Elastohydrodynamic Contact Stress under Different Screw Speeds
The central elastohydrodynamic contact stresses of the BSs under different helix angles and different screw speeds are shown in Figure 12. As shown in Figure 12a, the non-Hertzian solution is smaller than the Hertzian solution of the central elastohydrodynamic contact stress at contact point A. The central elastohydrodynamic contact stress decreases with the increase in the helix angle and the screw speed, respectively. As shown in Figure 12b, the non-Hertzian solution is greater than the Hertzian solution of the central elastohydrodynamic contact stress at contact point B. The central elastohydrodynamic contact stress increases as the helix angle increases, and decreases as the screw speed increases.
stress. The asperity contact stresses in the elliptical contact region of the BSs under different helix angles and different screw speeds are shown in Figure 13. The law of the asperity contact stress of points A and B is the same as that of the elastohydrodynamic contact stress of points A and B. Under the helix angle of 21.71°, the calculation accuracies of asperity contact stress for the non-Hertzian solution are 2.93% and 3.01% more accurate than that of the Hertzian solution at the screw-raceway contact region and the nut-raceway contact region, respectively.  The wear depth and wear rate of the BSs are mainly influenced by the asperity contact stress. The asperity contact stresses in the elliptical contact region of the BSs under different helix angles and different screw speeds are shown in Figure 13. The law of the asperity contact stress of points A and B is the same as that of the elastohydrodynamic contact stress of points A and B. Under the helix angle of 21.71 • , the calculation accuracies of asperity contact stress for the non-Hertzian solution are 2.93% and 3.01% more accurate than that of the Hertzian solution at the screw-raceway contact region and the nut-raceway contact region, respectively.

Conclusions
The central elastohydrodynamic contact stresses of the BSs under different helix angles and different screw speeds are shown in Figure 12. As shown in Figure 12a, the non-Hertzian solution is smaller than the Hertzian solution of the central elastohydrodynamic contact stress at contact point A. The central elastohydrodynamic contact stress decreases with the increase in the helix angle and the screw speed, respectively. As shown in Figure  12b, the non-Hertzian solution is greater than the Hertzian solution of the central elastohydrodynamic contact stress at contact point B. The central elastohydrodynamic contact stress increases as the helix angle increases, and decreases as the screw speed increases.
The wear depth and wear rate of the BSs are mainly influenced by the asperity contact stress. The asperity contact stresses in the elliptical contact region of the BSs under different helix angles and different screw speeds are shown in Figure 13. The law of the asperity contact stress of points A and B is the same as that of the elastohydrodynamic contact stress of points A and B. Under the helix angle of 21.71°, the calculation accuracies of asperity contact stress for the non-Hertzian solution are 2.93% and 3.01% more accurate than that of the Hertzian solution at the screw-raceway contact region and the nut-raceway contact region, respectively.

Conclusions
Due to the obvious asymmetry of the raceway contact region of the high-speed BSs, a non-Hertzian elastohydrodynamic contact stress calculation method based on the minimum excess principle is proposed, and the conclusions below can be drawn: (1) The normal contact stress error between the Hertzian solution and the non-Hertzian solution increases as the helix angle, and the contact stress error of the screw raceway contact region is always greater than that of the nut raceway contact region.
(2) The non-Hertzian elastohydrodynamic contact stress of the BSs show different trends in the contact regions of the screw raceway and nut raceway with the increase of the helix angle and screw speed, respectively. The non-Hertzian solution is smaller than the Hertzian solution of the central elastohydrodynamic contact stress at the screw raceway contact region. However, the non-Hertzian solution is greater than the Hertzian solution of the central elastohydrodynamic contact stress at the nut raceway contact region.
(3) The positioning accuracy of the DNBSs is obviously influenced by the wear of the raceway due to the asperity contact problem. Under the helix angle of 21.71 • , the calculation accuracies of asperity contact stress for the non-Hertzian solution are 2.93% and 3.01% more accurate than that of the Hertzian solution at the screw raceway contact region and the nut raceway contact region, respectively. Therefore, the non-Hertzian elastohydrodynamic contact stress calculation method is essential to precisely analyze the wear mechanism of the high-speed BSs.
Author Contributions: T.S. proposed and simulated the calculation method, designed the experiment, and wrote the paper. X.G. implemented the experiment. M.W. guided the general idea of this paper. Y.Z. analyzed the experimental data. All authors have read and agreed to the published version of the manuscript.

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

The Finite Difference Method
The finite difference method is an approximation method for solving the numerical solution for differential equations. The main principle of the finite difference method is to make a direct difference approximation to the differential term in the differential equation, thus transforming the differential equation into the algebraic equation for the solution. Due to the advantages of conceptual clarity, generality, and ease of implementation on computers, the finite difference method is widely used in the field of differential equation solving. Therefore, the Reynolds equation describing the elastohydrodynamic contact characteristics of the lubricant film of DNBSs under different operating conditions is determined by the finite difference method in this paper.