Wind Turbine Gearbox Diagnosis Based on Stator Current

Abstract

Early detection of faults in wind energy systems can reduce downtime, operating, and maintenance costs while increasing productivity. This paper proposes a method based on the analysis of generator stator current signals to detect faults in a wind turbine gearbox equipped with a doubly fed induction generator (DFIG). A localized parameter model was established to simulate the vibratory response of a two-stage gear system under healthy and faulty conditions. The simulation was performed in the MATLAB/Simulink environment. The results include a detailed analysis of the mechanical part of the gearbox, highlighting mesh stiffness, output speed, and accelerations. Additionally, the electrical part was evaluated based on the current supplied by the doubly fed induction generator. The results were presented in the case of healthy gears and in the presence of faults such as a broken or cracked tooth. Fast Fourier transform (FFT) analysis was employed to detect gear defects in the stator current signal. The presence of a crack or broken tooth in the gearbox induces modulation of the DFIG stator current signals according to the shaft frequencies corresponding to the faulty gear. These findings provide a preliminary basis for the detection and diagnosis of this type of failure.

1. Introduction

An intricate electromechanical device called a wind turbine (WT) system transforms wind energy into electrical energy [1]. Figure 1 is an illustration of the most common setup. The three-bladed main rotor is supported by the main bearing, which also supplies the planetary gear with torque. The primary rotor is attached to a plate that serves as the planetary gear’s input. Three planets make up the planetary gear, and the shafts connecting them to the plate. The planets exert torque on the sun by rolling across the stationary ring. The planetary gear’s output is the sun shaft. Moreover, the two-stage parallel gear is dived by the sun. The parallel gear contains three shafts: the fast shaft, which drives the generator, the intermediate shaft, which is connected to the solar shaft, and the slow shaft, which drives the solar shaft. The parallel gear has a countershaft installed inside of it [2].

By the end of 2019, more than 650 gigawatts of wind capacity were installed worldwide [3,4]. However, the rate at which turbine failures occur makes it difficult to increase investment in wind energy [5].

In 2011, a study was carried out to identify the most critical component of WT [6]. This study pointed out the gearbox as one of the critical components. This component is known to have high downtimes and maintenance costs due to its repair and maintenance procedures [7]. Although, the failure rate of WT gearboxes is relatively low (8.9%). Its downtime is as long as 25.4 days in one year [8]. As a result, a large percentage of the energy is lost due to a malfunction of the latter [9]. Gearboxes come in second place in the failure downtime due to their size and robust linkage with other components, making them more difficult to access, repair, or even replace [10]. The majority of gearbox failures are caused by gear and bearing failures. According to [11], gear failures made up around 59% of wind turbine failure modes. The common gear faults are caused by anomalies in the gear tooth, such as a chipped tooth, broken tooth, root crack, spalling, wear, pitting, and surface damage [12].

The list of faults described above can be detected by analyzing vibration signals. Vibration analysis is the most widely used approach for fault detection in rotating machinery, particularly in WTs [13,14]. However, the application of vibration analysis requires additional vibration sensors and data acquisition devices. These sensors and devices are inevitably subject to failure, which could cause additional problems to system reliability and additional operating and maintenance costs. All these reasons led operators and researchers to look for other methods that would be less expensive to install [15]. It is noticed that a WT is an electromechanical system where the coupling between the generator and the damaged component in the gearbox generates vibrations induced by a gearbox fault which modulates the electrical signals measured from the generator terminals. As a result, gearbox faults can be diagnosed using generator current signals analysis (GCSA) [16].

This paper proposes a new approach for diagnosing gear tooth faults based on generator stator current analysis. Several research papers have discussed this technique [17,18,19,20,21,22,23], but they have overlooked one important aspect: the gearbox’s dynamic modeling. In this work, a mathematical model is developed for a system consisting of a two-stage gearbox and a doubly fed induction generator (DFIG) to simulate the system under healthy and faulty conditions. The fast Fourier transform (FFT) analysis method detects faults in the stator current signal.

This paper is organized as follows: Section 2 discusses the modeling of the studied system including the modeling of two stages gearbox, mesh stiffness calculation for (healthy, crack and a broken tooth) faults and modeling of the doubly fed induction generator (DFIG). The simulation results and discussion are in Section 3. Concluding remarks are presented in Section 4.

2. Modeling of the Studied System

Diagnosing gearbox failures is a way to improve wind turbine reliability, prevent catastrophic failures, and reduce downtime and maintenance costs [24]. Several methods can be used to detect faults in gearboxes such as lubrication analysis, acoustic emission analysis and vibration analysis. However, the challenge lies in the application of electrical analysis [15].

Our system is made up of a two-stage gearbox (multiplier), a four gear ($Z_1$ = 100 teeth, $Z_2$ = 29, $Z_3$ = 90 and $Z_4$ = 36) with a multiplication ratio of 8.62 ($Z_1 Z_3/Z_2 Z_4$ = 8.62). The input shaft of this gearbox is connected to a turbine, while the output shaft is connected to a doubly fed induction generator (DFIG).

The primary aim of this study is to estimate the stator currents from the output of the doubly fed induction generator (DFIG) to diagnose faults in the wind turbine gearbox. To achieve this objective, a detailed modeling of the system’s homelands (see Figure 2) must be provided. This modeling includes an in-depth analysis of the system’s components and their interactions, which allows us to simulate the system’s behavior under different operating conditions. By accurately modeling the system, we can identify any deviations from normal operation and diagnose any potential faults in the gearbox.

2.1. Modeling of Two Stages Gearbox

A nonlinear dynamic model of 10 degrees of freedom (DOF) was established using the concentrated mass method and ignoring the axial constraints, as shown in Figure 3. Some phenomena, such as friction between gear tooth, gearbox housing, and backlash, are ignored by this model. The following notation was used in this study [25,26]:

• ${\mathrm{m}}_1,{\mathrm{m}}_2,{\mathrm{m}}_3,{\mathrm{m}}_4:$ mass of the gears 1, 2, 3 and 4, respectively;
• ${\mathrm{I}}_1,{\mathrm{I}}_2,{\mathrm{I}}_3,{\mathrm{I}}_4:$ mass moment of inertia of gears 1, 2, 3 and 4, respectively;
• ${\mathrm{T}}_{\mathrm{T}},{\mathrm{T}}_{\mathrm{M}}$: turbine torque, machine torque;
• ${\mathrm{R}}_{\mathrm{b}1},{\mathrm{R}}_{\mathrm{b}2},{\mathrm{R}}_{\mathrm{b}3},{\mathrm{R}}_{\mathrm{b}4}$: base circle radius of gears 1, 2, 3 and 4, respectively;
• ${\mathsf{\theta }}_{\mathrm{T}}$,${\mathsf{\theta }}_{\mathrm{M}}$: turbine/machine angular displacement;
• ${\mathsf{\theta }}_1$,${\mathsf{\theta }}_2,{\mathsf{\theta }}_3,{\mathsf{\theta }}_4$: angular displacement of gears 1, 2, 3 and 4, respectively;
• ${\mathrm{I}}_{\mathrm{T}}$, ${\mathrm{I}}_{\mathrm{M}}$: mass moment of inertia of turbine/machine;
• ${\mathrm{C}}_{\mathrm{t}1},{\mathrm{C}}_{\mathrm{t}2}$: time-varying damping coefficient of first and second stage;
• ${\mathrm{C}}_{\mathrm{y}1},{\mathrm{C}}_{\mathrm{y}2},{\mathrm{C}}_{\mathrm{y}3},{\mathrm{C}}_{\mathrm{y}4}$: the radial damping coefficient along the direction of the y-axis of each bearing 1, 2, 3 and 4, respectively;
• ${\mathrm{C}}_1,{\mathrm{C}}_2,{\mathrm{C}}_3$: the damping coefficient of input, intermediate and output shafts, respectively;
• ${\mathrm{K}}_{\mathrm{t}1},{\mathrm{K}}_{\mathrm{t}2}$: time-varying meshing stiffness of first and second stages;
• ${\mathrm{k}}_1,{\mathrm{k}}_2,{\mathrm{k}}_3$: torsional stiffness of input, intermediate and output shafts, respectively.

We focused on the model response in the y direction, representing the effect of time-varying mesh stiffness, as the model response in the x direction is transient. According to Newton’s second law and considering the previous assumptions, the torsional vibration dynamics equations of the two-stage spur gear are as follows:

The equations of motion in the y direction of gear and pinion of the first stage are:

$${\mathrm{m}}_1{\ddot{\mathrm{y}}}_1=-{\mathrm{c}}_{{\mathrm{y}}_1}{\dot{\mathrm{y}}}_1-{\mathrm{k}}_{{\mathrm{y}}_1}{\mathrm{y}}_1+\mathrm{k}{\mathrm{t}}_1\left({\mathrm{R}}_{\mathrm{b}1}{\mathsf{\theta }}_1-{\mathrm{R}}_{\mathrm{b}2}{\mathsf{\theta }}_2-{\mathrm{y}}_1+{\mathrm{y}}_2\right)+{\mathrm{C}}_{\mathrm{t}1}\left({\mathrm{R}}_{\mathrm{b}1}{\dot{\mathsf{\theta }}}_1-{\mathrm{R}}_{\mathrm{b}2}{\dot{\mathsf{\theta }}}_2-{\dot{\mathrm{y}}}_1+{\dot{\mathrm{y}}}_2\right)$$

(1)

$${\mathrm{m}}_2{\ddot{\mathrm{y}}}_2=-{\mathrm{c}}_{{\mathrm{y}}_2}{\dot{\mathrm{y}}}_2-{\mathrm{k}}_{{\mathrm{y}}_2}{\mathrm{y}}_2-\mathrm{k}{\mathrm{t}}_1\left({\mathrm{R}}_{\mathrm{b}1}{\mathsf{\theta }}_1-{\mathrm{R}}_{\mathrm{b}2}{\mathsf{\theta }}_2-{\mathrm{y}}_1+{\mathrm{y}}_2\right)-{\mathrm{C}}_{\mathrm{t}1}\left({\mathrm{R}}_{\mathrm{b}1}{\dot{\mathsf{\theta }}}_1-{\mathrm{R}}_{\mathrm{b}2}{\dot{\mathsf{\theta }}}_2-{\dot{\mathrm{y}}}_1+{\dot{\mathrm{y}}}_2\right)$$

(2)

The equations of motion in the y direction of gear and pinion of the second stage are:

$${\mathrm{m}}_3{\ddot{\mathrm{y}}}_3=-{\mathrm{c}}_{{\mathrm{y}}_3}{\dot{\mathrm{y}}}_3-{\mathrm{k}}_{{\mathrm{y}}_3}{\mathrm{y}}_3+\mathrm{k}{\mathrm{t}}_2\left({\mathrm{R}}_{\mathrm{b}3}{\mathsf{\theta }}_3-{\mathrm{R}}_{\mathrm{b}4}{\mathsf{\theta }}_4-{\mathrm{y}}_3+{\mathrm{y}}_4\right)+{\mathrm{C}}_{\mathrm{t}2}\left({\mathrm{R}}_{\mathrm{b}3}{\dot{\mathsf{\theta }}}_3-{\mathrm{R}}_{\mathrm{b}4}{\dot{\mathsf{\theta }}}_4-{\dot{\mathrm{y}}}_3+{\dot{\mathrm{y}}}_4\right)$$

(3)

$${\mathrm{m}}_4{\ddot{\mathrm{y}}}_4=-{\mathrm{c}}_{{\mathrm{y}}_4}{\dot{\mathrm{y}}}_4-{\mathrm{k}}_{{\mathrm{y}}_4}{\mathrm{y}}_4-\mathrm{k}{\mathrm{t}}_2\left({\mathrm{R}}_{\mathrm{b}3}{\mathsf{\theta }}_3-{\mathrm{R}}_{\mathrm{b}4}{\mathsf{\theta }}_4-{\mathrm{y}}_3+{\mathrm{y}}_4\right)-{\mathrm{C}}_{\mathrm{t}2}\left({\mathrm{R}}_{\mathrm{b}3}{\dot{\mathsf{\theta }}}_3-{\mathrm{R}}_{\mathrm{b}4}{\dot{\mathsf{\theta }}}_4-{\dot{\mathrm{y}}}_3+{\dot{\mathrm{y}}}_4\right)$$

(4)

The rotary motion equations of the turbine and machine are:

$${\mathrm{I}}_{\mathrm{T}}{\ddot{\mathsf{\theta }}}_{\mathrm{T}}={\mathrm{T}}_{\mathrm{T}}-{\mathrm{k}}_1\left({\mathsf{\theta }}_{\mathrm{T}}-{\mathsf{\theta }}_1\right)-{\mathrm{C}}_1\left({\dot{\mathsf{\theta }}}_{\mathrm{T}}-{\dot{\mathsf{\theta }}}_1\right)$$

(5)

$${\mathrm{I}}_{\mathrm{M}}{\ddot{\mathsf{\theta }}}_{\mathrm{M}}=-{\mathrm{T}}_{\mathrm{M}}+{\mathrm{k}}_3\left({\mathsf{\theta }}_4-{\mathsf{\theta }}_{\mathrm{M}}\right)+{\mathrm{C}}_3\left({\dot{\mathsf{\theta }}}_4-{\dot{\mathsf{\theta }}}_{\mathrm{M}}\right)$$

(6)

The rotary motion equations of gear and pinion of the first stage are:

$${\mathrm{I}}_1{\ddot{\mathsf{\theta }}}_1={\mathrm{k}}_1\left({\mathsf{\theta }}_{\mathrm{T}}-{\mathsf{\theta }}_1\right)+{\mathrm{C}}_1\left({\dot{\mathsf{\theta }}}_{\mathrm{T}}-{\dot{\mathsf{\theta }}}_1\right)-{\mathrm{R}}_{\mathrm{b}1}\left[{\mathrm{k}}_{\mathrm{t}1}\left({\mathrm{R}}_{\mathrm{b}1}{\mathsf{\theta }}_1-{\mathrm{R}}_{\mathrm{b}2}{\mathsf{\theta }}_2-{\mathrm{y}}_1+{\mathrm{y}}_2\right)\right]+{\mathrm{C}}_{\mathrm{t}1}\left({\mathrm{R}}_{\mathrm{b}1}{\dot{\mathsf{\theta }}}_1-{\mathrm{R}}_{\mathrm{b}2}{\dot{\mathsf{\theta }}}_2-{\dot{\mathrm{y}}}_1+{\dot{\mathrm{y}}}_2\right)$$

(7)

$${\mathrm{I}}_2{\ddot{\mathsf{\theta }}}_2=-{\mathrm{k}}_2\left({\mathsf{\theta }}_2-{\mathsf{\theta }}_3\right)-{\mathrm{C}}_2\left({\dot{\mathsf{\theta }}}_2-{\dot{\mathsf{\theta }}}_3\right)+{\mathrm{R}}_{\mathrm{b}2}\left[{\mathrm{k}}_{\mathrm{t}1}\left({\mathrm{R}}_{\mathrm{b}1}{\mathsf{\theta }}_1-{\mathrm{R}}_{\mathrm{b}2}{\mathsf{\theta }}_2-{\mathrm{y}}_1+{\mathrm{y}}_2\right)\right]+{\mathrm{C}}_{\mathrm{t}1}\left({\mathrm{R}}_{\mathrm{b}1}{\dot{\mathsf{\theta }}}_1-{\mathrm{R}}_{\mathrm{b}2}{\dot{\mathsf{\theta }}}_2-{\dot{\mathrm{y}}}_1+{\dot{\mathrm{y}}}_2\right)$$

(8)

The rotary motion equations of gear and pinion of the second stage are:

$${\mathrm{I}}_3{\ddot{\mathsf{\theta }}}_3={\mathrm{k}}_2\left({\mathsf{\theta }}_2-{\mathsf{\theta }}_3\right)+{\mathrm{C}}_2\left({\dot{\mathsf{\theta }}}_2-{\dot{\mathsf{\theta }}}_3\right)-{\mathrm{R}}_{\mathrm{b}3}\left[{\mathrm{k}}_{\mathrm{t}2}\left({\mathrm{R}}_{\mathrm{b}3}{\mathsf{\theta }}_3-{\mathrm{R}}_{\mathrm{b}4}{\mathsf{\theta }}_4-{\mathrm{y}}_3+{\mathrm{y}}_4\right)\right]+{\mathrm{C}}_{\mathrm{t}2}\left({\mathrm{R}}_{\mathrm{b}3}{\dot{\mathsf{\theta }}}_3-{\mathrm{R}}_{\mathrm{b}4}{\dot{\mathsf{\theta }}}_4-{\dot{\mathrm{y}}}_3+{\dot{\mathrm{y}}}_4\right)$$

(9)

$${\mathrm{I}}_4{\ddot{\mathsf{\theta }}}_4=-{\mathrm{k}}_3\left({\mathsf{\theta }}_4-{\mathsf{\theta }}_{\mathrm{M}}\right)-{\mathrm{C}}_3\left({\dot{\mathsf{\theta }}}_4-{\dot{\mathsf{\theta }}}_{\mathrm{M}}\right)+{\mathrm{R}}_{\mathrm{b}4}\left[{\mathrm{k}}_{\mathrm{t}2}\left({\mathrm{R}}_{\mathrm{b}3}{\mathsf{\theta }}_3-{\mathrm{R}}_{\mathrm{b}4}{\mathsf{\theta }}_4-{\mathrm{y}}_3+{\mathrm{y}}_4\right)\right]+ {\mathrm{C}}_{\mathrm{t}2}\left({\mathrm{R}}_{\mathrm{b}3}{\dot{\mathsf{\theta }}}_3-{\mathrm{R}}_{\mathrm{b}4}{\dot{\mathsf{\theta }}}_4-{\dot{\mathrm{y}}}_3+{\dot{\mathrm{y}}}_4\right)$$

(10)

The mesh damping $C_{ti}$ is assumed proportional to the mesh stiffness in this study, and their ratio ${\mu }_{ti}$ is calculated using the formula below [27]:

$${\mathsf{\mu }}_{\mathrm{t}\mathrm{i}}=2{\mathsf{\epsilon }}_{\mathrm{t}\mathrm{i}}\sqrt{\frac{1}{{\mathrm{K}}_{\mathrm{m}}\left(\frac{{\mathrm{R}}_{\mathrm{b}1}^2}{{\mathrm{I}}_1}+\frac{{\mathrm{R}}_{\mathrm{b}2}^2}{{\mathrm{I}}_2}\right)}}$$

(11)

where: ${\mathrm{k}}_{\mathrm{m}}$ is the mean value of the mesh stiffness.

In this study, the damping ratio, denoted by ${\epsilon }_{ti}$, is taken to be 0.10 [27]. The mesh damping coefficient is calculated last, using the formula below:

$${\mathrm{C}}_{\mathrm{t}\mathrm{i}}={\mathsf{\mu }}_{\mathrm{t}\mathrm{i}}{\mathrm{k}}_{\mathrm{t}\mathrm{i}}$$

(12)

The values of torsional stiffnesses used in [25] are adopted in this study: the torsional stiffnesses of the input shaft ${\mathrm{k}}_1$ = 1381.6 Nm/ rad, of the intermediate shaft ${\mathrm{k}}_2$ = 3555.446 Nm/rad and of the output shaft ${\mathrm{k}}_3$ = 4824.705 Nm/ rad.

Mesh Stiffness Formulation

The gear mesh stiffness is a time-varying parameter that reflects the gear mesh conditions as the number of teeth in contact and the line of contact of the engaged gear teeth vary. It depends on the tooth geometry, position of the contact point, gear tooth deflections, gear tooth profile errors, gear hub torsional deformation, and local tooth faults [28].

(A)

Calculation of the mesh stiffness of a healthy gear

The potential energy technique presented by [29] was used as the basis for the gear stiffness model used in this study. This approach treats the gear tooth as a non-uniform cantilever beam that is subjected to a force. The various energies that can be stored in the gear tooth structure are calculated using beam theory. These energies are then converted to appropriate stiffnesses [30].

The analytical expressions of the Hertzian contact stiffness ${\mathrm{k}}_{\mathrm{h}}$, the bending stiffness ${\mathrm{k}}_{\mathrm{b}}$, shear stiffness ${\mathrm{k}}_{\mathrm{s}}$ and axial compressive stiffness ${\mathrm{k}}_{\mathrm{a}}$ are given as follows [31,32]:

$${\mathrm{k}}_{\mathrm{h}}=\frac{\mathsf{\pi }\mathrm{E}\mathrm{L}}{4\left(1-{\mathrm{v}}^2\right)}$$

(13)

$$\frac{1}{{\mathrm{k}}_{\mathrm{a}}}={{\displaystyle \int }}_{-{\mathsf{\alpha }}_1}^{{\mathsf{\alpha }}_2}\frac{\left({\mathsf{\alpha }}_2-\mathsf{\alpha }\right) \cos \mathsf{\alpha }\mathrm{s}\mathrm{i}{n}^2{\mathsf{\alpha }}_1}{2\mathrm{E}\mathrm{L}[\sin \mathsf{\alpha }+\left({\mathsf{\alpha }}_2-\mathsf{\alpha }\right)\cos \mathsf{\alpha }]}\mathrm{d}\mathsf{\alpha }$$

(14)

$$\frac{1}{{\mathrm{k}}_{\mathrm{s}}}={{\displaystyle \int }}_{-{\mathsf{\alpha }}_1}^{{\mathsf{\alpha }}_2}\frac{1.2\left(1+\mathrm{v}\right)\left({\mathsf{\alpha }}_2-\mathsf{\alpha }\right)\cos \mathsf{\alpha }{\cos }^2{\mathsf{\alpha }}_1}{\mathrm{E}\mathrm{L}[\sin \mathsf{\alpha }+\left({\mathsf{\alpha }}_2-\mathsf{\alpha }\right)\cos \mathsf{\alpha }]}\mathrm{d}\mathsf{\alpha }$$

(15)

$$\frac{1}{{\mathrm{k}}_{\mathrm{b}}}={{\displaystyle \int }}_{-{\mathsf{\alpha }}_1}^{{\mathsf{\alpha }}_2}\frac{3{\left\{1+\cos {\mathsf{\alpha }}_1\left.\left[\left({\mathsf{\alpha }}_2-\mathsf{\alpha }\right)\sin \mathsf{\alpha }-\cos \mathsf{\alpha }\right]\right\}\right.}^2\left({\mathsf{\alpha }}_2-\mathsf{\alpha }\right)\cos \mathsf{\alpha }}{2\mathrm{E}\mathrm{L}{[\sin \mathsf{\alpha }+\left({\mathsf{\alpha }}_2-\mathsf{\alpha }\right)\cos \mathsf{\alpha }]}^3}\mathrm{d}\mathsf{\alpha }$$

(16)

where $\mathrm{E}$ is the Young modulus, $\mathrm{v}$ is the Poisson’s ratio, and $\mathrm{L}$ is the tooth face width; the other variables are shown in Figure 4.

Hence, the total effective for the single tooth pair meshing duration can be represented as follows [29,33]:

$${\mathrm{k}}_{\mathrm{t}}=\frac{1}{\frac{1}{{\mathrm{k}}_{\mathrm{h}}}+\frac{1}{{\mathrm{k}}_{\mathrm{b}1}}+\frac{1}{{\mathrm{k}}_{\mathrm{s}1}}+\frac{1}{{\mathrm{k}}_{\mathrm{a}1}}+\frac{1}{{\mathrm{k}}_{\mathrm{b}2}}+\frac{1}{{\mathrm{k}}_{\mathrm{s}2}}+\frac{1}{{\mathrm{k}}_{\mathrm{a}2}}}$$

(17)

The driving and driven gears are denoted, respectively, by subscripts 1 and 2.

For the double-tooth-pair meshing duration, the total effective mesh stiffness is the sum of the two pair’s stiffness, which can be expressed as [29,33]:

$${\mathrm{k}}_{\mathrm{t}}={{\displaystyle \sum }}_{\mathrm{i}=1}^2\frac{1}{\frac{1}{{\mathrm{K}}_{\mathrm{h},\mathrm{i}}}+\frac{1}{{\mathrm{K}}_{\mathrm{b}1,\mathrm{i}}}+\frac{1}{{\mathrm{K}}_{\mathrm{s}1,\mathrm{i}}}+\frac{1}{{\mathrm{K}}_{\mathrm{a}1,\mathrm{i}}}+\frac{1}{{\mathrm{K}}_{\mathrm{b}2,\mathrm{i}}}+\frac{1}{{\mathrm{K}}_{\mathrm{s}2,\mathrm{i}}}+\frac{1}{{\mathrm{K}}_{\mathrm{a}2,\mathrm{i}}}}$$

(18)

where i represents the i th pair of the meshing tooth.

(B)

Mesh stiffness calculation of gear with a cracked tooth

When a crack has been initiated at the base of one of the pinion gears, the bending and shear stiffness will change. This phenomenon will occur because when the crack is present, the effective moment of inertia of the surface and the cross-sectional area will change. So, for the single tooth mesh period, the total effective mesh stiffness is given by Equation (19), while Figure 5 shows a schematic of a pinion gear with a cracked tooth [34]. In this figure, ${\mathrm{q}}_1$ represents the length of the crack, ${\mathrm{h}}_{\mathrm{c}1}$ is the distance between the root of the crack and the tooth centerline, corresponding to point 𝑮 on the tooth profile and the angle ${\alpha }_g$, M represents the torque causing the bending effect of Fa, and v is the crack intersection angle.

$${\mathrm{K}}_{\mathrm{t}\_\mathrm{c}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{k}}=\frac{1}{\frac{1}{{\mathrm{k}}_{\mathrm{h}}}+\frac{1}{{\mathrm{k}}_{\mathrm{b}\_\mathrm{c}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{k}}}+\frac{1}{{\mathrm{k}}_{\mathrm{s}\_\mathrm{c}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{k} }}+\frac{1}{{\mathrm{k}}_{\mathrm{a}1}}+\frac{1}{\mathrm{k}{ }_{\mathrm{b}2}}+\frac{1}{{\mathrm{k}}_{\mathrm{s}2}}+\frac{1}{{\mathrm{k}}_{\mathrm{a}2}}}$$

(19)

${\mathrm{k}}_{\mathrm{b}\_\mathrm{c}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{k}}$ and ${\mathrm{k}}_{\mathrm{s}\_\mathrm{c}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{k}}$ are the cracked tooth’s bending and shear mesh stiffness, respectively. They are calculated as follows [34]:

$$\frac{1}{{\mathrm{k}}_{\mathrm{b}\_\mathrm{c}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{k}}}={{\displaystyle \int }}_{-{\mathsf{\alpha }}_1}^{{\mathsf{\alpha }}_2}\frac{12{\left\{1+\cos {\mathsf{\alpha }}_1\left.\left[\left({\mathsf{\alpha }}_2-\mathsf{\alpha }\right)\sin \mathsf{\alpha }-\cos \mathsf{\alpha }\right]\right\}\right.}^2\left({\mathsf{\alpha }}_2-\mathsf{\alpha }\right)\cos \mathsf{\alpha }}{\mathrm{E}\mathrm{L}[\sin {\mathsf{\alpha }}_2-\frac{{\mathrm{q}}_1}{{\mathrm{R}}_{\mathrm{b}1}}\sin \mathsf{\upsilon }+\sin \mathsf{\alpha }+\left({\mathsf{\alpha }}_2-\mathsf{\alpha }\right)\cos \mathsf{\alpha }]{ }^3}\mathrm{d}\mathsf{\alpha }$$

(20)

$$\frac{1}{{\mathrm{K}}_{\mathrm{s}\_\mathrm{c}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{k}}}={{\displaystyle \int }}_{-{\mathsf{\alpha }}_1}^{{\mathsf{\alpha }}_2}\frac{2.4\left(1+\mathrm{v}\right)\left({\mathsf{\alpha }}_2-\mathsf{\alpha }\right)\cos \mathsf{\alpha }{\left(\cos {\mathsf{\alpha }}_1\right)}^2}{\mathrm{E}\mathrm{L}[\sin {\mathsf{\alpha }}_2-\frac{{\mathrm{q}}_1}{{\mathrm{R}}_{\mathrm{b}1}}\sin \mathsf{\upsilon }+\sin \mathsf{\alpha }+\left({\mathsf{\alpha }}_2-\mathsf{\alpha }\right)\cos \mathsf{\alpha }] }\mathrm{d}\mathsf{\alpha }$$

(21)

(C)

Mesh stiffness calculation of gear with a broken tooth

Since the tooth will lose contact where the missing tooth was, only one pair of teeth will mesh during the first period of the double pair of teeth, contrary to the ideal scenario. As a result, the overall effective stiffness of the mesh will only be equal to the stiffness of a single pair of teeth. It will become:

$${\mathrm{K}}_{\mathrm{t}\_\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{k}}=\frac{1}{\frac{1}{{\mathrm{K}}_{\mathrm{h}}}+\frac{1}{{\mathrm{K}}_{\mathrm{b}1,1}}+\frac{1}{{\mathrm{K}}_{\mathrm{a}1,1}}+\frac{1}{{\mathrm{K}}_{\mathrm{a}2,1}}+\frac{1}{{\mathrm{K}}_{\mathrm{b}2,1}}}$$

(22)

2.2. Doubly Fed Induction Generator (DFIG) Modeling

The use of DFIG in wind energy systems offers several advantages. With the ability to control active and reactive power independently, wind turbines equipped with these generators can optimize their power output and improve efficiency. Additionally, the generators can operate at variable speeds, making them easier to integrate with the grid and providing reactive power support to maintain grid stability. Compared to other types of generators, asynchronous double-fed generators are cost-effective and have a higher power output, making them a popular and effective option for wind turbine manufacturers.

A DFIG is constructed similarly to a wound rotor induction generator and has a three-phase stator winding and a three-phase rotor winding. Slip rings are used to feed the latter. If neglecting the effects of the slotting, we assumed that the permeability of the iron parts is infinite and the flux density is radial in the air gap. The simplified and idealized DFIG model can be shown in Figure 6 [35,36].

The stator and rotor voltage are represented as follows [35,36]:

$$\left\{\begin{array}{l}{\mathrm{V}}_{\mathrm{sa}}={\mathrm{R}}_{\mathrm{s}}{\mathrm{I}}_{\mathrm{sa}}+\frac{\mathrm{d}}{\mathrm{dt}}{\mathsf{\Phi }}_{\mathrm{sa}} \\ {\mathrm{V}}_{\mathrm{sb}}={\mathrm{R}}_{\mathrm{s}}{\mathrm{I}}_{\mathrm{sb}}+\frac{\mathrm{d}}{\mathrm{dt}}{\mathsf{\Phi }}_{\mathrm{sb}}\\ {\mathrm{V}}_{\mathrm{sc}}={\mathrm{R}}_{\mathrm{s}}{\mathrm{I}}_{\mathrm{sc}}+\frac{\mathrm{d}}{\mathrm{dt}}{\mathsf{\Phi }}_{\mathrm{sc}}\end{array}\right.$$

(23)

$$\left\{\begin{array}{l}{\mathrm{V}}_{\mathrm{ra}}={\mathrm{R}}_{\mathrm{r}}{\mathrm{I}}_{\mathrm{ra}}+\frac{\mathrm{d}}{\mathrm{dt}}{\mathsf{\Phi }}_{\mathrm{ra}}\\ {\mathrm{V}}_{\mathrm{rb}}={\mathrm{R}}_{\mathrm{r}}{\mathrm{I}}_{\mathrm{rb}}+\frac{\mathrm{d}}{\mathrm{dt}}{\mathsf{\Phi }}_{\mathrm{rb}}\\ {\mathrm{V}}_{\mathrm{rc}}={\mathrm{R}}_{\mathrm{r}}{\mathrm{I}}_{\mathrm{rc}}+\frac{\mathrm{d}}{\mathrm{dt}}{\mathsf{\Phi }}_{\mathrm{rc}}\end{array}\right.$$

(24)

where the voltages applied to the stator are ${\mathrm{V}}_{\mathrm{s}\mathrm{a}}$, ${\mathrm{V}}_{\mathrm{s}\mathrm{b}}$ and ${\mathrm{V}}_{\mathrm{s}\mathrm{c}}$. The stator currents in phases A, B and C are ${\mathrm{I}}_{\mathrm{s}\mathrm{a}}$, ${\mathrm{I}}_{\mathrm{s}\mathrm{b}}$ and ${\mathrm{I}}_{\mathrm{s}\mathrm{c}}$; the rotor currents in phases a, b and c are ${\mathrm{I}}_{\mathrm{r}\mathrm{a}}$, ${\mathrm{I}}_{\mathrm{r}\mathrm{b}}$ and ${\mathrm{I}}_{\mathrm{r}\mathrm{c}}$; the rotor voltages are ${\mathrm{V}}_{\mathrm{r}\mathrm{a}}$, ${\mathrm{V}}_{\mathrm{r}\mathrm{b}}$ and ${\mathrm{V}}_{\mathrm{r}\mathrm{c}}$; and the stator and rotor fluxes are (${\mathrm{Փ}}_{\mathrm{s}\mathrm{a}}$,${\mathrm{Փ}}_{\mathrm{s}\mathrm{b}}\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{Փ}}_{\mathrm{s}\mathrm{c}}$), (${\mathrm{Փ}}_{\mathrm{r}\mathrm{a}}$,${\mathrm{Փ}}_{\mathrm{r}\mathrm{b}}$ and ${\mathrm{Փ}}_{\mathrm{r}\mathrm{c}}$).

The stator and rotor resistances of the DFIG are ${\mathrm{R}}_{\mathrm{s}}$ and ${\mathrm{R}}_{\mathrm{r}}$, respectively.

Applying the park transformation to the three phases of DFIG gives the general model of DFIG as follows [37,38,39,40]:

$$\left\{\begin{array}{l}{\mathrm{V}}_{\mathrm{s}\mathrm{d}}={\mathrm{R}}_{\mathrm{s}}{\mathrm{I}}_{\mathrm{s}\mathrm{d}}+\frac{\mathrm{d}}{\mathrm{d}\mathrm{t}}{\mathsf{\Phi }}_{\mathrm{s}\mathrm{d}}-{\mathsf{\omega }}_{\mathrm{s}}{\mathsf{\Phi }}_{\mathrm{s}\mathrm{q}}\\ {\mathrm{V}}_{\mathrm{s}\mathrm{q}}={\mathrm{R}}_{\mathrm{s}}{\mathrm{I}}_{\mathrm{s}\mathrm{q}}+\frac{\mathrm{d}}{\mathrm{d}\mathrm{t}}{\mathsf{\Phi }}_{\mathrm{s}\mathrm{q}}+{\mathsf{\omega }}_{\mathrm{s}}{\mathsf{\Phi }}_{\mathrm{s}\mathrm{d}}\\ {\mathrm{V}}_{\mathrm{r}\mathrm{d}}={\mathrm{R}}_{\mathrm{r}}{\mathrm{I}}_{\mathrm{r}\mathrm{d}}+\frac{\mathrm{d}}{\mathrm{d}\mathrm{t}}{\mathsf{\Phi }}_{\mathrm{r}\mathrm{d}}-{\mathsf{\omega }}_{\mathrm{r}}{\mathsf{\Phi }}_{\mathrm{r}\mathrm{q}}\\ {\mathrm{V}}_{\mathrm{r}\mathrm{q}}={\mathrm{R}}_{\mathrm{r}}{\mathrm{I}}_{\mathrm{r}\mathrm{q}}+\frac{\mathrm{d}}{\mathrm{d}\mathrm{t}}{\mathsf{\Phi }}_{\mathrm{r}\mathrm{q}}+{\mathsf{\omega }}_{\mathrm{r}}{\mathsf{\Phi }}_{\mathrm{r}\mathrm{d}}\end{array}\right.$$

(25)

where ${\mathrm{V}}_{\mathrm{s}\mathrm{d}},{\mathrm{V}}_{\mathrm{s}\mathrm{q}},{\mathrm{V}}_{\mathrm{r}\mathrm{d}},{\mathrm{V}}_{\mathrm{r}\mathrm{q}},{\mathrm{I}}_{\mathrm{s}\mathrm{d}},{\mathrm{I}}_{\mathrm{s}\mathrm{q}},{\mathrm{I}}_{\mathrm{r}\mathrm{d}},{\mathrm{I}}_{\mathrm{r}\mathrm{q}},{\mathrm{Փ}}_{\mathrm{s}\mathrm{d}},{\mathrm{Փ}}_{\mathrm{s}\mathrm{q}},{\mathrm{Փ}}_{\mathrm{r}\mathrm{d}} \mathrm{a}\mathrm{n}\mathrm{d} {\mathrm{Փ}}_{\mathrm{r}\mathrm{q}}$ represent, respectively, the components along the d and q axes of the stator and rotor voltages, currents and flux.

${\mathsf{\omega }}_{\mathrm{s}}$ and ${\mathsf{\omega }}_{\mathrm{r}}$ are the stator and rotor angular frequencies in $\mathrm{r}\mathrm{a}\mathrm{d}/\mathrm{s}$ and the relationship between them is:

$${\mathsf{\omega }}_{\mathrm{r}}+{\mathsf{\omega }}_{\mathrm{m}}={\mathsf{\omega }}_{\mathrm{s}}$$

(26)

where ${\omega }_m$ is the electrical angular frequency of the machine. In this case, the rotational speed of the gearbox output gear is defined by the following equation:

$${\mathsf{\omega }}_{\mathrm{m}}=\frac{\mathrm{d}{\mathsf{\theta }}_4}{\mathrm{d}\mathrm{t}}$$

(27)

Similarly, ${\mathrm{\Omega }}_m$ is the mechanical angular speed, linked to the electrical frequency by a pair of poles, p:

$${\mathsf{\omega }}_{\mathrm{m}}=p.{\Omega }_{\mathrm{m}}$$

(28)

The currents and fluxes for DFIG are related by the following equations:

$$\left\{\begin{array}{c}{\mathsf{\Phi }}_{\mathrm{sd}}={\mathrm{L}}_{\mathrm{s}}{\mathrm{I}}_{\mathrm{sd}}+\mathrm{M}{\mathrm{I}}_{\mathrm{rd}}\\ {\mathsf{\Phi }}_{\mathrm{sq}}={\mathrm{L}}_{\mathrm{s}}{\mathrm{I}}_{\mathrm{sq}}+\mathrm{M}{\mathrm{I}}_{\mathrm{rq}}\\ {\mathsf{\Phi }}_{\mathrm{rd}}={\mathrm{L}}_{\mathrm{r}}{\mathrm{I}}_{\mathrm{rd}}+\mathrm{M}{\mathrm{I}}_{\mathrm{sd}}\\ {\mathsf{\Phi }}_{\mathrm{rq}}={\mathrm{L}}_{\mathrm{r}}{\mathrm{I}}_{\mathrm{rq}}+\mathrm{M}{\mathrm{I}}_{\mathrm{sq}}\end{array}\right.$$

(29)

${\mathrm{L}}_{\mathrm{s}}$ and ${\mathrm{L}}_{\mathrm{r}}$ are the stator and rotor leakage inductances, respectively, with M serving as the mutual inductance.

We find Equation (30) by putting Equation (29) in (25) [41]:

$$\left\{\begin{array}{l}{\mathrm{V}}_{\mathrm{sd}}={\mathrm{R}}_{\mathrm{s}}{\mathrm{I}}_{\mathrm{sd}}+{\mathrm{L}}_{\mathrm{s}}\frac{{\mathrm{dI}}_{\mathrm{sd}}}{\mathrm{dt}}+\mathrm{M}\frac{{\mathrm{dI}}_{\mathrm{rd}}}{\mathrm{dt}}-{\mathsf{\omega }}_{\mathrm{s}}{\mathrm{L}}_{\mathrm{s}}{\mathrm{I}}_{\mathrm{sq}}-{\mathsf{\omega }}_{\mathrm{s}}{\mathrm{MI}}_{\mathrm{rq}} \\ {\mathrm{V}}_{\mathrm{sq}}={\mathrm{R}}_{\mathrm{s}}{\mathrm{I}}_{\mathrm{sq}}+{\mathrm{L}}_{\mathrm{s}}\frac{{\mathrm{dI}}_{\mathrm{sq}} }{\mathrm{dt}}+\mathrm{M}\frac{{\mathrm{dI}}_{\mathrm{rq}}}{\mathrm{dt}}+{\mathsf{\omega }}_{\mathrm{s}}{\mathrm{L}}_{\mathrm{s}}{\mathrm{I}}_{\mathrm{sd}}+{\mathsf{\omega }}_{\mathrm{s}}{\mathrm{MI}}_{\mathrm{rd}}\\ {\mathrm{V}}_{\mathrm{rd}}={\mathrm{R}}_{\mathrm{r}}{\mathrm{I}}_{\mathrm{rd}}+{\mathrm{L}}_{\mathrm{r}}\frac{{\mathrm{dI}}_{\mathrm{rd}}}{\mathrm{dt}}+\mathrm{M}\frac{{\mathrm{dI}}_{\mathrm{sd}}}{\mathrm{dt}}-{\mathsf{\omega }}_r{\mathrm{L}}_{\mathrm{r}}{\mathrm{I}}_{\mathrm{rq}}-{\mathsf{\omega }}_r{\mathrm{MI}}_{\mathrm{sq}}\\ {\mathrm{V}}_{\mathrm{rq}}={\mathrm{R}}_{\mathrm{r}}{\mathrm{I}}_{\mathrm{rq}}+{\mathrm{L}}_{\mathrm{r}}\frac{{\mathrm{dI}}_{\mathrm{rq}}}{\mathrm{dt}}+\mathrm{M}\frac{{\mathrm{dI}}_{\mathrm{sq}}}{\mathrm{dt}}+{\mathsf{\omega }}_r{\mathrm{L}}_{\mathrm{r}}{\mathrm{I}}_{\mathrm{rq}}+{\mathsf{\omega }}_r{\mathrm{MI}}_{\mathrm{sd}}\end{array}\right.$$

(30)

Equation (30) is written in the matrix form as follows:

$$\left[\begin{array}{c}{\mathrm{V}}_{\mathrm{sd}}\\ {\mathrm{V}}_{\mathrm{sq}}\\ {\mathrm{V}}_{\mathrm{rd}}\\ {\mathrm{V}}_{\mathrm{rq}}\end{array}\right]=\left[\begin{array}{cccc}{\mathrm{R}}_{\mathrm{s}}& -{\mathsf{\omega }}_{\mathrm{s}}{\mathrm{L}}_{\mathrm{s}}& 0& -{\mathsf{\omega }}_{\mathrm{s}}\mathrm{M}\\ {\mathsf{\omega }}_{\mathrm{s}}{\mathrm{L}}_{\mathrm{s}}& {\mathrm{R}}_{\mathrm{s}}& {\mathsf{\omega }}_{\mathrm{s}}\mathrm{M}& 0\\ 0& \left(-{\mathsf{\omega }}_{\mathrm{s}}+{\mathsf{\omega }}_{\mathrm{m}}\right)\mathrm{M}& {\mathrm{R}}_{\mathrm{r}}& \left(-{\mathsf{\omega }}_{\mathrm{s}}+{\mathsf{\omega }}_{\mathrm{m}}\right){\mathrm{L}}_{\mathrm{r}}\\ \left({\mathsf{\omega }}_{\mathrm{s}}-{\mathsf{\omega }}_{\mathrm{m}}\right)\mathrm{M}& 0& \left({\mathsf{\omega }}_{\mathrm{s}}-{\mathsf{\omega }}_{\mathrm{m}}\right){\mathrm{L}}_{\mathrm{r}}& {\mathrm{R}}_{\mathrm{r}}\end{array}\right]\left[\begin{array}{c}{\mathrm{I}}_{\mathrm{sd}}\\ {\mathrm{I}}_{\mathrm{sq}}\\ {\mathrm{I}}_{\mathrm{rd}}\\ {\mathrm{I}}_{\mathrm{rq}}\end{array}\right]+\left[\begin{array}{cccc}{\mathrm{L}}_{\mathrm{s}}& 0& \mathrm{M}& 0\\ 0& {\mathrm{L}}_{\mathrm{s}}& 0& \mathrm{M}\\ \mathrm{M}& 0& {\mathrm{L}}_{\mathrm{r}}& 0\\ 0& \mathrm{M}& 0& {\mathrm{L}}_{\mathrm{r}}\end{array}\right]\frac{\mathrm{d}}{\mathrm{d}\mathrm{t}}\left[\begin{array}{c}{\mathrm{I}}_{\mathrm{sd}}\\ {\mathrm{I}}_{\mathrm{sq}}\\ {\mathrm{I}}_{\mathrm{rd}}\\ {\mathrm{I}}_{\mathrm{rq}}\end{array}\right]$$

(31)

Starting from (31), we get the following equation:

$$\begin{gathered} \frac{\mathrm{d}}{\mathrm{dt}}\left[\begin{array}{c}{\mathrm{I}}_{\mathrm{sd}}\\ {\mathrm{I}}_{\mathrm{sq}}\\ {\mathrm{I}}_{\mathrm{rd}}\\ {\mathrm{I}}_{\mathrm{rq}}\end{array}\right]={\left[\begin{array}{cccc}{\mathrm{L}}_{\mathrm{s}}& 0& \mathrm{M}& 0\\ 0& {\mathrm{L}}_{\mathrm{s}}& 0& \mathrm{M}\\ \mathrm{M}& 0& {\mathrm{L}}_{\mathrm{r}}& 0\\ 0& \mathrm{M}& 0& {\mathrm{L}}_{\mathrm{r}}\end{array}\right]}^{-1}\left[\begin{array}{cccc}-{\mathrm{R}}_{\mathrm{s}}& {\mathsf{\omega }}_{\mathrm{s}}{\mathrm{L}}_{\mathrm{s}}& 0& {\mathsf{\omega }}_{\mathrm{s}}\mathrm{M}\\ -{\mathsf{\omega }}_{\mathrm{s}}{\mathrm{L}}_{\mathrm{s}}& -{\mathrm{R}}_{\mathrm{s}}& -{\mathsf{\omega }}_{\mathrm{s}}\mathrm{M}& 0\\ 0& \left({\mathsf{\omega }}_{\mathrm{s}}-{\mathsf{\omega }}_{\mathrm{m}}\right)\mathrm{M}& -{\mathrm{R}}_{\mathrm{r}}& \left({\mathsf{\omega }}_{\mathrm{s}}-{\mathsf{\omega }}_{\mathrm{m}}\right){\mathrm{L}}_{\mathrm{r}}\\ -\left({\mathsf{\omega }}_{\mathrm{s}}-{\mathsf{\omega }}_{\mathrm{m}}\right)\mathrm{M}& 0& -\left({\mathsf{\omega }}_{\mathrm{s}}-{\mathsf{\omega }}_{\mathrm{m}}\right){\mathrm{L}}_{\mathrm{r}}& -{\mathrm{R}}_{\mathrm{r}}\end{array}\right]\left[\begin{array}{c}{\mathrm{I}}_{\mathrm{sd}}\\ {\mathrm{I}}_{\mathrm{sq}}\\ {\mathrm{I}}_{\mathrm{rd}}\\ {\mathrm{I}}_{\mathrm{rq}}\end{array}\right] \\ +{\left[\begin{array}{cccc}{\mathrm{L}}_{\mathrm{s}}& 0& \mathrm{M}& 0\\ 0& {\mathrm{L}}_{\mathrm{s}}& 0& \mathrm{M}\\ \mathrm{M}& 0& {\mathrm{L}}_{\mathrm{r}}& 0\\ 0& \mathrm{M}& 0& {\mathrm{L}}_{\mathrm{r}}\end{array}\right]}^{-1}\left[\begin{array}{c}{\mathrm{V}}_{\mathrm{sd}}\\ {\mathrm{V}}_{\mathrm{sq}}\\ {\mathrm{V}}_{\mathrm{rd}}\\ {\mathrm{V}}_{\mathrm{rq}}\end{array}\right] \end{gathered}$$

(32)

After simplifying the previous equation, we get the following equation:

$$\frac{\mathrm{d}}{\mathrm{dt}}\left[\begin{array}{c}{\mathrm{I}}_{\mathrm{sd}}\\ {\mathrm{I}}_{\mathrm{sq}}\\ {\mathrm{I}}_{\mathrm{rd}}\\ {\mathrm{I}}_{\mathrm{rq}}\end{array}\right]=\left[\begin{array}{cccc}-{\mathsf{\lambda }}_2& {\mathsf{\lambda }}_1{\mathsf{\omega }}_{\mathrm{m}}+{\mathsf{\omega }}_{\mathrm{s}}& {\mathsf{\lambda }}_4& {\mathsf{\lambda }}_6{\mathsf{\omega }}_{\mathrm{m}}\\ -{\mathsf{\lambda }}_1{\mathsf{\omega }}_{\mathrm{m}}-{\mathsf{\omega }}_{\mathrm{s}}& -{\mathsf{\lambda }}_2& -{\mathsf{\lambda }}_6{\mathsf{\omega }}_{\mathrm{m}}& {\mathsf{\lambda }}_4\\ {\mathsf{\lambda }}_5& -{\mathsf{\lambda }}_7{\mathsf{\omega }}_{\mathrm{m}}& -{\mathsf{\lambda }}_3& -\frac{{\mathsf{\omega }}_{\mathrm{m}}}{\mathsf{\sigma }}+{\mathsf{\omega }}_{\mathrm{s}}\\ {\mathsf{\lambda }}_7{\mathsf{\omega }}_{\mathrm{m}}& {\mathsf{\lambda }}_5& \frac{{\mathsf{\omega }}_{\mathrm{m}}}{\mathsf{\sigma }}-{\mathsf{\omega }}_{\mathrm{s}}& -{\mathsf{\lambda }}_3\end{array}\right]\left[\begin{array}{c}{\mathrm{I}}_{\mathrm{sd}}\\ {\mathrm{I}}_{\mathrm{sq}}\\ {\mathrm{I}}_{\mathrm{rd}}\\ {\mathrm{I}}_{\mathrm{rq}}\end{array}\right]+\left[\begin{array}{cccc}{\mathsf{\gamma }}_1& 0& -{\mathsf{\gamma }}_3& 0\\ 0& {\mathsf{\gamma }}_1& 0& -{\mathsf{\gamma }}_3\\ -{\mathsf{\gamma }}_3& 0& {\mathsf{\gamma }}_2& 0\\ 0& -{\mathsf{\gamma }}_3& 0& {\mathsf{\gamma }}_2\end{array}\right]\left[\begin{array}{c}{\mathrm{V}}_{\mathrm{sd}}\\ {\mathrm{V}}_{\mathrm{sq}}\\ {\mathrm{V}}_{\mathrm{rd}}\\ {\mathrm{V}}_{\mathrm{rq}}\end{array}\right]$$

(33)

where:

$${\mathsf{\lambda } }_1=\frac{1-\mathsf{\sigma }}{\mathsf{\sigma }} ;{\mathsf{\lambda } }_2=\frac{{\mathrm{R}}_{\mathrm{s}}}{{\mathsf{\sigma }\mathrm{L}}_{\mathrm{s}}};{\mathsf{\lambda } }_3=\frac{{\mathrm{R}}_{\mathrm{r}}}{{\mathsf{\sigma }\mathrm{L}}_{\mathrm{r}}}{ \mathsf{\lambda } }_4=\frac{{\mathrm{R}}_{\mathrm{r}}\mathrm{M}}{{\mathsf{\sigma }\mathrm{L}}_{\mathrm{s}}{\mathrm{L}}_{\mathrm{r}}} ;{\mathsf{\lambda } }_5=\frac{{\mathrm{R}}_{\mathrm{S}}\mathrm{M}}{{\mathsf{\sigma }\mathrm{L}}_{\mathrm{s}}{\mathrm{L}}_{\mathrm{r}}};{ \mathsf{\lambda } }_6=\frac{\mathrm{M}}{{\mathsf{\sigma }\mathrm{L}}_{\mathrm{s}}};{ \mathsf{\lambda } }_7=\frac{\mathrm{M}}{{\mathsf{\sigma }\mathrm{L}}_{\mathrm{r}}}$$

$${\mathsf{\gamma }}_1=\frac{1}{{\mathsf{\sigma }\mathrm{L}}_{\mathrm{s}}} ;{\mathsf{\gamma }}_2=\frac{1}{{\mathsf{\sigma }\mathrm{L}}_{\mathrm{r}}} ;{\mathsf{\gamma }}_3=\frac{\mathrm{M}}{{\mathsf{\sigma }\mathrm{L}}_{\mathrm{s}}{\mathrm{L}}_{\mathrm{r}}}$$

$\mathsf{\sigma }$ is the leakage factor equal to $1-\frac{{\mathrm{M}}^2}{{\mathrm{L}}_{\mathrm{s}}{\mathrm{L}}_{\mathrm{r}}}$.

Finally, the state model describes the doubly fed induction machine given as follows:

$$\left\{\begin{array}{c}\dot{\mathrm{X}}=\mathrm{A}\mathrm{X}+\mathrm{B}\mathrm{U}\\ \mathrm{Y}=\mathrm{C}\mathrm{X}\end{array}\right.$$

(34)

where:

$$\mathrm{X}={\left[{\mathrm{I}}_{\mathrm{sd}}{ \mathrm{I}}_{\mathrm{sq}}{ \mathrm{I}}_{\mathrm{rd}}{ \mathrm{I}}_{\mathrm{rq}}\right]}^{\mathrm{T}}$$

(35)

$$\mathrm{U}={\left[{\mathrm{V}}_{\mathrm{sd}}{ \mathrm{V}}_{\mathrm{sq}}{ \mathrm{V}}_{\mathrm{rd}}{ \mathrm{V}}_{\mathrm{rq}}\right]}^{\mathrm{T}}$$

(36)

$$\mathrm{A}=\left[\begin{array}{cccc}-{\mathsf{\lambda }}_2& {\mathsf{\lambda }}_1{\mathsf{\omega }}_{\mathrm{m}}+{\mathsf{\omega }}_{\mathrm{s}}& {\mathsf{\lambda }}_4& {\mathsf{\lambda }}_6{\mathsf{\omega }}_{\mathrm{m}}\\ -{\mathsf{\lambda }}_1{\mathsf{\omega }}_{\mathrm{m}}-{\mathsf{\omega }}_{\mathrm{s}}& -{\mathsf{\lambda }}_2& -{\mathsf{\lambda }}_6{\mathsf{\omega }}_{\mathrm{m}}& {\mathsf{\lambda }}_4\\ {\mathsf{\lambda }}_5& -{\mathsf{\lambda }}_7{ \mathsf{\omega }}_{\mathrm{m}}& -{\mathsf{\lambda }}_3& -\frac{{\mathsf{\omega }}_{\mathrm{m}}}{\mathsf{\sigma }}+{\mathsf{\omega }}_{\mathrm{s}}\\ {\mathsf{\lambda }}_7{\mathsf{\omega }}_{\mathrm{m}}& {\mathsf{\lambda }}_5& \frac{{\mathsf{\omega }}_{\mathrm{m}}}{\mathsf{\sigma }}-{\mathsf{\omega }}_{\mathrm{s}}& -{\mathsf{\lambda }}_3\end{array}\right]$$

(37)

$$\mathrm{B}=\left[\begin{array}{cccc}{\mathsf{\gamma }}_1& 0& -{\mathsf{\gamma }}_3& 0\\ 0& {\mathsf{\gamma }}_1& 0& -{\mathsf{\gamma }}_3\\ -{\mathsf{\gamma }}_3& 0& {\mathsf{\gamma }}_2& 0\\ 0& -{\mathsf{\gamma }}_3& 0& {\mathsf{\gamma }}_2\end{array}\right]$$

(38)

$$\mathrm{C}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$

(39)

The electromagnetic torque generated by the DFIG is represented by the following equation [42]:

$${\mathrm{T}}_{\mathrm{em}}=\mathrm{PM}\left({\mathrm{I}}_{\mathrm{rd}}{\mathrm{I}}_{\mathrm{sq}}-{\mathrm{I}}_{\mathrm{rq}}{\mathrm{I}}_{\mathrm{sd}}\right)$$

(40)

3. Results and Discussion

To obtain digital solutions, we used Matlab/Simulink. Figure 7 shows the studied system’s Simulink model, and

Table 1. Two-stage gear drive system parameters (delivered by the manufacturer).

Symbol	Quantity	Value	Symbol	Quantity	Value
E	Young modulus	Pa	${\mathrm{m}}_4$	Second pinion mass	0.294 kg
$\mathrm{v}$	Poisson’s ratio	0.3	${\mathrm{I}}_1$	First gear inertia	$4.82\times {10}^{-3} \mathrm{k}\mathrm{g}·{\mathrm{m}}^2$
m	Modulus	1.5 mm	${\mathrm{I}}_2$	First pinion inertia	$4.76\times {10}^{-5} \mathrm{k}\mathrm{g}·{\mathrm{m}}^2$
${\mathrm{k}}_{\mathrm{b}}$	Bearing stiffness	$6.56\times {10}^8 \mathrm{N}/\mathrm{m}$	${\mathrm{I}}_3$	Second gear inertia	$3.89\times {10}^{-3} \mathrm{k}\mathrm{g}·{\mathrm{m}}^2$
${\mathrm{C}}_{\mathrm{b}}$	Bearing damping	$1.8\times {10}^3 \mathrm{N}\mathrm{s}/\mathrm{m}$	${\mathrm{I}}_4$	Second pinion inertia	$1.21\times {10}^{-4} \mathrm{k}\mathrm{g}·{\mathrm{m}}^2$
${\mathrm{m}}_1$	First gear mass	1.74 kg	${\mathrm{I}}_{\mathrm{T}}$	Turbine inertia	$3.3\times {10}^{-4} \mathrm{k}\mathrm{g}·{\mathrm{m}}^2$
${\mathrm{m}}_2$	First pinion mass	0.16 kg	${\mathrm{I}}_{\mathrm{M}}$	Machine inertia	$6.63\times {10}^{-3} \mathrm{k}\mathrm{g}·{\mathrm{m}}^2$
${\mathrm{m}}_3$	Second gear mass	1.79 kg	F	Turbine friction coefficient	$0.0016$

and

Table 2. DFIG parameters [43].

Symbol	Quantity	Value	Symbol	Quantity	Value
${\mathrm{P}}_{\mathrm{n}}$	Rated Power	$1.5 \mathrm{K}\mathrm{w}$	$\mathrm{R}\mathrm{s}$	Stator Resistance	1.18 Ω
${\mathrm{V}}_{\mathrm{s}}$	Rated Voltage	220/380 V	$\mathrm{R}\mathrm{r}$	Rotor Resistance	1.66 Ω
$f_{\mathrm{s}}$	Frequency	50 Hz	$\mathrm{L}\mathrm{s}$	Stator Inductance	0.20 H
${\mathrm{\Omega }}_{\mathrm{s}}$	Rated Speed	3000 rpm	$\mathrm{L}\mathrm{r}$	Rotor Inductance	$0.18 \mathrm{H}$
p	Number of Pole Pairs	1	$\mathrm{M}$	Mutual Inductance	0.17 H

list the system’s parameters.

The two-stage gearbox was implemented using Equations (1)–(10). This gearbox is driven by a constant turbine torque $T_T$ with a value of 61.21 Nm.

On the other hand, the doubly-fed induction generator (DFIG) was implemented using Equations (25)–(40). The DFIG trained at a fixed speed of 310 radians/s, powered directly by two perfect three-phase voltag sources. One is located at the stator with a frequency of $f_s=$ 50 Hz and an amplitude of $220\sqrt{2}$ V, and the other is located at the rotor with an amplitude of 10 V and a frequency equal to the rotor frequency ($f_r$ = s·$f_s)$.

A broken tooth fault is created at the level of the first gear in the first stage. The following Figure 8 and Figure 9 show the changes in the variation of meshing stiffness in the first and second stages, respectively, during the rotation of the gears. From these figures, we can see that the meshing stiffness of the gear system is approximately parabolic.

To explain and extract the characteristics of the first stage stiffness signal shown in Figure 8, a spectral analysis (FFT) has been presented in Figure 10 and Figure 11 for healthy and faulty gear.

The spectrum of meshing stiffness of the first stage (${Kt}_1$) is represented in Figure 10; this spectrum shows a multitude of peaks whose amplitude decreases very quickly. The largest amplitude line is observed at the fundamental meshing frequency $f_m$ = 572.6 Hz ($f_m$=$f_1$·$Z_1$ = 572.6 Hz), then the next in decreasing order of amplitude at 1145 Hz, then at 1718 Hz, then at 2290 Hz, etc. For the healthy gear, the signal spectrum makes it possible to read the meshing frequency and its harmonics (2$·f_m$ … n.$f_m$).

Figure 11 represents the spectrum of meshing stiffness for the first stage in the presence of broken tooth faults of a Gear one notices on the spectrum of the meshing stiffness of the lateral bands on either side of the meshing frequency ($f_m$ = $f_1 Z_1$ = 572.6 Hz) and its harmonics, the distances between the peaks is equal to 5.7 Hz. This frequency corresponds to the frequency of rotation of the Gear it is the frequency of repetition of the broken tooth.

Figure 12 represents the accelerations of $\ddot{y_1}$, $\ddot{y_2}$, $\ddot{y_3}$ and $\ddot{y_4}$, respectively. We notice multitudes of peaks separated by a distance of 0.1747 s between two adjacent peaks. It represents the period of rotation of the first shaft (${\mathrm{T}}_1=\frac{1}{f_1}$=$\frac{1}{5.726}$). It can be seen that with a broken tooth fault, the whole system will be influenced by this defect. Although the fault is on the first floor, it can be observed in the vibration responses of the entire system.

In Figure 13, the rotational speed of the gearbox output gear is presented. The speed remains steady at a value of 310 radians/s, which is consistent with the normal operating range of the system. However, the appearance of sharp peaks at specific times indicates the presence of faults in the system.

These peaks appear every 0.1747 s, corresponding to the period of repetition of the broken tooth.

Researchers using vibration analysis to detect gearbox problems will find these results a valuable resource. They offer interesting data and references that advance diagnostic and fault prevention techniques.

To clarify the effect of crack fault on the mesh stiffness, we presented the gear mesh stiffness for different crack lengths with the same crack angle v = 22.5° in Figure 14; we note a degradation of the value of total stiffness proportionally with the crack size due to the change of the bending stiffness.

In this party, we consider a crack fault at the gable root with a depth of $q_{1 }$ = 7.5 mm (50% tooth width) and an angle of v = 22.5° on the gear at stage 1. The accelerations $\ddot{y_1}$, $\ddot{y_2}$, $\ddot{y_3}$ and $\ddot{y_4}$ are depicted in Figure 15. Once again, we observe numerous peaks that are smaller compared to the case of a broken tooth. These peaks are separated by a distance of 0.1747 s, representing the period of rotation of the first shaft (period of repetition of the cracked tooth).

Figure 16 shows the angular rotation speed of the output shaft of the gearbox. A constant speed of 310 radians/s is observed. However, periodic peaks appear every 0.1747 s. These peaks are related to the rotation of the first wheel with a cracked tooth. In comparison to the peaks observed in the case of a broken tooth, these spikes are less pronounced.

Figure 17 represents the stator current of the DFIG. We can see that it has a sinusoidal format after a transient state.

To detect broken and cracked tooth faults in the gear system, a spectral analysis of phases currents A, B, and C are presented in Figure 18, Figure 19 and Figure 20, respectively. From these Figures, we can notice that the current spectra in healthy gear are composed of a fundamental $f_s=50 \mathrm{H}\mathrm{z}$ (network frequency) and for the case, broken, and cracked teeth are composed of a fundamental $f_s$ and a series of harmonics with frequencies ($f_s+K{f}_1$) and ($f_s-K{f}_1$), where $f_1$ is the rotation frequency of the first shaft corresponding to the defective wheel and K is the integer number.

A small difference, of 0.016 Hz, between the calculated theoretical frequency ($f_1$ = 5.726 Hz) and that of the spectrum ($f_1$ = 5.71 Hz) is quite normal.

As a comparison, we can notice that the spectrum of stator current in the case of a cracked tooth has the same characteristic as the spectrum of the case of a broken tooth with smaller amplitudes.

4. Conclusions

A method based on the analysis of stator current signals has been presented for fault detection of wind gearbox equipped with a DFIG system. This one was modeled to simulate the behavior in a healthy, cracked, and broken tooth fault of the first gear of the gearbox. The results showed that the proposed stator current signal analysis approach effectively detects and diagnoses various gear failures. The presence of a crack or broken fault in the gear tooth of the gearbox induces a modulation of the DFIG stator current signals by the frequencies of the shaft corresponding to the defect gear. These results can provide a preliminary basis for detecting and diagnosing this type of failure.

Looking forward, the authors plan to further test the proposed method in a practical setting using a dedicated test bed. By doing so, they hope to validate the effectiveness of the approach in a real-world scenario and optimize the method for industrial application. This study serves as a promising starting point for future research aimed at improving wind turbine gearbox reliability and reducing the overall cost of wind energy production.