# On Dynamic Mesh Force Evaluation of Spiral Bevel Gears

Authors: Xiaoyu Sun; Yongqiang Zhao; Ming Liu; Yanping Liu
Journal: Shock and Vibration
 (2019)

---

## Abstract

The mesh model and mesh stiffness representation are the two main factors affecting the calculation method and the results of the dynamic mesh force. Comparative studies considering the two factors are performed to explore appropriate approaches to estimate the dynamic meshing load on each contacting tooth flank of spiral bevel gears. First, a tooth pair mesh model is proposed to better describe the mesh characteristics of individual tooth pairs in contact. The mesh parameters including the mesh vector, transmission error, and mesh stiffness are compared with those of the extensively applied single-point mesh model of a gear pair. Dynamic results from the proposed model indicate that it can reveal a more realistic and pronounced dynamic behavior of each engaged tooth pair. Second, dynamic mesh force calculations from three different approaches are compared to further investigate the effect of mesh stiffness representations. One method uses the mesh stiffness estimated by the commonly used average slope approach, the second method applies the mesh stiffness evaluated by the local slope approach, and the third approach utilizes a quasistatically defined interpolation function indexed by mesh deflection and mesh position.


---

## Body

## 1. Introduction

Spiral bevel gears are commonly used in power transmission between intersecting shafts but often suffer from fatigue failures caused by excess dynamic loads. The estimation of dynamic loads carried by each pair of teeth in contact lays the root for analyses of bending and contact stresses of gear teeth, gear tooth surface wear (pitting and scoring), and lubrication performance between the mating tooth flanks, which are all effective ways to investigate the mechanism of gear failures. Consequently, a reasonably accurate prediction of dynamic mesh force (DMF) is the cornerstone for failure analysis.Two factors have a dominant influence on the calculation approaches and results of DMF: firstly the mesh model which determines the calculation of mesh parameters including effective mesh point coordinates, line of action vectors, mesh stiffness, and transmission errors and secondly the mesh stiffness representations.The determination of mesh parameters for spiral bevel gears is complicated not only because of the complex tooth flank geometries but also because the mesh properties are strongly influenced by load conditions. Nowadays, loaded tooth contact analysis (LTCA) based on the conventional finite element method [1–4] or finite element/contact mechanism method [5] is widely applied to obtain the mesh parameters of spiral bevel gears under the loaded condition. However, in the absence of the sophisticated commercial finite element software (e.g., ABAQUS and ANSYS) and a specialized LTCA tool, [6] the loaded mesh parameters were not available, and hence, experimental or simple analytical models which preclude exact mesh geometries were used in most early dynamic studies [7, 8]. Until the beginning of this century, Cheng and Lim first attempted to integrate the time-variant mesh parameters obtained from tooth contact analysis (TCA) [9] and LTCA [10] of gears with exact tooth flank geometry into a mesh model. The model was extensively applied in the subsequent studies on dynamics of the hypoid and bevel gear pair or geared system [10–17]. In addition to the time-dependent mesh characteristics, some other special factors affecting the dynamic responses were investigated, such as backlash nonlinearity [10–13], friction [13–15], mesh stiffness asymmetry [12], gyroscopic effect [13], assembly errors [13], and geometric eccentricity [13] as well as other driveline components such as the bearing [16] and elastic housing [17].In all contributions mentioned above, the single-point coupling mesh model where the gear mesh coupling is modeled by a single spring-damper unit located at the effective mesh point and acting along the instantaneous line of action was applied. Although the model is believed to be able to catch the overall effect of the mesh characteristics of a gear pair, it cannot reflect the mesh properties of each mesh interface when more than one tooth pair is engaged in meshing, and hence, the DMF for each couple of teeth in contact cannot be predicted directly. Karagiannis et al. [15] utilized the load distribution factor defined and calculated in quasistatic LTCA to estimate the dynamic contact load per tooth when load sharing occurs. This method may be feasible, but the results cannot accurately reflect the real features of the contact force under the dynamic condition. To overcome the deficiency of the single-point coupling mesh model, Wang [18] proposed a multipoint coupling mesh model in which each mesh interface is represented by the respective effective mesh point, line of action, and mesh stiffness. The mesh properties and dynamic responses from the model were then compared with those from a mesh model whose parameters are obtained from pitch cone design [19]. The differences between them stem from the fact that the mesh positions and line of action vectors which are time-variant for the multipoint coupling mesh model are constant for the mesh model based on the pitch cone design concept. The authors also believed that the former is capable of describing the exact mesh properties of loaded gears, while the latter is only valid for no-load and light-loaded conditions, as the mesh parameters of them are obtained under loaded and unloaded conditions, respectively. So far, the multipoint coupling mesh model has not been used as widely as the single-point coupling mesh model for dynamic analysis probably because most studies focus on the overall effect of mesh characteristics on dynamics of a gear pair rather than the detailed dynamic responses of each meshing tooth pair. Among recent studies, maybe the multipoint coupling mesh model was only applied by Wang et al. [20]. Their research proved that the model has the capability of predicting the DMF-time history for each tooth pair within a full tooth engagement cycle. Even so, because the model still uses the transmission error of a gear pair to determine the tooth contact and contact deformation of each engaged tooth pair, the results cannot properly reflect the real contact state of each pair of teeth when multiple pairs of teeth are in the zone of contact.The accuracy of DMF calculation results is, to a large extent, dependent on the accuracy of mesh stiffness. The deformation of an engaged gear pair or an individual contacting tooth pair is nonlinear with the applied load because of the nonlinearity of the local contact deformation with the contact force. It means that mesh stiffness should vary with mesh deformation to accurately describe the nonlinear force-deflection relationship. However, in previous studies, although the variation of mesh stiffness with the mesh position was considered, the change in mesh stiffness with deformation at a given mesh position was rarely taken into consideration in DMF calculation. There are generally two methods to calculate mesh stiffness: the commonly used one is called the average slope approach and the other, which is seldom applied but may be more appropriate for dynamic analysis, is called the local slope approach [21]. The mesh stiffness evaluated by the two approaches can both be applied in DMF estimation, but, to accurately calculate the DMF, an error-prone point which should be paid enough attention is that the calculation formulas of DMF are entirely different for different types of mesh stiffness applied. It is worth mentioning that not all dynamic models utilize prespecified mesh stiffness for the calculation of DMF. In researches on the dynamics of parallel axis gears, [22, 23], for example, the mesh stiffness is evaluated according to instantaneous mesh deflection at each time step of the dynamic simulation. This kind of model is referred to as the coupled dynamics and contact analysis model, as the dynamic and contact behaviors are mutually influential. The variation of mesh stiffness with deformation is realized in this kind of model. However, the model is currently not computationally feasible for spiral bevel gears because of the more complex and time-consuming TCA computation.This study aims to explore appropriate methods for more accurate prediction of the DMF of each meshing tooth pair of spiral bevel gears through comparative analyses considering the two factors mentioned above. First, an improved mesh model of a tooth pair (MMTP) which represents the variation of the mesh properties of a tooth pair throughout a full engagement cycle was proposed. The model applies the transmission error curve of a tooth pair instead of that of a gear pair to estimate mesh stiffness and identify the contact for each couple of teeth. This improvement enables the model to capture more realistic mesh characteristics of each contacting tooth pair than the multipoint coupling mesh model. Second, we compared the dynamic responses from the proposed mesh model and the widely used single-point coupling mesh model for a generic spiral bevel gear pair on the one hand to validate the proposed one and on the other hand to demonstrate the superiority of the proposed model in the evaluation of DMF per tooth pair. It is shown that the two models predict the similar variation tendencies of DMF response with the mesh frequency. Although some differences exist in the response amplitudes particularly in a higher mesh frequency range, the resonance frequencies estimated by the two models are almost the same. Besides, the time histories of DMF per tooth pair predicted by the two models are significantly different, but the results from the proposed model can provide much better insight into dynamic contact behavior of each contacting tooth pair. Third, we performed a comparative study on three calculation approaches of DMF to further investigate the effect of mesh stiffness representations on DMF calculations. One method uses the mesh stiffness estimated by the average slope approach, another one applies the mesh stiffness evaluated by the local slope approach, and the third method utilizes a quasistatically defined interpolation function indexed by the mesh deflection and the pinion roll angle. Weighing the accuracy and cost of the computation, we found that the local slope approach is preferred for DMF evaluation of spiral bevel gears.


## 2. Mesh Model

This study applies the mesh model with prespecified mesh parameters obtained through quasistatic LTCA which is performed with the aid of the commercial finite element analysis software ABAQUS [24]. For the development of finite element models, the technique without the application of the CAD software proposed by Litvin et al. [1] is adopted. In this way, not only are the nodes of the finite element mesh lying on tooth surfaces guaranteed to be the points of the exact tooth surfaces (which are calculated through the mathematical model of tooth surface geometry [25]) but also the inaccuracy of solid models due to the geometric approximation by the spline functions of the CAD software can be avoided. The procedure of using ABAQUS for LTCA had been introduced in detail by Zhou et al. [4] and Hou and Deng [26] and thus will not be described here. The correctness of the finite element model is verified by comparing the obtained kinematic transmission error (KTE) curve with that obtained by traditional TCA, as shown in Figure 1, where the two curves are approximately consistent. A pair of mismatched face-milled Gleason spiral bevel gears is selected to be the example for analysis. The contact pattern is located in the center of the tooth surface, and there is no contact in the extreme part including the tip, toe, and heel of the teeth. Table 1 lists the design parameters, with which the machine tool setting parameters used to derive tooth geometry are obtained by Gleason CAGE4Win.Figure 1
Angular transmission error curves.Table 1
Design parameters of the example spiral bevel gear pair.
Parameters
Pinion (left hand)
Gear (right hand)
Number of teeth
11
34
Helical angle (rad)
0.646
Pressure angle (rad)
0.349
Pitch angle (rad)
0.313
1.257
Pitch diameter (mm)
71.5
221
Pitch cone distance (mm)
116.139
Width of the tooth (mm)
35
Outer transverse module
6.5
Virtual contact ratio
1.897For the convenience of the following discussions, the coordinate systems (Figure2) used in this study are introduced first. The LTCA is performed in the fixed reference frame Oxyz, whose origin is at the intersection point of the pinion and gear rotation axes. Other coordinate systems are used for dynamic analysis. The coordinate systems Opxp0yp0zp0 and Ogxg0yg0zg0 are the local inertial reference frames of the pinion and gear, respectively, whose origins and z-axes coincide with the mass centroids and rolling axes of each body under the initial static equilibrium condition. The body-fixed coordinate systems Cpxp4yp4zp4 and Cgxg4yg4zg4 for the pinion and gear are also located with their origins at the mass centroids, and all the coordinate axes are assumed to coincide with the principal axes of the inertia of each body.Figure 2
Coordinate systems of a spiral bevel gear pair.


### 2.1. Mesh Model of a Gear Pair (MMGP)

In the single-point coupling mesh model (Figure3), the two mating gears are coupled by a single spring-damper element whose stiffness km, damping cm, acting point, and direction vary with respect to the angular position of the pinion θp. The spring stiffness km, which is the function of bending, shear, and contact deformations of all meshing tooth pairs in the zone of contact and the base rotation compliance, represents the effective mesh stiffness of a gear pair. The KTE ek due to tooth surface modification, machining error, and assembly misalignment also changes with the pinion roll angle θp, while the backlash b is often regarded as a constant for simplicity.Figure 3
Single-point coupling mesh model.The synthesis of the single-point coupling mesh model requires reducing the distributed contact forces obtained by LTCA to an equivalent net mesh forceFmFmx,Fmy,Fmz (Figure 4(a)). The direction of the net mesh force is defined as the line of action vector Lm which points from the pinion to the gear, and the point of action is the effective mesh point representing the mean contact position. The calculation method of contact parameters applied in previous studies [10–14] can be used here, while if LTCA is undertaken using ABAQUS, another similar but more convenient approach can be applied. As long as the history outputs of CFN, CMN, and XN for the contact pair are requested [24], the effective mesh force Fm,iFmx,i,Fmy,i,Fmz,i and the mesh moment Mm,iMmx,i,Mmy,i,Mmz,i due to contact pressure on each meshing tooth pair i and the mesh force acting point vector rm,ixm,i,ym,i,zm,i can be obtained directly at each specified time step. Thus, the magnitude of equivalent net mesh force Fnet is calculated by(1)Fnet=Fmx2+Fmy2+Fmz2,Fml=∑i=1NFml,i,l=x,y,z.Figure 4
Synthesis of the effective mesh parameters of the (a) single-point coupling mesh model and (b) multipoint coupling mesh model.
(a)
(b)The line of action vectorLmnx,ny,nz of the net mesh force is obtained by(2)nl=FmlFnet,l=x,y,z.The total contact-induced moments about each coordinate axis are as follows:(3)Mml=∑i=1NMml,i,l=x,y,z.The effective mesh point positionrmxm,ym,zm can be calculated from the following equation:(4)xm=∑i=1Nxm,iFm,i∑i=1NFm,i,zm=Mmy+FmzxmFmx,ym=Mmx+FmyzmFmz.In above-mentioned equations,N is the number of engaged tooth pairs at each pinion roll angle. In this work, the single-point coupling mesh model is named the mesh model of a gear pair (MMGP), as it reflects the overall mesh properties of a gear pair.


### 2.2. Mesh Model of a Tooth Pair (MMTP)

In the multipoint coupling mesh model (Figure5(b)) [18], there are a variable number of spring-damper units coupling the pinion and gear in each mesh position. The spring-damper element represents the effective mesh point, line of action, and mesh stiffness of each mesh interface, and the number of elements which depends on the gear rotational position and load torque indicates the number of tooth pairs in the zone of contact. This model describes the variation of the mesh characteristics of each tooth pair in the contact zone over a meshing period, and thereby, the mesh characteristics of the individual tooth pair can be considered the basic unit. Based on this point of view, the mesh model of a tooth pair (MMTP) (Figure 5(a)) which describes the mesh characteristic variation of a tooth pair over one tooth engagement cycle is proposed. The MMTP facilitates the development of models handling more complex situations, such as the model with nonidentical mesh characteristics of the tooth pair (e.g., one or several teeth have cracks, surface wear, or any other features influencing mesh properties).Figure 5
Illustrations of the (a) mesh model of a tooth pair and (b) multipoint coupling mesh model.
(a)
(b)The synthesis of the effective mesh point and line of action for the MMTP is the same as that for the multipoint coupling mesh model (Figure4(b)) [18]. However, as mentioned before, the effective mesh force Fm,iFmx,i,Fmy,i,Fmz,i due to contact pressure on each meshing tooth pair i and the corresponding acting point vector rm,ixm,i,ym,i,zm,i can be obtained directly from LTCA performed with ABAQUS. The direction of the effective mesh force is defined as the line of action for each contacting tooth pair, and the acting point vector is just the effective mesh point position vector. The calculation approach of mesh stiffness and DMF for the MMTP is modified to make the model describe the meshing state of each pair of teeth more accurately than the multipoint coupling mesh model. These will be elaborated in the following corresponding sections.


### 2.3. Comparison of Mesh Parameters

In addition to comparing the mesh properties embodied on the individual tooth pair and the gear pair, the effect of the applied load on their variations was investigated. To this end, we synthesized the mesh parameters for the MMTP and the MMGP at 50 Nm and 300 Nm torque loads. The mesh vectors, which represent the magnitude, direction, and point of action of the equivalent mesh force, at each mesh position throughout one tooth engagement cycle for the MMTP and over one mesh cycle for the MMGP are presented in Figures6(a) and 6(b), respectively. As can be seen, the magnitude of the equivalent mesh force for the MMTP (i.e., the effective mesh force on a single pair of meshing teeth) varies obviously with the mesh position, while only the tiny change can be observed in the size of the equivalent net mesh force for the MMGP (i.e., the combined mesh force of all engaged tooth pairs). The former is mainly attributed to the load sharing among the tooth pairs in contact simultaneously, whereas the latter is only caused by small changes in the directional rotation radius (which can be regarded as the arm of force) at different mesh positions. Figure 6 also shows that the direction of the equivalent mesh force changes little for both the MMTP and the MMGP.Figure 6
Mesh vectors of the MMTP (a) and the MMGP (b).
(a)
(b)Furthermore, the line of action directional cosines and effective mesh point coordinates of the two models over a full tooth (tooth pair 2) engagement cycle are compared in Figure7. It shows clearly that all mesh parameters of the tooth pair are altered significantly with the gear rotation relative to the corresponding changes in the mesh parameters of the gear pair. In addition, from both Figures 6 and 7, we can observe that the range of change in mesh parameters of the gear pair becomes small under the heavy-loaded condition (300 Nm), but no obvious change for the mesh parameters of the tooth pair was observed. A reasonable explanation is as follows: with the increase of load, the extended tooth contact appears because of the deformation of gear teeth. It is manifested in the advance engagement of the incoming tooth pair 3 and the deferred disengagement of the receding tooth pair 1 in Figure 7. This change increases the contact ratio and meanwhile decreases the amplitude of the variation of the mesh parameters for the gear pair. Note that it does not influence the mesh frequency of the gear pair because according to the compatibility condition, the tooth rotation due to deformation for each couple of teeth should be the same to maintain tooth contact and continuity in power transmission, and hence, the tooth pitch does not change. Overall, the results indicate that although the time-varying mesh characteristics exhibited on each pair of teeth are prominent, their overall effect on the gear pair is weakened, and this weakening will become more pronounced as the contact ratio increases.Figure 7
Comparison of time-varying line of action directional cosines (a–c) and effective mesh point coordinates (d–f) between the MMTP and the MMGP at 50 Nm and 300 Nm torque loads.
(a)
(b)
(c)
(d)
(e)
(f)


### 2.4. Remark on the Mesh Model Based on Quasistatic LTCA

From the above analysis, it is substantiated again that the mesh parameters of spiral bevel gears are load-dependent. Therefore, the mesh model based on quasistatic LTCA is superior to the mesh model relying on the geometric relationship in terms of describing the mesh properties of loaded gear pairs. However, as the mesh parameters are estimated at a constant applied torque, the mesh model is theoretically only suitable for the dynamic analysis of a gear pair subjected to this torque load, and it seems to be incompetent in the study of a gear pair subjected to a time-variant torque load. Nevertheless, it has been discovered that the mesh stiffness shows more prominent time-varying parametric excitation effects than the mesh vector [11, 13]. Therefore, a compromising way to calculate the DMF is applying the mesh vector predicted at the mean applied torque and the calculation approach using the mesh force interpolation function (as introduced in DMF Calculation) which takes the load effects on mesh stiffness into consideration.


## 2.1. Mesh Model of a Gear Pair (MMGP)

In the single-point coupling mesh model (Figure3), the two mating gears are coupled by a single spring-damper element whose stiffness km, damping cm, acting point, and direction vary with respect to the angular position of the pinion θp. The spring stiffness km, which is the function of bending, shear, and contact deformations of all meshing tooth pairs in the zone of contact and the base rotation compliance, represents the effective mesh stiffness of a gear pair. The KTE ek due to tooth surface modification, machining error, and assembly misalignment also changes with the pinion roll angle θp, while the backlash b is often regarded as a constant for simplicity.Figure 3
Single-point coupling mesh model.The synthesis of the single-point coupling mesh model requires reducing the distributed contact forces obtained by LTCA to an equivalent net mesh forceFmFmx,Fmy,Fmz (Figure 4(a)). The direction of the net mesh force is defined as the line of action vector Lm which points from the pinion to the gear, and the point of action is the effective mesh point representing the mean contact position. The calculation method of contact parameters applied in previous studies [10–14] can be used here, while if LTCA is undertaken using ABAQUS, another similar but more convenient approach can be applied. As long as the history outputs of CFN, CMN, and XN for the contact pair are requested [24], the effective mesh force Fm,iFmx,i,Fmy,i,Fmz,i and the mesh moment Mm,iMmx,i,Mmy,i,Mmz,i due to contact pressure on each meshing tooth pair i and the mesh force acting point vector rm,ixm,i,ym,i,zm,i can be obtained directly at each specified time step. Thus, the magnitude of equivalent net mesh force Fnet is calculated by(1)Fnet=Fmx2+Fmy2+Fmz2,Fml=∑i=1NFml,i,l=x,y,z.Figure 4
Synthesis of the effective mesh parameters of the (a) single-point coupling mesh model and (b) multipoint coupling mesh model.
(a)
(b)The line of action vectorLmnx,ny,nz of the net mesh force is obtained by(2)nl=FmlFnet,l=x,y,z.The total contact-induced moments about each coordinate axis are as follows:(3)Mml=∑i=1NMml,i,l=x,y,z.The effective mesh point positionrmxm,ym,zm can be calculated from the following equation:(4)xm=∑i=1Nxm,iFm,i∑i=1NFm,i,zm=Mmy+FmzxmFmx,ym=Mmx+FmyzmFmz.In above-mentioned equations,N is the number of engaged tooth pairs at each pinion roll angle. In this work, the single-point coupling mesh model is named the mesh model of a gear pair (MMGP), as it reflects the overall mesh properties of a gear pair.


## 2.2. Mesh Model of a Tooth Pair (MMTP)

In the multipoint coupling mesh model (Figure5(b)) [18], there are a variable number of spring-damper units coupling the pinion and gear in each mesh position. The spring-damper element represents the effective mesh point, line of action, and mesh stiffness of each mesh interface, and the number of elements which depends on the gear rotational position and load torque indicates the number of tooth pairs in the zone of contact. This model describes the variation of the mesh characteristics of each tooth pair in the contact zone over a meshing period, and thereby, the mesh characteristics of the individual tooth pair can be considered the basic unit. Based on this point of view, the mesh model of a tooth pair (MMTP) (Figure 5(a)) which describes the mesh characteristic variation of a tooth pair over one tooth engagement cycle is proposed. The MMTP facilitates the development of models handling more complex situations, such as the model with nonidentical mesh characteristics of the tooth pair (e.g., one or several teeth have cracks, surface wear, or any other features influencing mesh properties).Figure 5
Illustrations of the (a) mesh model of a tooth pair and (b) multipoint coupling mesh model.
(a)
(b)The synthesis of the effective mesh point and line of action for the MMTP is the same as that for the multipoint coupling mesh model (Figure4(b)) [18]. However, as mentioned before, the effective mesh force Fm,iFmx,i,Fmy,i,Fmz,i due to contact pressure on each meshing tooth pair i and the corresponding acting point vector rm,ixm,i,ym,i,zm,i can be obtained directly from LTCA performed with ABAQUS. The direction of the effective mesh force is defined as the line of action for each contacting tooth pair, and the acting point vector is just the effective mesh point position vector. The calculation approach of mesh stiffness and DMF for the MMTP is modified to make the model describe the meshing state of each pair of teeth more accurately than the multipoint coupling mesh model. These will be elaborated in the following corresponding sections.


## 2.3. Comparison of Mesh Parameters

In addition to comparing the mesh properties embodied on the individual tooth pair and the gear pair, the effect of the applied load on their variations was investigated. To this end, we synthesized the mesh parameters for the MMTP and the MMGP at 50 Nm and 300 Nm torque loads. The mesh vectors, which represent the magnitude, direction, and point of action of the equivalent mesh force, at each mesh position throughout one tooth engagement cycle for the MMTP and over one mesh cycle for the MMGP are presented in Figures6(a) and 6(b), respectively. As can be seen, the magnitude of the equivalent mesh force for the MMTP (i.e., the effective mesh force on a single pair of meshing teeth) varies obviously with the mesh position, while only the tiny change can be observed in the size of the equivalent net mesh force for the MMGP (i.e., the combined mesh force of all engaged tooth pairs). The former is mainly attributed to the load sharing among the tooth pairs in contact simultaneously, whereas the latter is only caused by small changes in the directional rotation radius (which can be regarded as the arm of force) at different mesh positions. Figure 6 also shows that the direction of the equivalent mesh force changes little for both the MMTP and the MMGP.Figure 6
Mesh vectors of the MMTP (a) and the MMGP (b).
(a)
(b)Furthermore, the line of action directional cosines and effective mesh point coordinates of the two models over a full tooth (tooth pair 2) engagement cycle are compared in Figure7. It shows clearly that all mesh parameters of the tooth pair are altered significantly with the gear rotation relative to the corresponding changes in the mesh parameters of the gear pair. In addition, from both Figures 6 and 7, we can observe that the range of change in mesh parameters of the gear pair becomes small under the heavy-loaded condition (300 Nm), but no obvious change for the mesh parameters of the tooth pair was observed. A reasonable explanation is as follows: with the increase of load, the extended tooth contact appears because of the deformation of gear teeth. It is manifested in the advance engagement of the incoming tooth pair 3 and the deferred disengagement of the receding tooth pair 1 in Figure 7. This change increases the contact ratio and meanwhile decreases the amplitude of the variation of the mesh parameters for the gear pair. Note that it does not influence the mesh frequency of the gear pair because according to the compatibility condition, the tooth rotation due to deformation for each couple of teeth should be the same to maintain tooth contact and continuity in power transmission, and hence, the tooth pitch does not change. Overall, the results indicate that although the time-varying mesh characteristics exhibited on each pair of teeth are prominent, their overall effect on the gear pair is weakened, and this weakening will become more pronounced as the contact ratio increases.Figure 7
Comparison of time-varying line of action directional cosines (a–c) and effective mesh point coordinates (d–f) between the MMTP and the MMGP at 50 Nm and 300 Nm torque loads.
(a)
(b)
(c)
(d)
(e)
(f)


## 2.4. Remark on the Mesh Model Based on Quasistatic LTCA

From the above analysis, it is substantiated again that the mesh parameters of spiral bevel gears are load-dependent. Therefore, the mesh model based on quasistatic LTCA is superior to the mesh model relying on the geometric relationship in terms of describing the mesh properties of loaded gear pairs. However, as the mesh parameters are estimated at a constant applied torque, the mesh model is theoretically only suitable for the dynamic analysis of a gear pair subjected to this torque load, and it seems to be incompetent in the study of a gear pair subjected to a time-variant torque load. Nevertheless, it has been discovered that the mesh stiffness shows more prominent time-varying parametric excitation effects than the mesh vector [11, 13]. Therefore, a compromising way to calculate the DMF is applying the mesh vector predicted at the mean applied torque and the calculation approach using the mesh force interpolation function (as introduced in DMF Calculation) which takes the load effects on mesh stiffness into consideration.


## 3. Mesh Stiffness Evaluation

The force-deflection curve (a typical one is illustrated in Figure8) describes the nonlinear relationship between the force and the deformation precisely. Therefore, using this curve to predict the DMF for a known mesh deformation could be the most accurate approach. However, it is by no means the most efficient method since the establishment of the curve requires a series of time-consuming LTCAs at different load levels. Usually, the DMF for a given torque load is predicted by using the mesh stiffness evaluated at the same load. Note that although at each rolling position the mesh stiffness is a constant, it varies with respect to mesh positions. Hence, the mesh stiffness has a spatially varying feature.Figure 8
A typical fitted force-deflection curve and illustration of the local tangent mesh stiffness (LTMS) and the average secant mesh stiffness (ASMS).In general, there are two calculation approaches for mesh stiffness, namely, the average slope approach and the local slope approach [21]. The corresponding calculated mesh stiffness is named the average secant mesh stiffness (ASMS) and local tangent mesh stiffness (LTMS) in this paper. As shown in Figure 8, they are the slope of the secant (blue dashed line) and the slope of the tangent (red dash-dotted line) of the force-deflection curve, respectively.


### 3.1. Average Secant Mesh Stiffness

The average secant mesh stiffness has been extensively used to calculate DMF for different types of gears, such as spur gears [27–29], helical gears [29, 30], and spiral bevel or hypoid gears [9, 10, 13–18, 20]. It is calculated through dividing the magnitude of mesh force Fm by the mesh deflection along the line of action δ:(5)kma=Fmδ.The difference between the translational KTEε0 and the loaded transmission error (LTE) ε can be attributed to the mesh deflection δ:(6)δ=ε0−ε.In practice, the KTE always exists for spiral bevel gears, even if there are no machining and assembly errors, because gear tooth flank profiles are often modified to be mismatched to avoid contact of extreme parts (top, heel, or toe) of the teeth and to reduce the sensitivity of transmission performance to installation and machining errors. Because of the mismatch, the basic mating equation of the contacting tooth flanks is satisfied only at one point of the path of contact. It has been proved that although the KTE is always small compared to the LTE, the KTE cannot be omitted in the calculation of mesh stiffness, especially for the light-loaded condition [13].The translational KTE and LTE,ε0 and ε, are the projections of the corresponding angular transmission errors, e0 and e, along the line of action. As y-axis of the coordinate system Oxyz is the gear-shaft rotation axis (Figure 2), at which the angular transmission error is expressed about, the translational KTE and LTE are calculated by(7)ε0=e0λy,ε=eλy,where λy=nxzm−nzxm is the directional rotation radius for the angular transmission error. Note that equations (5)–(7) are derived based on the assumption that the mesh vectors for unloaded and loaded conditions are identical. However, the results shown in Figure 7 indicate that the hypothesis deviates from the fact and thereby inevitably brings about calculation errors.The angular KTE and LTE,e0 and e, can be determined as(8)e0=Δθg0−NpNgΔθp0,e=Δθg−NpNgΔθp,where Δθp0 and Δθg0 (obtained by TCA) are the unloaded rotational displacements of the pinion and gear relative to the initial conjugated contact position where the transmission error equals zero, Δθp and Δθg (obtained by LTCA) are the corresponding loaded rotational displacements, and Np and Ng are the tooth numbers of the pinion and gear, respectively.Equations (5)–(8) can be directly applied to calculate the ASMS of a gear pair. In this case, Fm is the magnitude of the net mesh force of the gear pair and δ is the total deformation of all meshing tooth pairs in the zone of contact.For the MMTP, the ASMS of a single tooth pair needs to be estimated. The effective mesh force on each meshing tooth pairFm,i can be obtained directly through LTCA, while the troublesome issue is how to obtain the mesh deflection of a single tooth pair when multiple sets of teeth are in contact. In this work, the compatibility condition is applied to address this problem. It is based on the fact that the total tooth rotation Δθ (which is the sum of the composite deformation of the tooth and base ΔB caused by bending moment and shearing load, contact deformation ΔH, and profile separation ΔS divided by the directional rotation radius about the rolling axis λ) must be the same for all simultaneously engaged tooth pairs to maintain the tooth contact and continuity in power transmission [2]. If three pairs of teeth are in contact, the compatibility condition can be expressed by(9)ΔSi−1+ΔBi−1+ΔHi−1λi−1=ΔSi+ΔBi+ΔHiλi=ΔSi+1+ΔBi+1+ΔHi+1λi+1=Δθ.If the ASMS of a tooth pair is defined as the mesh force required to produce a unit total deformation (which contains the composite deformation of the tooth and baseΔB and the contact deformation ΔH), equation (9) can be rewritten in a form with the ASMS of each meshing tooth pair kma,i as follows:(10)ΔSi−1+Fm,i−1/kma,i−1λi−1=ΔSi+Fm,i/kma,iλi=ΔSi+1+Fm,i+1/kma,i+1λi+1=Δθ.


### 3.2. Local Tangent Mesh Stiffness

The local tangent mesh stiffness can be approximately calculated by the central difference method which requires the LTCA to be performed at three different torque loads: one at the nominal torque load, another above it, and the third below it:(11)kml=dFmdδδ=δT≈12FmT−FmT−ΔTδT−δT−ΔT+FmT+ΔT−FmTδT+ΔT−δT,where ΔT is a specified small change in torque load. As δ=ε0−ε, equation (11) can be transformed to(12)kml≈12FmT−FmT−ΔTεT−ΔT−εT+FmT+ΔT−FmTεT−εT+ΔT.Equation (12) can be directly used to calculate the LTMS of both a gear pair and a tooth pair, as the profile separation of each tooth pair ΔSi which influences the calculation of the ASMS of a tooth pair does not function in the calculation of the LTMS of a tooth pair. The reason is that the KTE ε0 which contains the profile separation does not appear in equation (12). In contrast with the calculation of the ASMS, for the evaluation of the LTMS, assuming the mesh vectors at the torque loads of T, T+ΔT, and T−ΔT are identical is reasonable when ΔT is small.


### 3.3. Example

Taking the evaluation of the mesh stiffness at a torque load of 300 Nm as an example, we demonstrate the general calculation approaches of the ASMS and LTMS. Figure9 shows the effective mesh force of three consecutively engaged tooth pairs and the net mesh force of the gear pair over one tooth engagement cycle. Note that the actual length of the period of one tooth engagement is longer than the one shown in Figure 9 since tooth pair 2 comes into contact at some pinion roll angle less than 0.171 radians and leaves the contact at some pinion roll angle greater than 1.107 radians. By assuming that, at the beginning and end of a tooth engagement cycle, the mesh force of each meshing tooth pair increases and decreases linearly, the pinion roll angles where tooth contact begins and ends can be approximately predicted through linear extrapolation. Using the data given in Figure 9, we can predict that tooth pair 2 engages at approximately 0.159 radians and disengages at nearly 1.110 radians, and tooth pair 3 engages at about 0.730 radians. Therefore, if the pitch error is not considered, the tooth meshing is repeated every 0.571 radians, and the actual length of one tooth engagement cycle is 0.951 radians. With the known roll angles at the beginning of engagement and disengagement, mesh parameters at the two points can be predicted. In this work, the LS-SVMlab toolbox [31] is utilized to estimate the functions of mesh parameters with respect to the rotation angle of the pinion tooth, and then the parameters at the starting positions of engagement and disengagement are evaluated by the estimated functions.Figure 9
Finite element calculation of the effective mesh force of the tooth pair and the net mesh force of the gear pair over one tooth engagement cycle.Figure1 shows the angular KTE and LTE of the example gear pair. The KTE of the gear pair is obtained by LTCA performed at a small torque load. This method is named the FEM-based TCA. The torque load should be carefully chosen (3 Nm in this work) to keep the teeth in contact but not to cause appreciable deformation. The KTEs of tooth pairs 1, 2, and 3, which are necessary for the calculation of the mesh deflection of the tooth pair, are acquired by traditional TCA [32] (to the authors’ knowledge, the complete KTE curve of a tooth pair cannot be obtained by the FEM-based TCA). The appropriateness of the selected torque for FEM-based TCA can be validated if the obtained KTE is approximately consistent with KTEs acquired by traditional TCA. It can be seen from Figure 1 that, in contact zone 1, the KTE (absolute value) of tooth pair 1 is less than that of tooth pair 2, which means in this region tooth pair 1 is in contact, whereas tooth pair 2 is not in contact in the no-loading condition. It is because the mating tooth surfaces of tooth pair 2 are separated from each other by the profile separation which can be measured by the difference in KTEs between tooth pairs 1 and 2. Likewise, for contact zones 2 and 3, only tooth pairs 2 and 3 are in contact, respectively, in the unloaded condition. That is why the KTEs of the gear pair in contact zones 1, 2, and 3 correspond to the KTEs of tooth pairs 1, 2, and 3, respectively.The total deformation of a gear pair is usually considered the difference between the KTE and the LTE of the gear pair. This difference, however, is not always the mesh deflection of a tooth pair caused by the load it bears. The schematic diagram for illustrating the deflection of the meshing tooth pair is shown in Figure10. As previously mentioned, at any rolling position in contact zone 1 (shown in Figure 1), only tooth pair 1 is in contact in the no-loading condition and the mating tooth surfaces of tooth pair 2 are separated by the profile separation, as shown in Figure 10. With the augment of torque load, the LTE of the gear pair increases gradually because of the mesh deflection of tooth pair 1. If the LTE of the gear pair at any rolling position in contact zone 1 is greater than the KTE of tooth pair 2, which means the profile separation of tooth pair 2 at this position has been closed because of tooth rotation, tooth pair 2 will be in contact and share the whole load with tooth pair 1. Therefore, as illustrated in Figure 10, part of the rotation of tooth pair 2 compensates the profile separation between the mating flanks and the remaining part is the result of mesh deflection. At any rolling position in contact zone 2, as tooth pair 2 is in contact initially, the tooth rotation of tooth pair 2 is wholly attributed to its deflection. In contact zone 3, however, the profile separation of tooth pair 2 must be closed before contact, and hence, the remaining part of the tooth rotation which excludes the angular profile separation is the rotational deformation of tooth pair 2. Based on above analysis, it can be speculated that, at the torque load of 300 Nm, tooth pairs 1 and 2 are in contact in region AB, only tooth pair 2 is in contact in region BC, and tooth pairs 2 and 3 are in contact in region CD. This speculation is proved to be true by the time histories of mesh force on each engaged tooth pair shown in Figure 9.Figure 10
Deflection of the meshing tooth pair at the meshing position in contact zone 1.Above discussions elaborate the meaning of equation (10). Namely, the tooth rotations Δθ for all engaged tooth pairs are the same and equal the rotation of the pinion and gear blanks which is measured by the difference between the angular KTE and LTE of the gear pair, while the mesh deflection along the line of action δi for the tooth pair i due to the load it shares should be calculated by subtracting the profile separation ΔSi from the projection of the tooth rotation Δθ along the line of action. Actually, the mesh deflection for a tooth pair equals the difference between the KTE of the tooth pair and the LTE of the gear pair, as shown in Figure 1. Therefore, in the MMTP, the KTE of a tooth pair is applied instead of the KTE of a gear pair in the estimation of mesh stiffness and the identification of the contact for each pair of teeth. It is the main difference from the multipoint coupling mesh model.Figure11 shows the mesh deflections along the line of action for the gear pair and each tooth pair. The tooth engagement cycle is divided into regions where two meshing tooth pairs are in contact and one meshing tooth pair is in contact (denoted by DTC and STC, respectively, in Figures 11–13). In the double-tooth contact zone, although the rotational deformation of the gear pair equals the rotational deformation of the tooth pair which is initially in contact in the unloaded condition, their translational deflections along the line of action are somewhat different. This is because in the double-tooth contact zone, the effective mesh point of the gear pair is different from that of the tooth pair, as shown in Figure 4, and thus, the corresponding lines of action of the net mesh force acting on the effective mesh point are different as well.Figure 11
Mesh deflection along the line of action of the gear pair and the tooth pair over one tooth engagement cycle at 300 Nm applied torque.Figure 12
The ASMS of the gear pair and the tooth pair over one tooth engagement cycle at 300 Nm applied torque.Figure 13
The LTMS of the gear pair and the tooth pair over one tooth engagement cycle at 300 Nm applied torque.Figure12 compares the ASMS calculations of the gear pair and the tooth pair. It can be seen that, in the double-tooth contact zone, the ASMS of the gear pair experiences a gradual increase and then a steady decrease. The variation significantly differs from that of unmodified spur gears whose mesh stiffness presents a sharp increase and decrease at the instants when the approaching tooth pair comes into contact and the receding tooth pair leaves the mesh, respectively [21]. Apart from that, with the contact ratio increase, the variation of the mesh stiffness of the gear pair tends to become smaller. This conclusion can be proved by the comparison of the ASMS calculations between this study and studies of Peng [13] and Wang et al. [20].An important fact shown in Figure12 is that, in the double-tooth contact zone, the ASMS of a gear pair does not equal the sum of the ASMS of each meshing tooth pair, if the KTE of a tooth pair is applied to estimate the mesh stiffness of a tooth pair, except for the position where the two meshing tooth pairs are in contact in the unloaded condition, such as point A in the graph. The reason is explained in detail in Appendix A. In many previous studies on mesh stiffness calculation of spur gears [33–36], the total mesh stiffness of the gear pair in the multiple-tooth engagement region is calculated by summation of the mesh stiffness of all tooth pairs in contact. This estimation method is probably inapplicable for spiral bevel or hypoid gears in terms of the ASMS evaluation, whereas it is proved to be feasible for the LTMS calculation.For the calculation of the LTMS expressed in equation (12), LTCA has to be performed at the torque loads of T, T+ΔT, and T−ΔT to obtain the corresponding LTEs εT, εT+ΔT, and εT−ΔT. Since the mesh parameters are all load-dependent, ΔT should be properly chosen to ensure that the effective mesh points and the directions of the line of action are almost the same for the three torque loads at each rolling position. However, it does not mean that the smaller the ΔT is, the better the result will be since the effect of numerical calculation errors on the calculation results will become more pronounced as ΔT decreases. In the present study, ΔT=10Nm can weigh the two factors mentioned above. Figure 13 compares the LTMS calculations of the gear pair and the tooth pair. It can be seen that, in the double-tooth contact zone, the LTMS of the gear pair approximately equals the sum of the LTMS of simultaneously contacting tooth pairs. The reason is mathematically clarified in Appendix A. Comparing calculations in Figures 12 and 13, we can observe that, for both the mesh stiffness of a gear pair and a tooth pair, the LTMS is larger than the ASMS over the entire tooth engagement cycle. It is consistent with the theoretical prediction from the concave force-deflection curve shown in Figure 8. It has been known that the nonlinear contact deformation due to the growth of the contact area as applied load increases results in the concave shape of the force-deflection curve. Note that the difference between the computed results of the ASMS and LTMS is the result of their definitions and calculation approaches, and we cannot make a conclusion about which prediction is more accurate. Both the ASMS and LTMS can be used to calculate DMF, but the corresponding calculation formulas are entirely different.


## 3.1. Average Secant Mesh Stiffness

The average secant mesh stiffness has been extensively used to calculate DMF for different types of gears, such as spur gears [27–29], helical gears [29, 30], and spiral bevel or hypoid gears [9, 10, 13–18, 20]. It is calculated through dividing the magnitude of mesh force Fm by the mesh deflection along the line of action δ:(5)kma=Fmδ.The difference between the translational KTEε0 and the loaded transmission error (LTE) ε can be attributed to the mesh deflection δ:(6)δ=ε0−ε.In practice, the KTE always exists for spiral bevel gears, even if there are no machining and assembly errors, because gear tooth flank profiles are often modified to be mismatched to avoid contact of extreme parts (top, heel, or toe) of the teeth and to reduce the sensitivity of transmission performance to installation and machining errors. Because of the mismatch, the basic mating equation of the contacting tooth flanks is satisfied only at one point of the path of contact. It has been proved that although the KTE is always small compared to the LTE, the KTE cannot be omitted in the calculation of mesh stiffness, especially for the light-loaded condition [13].The translational KTE and LTE,ε0 and ε, are the projections of the corresponding angular transmission errors, e0 and e, along the line of action. As y-axis of the coordinate system Oxyz is the gear-shaft rotation axis (Figure 2), at which the angular transmission error is expressed about, the translational KTE and LTE are calculated by(7)ε0=e0λy,ε=eλy,where λy=nxzm−nzxm is the directional rotation radius for the angular transmission error. Note that equations (5)–(7) are derived based on the assumption that the mesh vectors for unloaded and loaded conditions are identical. However, the results shown in Figure 7 indicate that the hypothesis deviates from the fact and thereby inevitably brings about calculation errors.The angular KTE and LTE,e0 and e, can be determined as(8)e0=Δθg0−NpNgΔθp0,e=Δθg−NpNgΔθp,where Δθp0 and Δθg0 (obtained by TCA) are the unloaded rotational displacements of the pinion and gear relative to the initial conjugated contact position where the transmission error equals zero, Δθp and Δθg (obtained by LTCA) are the corresponding loaded rotational displacements, and Np and Ng are the tooth numbers of the pinion and gear, respectively.Equations (5)–(8) can be directly applied to calculate the ASMS of a gear pair. In this case, Fm is the magnitude of the net mesh force of the gear pair and δ is the total deformation of all meshing tooth pairs in the zone of contact.For the MMTP, the ASMS of a single tooth pair needs to be estimated. The effective mesh force on each meshing tooth pairFm,i can be obtained directly through LTCA, while the troublesome issue is how to obtain the mesh deflection of a single tooth pair when multiple sets of teeth are in contact. In this work, the compatibility condition is applied to address this problem. It is based on the fact that the total tooth rotation Δθ (which is the sum of the composite deformation of the tooth and base ΔB caused by bending moment and shearing load, contact deformation ΔH, and profile separation ΔS divided by the directional rotation radius about the rolling axis λ) must be the same for all simultaneously engaged tooth pairs to maintain the tooth contact and continuity in power transmission [2]. If three pairs of teeth are in contact, the compatibility condition can be expressed by(9)ΔSi−1+ΔBi−1+ΔHi−1λi−1=ΔSi+ΔBi+ΔHiλi=ΔSi+1+ΔBi+1+ΔHi+1λi+1=Δθ.If the ASMS of a tooth pair is defined as the mesh force required to produce a unit total deformation (which contains the composite deformation of the tooth and baseΔB and the contact deformation ΔH), equation (9) can be rewritten in a form with the ASMS of each meshing tooth pair kma,i as follows:(10)ΔSi−1+Fm,i−1/kma,i−1λi−1=ΔSi+Fm,i/kma,iλi=ΔSi+1+Fm,i+1/kma,i+1λi+1=Δθ.


## 3.2. Local Tangent Mesh Stiffness

The local tangent mesh stiffness can be approximately calculated by the central difference method which requires the LTCA to be performed at three different torque loads: one at the nominal torque load, another above it, and the third below it:(11)kml=dFmdδδ=δT≈12FmT−FmT−ΔTδT−δT−ΔT+FmT+ΔT−FmTδT+ΔT−δT,where ΔT is a specified small change in torque load. As δ=ε0−ε, equation (11) can be transformed to(12)kml≈12FmT−FmT−ΔTεT−ΔT−εT+FmT+ΔT−FmTεT−εT+ΔT.Equation (12) can be directly used to calculate the LTMS of both a gear pair and a tooth pair, as the profile separation of each tooth pair ΔSi which influences the calculation of the ASMS of a tooth pair does not function in the calculation of the LTMS of a tooth pair. The reason is that the KTE ε0 which contains the profile separation does not appear in equation (12). In contrast with the calculation of the ASMS, for the evaluation of the LTMS, assuming the mesh vectors at the torque loads of T, T+ΔT, and T−ΔT are identical is reasonable when ΔT is small.


## 3.3. Example

Taking the evaluation of the mesh stiffness at a torque load of 300 Nm as an example, we demonstrate the general calculation approaches of the ASMS and LTMS. Figure9 shows the effective mesh force of three consecutively engaged tooth pairs and the net mesh force of the gear pair over one tooth engagement cycle. Note that the actual length of the period of one tooth engagement is longer than the one shown in Figure 9 since tooth pair 2 comes into contact at some pinion roll angle less than 0.171 radians and leaves the contact at some pinion roll angle greater than 1.107 radians. By assuming that, at the beginning and end of a tooth engagement cycle, the mesh force of each meshing tooth pair increases and decreases linearly, the pinion roll angles where tooth contact begins and ends can be approximately predicted through linear extrapolation. Using the data given in Figure 9, we can predict that tooth pair 2 engages at approximately 0.159 radians and disengages at nearly 1.110 radians, and tooth pair 3 engages at about 0.730 radians. Therefore, if the pitch error is not considered, the tooth meshing is repeated every 0.571 radians, and the actual length of one tooth engagement cycle is 0.951 radians. With the known roll angles at the beginning of engagement and disengagement, mesh parameters at the two points can be predicted. In this work, the LS-SVMlab toolbox [31] is utilized to estimate the functions of mesh parameters with respect to the rotation angle of the pinion tooth, and then the parameters at the starting positions of engagement and disengagement are evaluated by the estimated functions.Figure 9
Finite element calculation of the effective mesh force of the tooth pair and the net mesh force of the gear pair over one tooth engagement cycle.Figure1 shows the angular KTE and LTE of the example gear pair. The KTE of the gear pair is obtained by LTCA performed at a small torque load. This method is named the FEM-based TCA. The torque load should be carefully chosen (3 Nm in this work) to keep the teeth in contact but not to cause appreciable deformation. The KTEs of tooth pairs 1, 2, and 3, which are necessary for the calculation of the mesh deflection of the tooth pair, are acquired by traditional TCA [32] (to the authors’ knowledge, the complete KTE curve of a tooth pair cannot be obtained by the FEM-based TCA). The appropriateness of the selected torque for FEM-based TCA can be validated if the obtained KTE is approximately consistent with KTEs acquired by traditional TCA. It can be seen from Figure 1 that, in contact zone 1, the KTE (absolute value) of tooth pair 1 is less than that of tooth pair 2, which means in this region tooth pair 1 is in contact, whereas tooth pair 2 is not in contact in the no-loading condition. It is because the mating tooth surfaces of tooth pair 2 are separated from each other by the profile separation which can be measured by the difference in KTEs between tooth pairs 1 and 2. Likewise, for contact zones 2 and 3, only tooth pairs 2 and 3 are in contact, respectively, in the unloaded condition. That is why the KTEs of the gear pair in contact zones 1, 2, and 3 correspond to the KTEs of tooth pairs 1, 2, and 3, respectively.The total deformation of a gear pair is usually considered the difference between the KTE and the LTE of the gear pair. This difference, however, is not always the mesh deflection of a tooth pair caused by the load it bears. The schematic diagram for illustrating the deflection of the meshing tooth pair is shown in Figure10. As previously mentioned, at any rolling position in contact zone 1 (shown in Figure 1), only tooth pair 1 is in contact in the no-loading condition and the mating tooth surfaces of tooth pair 2 are separated by the profile separation, as shown in Figure 10. With the augment of torque load, the LTE of the gear pair increases gradually because of the mesh deflection of tooth pair 1. If the LTE of the gear pair at any rolling position in contact zone 1 is greater than the KTE of tooth pair 2, which means the profile separation of tooth pair 2 at this position has been closed because of tooth rotation, tooth pair 2 will be in contact and share the whole load with tooth pair 1. Therefore, as illustrated in Figure 10, part of the rotation of tooth pair 2 compensates the profile separation between the mating flanks and the remaining part is the result of mesh deflection. At any rolling position in contact zone 2, as tooth pair 2 is in contact initially, the tooth rotation of tooth pair 2 is wholly attributed to its deflection. In contact zone 3, however, the profile separation of tooth pair 2 must be closed before contact, and hence, the remaining part of the tooth rotation which excludes the angular profile separation is the rotational deformation of tooth pair 2. Based on above analysis, it can be speculated that, at the torque load of 300 Nm, tooth pairs 1 and 2 are in contact in region AB, only tooth pair 2 is in contact in region BC, and tooth pairs 2 and 3 are in contact in region CD. This speculation is proved to be true by the time histories of mesh force on each engaged tooth pair shown in Figure 9.Figure 10
Deflection of the meshing tooth pair at the meshing position in contact zone 1.Above discussions elaborate the meaning of equation (10). Namely, the tooth rotations Δθ for all engaged tooth pairs are the same and equal the rotation of the pinion and gear blanks which is measured by the difference between the angular KTE and LTE of the gear pair, while the mesh deflection along the line of action δi for the tooth pair i due to the load it shares should be calculated by subtracting the profile separation ΔSi from the projection of the tooth rotation Δθ along the line of action. Actually, the mesh deflection for a tooth pair equals the difference between the KTE of the tooth pair and the LTE of the gear pair, as shown in Figure 1. Therefore, in the MMTP, the KTE of a tooth pair is applied instead of the KTE of a gear pair in the estimation of mesh stiffness and the identification of the contact for each pair of teeth. It is the main difference from the multipoint coupling mesh model.Figure11 shows the mesh deflections along the line of action for the gear pair and each tooth pair. The tooth engagement cycle is divided into regions where two meshing tooth pairs are in contact and one meshing tooth pair is in contact (denoted by DTC and STC, respectively, in Figures 11–13). In the double-tooth contact zone, although the rotational deformation of the gear pair equals the rotational deformation of the tooth pair which is initially in contact in the unloaded condition, their translational deflections along the line of action are somewhat different. This is because in the double-tooth contact zone, the effective mesh point of the gear pair is different from that of the tooth pair, as shown in Figure 4, and thus, the corresponding lines of action of the net mesh force acting on the effective mesh point are different as well.Figure 11
Mesh deflection along the line of action of the gear pair and the tooth pair over one tooth engagement cycle at 300 Nm applied torque.Figure 12
The ASMS of the gear pair and the tooth pair over one tooth engagement cycle at 300 Nm applied torque.Figure 13
The LTMS of the gear pair and the tooth pair over one tooth engagement cycle at 300 Nm applied torque.Figure12 compares the ASMS calculations of the gear pair and the tooth pair. It can be seen that, in the double-tooth contact zone, the ASMS of the gear pair experiences a gradual increase and then a steady decrease. The variation significantly differs from that of unmodified spur gears whose mesh stiffness presents a sharp increase and decrease at the instants when the approaching tooth pair comes into contact and the receding tooth pair leaves the mesh, respectively [21]. Apart from that, with the contact ratio increase, the variation of the mesh stiffness of the gear pair tends to become smaller. This conclusion can be proved by the comparison of the ASMS calculations between this study and studies of Peng [13] and Wang et al. [20].An important fact shown in Figure12 is that, in the double-tooth contact zone, the ASMS of a gear pair does not equal the sum of the ASMS of each meshing tooth pair, if the KTE of a tooth pair is applied to estimate the mesh stiffness of a tooth pair, except for the position where the two meshing tooth pairs are in contact in the unloaded condition, such as point A in the graph. The reason is explained in detail in Appendix A. In many previous studies on mesh stiffness calculation of spur gears [33–36], the total mesh stiffness of the gear pair in the multiple-tooth engagement region is calculated by summation of the mesh stiffness of all tooth pairs in contact. This estimation method is probably inapplicable for spiral bevel or hypoid gears in terms of the ASMS evaluation, whereas it is proved to be feasible for the LTMS calculation.For the calculation of the LTMS expressed in equation (12), LTCA has to be performed at the torque loads of T, T+ΔT, and T−ΔT to obtain the corresponding LTEs εT, εT+ΔT, and εT−ΔT. Since the mesh parameters are all load-dependent, ΔT should be properly chosen to ensure that the effective mesh points and the directions of the line of action are almost the same for the three torque loads at each rolling position. However, it does not mean that the smaller the ΔT is, the better the result will be since the effect of numerical calculation errors on the calculation results will become more pronounced as ΔT decreases. In the present study, ΔT=10Nm can weigh the two factors mentioned above. Figure 13 compares the LTMS calculations of the gear pair and the tooth pair. It can be seen that, in the double-tooth contact zone, the LTMS of the gear pair approximately equals the sum of the LTMS of simultaneously contacting tooth pairs. The reason is mathematically clarified in Appendix A. Comparing calculations in Figures 12 and 13, we can observe that, for both the mesh stiffness of a gear pair and a tooth pair, the LTMS is larger than the ASMS over the entire tooth engagement cycle. It is consistent with the theoretical prediction from the concave force-deflection curve shown in Figure 8. It has been known that the nonlinear contact deformation due to the growth of the contact area as applied load increases results in the concave shape of the force-deflection curve. Note that the difference between the computed results of the ASMS and LTMS is the result of their definitions and calculation approaches, and we cannot make a conclusion about which prediction is more accurate. Both the ASMS and LTMS can be used to calculate DMF, but the corresponding calculation formulas are entirely different.


## 4. DMF Calculation

### 4.1. Calculation Approach Using ASMS

A common representation of DMFFmd utilizing the ASMS is [13, 14, 18, 20](13)Fmd=kmaεd−ε0+cmε˙d−ε˙0,εd−ε0>0,0,−b≤εd−ε0≤0,kmaεd−ε0+b+cmε˙d−ε˙0,εd−ε0<−b,where kma and cm are the ASMS and proportional mesh damping and εd and ε0 are the translation form of the dynamic transmission error and KTE along the line of action. The gear backlash b is expressed on the coast side here. It is noted that here the translational KTE ε0 should be determined by ε0=Np/NgΔθp0−Δθg0λ which is different from the traditional definition of the KTE (equation (8)).The DMF is composed of the elastic contact forceFm,a=kmaεd−ε0 and the viscous damping force Fc=cmε˙d−ε˙0. The elastic contact force Fm,a is illustrated by the blue dashed line in Figure 8. It can be observed that as long as the dynamic mesh deflection δd is not equal to the static mesh deflection δT caused by the nominal load, the calculated contact force will differ from the actual one Fm. The deviation can be acceptable if the dynamic mesh deflection fluctuates around the nominal static mesh deflection within a small range. However, as for moderate and large amplitude vibrations where a substantial difference exists between δd and δT, the validity of Fm,a may be questioned.


### 4.2. Calculation Approach Using LTMS

As shown in Figure8, for small mesh deflections relative to the nominal static deflection δT, the elastic contact force Fm,l (demonstrated by the red dash-dotted line) evaluated by using the LTMS is closer to the true value Fm than the contact force Fm,a. Therefore, in this case, it could be more appropriate to apply the LTMS to calculate DMF. The DMF Fmd considering backlash nonlinearity can be defined as(14)Fmd=Fm0+kmlεd−εl+cmε˙d−ε˙l,Fr>0,0,Fr≤0andεd−ε0≥−b,kmaεd−ε0+b+cmε˙d−ε˙0,Fr≤0andεd−ε0<−b,where kml is the LTMS and εl is the translational LTE along the line of action. The nominal contact force Fm0 can be acquired from LTCA. The reference force Fr which is used to determine the operational condition is given by(15)Fr=Fm0+kmlεd−εl.Since the calculation approach using the LTMS is only suitable for small mesh deflections about the nominal static deflection, this approach is only applied to calculate the DMF for the drive side, while the DMF for the coast side is still estimated by the approach using the ASMS. It has been found that mesh properties for the coast side have less effect on dynamic responses than those of the drive side, and their influences appear only in double-sided impact regions for light-loaded cases [12]. Therefore, if the evaluation of DMF for the drive side is of primary interest, it is reasonable to ignore the mesh characteristic asymmetry. In this study, the mesh parameters of the drive side are used to approximately evaluate the DMF of the coast side for simplicity.It is noteworthy that the representations of the elastic contact force in equations (13) and (14) are fundamentally different, and hence, it is incorrect to apply the ASMS in equation (14) and apply the LTMS in equation (13). For example, the use of the LTMS in equation (13) will give a false result Fm,e (demonstrated by the green dash-dotted line in Figure 8) which differs from the actual one Fm dramatically.


### 4.3. Calculation Approach Using Mesh Force Interpolation Function

As illustrated in Figure8, both the calculation approach using the ASMS and the one applying the LTMS will lose accuracy when the dynamic mesh deflection substantially differs from the nominal static mesh deflection because at any mesh position, the elastic contact is modeled as a linear spring, and therefore, the nonlinear force-deflection relationship is ignored. The effect of such simplification will be examined later by comparing the DMF calculated by the two approaches with the DMF interpolated from the force-deflection function. Since the mesh force also varies with the contact position of meshing tooth pairs, the interpolation function indexed by the nominal mesh position and the mesh deflection at this position is required. This method is also applied by Peng [13] and recently by Dai et al. [37] to the prediction of the total DMF of hypoid gear pairs and spur gear pairs, respectively. Here, the interpolation function for the mesh force of individual tooth pairs at each mesh position over one tooth engagement cycle is developed through a series of quasistatic LTCAs with the applied torque load ranging from 10 Nm to 3000 Nm. This range far exceeds the nominal load (50 Nm and 300 Nm in the present work) that the gear bears because we need to obtain the relationship between the mesh force and the mesh deflection at the beginning and end of the tooth meshing cycle under the nominal load. Looking back at Figure 6(a), we can observe that the mesh force at these two positions will be zero if the applied load is below the nominal load, or be very small if the applied load is just over the nominal load. Therefore, in order to obtain sufficient data for the development of the force-deflection relationship at these two positions, the range of the applied torque load should be wide enough. The force-deflection function at each mesh position over one tooth engagement cycle is graphically shown in Figure 14. Note that even when the load applied for LTCA is 3000 Nm, the data for the mesh deformation exceeding 0.004 mm at the position of tooth engagement and disengagement are not available. Likewise, because of load sharing, the data points marked in Figure 14 cannot be obtained by LTCA with the maximum applied load of 3000 Nm. These marked data are estimated by the LS-SVMlab toolbox [31] based on the data obtained from LTCA. Although the predicted data may not be very precise, they have a marginal effect on the steady-state dynamic response since the dynamic mesh deformation of a tooth pair in the region of tooth engagement and disengagement is also small, and thereby, the precise data from LTCA are used.Figure 14
Relationship between the tooth contact force and its corresponding deformation along the line of action at each mesh position.The interpolated DMFFmd is expressed as follows:(16)Fmd=fdεd−ε0,θt+cmε˙d−ε˙0,εd−ε0>0,0,−b≤εd−ε0≤0,fcεd−ε0+b,θt+cmε˙d−ε˙0,εd−ε0<−b,where fd and fc represent the interpolation functions for the drive and coast sides, which are the same in this work; εd−ε0 and εd−ε0+b express the dynamic mesh deformation along the line of action for the drive and coast sides, respectively; and θt is the mesh position of a tooth pair.


## 4.1. Calculation Approach Using ASMS

A common representation of DMFFmd utilizing the ASMS is [13, 14, 18, 20](13)Fmd=kmaεd−ε0+cmε˙d−ε˙0,εd−ε0>0,0,−b≤εd−ε0≤0,kmaεd−ε0+b+cmε˙d−ε˙0,εd−ε0<−b,where kma and cm are the ASMS and proportional mesh damping and εd and ε0 are the translation form of the dynamic transmission error and KTE along the line of action. The gear backlash b is expressed on the coast side here. It is noted that here the translational KTE ε0 should be determined by ε0=Np/NgΔθp0−Δθg0λ which is different from the traditional definition of the KTE (equation (8)).The DMF is composed of the elastic contact forceFm,a=kmaεd−ε0 and the viscous damping force Fc=cmε˙d−ε˙0. The elastic contact force Fm,a is illustrated by the blue dashed line in Figure 8. It can be observed that as long as the dynamic mesh deflection δd is not equal to the static mesh deflection δT caused by the nominal load, the calculated contact force will differ from the actual one Fm. The deviation can be acceptable if the dynamic mesh deflection fluctuates around the nominal static mesh deflection within a small range. However, as for moderate and large amplitude vibrations where a substantial difference exists between δd and δT, the validity of Fm,a may be questioned.


## 4.2. Calculation Approach Using LTMS

As shown in Figure8, for small mesh deflections relative to the nominal static deflection δT, the elastic contact force Fm,l (demonstrated by the red dash-dotted line) evaluated by using the LTMS is closer to the true value Fm than the contact force Fm,a. Therefore, in this case, it could be more appropriate to apply the LTMS to calculate DMF. The DMF Fmd considering backlash nonlinearity can be defined as(14)Fmd=Fm0+kmlεd−εl+cmε˙d−ε˙l,Fr>0,0,Fr≤0andεd−ε0≥−b,kmaεd−ε0+b+cmε˙d−ε˙0,Fr≤0andεd−ε0<−b,where kml is the LTMS and εl is the translational LTE along the line of action. The nominal contact force Fm0 can be acquired from LTCA. The reference force Fr which is used to determine the operational condition is given by(15)Fr=Fm0+kmlεd−εl.Since the calculation approach using the LTMS is only suitable for small mesh deflections about the nominal static deflection, this approach is only applied to calculate the DMF for the drive side, while the DMF for the coast side is still estimated by the approach using the ASMS. It has been found that mesh properties for the coast side have less effect on dynamic responses than those of the drive side, and their influences appear only in double-sided impact regions for light-loaded cases [12]. Therefore, if the evaluation of DMF for the drive side is of primary interest, it is reasonable to ignore the mesh characteristic asymmetry. In this study, the mesh parameters of the drive side are used to approximately evaluate the DMF of the coast side for simplicity.It is noteworthy that the representations of the elastic contact force in equations (13) and (14) are fundamentally different, and hence, it is incorrect to apply the ASMS in equation (14) and apply the LTMS in equation (13). For example, the use of the LTMS in equation (13) will give a false result Fm,e (demonstrated by the green dash-dotted line in Figure 8) which differs from the actual one Fm dramatically.


## 4.3. Calculation Approach Using Mesh Force Interpolation Function

As illustrated in Figure8, both the calculation approach using the ASMS and the one applying the LTMS will lose accuracy when the dynamic mesh deflection substantially differs from the nominal static mesh deflection because at any mesh position, the elastic contact is modeled as a linear spring, and therefore, the nonlinear force-deflection relationship is ignored. The effect of such simplification will be examined later by comparing the DMF calculated by the two approaches with the DMF interpolated from the force-deflection function. Since the mesh force also varies with the contact position of meshing tooth pairs, the interpolation function indexed by the nominal mesh position and the mesh deflection at this position is required. This method is also applied by Peng [13] and recently by Dai et al. [37] to the prediction of the total DMF of hypoid gear pairs and spur gear pairs, respectively. Here, the interpolation function for the mesh force of individual tooth pairs at each mesh position over one tooth engagement cycle is developed through a series of quasistatic LTCAs with the applied torque load ranging from 10 Nm to 3000 Nm. This range far exceeds the nominal load (50 Nm and 300 Nm in the present work) that the gear bears because we need to obtain the relationship between the mesh force and the mesh deflection at the beginning and end of the tooth meshing cycle under the nominal load. Looking back at Figure 6(a), we can observe that the mesh force at these two positions will be zero if the applied load is below the nominal load, or be very small if the applied load is just over the nominal load. Therefore, in order to obtain sufficient data for the development of the force-deflection relationship at these two positions, the range of the applied torque load should be wide enough. The force-deflection function at each mesh position over one tooth engagement cycle is graphically shown in Figure 14. Note that even when the load applied for LTCA is 3000 Nm, the data for the mesh deformation exceeding 0.004 mm at the position of tooth engagement and disengagement are not available. Likewise, because of load sharing, the data points marked in Figure 14 cannot be obtained by LTCA with the maximum applied load of 3000 Nm. These marked data are estimated by the LS-SVMlab toolbox [31] based on the data obtained from LTCA. Although the predicted data may not be very precise, they have a marginal effect on the steady-state dynamic response since the dynamic mesh deformation of a tooth pair in the region of tooth engagement and disengagement is also small, and thereby, the precise data from LTCA are used.Figure 14
Relationship between the tooth contact force and its corresponding deformation along the line of action at each mesh position.The interpolated DMFFmd is expressed as follows:(16)Fmd=fdεd−ε0,θt+cmε˙d−ε˙0,εd−ε0>0,0,−b≤εd−ε0≤0,fcεd−ε0+b,θt+cmε˙d−ε˙0,εd−ε0<−b,where fd and fc represent the interpolation functions for the drive and coast sides, which are the same in this work; εd−ε0 and εd−ε0+b express the dynamic mesh deformation along the line of action for the drive and coast sides, respectively; and θt is the mesh position of a tooth pair.


## 5. Gear Pair Dynamic Model

### 5.1. A 12-DOF Lumped-Parameter Model

The gear pair model is illustrated in Figure15(a). Each gear with its supporting shaft is modeled as a single rigid body. The gear-shaft bodies are each supported by two compliant rolling element bearings placed at arbitrary axial locations. In the current case, the pinion is overhung supported, while the gear is simply supported.Figure 15
(a) Gear pair model and (b) coordinate systems of each gear-shaft body.
(a)
(b)As shown in Figure15(b), the position of each gear-shaft body in the space is determined by its centroid position Cj and the attitude of the body relative to the centroid, which are defined by the translational coordinates xj,yj,zj and the Cardan angles αj,βj,γj, respectively. The rotation process of the body about the centroid described by the Cardan angles is as follows: the body firstly rotates about the xj1-axis with the angle αj from the original position indicated by the coordinate system Cjxj1yj1zj1 to the position represented by the coordinate system Cjxj2yj2zj2, and then it rotates about the yj2-axis with the angle βj reaching the location shown by the coordinate system Cjxj3yj3zj3. The last rotation is about the zj3-axis with the angle γj, and the body reaches the final instantaneous position represented by the coordinate system Cjxj4yj4zj4. Generally, αj and βj are small angles, while γj is the large rolling angle of the gear. In all the coordinates mentioned above and the following equations, the subscript j=p,g indicates that the quantities belong to the pinion or gear component. The derivation of the equations of motion is demonstrated in Appendix B. Here, the detailed equations of motion of the gear pair with 12 DOFs are given directly as follows:(17a)mpx¨p=−xp∑l=A,Bkpxtl−βp∑l=A,BkpxtlZpl−x˙p∑l=A,Bcpxtl−β˙p∑l=A,BcpxtlZpl+∑i=1Nnpx,iFm,i,(17b)mpy¨p=−yp∑l=A,Bkpxtl+αp∑l=A,BkpxtlZpl−y˙p∑l=A,Bcpxtl+α˙p∑l=A,BcpxtlZpl+∑i=1Nnpy,iFm,i,(17c)mpz¨p=−zp∑l=A,Bkpztl−z˙p∑l=A,Bcpztl+∑i=1Nnpz,iFm,i,(17d)Jpxα¨p+Jpzβ˙pγ˙p+Jpzβpγ¨p=−θpx∑l=A,Bkpxrl+∑l=A,BkpxtlZpl−αpZpl+yp−θ˙px∑l=A,Bcpxrl+∑l=A,BcpxtlZpl−α˙pZpl+y˙p+∑i=1Nλpx,iFm,i,(17e)Jpxβ¨p−Jpzα˙pγ˙p−Jpzαpγ¨p=−θpy∑l=A,Bkpxrl−∑l=A,BkpxtlZplβpZpl+xp−θ˙py∑l=A,Bcpxrl−∑l=A,BcpxtlZplβ˙pZpl+x˙p+∑i=1Nλpy,iFm,i,(17f)Jpzγ¨p=Tp+∑i=1Nλpz,iFm,i,(17g)mgx¨g=−xg∑l=A,Bkgxtl−βg∑l=A,BkgxtlZgl−x˙g∑l=A,Bcgxtl−β˙g∑l=A,BcgxtlZgl+∑i=1Nngx,iFm,i,(17h)mgy¨g=−yg∑l=A,Bkgxtl+αg∑l=A,BkgxtlZgl−y˙g∑l=A,Bcgxtl+α˙g∑l=A,BcgxtlZgl+∑i=1Nngy,iFm,i,(17i)mgz¨g=−zg∑l=A,Bkgztl−z˙g∑l=A,Bcgztl+∑i=1Nngz,iFm,i,(17j)Jgxα¨g+Jgzβ˙gγ˙g+Jgzβgγ¨g=−θgx∑l=A,Bkgxrl+∑l=A,BkgxtlZgl−αgZgl+yg−θ˙gx∑l=A,Bcgxrl+∑l=A,BcgxtlZgl−α˙gZgl+y˙g+∑i=1Nλgx,iFm,i,(17k)Jgxβ¨g−Jgzα˙gγ˙g−Jgzαgγ¨g=−θgy∑l=A,Bkgxrl−∑l=A,BkgxtlZglβgZgl+xg−θ˙gy∑l=A,Bcgxrl−∑l=A,BcgxtlZglβ˙gZgl+x˙g+∑i=1Nλgy,iFm,i,(17l)Jgzγ¨g=−Tg+∑i=1Nλgz,iFm,i,where N is the maximum number of tooth pairs that simultaneously participate in the meshing; Tp and Tg are the driving torque acting on the pinion body and the torque load applied to the gear body, respectively; Fm,i is the DMF of each tooth pair and can be calculated by equations (13), (14), or (16) according to the needs; and nj,injx,i,njy,i,njz,i is the contact force directional cosine vector of each meshing tooth pair i. The calculation of the directional rotation radius λjx,i, λjy,i, and λjz,i is given by(18)λjx,i,λjy,i,λjz,i=njz,iyjm,i−njy,izjm,i,njx,izjm,i−njz,ixjm,i,njy,ixjm,i−njx,iyjm,i.The evaluation of sliding friction is beyond the scope of this article, and thereby, the current dynamic model is formulated without consideration of sliding friction. However, it should be noted that friction is a significant external source of excitation influencing the dynamic behavior of gearing systems. Explanations for other symbols used above are given in AppendixB and nomenclature. The equations of motion given above are expressed based on the MMTP. If the MMGP is applied, the DMF and mesh parameters of the individual tooth pair should be replaced by those of the gear pair. The following description is also based on the MMTP. The related information about the MMGP is introduced in the study by Peng [13].Geometric modeling of spiral bevel gears, TCA, and LTCA are performed in the global fixed reference frameOxyz (as shown in Figure 2), so it is necessary to transform the mesh point position rm,i and the line of action vector Li to the local inertial reference frames of the pinion and gear Ojxj0yj0zj0 in which the dynamic analysis is performed. The relationship between rpm,ixpm,i,ypm,i,zpm,i and rm,ixm,i,ym,i,zm,i is(19)xpm,iypm,izpm,i1=MOpO′rm,i=00−100100100−OpO¯0001xm,iym,izm,i1=−zm,iym,ixm,i−OpO¯1.The directional cosine vector of the mesh force acting on the pinionnp,inpx,i,npy,i,npz,i is(20)npx,inpy,inpz,i=MOpO−Li=00−1010100−nx,i−ny,i−nz,i=nz,i−ny,i−nx,i.The relationship betweenrgm,ixgm,i,ygm,i,zgm,i and rm,ixm,i,ym,i,zm,i is(21)xgm,iygm,izgm,i1=MOgO′rm,i=00−1010000−10−OgO¯0001xm,iym,izm,i1=−zm,ixm,i−ym,i−OgO¯1.The directional cosine vector of the mesh force acting on the gearng,ingx,i,ngy,i,ngz,i is(22)ngx,ingy,ingz,i=MOgOLi=00−11000−10nx,iny,inz,i=−nz,inx,i−ny,i.It is noted that the direction of the mesh force acting on the pinion is opposite to the direction of the line of action vectorLi.


### 5.2. Dynamic Transmission Error Calculation

The translational dynamic transmission error vector is defined as the difference between the effective mesh point position vectors of the pinion and gear. The dynamic transmission error along the line of actionεd can be obtained by projecting the translational dynamic transmission error vector along the line of action vector.For the MMGP, the time-dependent mesh parameters are often treated as the function of the position on the mesh cycle (PMC) [29], while for the MMTP, the mesh parameters of each tooth pair should be treated as the function of the position on the one tooth engagement cycle (PTEC). For current formulation, all the teeth are assumed to be equally spaced and have the same surface topology. The PMC and PTEC are calculated according to the following equations:(23)PMC=θpz+θpz,0−floorθpz+θpz,0θpz,pθpz,p,(24)PTEC=θpz+θt,0−floorθpz+θt,0Nθpz,pNθpz,p,where θpz is the pinion roll angle which approximately equals γp as shown in equation (B.11); θpz,0 is the initial angular position of the pinion which can be arbitrarily selected; and θpz,p is the angular pitch. If equation (24) is applied to calculate the PTEC of the reference tooth pair, θt,0 is the initial angular position of the reference tooth pair which can also be chosen discretionarily, and N is the maximum number of tooth pairs that participate in the meshing simultaneously. The PTEC of other simultaneously engaged tooth pairs can be deduced by the PTEC of the reference tooth pair and the number of angular pitches between them. In both equations, the “floor” function rounds a number to the nearest smaller integer.On the pinion side, it is convenient to express the instantaneous position vector of the effective mesh point in the coordinate systemCpxp3yp3zp3:(25)rpm,i3=xpm,iPTECiip3+ypm,iPTECijp3+zpm,iPTECikp3.To consider the effect of large rotational anglesθpzorγp and θgzorγg on the dynamic transmission error, an imaginary coordinate system CgxgIygIzgI considering the instantaneous angular TE about the z-axis on the gear side θgz−θpz/R is established. The effective mesh point position vector of the gear expressed in CgxgIygIzgI is(26)rgm,iI=xgm,iPTECiigI+ygm,iPTECijgI+zgm,iPTECikgI.Transforming the coordinates to the global fixed reference frameOxyz, we obtain(27)rpm,iO=MOOp′MOpCp′rpm,i3,rgm,iO=MOOg′MOgCg′rgm,iI,where the transformation matrices for Cpxp3yp3zp3 to Oxyz are(28)MOOp′=001OpO¯0100−10000001,MOpCp′=10βpxp01−αpyp−βpαp1zp0001.The transformation matrices forCgxgIygIzgI to Oxyz are(29)MOOg′=010000−1−OgO¯−10000001,MOgCg′=1−θgz−θpz/Rβgxgθgz−θpz/R1−αgyg−βgαg1zg0001.The dynamic transmission error along the line of action for each pair of teeth in contact is calculated by(30)εd,i=rpm,iO−rgm,iO⋅Li.


## 5.1. A 12-DOF Lumped-Parameter Model

The gear pair model is illustrated in Figure15(a). Each gear with its supporting shaft is modeled as a single rigid body. The gear-shaft bodies are each supported by two compliant rolling element bearings placed at arbitrary axial locations. In the current case, the pinion is overhung supported, while the gear is simply supported.Figure 15
(a) Gear pair model and (b) coordinate systems of each gear-shaft body.
(a)
(b)As shown in Figure15(b), the position of each gear-shaft body in the space is determined by its centroid position Cj and the attitude of the body relative to the centroid, which are defined by the translational coordinates xj,yj,zj and the Cardan angles αj,βj,γj, respectively. The rotation process of the body about the centroid described by the Cardan angles is as follows: the body firstly rotates about the xj1-axis with the angle αj from the original position indicated by the coordinate system Cjxj1yj1zj1 to the position represented by the coordinate system Cjxj2yj2zj2, and then it rotates about the yj2-axis with the angle βj reaching the location shown by the coordinate system Cjxj3yj3zj3. The last rotation is about the zj3-axis with the angle γj, and the body reaches the final instantaneous position represented by the coordinate system Cjxj4yj4zj4. Generally, αj and βj are small angles, while γj is the large rolling angle of the gear. In all the coordinates mentioned above and the following equations, the subscript j=p,g indicates that the quantities belong to the pinion or gear component. The derivation of the equations of motion is demonstrated in Appendix B. Here, the detailed equations of motion of the gear pair with 12 DOFs are given directly as follows:(17a)mpx¨p=−xp∑l=A,Bkpxtl−βp∑l=A,BkpxtlZpl−x˙p∑l=A,Bcpxtl−β˙p∑l=A,BcpxtlZpl+∑i=1Nnpx,iFm,i,(17b)mpy¨p=−yp∑l=A,Bkpxtl+αp∑l=A,BkpxtlZpl−y˙p∑l=A,Bcpxtl+α˙p∑l=A,BcpxtlZpl+∑i=1Nnpy,iFm,i,(17c)mpz¨p=−zp∑l=A,Bkpztl−z˙p∑l=A,Bcpztl+∑i=1Nnpz,iFm,i,(17d)Jpxα¨p+Jpzβ˙pγ˙p+Jpzβpγ¨p=−θpx∑l=A,Bkpxrl+∑l=A,BkpxtlZpl−αpZpl+yp−θ˙px∑l=A,Bcpxrl+∑l=A,BcpxtlZpl−α˙pZpl+y˙p+∑i=1Nλpx,iFm,i,(17e)Jpxβ¨p−Jpzα˙pγ˙p−Jpzαpγ¨p=−θpy∑l=A,Bkpxrl−∑l=A,BkpxtlZplβpZpl+xp−θ˙py∑l=A,Bcpxrl−∑l=A,BcpxtlZplβ˙pZpl+x˙p+∑i=1Nλpy,iFm,i,(17f)Jpzγ¨p=Tp+∑i=1Nλpz,iFm,i,(17g)mgx¨g=−xg∑l=A,Bkgxtl−βg∑l=A,BkgxtlZgl−x˙g∑l=A,Bcgxtl−β˙g∑l=A,BcgxtlZgl+∑i=1Nngx,iFm,i,(17h)mgy¨g=−yg∑l=A,Bkgxtl+αg∑l=A,BkgxtlZgl−y˙g∑l=A,Bcgxtl+α˙g∑l=A,BcgxtlZgl+∑i=1Nngy,iFm,i,(17i)mgz¨g=−zg∑l=A,Bkgztl−z˙g∑l=A,Bcgztl+∑i=1Nngz,iFm,i,(17j)Jgxα¨g+Jgzβ˙gγ˙g+Jgzβgγ¨g=−θgx∑l=A,Bkgxrl+∑l=A,BkgxtlZgl−αgZgl+yg−θ˙gx∑l=A,Bcgxrl+∑l=A,BcgxtlZgl−α˙gZgl+y˙g+∑i=1Nλgx,iFm,i,(17k)Jgxβ¨g−Jgzα˙gγ˙g−Jgzαgγ¨g=−θgy∑l=A,Bkgxrl−∑l=A,BkgxtlZglβgZgl+xg−θ˙gy∑l=A,Bcgxrl−∑l=A,BcgxtlZglβ˙gZgl+x˙g+∑i=1Nλgy,iFm,i,(17l)Jgzγ¨g=−Tg+∑i=1Nλgz,iFm,i,where N is the maximum number of tooth pairs that simultaneously participate in the meshing; Tp and Tg are the driving torque acting on the pinion body and the torque load applied to the gear body, respectively; Fm,i is the DMF of each tooth pair and can be calculated by equations (13), (14), or (16) according to the needs; and nj,injx,i,njy,i,njz,i is the contact force directional cosine vector of each meshing tooth pair i. The calculation of the directional rotation radius λjx,i, λjy,i, and λjz,i is given by(18)λjx,i,λjy,i,λjz,i=njz,iyjm,i−njy,izjm,i,njx,izjm,i−njz,ixjm,i,njy,ixjm,i−njx,iyjm,i.The evaluation of sliding friction is beyond the scope of this article, and thereby, the current dynamic model is formulated without consideration of sliding friction. However, it should be noted that friction is a significant external source of excitation influencing the dynamic behavior of gearing systems. Explanations for other symbols used above are given in AppendixB and nomenclature. The equations of motion given above are expressed based on the MMTP. If the MMGP is applied, the DMF and mesh parameters of the individual tooth pair should be replaced by those of the gear pair. The following description is also based on the MMTP. The related information about the MMGP is introduced in the study by Peng [13].Geometric modeling of spiral bevel gears, TCA, and LTCA are performed in the global fixed reference frameOxyz (as shown in Figure 2), so it is necessary to transform the mesh point position rm,i and the line of action vector Li to the local inertial reference frames of the pinion and gear Ojxj0yj0zj0 in which the dynamic analysis is performed. The relationship between rpm,ixpm,i,ypm,i,zpm,i and rm,ixm,i,ym,i,zm,i is(19)xpm,iypm,izpm,i1=MOpO′rm,i=00−100100100−OpO¯0001xm,iym,izm,i1=−zm,iym,ixm,i−OpO¯1.The directional cosine vector of the mesh force acting on the pinionnp,inpx,i,npy,i,npz,i is(20)npx,inpy,inpz,i=MOpO−Li=00−1010100−nx,i−ny,i−nz,i=nz,i−ny,i−nx,i.The relationship betweenrgm,ixgm,i,ygm,i,zgm,i and rm,ixm,i,ym,i,zm,i is(21)xgm,iygm,izgm,i1=MOgO′rm,i=00−1010000−10−OgO¯0001xm,iym,izm,i1=−zm,ixm,i−ym,i−OgO¯1.The directional cosine vector of the mesh force acting on the gearng,ingx,i,ngy,i,ngz,i is(22)ngx,ingy,ingz,i=MOgOLi=00−11000−10nx,iny,inz,i=−nz,inx,i−ny,i.It is noted that the direction of the mesh force acting on the pinion is opposite to the direction of the line of action vectorLi.


## 5.2. Dynamic Transmission Error Calculation

The translational dynamic transmission error vector is defined as the difference between the effective mesh point position vectors of the pinion and gear. The dynamic transmission error along the line of actionεd can be obtained by projecting the translational dynamic transmission error vector along the line of action vector.For the MMGP, the time-dependent mesh parameters are often treated as the function of the position on the mesh cycle (PMC) [29], while for the MMTP, the mesh parameters of each tooth pair should be treated as the function of the position on the one tooth engagement cycle (PTEC). For current formulation, all the teeth are assumed to be equally spaced and have the same surface topology. The PMC and PTEC are calculated according to the following equations:(23)PMC=θpz+θpz,0−floorθpz+θpz,0θpz,pθpz,p,(24)PTEC=θpz+θt,0−floorθpz+θt,0Nθpz,pNθpz,p,where θpz is the pinion roll angle which approximately equals γp as shown in equation (B.11); θpz,0 is the initial angular position of the pinion which can be arbitrarily selected; and θpz,p is the angular pitch. If equation (24) is applied to calculate the PTEC of the reference tooth pair, θt,0 is the initial angular position of the reference tooth pair which can also be chosen discretionarily, and N is the maximum number of tooth pairs that participate in the meshing simultaneously. The PTEC of other simultaneously engaged tooth pairs can be deduced by the PTEC of the reference tooth pair and the number of angular pitches between them. In both equations, the “floor” function rounds a number to the nearest smaller integer.On the pinion side, it is convenient to express the instantaneous position vector of the effective mesh point in the coordinate systemCpxp3yp3zp3:(25)rpm,i3=xpm,iPTECiip3+ypm,iPTECijp3+zpm,iPTECikp3.To consider the effect of large rotational anglesθpzorγp and θgzorγg on the dynamic transmission error, an imaginary coordinate system CgxgIygIzgI considering the instantaneous angular TE about the z-axis on the gear side θgz−θpz/R is established. The effective mesh point position vector of the gear expressed in CgxgIygIzgI is(26)rgm,iI=xgm,iPTECiigI+ygm,iPTECijgI+zgm,iPTECikgI.Transforming the coordinates to the global fixed reference frameOxyz, we obtain(27)rpm,iO=MOOp′MOpCp′rpm,i3,rgm,iO=MOOg′MOgCg′rgm,iI,where the transformation matrices for Cpxp3yp3zp3 to Oxyz are(28)MOOp′=001OpO¯0100−10000001,MOpCp′=10βpxp01−αpyp−βpαp1zp0001.The transformation matrices forCgxgIygIzgI to Oxyz are(29)MOOg′=010000−1−OgO¯−10000001,MOgCg′=1−θgz−θpz/Rβgxgθgz−θpz/R1−αgyg−βgαg1zg0001.The dynamic transmission error along the line of action for each pair of teeth in contact is calculated by(30)εd,i=rpm,iO−rgm,iO⋅Li.


## 6. Numerical Simulation Analysis

The proposed equations of motion are solved by MATLAB’s ODE45 [38] which implements the 4/5th order Runge–Kutta integration routine with the adaptive step size. The second-order differential form of equations (17a)–(17l) must be rewritten in the state space (first-order) form before they can be solved by the solver. Step speed sweep simulations are performed to obtain the dynamic response in a prescribed range of mesh frequency. Numerical integration at each mesh frequency is performed up to the time when the steady-state response is reached. The variable values of the final state at this mesh frequency are applied as the initial conditions for dynamic analysis at the next one. Time-domain DMF response at a specified mesh frequency is computed directly from the solution of numerical integration. Frequency response is obtained by root mean square (RMS) values of the time-domain response over one tooth engagement period at each frequency. Note that time histories of DMF on each pair of teeth over one tooth engagement cycle are not the same at the frequencies where subharmonic or chaotic motion occurs, which results in different RMS values for different tooth pairs. In such cases, the average value is adopted for graphing.The system parameters are listed in Table2. The rated torque load of the gear pair is about 700 Nm. Considering the dynamic characteristics for medium- and heavy-load cases are similar to but different from those for the light-load case [13], we performed dynamic simulations at 50 Nm and 300 Nm torque loads to examine the modeling effects for the light-load and medium-to-heavy-load conditions, respectively.Table 2
System parameters.
Parameters
Pinion
Gear
Massm (kg)
1.82
6.71
Torsional moment of inertiaJz (kg·m2)
1.065E − 3
3.94E − 2
Bending moment of inertiaJx (kg·m2)
7.064E − 4
2E − 2
Centroid position of the gearL (mm)
97.507
−41.718
Bearing A axial positionLA (mm)
120
200
Bearing B axial positionLB (mm)
80
80
Translational bearing stiffnesskxt (N/m)
2.8E8
2.8E8
Axial bearing stiffnesskzt (N/m)
8E7
8E7
Tilting bearing stiffnesskxr (Nm/rad)
6E6
8E6
Mesh damping ratio
0.06
Support component damping ratio
0.02
Gear backlash (mm)
0.15


### 6.1. Effects of Mesh Model

To preclude the influence of mesh stiffness on the calculation results of DMF, the more often used ASMS is adopted for both the MMTP and the MMGP.Figure16 compares the DMF responses with respect to the mesh frequency from the two models. It can be seen that the simulation results from the two models are overall consistent, which partly verifies the validity of the MMTP proposed in this work. Differences exist in the response amplitudes particularly at higher mesh frequency, but the resonance frequencies estimated by the two models are approximately the same. In the lightly loaded condition (50 Nm), jump discontinuities in the vicinity of the resonance frequencies are predicted by both models. The tooth separation and impact occur in the frequency range between the upward and downward jump discontinuities. The jump frequencies near the resonance frequency of about 3000 Hz are different for speed sweep up and down simulations owing to the strong nonlinearity of the system. For the medium-to-heavy-loaded case (300 Nm), the nonlinear time-varying responses from both models are continuous curves, and the results for speed sweep up and down cases are just the same. The disappearance of nonlinear jumps means the load has become large enough to prevent the loss of contact and the backlash nonlinearity is ineffectual. Generally, the increase of load tends to weaken the effect of backlash nonlinearity not only because of the higher mean load which plays the role in keeping tooth surfaces in contact [11] but also because the variations of mesh stiffness, mesh vector, and LTE which all aggravate the degree of backlash nonlinearity are alleviated with the increase of load.Figure 16
RMS values of the DMF acting on the individual tooth pair evaluated at (a) 50 Nm and (b) 300 Nm torque loads.
(a)
(b)The essential differences between the MMTP and the MMGP can be explored further by comparing the time histories of DMF per tooth pair. Analysis results that can describe some typical phenomena of gear meshing are illustrated as follows.Figure17 compares the response histories at 50 Hz mesh frequency and 300 Nm torque load. As can be seen, globally, the DMF predicted by the MMGP waves smoothly around the static MF obtained through LTCA. By contrast, the DMF from the MMTP seems to fluctuate more frequently with dramatic oscillation (due to the tooth engaging-in impact) occurring at the start of tooth engagement which, however, cannot be observed in the results from the MMGP. The overall difference in the predicted DMFs can be attributed to the fact that the MMGP uses the mesh parameters which reflect the gentle overall effect of time-varying mesh characteristics on the entire gear pair, while the more strongly changed mesh parameters of a tooth pair are applied in the MMTP.Figure 17
Time histories of DMF per tooth pair at 50 Hz mesh frequency and 300 Nm torque load.It is known that tooth impact always occurs at the moment of the first contact of each tooth pair. This phenomenon can be explained by the short infinite high acceleration appearing in the acceleration graph (which is the second derivative of the motion or KTE graph) at the changeover point between two adjacent pairs of teeth [39]. The results shown in Figure 17 indicate that the MMTP is capable of describing tooth impact at the instant of engagement and disengagement, but the MMGP is incompetent in this respect. No impact occurs in the estimation from the MMGP because the DMF on each tooth pair in simultaneous contact can only be roughly estimated through multiplying the net mesh force by the quasistatically defined load distribution factor which gradually increases from zero in the engaging-in stage and decreases to zero when the tooth pair leaves the contact.Another interesting phenomenon shown in the graph is that the DMF from the MMTP oscillates above the static mesh force in the engaging-in stage but below the static mesh force in the engaging-out stage. To find the reason, we applied the MMTP to a two-DOF dynamic model which takes into account only the torsional motions of the pinion and gear. It is shown in Figure17 that the predicted DMF fluctuates around the static mesh force from LTCA. This seems to be reasonable because, in LTCA, the pinion and gear are restricted to rotating about their rolling axes as well. The results thereby indicate that the deviation between the DMF and the static mesh force is because of the change of the relative position of the pinion and gear due to the release of additional lateral degrees of freedom. For the MMGP, the mesh force acting on each pair of teeth obtained according to the principle of load distribution waves around the static mesh force. The estimation, however, may deviate from the actuality if the pose of the gear can change in any degree of freedom. The simulation results also indicate that the mesh parameters obtained by the current LTCA method are more suitable for the two-DOF dynamic model.At the light torque load of 50 Nm (as shown in Figure18), the deviation between the DMF and the static mesh force is reduced, and the engaging-in impact is less significant. This is because the change of relative position of the pinion and gear caused by the load is relatively small at the lightly loaded condition. Therefore, we can find that the transmission performance of a spiral bevel gear pair is susceptible to the relative position of the two gears, and a reasonable design of the external support structure plays an essential role in improving the stability of the gear transmission. Tooth flank modification is also necessary to reduce the sensitivity of transmission performance to the inevitable change of relative position caused by the load or installation errors.Figure 18
Time histories of DMF per tooth pair at the mesh frequency of 50 Hz and 50 Nm torque load.It is noteworthy that the prediction from the mesh model with predefined mesh parameters has an unavoidable deviation from the reality as the mesh parameters under the dynamic condition differ from those for the static case. Therefore, the impact force sometimes may be overpredicted because of the inaccurate initial contact position of a tooth pair. The inaccuracy leads to overestimation of the contact deformation at the start of tooth engagement. However, according to the contact conditions of gear teeth discussed before, the tooth pair will get into contact ahead of the quasistatically defined beginning of contact because the large deformation evaluated at the point of the predefined start of contact means the profile separation between the mating flanks of the extended contact area has been closed.Figure19 compares the DMF histories per tooth pair at the mesh frequency of 400 Hz where loss of contact occurs. As displayed, near resonant frequency, DMF varies drastically with the peak value much higher than static mesh force. The two mesh models provide almost identical results with a slight difference near the end of tooth engagement where transient loss of contact only appears in the results from the MMTP. It is worth noting that the differences between the calculations from the MMTP and MMGP usually appear near the beginning and end of tooth contact since in these two regions, the KTE curve of the tooth pair and the gear pair, which are separately applied in the two models, is distinct, as shown in Figure 1.Figure 19
Time histories of DMF per tooth pair at the mesh frequency of 400 Hz and 50 Nm torque load.


### 6.2. Effects of Mesh Stiffness

Here, the MMTP is applied for all analyses. The dynamic responses and time histories of DMF per tooth pair predicted by using the ASMS and LTMS are compared with the results from the mesh force interpolation function, which, in theory, are believed to be more accurate.Figure20 shows the comparison of the DMF responses at 300 Nm torque load. Differences occur in both response amplitudes and resonant frequencies and become noticeable in the range of high mesh frequency. Overall, both the approaches using the ASMS and LTMS tend to underestimate the effective values of DMF, and the predicted resonance frequencies consistently shift to the lower mesh frequency range in comparison with the results from the mesh force interpolation function. Generally, the results from the model applying the LTMS seem to be closer to the interpolated results. This finding can also be observed from the time histories of the individual DMF illustrated in Figure 21.Figure 20
RMS values of the DMF acting on each tooth pair evaluated by the approaches using the ASMS, LTMS, and mesh force interpolation function at 300 Nm toque load.Figure 21
Time histories of DMF per tooth pair estimated by the approaches using the ASMS, LTMS, and mesh force interpolation function at the mesh frequency of 50 Hz and 300 Nm torque load.As shown in Figure21, DMF predicted by the mesh force interpolation function changes most smoothly, and the tooth engaging-in impact is the minimum among the three approaches. It manifests that if the nonlinear relationship between the mesh force and the mesh deflection (i.e., the variation of mesh stiffness with mesh deflection) is considered, the fluctuation of DMF is actually moderate. From the graph, we can also observe that the DMF estimated with the LTMS varies smoothly relative to the DMF predicted by the ASMS, and the tooth engaging-in impact is less powerful. Although the two methods do not give substantially different results, the approach using the LTMS is preferred for vibration analysis of gears around the nominal static deformation. This is not only because the results predicted with the LTMS are more similar to the results from the mesh force interpolation function which precisely describes the force versus deformation behavior but also because the calculation approach using the LTMS shown in equation (14) better represents the relationship between the DMF and the nominal static mesh force when gears vibrate around the nominal static deformation.Similar to the results shown in Figure21, for the example gear pair, in the predefined mesh frequency range (0–6000 Hz), the impact rarely occurs when the tooth pair leaves the contact. It should be noted that whether impact occurs at the instant of tooth engagement and disengagement and the magnitude of the impact largely depends on the transmission error which can be modified through tooth surface modification.


## 6.1. Effects of Mesh Model

To preclude the influence of mesh stiffness on the calculation results of DMF, the more often used ASMS is adopted for both the MMTP and the MMGP.Figure16 compares the DMF responses with respect to the mesh frequency from the two models. It can be seen that the simulation results from the two models are overall consistent, which partly verifies the validity of the MMTP proposed in this work. Differences exist in the response amplitudes particularly at higher mesh frequency, but the resonance frequencies estimated by the two models are approximately the same. In the lightly loaded condition (50 Nm), jump discontinuities in the vicinity of the resonance frequencies are predicted by both models. The tooth separation and impact occur in the frequency range between the upward and downward jump discontinuities. The jump frequencies near the resonance frequency of about 3000 Hz are different for speed sweep up and down simulations owing to the strong nonlinearity of the system. For the medium-to-heavy-loaded case (300 Nm), the nonlinear time-varying responses from both models are continuous curves, and the results for speed sweep up and down cases are just the same. The disappearance of nonlinear jumps means the load has become large enough to prevent the loss of contact and the backlash nonlinearity is ineffectual. Generally, the increase of load tends to weaken the effect of backlash nonlinearity not only because of the higher mean load which plays the role in keeping tooth surfaces in contact [11] but also because the variations of mesh stiffness, mesh vector, and LTE which all aggravate the degree of backlash nonlinearity are alleviated with the increase of load.Figure 16
RMS values of the DMF acting on the individual tooth pair evaluated at (a) 50 Nm and (b) 300 Nm torque loads.
(a)
(b)The essential differences between the MMTP and the MMGP can be explored further by comparing the time histories of DMF per tooth pair. Analysis results that can describe some typical phenomena of gear meshing are illustrated as follows.Figure17 compares the response histories at 50 Hz mesh frequency and 300 Nm torque load. As can be seen, globally, the DMF predicted by the MMGP waves smoothly around the static MF obtained through LTCA. By contrast, the DMF from the MMTP seems to fluctuate more frequently with dramatic oscillation (due to the tooth engaging-in impact) occurring at the start of tooth engagement which, however, cannot be observed in the results from the MMGP. The overall difference in the predicted DMFs can be attributed to the fact that the MMGP uses the mesh parameters which reflect the gentle overall effect of time-varying mesh characteristics on the entire gear pair, while the more strongly changed mesh parameters of a tooth pair are applied in the MMTP.Figure 17
Time histories of DMF per tooth pair at 50 Hz mesh frequency and 300 Nm torque load.It is known that tooth impact always occurs at the moment of the first contact of each tooth pair. This phenomenon can be explained by the short infinite high acceleration appearing in the acceleration graph (which is the second derivative of the motion or KTE graph) at the changeover point between two adjacent pairs of teeth [39]. The results shown in Figure 17 indicate that the MMTP is capable of describing tooth impact at the instant of engagement and disengagement, but the MMGP is incompetent in this respect. No impact occurs in the estimation from the MMGP because the DMF on each tooth pair in simultaneous contact can only be roughly estimated through multiplying the net mesh force by the quasistatically defined load distribution factor which gradually increases from zero in the engaging-in stage and decreases to zero when the tooth pair leaves the contact.Another interesting phenomenon shown in the graph is that the DMF from the MMTP oscillates above the static mesh force in the engaging-in stage but below the static mesh force in the engaging-out stage. To find the reason, we applied the MMTP to a two-DOF dynamic model which takes into account only the torsional motions of the pinion and gear. It is shown in Figure17 that the predicted DMF fluctuates around the static mesh force from LTCA. This seems to be reasonable because, in LTCA, the pinion and gear are restricted to rotating about their rolling axes as well. The results thereby indicate that the deviation between the DMF and the static mesh force is because of the change of the relative position of the pinion and gear due to the release of additional lateral degrees of freedom. For the MMGP, the mesh force acting on each pair of teeth obtained according to the principle of load distribution waves around the static mesh force. The estimation, however, may deviate from the actuality if the pose of the gear can change in any degree of freedom. The simulation results also indicate that the mesh parameters obtained by the current LTCA method are more suitable for the two-DOF dynamic model.At the light torque load of 50 Nm (as shown in Figure18), the deviation between the DMF and the static mesh force is reduced, and the engaging-in impact is less significant. This is because the change of relative position of the pinion and gear caused by the load is relatively small at the lightly loaded condition. Therefore, we can find that the transmission performance of a spiral bevel gear pair is susceptible to the relative position of the two gears, and a reasonable design of the external support structure plays an essential role in improving the stability of the gear transmission. Tooth flank modification is also necessary to reduce the sensitivity of transmission performance to the inevitable change of relative position caused by the load or installation errors.Figure 18
Time histories of DMF per tooth pair at the mesh frequency of 50 Hz and 50 Nm torque load.It is noteworthy that the prediction from the mesh model with predefined mesh parameters has an unavoidable deviation from the reality as the mesh parameters under the dynamic condition differ from those for the static case. Therefore, the impact force sometimes may be overpredicted because of the inaccurate initial contact position of a tooth pair. The inaccuracy leads to overestimation of the contact deformation at the start of tooth engagement. However, according to the contact conditions of gear teeth discussed before, the tooth pair will get into contact ahead of the quasistatically defined beginning of contact because the large deformation evaluated at the point of the predefined start of contact means the profile separation between the mating flanks of the extended contact area has been closed.Figure19 compares the DMF histories per tooth pair at the mesh frequency of 400 Hz where loss of contact occurs. As displayed, near resonant frequency, DMF varies drastically with the peak value much higher than static mesh force. The two mesh models provide almost identical results with a slight difference near the end of tooth engagement where transient loss of contact only appears in the results from the MMTP. It is worth noting that the differences between the calculations from the MMTP and MMGP usually appear near the beginning and end of tooth contact since in these two regions, the KTE curve of the tooth pair and the gear pair, which are separately applied in the two models, is distinct, as shown in Figure 1.Figure 19
Time histories of DMF per tooth pair at the mesh frequency of 400 Hz and 50 Nm torque load.


## 6.2. Effects of Mesh Stiffness

Here, the MMTP is applied for all analyses. The dynamic responses and time histories of DMF per tooth pair predicted by using the ASMS and LTMS are compared with the results from the mesh force interpolation function, which, in theory, are believed to be more accurate.Figure20 shows the comparison of the DMF responses at 300 Nm torque load. Differences occur in both response amplitudes and resonant frequencies and become noticeable in the range of high mesh frequency. Overall, both the approaches using the ASMS and LTMS tend to underestimate the effective values of DMF, and the predicted resonance frequencies consistently shift to the lower mesh frequency range in comparison with the results from the mesh force interpolation function. Generally, the results from the model applying the LTMS seem to be closer to the interpolated results. This finding can also be observed from the time histories of the individual DMF illustrated in Figure 21.Figure 20
RMS values of the DMF acting on each tooth pair evaluated by the approaches using the ASMS, LTMS, and mesh force interpolation function at 300 Nm toque load.Figure 21
Time histories of DMF per tooth pair estimated by the approaches using the ASMS, LTMS, and mesh force interpolation function at the mesh frequency of 50 Hz and 300 Nm torque load.As shown in Figure21, DMF predicted by the mesh force interpolation function changes most smoothly, and the tooth engaging-in impact is the minimum among the three approaches. It manifests that if the nonlinear relationship between the mesh force and the mesh deflection (i.e., the variation of mesh stiffness with mesh deflection) is considered, the fluctuation of DMF is actually moderate. From the graph, we can also observe that the DMF estimated with the LTMS varies smoothly relative to the DMF predicted by the ASMS, and the tooth engaging-in impact is less powerful. Although the two methods do not give substantially different results, the approach using the LTMS is preferred for vibration analysis of gears around the nominal static deformation. This is not only because the results predicted with the LTMS are more similar to the results from the mesh force interpolation function which precisely describes the force versus deformation behavior but also because the calculation approach using the LTMS shown in equation (14) better represents the relationship between the DMF and the nominal static mesh force when gears vibrate around the nominal static deformation.Similar to the results shown in Figure21, for the example gear pair, in the predefined mesh frequency range (0–6000 Hz), the impact rarely occurs when the tooth pair leaves the contact. It should be noted that whether impact occurs at the instant of tooth engagement and disengagement and the magnitude of the impact largely depends on the transmission error which can be modified through tooth surface modification.


## 7. Conclusion

In this paper, two critical factors influencing the calculation of dynamic mesh force (DMF) for spiral bevel gears, i.e., the mesh model and the mesh stiffness representation, have been investigated. An improved mesh model of a tooth pair (MMTP) which better describes the mesh characteristics of each pair of teeth is proposed. As the difference between the kinematic transmission error (KTE) of a tooth pair and the loaded transmission error (LTE) of a gear pair represents the actual deformation of the tooth pair caused by the load it bears, instead of the KTE of a gear pair which is commonly used in the single-point and multipoint coupling mesh models, the KTE of a tooth pair is applied in the MMTP for the calculation of the mesh stiffness of individual tooth pairs, the identification of tooth contact, and the estimation of dynamic mesh deflection of each pair of teeth in contact. The mesh parameters and dynamic responses are compared between the proposed model and the mesh model of a gear pair (MMGP) (i.e., the single-point coupling mesh model). Overall, the two models predict the similar variation trend of dynamic responses with the mesh frequency. Although noticeable differences can be observed in the response amplitudes particularly in the range of high mesh frequency, the resonance frequencies estimated by the two models are almost the same. The time histories of dynamic mesh force (DMF) per tooth pair from the two models are significantly different since the variations of mesh parameters which act as the parametric excitations in the dynamic model are dramatic for the MMTP but moderate for the MMGP. Compared with previous mesh models, the MMTP is superior in simulating the dynamic contact behavior of each pair of teeth especially at the instant of tooth engagement and disengagement. The results from the MMTP also show that the DMF will deviate from the static mesh force if the pose of the gear can change in any degree of freedom. This phenomenon, however, has never been observed in previous studies where the single-point or multipoint coupling mesh model is applied.For the MMGP, the total mesh stiffness of all contacting tooth pairs is applied, while for the MMTP, the stiffness of the individual tooth pair is used. Both stiffnesses can be calculated by the average slope and local slope approaches. The average secant mesh stiffness (ASMS) derived from the average slope approach substantially differs from the local tangent mesh stiffness (LTMS) estimated by the local slope approach owing to their definitions and calculation methods, and we cannot make a conclusion about which calculation is more appropriate. Besides, because of the particularity of the tooth shape, the ASMS for the MMGP cannot be directly computed by the summation of the ASMS of each tooth pair, but the LTMS of a gear pair approximately equals the sum of the LTMS of all tooth pairs in contact. Both the ASMS and LTMS can be used to calculate DMF, but the calculation formulas are rather different. The misuse of mesh stiffness for a given method will lead to errors in the physical sense.DMF calculations from the method applying the ASMS and LTMS are compared with those from the mesh force interpolation function which accurately represents the nonlinear force-deformation relationship. Although the results given by the methods using the LTMS and ASMS are not substantially different, the former approach is recommended to calculate the DMF when simulating the vibration of gears about the nominal static deformation. One of the reasons is that the results predicted with the LTMS are closer to the results from the mesh force interpolation function. Another reason is that the model using the LTMS better represents the relationship between the DMF and the nominal static mesh force when gears vibrate around the nominal static deformation. The third reason is that the assumption, which believes the position of the effective mesh point does not change before and after the tooth deformation, is more reasonable for the calculation of LTMS. This is because the change of the effective mesh point position caused by the mean tooth deflection relative to the undeformed state is noticeable, but the change caused by the small deformation relative to the nominal state is negligible. The approach applying the mesh force interpolation function can provide the most accurate prediction. However, the development of the mesh force interpolation function requires multiple time-consuming loaded tooth contact analyses at different load levels, and the interpolation in solving dynamic problems needs more computational efforts than the other two approaches.



---

*Source: 5614574-2019-10-20.xml*