Dynamic behavior analysis method for the normal and faulty drivetrain of wind turbine under multi‐working conditions

Due to the random fluctuation of wind energy, the working conditions of the drivetrain have complex time‐varying characteristics. To deeply understand time‐varying dynamic behavior, a dynamic behavior method for normal and faulty drivetrain under multiple operating status is proposed. First, the working condition is divided into multiple working conditions by combining the control principle of the wind turbine with the K‐means clustering algorithm. Second, the torsional dynamics model of the drivetrain is established, solved, and verified. Finally, the dynamic behavior of normal and faulty drivetrain under multiple working conditions is analyzed. The dynamic behavior of drivetrain in normal and faulty states under different working conditions can be determined by the method proposed in this paper, which is affected by the system itself, external load, and the control strategy of the wind turbine. At the same time, the fault‐sensitive working condition range and the fault feature index under each working condition are determined. This research can provide an important theoretical basis for variable condition fault diagnosis.


| INTRODUCTION
The drivetrain of the wind turbine, which is responsible for energy and load transfer, is always in a complex, changeable, and harsh working environment. 1Due to the wide operating speed range, wide excitation frequency bandwidth, and the coupling effect of inertial force and aerodynamic load, the failure rate of drivetrain increases, accounting for 40%-60% of the wind turbine downtime.It has a serious impact on the efficiency and economic benefits of wind turbines. 2At present, wind turbines use fixed thresholds to monitor the operating status, which cannot meet the requirements of condition monitoring due to the complexity and diversity of the operating conditions of the drivetrain.Therefore, condition assessment and fault diagnosis need to consider multiple operating conditions, and set corresponding thresholds according to different operating conditions.Before that, it is necessary to study the dynamic behavior of the drivetrain under complex and variable working conditions.
Many researchers have studied the dynamic behavior of the wind turbine drivetrain.Tan et al. 3 proposed a rigid-flexible coupling dynamic modeling method and investigated the modal shape, modal energy, and amplitude-frequency responses of the drivetrain considering the electromagnetic stiffness.Dewangan et al. 4 developed a dynamic model to implement the amplitude modulation effect due to carrier rotation and the effect of gravity, and verified it by experiments.Li et al. 5 established an integrated drivetrain coupling analysis model that showed that the non-torque loads aggravate the nonuniform load sharing between planet gears, and the non-torque loads enlarge the bearing dynamic supporting forces.Park et al. 6 analyzed the multibody dynamics of the drivetrain considering the shaft bending effect and a variable gear mesh, and pointed out that the planetary gear stages are more sensitive to shaft bending and eccentricity.Yang et al. 7 analyzed the nonlinear dynamic response of cyclic load and random wind load on the drivetrain, proposed that the planetary gear can get rid of chaos and enter into stable periodic motion by changing the stiffness ratio properly.Guo et al. 8 developed the integrated finite element model of a 5-MW wind turbine and studied the load transmission along its transmitting route and mechanical energy distribution during dynamic response under random wind loads.Sun et al. 9 established the translation-torsion dynamic analysis model considering the tooth shape error factor, and analyzed the influence of each order error on the dynamic characteristics of the system.Zhao et al. 10 established the full flexible multibody dynamic model of the transmission system, studied the effect of the structural flexibility on the dynamic response, and pointed out that the housing flexibility is the most significant factor affecting the dynamic response.Wu et al. 11 established drivetrain dynamic model with gear ring flexibility and analyzed the influence of ring flexibility, number of ring support points, and load torque on ring deformation and meshing force.Xiang et al. 12 analyzed the influence of external load excitation frequency and meshing damping ratio on the vibration response and bifurcation characteristics of the drivetrain.These studies analyzed the dynamics of the normal drivetrain considering electromagnetic stiffness, shaft bending effect, time-varying meshing stiffness, and so on, but did not consider the operating condition of the drivetrain.Some researchers have studied the dynamics of the faulty drivetrain.Zhang et al. 13 studied the effect of wear fault on the dynamic characteristics of wind turbine drivetrain considering tooth contact temperature and showed that the wear fault under the influence of tooth contact temperature will make the chaotic motion occur in advance.He et al. 14 established a rigid-flexible coupling dynamics model to analyze the vibration characteristics of the floating sun gear and deduced the mechanism of modulating sidebands indicating different kinds of defects.Dadona et al. 15 considered the problem of partial contact loss along the tooth contact line due to the fault's geometry and proposed a method to simulate partial tooth face faults and validated it by experiments.Wang et al. 16 investigated the effect of variability and uncertainty of wind and wave conditions on the short-term fatigue damage of wind turbine drivetrain and pointed out that wind has a large effect on the fatigue damage while waves are relatively small.Lei et al. 17 established a dynamic model of healthy, cracked, and spalled drivetrain considering the vibration transmission path, and analyzed its dynamic characteristics.Doğan et al. 18 proposed a numerical fault detection model based on dynamic transmission error, and studied the influence of crack length on gear meshing stiffness.Bi et al. 19 analyzed the influence of different tooth root crack growth degrees on the vibration characteristics of highspeed gear system, and identified the tooth root crack failure mode by using the generalized back propagation neural network.These studies analyzed the dynamics of the faulty drivetrain considering cracks, wear, tooth spalling, and so on, but did not consider the fault dynamics under various operating conditions.
Therefore, some researchers have studied the dynamic behavior of the drivetrain under different working conditions.Qin et al. 20 established a dynamic external load calculation model to analyze the input load variation pattern and influencing factors of wind turbine gearbox under start-up, operation, and braking conditions.Tang et al. 21discretized the input torque curve into 10 torque conditions, and analyzed the transmission characteristics and dynamic response of the wind turbine gearbox under different working conditions.Boukhezzar et al. 22 divided the operating status of wind turbine into above-rated wind speed and below-rated wind speed, two working regions, and compared linear and nonlinear control strategies based on nonlinear dynamics with the accompanying wind speed.Rezamand et al. 23 employed the kernel fuzzy C-means to categorize different operating conditions and analyzed the dynamics of wind turbine drivetrain bearings to predict the remaining useful life of faulty bearings.Chen et al. 24 established a mechanical-electrical coupling model of wind turbine drivetrain with variable speed and load, and revealed the relationship between dynamic characteristics and electromagnetic characteristics under different operating conditions.Tian et al. 25 analyzed the vibration time-frequency characteristics under different wear severity and torque conditions, and verified it by experiments.However, the classification of working conditions in these studies is relatively simple, and it is difficult to accurately describe the dynamics of wind turbines under time-varying working conditions.Also, there are fewer studies on the multi-working condition dynamics characteristics of faulty drivetrain.
In this paper, a dynamic behavior analysis method for normal and faulty drivetrain under multiple working conditions is proposed, which can provide a theoretical basis for condition assessment and variable condition fault diagnosis.The overall flowchart is illustrated in Figure 1.The remainder is organized as follows.In Section 2, the working condition of the drivetrain is divided and the variation of working conditions feature parameters under different working conditions is analyzed.In Section 3, the torsional dynamics model considering external excitation and time-varying meshing stiffness is established.In the last section, the dynamics model for normal and fault drivetrain is solved and verified, and its dynamic behavior under multiple working conditions is analyzed.

| WORKING CONDITION CLASSIFICATION OF DRIVETRAIN
The method combining the wind turbine control principle and K-means clustering is used to classify the wind turbine operating status. 26First, the historical SCADA (supervisory control and data acquisition) data of wind turbines are cleaned to remove abnormal values such as power curtailment, outlier, and so on. 27Then, the operation status is divided into five stages (shutdown, start-up, maximum wind-energy tracking, constant speed, and constant power) by the wind turbine control principle.Considering the influence of environmental factors on the operation status of the wind turbine, the active power, 30-s average wind speed, hub speed, and nacelle temperature of SCADA parameters are selected as feature parameters of working conditions.Due to the maximum wind-energy tracking stage and constant speed stage accounting for a large proportion, these two stages are subdivided by K-means clustering according to the selected parameters.The Calinski-Harabas criterion is adopted to determine the number of clusters k, and the k for both stages is 4. Therefore, the working condition of the drivetrain is finally classified into nine categories (Figure 2).This working condition classification method classifies the drivetrain operating status as more complex through secondary division.This provides more accurate working conditions for condition evaluation, fault diagnosis, and other work.
The variation of the working conditions feature parameters under different working conditions is analyzed to clarify the relationship between the parameters and multiple working conditions.For the feature parameters of each working condition, 1000 sample points are selected for comparison and analysis (Figure 3).The 30-s average wind speed is directly related to the external excitation, and the wind speed ranges are clearly distinguished under different working conditions.The nacelle temperature is related to the external environment of the drivetrain and is weakly correlated to the working conditions compared to the other three parameters.The hub speed and active power are directly related to the working conditions.Working condition-1 (shutdown stage) includes the part where the unit stops and the hub starts rotating but does not start generating power, so it has nonzero values in it.The range of each parameter at different working conditions is shown in Table 1.
The working condition classification.Active power versus 30-s average wind speed and hub speed.

| Torsional dynamics model
To analyze the internal load transfer and dynamic behavior under multiple working conditions, the dynamics model of the drivetrain must be established first.Kahraman et al. proposed in 1944 that when the support stiffness is 10 times or more than the meshing stiffness, the torsion model can achieve higher accuracy. 28he gearbox consists of one-stage planetary gear and two-stage fixed shaft gears.The planetary gear is the lowspeed stage with planetary carrier input, sun gear output, and fixed ring gear.The two-stage fixed shaft gears are medium-and high-speed stages.The wind energy is transferred to the planetary carrier through the impeller, distributed by four planetary gears, then transferred from the sun gear to the medium-speed stage gear pair, and finally output by the high-speed stage.The torsional dynamics model established by the lumped mass method is shown in Figure 4.
The angular displacement and torque of each member are specified to be positive counterclockwise.The derived equations for the torsional dynamics of the drivetrain are as follows: n b c pi s = , , , , 1, 2, 3, and 4; m p is planetary gear mass; r c is the distance from the center of the planetary frame to the center of the planetary gear; r bi is base circle radius of each gear; K m , K l , K h and C m , C l , and C h are torsional stiffness and damping of the shaft at all levels; F spi , F rpi , F 12 , and F 34 are the meshing force of each gear pair; α sp is the angle of meshing between the sun gear and planetary gear; α rp is the angle of meshing between planetary gear and gear ring; T in and T out are input and output torque.
The relative gear mesh displacement along the action line of each gear pair is as follows: (2) The meshing force of each gear pair is as follows: where K t and C t are the time-varying meshing stiffness and damping of each gear pair, t spi rpi = , , 1, 2, 3, and 4.

| Stiffness and damping
Gear meshing stiffness is one of the main internal excitations of gear dynamics, and accurate evaluation of meshing stiffness is the prerequisite of gear dynamics analysis.When using the potential energy method, which simplified the gear tooth as a cantilever beam on the base circle to solve the meshing stiffness, the noncoincidence between the base circle and the root circle will lead to calculation errors.Therefore, the time-varying meshing stiffness is calculated based on the improved potential energy method, which simplified the gear tooth as a cantilever beam on the root circle to reduce calculation errors. 29n this paper, the root crack fault of the sun gear is taken as a case to study the dynamic behavior of the faulty drivetrain under multi-working conditions, thus the influence of the crack on meshing stiffness is considered.The sun gear tooth models of normal tooth and cracked tooth are shown in Figure 7.The total energy stored in the paired meshing gears includes hertz energy, bending energy, shear energy, and axial compression energy, which correspond to hertz stiffness, bending stiffness, shear stiffness, and radial compression stiffness respectively.The meshing stiffness during single tooth meshing can be expressed as When two pairs of gear teeth are engaged at the same time, the meshing stiffness can be expressed as 29 where k h is the hertzian stiffness, k b is the bending stiffness, k s is the shear stiffness, k a is the radial compression stiffness, 1 and 2 are the active and driven gears, and i indicates which pair of gear teeth is engaged.The formula for each meshing stiffness component (hertz stiffness, bending stiffness, shear stiffness, and radial compression stiffness) and the explanation of the symbols in Figure 5 are shown in the appendix.
The angle φ between the crack line and tooth center line is assumed to be constant and is set as 45°.This paper analyzed five levels of crack growth: 10%, 30%, and 50%.The level of crack growth is defined as The time-varying mesh stiffness of normal and faulty sun-planet pairs is shown in Figure 6.The X-axis is the rotation angular displacement of the planet gear θ p .It can be observed that as the size of the crack grows, the mesh stiffness reduces correspondingly.
Meshing damping is as follows 12 : where ξ is the damping ratio, generally taken as 0.005-0.075;m 1 and m 2 are the masses of the active and driven gears.

| External excitation
The wind speed data collected by the SCADA system is used as the actual wind speed, and the induced timevarying load is the external excitation.According to the aerodynamic Betz theory, the aerodynamic power of the wind turbine impeller is mainly related to the wind speed and rotor power coefficient 30 where P is aerodynamic power, ρ is air density, R is impeller radius, V w is wind speed, and C p is rotor power coefficient.
Rotor power coefficient C p , which describes wind energy utilization rate of wind turbine, is adjusted by the pitch angle, so the rotor power coefficient can be regarded as a nonlinear function of the pitch angle and speed ratio 30 where c 1 -c 9 are wind turbine parameter, β is pitch angle, γ is intermediate variable, λ is speed ratio, and ω is hub speed.
Input torque T in can be calculated from aerodynamic power P and hub speed ω as in (10)

| Model solution
In this paper, the 2-MW-WT2000D110H80 double-fed wind turbine of a wind farm is used as the research object.The main design parameters of the wind turbine are shown in Table 2.The gearbox parameters of the wind turbine main drivetrain are shown in Table 3.
The dynamic behavior under different working conditions is studied based on the normal operation data of wind turbine on January 4, 2021.First, the working conditions are identified based on historical SCADA data (Figure 7).Then, the stiffness, damping, and external excitation of the torsional dynamics model are calculated based on the drivetrain parameters.Finally, the fourth-order Runge-Kutta method is used to solve the dynamics equations of the drivetrain.

| Dynamics model validation
To verify the accuracy of the torsional dynamics model, a three-axis vibration acceleration sensor is installed on the The mesh stiffness beam model.(A) Normal gear.(B) Cracked gear.
F I G U R E 6 Time-varying meshing stiffness.
HUANG ET AL.
| 3625 housing at the input end of the gearbox (Figure 8).The Xdirection is the spindle axis, and the sampling rate is 500 Hz.The operating data of the sun gear under working condition-9 in the constant power stage is selected for verification.Spectral analysis is conducted on the output signal of the dynamic model (torsional acceleration of sun gear, and relative gear mesh displacement between the sun gear and planetary gear) and the real collected signal (Y-direction acceleration), as shown in Figure 9.
In the selected operating data, the maximum value of the hub speed is 13.62 rpm.Therefore, the theoretical rotational frequency of the sun gear is 1.324 Hz, and the theoretical meshing frequency of the sun gear and planetary gear is 19.59 Hz.From Figure 11, it can be seen that the rotational frequency and meshing frequency of the dynamics model output signal are 1.367 and 19.59 Hz, respectively, and the real rotational frequency and actual meshing frequency are 1.333 and 19.98 Hz, respectively.Compared the output signal of the dynamics model with the real collected signal, the error of the sun gear rotational frequency and the sun gearplanetary gear meshing frequency are within a reasonable range.
Meanwhile, the torque feedback value in SCADA data describes the gearbox output torque.Comparing the output torque from the torsional dynamics model and SCADA data (Figure 10), the overall trend is relatively consistent.

| Dynamic behavior of normal drivetrain under multi-working conditions
With the increase in wind speed, the input load keeps increasing.At the same time, the wind turbine controls the input load variation of the drivetrain by controlling the electromagnetic torque and adjusting the pitch angle according to the control strategy, so that the dynamic behavior of the drivetrain components has different characteristics under different working conditions.The variation law of relative gear meshing displacement and meshing frequency of each gear pair under each working condition at other stages is shown in Table 4.
Under working condition-2, the wind turbine starts to operate.The pitch angle is quickly optimized, and the input load is determined by the wind speed.The speed and load of each component of the drivetrain increase rapidly.
The maximum wind-energy tracking phase aims to track the rotor power coefficient.The electromagnetic torque of the generator is controlled to adjust the hub speed to follow wind speed.The input load is basically determined by the wind speed, but is influenced by the generator's electromagnetic torque control.The rotor power coefficient reaches the maximum value rapidly under working condition-3, keeps fluctuating at the maximum value under working condition-4, and the input load and hub speed increases.Under working condition-5, as the hub speed reaches the rated speed, the wind turbine controls the generator's electromagnetic torque to increase.As a result, the rotor power coefficient fluctuates sharply at the maximum value and the input load is affected.
In the constant speed stage, the hub speed reaches the rated speed.The wind turbine controls the generator's electromagnetic torque to absorb energy exceeding the rated wind speed and achieves constant hub speed operation.The input load is determined by the wind speed in conjunction with the generator's electromagnetic torque control and is influenced by the pitch angle control.From working condition-6 to working condition-8, the hub speed remains at the rated speed and is stable, but the rotor power coefficient keeps fluctuating, and the input load increases steadily.Meanwhile, under working condition-8, as the power reaches the rated power, the wind turbine starts to adjust the pitch angle, resulting in an increase in the fluctuation of the rotor power coefficient and an impact on the input load.
In the constant power phase, under working condition-9, the wind turbine continuously adjusts the pitch angle to control the power to maintain the rated power output.The input load is determined by the paddle pitch angle control and remains essentially constant.
Taking the sun gear as an example, the frequency domain diagram of the relative gear mesh displacement between the sun gear and the planetary gear under | 3627 different working conditions is shown in Figure 11.The frequency components, which are mainly composed of meshing frequency and its harmonic frequency, increase continuously during the start-up phase and the maximum wind-energy tracking phase.In the constant speed phase, they reach and fluctuate around the maximum value.The variation law is consistent with Table 4.

| Dynamic behavior of faulty drivetrain under multi-working conditions
The multi-working condition dynamics analysis of the faulty drivetrain can clarify the fault-sensitive working conditions and quantify the fault under different The displacement increases rapidly, but with a small amplitude.
The frequency amplitude is low and the range is small.

Maximum windenergy tracking phase 3
The displacement amplitude and range increase.
The frequency amplitude and range increase rapidly.

4
The displacement amplitude increases and the range increases rapidly.
The frequency amplitude and range increases.

5
The displacement amplitude and range increase continuously.
The frequency amplitude increases but the range decreases.

Constant speed phase 6
The displacement amplitude increases but the growth rate decreases and the range increases slightly.
The frequency amplitude increases slightly, and the range is large compared to other working conditions in this stage. 7 The displacement amplitude and range increase insignificantly.
The frequency amplitude remains essentially the same but the range is reduced. 8 The displacement amplitude increases and the growth rate improves, and the range increases rapidly.
The frequency amplitude remains essentially the same and the range is reduced again.

Constant power phase 9
The displacement amplitude no longer increases and the range decreases.
The frequency amplitude remains essentially constant and the range is minimal.
F I G U R E 11 Frequency domain diagram of the sun gear and planetary gear relative gear mesh displacements for each working condition.
conditions through the fault feature indexes, providing a theoretical basis for the fault diagnosis of the wind turbine drivetrain.In this section, the dynamic behavior with different crack levels under the same working condition, the fault-sensitive working conditions, and the fault feature indexes under different working conditions are analyzed, taking the sun gear tooth root crack fault as an example.

| Dynamic behavior with different crack levels under the same working conditions
Under the same working conditions, the dynamic behavior of drivetrain changes with the crack size.The time-domain and frequency-domain diagrams of the relative meshing displacement of sun-planetary gear for the normal and faulty states under working condition-9 are shown in Figure 12.In the normal state, the gear teeth are perfectly engaged and shock-free.However, when there is a root crack fault, the meshing of the cracked gear will produce larger displacement and impact.The fault feature appears and increases as the crack level rises.
In Figure 14A, the fault feature at the 30% crack level is slight but not obvious, while they are particularly pronounced at the 50% crack level.However, it is difficult to detect at the 10% crack level, which is also the reason why tooth root crack faults are difficult to detect at the early stage.The periodicity of the fault feature generation is related to the rotation speed of the sun gear.The cracked gear teeth mesh with each planetary gear once for each rotation of the sun gear.In Figure 14B, at each crack level, the fault feature frequency caused by the crack appears, and the sideband caused by the modulation of the fault feature frequency appears near the meshing frequency.With the increase of crack level, the amplitude of the fault feature frequency and the sideband increases continuously.Taking 50% crack as an example, the frequency domain of 0-100 Hz is shown in Figure 13.
The fault feature frequency of the sun gear tooth root crack fault can be calculated by where f crack is the fault feature frequency, f m is the meshing frequency, N is the number of planetary gears, and Z s is the number of sun gear teeth.
The theoretically calculated and simulated fault feature frequencies are 4.354 and 4.350 Hz, respectively, which are basically the same.When the fault occurs, the fault feature frequency and its harmonic frequency appear in the frequency domain.At the same time, the sideband caused by the fault feature frequency is generated near the meshing frequency.The sideband appears in the following locations:  In the time domain, the 50% crack-induced shock begins to appear at condition-4 and becomes more and more obvious with the increase in working conditions.In the frequency domain, the fault feature frequency can already be found at condition-2, and the amplitude continues to increase with the increase of condition.The fault feature frequency is related to the rotation speed of the sun gear.Similarly, it keeps increasing from condition-2 to condition-5.The growth rate keeps increasing from condition-2 to condition-4 and starts to decrease continuously at condition-5.And from condition-6 to condition-9, it fluctuates around the maximum and tends to be stable.The frequency comparison under different working conditions in normal state and 50% crack level faulty state is shown in Figure 15.
The shock caused by the 10% crack does not appear on the time domain in every working condition, but the fault feature frequency caused by the crack starts to appear at condition-4.The shock caused by 30% cracks begins to appear in the time domain at condition-7, and the fault feature frequency appears in the frequency domain at condition-3.

| Feature index analysis of tooth root crack fault under different working conditions
Taking eight typical feature indexes as examples, the changes of each feature index under different working conditions are analyzed to clarify the tooth root crack fault feature indexes under each working condition.The selected typical feature indexes and their equations are shown in Table 5.
To clarify the sensitivity of the feature indexes under different working conditions, the feature index percentage change (P F ) of fault state value relative to the normal state value is used to measure the change degree of the vibration signal.The equation is as follows: where F crack is the fault state value of the feature indexes and F is the normal state value of the feature indexes.The analysis of the feature indexes of different crack levels shows that the change degree of the feature indexes of each crack level varies under different working conditions.In general, the peak-to-peak value (F 3 ) has the best sensitivity among the eight feature indexes.Starting from the working condition-5, the clearance factor (F 5 ), crest factor (F 7 ), and impulse factor (F 8 ) also have good sensitivity.T A B L E 5 The selected typical feature indexes and their equations.

Label
Feature index Equation Label Feature index Equation Taking the working condition-9 as an example, the percentage changes of feature indexes for different crack levels are shown in Table 6.The percentage changes of the feature indexes all increase continuously with the growth of the crack.The percentage change of peak-topeak value (F 3 ) is large compared to other feature indexes, and the clearance factor (F 5 ), crest factor (F 7 ), and impulse factor (F 8 ) also have good sensitivity at a 50% crack level.
The variation of feature indexes for the different crack level faults coincides with their fault-sensitive working conditions.Taking the 50% crack level as an example, the percentage changes of each feature index under different working conditions are shown in Figure 16.Working condition-2 and condition-3 are affected by noise, and are difficult to be characterized by each feature index.The sensitivity of the peak-to-peak value (F 3 ) is well at condition-4, followed by the root mean square (F 2 ).The sensitivity of the clearance factor (F 5 ), crest factor (F 7 ), and impulse factor (F 8 ) are well at condition-5.The sensitivity of the peak-to-peak value (F 3 ) is well at conditions-6 to conditions-9, and the clearance factor (F 5 ), crest factor (F 7 ), and impulse factor (F 8 ) are second.

| CONCLUSIONS
In this paper, a dynamic characteristics analysis method of normal and faulty main drivetrain under multiple operating conditions is proposed.This research provides important theoretical values for structural design and fault diagnosis.Several conclusions are obtained: 1.The mechanism of correlation between working condition feature parameters and working conditions is clarified, and the load variation of the normal condition of the main drive system is revealed.The results show that the wind turbine regulates the input load variation according to the control strategy, so the drivetrain dynamic behavior has different variation characteristics under different working conditions.The input load of the drivetrain under working conditions 2-5 is determined by the wind speed, but working conditions 4 and 5 are influenced by the generator's electromagnetic torque control.The input load under working conditions 6-8 is determined by the wind speed combined with the electromagnetic torque control, and working condition-8 is influenced by the pitch angle control.The input load under working condition-9 is determined by the pitch angle control.2. The faulty sensitive working conditions for each crack level of the sun gear crack fault and the faulty feature indexes under each working condition are analyzed.
The results show that the peak-to-peak feature index has a good sensitivity to faults, and the variation of feature indexes for the different crack level faults coincides with their fault-sensitive working conditions.3. The proposed method can simulate more types of faults or multiple fault coupling, and clarify the faulty sensitive working condition range of each fault as well as the faulty feature index.To calculate the hertzian stiffness, the gear is assumed to be an isotropic elastomer and the elastic deformation of its contact region is approximated as a parabolic contact.Taking the planetary gear and gear ring mesh as an example, the planetary gear is equated to a cylinder and the gear ring is equated to a circular groove, as shown in Figure A1.The half-width of the contact zone is

AUTHOR CONTRIBUTIONS
where E is the elastic modulus, L is the tooth width, v is the Poisson ratio, and F is the contact force.The contact deformation is Expanding the quadratic term of Equation (A2) and taking the approximation Combining Equation (A4) with Equation (A1), the hertzian stiffness k h can be calculated as follows: For an external gear, if the gear is a standard spur gear with the pressure angle of 20°, the root circle is bigger if the tooth number is more than 41.It is smaller if the tooth number is less than 41.These two cases will be discussed separately.

Case 1:
The gear root circle is smaller than the base circle.
The gear tooth model for this case is shown in Figure A2, where the gear tooth profile follows an involute curve and the curve is simplified to the straight line (NN′, DD′).
The action force F which is along the action line can be decomposed into two orthogonal forces F r and According to the beam theory, the bending, shear, and axial compressive energies stored in a tooth can be calculated as Contact deformation the planetary gear-gear ring mesh.

F I G U R E A2
The gear tooth model when the root circle is smaller than the base circle.
where k b , k s , and k a are the bending, shear and axial compressive stiffness, d is the distance from the contact point to the gear root, h is the distance from the contact point to the symmetry line, G is shear modulus, A x and I x are the area and the area moment of inertia of the section where the distance to the tooth root is x.h, h x , d, I x , and A x can be expressed as follows: where α 1 is the contact angle, α 2 is the half tooth angle on the base circle, and α 3 is the approximated half tooth angle on the root circle.α 1 , α 2 , and α 3 can be expressed as follows: where θ is the angular displacement, Z 1 and Z 2 are the number of teeth of the meshing gear, and α 0 is the pressure angle.Substituting Equations (A6) and (A10) into (A7), (A8), and (A9), k b , k s , and k a can be expressed as  ( ) The gear root circle is bigger than the base circle.
The gear tooth model for this case is shown in Figure A3.Compared with Case 1, the equations of d and h x are changed and the half tooth angle on the root circle of Case 2 is α 4 The gear tooth model when the root circle is bigger than the base circle.
Similar to Case 1, the k b , k s , and k a of Case 2 can be expressed as where α 5 is the angle between F and F b when the contact point is on the root circle, which can be calculated by Equation (A21)

(A21)
For the single-tooth-pair meshing duration, the total mesh stiffness can be expressed as For the double-tooth-pair meshing duration, the total mesh stiffness can be expressed as where 1 and 2 are the active and driven gears and i indicates which pair of gear teeth is meshed.

Calculation of meshing stiffness when gears with cracked tooth.
The cracked tooth model is shown in Figure A4, assuming that the crack arises at the root of the tooth.Due to the continuous smooth and small curvature of the crack expansion path, the root crack is approximated as a straight line.The crack expands along an angle (v) to the tooth center line.
Also consider the two cases: the gear root circle is smaller or bigger than the base circle.
1. Case 1: The gear root circle is smaller than the base circle.
The gear tooth model for this case is shown in Figure A5.The starting point of the crack is M, expanding along the straight line to the center line intersection B and extending to point D. Due to the fact that hertz stiffness and compression stiffness are not affected by the crack, their expressions remain unchanged.However, the expressions of the bending stiffness and the shear stiffness will be affected because of the change in the tooth length and the tooth height caused by the crack.

Condition 1: When
The symbol h a represents the distance between the end of the crack and the gear tooth center line when the crack has not extended to the center line.The symbol h 0 represents half of the roof chordal tooth thickness.The angle α a corresponds to the force action point K.For the cracked tooth, the tooth section area and the area moment of inertia can be expressed as follows: where h a and h x can be expressed as follows: Substituting Equations (A24) and (A25) into (A7) and (A8), k b and k s can be expressed as  The tooth section area and the area moment of inertia can be expressed as follows: Substituting Equations (A30) and (A31) into (A7) and (A8), k b and k s can be expressed as | 3637 The symbol h c represents the distance between the end of the crack and the gear tooth center line when the crack has extended to the center line and in another direction.The angle α c corresponds to the force action point E.
The tooth section area and the area moment of inertia can be expressed as follows: where h c can be expressed as follows: Substituting Equations (A34) and (A35) into (A7) and (A8), k b and k s can be expressed as  In this condition, if the section height h x is smaller than h c , the corresponding tooth part can be ignored.If h x is bigger than h c , the expressions of the tooth section area and the area moment of inertia are the same as Equations (A34) and (A35).Therefore, k b and k s can be expressed as The cracked gear tooth model when the root circle is bigger than the base circle.The gear tooth model for this case is shown in Figure A6.Similar to Case 1, four conditions are considered in Case 2. In each condition, the expressions of the tooth section area and the area moment of inertia are the same as Case 1. But, the expression of the tooth height h x changes as < .

3
Comparison of feature parameters under different working conditions.(A) Thirty-second average wind speed.(B) Nacelle temperature.(C) Hub speed.(D) Active power.DRIVETRAIN T A B L E 2 Main design parameters of the wind turbine.−25 to +50 (normal temperature); −40 to +50 (low temperature) Operating environment temperature −15 to +40 (normal temperature); −30 to +40 (low temperature) <2500 (plain type); <4000 (highland type) T A B L E 3 The gearbox parameters of the wind turbine drivetrain.

9
Comparison of the output signal of dynamic model spectrum map with real collected signal spectrum map.(A) Torsional acceleration of sun gear.(B) Relative gear mesh displacement between the sun gear and planetary gear.(C) Y-direction acceleration.F I G U R E 10 Comparison of supervisory control and data acquisition torque feedback values with simulated output torque.HUANG ET AL.
c k (m and n are integers).

4. 4 . 2 |
Sensitive working condition analysis of tooth root cracking fault Sensitive working condition analysis can clarify which working condition has the more obvious fault feature and can improve the accuracy of fault diagnosis.Also taking (A) F I G U R E Time-domain and frequency-domain diagram of the relative gear meshing displacement for the normal and faulty states.(A) Time domain.(B) Frequency domain.HUANG ET AL. | 3629 F I G U R E 13 Frequency domain of the relative sun-planet mesh displacement at 50% crack level (0-100 Hz).U R E 14 Time domain and frequency domain of the relative sun-planet meshing displacement at 50% crack level under different working conditions.(A) Time domain.(B) Frequency 50% tooth root as example, time and frequency domain of the relative sun-planet meshing displacement under different working conditions are shown in Figure 14.

F
I G U R E 15 Frequency comparison of normal state and 50% crack level faulty state under different working conditions.(A) Normal state frequency domain diagram.(B) Fifty percent crack level fault state frequency domain diagram.

F
I G U R E A4The cracked tooth model.F I G U R E A5The cracked gear tooth model when root circle is smaller than base circle.

2 . 2 :
Case The gear root circle is bigger than the base circle.
procedure as in Case 1, the expression of the k b and k s for the four conditions are derived and listed as follows: 1 + )( − )cos cos sin + ( − )cos − sin . 1 + )( − )cos cos sin + ( − )cos − sin .
Range of each parameter at different working conditions.
Variation of each gear pair under different working conditions.
T A B L E 4 The percentage changes of feature indexes for different crack levels under working condition-9.
review and editing: Yuhao Huang, Huanguo Chen, and Juchuan Dai.Project administration: Huanguo Chen.Funding acquisition: Huanguo Chen.All authors have read and agreed to the published version of the manuscript.T A B L E 6