Dynamic Analysis of Planetary Gear with Root Crack in Sun Gear Based on Improved Time-Varying Mesh Stiffness

The time-varying mesh stiffness (TVMS) is the crucial parameter of the dynamic model of the gear system. Accurate calculation of TVMS is essential for effective fault diagnosis of the gear system. A mesh stiffness improved method considering both the gear tooth transition curve and deformation of the gear body is presented in this paper, and the stiffness calculation expressions under healthy and crack states are given, respectively. Based on the lumped parameter method, the dynamic model of planetary gear with crack in sun gear is constructed, and the vibration response is solved. The simulation results show when the tooth root cracks appear, the vibration response of the tooth has obvious shock response characteristics. The characteristic frequency and frequency multiplication of sun gear fault can be found obviously by envelope analysis. The simulation signal and the test signal obtained by the drivetrain dynamics simulator gearbox test rig are compared and verified. The comparison and verification results show that the proposed TVMS calculation method and the dynamic model are accurate, which can provide a certain theoretical basis for the fault diagnosis of planetary gear with crack fault. The test results are consistent with the numerical analysis results, which provides a theoretical basis for the fault diagnosis of tooth crack.


Introduction
The planetary gear transmission system has the advantages of high transmission efficiency and compactness performance, etc. Therefore, it is widely used in aeronautics and astronautics, wind turbines, construction machinery, and so on; furthermore, it is a crucial component for transferring torque and motion. Due to the change of the pairs and position of gear teeth during the meshing, the mesh stiffness of gear changes constantly. The time-varying characteristic of mesh stiffness is a crucial parameter excitation source, which has a significant influence on the dynamic response characteristic of the gear transmission system. Therefore, how to accurately calculate the tooth mesh stiffness has become a focus of gear dynamics research.
So far, scholars have carried out some research on the calculation methods of mesh stiffness in healthy gear. In the early 1940s, Weber and Banascheck [1] calculated the displacement of a cylindrical To sum up, most of the literature decompose the mesh stiffness of gears into compression stiffness, bending stiffness, shear stiffness, and Hertz contact stiffness. Then these kinds of stiffness are synthesized, and the total stiffness is obtained. In these studies, especially for planetary gears, a few scholars consider the deformation stiffness of the gear body, and a few scholars consider the influence of tooth root transition on the stiffness solution. However, in the process of solving the stiffness of planetary gear, few scholars consider the influence factors of the fillet transition part of the tooth root and gear body at the same time. Therefore, in order to make up for this defect, on the one hand, the influence factors of the root fillet transition parts should be considered, and on the other hand, the influence of the gear body deformation should also be considered. On this basis, the formulas for calculating the mesh stiffness of gear teeth under healthy and fault conditions are derived respectively.
The structure of this paper is as follows: In Section 2, the gear meshing stiffness is evaluated. The meshing stiffness of healthy gears and cracked gears are resolved in Sections 2.1 and 2.2, respectively. In Section 3, a dynamic model with eight degrees of freedom is established. In Section 4, the simulation and experimental results are analyzed and discussed, and finally, conclusions are given in Section 5.

Stiffness Calculation
The gear tooth number engaged in meshing and meshing position change during the meshing process; as a result, gear mesh stiffness demonstrates strong time-varying property. An illustration of the meshing process is shown in Figure 1. TVMS is one of the main excitation sources of gear system vibration. When the tooth cracks, the damaged gear tooth is weaker than the health one. Thus, the TVMS of the gear will be decreased, which will lead to the change of the vibration characteristics of the gear system. If the stiffness values of different crack sizes can be calculated quantitatively, the vibration signal at the meshing point of each pair of gears can be obtained through solution of dynamic equation. By processing the vibration signal, we can achieve the purpose of monitoring the state of the system. Therefore, an accurate solution of mesh stiffness is conducive to dynamic simulation analysis and accurate dynamic vibration characteristics of gears.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 21 established an analytical expression for calculating the gear mesh stiffness and built a six-DOF gear dynamics model considering the mesh stiffness and sliding friction factors. To sum up, most of the literature decompose the mesh stiffness of gears into compression stiffness, bending stiffness, shear stiffness, and Hertz contact stiffness. Then these kinds of stiffness are synthesized, and the total stiffness is obtained. In these studies, especially for planetary gears, a few scholars consider the deformation stiffness of the gear body, and a few scholars consider the influence of tooth root transition on the stiffness solution. However, in the process of solving the stiffness of planetary gear, few scholars consider the influence factors of the fillet transition part of the tooth root and gear body at the same time. Therefore, in order to make up for this defect, on the one hand, the influence factors of the root fillet transition parts should be considered, and on the other hand, the influence of the gear body deformation should also be considered. On this basis, the formulas for calculating the mesh stiffness of gear teeth under healthy and fault conditions are derived respectively.
The structure of this paper is as follows: In Section 2, the gear meshing stiffness is evaluated. The meshing stiffness of healthy gears and cracked gears are resolved in Sections 2.1 and 2.2, respectively. In Section 3, a dynamic model with eight degrees of freedom is established. In Section 4, the simulation and experimental results are analyzed and discussed, and finally, conclusions are given in Section 5.

Stiffness Calculation
The gear tooth number engaged in meshing and meshing position change during the meshing process; as a result, gear mesh stiffness demonstrates strong time-varying property. An illustration of the meshing process is shown in Figure 1. TVMS is one of the main excitation sources of gear system vibration. When the tooth cracks, the damaged gear tooth is weaker than the health one. Thus, the TVMS of the gear will be decreased, which will lead to the change of the vibration characteristics of the gear system. If the stiffness values of different crack sizes can be calculated quantitatively, the vibration signal at the meshing point of each pair of gears can be obtained through solution of dynamic equation. By processing the vibration signal, we can achieve the purpose of monitoring the state of the system. Therefore, an accurate solution of mesh stiffness is conducive to dynamic simulation analysis and accurate dynamic vibration characteristics of gears.

Mesh Stiffness Calculation of Healthy Gears
The gear tooth is simplified as cantilever beams in this section. For a standard external meshing tooth, when the number of teeth is 42, the base circle and the root circle of the gear coincide. When the number of teeth is less than 42, the root circle is smaller than the base circle. Otherwise, it is larger than the base circle. Therefore, the next two cases are analyzed.

Mesh Stiffness Calculation of Healthy Gears
The gear tooth is simplified as cantilever beams in this section. For a standard external meshing tooth, when the number of teeth is 42, the base circle and the root circle of the gear coincide. When the number of teeth is less than 42, the root circle is smaller than the base circle. Otherwise, it is larger than the base circle. Therefore, the next two cases are analyzed. AB is an additive curve, BC is the involved curve of gear teeth, and CD is the transition curve of the tooth profile, which is related to the tool path, as shown in Figure 3. According to [24], the gear tooth root transition curve equation can be expressed as: where R is the pitch circle radius, a1 is the distance between the tool fillet center and the meshing line, rp is the tool radius, b1 is the distance between the tool fillet center and the cogging centerline, γ is the variable (α0 ≤ γ ≤ π/2), ha* and c* are the addendum height coefficient and the addendum clearance coefficient, respectively, and α0 is the pressure angle of the pitch circle.  AB is an additive curve, BC is the involved curve of gear teeth, and CD is the transition curve of the tooth profile, which is related to the tool path, as shown in Figure 3. When the number of teeth is less than 42, the tooth beam model is shown in Figure 2. AB is an additive curve, BC is the involved curve of gear teeth, and CD is the transition curve of the tooth profile, which is related to the tool path, as shown in Figure 3. According to [24], the gear tooth root transition curve equation can be expressed as: where R is the pitch circle radius, a1 is the distance between the tool fillet center and the meshing line, rp is the tool radius, b1 is the distance between the tool fillet center and the cogging centerline, γ is the variable (α0 ≤ γ ≤ π/2), ha* and c* are the addendum height coefficient and the addendum clearance coefficient, respectively, and α0 is the pressure angle of the pitch circle.  According to [24], the gear tooth root transition curve equation can be expressed as: where R is the pitch circle radius, a 1 is the distance between the tool fillet center and the meshing line, r p is the tool radius, b 1 is the distance between the tool fillet center and the cogging centerline, γ is the variable (α 0 ≤ γ ≤ π/2), h a * and c* are the addendum height coefficient and the addendum clearance coefficient, respectively, and α 0 is the pressure angle of the pitch circle.
Appl. Sci. 2020, 10, 8379 5 of 20 Four parts of potential energy will be stored during gear meshing, including Hertz energy U h , bending potential energy U b , radial compression denaturation energy U a , and shear deformation energy U b . Based on the potential energy mentioned above, Hertz stiffness k h , bending stiffness k b , radial compression stiffness k a, and shear stiffness k s can be calculated respectively. According to the mechanics of materials and elasticity, the following conclusions can be obtained: where F is the interaction force at the mesh point, which can be decomposed into radial force F a and tangential force F b , E is the modulus of elasticity, G is the modulus of shear, Ix and Ix 1 are the moments of inertia of the gear section at the distance of x and x 1 from the base circle, A x and A x1 are the cross-sectional areas, d is the distance between the meshing point and the base circle, h is the distance between the meshing point and the gear pair line, R b and R r is the base circle radius and the root circle radius respectively. Combined with Figure 2 and the characteristics of the involute curve, it can be concluded that: where α 2 = π 2N + tan α 0 − α 0 , α 3 = arcsin(R b sin α 2 /R r ). The calculation formulas of the base circle radius and root circle radius are as follows: where m is the modulus, z is the number of gear teeth, α is the tooth pressure angle, h a * and c* is the coefficient of addendum height and clearance respectively. For standard gears, h a * = 1, c* = 0.25, α = 20 • . According to Formulas (18) and (19), the radius of the root circle and base circle are almost equal when the number of teeth is 42.
Based on the Hertz theory, the Hertz stiffness of a pair of mesh teeth of the same material on the meshing line is independent of the contact position, thus: where E is the elastic modulus, L is the gear axial thickness, v is the Poisson ratio. By substituting Equations (10), (12), (14), (15), (18) and (19) into Equation (7), the bending stiffness can be obtained: Similarly, shear stiffness and radial compression stiffness can be obtained as follows:

The Diameter of the Root Circle Is Larger than the Base Circle
When the number of teeth is greater than 42, the tooth beam model is shown in Figure 4.
where m is the modulus, z is the number of gear teeth, α is the tooth pressure angle, ha* and c* is the coefficient of addendum height and clearance respectively. For standard gears, ha* = 1, c* = 0.25, α = 20°. According to Formulas (18) and (19), the radius of the root circle and base circle are almost equal when the number of teeth is 42. Based on the Hertz theory, the Hertz stiffness of a pair of mesh teeth of the same material on the meshing line is independent of the contact position, thus: where E is the elastic modulus, L is the gear axial thickness, v is the Poisson ratio. By substituting Equations (10), (12), (14), (15), (18) and (19) into Equation (7), the bending stiffness can be obtained: Similarly, shear stiffness and radial compression stiffness can be obtained as follows:

The Diameter of the Root Circle Is Larger than the Base Circle
When the number of teeth is greater than 42, the tooth beam model is shown in Figure 4. When the tooth root circle is larger than the base circle, the bending potential energy, shear potential energy, and axial compression potential energy formula should be modified as follows: When the tooth root circle is larger than the base circle, the bending potential energy, shear potential energy, and axial compression potential energy formula should be modified as follows: Similarly, it can be concluded that:

Calculation of Comprehensive Mesh Stiffness
For the above two cases, since the gear is not a rigid body absolutely, the flexible deformation of the gear matrix should be considered when calculating the comprehensive stiffness, and the corresponding flexible stiffness can be expressed as [4]: where L is the gear tooth width, µ f and S f are shown in Figure 5. The coefficients L*, M*, P* and Q* can be calculated from the following polynomials.
where, X * i can be the coefficient L*, M*, P* and Q*. h fi = r f /r int . The definition of θ f is given in Figure 5. The parameter values in the Formula (31) are shown in reference [25].
Similarly, it can be concluded that:

Calculation of Comprehensive Mesh Stiffness
For the above two cases, since the gear is not a rigid body absolutely, the flexible deformation of the gear matrix should be considered when calculating the comprehensive stiffness, and the corresponding flexible stiffness can be expressed as [4]: where L is the gear tooth width, μf and Sf are shown in Figure 5. The coefficients L*, M*, P* and Q* can be calculated from the following polynomials.
( ) where, * i X can be the coefficient L*, M*, P* and Q*. hfi = rf/rint. The definition of θf is given in Figure   5. The parameter values in the Formula (31) are shown in reference [25]. After considering four kinds of potential energy and flexible deformation of gear matrix, the comprehensive mesh stiffness of a pair of gear tooth can be expressed as follows: After considering four kinds of potential energy and flexible deformation of gear matrix, the comprehensive mesh stiffness of a pair of gear tooth can be expressed as follows: Appl. Sci. 2020, 10, 8379 8 of 20 The total potential energy of a pair of mesh gears can be expressed as follows: where the subscripts 1 and 2 represent the driving gear and driven gear, respectively. For example, k b1 indicates the bending stiffness of the driving gear, and other symbols follow similar naming rules.

Mesh Stiffness Calculation of Crack Gear
Due to the influence of overload, fatigue, and other factors on the tooth root, after the planetary gear transmission system works for a long time, the roots of the gear teeth will fatigue cracks occur. Figure 6 is a schematic diagram of a gear tooth crack. It is assumed that the root crack of a gear tooth is a through type crack and its propagation path is a straight line.
The total potential energy of a pair of mesh gears can be expressed as follows: where the subscripts 1 and 2 represent the driving gear and driven gear, respectively. For example, kb1 indicates the bending stiffness of the driving gear, and other symbols follow similar naming rules.

Mesh Stiffness Calculation of Crack Gear
Due to the influence of overload, fatigue, and other factors on the tooth root, after the planetary gear transmission system works for a long time, the roots of the gear teeth will fatigue cracks occur. Figure 6 is a schematic diagram of a gear tooth crack. It is assumed that the root crack of a gear tooth is a through type crack and its propagation path is a straight line.  Figure 7 shows the crack schematic diagram of gear teeth when the root circle diameter is smaller than the base circle diameter. q1 is the initial crack length, which does not reach the centerline of the tooth, and q2 is that the crack size exceeds the centerline of the tooth. v is the angle between the crack line and the centerline of the tooth.
As in the case of healthy gear, the combined mesh stiffness consists of Hertz stiffness, bending stiffness, shear stiffness, and axial compression stiffness. The Hertz stiffness does not change when the gear teeth crack, but the bending stiffness, shear stiffness, and axial compression stiffness will change with the change of the height. In order to deduce these three kinds of stiffness, we need to consider four cases, case 1 and 2 are that the crack does not pass through the centerline of the tooth, and case 3 and 4 are that the crack passes through the centerline of the tooth.  Figure 7 shows the crack schematic diagram of gear teeth when the root circle diameter is smaller than the base circle diameter. q 1 is the initial crack length, which does not reach the centerline of the tooth, and q 2 is that the crack size exceeds the centerline of the tooth. v is the angle between the crack line and the centerline of the tooth. With the continuous expansion of cracks, the cross-sectional area Ax and moment of inertia Ix can be refined according to the specific situations. There are four situations. Each situation will be discussed in detail one after another.
Situation 1: When ha ≥ h0 and α1 ≥ αa As in the case of healthy gear, the combined mesh stiffness consists of Hertz stiffness, bending stiffness, shear stiffness, and axial compression stiffness. The Hertz stiffness does not change when the gear teeth crack, but the bending stiffness, shear stiffness, and axial compression stiffness will change with the change of the height. In order to deduce these three kinds of stiffness, we need to consider four cases, case 1 and 2 are that the crack does not pass through the centerline of the tooth, and case 3 and 4 are that the crack passes through the centerline of the tooth.
With the continuous expansion of cracks, the cross-sectional area Ax and moment of inertia I x can be refined according to the specific situations. There are four situations. Each situation will be discussed in detail one after another.
Situation 1: When h a ≥ h 0 and α 1 ≥ α a Situation 2: When h a < h 0 or when h a ≥ h 0 and α 1 ≤ α a where h a = R b sin α 2 − q 1 sin ν Situation 3: When h c < h 0 or when h c ≥ h 0 and α 1 ≤ α c Situation 4: When h c < h 0 and α 1 ≥ α c According to the derivation method of TVMS, when the root circle of the healthy gear tooth is smaller than the base circle, the TVMS value under crack state can be derived. Figure 8 shows the crack schematic diagram of gear teeth when the base circle diameter is smaller than the root circle diameter.

The Diameter of the Root Circle Is Larger than That of the Base Circle
According to the derivation method of TVMS, when the root circle of the healthy gear tooth is smaller than the base circle, the TVMS value under crack state can be derived.
2.2.2. The Diameter of the Root Circle Is Larger than That of the Base Circle Figure 8 shows the crack schematic diagram of gear teeth when the base circle diameter is smaller than the root circle diameter. With the continuous expansion of cracks, the cross-sectional area Ax and moment of inertia Ix can be refined according to the specific situations. There are four situations. Each situation will be discussed in detail one after another.
Situation 1: When ha ≥ h0 and α1 ≥ αa Situation 4: When hc < h0 and α1 ≥ αc With the continuous expansion of cracks, the cross-sectional area A x and moment of inertia I x can be refined according to the specific situations. There are four situations. Each situation will be discussed in detail one after another.
Situation 1: When h a ≥ h 0 and α 1 ≥ α a Situation 2: When h a < h 0 or when h a ≥ h 0 and α 1 ≤ α a where h a = R b sin α 2 − q 1 sin ν Situation 3: When h c < h 0 or when h c ≥ h 0 and α 1 ≤ α c Situation 4: When h c < h 0 and α 1 ≥ α c According to the derivation method of TVMS when the root circle of a healthy tooth is larger than the base circle, the TVMS value under crack state can be derived.

Dynamic Equation Modeling of the Planetary Gear Train
A pure torsional dynamic model of planetary gear train is established considering nonlinear factors such as TVMS and meshing clearance. The structural diagram of the 2K-H planetary gear set is shown in Figure 9 and the pure torsional vibration model is shown in Figure 10.
( ) According to the derivation method of TVMS when the root circle of a healthy tooth is larger than the base circle, the TVMS value under crack state can be derived.

Dynamic Equation Modeling of the Planetary Gear Train
A pure torsional dynamic model of planetary gear train is established considering nonlinear factors such as TVMS and meshing clearance. The structural diagram of the 2K-H planetary gear set is shown in Figure 9 and the pure torsional vibration model is shown in Figure 10.  The connecting shaft of the sun gear is the power input shaft and the cage where the planetary gear is located in the output. Assuming that the number of sun gear is Zs, the number of teeth of planetary gear is Zp, the number of teeth of inner gear ring is Zr, and the transmission ratio of the first stage 2K-H planetary gear train is isc.
The planetary gear set is different from the traditional fixed shaft gear, there are many meshing points in the transmission process, and the meshing position is constantly changing. The planetary gear has meshed with the ring and the sun gear simultaneously, and there is mutual coupling and restraining effect in the meshing process, so that the mesh frequency component may be affected. Due to the rotation and revolution of the planetary gear, the vibration impact direction caused by ( ) According to the derivation method of TVMS when the root circle of a healthy tooth is larger than the base circle, the TVMS value under crack state can be derived.

Dynamic Equation Modeling of the Planetary Gear Train
A pure torsional dynamic model of planetary gear train is established considering nonlinear factors such as TVMS and meshing clearance. The structural diagram of the 2K-H planetary gear set is shown in Figure 9 and the pure torsional vibration model is shown in Figure 10.  The connecting shaft of the sun gear is the power input shaft and the cage where the planetary gear is located in the output. Assuming that the number of sun gear is Zs, the number of teeth of planetary gear is Zp, the number of teeth of inner gear ring is Zr, and the transmission ratio of the first stage 2K-H planetary gear train is isc.
The planetary gear set is different from the traditional fixed shaft gear, there are many meshing points in the transmission process, and the meshing position is constantly changing. The planetary gear has meshed with the ring and the sun gear simultaneously, and there is mutual coupling and restraining effect in the meshing process, so that the mesh frequency component may be affected. Due to the rotation and revolution of the planetary gear, the vibration impact direction caused by The connecting shaft of the sun gear is the power input shaft and the cage where the planetary gear is located in the output. Assuming that the number of sun gear is Z s , the number of teeth of planetary gear is Z p , the number of teeth of inner gear ring is Z r , and the transmission ratio of the first stage 2K-H planetary gear train is i sc .
The planetary gear set is different from the traditional fixed shaft gear, there are many meshing points in the transmission process, and the meshing position is constantly changing. The planetary gear has meshed with the ring and the sun gear simultaneously, and there is mutual coupling and restraining effect in the meshing process, so that the mesh frequency component may be affected. Due to the rotation and revolution of the planetary gear, the vibration impact direction caused by planetary gear damage will change constantly. Because the position of the sensor on the planetary gear housing is fixed, the angle between the impact direction and the sensor will constantly change, which makes the fault frequency more difficult to obtain. These characteristics make the traditional signal processing and pattern recognition methods not competent for planetary gear train fault diagnosis and prediction.
For fault mechanism, the dynamic model of the planetary gear train is established to study the influence of damage mode on dynamic response results.
The dynamic model includes a lumped parameter model and a distributed parameter model. The lumped parameter model can also be called a multi-degree of freedom spring-mass vibration system. There are two types of coupled planetary vibration models: the pure torsion model and translational torsion model. Among them, the core parameters of the pure torsion model include circumferential torsional stiffness, the moment of inertia of center construction, and mass of planetary gear. The pure torsion model has the advantages of fewer parameters, simpler model and less computation. The translational torsional coupled vibration model considers not only torsion but also translational motion in two orthogonal directions in the vertical plane of the axis. The translational torsional coupled vibration model has more parameters. Therefore, the model is more complex and challenging to solve. When the support stiffness is 10 times or more than the mesh stiffness, the pure torsional stiffness model can also achieve high calculation accuracy [26].
The meshing force of gear is approximately elastic meshing force P and viscous meshing force D. Where, K spi and K rpi denote the mesh stiffness between the sun gear and the i-th planet gear, the ring gear, and the i-th planet gear respectively. C spi and C rpi represent the meshing damping.  (x, b) is a nonlinear clearance function, which can be specifically expressed as follows: where b is the gap constant. The dynamic simulation model of the planetary gear set can be shown as follows: The above equation is a system of positive semidefinite nonlinear second order differential equations with N + 2 degrees of freedom. The angular displacement of the rigid body in the system is transformed into relative linear displacement. The relative displacements between the sun gear and the planet gear and between the sun gear and the carrier are defined as follows: x spi = θ s r bs − θ pi r bpi − θ c r c cos α x sc = θ s r bs − 2θ c r c cos α Substituting Formulas (52) and (54) into the Formula (51), we can get the following results: Substituting Formulas (54) and (55) into the Formula (53), we can get the following results: Table 1 shows the parameters of the planetary gearbox. Obviously, the diameter of the root circle of both the sun gear and the planet gear is smaller than the diameter of the base circle. Based on the above-mentioned improved method in Section 2, the TVMS of a pair of the sun-planet gears is calculated and shown in Figure 11 for the different health conditions. It shows the period that was experienced from the crack tooth starting to mesh to the ending of mesh. As the size of the crack grows, the mesh stiffness reduces correspondingly, which is consistent with the change of actual gear mesh stiffness. Correct stiffness calculation is an integral part of gear fault diagnosis. For each crack level, the corresponding vibration signal can be obtained through simulation calculation. By analyzing these vibration signals, the index of fault severity can be obtained, and then the fault severity of the gear system can be detected.

Simulation Analysis
The fourth-order Runge-Kutta method is used to numerically simulate the sun planetary vibration signal of the planetary gear in this section. Before obtaining the dynamic response, first to calculate the corresponding fault characteristic information under the set working conditions to

Simulation Analysis
The fourth-order Runge-Kutta method is used to numerically simulate the sun planetary vibration signal of the planetary gear in this section. Before obtaining the dynamic response, first to calculate the corresponding fault characteristic information under the set working conditions to verify the correctness of the solution stiffness and the establishment of the model. The following is the calculation process of fault characteristic information: Sun gear rotation frequency: Planetary carrier rotation frequency: Planetary gear rotation frequency: Gear mesh frequency: Among them, Z 1 is the number of sun gear teeth, Z 2 is the number of planetary gear teeth, Z 3 is the number of ring gear teeth, and n i is the input speed, that is, the sun gear speed. When n i = 1200 r/min, the calculated sun gear rotation frequency f s = 20 Hz, the rotation frequency of planet carrier is The relevant parameters of planetary gear system dynamics modeling and simulation are shown in the following Table 1. The main purpose of this paper is to verify the effectiveness of the TVMS solution method and the establishment of the dynamic model. To verify the correctness of the model, the vibration response of the planetary gear and sun gear under the healthy state and the cracked state is solved. The timedomain waveforms of the TVMS of the sun gear in healthy and faulty states are obtained, as shown in the following Figure 12. The main purpose of this paper is to verify the effectiveness of the TVMS solution method and the establishment of the dynamic model. To verify the correctness of the model, the vibration response of the planetary gear and sun gear under the healthy state and the cracked state is solved. The time-domain waveforms of the TVMS of the sun gear in healthy and faulty states are obtained, as shown in the following Figure 12. It can be seen from Figure 12 that as the degree of failure increases, the TVMS of the gear will gradually decrease. The mesh stiffness values in health and fault state were substituted into the simulation model, and the dynamic response results were shown in Figures 13 and 14. The envelope spectrum of the vibration signal is shown in Figures 15 and 16. By comparing Figures 13 and 14, it can be found that when the sun gear cracks, periodic shock signals appear in the vibration signal. Figure 15 shows the envelope spectrum in a healthy condition. The gear mesh frequency fm = 437.5 Hz and its multiplication can be clearly found. Figure 16 shows the envelope spectrum in a crack condition, from which the mesh frequency fm = 437.5 Hz and its multiplication can be found clearly. Simultaneously, a large number of sidebands appear near the mesh frequency in the vibration signal spectrum. In order to facilitate the analysis, the low-frequency band was amplified, and it was found that the sun gear fault characteristic It can be seen from Figure 12 that as the degree of failure increases, the TVMS of the gear will gradually decrease. The mesh stiffness values in health and fault state were substituted into the simulation model, and the dynamic response results were shown in Figures 13 and 14. The envelope spectrum of the vibration signal is shown in Figures 15 and 16. By comparing Figures 13 and 14, it can be found that when the sun gear cracks, periodic shock signals appear in the vibration signal. Figure 15 shows the envelope spectrum in a healthy condition. The gear mesh frequency f m = 437.5 Hz and its multiplication can be clearly found.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 16 of 21 frequency (15.63 Hz) and its multiplication frequency existed (31.26 Hz and 46.9 Hz). The spectrum structure obtained by simulation analysis is also consistent with the theoretical analysis results.

Experimental Analysis
The fault simulation tests were carried out on the two stage planetary gearbox in healthy and crack condition of the sun gear. The fault simulation test bench is shown in Figure 17. The planetary gearbox contains a single planetary row, the sun gear is the input, and the planet carrier is the output. The specific parameters of the test bench are shown in Table 1. The gears used in the test are all standard gears. In addition, we also assume that the gear has no manufacturing errors. In the test, the vibration sensor is install in the vertical direction of the case of the planetary gearbox. The input shaft speed of the gearbox, that is, the speed of the sun gear is 1200 r/min. The sampling frequency was 20,480 Hz. Under the same load conditions, the vibration signal data of the planetary gearbox were collected in the healthy and crack state of the sun gear. The root crack implantation method is wire cut electrical discharge machining (WEDM) in the experiments.  Figure 16 shows the envelope spectrum in a crack condition, from which the mesh frequency f m = 437.5 Hz and its multiplication can be found clearly. Simultaneously, a large number of sidebands appear near the mesh frequency in the vibration signal spectrum. In order to facilitate the analysis, the low-frequency band was amplified, and it was found that the sun gear fault characteristic frequency (15.63 Hz) and its multiplication frequency existed (31.26 Hz and 46.9 Hz). The spectrum structure obtained by simulation analysis is also consistent with the theoretical analysis results.

Experimental Analysis
The fault simulation tests were carried out on the two stage planetary gearbox in healthy and crack condition of the sun gear. The fault simulation test bench is shown in Figure 17. The planetary gearbox contains a single planetary row, the sun gear is the input, and the planet carrier is the output. The specific parameters of the test bench are shown in Table 1. The gears used in the test are all standard gears. In addition, we also assume that the gear has no manufacturing errors. In the test, the vibration sensor is install in the vertical direction of the case of the planetary gearbox. The input shaft speed of the gearbox, that is, the speed of the sun gear is 1200 r/min. The sampling frequency was 20,480 Hz. Under the same load conditions, the vibration signal data of the planetary gearbox were collected in the healthy and crack state of the sun gear. The root crack implantation method is wire cut electrical discharge machining (WEDM) in the experiments.

Experimental Analysis
The fault simulation tests were carried out on the two stage planetary gearbox in healthy and crack condition of the sun gear. The fault simulation test bench is shown in Figure 17. The planetary gearbox contains a single planetary row, the sun gear is the input, and the planet carrier is the output. The specific parameters of the test bench are shown in Table 1. The gears used in the test are all standard gears. In addition, we also assume that the gear has no manufacturing errors. In the test, the vibration sensor is install in the vertical direction of the case of the planetary gearbox. The input shaft speed of the gearbox, that is, the speed of the sun gear is 1200 r/min. The sampling frequency was 20,480 Hz. Under the same load conditions, the vibration signal data of the planetary gearbox were collected in the healthy and crack state of the sun gear. The root crack implantation method is wire cut electrical discharge machining (WEDM) in the experiments. The voltage signal obtained based on the sun gear acceleration sensor is shown in Figure 18. The mesh frequency of planetary gears can be easily found as well in Figure 19. Figure 20 shows the timedomain waveform of the experimental data of the crack condition of the sun gear. In the case of sun gear tooth cracks, there will be some spikes in the vibration signal. In addition, the effect of amplitude modulation is not obvious, because amplitude modulation and fault signal together cause the complexity of the vibration signal. Similar to the simulation signal verification, the envelope The voltage signal obtained based on the sun gear acceleration sensor is shown in Figure 18. The mesh frequency of planetary gears can be easily found as well in Figure 19. Figure 20 shows the time-domain waveform of the experimental data of the crack condition of the sun gear. In the case of sun gear tooth cracks, there will be some spikes in the vibration signal. In addition, the effect of amplitude modulation is not obvious, because amplitude modulation and fault signal together cause the complexity of the vibration signal. Similar to the simulation signal verification, the envelope spectrum analysis of the acquired time-domain signal is performed, and the obtained results are shown in the following Figure 21: Appl. Sci. 2020, 10, x FOR PEER REVIEW 18 of 21 spectrum analysis of the acquired time-domain signal is performed, and the obtained results are shown in the following Figure 21: Figure 18. Test signal in the healthy condition. Figure 18. Test signal in the healthy condition. Figure 18. Test signal in the healthy condition.    The corresponding gear mesh frequency fm = 437.5 Hz and its multiplier can be found from the spectrogram of the sun gear tooth root crack' experimental signal. In the low-frequency band, the corresponding sun gear fault characteristic frequency 15.63 Hz can also be found clearly. Therefore, through the analysis of the experimental signal, the correctness of the TVMS solution method and the