EHL Analysis of Spiral Bevel Gear Pairs considering the Contact Point Migration due to Deformation under Load

(e 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. (e 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. (en, the osculating paraboloids at this point are extracted from the gear tooth surface geometry. (e 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. (e 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.


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. erefore, 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. is 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. e 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. e 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. en, the geometric and kinematic parameters for EHL simulation are determined according to the differential geometry theory. e 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.

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. e 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. e 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.
erefore, the center of the contact pressure distribution is utilized to determine the LCP in the present work. e 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.
e 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]. e 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. e global fixed reference frame Oxyz 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. O 1 x 1 y 1 z 1 and O 2 x 2 y 2 z 2 are the local fixed reference frames used to describe the rotational positions of the gear and the pinion, respectively. e gear and the pinion rotate about the x-axis of their local fixed coordinates. c is the angle between the gear and the pinion axis. O g x g y g z g and O p x p y p z p are the body-fixed coordinate systems of the gear and the pinion in which their tooth surfaces are represented. e symbol A b c will be often used in this article where the subscript c indicates the coordinate system in which A is represented and the superscript b indicates the object to which A belongs. e procedures of using ABAQUS for LTCA have been introduced in detail by Zhou [14] and Sun et al. [15]. e 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]). e 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 P m is the EMP at a given meshing position and its position vector represented in the fixed coordinate system Oxyz is r m (x m , y m , z m ).  acting on the EMP. P m can first be projected onto the projection plane P 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 (P p P ) and the gear (P g P ) can be calculated by e mathematical model of the tooth flank geometry depends on how the gear teeth are machined but can be generally represented as [16] where r t (s, θ) indicates the cutting surface formed by the rotation of the cutting edge about the axis of the cutter head, s is the parameter describing the cutting edge geometry, θ is the cutter phase angle, and M bt (ϕ) 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. e form of the matrix M bt (ϕ) hinges on the motion of the machine tool. Readers can refer to relevant literatures (e.g., [16,17]) for the definition of the matrix M bt (ϕ). e cutting surface will generate a family of surfaces r b (s, θ, ϕ) 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 r b (s, θ, ϕ). e points on r b (s, θ, ϕ) which simultaneously satisfy the equation of meshing f(s, θ, ϕ) are those on the envelope surface (namely, the tooth surface). In the equation of meshing, n t denotes the unit normal vector of the cutting surface and v tb 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, by f g s g , θ g , ϕ g � n t,g · v tb,g � 0, y g P � x g g s g , θ g , ϕ g , x g P � �������������������������� y g g s g , θ g , ϕ g 2 + z g g s g , θ g , ϕ g 2 ,  where x i i , y i i , and z i i (i � p, g) are the three components of the position vector of the LCP r i i (x i i , y i i , z i i ) represented in the body-fixed frames of the pinion and the gear (O i x i y i z i ).
e 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: As mentioned previously, the LCP is different from the theoretical unloaded contact point calculated based on the principle of local conjugation [16]. at 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 (P p ) and gear (P g ) tooth surfaces do not coincide and their normals (n p and n g ) 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. e 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. e rotation of the gear and pinion teeth is one of the manifestations of the total deformation [22]. ey are, respectively, simulated by the rotation angles of the gear (Δθ g ) and the pinion (Δθ p ) around their rotational shafts. According to the calculation process of the LCP described before, it is obvious that P p and P g will coincide after the rotation as illustrated in Figure 5(b). e gear and pinion tooth surfaces are now represented in the coordinate system Oxyz as follows: and at the LCP, there is Usually, c � π/2, and thus the transformation matrices in equations (5) and (6) are defined as en, the rotation angle of the gear tooth (Δθ g ) and the pinion tooth (Δθ p ) 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 (n g and n p ) 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 Δc g and Δc p , respectively ( Figure 5(c)). Taking the pinion as an example, the procedure of eliminating the angle Δc p is described. A local coordinate system where the three components of the direction vector L m (n mx , n my , n mz ) of the effective mesh force F m (F mx , F my , F mz ) 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 Oxyz to where x, y, and z are the unit vectors along the axes of the coordinate system Oxyz and x p , y p , and z p are the three components of the position vector of P p represented in Oxyz. en, the pinion tooth surface along with the co e 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. e tooth surfaces of the gear and the pinion should be represented in the same coordinate system, e.g., Oxyz, 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 r i L′ � r i L (i � p, g). erefore, the tooth surfaces of the gear and the pinion in Oxyz 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 Oxyz are r p (θ p , ϕ p ) and r g (θ g , ϕ g ), and the corresponding normal vectors of tooth surfaces at the LCP are n p (θ p , ϕ p ) and n g (θ g , ϕ g ). e 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.
e smaller the value of the induced normal curvature, the closer the two surfaces. erefore, 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.
To obtain the minimum value of the induced normal curvature, the principal curvatures and principal directions at the LCP can be calculated firstly. e normal curvature in the arbitrary direction is given by g x z y p y 1 (y g )  Mathematical Problems in Engineering where Φ 1 and Φ 2 are the first and second fundamental forms of the surfaces, the coefficients (E, F, and G) of the first fundamental form are given as and the coefficients (L, M, and N) of the second fundamental form are defined as where the unit normal to the surface is calculated by Equation (15) can be rewritten as follows: where μ � dθ/dϕ 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 dk n /dμ � 0, one can obtain the relationship between the principal curvatures and principal directions as e equation of principal curvatures can be derived by eliminating dθ/dϕ from equation (20): e calculated principal curvatures can be substituted into equation (19) to evaluate the corresponding principal directions. After getting the principal curvatures of both the gear (k g 1 and k g 2 ) and pinion (k p 1 and k p 2 ) tooth surfaces, as well as the corresponding principal directions (e g 1 , e g 2 , e p 1 , and e p 2 ), 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 w and u, respectively. w has an unknown angle α from the principal direction of the gear e g 1 . β is the angle between e g 1 and e p 1 and can be calculated by Utilizing the Euler equation, the gear and pinion normal curvatures along w are given, respectively, as Hence, the induced normal curvature along w is defined as k r w � k g w − k p w � k g 1 cos 2 α + k g 2 sin 2 α − k p 1 cos 2 (α + β) − k p 2 sin 2 (α + β).

(25)
Taking the extremum of the induced normal curvature dk r w /dα � 0, one can obtain two values of α. One of them is the angle between e g 1 and w and the other is the angle between e g 1 and u. e normal curvature of the gear (k g w and k g u ) and the pinion (k p w and k p u ) in the direction of w and u ω g ω p y 1 can then be obtained by substituting α into equations (23) and (24). e radii of curvature in the direction of w and u 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 as where ω g and ω p are the angular velocity vectors of the gear and pinion, respectively. e velocity vectors v g and v p can be resolved along the unit normal of the tooth surfaces and in the tangential plane T along the major and minor axes of the elliptical contact footprint as shown in Figure 7. ese tangential components, which are used to calculate the entrainment and sliding velocities, can be obtained by the vector dot product: erefore, the speed of entraining motion along the major and minor axes of the contact ellipse is defined as e entrainment and sliding velocities can then be calculated, respectively, by

EHL Basic Dimensionless Equations
e Reynolds equation describing the formation of hydrodynamic pressure for elliptical contacts with an angled entrainment flow in the steady-state condition is with the following nondimensional variables: Here, a and b 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]; p H is the maximum Hertzian pressure, and for elliptical contacts, p H � 3w/(2πab) where w is the load carried by a pair of meshing teeth; η 0 and ρ 0 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. e dimensionless coefficients for Poiseuille flow ε x and ε y are given by where k is the contact ellipticity and defined by k � b/a. e dimensionless boundary conditions of the Reynolds equation are P(X 0 , Y) � P(X, Y 0 ) � 0 at the inlet and P(X e , Y) � zP(X e , Y)/zX � P(X, Y e ) � zP(X, Y e )/zY � 0 at the outlet. e cavitation condition imposes P(X, Y) ≥ 0. e non-Newtonian behaviour of the lubricant is not taken into account in the present work for simplicity. is 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 Figure 7: Velocity relation.
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. erefore, the film thickness expression for the smooth elliptical contact is applied: where H c0 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 aX 2 /2R e u and kbY 2 /2R e w represent the undeformed contact geometry where R e u and R e w are, respectively, the equivalent radii along the minor (u) and major (w) axes of the contact ellipse. Applying the radii of curvature in the direction of u and w (calculated by equation (26)), one can obtain R e u and R e w by 1 It is important to note that "+" is used before 1/R p u as the curvature directions of the gear and pinion tooth surfaces along u are opposite, whereas "− " is used before 1/R p w since the curvature directions of the two surfaces along w are the same. e localized elastic contact deformation V(X, Y) is calculated by the Boussinesq integral: where E' 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 piezoviscous 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]. erefore, the Roelands equation is adopted here to describe the pressure-viscosity relation for simplicity: where Z is the pressure-viscosity index, c p � 1.96 MPa is a constant, η 0 denotes the atmospheric viscosity of lubricants, and α is the pressure-viscosity coefficient. e variation of lubricant density with pressure is described by the Dowson-Higginson pressure-density relationship: ρ � 1 + 0.6 × 10 − 9 p H P 1 + 1.7 × 10 − 9 p H P .
(37) e 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 p H is Pa.

Numerical Method
e 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. e discrete Reynolds equation is solved by Gauss-Seidel line relaxation. e 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). e domain of computation is defined as − 2.5 ≤ X ≤ 1.5 and − 2 ≤ Y ≤ 2. Four layers of the grid are used, with the finest grid consisting of 257 × 257 nodes equally spaced. e next coarser grid has (128 + 1) × (128 + 1) nodes, etc. e finest grid has the dimensionless mesh spacing of ΔX � ΔY � 0.015625, which meets the requirement of most engineering applications. e first target grid on which the full multigrid process is started should be sufficiently fine [28], and thus the grid having (64 + 1) × (64 + 1) nodes is selected. e pressure and load convergence criteria are used to determine whether the result converges on the target grid: where P n i,j and P n− 1 i,j , respectively, denote the new and old approximation of pressure in each grid point and the limit error ε p and ε w are set to be 10 − 5 and 10 − 3 , respectively. Figure 8 shows the flowchart of the present numerical procedure. Note that the normal approach H c0 will also be adjusted on the coarsest grid according to the multigrid method [28], which is not illustrated in Figure 8.

Results and Discussion
A pair of mismatched spiral bevel gears manufactured by the Gleason face milling approach is treated. e pinion is the driving member. e 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. e 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. e 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. e 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. e 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. e 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.
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. e 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. is 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. 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. ey 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.
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. erefore, 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. e 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. ough 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). e half-contact width a and half-contact length b of the Hertzian contact  Hertzian contact pressure caused by the difference in the tooth surface geometry between the LCP and the TCP is not very large. is 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. is 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    Ends of the contact ellipse  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. is indicates that owing to the deformation, the actual contact point moves to a position that is favorable for lubrication.

Conclusion
e 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. e 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.
erefore, 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. us, obtaining the lubrication performance at the loaded contact points under the actual working conditions will help to guide the tooth surface modification.