Dynamics of lubricated spiral bevel gears under different contact paths

To assess the meshing quality of spiral bevel gears, the static meshing characteristics are usually checked under different contact paths to simulate the deviation in the footprint from the design point to the heel or toe of the gear flank caused by the assembly error of two gear axes. However, the effect of the contact path on gear dynamics under lubricated conditions has not been reported. In addition, most studies regarding spiral bevel gears disregard the lubricated condition because of the complicated solutions of mixed elastohydrodynamic lubrication (EHL). Hence, an analytical friction model with a highly efficient solution, whose friction coefficient and film thickness predictions agree well with the results from a well-validated mixed EHL model for spiral bevel gears, is established in the present study to facilitate the study of the dynamics of lubricated spiral bevel gears. The obtained results reveal the significant effect of the contact path on the dynamic response and meshing efficiency of gear systems. Finally, a comparison of the numerical transmission efficiency under different contact paths with experimental measurements indicates good agreement.


Introduction
Dynamics, which interrelates noise, durability, and vibration problems, is believed to be an important indicator in gear design owing to the mutual effect of dynamics, tribology, and fatigue problems. Mesh forces may increase significantly under dynamic conditions, and they are transmitted through the shaft and bearing into the gear housing, resulting in excessive structure vibration. Moreover, the fatigue life of the two interaction surfaces is significantly affected by the fluctuating load generated by vibration. Owing to mounting errors or deformations of the bearing supporting system, the tooth surface contact area will differ from the designed contact path during actual operations. Hence, the contact path is typically moved to the heel and toe of the gear flank to verify the static contact quality [1]. However, unlike spur gears, the contact geometry, kinematics, and mesh stiffness, believed to be important excitations for gear dynamics [2], are sensitive to the contact paths owing to the complicated spatial surface of gear flanks in spiral gears. Consequently, investigations into the effect of contact path on the dynamics and meshing efficiency of spiral bevel gears can provide a full assessment of their transmission quality.
The dynamics of gears has been extensively investigated previously, particularly for parallel axis transmission, which focuses on various effect factors, such as time-variant parameters [2,3], lubrication [4,5],  [6,7], tooth profiles [8], and assembling errors [9]. Although numerous studies regarding gear dynamics have been published, studies regarding the dynamics of spiral bevel gears are limited owing to the complicated meshing geometry and kinematics. Donley et al. [10] proposed a dynamic hypoid gear model, in which the line-of-action and mesh position were assumed to be invariant. Furthermore, nonlinear dynamic behaviors of spiral bevel and hypoid gears have been simulated [11,12], where time-variant parameters were involved. Based on the proposed dynamic model, the effects of the drive and coast sides (asymmetry of mesh stiffness nonlinearity) on spiral bevel and hypoid gear dynamics were investigated [13]. In Refs. [11][12][13], a torsional dynamic model (two-DOF) was reduced to a one-DOF model that disregarded the bearing support and gear flank friction. Furthermore, multi-DOF models of bevel and hypoid gear systems have been proposed [14,15], and the dynamic responses to the bearing stiffness and torque load were investigated. To obtain more detailed dynamic characteristics for each meshing pair, a multipoint hypoid gear mesh model based on tooth contact analysis (TCA) was established in Ref. [16]. The aforementioned dynamic models were assumed to be dry instead of the lubricated condition of the meshing tooth pair. The dynamics of lubricated spiral bevel gears were analyzed [17] based on a torsional dynamic model, and the results were compared with those from a one-DOF model developed by Ref. [11]. Mohammadpour et al. [18][19][20][21] proposed a multiphysics tribo-dynamic model considering mixed lubrication and bearing supports to investigate the transmission efficiency and other dynamic behaviors. Yavuz et al. [22] investigated the dynamic mesh force in the frequency domain under different backlash and bearing stiffnesses. The shafts and their flexibilities were numerically simulated using Timoshenko beam finite elements, but the mesh line-of-action and position were equivalently treated as invariant. Alves et al. [23] proposed a static and dynamic model for spiral bevel gears to investigate the tooth flank contact pressure under dynamic and static conditions. Friction was omitted in the abovementioned studies [22,23]. As mentioned above, most studies focused on the effect of nonlinear time-varying mesh parameters, backlash nonlinearity, load, etc. on dynamic responses, whereas lubricated conditions were disregarded. Only a few reports regarding the effect of assembly errors on elastohydrodynamic lubrication [24] and the effect of contact path on contact fatigue [25] under static conditions in spiral bevels have been published. The conclusions indicated that the contact path affects the lubrication characteristics and fatigue life significantly. However, the effects of the contact path on the dynamics and efficiency of a lubricated spiral bevel gear have not been reported. Therefore, the investigation into the effect of the contact path on dynamics will benefit future studies pertaining to lubrication and fatigue life under nonlinear dynamic conditions. Hence, an eight-DOF dynamic model was developed in the present study based on a TCA model and an analytical friction model to simulate the nonlinear dynamics and meshing efficiency of spiral bevel gears under different meshing paths. The analytical friction model was demonstrated to be reasonable by comparing the present friction model with a previously published mixed elastohydrodynamic lubrication (EHL) model of spiral bevel gears. Finally, the meshing efficiency was calculated and compared with the numerical results.

Assembling parameters for different contact paths
The aim of this study is to reveal the effect of contact path on dynamic responses; a schematic illustration of the contact path is shown in Fig. 1. Unlike involute spur gears, the contact path and surface parameters of spiral bevel gears are difficult to obtain analytically. Therefore, before modeling the dynamics of spiral bevel gears, a TCA model is required to determine the contact path and relevant contact parameters, such as the principal directions, principal curvatures, contact radii, entraining and sliding vectors, contact load, and static transmission error at transient meshing positions. The TCA model was programmed as a computer package using Formula Translation (FORTRAN), and the methodology has been described in Refs. [26,27]. The derivations of tooth contact parameters are laborious; therefore, this study focuses on the effect of contact path on the dynamics and meshing efficiency of spiral bevel gears. Hence, the determination of the contact path is provided briefly below for clarity. To obtain the different contact paths, the gear and pinion were first assembled at the designed point ( Fig. 1) using the assembling parameters [27,28], which included the pinion axial, vertical offset, and gear axial adjustment, denoted as ΔH, ΔV, and ΔJ, respectively. The initial point was determined by the axial and radial projections gr L and gr R , respectively. Subsequently, the mesh parameters for the different contact paths were computed using the TCA model. Figure 2 shows the contact relationship between the pinion and gear, in which O and  O are the intersection points between the pinion axis p p and gear axis g p (unit vector) before and after the adjustment, respectively, whereas points p O and g O denote the predesigned crossing points of the two axes. As shown in Fig. 2 where bi R , i n , and i t (i = p, g) are the position vector, unit normal vector, and surface unit tangential vector at a transient meshing position, respectively.   is the two-axis angle (shaft angle) between p p and g p , and denotes the transformation matrix from It is assumed that Eq. (2) is satisfied when normal vectors (g) p n and g n rotate about (g) p p and g p with angle  p and  g [28], respectively.
In addition, the vectors in system g O S will be updated owing to the rotation angles of the pinion and gear, which are expressed as

M p
is the rotational transform matrix of the pinion with respect to vector (g) p p with angle

M p
represents the rotational transform matrix of the gear. Furthermore, the conjugated points must satisfy the conjugation theory of a space curved surface as follows [27]: where s V is the relative sliding velocity of two conjugated surfaces.
When the initial running position (designed point) is determined, the mating gear and pinion are assembled in the target position through adjustments H , J , and V , as depicted in Fig. 1. The adjustments can be computed as follows [25,27]: If H , J , and V are calculated, the pinion and gear can be assembled at the expected contact point based on the corresponding adjustment values. Generally, H and V are sufficient for mating the pinion and gear on the designed point, i.e., J can be set as zero.
After the pinion and gear are assembled, the contact parameters can be obtained using the TCA model [25] under different contact paths in a mesh cycle. In fact, the contact parameters are dependent on the machining settings during the machining process, particularly the relative kinematics between the cutter and gear blank [26]. Relevant descriptions of the contact geometries and surface parameters are available in a previous study [25].

Dynamic model
The geared system adopted in the present study comprised a spiral bevel gear pair and tapered roller bearings, as illustrated in Fig. 3. If the flexibility of the shaft is considered, then a finite element method (FEM) can generally be used to model the gear shafts [22]. It is well known that the FEM is time consuming. In fact, the bending effect of a shaft on the system dynamics is limited, as indicated experimentally (Fujii et al. [29]) and theoretically (Gosselin [30]) for a similar dynamic system. Hence, the deformation of the shaft was not considered in the present study. A threedimensional (3D) dynamic model under different contact paths in the spiral bevel gears is illustrated in Fig. 4. The transmission model of the pinion and gear was discretized in terms of the time-varying mesh stiffness m ( ) k t , mesh damping m ( ) c t , gear backlash 2b, and kinematic transmission error m ( ) e t along the line-of-action direction. As shown in Fig. 4, the translational displacements, which can be defined as   ( , , , ) , were considered; furthermore, the subscript  , i p g refers to the pinion and gear, respectively. It is noteworthy that the dynamic model is described in the global coordinate system ( , , , ) where p R and g R are the contact radii. Owing to the change in the contact path, the contact radii are variant and can be computed as follows: It is noteworthy that  T p p p ( , , ) x y y n and  g g ( , , x y n T g ) y denote nonlinear displacements along the line-of-action due to the lateral and axial motions of the pinion and gear axis, respectively. Using the backlash nonlinear, the dynamic mesh force m F can be expressed as In Eq. (10), m ( ) k t is the mesh stiffness that can be calculated using the loaded tooth contact analysis (LTCA) model. LTCA is typically developed based on a finite element (FE) model or FE-based models [31,32]. However, the FE model is extremely time consuming [33]. In this study, an efficient LTCA model proposed by Sheveleva et al. [34] was adopted, and detailed explanations of this model are available in Ref. [34].
Displacements i x , i y , and i z (  , i p g ) are axial and lateral motions that correspond to the deflections of the supporting bearings. The tapered roller bearing is shown in Fig. 5. The method for calculating the load and stiffness calculation is mature [35]. For conciseness, only a brief introduction of the bearing load is presented herein. The bearing loads caused by the axial and radial displacements are expressed in the integral form as follows [35]: where n is a constant, i.e., n = 10/9 for a line contact; Z is the number of tapered rollers; n k is the nonlinear stiffness due to the assembly of the inner ring, outer ring, and roller elements, and it is related to the material properties and bearing geometry; represents the maximum bearing deflection in the direction of the resultant force vector;  l is the half-loaded area angle;  1 denotes the bearing contact angle. When the bearing load is attained, the bearing supporting stiffness is calculated.  Methods to calculate the gear mesh and bearing forces have been developed; therefore, the differential equation governing the dynamics of the spiral bevel gear system is expressed as where    p g p p p g g g ( , , , , , , , ) where p I and g I denote the rotational inertia of the pinion and gear about its axis, respectively; p m and g m are the masses of the pinion and gear, respectively. The stiffness matrix K includes the mesh stiffness and bearing stiffness. The damping matrix C is expressed as   2 m C K, where  is the damping ratio, which can be obtained from Refs. [17,18]. F is the force vector that includes external excitations and internal forces. The external excitation force is the torque fluctuation, and the internal excitation force is a result of the time-varying spatial vector, transmission error, backlash, and friction torque [17].
Matrices K and C , and vector F will not be expanded comprehensively herein for brevity, as they have been derived previously [36]. It is noteworthy that was assumed in Refs. [18,21]; subsequently, the dynamic model was reduced as a seven-DOF system. However, the rate of change of gear teeth contact radii may result in more complicated dynamic responses, such as severe tooth separations, particularly at higher speeds [17]. Hence, the rate of change of the contact radii was considered in the present study. To improve the computational efficiency when using Eq. (13), the normalization was performed in this study as follows: where  n is the reference frequency, which is often selected as the resonant frequency. Based on Eq. (16), the equation of motion is rewritten as In Eq. (18), i C , i K , and i F (i = 1, 8) are the corresponding elements in matrices C , K , and F , respectively. Subsequently, Eq. (17) can be solved using the Runge-Kutta method.

Gear friction model
The excitation in the torsional direction comprises the applied torques p T and g T as well as the friction torques pf T and gf T of the pinion and gear owing to gear flank friction, respectively. When the film in the conjugated gear flank is thin, mixed lubrication occurs, and the mesh load is supported by asperity contact and a film simultaneously. The authors have previously investigated the friction characteristics of spiral bevel gears under different contact paths [25] using a mixed EHL model that can accommodate 3D surface roughness. However, the computations of the governing equation of the mixed EHL model are time consuming. To reduce the solving burden, the friction coefficient was predicted using an analytical method, and it will be compared to the results from the mixed EHL model [25] in later discussions to demonstrate the feasibility of the proposed analytical friction model.
A mixed lubrication condition was considered. The friction force f action on the gear flank comprised viscous shear friction v f and boundary friction b f , expressed as follows: To calculate the boundary friction b f , a Gaussian asperity contact model [18,36] was used in the present study. The boundary friction force can be calculated using the boundary friction coefficient [25,36] (20) where  denotes the coefficient of dry or boundary contact, generally assumed to be constant [25,36]. In this case,  was set to 0.13. According to Ref. [37], the load shared by asperities a W and the asperity contact area a A can be expressed as As suggested by Greenwood and Tripp [37], the roughness parameter    where    c h is the film thickness ratio,  is the composite root mean square roughness, and c h is the film thickness. The c h was calculated using an analytical film thickness formula for elliptical point contacts considering the oblique entraining angle [38,39], which was originally obtained under light load conditions [38]. However, Wang et al. [40] and Jalali-Vahid et al. [41] discovered that the curve-fitting formula by Chittenden et al. [38] can yield reasonable predictions of the film thickness compared with numerical results under a heavy-load operating environment with arbitrary entrainment. The curvefitting formula is expressed as follows: For viscous stress  , a viscoelastic non-Newtonian fluid model (Bair and Winer [42]) can be used as follows: where the lubricant viscosity  is assumed to be a function of pressure, and a typical relationship is    p e [25], which has been justified to be suitable experimentally by He et al. [43] for computing the shear force in a wide range of loads. The limiting shear elastic modulus  G and the limiting shear stress The viscous shear stress in the contact zone is related to the contact pressure. In the present study, contact pressure was discretized using a Hertzian contact model [39], which has been demonstrated as a reasonable assumption for spiral bevel gears [45].
Once the central film thickness and sliding velocity www.Springer.com/journal/40544 | Friction vector are provided, the shear rate of the lubricant at the center of the mesh can be computed. The shear rate can be expressed as a linear relationship, as widely used in Refs. [39,44], which can be expressed as Solving Eq. (30), the average viscous shear stress,  , can be obtained by averaging the local shear in the elliptical contact zone. Subsequently, the viscous friction is obtained as follows: Before calculating the frictional torque, the moment arms pf R and gf R applied to the pinion and gear must be computed. The sign of friction is determined by the direction of the sliding velocity. The friction torque may assist or resist the motion of the pinion and gear.
Hence, it is necessary to compute the moment arms pf R and gf R while considering the sign of the relative sliding velocity, as follows: Subsequently, the total frictional torques pf T and gf T are expressed as where   1, , k N is the k-th meshing gear pair that is determined using the TCA model. The friction coefficient  ( ) k for each conjugated gear pair k is computed using Eq. (19).
Based on the friction model, the instantaneous efficiency of the spiral bevel gear can be estimated as It is noteworthy that the rolling friction loss is considered, and the rolling friction force ro F is calculated as [46,47]   The thermal reduction factor T C is defined as [45,46]

Numerical result analysis
The parameters of the spiral bevel gears and assembled bearings are listed in Table 1. Additionally, three different contact paths are depicted in Fig. 1. The width of the gear flank is w B , and design points 1, 2, and 3 are located at the pitch cone; their coordinates ( gr L , gr R ) are (40.01 mm, 117.43 mm), (36.54 mm, 107.25 mm), and (33.08 mm, 97.08 mm), respectively. The contact paths through points 1, 2, and 3 are referred to as the heel, middle, and toe contacts. The input torque acting on the pinion was set as 200 N·m. The flowchart of the methodology of the dynamics of a spiral bevel gear under different contact paths is summarized in Fig. 6. As shown in Fig. 6, the TCA analysis involves complicated numerical processes for attaining the assembling and meshing parameters under different contact paths.
The three types of tooth contact trajectories are plotted in Fig. 7, and the corresponding assembling adjustments, obtained using the methods described in Section 2.1, are listed in Table 2. Under different contact paths, the relevant parameters for the dynamic model were calculated using the TCA model.   static transmission error (kinematic error) from the meshing-in to the meshing-out point. It is clear that the mesh stiffness m ( ) k t was relatively large for the heel contact, and the stiffness was affected by the contact ratio. The static transmission error m ( ) e t depended on the microgeometry and manufacturing, and it appeared as a sinusoidal-like form, as shown in Fig. 8. The transmission error was significant at the toe contact. Figure 9 summarizes the pinion and gear contact radii, p R and g R . The results show that the variation in the contact radii was limited. Therefore, the assumptions of constant contact radii and invariant rate of change of the contact radii can be reasonable at low speeds. Figure 10 shows the curvature radii along the minor and major axes of the contact ellipse, which are related to friction calculations. The frictional moment arms of the pinion and gear are shown in Fig. 11, and it is clear that the sign of the arms changed at design points 1, 2, and 3. To incorporate these time-variant parameters into a dynamic model, Fourier series functions with respect to the pinion rotational angle were applied in the present study to simulate the periodical parameters [17] during the meshing of spiral bevel gears.
The gear materials, lubricant, and roughness parameters for the present simulations were based on those in Ref. [25]. Figure 12 presents the maximum and minimum amplitudes of the DTE during different speeds for the heel, middle, and toe contacts. During the speed sweep, the critical resonance regions occurred at approximately 10,400 rpm for the toe and middle contacts and 11,000 rpm for the heel contact. In the resonance region, the amplitudes of the DTE of the middle and heel contacts fluctuated in a range larger      | https://mc03.manuscriptcentral.com/friction than that of the toe contact. Except for the resonance, the toe contact exhibited a large DTE. A clear jump phenomenon was observed, as was discovered in Refs. [11,18], particularly for middle and heel contacts. The time histories of the dynamic transmission error for the toe, middle, and heel contacts under the critical resonance speed are plotted in Fig. 13, depicting that the contact paths primarily affected the values of the minimum DTE instead of the maximum DTE at the resonance regions. The dynamic mesh force amplitudes at different speeds are illustrated in Fig. 14. The responses of the dynamic mesh force with respect to the pinion speed were similar to the dynamic transmission error. In the vicinity of resonance, the minimum force became zero, indicating the occurrence of teeth separation, resulting in contact loss. In addition, in the frequency region, the heel contact occupied a wider speed range, where separation occurred, compared with the case of middle and toe contacts. The periods of responses of the dynamic mesh force and its corresponding maximum Hertzian contact pressure are summarized in Figs. 15 and 16. As shown in Fig. 15, the dynamic mesh force of the heel contact was the greatest,     whereas the force was the minimum for the toe contact. However, as shown in Fig. 16, the maximum Hertzian contact pressure h p was high for the toe contact compared with those of the heel and middle contacts, although the meshing force was relatively low for the toe contact. This was because the surface geometries were different under different contact paths, as indicated in Fig. 10 Fig. 17. The maximum octahedral stresses were 1.76, 2.27, and 2.44 GPa under the heel, middle, and toe contacts, respectively. Despite the relatively small contact force for the toe contact, as shown in Fig. 15, conspicuous surface stress concentrations were observed owing to intermittent asperity contacts, which directly caused premature surface micropitting [48,49]. The stress solution was obtained from a mixed EHL model and an octahedral stress equation, which have been described in our previous study [25]. For brevity, the formulae of the mixed EHL model and stress are omitted herein, and readers can refer to Ref. [25] for details. Additionally, the higher Hertzian contact pressure generated larger stress distributions and stress-affected volumes, which dominated the contact fatigue life [25]. The radial and axial displacements of the pinion and gear under different contact paths during a speed sweep are shown in Figs. 18 and 19, respectively. For the pinion, the radial displacement was the resultant displacement of p x and p y , and the axial displacement was p z . For the gear, g y and g z represent the radial displacement, and g x represents the axial displacement. The radial and axial displacements of the pinion exhibited a trend similar to that of the dynamic transmission error. In a wide speed range, the amplitude of the radial displacement response of the pinion was greater than that of the gear. However, for the toe contact of the gear, a significant discontinuity in radial displacement was discovered at 8,800 r/min, and the amplitude was approximately 100 μm, which was much larger than the radial displacement of the pinion. In addition, the tendency of the gear axial displacement with respect to speed differed from that of the pinion, as shown in Fig. 19. Compared with the middle and heel contacts, the axial displacement amplitude of the toe contact fluctuated in a wide range, and the maximum displacement was large. Analyses of Figs. 12, 14, 18, and 19 show that the responses of the mesh force and DTE were similar to those of the axial and radial displacements of the pinion. It can be concluded that the dynamic mesh force and dynamic transmission error under different contact paths were primarily affected by the pinion displacements. In addition, the vibration of the gear was severe under the toe contact path.
The lateral and axial displacements of the shaft resulted in structural excitations that transmitted to the differential housing through bearings. A case study of bearings A nd C was performed, and the variation in the transmitted force through the supporting bearings in the axial and lateral directions are depicted in Figs. 20 and 21, respectively. For bearing A, the results were generally similar to the trends of the DTE and dynamic mesh force variation. For bearing C, the axial and radial bearing forces under toe contact were extremely high at approximately 8,800 r/min, consistent with the variation in the gear lateral displacement, as depicted in Fig. 18. Furthermore, it was discovered that the bearing force under the toe contact was greater than those under the middle and heel contacts apart from the resonance regions. Additionally, it was observed that the axial bearing force was much lower than the lateral bearing force, particularly in the resonance region.
The meshing efficiency of spiral bevel gears is related to the friction power loss; therefore, an accurate friction model is required for predicting the instantaneous   meshing efficiency. Only a few studies have focused on friction in spiral bevel or hypoid gears, such as those from Xu and Kahraman [46], Kolivand at al. [47], and Paouris et al. [39]. An analytical method of the friction model was used in Refs. [18,19,39]; however, it has not been validated for the application of spiral bevel or hypoid gears. Xu et al. [46,47] investigated the efficiency of hypoid gears, whereas the contact was assumed to be a line contact. Xu and Kahraman [46] proposed a fitting formula for the friction coefficient based on a significant amount of mixed EHL (line-contact model) analyses; it was expressed as a function of the maximum Hertzian contact pressure h p , slid-to-roll ratio SRR, entraining velocity  . To indicate the effect of the line-contact assumption on friction predictions, the results obtained using the method from Xu and Kahraman [46] were compared to those obtained from the mixed EHL model of spiral bevel gears [25]. The reliability of the mixed EHL model applied in spiral bevel gears was validated in Ref. [50]. In addition, the predictions of the present analytical friction model were compared with the results from the mixed EHL model. The friction coefficient predictions from different friction models under different contact paths are plotted in Fig. 22. It is noteworthy that the applied rotational www.Springer.com/journal/40544 | Friction speed and torque of the pinion were 3,000 r/min and 190 N·m, respectively. It was observed that the friction coefficient from the mixed EHL model [25] first increased and subsequently decreased, reaching the maximum at the pitch cone. Similar results have been reported in Refs. [51,52], where a relatively realistic lubrication model (the entrainment angle was considered) of a spiral bevel gear was employed. The friction coefficient of the toe contact was relatively high compared with those of the middle and heel contacts. As shown in Fig. 22, the friction model with a line-contact assumption proposed by Xu and Kahraman [46] indicated a relatively large prediction error around the pitch cone owing to the negligence of the entrainment angle. This indicates that the simplification of the line contact was reasonable for the friction analysis of spiral bevel gears apart from the neighboring pitch cone. The friction coefficient of the present analytical model was consistent with the results of the mixed EHL model for the toe, middle, and heel contacts. To further demonstrate the analytical model, the center film thickness was analyzed, as shown in Fig. 23. It was clear that the film thickness from the analytical model agreed well with the mixed EHL predictions. The static meshing efficiency achieved by the proposed model for the toe, middle, and heel contacts is plotted in Fig. 24. The maximum efficiency was reached in the vicinity of the pitch cone where the sliding velocity was the minimum [25].
Once the friction model was developed, the instantaneous meshing efficiency can be analyzed using the tribo-dynamic model. Figure 25 shows the averaged meshing efficiency over a wide speed    | https://mc03.manuscriptcentral.com/friction range. It was observed that the efficiency increased with the pinion speed when the pinion speed was less than 6,000 rpm. In the resonance regions, the efficiency fluctuated significantly owing to the tooth separations, thereby resulting in the disappearance of friction loss. Furthermore, it was evident that the efficiency of the toe contact was higher than those of the middle and heel contacts. Figure 26 shows the history of the meshing efficiency and the dynamic friction coefficient in a mesh cycle for the case where the rotational speed and torque of the pinion were 3,000 r/min and 190 N·m, respectively. Compared with Fig. 24, the dynamic meshing efficiency was lower than the static efficiency, as expected, owing to the power loss in vibration of the shaft in the spiral bevel gears along the lateral and axial directions. Although the difference in the friction coefficient was limited for different contact paths, the minimum instantaneous efficiencies were 89.1%, 89.5%, and 91.6% for the heel, middle, and toe contacts, respectively. This was because the sliding velocity was relatively high for the heel contact [25] owing to the large rotational radii, as illustrated in Fig. 9.

Experimental results
The friction, which is related to the transmission efficiency, was introduced to the dynamic model under different contact paths. Hence, the transmission efficiency was tested to verify the methodology used in the present study. Transmission efficiency tests were performed using a gear transmission system test rig, as shown in Fig. 27, to validate the dynamic model coupled with friction. The parameters of the tested gear pair are shown in Table 3, and the parameters of the assembled bearings in the test rig were the same as those listed in Table 1. The assembly adjustments for the toe, middle, and heel contacts, obtained using the methods described in Section 2.1, are listed in Table 4. In the experiment, Mobil gear oil 600XP150 was used as the lubricant. The parameters of the gear materials, lubricant, and root mean square (RMS)     , where the torque and angular speed were measured based on the torque sensor and angular encoder at a sampling frequency of 1,000 Hz.
The maximum mechanical speed of the output angular encoder (mounted on the driven side) and input angular encoder (mounted on the driving side) were 1,000 and 3,000 r/min, respectively. The maximum input and output torques of the motor were 96 and 236 N·m, respectively. It is noteworthy that the shaft speeds were measured using an angular encoder integrated in a motor with a wide speed range of 0-6,000 r/min, and they were not affected by the protective speed of the output angular encoder (1,000 r/min). In a smaller torque range, the effect of torque on efficiency was limited compared with that of speed. Hence, the efficiency was tested in a pinion speed range of 10-1,500 r/min with a load of 60 N·m acting on the gear, and the results are summarized in Fig. 28. As shown in Fig. 28(a), the measured efficiency increased with the speed, and the efficiency from large to small was that of the toe, middle, and heel contacts, coinciding with the trend of the numerical results. In addition, the numerical predictions agreed well with the tests at different speeds and contact paths; however, the former appeared slightly larger than the latter. This deviation may be a result of the subtraction error of the internal friction caused by the motor, bearing, and shafting, particularly at 10 r/min. The deviations between the experimental and numerical results were significant because of the effect of internal friction loss.

Conclusions
The static meshing quality in spiral bevel gears is generally verified under different contact paths; however, the dynamic characteristics under different contact paths have not been reported. Hence, the effects of contact paths on the dynamic response and meshing efficiency of a lubricated spiral bevel gear pair were analyzed based on the combination of an eight-DOF dynamic model, a TCA model, and an analytical friction model. The friction model was validated through a comparison between the present analytical results and the predictions of a mixed EHL model proposed previously in terms of the friction coefficient and film thickness. Based on the presented results, the following conclusions were obtained: 1) The effects of contact paths on gear dynamics revealed a complicated nonlinear response in the vicinity of resonance, where the amplitudes of DTE of the middle and heel contacts exhibited significant jump discontinuities. Except for resonance, the DTE amplitudes, dynamic meshing force, and lateral and axial bearing forces of the toe contact fluctuated significantly during a wide speed sweep.
2) At resonance, the dynamic meshing force was small for the toe contact. However, the maximum Hertzian contact pressure was higher than those of the middle and heel contacts owing to the effect of contact geometry, causing high surface stress concentrations, which were closely related to surface micropitting and contact fatigue.
3) The friction coefficient and film thickness from the present analytical model agreed well with the results from a mixed EHL model of spiral bevel gears proposed previously. In addition, the line contact assumption for the conjugation of the spiral bevel gear appeared unreasonable owing to the significant prediction error of the friction coefficient at the neighbor of the pitch cone.
4) The dynamic efficiency was lower than the quasistatic efficiency, as expected, owing to the energy loss caused by the vibration of the gear shaft. At resonance, the efficiency fluctuated because of the tooth separations. The contact radii of the toe contact were relatively small, and correspondingly, the sliding velocity was relatively low, resulting in a high meshing efficiency for the toe contact. 5) A comparison of the numerical transmission efficiencies under different contact paths with the experimental measurements indicated good agreement. The tested efficiency was slightly smaller than the predicted values owing to the effect of the internal friction loss.