# Experimental and Theoretical Study of Gear Dynamical Transmission Characteristic Considering Measured Manufacturing Errors

Authors: Fang Guo; Zongde Fang
Journal: Shock and Vibration
 (2018)

---

## Abstract

In the research of gear transmission, the vibration and noise problem has received many concerns all the times. Scholars use tooth modification technique to improve the meshing state of gearings in order to reduce the vibration and noise. However, few of researchers consider the influence of measured manufacturing errors when they do the study of tooth modification. In order to investigate the efficiency of the tooth modification in the actual project, this paper proposes a dynamic model of a helical gear pair including tooth modification and measured manufacturing errors to do a deterministic analysis on the dynamical transmission performance. In this analysis, based on the measured tooth deviation, a real tooth surface (including modification and measured tooth profile error) is fitted by a bicubic B-spline. With the tooth contact analysis (TCA) and loaded tooth contact analysis (LTCA) on the real tooth surface, the loaded transmission error, tooth surface elastic deformation, and load distribution can be determined. Based on the results, the time-varying mesh stiffness and gear mesh impact are computed. Taking the loaded transmission error, measured cumulative pitch error, eccentricity error, time-varying mesh stiffness, and gear mesh impact as the internal excitations, this paper establishes a 12-degree-of-freedom (DOF) dynamic model of a helical gear pair and uses the Fourier series method to solve it. In two situations of low speed and high speed, the gear system dynamic response is analyzed in the time and frequency domains. In addition, an experiment is performed to validate the simulation results. The study shows that the proposed technique is useful and reliable for predicting the dynamic response of a gear system.


---

## Body

## 1. Introduction

In the gear machining process, the design deviation, manufacturing error, and fixing error of the machine cutting tool can cause a change in conjugate conditions and produce a tooth profile error, pitch error, eccentric error, etc. These errors affect the contact characteristics of gears, worsen the stability and bearing capability of gear transmission, and cause enormous noise and vibration. Thus, many scholars pay close attention to the prediction and control of gear vibration and noise in their transmission designs [1–8]. After an extensive study, some researchers have used the tooth profile modification to compensate the elastic deformation of gear teeth to reduce the gear vibration [9–14]. However, there are errors remained during the theoretical profile modification in the machining process. Therefore, the accuracy and efficiency of the tooth profile modification must be checked to avoid bad working conditions of the gear trains.In the related studies, there are some cases which address the gear errors in the study of gear dynamics analysis. Case one: the effect of errors is ignored. Case two: although the errors are considered, the range of every error is estimated in theory, and errors supposedly change according to a particular curve. For example, Hu et al. [15] researched the effect of tooth modification of a high-speed gear system and set the variation of the tooth profile error and pitch error as a sine function. Kubo and Kiyono [16] analyzed the effect of tooth form errors on the helical gear dynamics, and they assume a composite tooth form error with four types of shape. Wei et al. [17] established a coupled nonlinear dynamics model which is used for planetary gear transmission systems in wind turbines and incorporates the effects of the time-varying mesh stiffness, dynamic transmission error, gear mesh impact, and input varying load; thus, the transmission error changes as a sine wave with time. Case three: the errors are shown as random processes. Chen et al. [18] used a stochastic volatility model to simulate the random wind velocity and analyzed the effect on the dynamic characteristics of wind turbine gear transmission system. Wang and Zhang [19] developed a dynamic and stochastic simulation model to analyze the vibration of the gear transmission system; they considered the gear transmission errors as the displacement excitation which was decomposed into harmonic and random components. They used a second-order Markov process with time-variant parameters considering the influence of the rotational speed to simulate a random component. Chen et al. [20] proposed a rotational nonlinear dynamic model of a planetary gear transmission system, which included the time-varying mesh stiffness and synthetic mesh error with random fluctuation; the synthetic mesh error was resolved into a sine signal and white noise. Nevertheless, the most related studies did not consider the low-frequency error (such as eccentricity error, cumulative pitch error, etc.) and hardly showed the real manufacturing errors of gear, so the simulation results might produce an obviously different dynamic response from the experiment.In this paper, the main focus is to investigate the impact of the measured tooth profile error and cumulative pitch error on the dynamic response of gear system with tooth modification. The main objective is to do dynamics analysis of a helical gear pair based on the real tooth surface (including tooth modification and tooth profile error) and the measured cumulative pitch error, so that the transmission performance of the machined gear system can be predicted. The internal excitations, including time-varying mesh stiffness, gear mesh impact, loaded transmission error, measured cumulative pitch error, and eccentric errors, will be considered in the analysis. At the end, the proposed method is compared with the traditional method and verified by a noise and vibration experiment.


## 2. Structure of Real Tooth Surface

A pair of helical gears is used as an example, which parameters are shown in Table1. The tooth surface is finish machined by grinding, and the processing precision is level 6 (adopting the standard of ISO 1328-1:1997). The pinion (drive wheel) has profile and lead modification, and the gear (driven wheel) is a standard involute gear. The theoretical modification curve is shown in Figure 1. Figure 2 is the actual modification curve measured by a Klingelnberg P100 gear measuring device. The gear’s measured surface deviation from the theoretical involute is displayed in Figure 3. The tooth surface deviation of the gear is notably small; thus, it is considered as the theoretical involute in the following calculation, and only the real tooth surface of the pinion must be structured.Table 1
Parameters of pair of helical gears.
Parameters
Pinion
Gear
Number of teeth,Z
19
47
Handedness
Right
Left
Normal modulus,mn (mm)
6
6
Normal pressure angle,α (deg.)
20
20
Helix angle,β (deg.)
–9.911
9.911
Tooth width,b (mm)
79
75
Load torque,Tg (N·m)
830
Input speed, rpm
50.8/2104Figure 1
Theoretical modification curve: (a) profile modification curve along meshing line direction; (b) lead modification curve.
(a)
(b)Figure 2
Measured modification curve of pinion: (a) profile modification curve along meshing line direction; (b) lead modification curve. (In actual processing, only one side of the tooth, which is the right flank, is modified, as shown in Figure 2. Three teeth are measured, namely tooth 1, tooth 8, and tooth 14. This paper selects tooth 8 to structure the real tooth surface.)
(a)
(b)Figure 3
Tooth surface deviation of gear: (a) profile deviation along meshing line direction; (b) lead deviation.
(a)
(b)


### 2.1. Fitting of Surface Deviation

Some scholars have applied the bicubic B-spline method to many studies. Gardner and Gardner [21] used it to model a magnetohydrodynamic duct flow and obtained excellent consistency between the numerical result and the analytic solution even for large values of the applied magnetic field. Na [22] used the method to fit a complex curved surface, and the fitting surface was C2 conforming. This paper also uses a bicubic B-spline to form the deviation surface.Before fitting the surface, some discrete points must be set on the measured modification curve in Figure2. To improve the fitting accuracy, we select more discrete points at both ends of the measured modification curve than in the middle. Finally, 42 points are selected in the meshing line direction and 47 points are selected in the lead direction. Then, the deviations of 42 points in the meshing line direction must be converted into the deviations along the tooth profile direction. According to the meshing relationship shown in Figure 4, the relationship between the tooth depth and the length of the meshing line is as follows:(1)Yprofile=Ymeshing2+rbp2,where Yprofile is the tooth depth, Ymeshing is the length of the meshing line, and rbp is the base radius of the pinion.Figure 4
Relationship between tooth depth and length of meshing line.As shown in Equation (2), the 42 × 47 discrete points are arranged in a matrix.(2)P=P00P01⋯P0nP10P11⋯P1n⋮⋮⋱⋮Pm0Pm1⋯Pmn,m=41,n=46,Based on the 42 × 47 discrete points, we fit the cubic B-spline curves in the directions of lead and tooth depth. The control vertices of the cubic B-spline curve are computed by using the method of reverse vertices, which is presented in Appendix [22]. Then, the cubic B-spline curve equations generated by these vertices are as follows:The lead direction,(3)Li,js=1ss2s3⋅161410−30303−630−13−31Vi,jVi,j+1Vi,j+2Vi,j+3,v=x−xixi+1−xi,x∈xi,xi+1,i=0,1,…,m;j=−1,0,1,…,n−2.The direction of tooth depth,(4)Ii,jv,s=1vv2v3⋅161410−30303−630−13−31Wsi,jWsi+1,jWsi+2,jWsi+3,j,v=x−xixi+1−xi,x∈xi,xi+1,i=−1,0,1,…,m−2;j=−1,0,1,…,n−2,where Vi, j and W(s)i, j are the vertex of the cubic B-spline curve in the direction of lead and tooth depth.If the expressions ofs and v are plugged into (4), the surface deviation δi,j (x,y) can be obtained, which is on the rotary projection of the tooth surface. When i changes from −1 to m−2 and j changes from −1 to n−2, we can obtain m × n patches constituting the deviation surface, shown in Figure 5.Figure 5
Deviation from theoretical tooth surface.


### 2.2. Superposition of Deviation Surface and Theoretical Tooth Surface

Figure6 shows the coordinates of the skew parabolic rack-cutter applied to the machined pinion. Here, M is a point on the surface of the rack-cutter, which coordinate is (u, −acu2, 0) in the coordinate system Sb. Parameters u (Figure 6(a)) and l = OcOa¯ (Figure 6(b)) are the surface parameters of the skew rack-cutter [23]. dp is the location parameter of the parabolic pole. β is the helix angle. The tooth surface of the pinion is enveloped by the rack-cutter surface during processing. According to the meshing principle and enveloping principle [24], the position vector and normal vector of the theoretical tooth surface can be derived and denoted as Rt (u, l) and nt (u, l). The real tooth surface is constituted by superposing the deviation surface and theoretical tooth surface. Then, the position vector and normal vector of the real tooth surface are [25](5)Rru,l=δ′u,l⋅ntu,l+Rtu,l,(6)nru,l=∂Rru,l∂u×∂Rru,l∂l=∂Rtu,l∂u+∂δ′u,l∂u⋅ntu,l+∂ntu,l∂u⋅δ′u,l×∂Rtu,l∂l+∂δ′u,l∂l⋅ntu,l+∂ntu,l∂l⋅δ′u,l,where δʹ (u, l) is the deviation of the tooth surface, which is converted from the deviation of the rotary projection of the tooth surface δi,j (x, y). The conversion can be realized by the relation between (u, l) and (x, y) shown in the following equation:(7)x=Rtx2u,l+Rty2u,l,y=Rtzu,l,where Rtx (u, l), Rty (u, l), and Rtz (u, l) are the coordinate components of Rt (u, l). Next, the partial derivatives ∂δ′u,l/∂u and ∂δ′u,l/∂l in (6) can be computed using the following equation:(8)∂δ′u,l∂u=∂δx,y∂x∂x∂u+∂δx,y∂y∂y∂u,∂δ′u,l∂l=∂δx,y∂x∂x∂l+∂δx,y∂y∂y∂l.Figure 6
Coordinates of skew rack-cutter. (a) Parabolic profile of rack-cutter in normal section. (b) For derivation of rack-cutter surface.
(a)
(b)


## 2.1. Fitting of Surface Deviation

Some scholars have applied the bicubic B-spline method to many studies. Gardner and Gardner [21] used it to model a magnetohydrodynamic duct flow and obtained excellent consistency between the numerical result and the analytic solution even for large values of the applied magnetic field. Na [22] used the method to fit a complex curved surface, and the fitting surface was C2 conforming. This paper also uses a bicubic B-spline to form the deviation surface.Before fitting the surface, some discrete points must be set on the measured modification curve in Figure2. To improve the fitting accuracy, we select more discrete points at both ends of the measured modification curve than in the middle. Finally, 42 points are selected in the meshing line direction and 47 points are selected in the lead direction. Then, the deviations of 42 points in the meshing line direction must be converted into the deviations along the tooth profile direction. According to the meshing relationship shown in Figure 4, the relationship between the tooth depth and the length of the meshing line is as follows:(1)Yprofile=Ymeshing2+rbp2,where Yprofile is the tooth depth, Ymeshing is the length of the meshing line, and rbp is the base radius of the pinion.Figure 4
Relationship between tooth depth and length of meshing line.As shown in Equation (2), the 42 × 47 discrete points are arranged in a matrix.(2)P=P00P01⋯P0nP10P11⋯P1n⋮⋮⋱⋮Pm0Pm1⋯Pmn,m=41,n=46,Based on the 42 × 47 discrete points, we fit the cubic B-spline curves in the directions of lead and tooth depth. The control vertices of the cubic B-spline curve are computed by using the method of reverse vertices, which is presented in Appendix [22]. Then, the cubic B-spline curve equations generated by these vertices are as follows:The lead direction,(3)Li,js=1ss2s3⋅161410−30303−630−13−31Vi,jVi,j+1Vi,j+2Vi,j+3,v=x−xixi+1−xi,x∈xi,xi+1,i=0,1,…,m;j=−1,0,1,…,n−2.The direction of tooth depth,(4)Ii,jv,s=1vv2v3⋅161410−30303−630−13−31Wsi,jWsi+1,jWsi+2,jWsi+3,j,v=x−xixi+1−xi,x∈xi,xi+1,i=−1,0,1,…,m−2;j=−1,0,1,…,n−2,where Vi, j and W(s)i, j are the vertex of the cubic B-spline curve in the direction of lead and tooth depth.If the expressions ofs and v are plugged into (4), the surface deviation δi,j (x,y) can be obtained, which is on the rotary projection of the tooth surface. When i changes from −1 to m−2 and j changes from −1 to n−2, we can obtain m × n patches constituting the deviation surface, shown in Figure 5.Figure 5
Deviation from theoretical tooth surface.


## 2.2. Superposition of Deviation Surface and Theoretical Tooth Surface

Figure6 shows the coordinates of the skew parabolic rack-cutter applied to the machined pinion. Here, M is a point on the surface of the rack-cutter, which coordinate is (u, −acu2, 0) in the coordinate system Sb. Parameters u (Figure 6(a)) and l = OcOa¯ (Figure 6(b)) are the surface parameters of the skew rack-cutter [23]. dp is the location parameter of the parabolic pole. β is the helix angle. The tooth surface of the pinion is enveloped by the rack-cutter surface during processing. According to the meshing principle and enveloping principle [24], the position vector and normal vector of the theoretical tooth surface can be derived and denoted as Rt (u, l) and nt (u, l). The real tooth surface is constituted by superposing the deviation surface and theoretical tooth surface. Then, the position vector and normal vector of the real tooth surface are [25](5)Rru,l=δ′u,l⋅ntu,l+Rtu,l,(6)nru,l=∂Rru,l∂u×∂Rru,l∂l=∂Rtu,l∂u+∂δ′u,l∂u⋅ntu,l+∂ntu,l∂u⋅δ′u,l×∂Rtu,l∂l+∂δ′u,l∂l⋅ntu,l+∂ntu,l∂l⋅δ′u,l,where δʹ (u, l) is the deviation of the tooth surface, which is converted from the deviation of the rotary projection of the tooth surface δi,j (x, y). The conversion can be realized by the relation between (u, l) and (x, y) shown in the following equation:(7)x=Rtx2u,l+Rty2u,l,y=Rtzu,l,where Rtx (u, l), Rty (u, l), and Rtz (u, l) are the coordinate components of Rt (u, l). Next, the partial derivatives ∂δ′u,l/∂u and ∂δ′u,l/∂l in (6) can be computed using the following equation:(8)∂δ′u,l∂u=∂δx,y∂x∂x∂u+∂δx,y∂y∂y∂u,∂δ′u,l∂l=∂δx,y∂x∂x∂l+∂δx,y∂y∂y∂l.Figure 6
Coordinates of skew rack-cutter. (a) Parabolic profile of rack-cutter in normal section. (b) For derivation of rack-cutter surface.
(a)
(b)


## 3. Internal Excitations

### 3.1. Time-Varying Mesh Stiffness

Using the TCA [24] and LTCA [26] on the real tooth surface of the pinion and the involute tooth surface of the gear, the time-varying mesh stiffness was computed. Using the TCA, the transmission error and backlash between teeth were obtained at zero load. Thus, the LTCA uses the finite element method and nonlinear programming method to obtain the comprehensive contact deformations Δn in the normal direction at n contact points of the tooth surface in a meshing period (the comprehensive contact deformation Δn simultaneously considers the effect of other meshed teeth). Then, the comprehensive mesh stiffness km is calculated by the following equation:(9)km=FnΔn,Fn=Ft⋅cosβb,Ft=Tgrbg,where Ft is the tangent force, Tg is the load torque, rbg is the base radius of the gear, Fn is the normal force, and βb is the base helix angle. Finally, the comprehensive mesh stiffness km consists of n discrete values in a meshing period, so it must be interpolated by the piecewise cubic-Hermite interpolation function and expressed as a periodic function by the Fourier series. Based on this idea, the comprehensive time-varying mesh stiffness was computed at 830 N·m of load torque as shown in Figure 7.Figure 7
Time-varying mesh stiffness of gear pair at low speed (50.8 rpm) (the time-varying mesh stiffness at high speed (2104 r/min) is identical to that in this figure).


### 3.2. Gear Meshing Impact

The corner mesh impact mainly includes the meshing-in and meshing-out impacts. Seireg and Houser [27] have experimentally verified that the meshing-in impact had an obviously greater effect of than the meshing-out impact on the gear transmission behavior. Thus, we only consider the meshing-in impact.Wu et al. [28] and Wang et al. [29] solved the meshing-in impact of unmodified and modified gear pairs, respectively. This paper refers to the theory and method of these papers to calculate the meshing-in impact of the practical modified gear pair. The specific solution steps are shown below.After the modification, in the initial phase of the meshing process, the meshing points obtained by TCA do not bear load in the LTCA. Thus, the actually original position of the meshing-in contact is nearly the first point which the LTCA-calculated load sharing factor between teeth is not zero as shown in Figure8. Figure 9(a) depicts the principle of the shock velocity. In the diagram, N1N2 is the theoretical meshing line, ignoring the manufacturing errors and the deformation of tooth; M is the actual original position of the meshing-in contact; M1M2 is the instantaneous meshing line through point M; vp and vg are the tangential speeds of the pinion and gear; vnp and vng are the normal speeds of the pinion and gear; rbp and rbg are the theoretical base radii of the pinion and gear; and r′bg is the actual base radius of the gear. The shock velocity with vs can be computed using the following equation:(10)vp=wp×RM,vg=wg×ρM,vnp=vp⋅nM,vng=vg⋅nM,vs=vnp−vng,where RM, nM are the position vector and normal vector of M, respectively, which can be obtained by the TCA; ρM is the vector from Og to M, with ρM = RM − E ; E is the vector from Of to Og; wp, wg are the angular speeds of pinion and gear.Figure 8
Approximate location of original meshing point (the contact positions on the curve appear in the process from meshing-in to meshing-out for a single tooth).Figure 9
Sketch map of meshing-in impact (on end face): (a) shock velocity principle, (b) dynamic model of meshing-in impact.
(a)
(b)Because of the gear mass and inertia, if the shock speedvs emerges between the pinion and the gear, the impact must be generated at the meshing point M. As shown in Figure 9(b), mep and meg are the equivalent masses of particles for the pinion and gear:(11)mep=Iprbp2,meg=Igr′bg2.According to the theory of impact dynamics, the impact kinetic energy and elastic potential energy are:(12)Ek=12mevs2,me=11/mep+1/meg=IpIgIpr′bg2+Igrbp2,(13)Ep=∫0WsKsxdx=12KsWs2,Fs=Ks⋅Ws,where me is the total equivalent mass of the gear pair; Fs is the maximal meshing-in impact; Ws is the deformation caused by Fs; and Ks is the stiffness at the original meshing point M, which value is the single tooth stiffness at M. Based on the law of conservation of energy, the relation can be obtained as follows:(14)Ek=Ep.By solving Equations (12)–(14), the maximal impact force can be computed as follows:(15)Fs=IpIgIpr′bg2+Igrbp2Ksvs21/2≈IpIgIprbg2+Igrbp2Ksvs21/2.Because the difference betweenr′bg and rbg is small, this paper uses an approximation solution for the impact force by substituting rbg for r′bg [28].We assume that the impact force is a half sine wave; then, it can be expressed as(16)fst=Fssinωst,0≤t≤ts,where ωs is the angular frequency, with ωs = π/ts; ts is the working time, which is defined as the process of changing the shock speed from versus to 0. According to the momentum theorem, the following relation can be obtained:(17)∫0tsfstdt=meΔv=mevs.Using (18), the working time ts [30] can be solved:(18)ts=π2mevsFs.Thus, the impact force is calculated at 830 N·m of load torque as shown in Figure10.Figure 10
Meshing-in impact force at (a) low speed (50.8 rpm) and (b) high speed (2104 rpm).
(a)
(b)


### 3.3. Error Excitation

We consider the cumulative pitch error, eccentric error, and loaded transmission error as the error excitation. Figures11(a) and 11(b) show the measured cumulative pitch errors of the pinion and gear, so the relative cumulative pitch error of the gear pair is equal to the difference between the gear and pinion. Because the installed eccentricity of the gear system is difficult to be measured, we define the eccentric errors of the pinion and gear as 15 μm and 20 μm based on the radial run-out tolerance Fr (adopting the standard of ISO 1328-2:1997). Figure 12 shows the relationship of the eccentric error and the line of action on the end face for the pinion and gear. Ei (i=p,g) is the eccentric error of the pinion or gear, and Φi is the angle between Ei and the line of action on the end face. As shown in Figure 12, Φi can be obtained as follows:(19)Φi=π2−αt+ωit+γi−φ,i=p,g,where αt is the operating pressure angle; φ is the position angle of the center line; ωi is the angular velocity of the pinion or gear; and γi is the initial phase of the eccentric error. Thus, the equivalent error of Ei projected onto the meshing line direction is(20)eEi=Ei⋅cosΦi⋅cosβb,i=p,g,where βb is the base helix angle. Finally, the relative cumulative pitch error projected onto the meshing line direction, eccentric error, and loaded transmission error obtained by the LTCA are superimposed, and the result is the displacement excitation of the gear system, which value in one period is shown in Figure 13.Figure 11
Measured cumulative pitch error: (a) pinion; (b) gear.
(a)
(b)Figure 12
Relationship of eccentric error and line of action on end face.Figure 13
Displacement excitation of gear pair in low-speed situation (50.8 rpm) (the value of displacement excitation in the high-speed situation (2104 rpm) is the same as that in this figure).


## 3.1. Time-Varying Mesh Stiffness

Using the TCA [24] and LTCA [26] on the real tooth surface of the pinion and the involute tooth surface of the gear, the time-varying mesh stiffness was computed. Using the TCA, the transmission error and backlash between teeth were obtained at zero load. Thus, the LTCA uses the finite element method and nonlinear programming method to obtain the comprehensive contact deformations Δn in the normal direction at n contact points of the tooth surface in a meshing period (the comprehensive contact deformation Δn simultaneously considers the effect of other meshed teeth). Then, the comprehensive mesh stiffness km is calculated by the following equation:(9)km=FnΔn,Fn=Ft⋅cosβb,Ft=Tgrbg,where Ft is the tangent force, Tg is the load torque, rbg is the base radius of the gear, Fn is the normal force, and βb is the base helix angle. Finally, the comprehensive mesh stiffness km consists of n discrete values in a meshing period, so it must be interpolated by the piecewise cubic-Hermite interpolation function and expressed as a periodic function by the Fourier series. Based on this idea, the comprehensive time-varying mesh stiffness was computed at 830 N·m of load torque as shown in Figure 7.Figure 7
Time-varying mesh stiffness of gear pair at low speed (50.8 rpm) (the time-varying mesh stiffness at high speed (2104 r/min) is identical to that in this figure).


## 3.2. Gear Meshing Impact

The corner mesh impact mainly includes the meshing-in and meshing-out impacts. Seireg and Houser [27] have experimentally verified that the meshing-in impact had an obviously greater effect of than the meshing-out impact on the gear transmission behavior. Thus, we only consider the meshing-in impact.Wu et al. [28] and Wang et al. [29] solved the meshing-in impact of unmodified and modified gear pairs, respectively. This paper refers to the theory and method of these papers to calculate the meshing-in impact of the practical modified gear pair. The specific solution steps are shown below.After the modification, in the initial phase of the meshing process, the meshing points obtained by TCA do not bear load in the LTCA. Thus, the actually original position of the meshing-in contact is nearly the first point which the LTCA-calculated load sharing factor between teeth is not zero as shown in Figure8. Figure 9(a) depicts the principle of the shock velocity. In the diagram, N1N2 is the theoretical meshing line, ignoring the manufacturing errors and the deformation of tooth; M is the actual original position of the meshing-in contact; M1M2 is the instantaneous meshing line through point M; vp and vg are the tangential speeds of the pinion and gear; vnp and vng are the normal speeds of the pinion and gear; rbp and rbg are the theoretical base radii of the pinion and gear; and r′bg is the actual base radius of the gear. The shock velocity with vs can be computed using the following equation:(10)vp=wp×RM,vg=wg×ρM,vnp=vp⋅nM,vng=vg⋅nM,vs=vnp−vng,where RM, nM are the position vector and normal vector of M, respectively, which can be obtained by the TCA; ρM is the vector from Og to M, with ρM = RM − E ; E is the vector from Of to Og; wp, wg are the angular speeds of pinion and gear.Figure 8
Approximate location of original meshing point (the contact positions on the curve appear in the process from meshing-in to meshing-out for a single tooth).Figure 9
Sketch map of meshing-in impact (on end face): (a) shock velocity principle, (b) dynamic model of meshing-in impact.
(a)
(b)Because of the gear mass and inertia, if the shock speedvs emerges between the pinion and the gear, the impact must be generated at the meshing point M. As shown in Figure 9(b), mep and meg are the equivalent masses of particles for the pinion and gear:(11)mep=Iprbp2,meg=Igr′bg2.According to the theory of impact dynamics, the impact kinetic energy and elastic potential energy are:(12)Ek=12mevs2,me=11/mep+1/meg=IpIgIpr′bg2+Igrbp2,(13)Ep=∫0WsKsxdx=12KsWs2,Fs=Ks⋅Ws,where me is the total equivalent mass of the gear pair; Fs is the maximal meshing-in impact; Ws is the deformation caused by Fs; and Ks is the stiffness at the original meshing point M, which value is the single tooth stiffness at M. Based on the law of conservation of energy, the relation can be obtained as follows:(14)Ek=Ep.By solving Equations (12)–(14), the maximal impact force can be computed as follows:(15)Fs=IpIgIpr′bg2+Igrbp2Ksvs21/2≈IpIgIprbg2+Igrbp2Ksvs21/2.Because the difference betweenr′bg and rbg is small, this paper uses an approximation solution for the impact force by substituting rbg for r′bg [28].We assume that the impact force is a half sine wave; then, it can be expressed as(16)fst=Fssinωst,0≤t≤ts,where ωs is the angular frequency, with ωs = π/ts; ts is the working time, which is defined as the process of changing the shock speed from versus to 0. According to the momentum theorem, the following relation can be obtained:(17)∫0tsfstdt=meΔv=mevs.Using (18), the working time ts [30] can be solved:(18)ts=π2mevsFs.Thus, the impact force is calculated at 830 N·m of load torque as shown in Figure10.Figure 10
Meshing-in impact force at (a) low speed (50.8 rpm) and (b) high speed (2104 rpm).
(a)
(b)


## 3.3. Error Excitation

We consider the cumulative pitch error, eccentric error, and loaded transmission error as the error excitation. Figures11(a) and 11(b) show the measured cumulative pitch errors of the pinion and gear, so the relative cumulative pitch error of the gear pair is equal to the difference between the gear and pinion. Because the installed eccentricity of the gear system is difficult to be measured, we define the eccentric errors of the pinion and gear as 15 μm and 20 μm based on the radial run-out tolerance Fr (adopting the standard of ISO 1328-2:1997). Figure 12 shows the relationship of the eccentric error and the line of action on the end face for the pinion and gear. Ei (i=p,g) is the eccentric error of the pinion or gear, and Φi is the angle between Ei and the line of action on the end face. As shown in Figure 12, Φi can be obtained as follows:(19)Φi=π2−αt+ωit+γi−φ,i=p,g,where αt is the operating pressure angle; φ is the position angle of the center line; ωi is the angular velocity of the pinion or gear; and γi is the initial phase of the eccentric error. Thus, the equivalent error of Ei projected onto the meshing line direction is(20)eEi=Ei⋅cosΦi⋅cosβb,i=p,g,where βb is the base helix angle. Finally, the relative cumulative pitch error projected onto the meshing line direction, eccentric error, and loaded transmission error obtained by the LTCA are superimposed, and the result is the displacement excitation of the gear system, which value in one period is shown in Figure 13.Figure 11
Measured cumulative pitch error: (a) pinion; (b) gear.
(a)
(b)Figure 12
Relationship of eccentric error and line of action on end face.Figure 13
Displacement excitation of gear pair in low-speed situation (50.8 rpm) (the value of displacement excitation in the high-speed situation (2104 rpm) is the same as that in this figure).


## 4. Dynamic Model and Solution Method

### 4.1. Dynamic Model

Figure14 depicts a 12-DOF dynamic model of a helical gear pair. In this picture, the pinion is left-handed and rotates about the z-axis counterclockwise. Both gears are treated as rigid discs, which are connected to each other by a time-varying mesh spring km in the plane of action along the normal direction. The displacement excitation e is shown in Figure 13 on the line of action.Figure 14
12-DOF dynamic model of helical gear pair.Both gears have radial motions in thex and y directions and axial motion in the z direction. Moreover, they have the swinging motion around the x-axis and y-axis and torsional motion around the z-axis. Thus, the equations of motion of the gear pair are [31](21)mp⋅x¨p+kxp⋅xp+cxp⋅x˙p+km⋅δ+cm⋅δ˙⋅cosβb⋅sinψ=−fs⋅cosβb⋅sinψ,mp⋅y¨p+kyp⋅yp+cyp⋅y˙p+km⋅δ+cm⋅δ˙⋅cosβb⋅cosψ=−fs⋅cosβb⋅cosψ,mp⋅z¨p+kzp⋅zp+czp⋅z˙p−km⋅δ+cm⋅δ˙⋅sinβb=fs⋅sinβb,Ixprbp2⋅u¨xp+kuxp⋅uxp+cuxp⋅u˙xp+km⋅δ+cm⋅δ˙⋅sinβb⋅sinψ=−fs⋅sinβb⋅sinψ,Iyprbp2⋅u¨yp+kuyp⋅uyp+cuyp⋅u˙yp+km⋅δ+cm⋅δ˙⋅sinβb⋅cosψ=−fs⋅sinβb⋅cosψ,Izprbp2⋅u¨zp+kuzp⋅uzp+cuzp⋅u˙zp+km⋅δ+cm⋅δ˙⋅cosβb=−fs⋅cosβb+Tprbp,mg⋅x¨g+kxg⋅xg+cxg⋅x˙g−km⋅δ+cm⋅δ˙⋅cosβb⋅sinψ=fs⋅cosβb⋅sinψ,mg⋅y¨g+kyg⋅yg+cyg⋅y˙g−km⋅δ+cm⋅δ˙⋅cosβb⋅cosψ=fs⋅cosβb⋅cosψ,mg⋅z¨g+kzg⋅zg+czg⋅z˙g+km⋅δ+cm⋅δ˙⋅sinβb=−fs⋅sinβb,Ixgrbg2⋅u¨xg+kuxg⋅uxg+cuxg⋅u˙xg+km⋅δ+cm⋅δ˙⋅sinβb⋅sinψ=−fs⋅sinβb⋅sinψ,Iygrbg2⋅u¨yg+kuyg⋅uyg+cuyg⋅u˙yg+km⋅δ+cm⋅δ˙⋅sinβb⋅cosψ=−fs⋅sinβb⋅cosψ,Izgrbg2⋅u¨zg+kuzg⋅uzg+cuzg⋅u˙zg+km⋅δ+cm⋅δ˙⋅cosβb=−fs⋅cosβb+Tgrbg,where mi is the mass of gear i (i=p,g); Ixi, Iyi, and Izi are the moments of inertia around the x-, y- and z-axes for gear i; rbi is the base radius of gear i; βb is the base helix angle; fs is the impact force; Tp is the driving torque; and Tg is the load torque. km and cm are the time-varying mesh stiffness and mesh damping of the engagement pair:(22)cm=2ξk¯mIzpIzgIzprbg2+Izgrbp2,where ξ is the mesh damping ratio, which value is 0.07 in this paper; k¯m is the mean mesh stiffness of the engagement pair. kxi, kyi, kzi, kuxi, kuyi, kuzi and cxi, cyi, czi, cuxi, cuyi, cuzi are the support stiffness and support damping in the directions of the 6-DOF, respectively, for gear i (i=p,g). In Equation (21), δ is the relative displacement between the gear pair in the direction of the meshing line, and it is expressed as follows:(23)δ=xp−xg⋅sinψ+yp−yg⋅cosψ+uzp+uzg⋅cosβb+uxp+uxg⋅sinψ+uyp+uyg⋅cosψ−zp+zg⋅sinβb−e,where ψ is the angle between the plane of action and the positive y-axis, which is defined as(24)ψ=αt−φ,counterclockwiserotationforpinion,αt+φ,clockwiserotationforpinion,where αt is the operating pressure angle, and φ is the position angle of the center line in Figure 14. Equation (21) is written in the matrix form as(25)M⋅q¨+Kb+Km⋅q+Cb+Cm⋅q˙=km⋅e⋅VT+cm⋅e˙⋅VT−fs⋅VT+P,M=diagmpmpmpIxprbp2Iyprbp2Izprbp2mgmgmgIxgrbg2Iygrbg2Izgrbg2,q=xpypzpuxpuypuzpxgygzguxguyguzgT,V=sinψ⋅cosβbcosψ⋅cosβb−sinβbsinψ⋅sinβbcosψ⋅sinβbcosβb−sinψ⋅cosβb−cosψ⋅cosβbsinβbsinψ⋅sinβbcosψ⋅sinβbcosβb,Km=km⋅VT⋅V,Cm=cm⋅VT⋅V,Kb=diagkxpkypkzpkuxpkuypkuzpkxgkygkzgkuxgkuygkuzg,Cb=diagcxpcypczpcuxpcuypcuzpcxgcygczgcuxgcuygcuzg,Cb=diagcxpcypczpcuxpcuypcuzpcxgcygczgcuxgcuygcuzg,P=00000Tprbp00000−TgrbgT,where M is the mass matrix; Kb and Cb are the supporting stiffness matrix and supporting damping matrix; Km and Cm are the meshing stiffness matrix and meshing damping matrix; and V is the projection vector of displacement in the directions of the 12-DOF, which convert to the meshing line direction.


### 4.2. Solution Method

There are two main types of methods to solve the gear system dynamic equation: numerical integration methods and parse methods. Numerical integration methods include the Gill method, Runge-Kutta method, Ritz method, etc. Parse methods include the harmonic balance method, Fourier series method, mode superposition method, etc. Considering of the low-frequency error, the computation period of the dynamic equation hasn meshing periods, and n is equal to the smallest common multiple of teeth number of the pinion and gear. However, distinctive features of the parse method are the quick calculation speed and the ability to directly obtain the steady-state solution. Thus, this paper uses the method of Fourier series to solve the dynamic equations.The general idea is that the response and exciting force are expanded as Fourier series; the coefficient of each Fourier component of the response is ordered to be equal to that of the exciting force; finally, the stable solution is obtained.First, letq, km, and e be expressed as follows:(26)q=q0+Δq,km=km0+Δkm,e=e0+Δe,where q0, km0, and e0 are the mean value of the dynamic transmission error, time-varying mesh stiffness, and error excitation; Δq, Δkm, and Δe are the fluctuating components. By substituting Equation (26) into (25), we obtain a new equation:(27)M⋅Δq¨+C⋅Δq˙+K0⋅q0+ΔK⋅q0+K0⋅Δq+ΔK⋅Δq=km0⋅e0⋅VT+km0⋅Δe⋅VT+Δkm⋅e0⋅VT+Δkm⋅Δe⋅VT+cm⋅Δe˙⋅VT−fs⋅VT+P,where K0=km0⋅VT⋅V+Kb, ΔK=Δkm⋅VT⋅V, C=cm⋅VT⋅V+Cb.Equation (27) is expressed in terms of a time-varying part and a time-invariant part as follows:(28)M⋅Δq¨+C⋅Δq˙+K0⋅Δq=km0⋅Δe⋅VT+cm⋅Δe˙⋅VT−ΔK⋅q0+Δq+Δkm⋅e⋅VT−fs⋅VT,K0⋅q0=km0⋅e0⋅VT+P,For the termΔK·(q0+Δq) on the right side of Equation (28), q0 can be computed using Equation (28); Δq is unknown, so it is approximated by the fluctuating component of the loaded transmission error obtained by the LTCA, which is denoted as Δqs. Let F be equal to the right side of Equation (28), that is,(29)M⋅Δq¨+C⋅Δq˙+K0⋅Δq≈km0⋅Δe⋅VT+cm⋅Δe˙⋅VT−ΔK⋅q0+Δqs+Δkm⋅e⋅VT−fs⋅VT=F,where F is the exciting force of the gear system. So, the dynamic equation (25) has been converted to the constant coefficient differential equation (29).Second, the dynamic responseΔq and exciting force F are expanded into n order Fourier series:(30)Δq=Δq0+∑i=1nAi⋅cosiωt+Bi⋅siniωt,F=F0+∑i=1nCi⋅cosiωt+Di⋅siniωt,where ω is the fundamental frequency, with ω=2πnp/60zg; np is the speed of the pinion; and zg is the number of gear teeth. The deciding rule of n order is that the exciting force F that is expanded as a Fourier series should be consistent with the one before expansion; here, n = 5000.Finally, by plugging theith-order (i = 1, 2, …, n) expressions of the response and the exciting force into Equation (29), the linear algebraic equations can be obtained.(31)−ωi2⋅M+K0ωi⋅C−ωi⋅C−ωi2⋅M+K0AiBi=CiDi,where ωi is the frequency of the ith-order harmonic, and ωi = iω. When Ai and Bi are calculated by Equation (31), the responses Δq and q can be obtained using Equations (30) and (26).


### 4.3. Comparison of Solution Methods

To verify the accuracy and calculation speed of the Fourier series method, a simple example is solved using this method and the Runge–Kutta method. The basic parameters and dynamic equation of the example are shown in Table2 and the following equation:(32)Iprbp2u¨p+cm⋅u˙p+u˙g⋅cosβb−e˙+km⋅up+ug⋅cosβb−e⋅cosβb=Tprbp,Igrbg2u¨g+cm⋅u˙p+u˙g⋅cosβb−e˙+km⋅up+ug⋅cosβb−e⋅cosβb=Tgrbg,where e is the sum of the eccentric errors of the pinion and gear. This paper calculates the relative rotational displacement between the pinion and gear in two cases: case one: the amplitude of the eccentric error of the pinion and gear is 0 μm; case two: both amplitudes are 10 μm. The results of two cases are shown in Figure 15.Table 2
Parameters of simple example.
Parameters
Pinion
Gear
Number of teeth,Z
20
40
Handedness
Left
Right
Normal modulus,mn (mm)
2.5
2.5
Normal pressure angle,α (deg.)
20
20
Helix angle,β (deg.)
25
-25
Tooth width,b (mm)
40
40
Load torque,Tg (N·m)
—
500
Input speed, rpm
2000
—
Moment of inertia,I (kg·m2)
5.97 × 10−4
5.31 × 10−3Figure 15
Relative rotational displacement in (a) case one and (b) case two.
(a)
(b)Figure15 shows that only the amplitude of the mesh frequency wave slightly varies between the results of Fourier series method and Runge–Kutta method, which illustrate that the gear system dynamic equation can be solved using the Fourier series method. In terms of computing time, there is no difference between these two methods in case one; in case two, the computing time of the Fourier series method is approximately 5 s, whereas, the computing time of the Runge–Kutta method is approximately 60 s. If the example in Table 1 is solved by the Runge–Kutta method, the computational cycle will be approximately 20 times longer than that of the example in this section, so the computing time will also significantly increase. To improve the calculation efficiency, this paper adopts the Fourier series method to solve the system response in the following analysis.


## 4.1. Dynamic Model

Figure14 depicts a 12-DOF dynamic model of a helical gear pair. In this picture, the pinion is left-handed and rotates about the z-axis counterclockwise. Both gears are treated as rigid discs, which are connected to each other by a time-varying mesh spring km in the plane of action along the normal direction. The displacement excitation e is shown in Figure 13 on the line of action.Figure 14
12-DOF dynamic model of helical gear pair.Both gears have radial motions in thex and y directions and axial motion in the z direction. Moreover, they have the swinging motion around the x-axis and y-axis and torsional motion around the z-axis. Thus, the equations of motion of the gear pair are [31](21)mp⋅x¨p+kxp⋅xp+cxp⋅x˙p+km⋅δ+cm⋅δ˙⋅cosβb⋅sinψ=−fs⋅cosβb⋅sinψ,mp⋅y¨p+kyp⋅yp+cyp⋅y˙p+km⋅δ+cm⋅δ˙⋅cosβb⋅cosψ=−fs⋅cosβb⋅cosψ,mp⋅z¨p+kzp⋅zp+czp⋅z˙p−km⋅δ+cm⋅δ˙⋅sinβb=fs⋅sinβb,Ixprbp2⋅u¨xp+kuxp⋅uxp+cuxp⋅u˙xp+km⋅δ+cm⋅δ˙⋅sinβb⋅sinψ=−fs⋅sinβb⋅sinψ,Iyprbp2⋅u¨yp+kuyp⋅uyp+cuyp⋅u˙yp+km⋅δ+cm⋅δ˙⋅sinβb⋅cosψ=−fs⋅sinβb⋅cosψ,Izprbp2⋅u¨zp+kuzp⋅uzp+cuzp⋅u˙zp+km⋅δ+cm⋅δ˙⋅cosβb=−fs⋅cosβb+Tprbp,mg⋅x¨g+kxg⋅xg+cxg⋅x˙g−km⋅δ+cm⋅δ˙⋅cosβb⋅sinψ=fs⋅cosβb⋅sinψ,mg⋅y¨g+kyg⋅yg+cyg⋅y˙g−km⋅δ+cm⋅δ˙⋅cosβb⋅cosψ=fs⋅cosβb⋅cosψ,mg⋅z¨g+kzg⋅zg+czg⋅z˙g+km⋅δ+cm⋅δ˙⋅sinβb=−fs⋅sinβb,Ixgrbg2⋅u¨xg+kuxg⋅uxg+cuxg⋅u˙xg+km⋅δ+cm⋅δ˙⋅sinβb⋅sinψ=−fs⋅sinβb⋅sinψ,Iygrbg2⋅u¨yg+kuyg⋅uyg+cuyg⋅u˙yg+km⋅δ+cm⋅δ˙⋅sinβb⋅cosψ=−fs⋅sinβb⋅cosψ,Izgrbg2⋅u¨zg+kuzg⋅uzg+cuzg⋅u˙zg+km⋅δ+cm⋅δ˙⋅cosβb=−fs⋅cosβb+Tgrbg,where mi is the mass of gear i (i=p,g); Ixi, Iyi, and Izi are the moments of inertia around the x-, y- and z-axes for gear i; rbi is the base radius of gear i; βb is the base helix angle; fs is the impact force; Tp is the driving torque; and Tg is the load torque. km and cm are the time-varying mesh stiffness and mesh damping of the engagement pair:(22)cm=2ξk¯mIzpIzgIzprbg2+Izgrbp2,where ξ is the mesh damping ratio, which value is 0.07 in this paper; k¯m is the mean mesh stiffness of the engagement pair. kxi, kyi, kzi, kuxi, kuyi, kuzi and cxi, cyi, czi, cuxi, cuyi, cuzi are the support stiffness and support damping in the directions of the 6-DOF, respectively, for gear i (i=p,g). In Equation (21), δ is the relative displacement between the gear pair in the direction of the meshing line, and it is expressed as follows:(23)δ=xp−xg⋅sinψ+yp−yg⋅cosψ+uzp+uzg⋅cosβb+uxp+uxg⋅sinψ+uyp+uyg⋅cosψ−zp+zg⋅sinβb−e,where ψ is the angle between the plane of action and the positive y-axis, which is defined as(24)ψ=αt−φ,counterclockwiserotationforpinion,αt+φ,clockwiserotationforpinion,where αt is the operating pressure angle, and φ is the position angle of the center line in Figure 14. Equation (21) is written in the matrix form as(25)M⋅q¨+Kb+Km⋅q+Cb+Cm⋅q˙=km⋅e⋅VT+cm⋅e˙⋅VT−fs⋅VT+P,M=diagmpmpmpIxprbp2Iyprbp2Izprbp2mgmgmgIxgrbg2Iygrbg2Izgrbg2,q=xpypzpuxpuypuzpxgygzguxguyguzgT,V=sinψ⋅cosβbcosψ⋅cosβb−sinβbsinψ⋅sinβbcosψ⋅sinβbcosβb−sinψ⋅cosβb−cosψ⋅cosβbsinβbsinψ⋅sinβbcosψ⋅sinβbcosβb,Km=km⋅VT⋅V,Cm=cm⋅VT⋅V,Kb=diagkxpkypkzpkuxpkuypkuzpkxgkygkzgkuxgkuygkuzg,Cb=diagcxpcypczpcuxpcuypcuzpcxgcygczgcuxgcuygcuzg,Cb=diagcxpcypczpcuxpcuypcuzpcxgcygczgcuxgcuygcuzg,P=00000Tprbp00000−TgrbgT,where M is the mass matrix; Kb and Cb are the supporting stiffness matrix and supporting damping matrix; Km and Cm are the meshing stiffness matrix and meshing damping matrix; and V is the projection vector of displacement in the directions of the 12-DOF, which convert to the meshing line direction.


## 4.2. Solution Method

There are two main types of methods to solve the gear system dynamic equation: numerical integration methods and parse methods. Numerical integration methods include the Gill method, Runge-Kutta method, Ritz method, etc. Parse methods include the harmonic balance method, Fourier series method, mode superposition method, etc. Considering of the low-frequency error, the computation period of the dynamic equation hasn meshing periods, and n is equal to the smallest common multiple of teeth number of the pinion and gear. However, distinctive features of the parse method are the quick calculation speed and the ability to directly obtain the steady-state solution. Thus, this paper uses the method of Fourier series to solve the dynamic equations.The general idea is that the response and exciting force are expanded as Fourier series; the coefficient of each Fourier component of the response is ordered to be equal to that of the exciting force; finally, the stable solution is obtained.First, letq, km, and e be expressed as follows:(26)q=q0+Δq,km=km0+Δkm,e=e0+Δe,where q0, km0, and e0 are the mean value of the dynamic transmission error, time-varying mesh stiffness, and error excitation; Δq, Δkm, and Δe are the fluctuating components. By substituting Equation (26) into (25), we obtain a new equation:(27)M⋅Δq¨+C⋅Δq˙+K0⋅q0+ΔK⋅q0+K0⋅Δq+ΔK⋅Δq=km0⋅e0⋅VT+km0⋅Δe⋅VT+Δkm⋅e0⋅VT+Δkm⋅Δe⋅VT+cm⋅Δe˙⋅VT−fs⋅VT+P,where K0=km0⋅VT⋅V+Kb, ΔK=Δkm⋅VT⋅V, C=cm⋅VT⋅V+Cb.Equation (27) is expressed in terms of a time-varying part and a time-invariant part as follows:(28)M⋅Δq¨+C⋅Δq˙+K0⋅Δq=km0⋅Δe⋅VT+cm⋅Δe˙⋅VT−ΔK⋅q0+Δq+Δkm⋅e⋅VT−fs⋅VT,K0⋅q0=km0⋅e0⋅VT+P,For the termΔK·(q0+Δq) on the right side of Equation (28), q0 can be computed using Equation (28); Δq is unknown, so it is approximated by the fluctuating component of the loaded transmission error obtained by the LTCA, which is denoted as Δqs. Let F be equal to the right side of Equation (28), that is,(29)M⋅Δq¨+C⋅Δq˙+K0⋅Δq≈km0⋅Δe⋅VT+cm⋅Δe˙⋅VT−ΔK⋅q0+Δqs+Δkm⋅e⋅VT−fs⋅VT=F,where F is the exciting force of the gear system. So, the dynamic equation (25) has been converted to the constant coefficient differential equation (29).Second, the dynamic responseΔq and exciting force F are expanded into n order Fourier series:(30)Δq=Δq0+∑i=1nAi⋅cosiωt+Bi⋅siniωt,F=F0+∑i=1nCi⋅cosiωt+Di⋅siniωt,where ω is the fundamental frequency, with ω=2πnp/60zg; np is the speed of the pinion; and zg is the number of gear teeth. The deciding rule of n order is that the exciting force F that is expanded as a Fourier series should be consistent with the one before expansion; here, n = 5000.Finally, by plugging theith-order (i = 1, 2, …, n) expressions of the response and the exciting force into Equation (29), the linear algebraic equations can be obtained.(31)−ωi2⋅M+K0ωi⋅C−ωi⋅C−ωi2⋅M+K0AiBi=CiDi,where ωi is the frequency of the ith-order harmonic, and ωi = iω. When Ai and Bi are calculated by Equation (31), the responses Δq and q can be obtained using Equations (30) and (26).


## 4.3. Comparison of Solution Methods

To verify the accuracy and calculation speed of the Fourier series method, a simple example is solved using this method and the Runge–Kutta method. The basic parameters and dynamic equation of the example are shown in Table2 and the following equation:(32)Iprbp2u¨p+cm⋅u˙p+u˙g⋅cosβb−e˙+km⋅up+ug⋅cosβb−e⋅cosβb=Tprbp,Igrbg2u¨g+cm⋅u˙p+u˙g⋅cosβb−e˙+km⋅up+ug⋅cosβb−e⋅cosβb=Tgrbg,where e is the sum of the eccentric errors of the pinion and gear. This paper calculates the relative rotational displacement between the pinion and gear in two cases: case one: the amplitude of the eccentric error of the pinion and gear is 0 μm; case two: both amplitudes are 10 μm. The results of two cases are shown in Figure 15.Table 2
Parameters of simple example.
Parameters
Pinion
Gear
Number of teeth,Z
20
40
Handedness
Left
Right
Normal modulus,mn (mm)
2.5
2.5
Normal pressure angle,α (deg.)
20
20
Helix angle,β (deg.)
25
-25
Tooth width,b (mm)
40
40
Load torque,Tg (N·m)
—
500
Input speed, rpm
2000
—
Moment of inertia,I (kg·m2)
5.97 × 10−4
5.31 × 10−3Figure 15
Relative rotational displacement in (a) case one and (b) case two.
(a)
(b)Figure15 shows that only the amplitude of the mesh frequency wave slightly varies between the results of Fourier series method and Runge–Kutta method, which illustrate that the gear system dynamic equation can be solved using the Fourier series method. In terms of computing time, there is no difference between these two methods in case one; in case two, the computing time of the Fourier series method is approximately 5 s, whereas, the computing time of the Runge–Kutta method is approximately 60 s. If the example in Table 1 is solved by the Runge–Kutta method, the computational cycle will be approximately 20 times longer than that of the example in this section, so the computing time will also significantly increase. To improve the calculation efficiency, this paper adopts the Fourier series method to solve the system response in the following analysis.


## 5. Example Analysis and Experimental Validation

### 5.1. Example Analysis

This paper provides two calculation examples using the proposed method. The example at low speed (50.8 rpm) shows the computation of the quasistatic transmission error, and the example at a high speed (2104 rpm) shows the analysis of the vibration response of the gear system. The basic parameters of the examples are shown in Table1; other system parameters are listed in Table 3.Table 3
Parameters of example system.
Gear parameters
Shaft parameters
m (kg)
I
x (kg·m2)
I
y (kg·m2)
I
z (kg·m2)
Maximum diameter (mm)
Length (mm)
Pinion
8.03
1.62 × 10−2
1.62 × 10−2
1.13 × 10−2
55
310
Gear
38.99
2.17 × 10−1
2.17 × 10−1
3.89 × 10−1
70
315
Bearing parameters
Stiffness
k
x (N/m)
k
y (N/m)
k
z (N/m)
k
ux (N/m)
k
uy (N/m)
k
uz (N/m)
2 × 108
2 × 108
1 × 108
1 × 107
1 × 107
0
Damping
c
x (N·s/m)
c
y (N·s/m)
c
z (N·s/m)
c
ux (N·s/m)
c
uy (N·s/m)
c
uz (N·s/m)
Pinion
4.01 × 102
4.01 × 102
2.83 × 102
40.2
40.2
0
Gear
8.83 × 102
8.83 × 102
6.24 × 102
1.47 × 102
1.47 × 102
0
The mass and inertia of pinion and gear include the mass and inertia of the corresponding shaft.Based on the above examples, the equations of motion in Section4.1 are solved, and the relative vibratory displacement along the line of action on the end face can be obtained. The specific results in the time and frequency domains are shown in Figures 18–25.


### 5.2. Example Calculation by Using the Traditional Method

In order to validate the effectiveness and feasibility of the proposed method, we compare it with the traditional method. The novelty of the proposed method lies in considering the measured manufacturing error in the study of tooth modification. So, in this section, we give two examples by using the traditional method: one is without regard to the manufacturing errors and the other is considering the given value and variation of the manufacturing errors according to the manufacturing tolerance.Case 1. Ignoring the effect of the manufacturing errors.In this case, the gears only have the tooth modification and eccentric errors. The value of modification is the theoretical modification curve in Figure1. The eccentric errors are consistent with the above examples in Section 5.1. We do TCA and LTCA on the theoretical modification tooth surface fitted by a bicubic B-spline and compute the time-varying mesh stiffness and meshing-in impact based on the results of TCA and LTCA. Finally, taking the time-varying mesh stiffness, meshing-in impact, loaded transmission error and eccentric errors as the internal excitations, we do dynamic analysis for the gear system.Case 2. Considering the given value and variation of the manufacturing errors.On the basis of Case one, we increase tooth profile error and cumulative pitch error to the gear system. Refer to the literature [15], we set the variation of the tooth profile error and cumulative pitch error as a harmonic function, shown in the following equation:(33)EMEt=Efα⋅sinωmt+φfα+Efp2⋅sinωr2t+φfp2−Efp1⋅sinωr1t+φfp1,where Efα is the amplitude of the tooth profile error of pinion; Efp1 and Efp2 are the amplitudes of the cumulative pitch errors of pinion and gear. Referring to the standard of ISO 1328-1:1997, we find out Efα = 13 μm, Efp1 = 28 μm, and Efp2 = 36 μm; ωm is the mesh frequency; ωr1 and ωr2 are the rotation frequencies of pinion and gear; and φi (i = fα, fp1, fp2) are the phase angles, which are set as φfα = 0, φfp1 = 3π/2, and φfp2 = π/2.In the process of dynamic analysis, we add the manufacturing errors to the equations of motion directly like (22) in reference [15], whereas the proposed method brings the tooth profile error to the loaded transmission error by TCA and LTCA.All the two cases are calculated at 50.8 rpm and 2104 rpm. The basic parameters of gears, shafts, and bearings are consistent with the above examples in Section5.1. The results are compared and analyzed in Section 5.4.


### 5.3. Experimental Validation

A test system in Figure16 was assembled to measure the system dynamic response. The test system can produce a circulating power flow through a close loop between the test and reaction gears. The test and reaction gears were connected by flexible shafts and elastic couplings, respectively, which were used to transmit a torque to the loop. Circular gratings were installed in the axle end of each test gear to measure the circumferential turning angle, which precision was ±5″. The grating signal was collected using a data acquisition card and an M + P acquisition system as shown in Figure 17.Figure 16
Gear transmission test system to implement proposed measurement scheme.Figure 17
Signal acquisition system.The turning angles of the test gears were converted from the grating signal, and the gear transmission error along the line of action on the end face was calculated using the following equation:(34)Δx=ϕg⋅rbg−ϕp⋅rbp,where ϕp and ϕg were the actual turning angles of the pinion and gear; rbp and rbg were the base radii. The specific results in the time and frequency domains are shown in Figures 18–25.Figure 18
Quasistatic transmission error.Figure 19
Spectrum and partial enlarged view of quasistatic transmission error: (a) experiment; (b) the measured manufacturing errors; (c) no manufacturing errors; (d) the given manufacturing errors.
(a)
(b)
(c)
(d)Figure 20
Mesh frequency component of quasistatic transmission error.Figure 21
Spectrum and partial enlarged view of gear frequency component of quasistatic transmission error: (a) experiment; (b) simulation by using the proposed method.
(a)
(b)Figure 22
Dynamic transmission error.Figure 23
Spectrum and partial enlarged view of dynamic transmission error: (a) experiment; (b) the measured manufacturing errors; (c) no manufacturing errors; (d) the given manufacturing errors.
(a)
(b)
(c)
(d)Figure 24
Mesh frequency component of dynamic transmission error.Figure 25
Spectrum and partial enlarged view of gear frequency component of dynamic transmission error: (a) experiment; (b) simulation by using the proposed method.
(a)
(b)


### 5.4. Comparison and Analysis

#### 5.4.1. Result Analysis at 50.8 rpm

Figures18 and 19 show the quasistatic transmission error of the helical gear pair and its spectrum for the above three examples and the experiment. The amplitudes corresponding to the fundamental frequencies in Figure 19 are summarized in Table 4, where fzge/zgs is the rotation frequency of gear; fzpe/zps is the rotation frequency of pinion; and fme/ms is the mesh frequency.Table 4
The amplitudes corresponding to the fundamental frequencies in experimental and simulation results at 50.8 rpm.
Experimental result
f
zge: 0.197 Hz (μm)
f
zpe: 0.786 Hz (μm)
f
me: 15.92 Hz (μm)
25.21
29.6
3.62
f
zgs: 0.344 Hz (μm)
f
zps: 0.848 Hz (μm)
f
ms: 16.1 Hz (μm)
Simulation result
1
24.91
16.12
2.42
2
20
15.2
0.90
3
49.63
37.77
12.86
The simulation result 1 is corresponding to the case of the measured manufacturing errors, 2 is corresponding to the case of no manufacturing errors, and 3 is corresponding to the case of the given manufacturing errors.Comparing the results in Figures18 and 19, we can see that the transmission error curves in both experimental and simulation results include the fluctuations of the rotation frequency and mesh frequency, and the values of these fundamental frequencies in the simulation results are largely consistent with ones in the experimental results.However, combining Table4, some other conclusions can be drawn:(1)
For the case of the given manufacturing errors, the vibration amplitudes of the fundamental frequencies are clearly larger than the experimental result, which indicates that the amplitudes of the given cumulative pitch error and tooth profile error are so large as to be not practical. This case’s results are certainly not desirable.(2)
For the case of no manufacturing errors, the amplitude corresponding to mesh frequency is far less than the experimental result. So, it is also not advisable to predict the transmission performance of the machined gear system on the basis of the result of this case.(3)
Only in the case of the measured manufacturing errors, the amplitudes of the base frequency components are close to the ones in the experimental results. Having some differences between the two results is because the internal and external excitations in the experiment are more complex than the simulation. By and large, the proposed method can serve as an estimate when the gear transmission experiment cannot be carried out.Through the above analysis, we learn that the proposed method is viable and practical. So, using a Chebyshev high-pass filter, we handle the results of the experiment and the simulation calculated by using the proposed method. The new curves and their FFT spectra after filtering are shown in Figures20 and 21. The displacement curves for both the experimental and simulation results in Figure 20 have the modulation envelope. Table 5 summarizes the sidebands from the experiment and simulation according to Figure 21. In this table, fme − fzpe + fzge, 2fme − fzpe − fzge, fms − fzps, fms − fzgs, 2fms − fzps, and 2fms − fzgs are the lower sidebands; fme + 2fzge, fme + fzpe, 2fme + 2fzge, 2fme + fzpe, 2fme + 2fzpe, fms + fzgs, fms + fzps, 2fms + fzgs, and 2fms + fzps are the upper sidebands. According to the results, the mesh frequency and its harmonic components are all modulated by the rotation frequency signal of the pinion and gear. However, because of the more complex low-frequency excitations in experiment, the distribution of the sidebands is symmetrical for the simulation, but not for the experiment.Table 5
Sidebands in experimental and simulation results at 50.8 rpm.
Experimental result
Base frequency
Sideband
f
me (Hz)
f
me − fzpe + fzge (Hz)
f
me + 2fzge (Hz)
f
me + fzpe (Hz)
15.92
15.33
16.31
16.7
Second frequency
Sideband
2fme (Hz)
2fme − fzpe − fzge (Hz)
2fme + 2fzge (Hz)
2fme + fzpe (Hz)
2fme + 2fzpe (Hz)
31.83
30.85
32.22
32.62
33.4
Simulation result by using the proposed method
Base frequency
Sideband
f
ms (Hz)
f
ms − fzps (Hz)
f
ms − fzgs (Hz)
f
ms + fzgs (Hz)
f
ms + fzps (Hz)
16.1
15.26
15.76
16.45
16.95
Second frequency
Sideband
2fms (Hz)
2fms − fzps (Hz)
2fms − fzgs (Hz)
2fms + fzgs (Hz)
2fms + fzps (Hz)
32.2
31.35
31.85
32.54
33.06


#### 5.4.2. Result Analysis at 2104 rpm

At 2104 rpm, we do the vibration test and the simulations in three cases mentioned in Section5.4.1. Figures 22 and 23 show the dynamic transmission error and its spectrum for the experiment and three simulations. Table 6 summarizes the vibration amplitudes of the main frequencies for four kinds of results, where the meanings of fzge/zgs, fzpe/zps, and fme/ms are same as Table 4.Table 6
The amplitudes corresponding to the fundamental frequencies in experimental and simulation results at 2104 rpm.
Experimental result
f
zge: 14.22 Hz (μm)
f
zpe: 34.53 Hz (μm)
f
me: 666.4 Hz (μm)
45.53
12.86
15.72
f
zgs: 14.23 Hz (μm)
f
zps: 35.08 Hz (μm)
f
ms: 666.2 Hz (μm)
Simulation result
1
24.66
15.93
5.28
2
20
15.21
1.57
3
49.63
37.78
14.47
The simulation result 1 is corresponding to the case of the measured manufacturing errors, 2 is corresponding to the case of no manufacturing errors, and 3 is corresponding to the case of the given manufacturing errors.Comparing with the results calculated at 50.8 rpm, we can see that the vibration amplitudes of mesh frequency and rotation frequencies opinion and gear increase a lot for the experiment. However, to the simulations of three cases, only the vibration amplitude of mesh frequency grows a little. Next, we analyze these phenomena:(1)
For the simulations of three cases, the cumulative pitch error at 2104 rpm is same as one at 50.8 rpm, so the vibration of rotation frequencies is almost unchanged. With the increase of speed, the meshing-in impact force grows, which leads to the vibration of mesh frequency intensified.(2)
For the experiment, because of speed rise, the amplitude corresponding to rotation frequency of gear nearly doubles the one at 50.8 rpm, which may be caused by low-frequency excitations other than the cumulative pitch error. About the vibration of mesh frequency, the extent of increase is much larger than one in the simulations, which may result from meshing-in impact and other new high-frequency excitations.Based on the above analysis, we can draw some conclusions as follows:(1)
For the case of the given measured manufacturing errors, although the amplitudes corresponding to mesh frequency and rotation frequency of gear are close to the ones in the experimental results, the simulation result is still undesirable. This is because the vibration intensified in the experiment is not caused by the cumulative pitch error and tooth profile error. The simulation of this case is distorted.(2)
For the case of no manufacturing errors, the vibration of high frequency is still far less than one in the experimental result. So, it is also totally impracticable.(3)
About the case of the measured manufacturing errors, the waveform and trend of the dynamic transmission error are roughly identical to ones in the experimental results. In a future study, if we attempt to consider other factors such as the bearing dynamic force, dynamic load, and assemble error of gear system, the simulation results would approximate the experiment results more and more. So the proposed method is feasible and worth studying.We analyze the high-frequency components for the experiment and the simulation computed by using the proposed method further on. Figures24 and 25 show the relative vibration displacement and its FFT spectra after filtering out the low-frequency components. The sidebands around the meshing frequency in the experimental results are more complex than those in the simulation results. It suggests that some new low-frequency excitations come into being in the experiment. Table 7 shows the corresponding sidebands in the experiment and simulation.Table 7
Sidebands in experimental and simulation results at 2104 rpm.
Experimental result
Base frequency
Sideband
f
me (Hz)
f
me − fzpe (Hz)
f
me + fzge (Hz)
f
me + fzpe (Hz)
666.4
631.7
680.4
700.7
Second frequency
Sideband
2fme (Hz)
2fme − fzpe (Hz)
2fme − fzge (Hz)
2fme + fzge (Hz)
2fme + fzpe (Hz)
1332
1298
1318
1347
1367
Simulation result by using the proposed method
Base frequency
Sideband
f
ms (Hz)
f
ms − fzps (Hz)
f
ms − fzgs (Hz)
f
ms + fzgs (Hz)
f
ms + fzps (Hz)
666.2
630.8
652.2
680.7
701.5
Second frequency
Sideband
2fms (Hz)
2fms − fzps (Hz)
2fms − fzgs (Hz)
2fms + fzgs (Hz)
2fms + fzps (Hz)
1332
1297
1318
1347
1368


## 5.1. Example Analysis

This paper provides two calculation examples using the proposed method. The example at low speed (50.8 rpm) shows the computation of the quasistatic transmission error, and the example at a high speed (2104 rpm) shows the analysis of the vibration response of the gear system. The basic parameters of the examples are shown in Table1; other system parameters are listed in Table 3.Table 3
Parameters of example system.
Gear parameters
Shaft parameters
m (kg)
I
x (kg·m2)
I
y (kg·m2)
I
z (kg·m2)
Maximum diameter (mm)
Length (mm)
Pinion
8.03
1.62 × 10−2
1.62 × 10−2
1.13 × 10−2
55
310
Gear
38.99
2.17 × 10−1
2.17 × 10−1
3.89 × 10−1
70
315
Bearing parameters
Stiffness
k
x (N/m)
k
y (N/m)
k
z (N/m)
k
ux (N/m)
k
uy (N/m)
k
uz (N/m)
2 × 108
2 × 108
1 × 108
1 × 107
1 × 107
0
Damping
c
x (N·s/m)
c
y (N·s/m)
c
z (N·s/m)
c
ux (N·s/m)
c
uy (N·s/m)
c
uz (N·s/m)
Pinion
4.01 × 102
4.01 × 102
2.83 × 102
40.2
40.2
0
Gear
8.83 × 102
8.83 × 102
6.24 × 102
1.47 × 102
1.47 × 102
0
The mass and inertia of pinion and gear include the mass and inertia of the corresponding shaft.Based on the above examples, the equations of motion in Section4.1 are solved, and the relative vibratory displacement along the line of action on the end face can be obtained. The specific results in the time and frequency domains are shown in Figures 18–25.


## 5.2. Example Calculation by Using the Traditional Method

In order to validate the effectiveness and feasibility of the proposed method, we compare it with the traditional method. The novelty of the proposed method lies in considering the measured manufacturing error in the study of tooth modification. So, in this section, we give two examples by using the traditional method: one is without regard to the manufacturing errors and the other is considering the given value and variation of the manufacturing errors according to the manufacturing tolerance.Case 1. Ignoring the effect of the manufacturing errors.In this case, the gears only have the tooth modification and eccentric errors. The value of modification is the theoretical modification curve in Figure1. The eccentric errors are consistent with the above examples in Section 5.1. We do TCA and LTCA on the theoretical modification tooth surface fitted by a bicubic B-spline and compute the time-varying mesh stiffness and meshing-in impact based on the results of TCA and LTCA. Finally, taking the time-varying mesh stiffness, meshing-in impact, loaded transmission error and eccentric errors as the internal excitations, we do dynamic analysis for the gear system.Case 2. Considering the given value and variation of the manufacturing errors.On the basis of Case one, we increase tooth profile error and cumulative pitch error to the gear system. Refer to the literature [15], we set the variation of the tooth profile error and cumulative pitch error as a harmonic function, shown in the following equation:(33)EMEt=Efα⋅sinωmt+φfα+Efp2⋅sinωr2t+φfp2−Efp1⋅sinωr1t+φfp1,where Efα is the amplitude of the tooth profile error of pinion; Efp1 and Efp2 are the amplitudes of the cumulative pitch errors of pinion and gear. Referring to the standard of ISO 1328-1:1997, we find out Efα = 13 μm, Efp1 = 28 μm, and Efp2 = 36 μm; ωm is the mesh frequency; ωr1 and ωr2 are the rotation frequencies of pinion and gear; and φi (i = fα, fp1, fp2) are the phase angles, which are set as φfα = 0, φfp1 = 3π/2, and φfp2 = π/2.In the process of dynamic analysis, we add the manufacturing errors to the equations of motion directly like (22) in reference [15], whereas the proposed method brings the tooth profile error to the loaded transmission error by TCA and LTCA.All the two cases are calculated at 50.8 rpm and 2104 rpm. The basic parameters of gears, shafts, and bearings are consistent with the above examples in Section5.1. The results are compared and analyzed in Section 5.4.


## 5.3. Experimental Validation

A test system in Figure16 was assembled to measure the system dynamic response. The test system can produce a circulating power flow through a close loop between the test and reaction gears. The test and reaction gears were connected by flexible shafts and elastic couplings, respectively, which were used to transmit a torque to the loop. Circular gratings were installed in the axle end of each test gear to measure the circumferential turning angle, which precision was ±5″. The grating signal was collected using a data acquisition card and an M + P acquisition system as shown in Figure 17.Figure 16
Gear transmission test system to implement proposed measurement scheme.Figure 17
Signal acquisition system.The turning angles of the test gears were converted from the grating signal, and the gear transmission error along the line of action on the end face was calculated using the following equation:(34)Δx=ϕg⋅rbg−ϕp⋅rbp,where ϕp and ϕg were the actual turning angles of the pinion and gear; rbp and rbg were the base radii. The specific results in the time and frequency domains are shown in Figures 18–25.Figure 18
Quasistatic transmission error.Figure 19
Spectrum and partial enlarged view of quasistatic transmission error: (a) experiment; (b) the measured manufacturing errors; (c) no manufacturing errors; (d) the given manufacturing errors.
(a)
(b)
(c)
(d)Figure 20
Mesh frequency component of quasistatic transmission error.Figure 21
Spectrum and partial enlarged view of gear frequency component of quasistatic transmission error: (a) experiment; (b) simulation by using the proposed method.
(a)
(b)Figure 22
Dynamic transmission error.Figure 23
Spectrum and partial enlarged view of dynamic transmission error: (a) experiment; (b) the measured manufacturing errors; (c) no manufacturing errors; (d) the given manufacturing errors.
(a)
(b)
(c)
(d)Figure 24
Mesh frequency component of dynamic transmission error.Figure 25
Spectrum and partial enlarged view of gear frequency component of dynamic transmission error: (a) experiment; (b) simulation by using the proposed method.
(a)
(b)


## 5.4. Comparison and Analysis

### 5.4.1. Result Analysis at 50.8 rpm

Figures18 and 19 show the quasistatic transmission error of the helical gear pair and its spectrum for the above three examples and the experiment. The amplitudes corresponding to the fundamental frequencies in Figure 19 are summarized in Table 4, where fzge/zgs is the rotation frequency of gear; fzpe/zps is the rotation frequency of pinion; and fme/ms is the mesh frequency.Table 4
The amplitudes corresponding to the fundamental frequencies in experimental and simulation results at 50.8 rpm.
Experimental result
f
zge: 0.197 Hz (μm)
f
zpe: 0.786 Hz (μm)
f
me: 15.92 Hz (μm)
25.21
29.6
3.62
f
zgs: 0.344 Hz (μm)
f
zps: 0.848 Hz (μm)
f
ms: 16.1 Hz (μm)
Simulation result
1
24.91
16.12
2.42
2
20
15.2
0.90
3
49.63
37.77
12.86
The simulation result 1 is corresponding to the case of the measured manufacturing errors, 2 is corresponding to the case of no manufacturing errors, and 3 is corresponding to the case of the given manufacturing errors.Comparing the results in Figures18 and 19, we can see that the transmission error curves in both experimental and simulation results include the fluctuations of the rotation frequency and mesh frequency, and the values of these fundamental frequencies in the simulation results are largely consistent with ones in the experimental results.However, combining Table4, some other conclusions can be drawn:(1)
For the case of the given manufacturing errors, the vibration amplitudes of the fundamental frequencies are clearly larger than the experimental result, which indicates that the amplitudes of the given cumulative pitch error and tooth profile error are so large as to be not practical. This case’s results are certainly not desirable.(2)
For the case of no manufacturing errors, the amplitude corresponding to mesh frequency is far less than the experimental result. So, it is also not advisable to predict the transmission performance of the machined gear system on the basis of the result of this case.(3)
Only in the case of the measured manufacturing errors, the amplitudes of the base frequency components are close to the ones in the experimental results. Having some differences between the two results is because the internal and external excitations in the experiment are more complex than the simulation. By and large, the proposed method can serve as an estimate when the gear transmission experiment cannot be carried out.Through the above analysis, we learn that the proposed method is viable and practical. So, using a Chebyshev high-pass filter, we handle the results of the experiment and the simulation calculated by using the proposed method. The new curves and their FFT spectra after filtering are shown in Figures20 and 21. The displacement curves for both the experimental and simulation results in Figure 20 have the modulation envelope. Table 5 summarizes the sidebands from the experiment and simulation according to Figure 21. In this table, fme − fzpe + fzge, 2fme − fzpe − fzge, fms − fzps, fms − fzgs, 2fms − fzps, and 2fms − fzgs are the lower sidebands; fme + 2fzge, fme + fzpe, 2fme + 2fzge, 2fme + fzpe, 2fme + 2fzpe, fms + fzgs, fms + fzps, 2fms + fzgs, and 2fms + fzps are the upper sidebands. According to the results, the mesh frequency and its harmonic components are all modulated by the rotation frequency signal of the pinion and gear. However, because of the more complex low-frequency excitations in experiment, the distribution of the sidebands is symmetrical for the simulation, but not for the experiment.Table 5
Sidebands in experimental and simulation results at 50.8 rpm.
Experimental result
Base frequency
Sideband
f
me (Hz)
f
me − fzpe + fzge (Hz)
f
me + 2fzge (Hz)
f
me + fzpe (Hz)
15.92
15.33
16.31
16.7
Second frequency
Sideband
2fme (Hz)
2fme − fzpe − fzge (Hz)
2fme + 2fzge (Hz)
2fme + fzpe (Hz)
2fme + 2fzpe (Hz)
31.83
30.85
32.22
32.62
33.4
Simulation result by using the proposed method
Base frequency
Sideband
f
ms (Hz)
f
ms − fzps (Hz)
f
ms − fzgs (Hz)
f
ms + fzgs (Hz)
f
ms + fzps (Hz)
16.1
15.26
15.76
16.45
16.95
Second frequency
Sideband
2fms (Hz)
2fms − fzps (Hz)
2fms − fzgs (Hz)
2fms + fzgs (Hz)
2fms + fzps (Hz)
32.2
31.35
31.85
32.54
33.06


### 5.4.2. Result Analysis at 2104 rpm

At 2104 rpm, we do the vibration test and the simulations in three cases mentioned in Section5.4.1. Figures 22 and 23 show the dynamic transmission error and its spectrum for the experiment and three simulations. Table 6 summarizes the vibration amplitudes of the main frequencies for four kinds of results, where the meanings of fzge/zgs, fzpe/zps, and fme/ms are same as Table 4.Table 6
The amplitudes corresponding to the fundamental frequencies in experimental and simulation results at 2104 rpm.
Experimental result
f
zge: 14.22 Hz (μm)
f
zpe: 34.53 Hz (μm)
f
me: 666.4 Hz (μm)
45.53
12.86
15.72
f
zgs: 14.23 Hz (μm)
f
zps: 35.08 Hz (μm)
f
ms: 666.2 Hz (μm)
Simulation result
1
24.66
15.93
5.28
2
20
15.21
1.57
3
49.63
37.78
14.47
The simulation result 1 is corresponding to the case of the measured manufacturing errors, 2 is corresponding to the case of no manufacturing errors, and 3 is corresponding to the case of the given manufacturing errors.Comparing with the results calculated at 50.8 rpm, we can see that the vibration amplitudes of mesh frequency and rotation frequencies opinion and gear increase a lot for the experiment. However, to the simulations of three cases, only the vibration amplitude of mesh frequency grows a little. Next, we analyze these phenomena:(1)
For the simulations of three cases, the cumulative pitch error at 2104 rpm is same as one at 50.8 rpm, so the vibration of rotation frequencies is almost unchanged. With the increase of speed, the meshing-in impact force grows, which leads to the vibration of mesh frequency intensified.(2)
For the experiment, because of speed rise, the amplitude corresponding to rotation frequency of gear nearly doubles the one at 50.8 rpm, which may be caused by low-frequency excitations other than the cumulative pitch error. About the vibration of mesh frequency, the extent of increase is much larger than one in the simulations, which may result from meshing-in impact and other new high-frequency excitations.Based on the above analysis, we can draw some conclusions as follows:(1)
For the case of the given measured manufacturing errors, although the amplitudes corresponding to mesh frequency and rotation frequency of gear are close to the ones in the experimental results, the simulation result is still undesirable. This is because the vibration intensified in the experiment is not caused by the cumulative pitch error and tooth profile error. The simulation of this case is distorted.(2)
For the case of no manufacturing errors, the vibration of high frequency is still far less than one in the experimental result. So, it is also totally impracticable.(3)
About the case of the measured manufacturing errors, the waveform and trend of the dynamic transmission error are roughly identical to ones in the experimental results. In a future study, if we attempt to consider other factors such as the bearing dynamic force, dynamic load, and assemble error of gear system, the simulation results would approximate the experiment results more and more. So the proposed method is feasible and worth studying.We analyze the high-frequency components for the experiment and the simulation computed by using the proposed method further on. Figures24 and 25 show the relative vibration displacement and its FFT spectra after filtering out the low-frequency components. The sidebands around the meshing frequency in the experimental results are more complex than those in the simulation results. It suggests that some new low-frequency excitations come into being in the experiment. Table 7 shows the corresponding sidebands in the experiment and simulation.Table 7
Sidebands in experimental and simulation results at 2104 rpm.
Experimental result
Base frequency
Sideband
f
me (Hz)
f
me − fzpe (Hz)
f
me + fzge (Hz)
f
me + fzpe (Hz)
666.4
631.7
680.4
700.7
Second frequency
Sideband
2fme (Hz)
2fme − fzpe (Hz)
2fme − fzge (Hz)
2fme + fzge (Hz)
2fme + fzpe (Hz)
1332
1298
1318
1347
1367
Simulation result by using the proposed method
Base frequency
Sideband
f
ms (Hz)
f
ms − fzps (Hz)
f
ms − fzgs (Hz)
f
ms + fzgs (Hz)
f
ms + fzps (Hz)
666.2
630.8
652.2
680.7
701.5
Second frequency
Sideband
2fms (Hz)
2fms − fzps (Hz)
2fms − fzgs (Hz)
2fms + fzgs (Hz)
2fms + fzps (Hz)
1332
1297
1318
1347
1368


## 5.4.1. Result Analysis at 50.8 rpm

Figures18 and 19 show the quasistatic transmission error of the helical gear pair and its spectrum for the above three examples and the experiment. The amplitudes corresponding to the fundamental frequencies in Figure 19 are summarized in Table 4, where fzge/zgs is the rotation frequency of gear; fzpe/zps is the rotation frequency of pinion; and fme/ms is the mesh frequency.Table 4
The amplitudes corresponding to the fundamental frequencies in experimental and simulation results at 50.8 rpm.
Experimental result
f
zge: 0.197 Hz (μm)
f
zpe: 0.786 Hz (μm)
f
me: 15.92 Hz (μm)
25.21
29.6
3.62
f
zgs: 0.344 Hz (μm)
f
zps: 0.848 Hz (μm)
f
ms: 16.1 Hz (μm)
Simulation result
1
24.91
16.12
2.42
2
20
15.2
0.90
3
49.63
37.77
12.86
The simulation result 1 is corresponding to the case of the measured manufacturing errors, 2 is corresponding to the case of no manufacturing errors, and 3 is corresponding to the case of the given manufacturing errors.Comparing the results in Figures18 and 19, we can see that the transmission error curves in both experimental and simulation results include the fluctuations of the rotation frequency and mesh frequency, and the values of these fundamental frequencies in the simulation results are largely consistent with ones in the experimental results.However, combining Table4, some other conclusions can be drawn:(1)
For the case of the given manufacturing errors, the vibration amplitudes of the fundamental frequencies are clearly larger than the experimental result, which indicates that the amplitudes of the given cumulative pitch error and tooth profile error are so large as to be not practical. This case’s results are certainly not desirable.(2)
For the case of no manufacturing errors, the amplitude corresponding to mesh frequency is far less than the experimental result. So, it is also not advisable to predict the transmission performance of the machined gear system on the basis of the result of this case.(3)
Only in the case of the measured manufacturing errors, the amplitudes of the base frequency components are close to the ones in the experimental results. Having some differences between the two results is because the internal and external excitations in the experiment are more complex than the simulation. By and large, the proposed method can serve as an estimate when the gear transmission experiment cannot be carried out.Through the above analysis, we learn that the proposed method is viable and practical. So, using a Chebyshev high-pass filter, we handle the results of the experiment and the simulation calculated by using the proposed method. The new curves and their FFT spectra after filtering are shown in Figures20 and 21. The displacement curves for both the experimental and simulation results in Figure 20 have the modulation envelope. Table 5 summarizes the sidebands from the experiment and simulation according to Figure 21. In this table, fme − fzpe + fzge, 2fme − fzpe − fzge, fms − fzps, fms − fzgs, 2fms − fzps, and 2fms − fzgs are the lower sidebands; fme + 2fzge, fme + fzpe, 2fme + 2fzge, 2fme + fzpe, 2fme + 2fzpe, fms + fzgs, fms + fzps, 2fms + fzgs, and 2fms + fzps are the upper sidebands. According to the results, the mesh frequency and its harmonic components are all modulated by the rotation frequency signal of the pinion and gear. However, because of the more complex low-frequency excitations in experiment, the distribution of the sidebands is symmetrical for the simulation, but not for the experiment.Table 5
Sidebands in experimental and simulation results at 50.8 rpm.
Experimental result
Base frequency
Sideband
f
me (Hz)
f
me − fzpe + fzge (Hz)
f
me + 2fzge (Hz)
f
me + fzpe (Hz)
15.92
15.33
16.31
16.7
Second frequency
Sideband
2fme (Hz)
2fme − fzpe − fzge (Hz)
2fme + 2fzge (Hz)
2fme + fzpe (Hz)
2fme + 2fzpe (Hz)
31.83
30.85
32.22
32.62
33.4
Simulation result by using the proposed method
Base frequency
Sideband
f
ms (Hz)
f
ms − fzps (Hz)
f
ms − fzgs (Hz)
f
ms + fzgs (Hz)
f
ms + fzps (Hz)
16.1
15.26
15.76
16.45
16.95
Second frequency
Sideband
2fms (Hz)
2fms − fzps (Hz)
2fms − fzgs (Hz)
2fms + fzgs (Hz)
2fms + fzps (Hz)
32.2
31.35
31.85
32.54
33.06


## 5.4.2. Result Analysis at 2104 rpm

At 2104 rpm, we do the vibration test and the simulations in three cases mentioned in Section5.4.1. Figures 22 and 23 show the dynamic transmission error and its spectrum for the experiment and three simulations. Table 6 summarizes the vibration amplitudes of the main frequencies for four kinds of results, where the meanings of fzge/zgs, fzpe/zps, and fme/ms are same as Table 4.Table 6
The amplitudes corresponding to the fundamental frequencies in experimental and simulation results at 2104 rpm.
Experimental result
f
zge: 14.22 Hz (μm)
f
zpe: 34.53 Hz (μm)
f
me: 666.4 Hz (μm)
45.53
12.86
15.72
f
zgs: 14.23 Hz (μm)
f
zps: 35.08 Hz (μm)
f
ms: 666.2 Hz (μm)
Simulation result
1
24.66
15.93
5.28
2
20
15.21
1.57
3
49.63
37.78
14.47
The simulation result 1 is corresponding to the case of the measured manufacturing errors, 2 is corresponding to the case of no manufacturing errors, and 3 is corresponding to the case of the given manufacturing errors.Comparing with the results calculated at 50.8 rpm, we can see that the vibration amplitudes of mesh frequency and rotation frequencies opinion and gear increase a lot for the experiment. However, to the simulations of three cases, only the vibration amplitude of mesh frequency grows a little. Next, we analyze these phenomena:(1)
For the simulations of three cases, the cumulative pitch error at 2104 rpm is same as one at 50.8 rpm, so the vibration of rotation frequencies is almost unchanged. With the increase of speed, the meshing-in impact force grows, which leads to the vibration of mesh frequency intensified.(2)
For the experiment, because of speed rise, the amplitude corresponding to rotation frequency of gear nearly doubles the one at 50.8 rpm, which may be caused by low-frequency excitations other than the cumulative pitch error. About the vibration of mesh frequency, the extent of increase is much larger than one in the simulations, which may result from meshing-in impact and other new high-frequency excitations.Based on the above analysis, we can draw some conclusions as follows:(1)
For the case of the given measured manufacturing errors, although the amplitudes corresponding to mesh frequency and rotation frequency of gear are close to the ones in the experimental results, the simulation result is still undesirable. This is because the vibration intensified in the experiment is not caused by the cumulative pitch error and tooth profile error. The simulation of this case is distorted.(2)
For the case of no manufacturing errors, the vibration of high frequency is still far less than one in the experimental result. So, it is also totally impracticable.(3)
About the case of the measured manufacturing errors, the waveform and trend of the dynamic transmission error are roughly identical to ones in the experimental results. In a future study, if we attempt to consider other factors such as the bearing dynamic force, dynamic load, and assemble error of gear system, the simulation results would approximate the experiment results more and more. So the proposed method is feasible and worth studying.We analyze the high-frequency components for the experiment and the simulation computed by using the proposed method further on. Figures24 and 25 show the relative vibration displacement and its FFT spectra after filtering out the low-frequency components. The sidebands around the meshing frequency in the experimental results are more complex than those in the simulation results. It suggests that some new low-frequency excitations come into being in the experiment. Table 7 shows the corresponding sidebands in the experiment and simulation.Table 7
Sidebands in experimental and simulation results at 2104 rpm.
Experimental result
Base frequency
Sideband
f
me (Hz)
f
me − fzpe (Hz)
f
me + fzge (Hz)
f
me + fzpe (Hz)
666.4
631.7
680.4
700.7
Second frequency
Sideband
2fme (Hz)
2fme − fzpe (Hz)
2fme − fzge (Hz)
2fme + fzge (Hz)
2fme + fzpe (Hz)
1332
1298
1318
1347
1367
Simulation result by using the proposed method
Base frequency
Sideband
f
ms (Hz)
f
ms − fzps (Hz)
f
ms − fzgs (Hz)
f
ms + fzgs (Hz)
f
ms + fzps (Hz)
666.2
630.8
652.2
680.7
701.5
Second frequency
Sideband
2fms (Hz)
2fms − fzps (Hz)
2fms − fzgs (Hz)
2fms + fzgs (Hz)
2fms + fzps (Hz)
1332
1297
1318
1347
1368


## 6. Summary

At 50.8 rpm and 2104 rpm, we simulate a helical gear system by using the proposed method and the traditional method and perform a test. Through comparative analysis of the whole results, it could be seen that the proposed method in this paper is significant and offers a fertile field for study.


## 7. Conclusions

This paper has incorporated the measured tooth profile error and cumulative pitch error into a gear dynamics model and calculated the time-varying mesh stiffness and meshing-in impact using the TCA and LTCA on an actual tooth surface. Because the period of the relative pitch error between pinion and gear is long, the fast and effective Fourier series method is used. Considering a pair of helical gears as an example, this paper does simulation and a verification test at 50.8 rpm and 2104 rpm, and the simulation results are compared with ones computed by using the traditional method. Based on the analysis results, the following conclusions can be drawn:(1)
Because of the quasistatic analysis based on the actual fitting tooth surface, the vibration amplitude of mesh frequency for the proposed method is closest to the corresponding experimental result among the three cases of simulation.(2)
At 2104 rpm, the waveform and trend of the dynamic transmission error for the proposed method are roughly identical to ones in the experimental results. However, the amplitudes of the fundamental frequency components are smaller than that in the experimental results. It shows that the proposed method offers a fertile field for study in order to be pretty approximate to the experiment results at a high speed.(3)
The traditional methods of either no manufacturing errors or the given manufacturing errors are too distorted to predict the actual transmission performance of gear system.(4)
In contrast, the proposed analysis method is reasonable, feasible and interesting to study.



---

*Source: 9645453-2018-11-04.xml*