Abstract
The actual contact point of a spiral bevel gear pair deviates from the theoretical contact point due to the gear deformation caused by the load. However, changes in meshing characteristics due to the migration of contact points are often ignored in previous studies on the elastohydrodynamic lubrication (EHL) analysis of spiral bevel gears. The purpose of this article is to analyze the impact of contact point migration on the results of EHL analysis. Loaded tooth contact analysis (LTCA) based on the finite element method is applied to determine the loaded contact point of the meshing tooth pair. Then, the osculating paraboloids at this point are extracted from the gear tooth surface geometry. The geometric and kinematic parameters for EHL simulation are determined according to the differential geometry theory. Numerical solutions to the Newtonian isothermal EHL of a spiral bevel gear pair at the migrated and theoretical contact points are compared to quantify the error involved in neglecting the contact point adjustment. The results show that under heavy-loaded conditions, the actual contact point of the deformed gear pair at a given pinion (gear) roll angle is different from the theoretical contact point considerably, and so do the meshing parameters. EHL analysis of spiral bevel gears under significant load using theoretical meshing parameters will result in obvious errors, especially in the prediction of film thickness.
1. Introduction
Spiral bevel gears, the most common form of gearing used for transmitting significant power between nonparallel shafts, have the potential to wear, pitting, and scuffing due to lubrication failure. Therefore, the lubrication behaviour of these gears has long been the concern of many research studies. However, fruitful results have only been yielded after the 21st century thanks to the significant improvement in tooth contact analysis technology, elastohydrodynamic lubrication (EHL) theory, numerical algorithms, and computer performance.
Early studies by Simon [1] and Chao [2] analyzed the lubricating properties of spiral bevel and hypoid gears assuming smooth surface due to the insufficient computational resource and techniques at that time. Recently, Xu and Kahraman [3] predicted the transmission efficiency of hypoid gear pairs, where the friction coefficient is calculated by applying a line contact EHL model which is only reasonable for the theoretical conjugate tooth surfaces. In applications, however, it is imperative to modify the tooth profiles of the mating gears to make them mismatched, thereby being capable of tolerating a certain amount of misalignments caused by the deflection of the supporting shaft and housing, as well as installation and machining errors. Due to the mismatch, an elliptical contact footprint, instead of the theoretical line contact one, is formed between each pair of teeth in contact. Kolivand et al. [4] regressed the friction coefficient formulae for the entire lubrication range, while the line contact assumption was still adopted. Simon [5] reported a thermo-EHL analysis of spiral bevel and hypoid gears where although a more realistic elliptical point contact model is applied, the entrainment flow vector is assumed to be along the minor axis of the elliptical footprint. Similarly, Zhang et al. [6] also predicted the tooth surface friction coefficient and wear through thermo-EHL analysis with the assumption that the rolling velocity direction coincides with the minor axis of the contact ellipse. Mohammadpour et al. [7] provided an isothermal Newtonian EHL solution of hypoid gear pairs at high loads, considering angled entrainment flow of the lubricant. This work was then extended to take non-Newtonian shear of thin films into account [8]. Lately, the unified lubrication model covering all the lubrication regions [9] was recognized by more and more researchers in the field. Pu et al. [10] used it to conduct the mixed EHL analysis of spiral bevel and hypoid gears with machined roughness. Cao et al. [11] also employed this model to study the effect of contact path on the lubrication performance, friction, and fatigue life. More recently, Gan et al. [12] further applied this model to investigate the temperature behaviour under mixed lubrication condition.
The data required for EHL analysis include the load acting on the meshing tooth pair and the geometric and kinematic meshing parameters at the contact points. In studies mentioned above, the load sharing among simultaneously engaged tooth pairs was evaluated using the loaded tooth contact analysis (LTCA), and the meshing parameters were calculated at the theoretical contact point (TCP) obtained through tooth contact analysis (TCA) which is, however, only valid under the no-load and light-loaded conditions. Many researchers such as Wang [13], Zhou et al. [14], and Sun et al. [15] have revealed that the contact path (which consists of all contact points over a tooth engagement cycle) of the mating spiral bevel gear pair under loaded condition is different from the theoretical one due to deformation of gear teeth and the body caused by the load. It turns out that the contact path changed by adjusting the mounting position of mating gears has a nonnegligible impact on lubrication performance [11]. However, changes in the contact path due to the load-induced deformation are usually ignored in EHL analysis.
The present study proposes a method for calculating the meshing parameters at the loaded contact point (LCP) which is synthesized from the load distribution results obtained by LTCA. Once the LCP is obtained, the osculating paraboloids at that point are extracted from the geometry of pinion and gear tooth surfaces. Then, the geometric and kinematic parameters for EHL simulation are determined according to the differential geometry theory. The flowchart of the proposed methodology is shown in Figure 1. Numerical solutions to the Newtonian isothermal EHL of a Gleason-type spiral bevel gear pair at the loaded and theoretical (unloaded) contact points are obtained and compared to quantify the error involved in neglecting the contact point adjustment due to deformation caused by the load.
2. Determination of Meshing Parameters
TCA algorithms reported in the literatures [16, 17] and many others can be applied to determine the TCPs of spiral bevel gear pairs under no-load and light-loaded conditions where the deformation caused by the load can be ignored. At high loads, the obtainment of LCPs needs LTCA. The LCP, in the present work, is defined as the contact point of the undeformed mating tooth surfaces of the pinion and the gear as shown in Figure 2. It is known that, for a pair of mismatched spiral bevel gears, an elliptical contact footprint is formed between the contacting tooth surfaces due to elastic contact deformation. The geometric center of the contact ellipse can be viewed as the projection point of the LCP on the deformed tooth flanks (Figure 2). Furthermore, the center of the contact pressure distribution, which is readily available from LTCA, is close to the center of the contact ellipse. Therefore, the center of the contact pressure distribution is utilized to determine the LCP in the present work. The center of the contact pressure is sometimes defined as the acting point of the effective meshing force due to contact pressure and named the “effective meshing point” (EMP) [13, 15], which will also be used in the following description.
The EMPs at each meshing position over a tooth engagement cycle (from a pair of teeth coming into contact to exit contact) can be obtained through quasistatic LTCA based on either a conventional finite element method [14] or a finite element/contact mechanism method [18]. For the second method, readers may refer to references [13, 19], where a specialized LTCA tool [20] is utilized. In this work, the first method is adopted, and LTCA is performed with the aid of the commercial FEA software ABAQUS [21].
The coordinate systems shown in Figure 3 are established to describe the relative kinematics of a spiral bevel gear pair. It should be noted that the proposed method is also applicable to the hypoid gear pair having a shaft offset. The global fixed reference frame whose origin is fixed to the intersection point of the pinion and gear rotation axes is established automatically in the simulation environment of ABAQUS, and the results from LTCA are given in this coordinate system. and are the local fixed reference frames used to describe the rotational positions of the gear and the pinion, respectively. The gear and the pinion rotate about the -axis of their local fixed coordinates. is the angle between the gear and the pinion axis. and are the body-fixed coordinate systems of the gear and the pinion in which their tooth surfaces are represented. The symbol will be often used in this article where the subscript indicates the coordinate system in which is represented and the superscript indicates the object to which belongs.
The procedures of using ABAQUS for LTCA have been introduced in detail by Zhou [14] and Sun et al. [15]. The position vector of EMP can be obtained conveniently as long as the history output XN (the center of the total force due to contact pressure [21]) of the contacting tooth pair is requested. Besides, the effective meshing force vectors that will be used to determine the position of tooth surfaces after deformation can be acquired from the history output CFN (the total force due to contact pressure [21]).
The next focus of the present study is the calculation of the coordinates of LCPs corresponding to the EMPs. As shown in Figure 4(a), it is assumed that is the EMP at a given meshing position and its position vector represented in the fixed coordinate system is . is the vector of effective mesh force acting on the EMP. can first be projected onto the projection plane of the gear and pinion teeth as shown in Figure 4(c). As illustrated in Figure 4(b), the coordinates of the projected points for the pinion () and the gear () can be calculated by
(a)
(b)
(c)
The mathematical model of the tooth flank geometry depends on how the gear teeth are machined but can be generally represented as [16]where indicates the cutting surface formed by the rotation of the cutting edge about the axis of the cutter head, is the parameter describing the cutting edge geometry, is the cutter phase angle, and is the transformation matrix representing the coordinate transformation from the coordinate system fixed to the cutter head to the coordinate system fixed to the gear blank. The form of the matrix hinges on the motion of the machine tool. Readers can refer to relevant literatures (e.g., [16, 17]) for the definition of the matrix . The cutting surface will generate a family of surfaces in the coordinate system rigidly connected to the gear blank. is the parameter of the family of surfaces and can either be the blank phase angle or the cradle angle. According to the enveloping theorem, the tooth surface of the machined gear is the envelope surface of the family of surfaces . The points on which simultaneously satisfy the equation of meshing are those on the envelope surface (namely, the tooth surface). In the equation of meshing, denotes the unit normal vector of the cutting surface and is the velocity in relative motion between the cutting surface and the tooth surface.
Applying equation (2) and the coordinate relationship between the tooth surface point and the projection point (as depicted in Figures 3 and 4(c)), the LCPs of the pinion and the gear can be calculated, respectively, bywhere , , and are the three components of the position vector of the LCP represented in the body-fixed frames of the pinion and the gear ().
The osculating paraboloids at LCP are extracted from the tooth surfaces of the pinion and the gear by the following method. In order to facilitate the following description, the equation of the gear tooth surface can be represented in a two-parameter form as follows:where is derived from the equation of meshing .
As mentioned previously, the LCP is different from the theoretical unloaded contact point calculated based on the principle of local conjugation [16]. That is to say, the contact conditions are not satisfied at the LCP for the undeformed gear set. As illustrated in Figure 5(a), the LCPs on the pinion () and gear () tooth surfaces do not coincide and their normals ( and ) are not collinear. If the gear and pinion tooth surfaces are in contact at the LCP, it is most likely that the relative position of the gear and pinion tooth flanks changes due to the deformation of the gear set caused by the load. The total deformation of the gear set is mainly composed of the tooth surface contact deformation, the bending and shear deflection of the gear teeth, and the bending and torsional deformation of the gear body and the shaft, as well as the deformation of housing. The rotation of the gear and pinion teeth is one of the manifestations of the total deformation [22]. They are, respectively, simulated by the rotation angles of the gear () and the pinion () around their rotational shafts. According to the calculation process of the LCP described before, it is obvious that and will coincide after the rotation as illustrated in Figure 5(b). The gear and pinion tooth surfaces are now represented in the coordinate system as follows:and at the LCP, there is
(a)
(b)
(c)
Usually, , and thus the transformation matrices in equations (5) and (6) are defined as
Then, the rotation angle of the gear tooth () and the pinion tooth () are determined by simultaneously solving equations (5)–(7) at the LCP.
At this time, the tooth surfaces of the gear and the pinion have coinciding position vectors at the LCP, while the contact condition is still not satisfied because the surface normal vectors at the LCP ( and ) are not collinear as illustrated in Figure 5(b). It is difficult to describe the changing process of the tooth surface position due to the composite deformation of the gear set. Nevertheless, it is reasonable to assume that the deformation makes the normal vectors of the gear and pinion tooth surfaces at the LCP coincide with the direction of the effective mesh force obtained by LTCA. Based on this assumption, the next step is to eliminate the angle between the tooth surface normal vector and the direction vector of the effective mesh force. It is assumed that the to-be-eliminated angles for the gear and the pinion are and , respectively (Figure 5(c)). Taking the pinion as an example, the procedure of eliminating the angle is described. A local coordinate system which is rigidly attached to the pinion tooth surface at the point is established (Figure 5(c)). The unit vectors along the axes of are defined aswhere the three components of the direction vector of the effective mesh force are obtained by
Before rotating, the pinion tooth surface should be represented in the local body-fixed coordinate system :and the coordinate transformation matrix from to iswhere , , and are the unit vectors along the axes of the coordinate system and , , and are the three components of the position vector of represented in . Then, the pinion tooth surface along with the coordinate system rotates about the axis to the position shown by . The rotation matrix is
The normal vector of the pinion tooth surface at the LCP is now parallel with the direction vector of the effective mesh force. A similar approach can be used to rotate the gear tooth flank to make its normal vector collinear with the direction vector of the effective mesh force as well. At this point, the tooth surfaces of the gear and the pinion meet the contact conditions at the LCP, and locally, they are osculating paraboloids at the LCP.
The tooth surfaces of the gear and the pinion should be represented in the same coordinate system, e.g., , and then the geometric characteristic parameters at the LCP can be determined according to the differential geometry principle [23]. Note that the representation of the tooth surface remains unchanged in the body-fixed coordinate system, and thereby . Therefore, the tooth surfaces of the gear and the pinion in are defined by
As depicted in Figure 6(a), it is assumed that the position vectors of the LCP for the pinion and gear tooth surfaces represented in are and , and the corresponding normal vectors of tooth surfaces at the LCP are and . The induced normal curvature of two conjugate surfaces, which is defined as the difference between the normal curvatures of the two surfaces in arbitrary direction, describes the closeness of the two surfaces along the specified direction near the contact point. The smaller the value of the induced normal curvature, the closer the two surfaces. Therefore, the direction in which the induced normal curvatures are minimum and maximum can be, respectively, defined as the direction of the major and minor axes of the instantaneous contact ellipse.
(a)
(b)
To obtain the minimum value of the induced normal curvature, the principal curvatures and principal directions at the LCP can be calculated firstly. The normal curvature in the arbitrary direction is given bywhere and are the first and second fundamental forms of the surfaces, the coefficients (, , and ) of the first fundamental form are given asand the coefficients (, , and ) of the second fundamental form are defined aswhere the unit normal to the surface is calculated by
Equation (15) can be rewritten as follows:where indicates the tangential direction of a surface corresponding to the normal curvature. Since the principal curvature is the maximum or minimum of the normal curvature, taking the derivative of equation (19) respect to and letting , one can obtain the relationship between the principal curvatures and principal directions as
The equation of principal curvatures can be derived by eliminating from equation (20):
The calculated principal curvatures can be substituted into equation (19) to evaluate the corresponding principal directions. After getting the principal curvatures of both the gear ( and ) and pinion ( and ) tooth surfaces, as well as the corresponding principal directions (, , , and ), the minimum value of the induced normal curvature can be calculated as follows: as shown in Figure 6(b), assume that the unit vectors in the direction of the major and minor axes of the instantaneous contact ellipse are and , respectively. has an unknown angle from the principal direction of the gear . is the angle between and and can be calculated by
Utilizing the Euler equation, the gear and pinion normal curvatures along are given, respectively, as
Hence, the induced normal curvature along is defined as
Taking the extremum of the induced normal curvature , one can obtain two values of . One of them is the angle between and and the other is the angle between and . The normal curvature of the gear ( and ) and the pinion ( and ) in the direction of and can then be obtained by substituting into equations (23) and (24). The radii of curvature in the direction of and is the inverse of the corresponding normal curvature values:
In addition to the geometrical properties of the gear and pinion tooth surfaces at the LCP, the implementation of EHL analysis also requires the motion characteristics. As illustrated in Figure 6(a), the velocity vectors of the gear and pinion at the ACP are given aswhere and are the angular velocity vectors of the gear and pinion, respectively. The velocity vectors and can be resolved along the unit normal of the tooth surfaces and in the tangential plane along the major and minor axes of the elliptical contact footprint as shown in Figure 7. These tangential components, which are used to calculate the entrainment and sliding velocities, can be obtained by the vector dot product:
Therefore, the speed of entraining motion along the major and minor axes of the contact ellipse is defined as
The entrainment and sliding velocities can then be calculated, respectively, by
3. EHL Basic Dimensionless Equations
The Reynolds equation describing the formation of hydrodynamic pressure for elliptical contacts with an angled entrainment flow in the steady-state condition iswith the following nondimensional variables: , , , , , and . Here, and are, respectively, the half-contact width (along the minor axis) and the half-contact length (along the major axis) of the Hertzian contact ellipse and can be calculated according to ref. [24]; is the maximum Hertzian pressure, and for elliptical contacts, where is the load carried by a pair of meshing teeth; and denote the ambient viscosity and density of the lubricant, respectively; is the angle of entrainment velocity relative to the minor axis of the contact ellipse. The dimensionless coefficients for Poiseuille flow and are given bywhere is the contact ellipticity and defined by . The dimensionless boundary conditions of the Reynolds equation are at the inlet and at the outlet. The cavitation condition imposes .
The non-Newtonian behaviour of the lubricant is not taken into account in the present work for simplicity. This simplification does not have a significant impact on the findings of this research because Newtonian models are suitable for the prediction of EHL film thickness, which are the index for evaluating lubrication performance used in this work, even under heavy-loaded conditions. It has been found in many previous studies that the entraining action that governs the formation of the lubricant film occurs in the inlet zone of the EHL contact where the pressure is relatively low, and thus the non-Newtonian behaviour of the lubricant is not apparent.
For the purpose of investigating the impacts of changes in meshing parameters due to the contact point migration on analysis results, smooth surfaces are assumed to exclude the influence of surface topography. Therefore, the film thickness expression for the smooth elliptical contact is applied:where is the normal approach of the gear and pinion tooth surfaces and needs to be adjusted during the iterative process to make the calculation result meet the load balance equation and and represent the undeformed contact geometry where and are, respectively, the equivalent radii along the minor () and major () axes of the contact ellipse. Applying the radii of curvature in the direction of and (calculated by equation (26)), one can obtain and by
It is important to note that “” is used before as the curvature directions of the gear and pinion tooth surfaces along are opposite, whereas “” is used before since the curvature directions of the two surfaces along are the same.
The localized elastic contact deformation is calculated by the Boussinesq integral:where is the reduced Young’s modulus of the two surfaces.
Isothermal condition is assumed in this work, so the viscosity and density of lubricants are only influenced by pressure. Although Doolittle’s free-volume model has been proven to be the most accurate description of the real piezo-viscous behaviour at the typical EHL pressures [25], it is not convenient to use as several Doolittle–Tait parameters must be experimentally measured. In contrast, the Roelands equation is relatively easy to use and can give reasonable calculation results close to those from the free-volume model [26]. Therefore, the Roelands equation is adopted here to describe the pressure-viscosity relation for simplicity:where is the pressure-viscosity index, is a constant, denotes the atmospheric viscosity of lubricants, and is the pressure-viscosity coefficient. The variation of lubricant density with pressure is described by the Dowson–Higginson pressure-density relationship:
The lubricant reaction force defined by the integral of the pressure over the whole solution domain should be balanced by the load acting on each pair of meshing teeth:
In the above equations the unit of maximum pressure is Pa.
4. Numerical Method
The finite difference method is applied in the discretization of the Reynolds equation (Appendix A). In order to ensure the solution convergence and numerical stability under heavy-loaded conditions, the semisystem approach described in ref. [27] is adopted. In accordance with this method, the coefficient matrix of the discretized Reynolds equation is constructed utilizing both the pressure flow and the entraining flow terms, and thereby the diagonal dominance of the coefficient matrix can be guaranteed even when the pressure flow terms almost vanish under extremely heavy load. The discrete Reynolds equation is solved by Gauss–Seidel line relaxation. The full multigrid algorithm [28] is applied to accelerate the solution process and improve numerical accuracy. Surface contact deformation is calculated using the multilevel multi-integration technique [28] (Appendix B). The domain of computation is defined as and . Four layers of the grid are used, with the finest grid consisting of nodes equally spaced. The next coarser grid has nodes, etc. The finest grid has the dimensionless mesh spacing of , which meets the requirement of most engineering applications. The first target grid on which the full multigrid process is started should be sufficiently fine [28], and thus the grid having nodes is selected. The pressure and load convergence criteria are used to determine whether the result converges on the target grid:where and , respectively, denote the new and old approximation of pressure in each grid point and the limit error and are set to be and , respectively.
Figure 8 shows the flowchart of the present numerical procedure. Note that the normal approach will also be adjusted on the coarsest grid according to the multigrid method [28], which is not illustrated in Figure 8.
5. Results and Discussion
A pair of mismatched spiral bevel gears manufactured by the Gleason face milling approach is treated. The pinion is the driving member. The drive sides of the gear pair, namely, the concave side of the pinion tooth and the mating convex side of the gear tooth, are the objects of interest. The basic design parameters are listed in Table 1, with which the machine tool settings used to derive tooth geometry are obtained by the Gleason CAGE4win. The material properties of the gear pair and the lubricant properties are shown in Table 2.
A series of quasistatic isothermal EHL solutions are given at each meshing position over a tooth engagement cycle. The meshing position represented by the pinion roll angle commences at the beginning of the tooth engagement. To investigate the effect of torque loads on contact geometry and lubrication performance, we performed EHL simulations under two typical torque loads of 50 Nm and 300 Nm, respectively, representing lightly loaded and heavily loaded conditions. The rated speed of the driving pinion is 500 r/min. TCA and LTCA were, respectively, performed to obtain the theoretical unloaded contact path and loaded contact paths at 50 Nm and 300 Nm torque load. LTCA also provides variations of contact load (shown in Figure 9), the coordinates of EMPs, and the direction cosine of the effective mesh force vector (listed in Tables 3 and 4) at each meshing position over the tooth engagement. The support of the example gear set is assumed to be rigid in LTCA, so the deformations of bearing and housing are not considered.
Figures 10(a) and 10(b), respectively, compare the contact paths (composed of the EMPs at each meshing position in the projection plane) at 50 Nm and 300 Nm torque loads with the theoretical contact path under no-load condition on the concave side of the pinion tooth and the convex side of the gear tooth. It is shown that the EMPs for both the pinion and the gear move from the heel to the toe of teeth throughout the engagement. Meanwhile, for the pinion, the EMP moves from the tooth root to the tooth top, but for the gear, on the contrary, it shifts from the tooth top to the tooth root. It can also be observed that the increase in load results in the extension of the EMP footprint toward the heel and the toe of teeth. As shown in Figure 10, the contact paths at 50 Nm and 300 Nm torque loads are longer than the theoretical unloaded contact path, and among them, the contact path at 300 Nm is the longest. It indicates that under loaded conditions, the tooth pair comes into meshing earlier and secedes meshing later than the unloaded case due to the deformation.
(a)
(b)
Utilizing the data in Tables 3 and 4, one can calculate the geometrical and kinematical parameters for EHL analysis according to the previously described approaches. Figure 11 displays the major axis direction of the contact ellipse at the LCP of each meshing position calculated by the proposed method. The contact pattern, which is the envelope of the ends of the contact ellipse predicted from LTCA shown in Figure 12, is also drawn in Figure 11 for comparison. It can be seen that the obtained major axis directions of each elliptical contact footprint in the contact pattern are in good agreement with those from LTCA. This indicates that the proposed method is capable of evaluating the contact geometry at the LCP while ensuring that the calculated contact footprints are consistent with the footprints predicted from LTCA.
(a)
(b)
Figures 13 and 14, respectively, show the variations of the radii of curvature and entraining velocities in the direction of the major and minor axes of the contact ellipse over one tooth engagement. They are the input geometric and kinematic meshing parameters of EHL analysis. It can be seen that since the LCPs do not significantly deviate from the TCPs at light load, as shown in Figure 10, the variation trend of meshing parameters over one tooth engagement under light-loaded conditions is similar to their theoretical change trend. In contrast, at heavy loads, a more noticeable difference can be found. Besides, in different meshing positions shown in Figure 8, the deviation direction of the LCP from the TCP is not consistent, so no obvious change rule of meshing parameters with load can be found from Figures 13 and 14.
(a)
(b)
(c)
(d)
(a)
(b)
Both the meshing parameters obtained at the LCPs and the corresponding values calculated at the TCPs are used in the EHL simulation to observe the error involved by ignoring the adjustment of contact points due to deformation. Note that in the areas where the meshing parameters are only presented under loaded conditions as shown in Figures 13 and 14, if the gear pair is considered to be a rigid body, the tooth pair of interest is out of contact, while the adjacent tooth pair is in contact. Therefore, only the calculation results in regions where the tooth pair of interest is in contact under both loaded and unloaded conditions are compared. Figure 15 shows the comparison of the central and minimum film thickness. It demonstrates that, at 50 Nm torque load, the results without considering the effect of deformation on the contact point are similar to those calculated at the LCPs. The maximum differences, which occur at the pinion roll angle of 0.711 rad, are about 6.4% for the central film thickness and 6.5% for the minimum film thickness. However, these differences become significant under the torque load of 300 Nm, reaching to 13.2% for the central film thickness and 15.2% for the minimum film thickness. Figure 16 compares the pressure distributions and film thickness profiles along the major and minor axes of the contact ellipse at the pinion roll angle of 0.711 rad. Though there is a clear difference between the film thickness predicted at the LCP and the TCP, the difference in the corresponding primary pressure peak is relatively small. Recall that although the gear tooth contact geometry at the LCP is different from that at the TCP, the loads carried by the meshing tooth pair at these two positions are the same (751.2 N and 4416 N at the applied torque of 50 Nm and 300 Nm, respectively, shown in Figure 9). The half-contact width and half-contact length of the Hertzian contact ellipse in the four cases are 0.171 mm and 3.220 mm (LCP, 50 Nm), 0.168 mm and 3.014 mm (TCP, 50 Nm), 0.313 mm and 6.166 mm (LCP, 300 Nm), and 0.304 mm and 5.758 mm (TCP, 300 Nm). Hence, the corresponding maximum Hertzian pressures are 0.65 GPa (LCP, 50 Nm), 0.7 GPa (TCP, 50 Nm), 1.1 GPa (LCP, 300 Nm), and 1.2 GPa (TCP, 300 Nm). It can be found that the difference in the maximum Hertzian contact pressure caused by the difference in the tooth surface geometry between the LCP and the TCP is not very large. This is why the primary pressure peak calculated at the LCP does not differ from that obtained at the TCP significantly. However, in addition to the tooth contact geometry, the film thickness also largely depends on the entrainment velocity in the direction of the minor axis of the contact ellipse, which has been found in many previous studies. As can be seen in Figure 14, the entrainment velocity along the minor axis of the contact ellipse calculated at the LCP is higher than the one obtained at the TCP. This may be the dominant reason why the film thickness calculated at the LCP is greater than the oil film thickness calculated at the TCP. Overall, the results indicate that the contact point migration, which may not significantly affect pressure distribution, has nonnegligible effect on the prediction of film thickness. An interesting result is that if the effect of gear deformation on the position of the contact point is taken into account, the lubrication performance of the example gear pair is improved to some extent after loading. For example, at the pinion roll angle of 0.711 rad where the load carried by the meshing tooth pair reaches to the maximum value as shown in Figure 9, the film thickness predicted at the LCP is larger than the value evaluated at the TCP. This indicates that owing to the deformation, the actual contact point moves to a position that is favorable for lubrication.
(a)
(b)
(a)
(b)
(c)
(d)
6. Conclusion
The present study provides a solution of the isothermal elliptical contact EHL for spiral bevel gears, considering the migration of contact points of the meshing tooth pair due to load-induced deformation. The loaded contact points at each meshing position during the tooth engagement are determined by the loaded tooth contact analysis. A method for calculating the geometric and kinematic meshing parameters at the loaded contact points is proposed.
At high loads, the loaded contact points of the deformed gear pair deviate significantly from the theoretical contact points of the corresponding rigid body. If the effect of deformation on the position of the contact point is not taken into account, it may cause a significant error in the EHL simulation results, especially the prediction of film thickness. Therefore, under tremendous loads, it is recommended to use the meshing parameters at the loaded contact points for EHL simulation. Changes in contact points due to gear deformation may improve or worsen lubrication performance. Thus, obtaining the lubrication performance at the loaded contact points under the actual working conditions will help to guide the tooth surface modification.
Appendix
A. Discrete Reynolds Equation
The Reynolds equation (equation (31)) can be transformed into the following discrete differential equation at each point of the grid:
The Poiseuille terms in Reynolds equation are discretized by the second-order central differential scheme aswhere
The coefficients of the Poiseuille terms in equation (A.1) can be obtained as follows:
To maintain numerical stability under heavy loads over a wide range of lubrication conditions, the first-order upstream discretization is selected for the Couette terms. Note that for thick film conditions with high rolling speeds, it is more appropriate to apply the second-order backward differential scheme as it provides higher accuracy due to its smaller truncation errors. The separated form of the Couette flow with respect to is
The discrete form of the density derivative term iswhere and .
The discrete form of the film thickness derivative term iswhere , , and .
The coefficients of the Couette term with respect to in equation (A.1) can thus be given by
Similarly, the separated form of Couette flow with respect to is
The discrete form of the density derivative term is
The discrete form of the film thickness derivative term iswhere , , and .
The coefficients of the Couette term with respect to in equation (A.1) can be obtained as follows:
The total coefficients for equation (A.1) are
In the iterative process, , , , and are all known and calculated by using current nodal pressure values.
B. Discrete Contact Deformation Equation
The discrete form of equation (35) readswhere denotes the pressure acting on the point () and is the influence coefficient which means the deflection at the point () due to the unit force applied at the point (). The coefficient is defined by
Its analytical solution iswhere
Abbreviations
| LCP: | Loaded contact point |
| EHL: | Elastohydrodynamic lubrication |
| EMP: | Effective meshing point |
| LTCA: | Loaded tooth contact analysis |
| TCA: | Tooth contact analysis |
| TCP: | Theoretical contact point. |
| : | Contact semiminor half-width |
| : | Contact semimajor half-width |
| , : | Principal directions of the gear |
| , : | Principal directions of the pinion |
| : | Reduced Young’s modulus |
| : | Vector of effective mesh force |
| : | Dimensionless film thickness |
| : | Normal approach of the gear and pinion tooth surfaces |
| : | Contact ellipticity |
| , : | Principal curvatures of the gear tooth surface |
| , : | Principal curvatures of the pinion tooth surface |
| , : | Normal curvature of the gear tooth surface in the direction of and |
| , : | Normal curvature of the pinion tooth surface in the direction of and |
| : | Direction vector of effective mesh force |
| : | Unit normal vector of the cutting surface |
| : | Unit normal vector of the gear tooth surface |
| : | Unit normal vector of the pinion tooth surface |
| : | Maximum Hertzian pressure |
| : | Actual contact position of the gear |
| : | Effective mesh point |
| : | Actual contact position of the pinion |
| : | Theoretical contact position under no-load condition |
| : | Projected point of the effective mesh point for the gear in plane |
| : | Projected point of the effective mesh point for the pinion in plane |
| : | Position vector of the effective mesh point |
| : | Position vector of the point on the cutting surface |
| : | Position vector of the point on the gear tooth surface in the body-fixed frame |
| : | Position vector of the point on the pinion tooth surface in the body-fixed frame |
| : | Position vector of the point on the gear tooth surface in the global coordinate system |
| : | Position vector of the point on the pinion tooth surface in the global coordinate system |
| , : | Radii of curvature of the gear tooth surface in the direction of and |
| , : | Radii of curvature of the pinion tooth surface in the direction of and |
| , : | Equivalent radii in the direction of and |
| , : | Roughness of gear and pinion surfaces |
| : | Cutting edge parameter |
| : | Unit vector in the direction of the minor axis of the contact ellipse |
| : | Entrainment velocity |
| , : | Speed of entraining motion along and |
| , : | Velocity vectors of the gear and the pinion at an actual contact position |
| , : | Speed of the gear tooth surface at an actual contact position along and |
| , : | Speed of the pinion tooth surface at an actual contact position along and |
| : | Sliding velocity |
| : | Velocity in relative motion between the cutting surface and the tooth surface |
| : | Dimensionless elastic contact deformation |
| : | Load carried by the meshing tooth pair |
| : | Unit vector in the direction of the major axis of the contact ellipse |
| : | Dimensionless load |
| : | Pressure-viscosity index. |
| : | Pressure-viscosity coefficient |
| : | Angle between the gear tooth surface normal vector and the direction vector of effective mesh force |
| : | Angle between the pinion tooth surface normal vector and the direction vector of effective mesh force |
| , : | Dimensionless coefficients for Poiseuille flow |
| : | Limit error of the pressure convergence criterion |
| : | Limit error of the load convergence criterion |
| : | Lubricant viscosity under ambient condition |
| : | Dimensionless lubricant viscosity |
| : | Entrainment angle |
| : | Rotation angle of gear teeth |
| : | Rotation angle of pinion teeth |
| : | Lubricant density under ambient condition |
| : | Dimensionless lubricant density |
| , : | Angular velocity vectors of the gear and the pinion. |
Data Availability
The data used and analyzed to support the findings of this study are available from both the first author and the corresponding author upon request.
Conflicts of Interest
The authors declare no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Acknowledgments
The authors acknowledge Professor Tun Liu of the Harbin Institute of Technology and Professor Jishou Ruan of Nankai University for their valuable discussions during the course of this study. This work was supported by the National Natural Science Foundation of China (no. 51505100); the Science Foundation for Young Scientists of Heilongjiang Province (no. QC2015046); the open project of the National Key Laboratory of Helicopter Transmission Technology, Nanjing University of Aeronautics and Astronautics (no. HTLO19G08).