Impact Statement
Classical methods for damage detection in gearboxes have proven to be successful in scenarios with steady operating conditions. Unfortunately, when considering wind turbine gearboxes, standard data analysis techniques are strongly challenged by the presence of wind fluctuations. To cope with the stochastic nature of the arising sensor signals due to wind turbulence, we introduce a damage detection approach for gearboxes based on a multi-resolution analysis of the sensor data.
1. Introduction
The condition monitoring of wind turbine gearboxes presents a challenging scenario that is often not amenable to classical data analysis techniques due to the existing operating conditions. Note that wind turbines are one of the main sources of renewable energy (Redl et al., Reference Redl, Hein, Buck, Graichen and Jones2021) nowadays and gearbox-related defects are among those problems with the highest amount of downtime hours per failure (Dao et al., Reference Dao, Kazemtabrizi and Crabtree2019). Therefore, damage detection in gearboxes is a significant aspect for keeping wind energy production high. To cope with the difficulties in the sensor signals, arising due to the prevailing wind fluctuations, we propose an approach for damage detection in gearboxes based on the analysis of sensor data with the multi-resolution dynamic mode decomposition (mrDMD).
Sensor data analysis is widely used to determine the health condition of mechanical devices. Analyzing the data generated by sensors installed on these devices, it is possible to continually monitor their integrity and detect anomalous behaviors, damages, and faults when they arise. Optimizing the monitoring process is of great interest because detecting damages in their early stage empowers companies to prevent more severe problems that could lead to a significant loss in terms of energy production and money.
Damage detection in gearbox vibration signals is a relatively easy task in steady operating conditions, as for example in a laboratory, see Doebling et al. (Reference Doebling, Farrar, Prime and Shevitz1996) for an overview of established methods. Unfortunately, when we deal with wind turbines, steady operating conditions represent an unrealistic scenario that does not taken into account weather condition, temperature variation, and, most importantly, wind turbulence, all of which results in varying load conditions on the blades of the wind turbine. These factors strongly affect gearbox vibration signals, introducing a stochastic component that makes them nondeterministic and far from smooth. In this scenario, damage detection is a much more difficult task since the data analysis strategies have not only to deal with signal anomalies induced by damage, but also by wind turbulence.
In this work, we consider the gearbox model and the scenario developed in Kahraman and Singh (Reference Kahraman and Singh1991) and further refined in Antoniadou et al. (Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015). The numerically produced signals include a stochastic component simulating the varying load condition caused by wind turbulence. Several kinds of damage can appear in gearboxes, here we focus on the crack of a gear’s tooth. The crack of a gear’s tooth is a bending fatigue failure. Bending fatigue failures usually occur in three progressive stages: crack initiation, crack propagation, and tooth fracture or breakage (Errichello, Reference Errichello2002). According to Ma et al. (Reference Ma, Song, Pang and Wen2014), “tooth fracture is the most serious fault in gearbox and it may cause complete failure of the gear.” Thus, the crack of a gear’s tooth is a damage in its early stage because the continuous interaction of the cracked tooth with the other gear leads to a progressive propagation of the crack until the breakage of the tooth occurs, leading to more severe damage that may even cause the total failure of the gear.
Historically speaking, fault detection in gearboxes was first performed using methods that focus on analyzing the spectral properties of sensor signals (Kumar et al., Reference Kumar, Ismail, Zhao, Noori, Yadav, Chen, Singh, Altabey, Silik, Kumar, Kumar and Balodi2021), such as the fast Fourier transform (FFT) (Rao et al., Reference Rao, Kim and Hwang2011) and cepstrum analysis (Badaoui et al., Reference Badaoui, Guillet and Danière2004). However, in the scenario we consider, these methods are strongly challenged by weather conditions, such as wind turbulence and temperature variation, which induce a stochastic component to the signals studied, changing their spectra nondeterministically over time. This stochastic change of the signals’ spectral properties strongly affects the effectiveness and reliability of spectral methods (Kumar et al., Reference Kumar, Ismail, Zhao, Noori, Yadav, Chen, Singh, Altabey, Silik, Kumar, Kumar and Balodi2021). Empirical mode decomposition (EMD) (Huang et al., Reference Huang, Shen, Long, Wu, Shih, Zheng, Yen, Tung and Liu1998) was proposed as an alternative to overcome the issues arising when studying spectral characteristics of signals. In the case of a damaged gearbox, EMD separates the time-domain signals in several modes that may contain damage features and may be used to detect faults. Unfortunately, EMD is affected by the so-called “mode mixing problem.” The problem consists of the fact that the modes produced by EMD either contain information of the signal relative to widely disparate time-scales or similar time-scale information resides in different modes (Wu and Huang, Reference Wu and Huang2009; Xu et al., Reference Xu, Yang, Wang, Bin, Chen and Zhao2016). The main implications are that the physical meaning of each mode is unclear and that a proper interpretation of the information they contain is not a trivial task. This does not mean that such modes cannot be studied or selected according to some specific strategy in order to extract meaningful insights about the data, it rather means that there is no principled way to do that.
In this work, we propose mrDMD (Kutz et al., Reference Kutz, Fu and Brunton2016b), a variation of dynamic mode decomposition (DMD) (Schmid, Reference Schmid2010), as an alternative to EMD that incorporates all its advantages and overcomes the mode mixing problem as well. mrDMD is a data-driven algorithm that has been used to perform different tasks, such as: spatio-temporal filtering of video or multi-scale separation of complex weather systems (Kutz et al., Reference Kutz, Fu and Brunton2016b). We propose it as a novel approach to detect the crack of a tooth in a wind turbine gearbox by analyzing gearbox vibration signals generated under varying load conditions. mrDMD is able to decompose the signal in spatio-temporal modes that capture its geometrical characteristics in the time-domain at different time-scales. The difference to EMD is that it associates to each of its modes a complex number that could either represent a frequency or an energy content, according to the point of view one wants to consider, and thanks to this additional information high/low-frequency or energy content separation of the modes is possible. On the one hand, low-frequency modes represent those structures that vary slowly over time, which therefore are not affected by the nondeterministic variation caused by wind turbulence that mainly influences the signal at smaller time-scales. On the other hand, high-frequency modes represent those structures that evolve at higher frequencies and are strongly affected by wind turbulence and sudden changes in the gear stiffness caused by the cracked tooth we want to identify. We will see that computing the instantaneous amplitude of high-frequency modes, via the Hilbert transform (HT), will enable us to easily identify the anomalies arising in vibration signals from the crack of a gear’s tooth.
We aim to show the potential of mrDMD as a fault detection method, in a context in which the analyzed vibration signals are affected by stochastic perturbations. We study the proposed approach on simulated signals representing the vibration response of a simple gearbox model with two spur gears and experimental data under controlled situations from a benchmark study.
The next sections are organized as follows: Section 2 presents an overview of related work on the topic of early damage detection in wind turbine gearboxes under varying load condition. Section 3 introduces our mrDMD-based strategy for tooth damage identification. In Section 4, the gearbox model is introduced and numerical simulations of acceleration signals representing the vibration response of the modeled gearbox are considered. Finally, Section 5 demonstrates the performance of our method, and shows the issues experienced by FFT as a spectral method, an approach based on time synchronous averaging (TSA), and an EMD-based approach due to the presence of wind fluctuations.
2. Related Work
Generally speaking, condition monitoring methods consist of damage detection methods, with the difference that one knows a priori about the frequency bands related to damage for the different monitored components. Reviews of methods used for damage detection in wind turbines can be found in Doebling et al. (Reference Doebling, Farrar, Prime and Shevitz1996), Kumar et al. (Reference Kumar, Ismail, Zhao, Noori, Yadav, Chen, Singh, Altabey, Silik, Kumar, Kumar and Balodi2021), Salameh et al. (Reference Salameh, Cauet, Etien, Sakout and Rambault2018), and Sharma (Reference Sharma2021). The most classical techniques used for condition monitoring, and more in general for damage detection, mainly focus on analyzing the spectral characteristics of vibration signals, such as the spectral kurtosis (Antoni, Reference Antoni2006), Fourier analysis and modulation sidebands (Inalpolat and Kahraman, Reference Inalpolat and Kahraman2010), and cepstrum analysis (Badaoui et al., Reference Badaoui, Guillet and Danière2004). These methods have proven their capability to provide good results in identifying damages under simplified conditions, such as steady loading of the gearboxes. However, when considering time-varying load conditions, the spectral characteristics of gearbox vibration signals change nondeterministically with time. This stochastic behavior cannot be handled using Fourier analysis as the Fourier transform expands a signal as a linear combination of wave functions with constant frequency over time. To overcome this issue, many time–frequency methods have been employed in damage detection, such as: the Wigner–Ville distribution (Staszewski et al., Reference Staszewski, Worden and Tomlinson1997), the Vold–Kalman filter (Feng et al., Reference Feng, Zhu and Zhang2019), wavelet analysis (Staszewski and Tomlinson, Reference Staszewski and Tomlinson1994), and cyclostationary analysis (Antoni and Randall, Reference Antoni and Randall2002; Xin et al., Reference Xin, Hamzaoui and Antoni2020). Probably, the most popular and successful of these methods is wavelet analysis (Peng and Chu, Reference Peng and Chu2004; Hu et al., Reference Hu, Tu, Li and Meng2018). Nonetheless, in the context we are considering, such a method has the disadvantage of relying on a fixed set of basis functions. This fact has a detrimental effect on its effectiveness and reliability in processing signals affected by perturbations of stochastic nature and that have spectral properties that change randomly with time. An additional approach to perform fault detection is given by entropy analysis (Sharma and Parey, Reference Sharma and Parey2016). However, many entropy features are available for different types of faults, and it is not clear what are the right features for a given fault. Moreover, under varying load condition, the stochastic behavior of the signal affects the performance of entropies which may show insignificant results, hence making them nonreliable (Sharma and Parey, Reference Sharma and Parey2016; Sharma, Reference Sharma2021). Another major contribution to the gear diagnostics field is TSA (Randall, Reference Randall2021). TSA is used to extract signal components associated with certain frequencies, as the meshing frequency, that may highlight damage features. Unfortunately, also this technique rely on being able to analyze signals recorded under steady speed and load conditions and it is strongly affected by the fluctuating operating conditions caused by wind turbulence.
Other strategies for damage detection, which try to overcome the problems experienced with TSA, entropy analysis and time–frequency methods, are based on the EMD (Huang et al., Reference Huang, Shen, Long, Wu, Shih, Zheng, Yen, Tung and Liu1998). The main advantage of using EMD with respect to entropy analysis and time–frequency methods is that it decomposes the signal into several components, or modes, which may contain meaningful information related to damage. After decomposing the vibration signal, the instantaneous frequency and amplitude of each component can be estimated, most commonly by applying the HT, in order to extract features of the signal that allow a better determination of whether there is damage or not. The decomposition’s procedure is based on local geometric characteristics of the signal in the time-domain, and it acts at different time-scales. It is completely data driven and there is no need of an a priori choice of a set of functions or a mother wavelet. Moreover, EMD-based methods for machinery damage identification do not rely on any kind of a priori chosen window function or any assumption about the regularity of the signal for any time-span (Junsheng et al., Reference Junsheng, Dejie and Yu2007; Feng et al., Reference Feng, Liang, Zhang and Hou2012; Antoniadou et al., Reference Antoniadou, Worden, Manson, Dervilis, Taylor and Farrar2013). As a consequence of that, these methods lead, in some cases, to a better estimation of even small variations in the signals. In general, EMD-based methods have the promise of high-quality results at a low computational cost. In Antoniadou et al. (Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015), an EMD-based strategy is developed to identify the crack of a gear’s tooth in a scenario that, similar to this work, takes into account the effects of the wind turbulence.
Unfortunately, EMD-based strategies have some limitations too. The major problem is the mode mixing problem studied in Wu and Huang (Reference Wu and Huang2009) and further analyzed in the wind turbine damage detection context in Antoniadou et al. (Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015). The mode mixing problem consists in the fact that each of the modes obtained after the signal’s decomposition, either contains information of the signal relative to widely disparate time-scales, or similar time-scale information resides in different modes (Wu and Huang, Reference Wu and Huang2009; Xu et al., Reference Xu, Yang, Wang, Bin, Chen and Zhao2016). In recent years, work has been carried out to provide a primary theoretical framework for the development of EMD. Despite the steps forward that give us a theoretical understanding of the algorithm when it is employed in a simplified scenario, for example, simple signals with only a pure oscillation component (Ge et al., Reference Ge, Chen, Yu, Chen and An2018), there is still no full mathematical justification for the general application of the EMD procedure, as it was already pointed out in Antoniadou et al. (Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015). This further affects the interpretability and reliability of the method.
In this work, we propose mrDMD as an alternative to EMD. An mrDMD-based procedure for damage detection in wind turbines was already developed in Dang et al. (Reference Dang and Lv2018). However, in Dang et al. (Reference Dang and Lv2018), the authors focused on the rolling bearing fault rather than the crack of a tooth, and they did not consider turbulent wind conditions as we do here. Moreover, there is a fundamental difference between the method we develop and the method they proposed: they analyze the spectral characteristics of the signals’ structures associated with the slow modes computed via mrDMD, while here, we focus on studying the information contained in high-frequency modes in their time-domain, avoiding all the issues related to examining the spectra of signals generated under turbulent wind conditions. This difference between the two approaches can be related to the different nature of the damages analyzed. On the one hand, we have the crack of a tooth that primarily affects gear mesh frequency and high order harmonics of the signals (del Rincon et al., Reference del Rincon, Viadero, Iglesias, de Juan, Garcia and Sancibrian2012; Antoniadou et al., Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015). On the other hand, the rolling bearing fault affects different frequency bands, including the low-frequencies (Hu et al., Reference Hu, Xiang, Xu and Lin2019).
3. Multi-resolution Dynamic Mode Decomposition for Damage Detection
The DMD algorithm was developed in the fluid dynamic community for the purpose of feature extraction (Schmid, Reference Schmid2010). Since its creation numerous variants of the original algorithm have been developed (Erichson et al., Reference Erichson, Brunton and Kutz2016; Arbabi and Mezić, Reference Arbabi and Mezić2017; Clainche and Vega, Reference Clainche and Vega2017; Hemati et al., Reference Hemati, Rowley, Deem and Cattafesta2017; Matsumoto and Indinger, Reference Matsumoto and Indinger2017). Moreover, DMD has been applied in several contexts, such as finance (Mann and Kutz, Reference Mann and Kutz2016), epidemiology (Proctor and Eckhoff, Reference Proctor and Eckhoff2015), highway traffic forecasting (Avila and Mezić, Reference Avila and Mezić2020), and biomedical informatics (Ingabire et al., Reference Ingabire, Wu, Amos, He, Peng, Wang, Li, Chen, Feng, Rao and Ren2021). It is important to mention that DMD has a solid theoretical interpretation through Koopman operator theory (Korda and Mezić, Reference Korda and Mezić2017), which makes the method understandable and interpretable. We now introduce the basics of DMD.
3.1. Dynamic mode decomposition
Let
$ \mathbf{y}\left({t}_i\right)\in {\mathrm{\mathbb{R}}}^n $
be a column vector containing data from sensors, evaluated at
$ n $
different points, at time
$ {t}_i\in {\mathrm{\mathbb{R}}}^{+} $
. Assuming we have vectors representing
$ T+1 $
consecutive equispaced time steps, we arrange them in two matrices:
$$ \mathbf{Y}\hskip0.35em =\hskip0.35em \left[\mathbf{y}\left({t}_0\right),\mathbf{y}\left({t}_1\right),\dots, \mathbf{y}\left({t}_{T-1}\right)\right]\hskip1em \mathrm{and}\hskip1em \overline{\mathbf{Y}}\hskip0.35em =\hskip0.35em \left[\mathbf{y}\left({t}_1\right),\mathbf{y}\left({t}_2\right),\dots, \mathbf{y}\left({t}_T\right)\right]. $$
These matrices are the so-called snapshots matrices and each one of their columns represents sensor values at a certain time step. The columns are arranged chronologically and those of
$ \overline{\mathbf{Y}} $
are shifted one time step forward in the future with respect to those of
$ \mathbf{Y} $
. Notice that, in this work, we always assume that the snapshots we analyze are taken at equispaced intervals in time, that is,
$ \Delta t\hskip0.35em =\hskip0.35em {t}_{i+1}-{t}_i $
for each
$ i\hskip0.35em =\hskip0.35em 0,1,\dots, T-1 $
with
$ \Delta t\in {\mathrm{\mathbb{R}}}^{+} $
.
DMD relies on the assumption that there exists a matrix
$ \mathbf{A}\in {\mathrm{\mathbb{R}}}^{n,n} $
such that
$$ \mathbf{y}\left({t}_{i+1}\right)\hskip0.35em =\hskip0.35em \mathbf{Ay}\left({t}_i\right)\hskip1em \mathrm{for}\;i\hskip0.35em =\hskip0.35em 0,1,\dots, T-1, $$
which written in a matrix form yields
$$ \overline{\mathbf{Y}}\hskip0.35em =\hskip0.35em \mathbf{AY}. $$
The algorithm’s goal is to find the matrix
$ A $
that best maps
$ \mathbf{y}\left({t}_i\right) $
into
$ \mathbf{y}\left({t}_{i+1}\right) $
, for
$ i\hskip0.35em =\hskip0.35em 0,1,\dots, \hskip0.35em T-1 $
, by solving the following minimization problem in the Frobenius norm
$$ \underset{\mathbf{A}}{\min}\parallel \overline{\mathbf{Y}}-\mathbf{AY}{\parallel}_F. $$
Thus, the matrix
$ A $
, solution of equation (4), is the matrix that best advances the sensor values in time, therefore best maps
$ \mathbf{Y} $
into
$ \overline{\mathbf{Y}} $
. The solution of the minimization problem (equation 4) is given by the matrix
$$ \mathbf{A}\hskip0.35em =\hskip0.35em \overline{\mathbf{Y}}{\mathbf{Y}}^{\dagger }, $$
where
$ {\mathbf{Y}}^{\dagger } $
is the Moore–Penrose inverse of the matrix
$ \mathbf{Y} $
. However, for computational reasons connected to the data’s dimension, computing the matrix
$ \mathbf{A} $
using equation (5) becomes computationally difficult and unstable (Tu et al., Reference Tu, Rowley, Luchtenburg, Brunton and Kutz2014). Instead, a low-rank approximation
$ \tilde{\mathbf{A}} $
of
$ \mathbf{A} $
is computed. The eigenvalues and eigenvectors of
$ \tilde{\mathbf{A}} $
are used to obtain an approximation of the eigenvalues and eigenvectors of
$ \mathbf{A} $
, called DMD eigenvalues and DMD modes, respectively.
Algorithm 1. Dynamic mode decomposition (DMD).
1: Arrange the data into matrices
$ \overline{\mathbf{Y}} $
and
$ \mathbf{Y} $
as in equation (1).
2: Compute the reduced singular value decomposition (SVD) of
$ \mathbf{Y} $
, that is,
$ \mathbf{Y}\approx {\mathbf{U}}_r{\boldsymbol{\Sigma}}_r{{\mathbf{V}}^{\ast}}_r $
. Here,
$ {\mathbf{U}}_r $
,
$ {\boldsymbol{\Sigma}}_r $
, and
$ {{\mathbf{V}}^{\ast}}_r $
are the rank-
$ r $
truncation of the matrices
$ \mathbf{U},\hskip0.35em \boldsymbol{\Sigma} $
and
$ {\mathbf{V}}^{\ast} $
computed via SVD and such that
$ \mathbf{Y}\hskip0.35em =\hskip0.35em \mathbf{U}\boldsymbol{\Sigma } {\mathbf{V}}^{\ast} $
. Further,
$ {\mathbf{V}}^{\ast} $
is the conjugate transpose of
$ \mathbf{V} $
and
$ r $
is a parameter to be chosen.
3: Construct the matrix
$ \tilde{\mathbf{A}}:= {{\mathbf{U}}^{\ast}}_r\overline{\mathbf{Y}}{\mathbf{V}}_r{\boldsymbol{\Sigma}}_r^{-1} $
.
4: Compute eigenvalues and eigenvectors of
$ \tilde{\mathbf{A}} $
, solving
$ \tilde{\mathbf{A}}\mathbf{W}\hskip0.35em =\hskip0.35em \boldsymbol{\Lambda} \mathbf{W} $
. With
$ \boldsymbol{\Lambda} $
a diagonal matrix of DMD eigenvalues and
$ \mathbf{W} $
a matrix where the columns are eigenvectors of
$ \tilde{\mathbf{A}} $
.
5: The DMD modes corresponding to DMD eigenvalues
$ \boldsymbol{\Lambda} $
are given by the columns of the matrix
$ \boldsymbol{\Phi} \hskip0.35em =\hskip0.35em \overline{\mathbf{Y}}{\mathbf{V}}_r{\mathbf{\sum}}_r^{-1}\mathbf{W} $
.
6: The reconstruction of the dynamics is given by the linear evolution
$$ \hat{\mathbf{y}}\left({t}_i\right)\hskip0.35em =\hskip0.35em \sum \limits_{k\hskip0.35em =\hskip0.35em 1}^M{b}_k{\boldsymbol{\phi}}_k\exp \left({\omega}_k{t}_i\right),\hskip1em i\hskip0.35em =\hskip0.35em 0,\dots, \hskip0.35em T, $$
where
$ M $
is the number of computed DMD modes, which depends on the size of
$ \tilde{\mathbf{A}} $
,
$ {\phi}_k $
is the
$ k $
-th DMD mode,
$ {b}_k $
is the
$ k $
-th entry of the vector
$ \mathbf{b}\hskip0.35em =\hskip0.35em {\boldsymbol{\Phi}}^{\dagger}\mathbf{y}\left({t}_0\right) $
, and
$ \hat{\mathbf{y}}\left({t}_i\right) $
is the reconstruction of the snapshot
$ \mathbf{y}\left({t}_i\right) $
. Moreover,
$ {\omega}_k $
is defined as the frequency value associated with the
$ k $
-th DMD mode and obtained via the DMD eigenvalue
$ {\lambda}_k $
:
$$ {\omega}_k\hskip0.35em =\hskip0.35em \log \left({\lambda}_k\right)/\Delta t. $$
Here,
$ \Delta t $
represents the time difference between two consecutive measurements. Since in this work we only consider equispaced measurements we have that
$ \Delta t\hskip0.35em =\hskip0.35em {t}_{i+1}-{t}_i $
for
$ i\hskip0.35em =\hskip0.35em 0,1,\dots, T-1 $
.
Algorithm 1 shows the main numerical steps of the DMD. The first important observation is that the temporal evolution of the analyzed signal can be reconstructed as a weighted sum of the DMD modes, as it can be seen in equation (6). Thus, DMD modes are signal’s components that encode meaningful spatio-temporal information. Notice that, since the DMD eigenvalue
$ {\lambda}_k $
is a complex number, in equation (6) the weight
$ \exp \left({\omega}_k{t}_i\right) $
, associated to the
$ k $
-th DMD mode, evolves periodically with the index “
$ i $
.” Consequently, each DMD mode contribution to the dynamic’s reconstruction changes periodically over time with oscillation frequency and growth/decay rate determined by its associated DMD eigenvalue. Furthermore, we can associate DMD modes with a concept of speed by looking at the frequencies associated with them. Given two DMD modes
$ {\boldsymbol{\phi}}_i,\ {\boldsymbol{\phi}}_j,\ i\ne j $
, we say that
$ {\boldsymbol{\phi}}_i $
is faster than
$ {\boldsymbol{\phi}}_j $
(or equivalently
$ {\boldsymbol{\phi}}_j $
is slower than
$ {\boldsymbol{\phi}}_i $
) if
$ \mid \exp \left({\omega}_i\right)\mid >\mid \exp \left({\omega}_j\right)\mid $
. On the one hand, slow modes represent those signal’s components whose contribution to the signal’s reconstruction vary slowly over time (Kutz et al., Reference Kutz, Fu and Brunton2016b). On the other hand, fast modes represent those structures whose contribution to the signal’s reconstruction vary with higher frequency.
Notice that, DMD operates at a single time-resolution, that is, DMD modes and eigenvalues fully characterize the signal’s evolution over all the given sampling window, as shown in equation (6). However, we aim to study signals that are locally affected by stochastic perturbations, and that can show a very different behavior at different time-scales. Thus, we aim to use an improved DMD that enables us to sift out information at different time-scales.
3.2. Multi-resolution dynamic mode decomposition
mrDMD (Kutz et al., Reference Kutz, Fu and Brunton2016b) is a variation of DMD, and it produces modes that extract spatial features at different time-scales. mrDMD basically consists of iteratively applying the DMD algorithm at different time-ranges. To be more precise, it starts by analyzing the largest sampling window, which consists of all the available snapshots. After that it computes DMD modes, identifies the slow ones and extracts them. Subsequently, it reduces the duration of the observation window by half, determines again slower modes in each half and extracts those modes from the residual signal. mrDMD iterates this process until the desired resolution or level of decomposition is reached. The final residual signal will be composed of fast modes representing spatial structures that characterize the signal locally, with a time-resolution determined by the number of iterations performed in the mrDMD algorithm. We now shortly recapture the formulation of mrDMD from Kutz et al. (Reference Kutz, Fu and Brunton2016b), where further details can be found.
In its first iteration, mrDMD reconstructs the time-domain signal as follows:
$$ \hat{\mathbf{y}}\left({t}_i\right)\hskip0.35em =\hskip0.35em \sum \limits_{k\hskip0.35em =\hskip0.35em 1}^M{b}_k{\boldsymbol{\phi}}_k\exp \left({\omega}_k{t}_i\right)\hskip0.35em =\hskip0.35em \underset{\mathrm{slow}\ \mathrm{modes}}{\underbrace{\sum \limits_{k\hskip0.35em =\hskip0.35em 1}^{m_1}{b}_k{\boldsymbol{\phi}}_k^{(1)}\exp \left({\omega}_k{t}_i\right)}}+\underset{\mathrm{fast}\ \mathrm{modes}}{\underbrace{\sum \limits_{k\hskip0.35em =\hskip0.35em {m}_1+1}^M{b}_k{\boldsymbol{\phi}}_k^{(1)}\exp \left({\omega}_k{t}_i\right)}}, $$
where
$ {\phi}_k^{(1)} $
represent the modes computed using the entire set of
$ T+1 $
snapshots,
$ M $
is the total number of modes, and
$ {m}_1 $
the number of slow modes at the first iteration level. The first sum in the right-hand side of equation (8) represents the slow-modes dynamics, whereas the second sum is everything else. At the second iteration level, DMD is performed after a split of the second sum
$$ {\mathbf{Y}}_{\left(T+1\right)/2}\hskip0.35em =\hskip0.35em {\mathbf{Y}}_{\left(T+1\right)/2}^{(1)}+{\mathbf{Y}}_{\left(T+1\right)/2}^{(2)}, $$
where the columns of the matrices
$ {\mathbf{Y}}_{\left(T+1\right)/2}^{(1)},{\mathbf{Y}}_{\left(T+1\right)/2}^{(2)} $
represent the fast modes reconstruction of the first
$ \left(T+1\right)/2 $
and last
$ \left(T+1\right)/2 $
snapshots, respectively. The iteration process continues by recursively removing slow-frequency components computed separately on each half of the snapshots. One builds the new matrices
$ {\mathbf{Y}}_{\left(T+1\right)/2},{\mathbf{Y}}_{\left(T+1\right)/4},{\mathbf{Y}}_{\left(T+1\right)/8},\dots $
until a suitable multi-resolution decomposition has been achieved.
The representation (equation 8) can be made more precise. Specifically, one must account for the number of levels
$ (L) $
of the decomposition, the number of time bins
$ (J) $
for each level, and the number of modes retained at each level
$ \left({m}_L\right) $
. Thus, the solution is parametrized by the following three indices:
$$ l\hskip0.35em =\hskip0.35em 1,\hskip0.35em 2,\dots, \hskip0.35em L:\mathrm{number}\ \mathrm{of}\ \mathrm{decomposition}\ \mathrm{levels}, $$
$$ j\hskip0.35em =\hskip0.35em 1,\hskip0.35em 2,\hskip0.35em \dots, \hskip0.35em J:\mathrm{number}\ \mathrm{of}\ \mathrm{time}\ \mathrm{bins}\;\mathrm{per}\;\mathrm{level}\;\left(J\hskip0.35em =\hskip0.35em {2}^{\left(l-1\right)}\right), $$
$$ k\hskip0.35em =1,\hskip0.35em 2,\dots, \hskip0.35em {m}_L:\mathrm{number}\ \mathrm{of}\ \mathrm{modes}\ \mathrm{extracted}\;\mathrm{at}\;\mathrm{level}\;L. $$
To formally determine the reconstructed snapshots
$ \hat{y}\left({t}_i\right) $
, the following indicator function is defined as
$$ {f}^{l,j}(t)\hskip0.35em := \hskip0.35em \left\{\begin{array}{rr}1,& \mathrm{for}\;t\in \left[{t}_j^l,{t}_{j+1}^l\right]\\ {}0,& \mathrm{elsewhere}\end{array}\right\},\mathrm{with}\hskip0.35em j\hskip0.35em =\hskip0.35em 1,\hskip0.35em 2,\dots, \hskip0.35em J\hskip0.35em \mathrm{and}\hskip0.35em J\hskip0.35em =\hskip0.35em {2}^{\left(l-1\right)}, $$
where
$ {t}_j^l,\ {t}_{j+1}^l\in {\mathrm{\mathbb{R}}}^{+} $
determine the interval of the
$ j $
-th time-bin at the
$ l $
-th decomposition level. The above indicator function is only nonzero in the interval, or time-bin, associated with the value of
$ j $
. The three indices and the indicator function in equation (11) are used to formally write equation (8) as
$$ \hat{\mathbf{y}}\left({t}_i\right)\hskip0.35em =\hskip0.35em \sum \limits_{l\hskip0.35em =\hskip0.35em 1}^L\sum \limits_{j\hskip0.35em =\hskip0.35em 1}^J\sum \limits_{k\hskip0.35em =\hskip0.35em 1}^{m_L}{f}^{\left(l,j\right)}\left({t}_i\right){b}_k^{\left(l,j\right)}{\boldsymbol{\phi}}_k^{\left(l,j\right)}\exp \left({\omega}_k^{\left(l,j\right)}{t}_i\right). $$
This concise definition of the mrDMD solution includes the information on the level, time-bin location, and the number of modes extracted. Figure 1 demonstrates mrDMD in terms of equation (12). In particular, each mode is represented in its respective time-bin and level.
Illustration of the multi-resolution dynamic mode decomposition (mrDMD) hierarchy. Represented are the modes
$ {\boldsymbol{\phi}}_k^{\left(l,j\right)} $
and their position in the decomposition structure. The triplet of integer values
$ l $
,
$ j $
, and
$ k $
uniquely expresses the level, bin, and mode of the decomposition.
3.2.1. Choice of parameters
The parameters in equation (10) play a relevant role in determining the effectiveness of the mrDMD algorithm in a specific application. Notice that, in the setting we consider, if the number of decomposition levels
$ L $
is determined, the number of time bins considered at each level is known. Thus, only the number of decomposition levels
$ L $
and the slow modes to consider at each level must be chosen.
To our knowledge, no principled method for choosing the parameters has been formulated so far. We follow Kutz et al. (Reference Kutz, Fu and Brunton2016b) and define general guidelines. First, notice that the choice of the parameter
$ L $
is related to the number of sampling points or snapshots available. For instance, if we know a priori that we want to compute
$ k\in \mathrm{\mathbb{N}} $
decomposition levels, and that in each time-bin at the
$ k $
-th decomposition level a certain amount
$ m\in \mathrm{\mathbb{N}} $
of sampling points is required, then
$ m{2}^{k-1} $
sampling points are needed at the first level in order to be able to reach the
$ k $
-th level of decomposition with the desired amount of snapshots in each time-bin. Thus, the number of available sampling points determines the number of decomposition levels,
$ L $
, that can be attained with the mrDMD. Another important factor to consider in the choice of the parameter
$ L $
is that the higher the decomposition level is, the higher is the frequency associated with the signal components represented by the modes computed at that level (Kutz et al., Reference Kutz, Fu and Brunton2016b). Thus, if in a given time window we are only interested in the slow modes, a high number of decomposition levels and a high snapshots’ sampling rate in that time window are unnecessary since there is no intention of extracting high-frequency temporal content. On the contrary, if we aim to analyze high-frequency components of the signal, we compute a higher number of decomposition levels so that the modes computed at the last level represent signal’s components related to higher frequencies.
Regarding the strategy to select the slow modes at each time-bin we can consider a simple threshold metric. Concretely, we define the threshold to choose DMD eigenvalues (and associated DMD modes) whose temporal behavior allows for a single wavelength or less to fit into the sampling window. As pointed out in Kutz et al. (Reference Kutz, Fu and Brunton2016b), given a specific application there could be more optimal threshold values that could be derived from the knowledge of the system analyzed. For a more detailed analysis on how the choice of these parameters affects the analysis performed by the mrDMD, see Kutz et al. (Reference Kutz, Fu and Brunton2016b).
3.3. mrDMD-based approach for damage detection
We now describe a mrDMD-based strategy to extract features that highlight the presence of a cracked tooth in a gearbox by analyzing signals representing its vibration response.
First step
Given snapshots
$ y\left({t}_i\right)\in {\mathrm{\mathbb{R}}}^1 $
for
$ i\hskip0.35em =\hskip0.35em 0,\dots, T $
, representing the temporal evolution of a sensor signal, the first step consists of constructing time-delay snapshots as follows
$$ \tilde{\mathbf{y}}\left({t}_i\right)\hskip0.35em =\hskip0.35em \left[y\left({t}_i\right),\dots, y\left({t}_{i+d}\right)\right], $$
and arranging them into matrices
$$ \overline{\mathbf{Y}}\hskip0.35em =\hskip0.35em \left[\tilde{\mathbf{y}}\left({t}_1\right)\dots, \tilde{\mathbf{y}}\left({t}_m\right)\right]\hskip1em \mathrm{and}\hskip1em Y\hskip0.35em =\hskip0.35em \left[\tilde{\mathbf{y}}\left({t}_0\right),\dots, \tilde{\mathbf{y}}\left({t}_{m-1}\right)\right]. $$
Here,
$ d $
is the length of the delay, and it has to be chosen a priori, and
$ m $
has to be chosen consequently according to the data availability.
Time-delay embedding is often applied when DMD-based methods are used to process univariate signals (Brunton et al., Reference Brunton, Brunton, Proctor, Kaiser and Kutz2017; Tirunagari et al., Reference Tirunagari, Kouchaki, Poh, Bober and Windridge2017; Clainche and Vega, Reference Clainche and Vega2018). In this work, the reasons why the time-delay embedding is used are connected to the fact that the spatial resolution of the signals we study is much lower than the temporal resolution, that is, we have several one-dimensional snapshots
$ y\left({t}_i\right)\in {\mathrm{\mathbb{R}}}^1 $
. It was first observed in Tirunagari et al. (Reference Tirunagari, Kouchaki, Poh, Bober and Windridge2017) that DMD is incapable of accurately representing even very simple one-dimensional signals. For instance, if we collect one-dimensional measurements representing the temporal evolution of a single sine wave, DMD reconstructs the dynamics using only a single real eigenvalue, which does not capture periodic oscillations (Kutz et al., Reference Kutz, Brunton, Brunton and Proctor2016a, Chapter 7). Fortunately, this issue can be addressed by considering the time delay embedding. See Kutz et al. (Reference Kutz, Brunton, Brunton and Proctor2016a) for a detailed and technical discussion on the motivations and implications of using DMD with the time-delay embedding.
Second step
The second step consists of applying the mrDMD algorithm to the time-delay snapshots to obtain the slow-modes reconstruction of the signal in the time-delay coordinates domain. The reconstruction of the time-delay snapshot
$ \tilde{\mathbf{y}}\left({t}_i\right) $
at the
$ i $
-th time step is given by equation (12) and is represented by
$ \hat{\mathbf{y}}\left({t}_i\right) $
. Since the information related to damage is assumed to be in high-frequency structures, as we will see in the next section, the general guideline is to choose the parameter
$ L $
large enough to enable the modes computed at the last decomposition level to represent the frequencies of interest.
Third step
In the third and final step, we want to obtain information about the signal’s geometrical structures related to perturbations caused by the damage we are considering, which is assumed to affect high frequencies of the signal. To do that, we first pick a time-range by selecting a time-delay snapshot, that is, the
$ i $
-th snapshot contains the evolution of the signal from time step
$ i $
to time step
$ i+d $
. After that, we subtract from the original time-delay snapshot,
$ \tilde{\mathbf{y}}\left({t}_i\right) $
, the slow modes reconstruction,
$ \hat{\mathbf{y}}\left({t}_i\right) $
, obtained using the mrDMD algorithm. After the subtraction, we are left with a residual, that is,
$$ {\mathbf{r}}_i\hskip0.35em =\hskip0.35em \tilde{\mathbf{y}}\left({t}_i\right)-\hat{\mathbf{y}}\left({t}_i\right). $$
The residual
$ {\mathbf{r}}_i $
gives us the information contained in the fastest modes computed at the deepest decomposition level for each time-bin. In the residual, anomalies related to damage are emphasized, making it possible to visually determine whether there is a cracked tooth or not in the time-range we are considering.
Algorithm 2 summarizes the main steps of the numerical procedure we just presented.
Algorithm 2. mrDMD-based approach for damage detection.
Input snapshots
$ y\left({t}_i\right)\in {\mathrm{\mathbb{R}}}^1 $
for
$ i\hskip0.35em =\hskip0.35em 0,\dots, T $
, delay’s length
$ d $
, number of decomposition levels
$ L $
.
Output Residuals
$ {\mathbf{r}}_i $
.
1: Construct time-delay snapshots equation (13) and arrange them into matrices equation (14).
2: Apply the mrDMD algorithm to the time-delay snapshots to obtain the slow-modes reconstruction of the signal in the time-delay coordinates domain. The reconstruction of the time-delay snapshot
$ \tilde{\mathbf{y}}\left({t}_i\right) $
at the
$ i $
-th time step is represented by
$ \hat{\mathbf{y}}\left({t}_i\right) $
.
3: Compute the residuals
$ {\mathbf{r}}_i\hskip0.35em =\hskip0.35em \tilde{\mathbf{y}}\left({t}_i\right)-\hat{\mathbf{y}}\left({t}_i\right) $
.
3.4. Computational cost
Algorithm 2 clearly shows that the computational cost of the proposed procedure is determined by the computational cost of the mrDMD algorithm performed in the second step, as the other steps consist of storing data and subtracting vectors. From work carried out in Kutz et al. (Reference Kutz, Fu and Brunton2016b), we know that computational efforts required by the mrDMD are dominated by the SVD of the matrix
$ \mathbf{Y} $
in the second step of Algorithm 1 and by the number of decomposition levels
$ L $
. In particular, if at a given decomposition level and time-bin we have
$ T $
time-delay snapshots with delay’s length
$ d $
, the computational cost of the SVD is
$ Q\hskip0.35em =\hskip0.35em O\left({T}^2d+{d}^3\right) $
floating point operations. Since in the mrDMD procedure, the SVD is performed at each time-bin of each decomposition level, we have that the overall cost of the mrDMD is
$ O\left({2}^LQ\right) $
. See Kutz et al. (Reference Kutz, Fu and Brunton2016b) for more details.
4. Gearbox Model
Data representing the response of gearboxes in presence of damages, such as cracks, can be either experimental using vibration measurements or numerical using simulation models. A lot of work has been conducted to analyze experimentally measured vibration signals in order to identify damages in vibrating structures (Antoniadou et al., Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015; Mohamed et al., Reference Mohamed, Sassi and Paurobally2018). The main advantage of experimental data is that they reflect the behavior of a real system. However, such data are expensive to obtain in terms of time and money, especially when repeated measurements have to be performed for different damage scenarios (Mohamed et al., Reference Mohamed, Sassi and Paurobally2018). Dynamic modeling and simulation of gearbox vibration signals can overcome these issues and can be a good alternative for studying the dynamic behavior of a gear system more simply and economically, in particular for our goal of developing a data analysis approach for early damage detection. Dynamic modeling and simulation also have the advantage of increasing the understanding of the system’s behavior before the initiation of a measurement campaign.
4.1. Theoretical setting and mathematical model
We consider a generic geared-rotor bearing system, shown in Figure 2, consisting of a rigid gearbox containing a spur gear pair mounted on flexible shafts, supported by rolling element bearings. Such a model has been previously used in other works (Kahraman and Singh, Reference Kahraman and Singh1991; Antoniadou et al., Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015) and its reliability has been shown in Kahraman and Singh (Reference Kahraman and Singh1991). We now introduce it briefly.
A generic gear rotor bearing system.
The vibration response of the gearbox model we consider can be modeled by the following dimensionless equation of motion (Antoniadou et al., Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015):
$$ {\ddot{x}}(t)+2z\dot{x}(t)+K(t)B\left(x(t)\right)\hskip0.35em =\hskip0.35em {F}_m+{F}_{te}(t)+{F}_{var}(t). $$
The coordinate
$ x(t) $
represents the difference between the dynamic transmission error and the static transmission error. The quantity
$ z $
is a dimensionless linear viscous dumping parameter.
$ {F}_m $
represents the mean force excitations, while
$ {F}_{te} $
and
$ {F}_{var} $
pertain to internal excitations related to the static transmission error and external excitations related to wind turbulence, respectively. The quantity
$ {F}_{var} $
is derived from the fluctuating component of external input torque that simulates the wind turbulence’s effect. The backlash is simulated by the function.
$$ B\left(x(t)\right):= \left\{\begin{array}{ll}x(t)-1,& \hskip4.5em \mathrm{for}\;x(t)\ge 1\\ {}0,& \mathrm{for}-1<x(t)<1\\ {}x(t)+1,& \hskip2.72em \mathrm{for}\;x(t)\le -1\end{array}\right\}, $$
while the time varying mesh stiffness is incorporated via the periodic function
$ K(t) $
. For a more detailed explanation of the specific structure of the functions mentioned and, more in general, for a deeper understanding of the gearbox model, see the Supplementary Appendix A, and Antoniadou et al. (Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015) and Kahraman and Singh (Reference Kahraman and Singh1991). Note that, this is a very simple wind turbine gearbox model, for example, it ignores the bearing vibrations and wind turbine gearboxes usually consist of three gear stages with one or two being planetary stages. According to Antoniadou et al. (Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015), it is sufficient for studies such as ours, since “more gear stages in a vibration signal would mean more frequency components at different frequency bands. Damage at a specific gear stage would therefore be shown in the vibration signal associated with the meshing frequency and its harmonics of the gear stage examined.” What is important in Antoniadou et al. (Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015) and in our study, is to investigate the influence of the varying load conditions on the vibration signals, which can be achieved.
4.2. Including damage and wind turbulence
4.2.1. Damage simulation
The damage we consider in this work is the crack of a gear’s tooth that, like other types of tooth damage, causes a reduction in the gear’s stiffness. This phenomenon can be represented by a periodic magnitude change in the mesh stiffness function
$ K(t) $
(Stander and Heyns, Reference Stander and Heyns2005). Thus, an amplitude change is periodically applied to the mesh stiffness function to simulate the cracked tooth. In the simulations used in this work, the crack of a gear’s tooth was modeled by decreasing the dimensionless mesh stiffness function by 13% of its nominal stiffness for 5 degrees of the shaft rotation, periodically, for every rotation of the damaged gear. Similar modeling was used in previous studies, see Antoniadou et al. (Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015) and Mohamed et al. (Reference Mohamed, Sassi and Paurobally2018) for a more detailed explanation. Recapturing results reported in Antoniadou et al. (Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015), we notice that in the model we consider, one should expect damage features to be evident at high-frequency components of the simulation signals because it is in the harmonics of the meshing frequency of the damaged gear pair that damage features occur. This is a phenomenon that was also observed in the experimental data used in Antoniadou et al. (Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015).
4.2.2. Wind turbulence’s simulation
The representation of the wind as a smooth flow is not realistic because it does not take into account all its irregular and stochastic fluctuations. Due to wind turbulence, wind turbines experience transient and time-varying load conditions. In order to build a realistic model, we need to be able to consider and simulate these phenomena. To do that, it is possible to use a series of wind turbine aerodynamics codes, developed by the National Renewable Energy Laboratory (Golden, CO, USA) (Bir, Reference Bir2005). Specifically, the FAST design code has been used to simulate the turbulent wind conditions, associated with wind speed of 5 and 13 m/s. See the Supplementary Appendix A for a deeper understanding of how the different wind conditions have been simulated.
4.3. Numerical data
The vibration response’s simulation of the gearbox model we considered in the previous section is given by numerical approximation of the function
$ {\ddot{x}}(t) $
obtained from the numerical solution of the dimensionless equation of motion (equation 16) (Kahraman and Singh, Reference Kahraman and Singh1991; Antoniadou et al., Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015; Mohamed et al., Reference Mohamed, Sassi and Paurobally2018). The numerical solution of equation (16) has been computed with the MATLAB
$ ode45 $
differential equation solver, with a fixed time step of 0.015. We originally intended to reproduce and use the simulated scenario as given in Antoniadou et al. (Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015). Unfortunately, the data are not publicly available and could not be provided. Moreover, using their stated model’s parameters we could not obtain similar signals to those illustrated in that work. Thus, we choose the model parameters as to get data of at least qualitatively similar behavior. The details on how the numerical model has been built are described in the Supplementary Appendix A, where also all the parameters used for the simulations are reported in Table 3 of Supplementary Appendix A.
Simulations can represent different scenarios. On the one hand, we can produce numerical signals where the fluctuations of the input torque, due to wind turbulence, are not taken into account, that is,
$ {F}_{var}(t)\hskip0.35em =\hskip0.35em 0 $
(Figure 3). These simulations represent an improbable scenario in which the load on the wind turbine blades is constant. Such signals can be helpful to test different numerical strategies in a simplified context, and we will refer to them as signals in steady load conditions. On the other hand, we can produce simulations representing the more realistic scenario where the wind is not considered as a uniform constant flow and its stochastic behavior is included in the model. These simulations represent a situation in which the load on the wind turbine blades fluctuates and has a stochastic behavior that determines random variations in the input torque. We will say that simulations generated in this scenario have varying load conditions. Figures 4 and 5 show the acceleration signals of simulations with varying load conditions produced considering two different wind speeds.
Simulation of an acceleration signal under a steady load condition and two different gearbox health scenarios. (a) Undamaged case: healthy gearbox. (b) Damaged case: one gear of the gearbox has a cracked tooth.
Simulation of an acceleration signal under varying load condition, with wind speed of 5 m/s, and two different gearbox health scenarios. (a) Undamaged case: healthy gearbox. (b) Damaged case: one gear of the gearbox has a cracked tooth.
Simulation of an acceleration signal under varying load condition, with wind speed of 13 m/s, and two different gearbox health scenarios. (a) Undamaged case: healthy gearbox. (b) Damaged case: one gear of the gearbox has a cracked tooth.
Independently of the scenario we consider, with steady or varying load conditions, we can produce signals representing a situation in which a perfectly healthy gearbox is modeled (Figures 3a, 4a, and 5a) or we can simulate a situation in which one of the gears is damaged (Figures 3b, 4b, and 5b). The simulated damage consists of a cracked tooth in one of the two gears. Recapturing some results from Antoniadou et al. (Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015), we notice that the signals simulated taking into account the wind turbulence are less smooth and show a more chaotic behavior than the signals where the load condition is steady. This is motivated by the fact that when wind turbulence is incorporated into the model it introduces a stochastic component that perturbs the acceleration signals. The higher the simulated wind condition, the noisier the signal. Note that, independently of the wind condition considered, the damaged and healthy signals do not show any relevant difference and look almost identical. This is because only high-frequency components are affected by damage and just for a short time. This fact underlines that it would be very difficult to address the damage identification task without relying on a method that can extract only the relevant information from the signals.
4.3.1. Dataset
In the following, we analyze the vibration response data of the gearbox from three gear revolutions. We use the same two different health conditions of healthy and with cracked tooth and the same two wind speeds of 5 m/s (Figure 6) and 13 m/s (Figure 7). The simulated signals are composed of 40,201 one-dimensional snapshots. In the simulations that include the damage, the presence of the cracked tooth affects the signals once per rotation at
$ {67}^{\circ } $
,
$ {427}^{\circ } $
, and
$ {787}^{\circ } $
. The MATLAB code to perform the simulations is available at https://github.com/Fraunhofer-SCAI/mrDMD-damage-detection and the simulated signal data is available at https://doi.org/10.24406/fordatis/191.
Simulation of an acceleration signal for three gear rotations under varying load condition, with wind speed of 5 m/s, and two different gearbox health scenarios. (a) Undamaged case: healthy gearbox. (b) Damaged case: one gear of the gearbox has a cracked tooth.
Simulation of an acceleration signal for three gear rotations under varying load condition, with wind speed of 13 m/s, and two different gearbox health scenarios. (a) Undamaged case: healthy gearbox. (b) Damaged case: one gear of the gearbox has a cracked tooth.
5. Results
5.1. Results with simulation data
In this section, we study the numerical strategy for early damage detection based on mrDMD proposed in Section 3. We compare our method with classical approaches used in this context, FFT as a frequency-domain approach, and TSA and EMD as a time-domain approaches.
5.1.1. Benchmark methods
5.1.1.1. Fast Fourier transform
We now apply the FFT (Rao et al., Reference Rao, Kim and Hwang2011) to the vibration signals shown in Figures 6 and 7. Figure 8a,b shows the Fourier spectra of the acceleration signals we have considered, representing the damaged and undamaged scenarios in varying load conditions associated with wind speed of 5 and 13 m/s, respectively. Independently of the wind condition considered, the differences between the healthy and cracked cases are very hard to identify just looking at the signals’ spectra. To be more specific, in Figure 8a,b, we highlight the meshing frequency and its harmonics, which are damage indicators because they relate to the frequency at which the damaged tooth has been excited (Li et al., Reference Li, Chen, Huangfu, Ma, Zhao and Yu2019). Unfortunately, also in the highlighted frequencies, there is no significant or obvious magnitude change that allows us to identify the presence of damage.
Fourier spectra of the dimensionless acceleration signals of the healthy and cracked tooth cases for varying load condition associated with wind speed of 5 m/s (a) and 13 m/s (b). The arrow-shaped cursors point to the mesh frequency,
$ {\Omega}_{mesh} $
, and its harmonics.
This shows us an arduous difficulty that frequency methods have to overcome in this scenario: when the effects of the wind turbulence are included, it is very challenging to distinguish between the spectral characteristics of the signals associated with damage and those associated with wind turbulence. In addition to that, we have that due to wind turbulence, the spectra of the signals change randomly with time, making the damage identification task even more difficult. The shown results for the comparison between the healthy and damaged signals spectra are rather basic. But note, we do not aim to show that conventional gear monitoring techniques, for example, cepstrum or kurtosis analysis, would be ineffective. These results mainly show that the operating conditions we are considering may have detrimental effects on the effectiveness of gear diagnostics techniques based on the signal’s spectral analysis.
5.1.1.2. Time synchronous averaging and Hilbert transform
TSA is a well-established procedure in rotating machines diagnostics (Randall, Reference Randall2021). One of the primary purposes of TSA analysis is to extract periodic waveforms from signals. For instance, it may be used to extract a periodic signal, such as the tooth meshing vibration of a single gear, from the whole machine’s vibration in order to detect damages by analyzing only the extracted signal’s component. Following along McFadden (Reference McFadden1986), we use TSA to demodulate the signals component associated with the gear mesh. After that, we use the HT to demodulate the extracted TSA signal in order to detect local variations caused by damage and not visible to the eye. See McFadden (Reference McFadden1986), for a more detailed description of the TSA-based strategy we are going to review.
Time synchronous averaging.
Given a signal
$ y(t)\in \mathrm{\mathbb{R}} $
, its time synchronous average
$ {y}_a(t) $
is defined as
$$ {y}_a(t)\hskip0.35em =\hskip0.35em \frac{1}{N}\sum \limits_{n\hskip0.35em =\hskip0.35em 0}^{N-1}y\left(t+{nT}_r\right), $$
where
$ N\in \mathrm{\mathbb{N}} $
and
$ {T}_r $
is the so-called periodic time. Equation (18) can be written as the convolution of
$ y(t) $
with a sequence of
$ N $
delta functions displaced by integer multiples of the periodic time
$ {T}_r $
, that is,
$$ {y}_a(t)\hskip0.35em =\hskip0.35em c(t)\ast y(t), $$
with
$$ c(t)\hskip0.35em =\hskip0.35em \frac{1}{N}\sum \limits_{n\hskip0.35em =\hskip0.35em 0}^{N-1}\delta \left(t+{nT}_r\right). $$
It can be shown that the desired signal component can be demodulated through TSA if the distance between the multiples of the periodic time is equal to its period. The main reason behind this fact is that equation (18) is equivalent in the frequency domain to the multiplication of the Fourier transform of
$ y(t) $
by a comb filter, removing all the components which fall outside the fundamental and harmonic frequencies of the desired signal. For more details, see McFadden (Reference McFadden1987).
Hilbert transform.
The HT of a signal
$ y(t) $
is defined as
$$ \hat{y}(t):= H\left[y(t)\right]\hskip0.35em =\hskip0.35em \frac{1}{\pi }{\int}_{-\infty}^{\infty}\frac{y(s)}{t-s} ds\hskip0.35em =\hskip0.35em y(t)\ast \frac{1}{\pi t}, $$
where
$ \ast $
denotes the convolution operator. Given the HT of
$ y(t) $
we can define the analytic signal
$$ z(t):= y(t)+i\hat{y}(t)\hskip0.35em =\hskip0.35em a(t){e}^{i\theta (t)}, $$
where
$ i:= \sqrt{-1} $
is the imaginary unit,
$ a(t):= \sqrt{y^2(t)+{\hat{y}}^2(t)} $
is the instantaneous amplitude, and
$ \theta (t)\hskip0.35em =\hskip0.35em \arctan \left(\frac{\hat{y}(t)}{y(t)}\right) $
is the instantaneous phase. The instantaneous frequency is defined as
$ \omega := \frac{\theta }{dt} $
. The HT does not change the domain of the variable. Indeed, the HT of a time-dependent signal is also a function of time.
We compute the HT using the “scipy.signal.function” in python.
5.1.1.3. Application of a TSA-based strategy
Figures 9a and 10a show the signals average obtained applying the TSA to acceleration signals in Figures 6 and 7, respectively. Since, our aim was to use TSA to extract the vibration signal associated with the meshing frequency, we initially set
$ {T}_r\hskip0.35em =\hskip0.35em \frac{1}{\Omega_{mesh}}\hskip0.35em =\hskip0.35em 2 $
, so that the distance between the multiples of the periodic time
$ {T}_r $
is equal to the meshing period. After that, we computed the instantaneous amplitude (Figures 9b and 10b) and phase (Figures 9c and 10c) of the TSA signals in order to detect small magnitude variations and highlight the presence of damage. As can be seen in Figures 9b,c and 10b,c, independently of the wind condition, both the instantaneous phase and amplitude associated with the healthy and cracked case show a very similar behavior. Thus, comparing the instantaneous phase and amplitude computed in the different health scenarios, we cannot emphasize any feature that indicates the presence of damage. The only exception is in Figure 10c, which shows that in the instantaneous phase associated with the cracked tooth, there is a sudden phase lag between
$ {900}^{\circ } $
and
$ 1,{000}^{\circ } $
of the gear rotation. Such a phase lag marks a substantial difference between the phase associated with healthy and cracked health conditions, and may be considered a damage indicator. However, Figure 10c shows us that phase lags also arise in the instantaneous phase associated with the healthy signal and therefore may be misleading to consider only them as damage indicators. Furthermore, shifting the analysis focus only on the instantaneous amplitude and phase associated with the cracked scenario, there is no particular anomaly arising at the angles at which the cracked tooth interacts with the other gear (
$ {67}^{\circ } $
,
$ {427}^{\circ } $
, and
$ {787}^{\circ } $
) that provides us with evidence of damage.
Results of TSA-based analysis, wind speed of 5 m/s. (a) TSA signals computed applying the TSA to the simulated acceleration signals. Instantaneous amplitude (b) and phase (c) of the computed TSA signals.
Results of TSA-based analysis, wind speed of 13 m/s. (a) TSA signals computed applying the TSA to the simulated acceleration signals. Instantaneous amplitude (b) and phase (c) of the computed TSA signals.
Notice that, we performed the above TSA-based analysis also using different values for the periodic time
$ {T}_r $
. Specifically, we set
$ {T}_r\hskip0.35em =\hskip0.35em \frac{1}{m\Omega} $
with
$ m\in \left\{1,2,\dots, 10\right\} $
, so that the distance between the multiples of the periodic time
$ {T}_r $
is equal to the periods of the meshing frequency and its harmonics. Unfortunately, independently of the parameter choice, the results were qualitatively similar to those attained by setting
$ {T}_r\hskip0.35em =\hskip0.35em \frac{1}{\Omega_{mesh}}\hskip0.35em =\hskip0.35em 2 $
. For the TSA approach, we used the MATLAB function “tsa,” which sets automatically the parameter
$ N $
in equation (18) as the maximal amount of meshing periods in the time-span considered to analyze the signal.
In this work, we investigated TSA in combination with the HT, but there are other TSA-based techniques (Randall, Reference Randall2021), for example, relying on kurtosis analysis instead of the HT, which might give different results. Here, we mainly want to show that TSA as other classical approaches may not always be effective in the operating conditions we are considering.
5.1.1.4. Empirical mode decomposition and Hilbert transform
The EMD-based strategy we present here, consists of applying EMD together with HT to detect damages, and was introduced in the wind turbine damage detection context to overcome issues experienced by studying the signals’ spectral characteristics. See Antoniadou et al. (Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015) for a more exhaustive analysis of the EMD-based fault detection method we are going to review.
Empirical mode decomposition. The essence of EMD is to decompose the signal into oscillatory functions, also called intrinsic mode functions (IMFs). Each IMF represents characteristics of the signal associated with a certain frequency band and are such that
$$ y(t)\hskip0.35em =\hskip0.35em \sum \limits_{k\hskip0.35em =\hskip0.35em 1}^N{C}_k(t)+{r}_N(t), $$
where
$ y(t)\in \mathrm{\mathbb{R}} $
is a given sensor signal,
$ {C}_k(t) $
is the
$ k $
-th IMF, and
$ {r}_N(t) $
is the rest of the approximation performed by the IMFs. To identify the desired information in the IMFs there is not a principled approach, but it is important to know that the first IMFs represent structures in the signal associated with the highest frequencies, and the ones that come after depict lower frequency components (Huang et al., Reference Huang, Shen, Long, Wu, Shih, Zheng, Yen, Tung and Liu1998). Generally, the first IMFs contain damage indicators, because they represent the highest frequency structures in the signal, which are those affected by the damage we are considering. For further details of the EMD algorithm, see Huang et al. (Reference Huang, Shen, Long, Wu, Shih, Zheng, Yen, Tung and Liu1998). We use the PyEMD python packageFootnote
1 for the following analysis.
5.1.1.5. Application of an EMD-based strategy
Figure 11 represents the IMFs obtained applying EMD on the simulated vibration signals in Figure 6. The signals in Figure 6 represent the vibration response of the modeled gearbox under varying load condition associated with a wind speed of 5 m/s.
Intrinsic mode functions of simulations under varying load condition, with wind speed of 5 m/s, and two different gearbox health scenarios. (a) Undamaged case: healthy gearbox. (b) Damaged case: one gear of the gearbox has a cracked tooth.
Looking at the IMFs in Figure 11 one can observe that in the damaged case EMD has produced more IMFs than in the healthy case. This happens because EMD perceived the effects of damage as more signal components, and therefore it added IMFs to represent them. Another explanation is that the additional IMF is related to the mode mixing problem: some IMFs include damage features that are probably already contained in another IMF.
Analyzing the IMFs in more detail, it can be seen that the IMF 1 in Figure 11b contains damage features. The damage features appear as increases in magnitude or pulses in the IMF, arising at the angles in which the cracked tooth affects the signal. Looking at the instantaneous amplitude of the first IMF (Figure 12), it can be clearly seen that at the degrees where the cracked tooth interacts with the other gear, there are sudden increases of magnitude. Such increases of magnitude allow us to identify the presence of damage and, therefore, to perform the damage identification task. Unfortunately, from Figure 12, it is also clear that together with those related to damage, there are other pulses in the IMF 1, probably related to the varying wind condition, which make the identification of damage features a very challenging task. In particular, without previous knowledge about where the cracked tooth affected the signal, we would not have been able to identify it. A possible explanation for such an event is that the mode mixing problem occurred. Thus, information related to the frequency bands affected by damage and by wind turbulence is represented in the same IMF and cannot be distinguished trivially. As a consequence of the mode mixing problem, we could also have that different IMFs contain the same frequency band information. For instance, using acceleration signals representing a scenario related to those we consider in this work, we found that sudden increases of magnitude caused by a cracked tooth could be identified in different IMFs. Thus, information related to the frequency band affected by the damage was contained in different IMFs. Such a problem was also experienced in Antoniadou et al. (Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015), where they use EMD to study signals similar to those analyzed here. Analyzing the simulated data generated with the higher wind speed of 13 m/s we could not successfully perform the damage identification task. Specifically, applying EMD to signals in Figure 7b, no IMFs could visually highlight the different nature of the effects of the wind and of the damage on the signal.
Instantaneous amplitude of IMF 1 computed in the damaged case with wind condition 5 m/s. The pulses related to the presence of damage are marked with red squares.
After this brief analysis, it is clear that one of the disadvantages of EMD is that information about the damage’s effects in the signal is not localized: EMD does not tell us which frequency-band each IMF exactly represents. Thus, we do not know in which IMF damage features can be detected. The second main weakness of EMD is the mode mixing problem, which affects the interpretability of the method and its capability to perform the fault detection task effectively. As we will see, the numerical procedure we propose overcomes these problems thanks to the characteristics of the mrDMD algorithm.
5.1.2. mrDMD-based approach for damage detection
For our numerical experiments with mrDMD, we use the simulated gearbox vibration signals shown in Figures 6 and 7 that represent the vibration response of the modeled gearbox in the damaged and undamaged scenarios under the two different wind conditions we are considering. The same data was used to test the EMD-based method previously. The length of the delay
$ d $
has been fixed to
$ d\hskip0.35em =\hskip0.35em 32,000 $
. Thus, looking at the residual
$ {\mathbf{r}}_0 $
we are analyzing a portion of the signal representing its temporal evolution for
$ {860}^{\circ } $
of the shaft rotation. The number of the decomposition levels
$ L $
in the mrDMD algorithm has been set to
$ L\hskip0.35em =\hskip0.35em 11 $
. We use the PyDMD python packageFootnote
2, where we use the previously mentioned threshold strategy for the choice of parameters.
5.1.2.1. mrDMD-residual analysis on simulation data
Figures 13 and 14 represent the results of the mrDMD-based approach that we proposed in Section 3. The first observation we make, comparing the residuals obtained from simulations produced with different wind conditions in the undamaged case (Figures 13a and 14a), is that when the wind speed is higher the high-frequency signals’ structures represented by the residuals show a more irregular behavior. This is because high-frequency structures are affected by wind turbulence, which is more relevant and volatile for higher wind speed conditions. The second observation concerns simulations representing the damaged and undamaged cases in the different wind conditions. Considering only the residual
$ {\mathbf{r}}_0 $
, the effects of the damaged tooth are clearly visible as sudden increases in magnitude arising at the angles in which the cracked tooth affects the signal (Figures 13b and 14b). The effects of the cracked tooth on the residuals can be further analyzed by considering their instantaneous amplitudes (Figures 15 and 16) computed through the HT. Independently of the wind condition considered and the residuals’ mean breadth, looking at the instantaneous amplitudes, the different characteristics of the residuals associated with the damaged and undamaged cases are even more evident. The instantaneous amplitudes related to both the healthy and damaged scenarios present a spiky structure. However, in the undamaged cases (Figures 15a and 16a), the spikes are of similar magnitude over the time-span analyzed, whereas in damaged scenarios (Figures 15b and 16b), the magnitudes of the spikes are much higher at the angles in which the cracked tooth affects the signal, making the damage identification even easier than it was considering only the residuals.
Residual of a simulation under varying load conditions with wind speed of 5 m/s, computed calculating
$ L=11 $
decomposition levels. (a) Undamaged case: healthy gearbox. (b) Damaged case: one gear of the gearbox has a cracked tooth.
Residual of a simulation under varying load conditions with wind speed of 13 m/s, computed calculating
$ L=11 $
decomposition levels. (a) Undamaged case: healthy gearbox. (b) Damaged case: one gear of the gearbox has a cracked tooth.
Instantaneous amplitude of residual
$ {\mathbf{r}}_0 $
(
$ L=11 $
), wind speed of 5 m/s. (a) Undamaged case: healthy gearbox. (b) Damaged case: one gear of the gearbox has a cracked tooth.
Instantaneous amplitude of residual
$ {\mathbf{r}}_0 $
(
$ L=11 $
), wind speed of 13 m/s. (a) Undamaged case: healthy gearbox. (b) Damaged case: one gear of the gearbox has a cracked tooth.
It can be clearly seen that the numerical procedure we propose can be effectively employed to highlight the effects of the cracked tooth on the signal, distinguishing them from the effects of the varying load condition, and enabling us to identify the damage independently of the wind condition. Moreover, a consistent advantage of the proposed strategy with respect to EMD is that the information related to damage is localized. With the proposed strategy, all the information required is contained in the residual. Note that, due to the mode mixing problem, such a residual strategy is not feasible with EMD. To effectively employ EMD to detect damages, a user search among the produced IMFs is needed to identify those structures containing information related to damage.
We remark that the choice of the number of decomposition levels
$ L $
to compute strongly affects the effectiveness of our analysis. Specifically, if the number of calculated decomposition levels is not high enough, the residual would consist of structures associated with a frequency band that would also include those low-frequencies not affected by the damage. As we will see shortly, the presence of the cracked tooth has a negligible effect on the structures associated with lower frequencies. Therefore, the residual would also consist of components that do not provide any meaningful insight for the damage identification process, which would dilute the valuable information. Consequently, the residual would not be able to highlight the presence of damage. To observe such a phenomenon, see Figure 17, which represents the residuals computed applying the proposed procedure on vibration signals in Figure 7 associated with wind speed of 13 m/s, considering a delay
$ d\hskip0.35em =\hskip0.35em 32,000 $
and calculating
$ L\hskip0.35em =\hskip0.35em 7 $
decomposition levels. It can be clearly seen that, in this case,
$ L\hskip0.35em =\hskip0.35em 7 $
decomposition levels are not enough to encode in the residual only the relevant information associated with the highest frequencies and related to damage. There is no obvious difference between the two residuals in Figure 17. Thus, they cannot be employed to perform damage detection.
Residual of a simulation under varying load conditions, with wind speed of 13 m/s, and two different gearbox health scenarios. (a) Undamaged case: healthy gearbox. (b) Damaged case: one gear of the gearbox has a cracked tooth. In each health scenario, the residual has been computed calculating
$ L=7 $
decomposition levels. Using these lower frequencies the damage cannot be detected.
To correctly choose the number of decomposition levels to compute and successfully apply the proposed procedure for damage detection, we need information about the frequencies affected by the damage we want to identify. Fortunately, in condition monitoring, the damage we want to identify and the frequencies affected by it can often be assumed to be known (Antoniadou et al., Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015).
5.1.2.2. Modal analysis with mrDMD on simulation data
We now investigate more closely the modes’ properties and how the proposed procedure works. To this aim, we analyze the modes computed via mrDMD, calculating
$ L\hskip0.35em =\hskip0.35em 11 $
decomposition levels, and used to obtain the residuals
$ {\mathbf{r}}_0 $
in Figure 14. We mainly focus on the scenario associated with the wind speed of 13 m/s because signals generated employing the same gearbox model we use and considering a similar wind condition were already simulated and studied in Antoniadou et al. (Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015). Thus, we can exploit the knowledge developed there about the signals’ physical characteristics to better interpret and understand the properties of the modes generated by mrDMD.
Recall that we can compute the residual
$ {\mathbf{r}}_0 $
associated with a signal by subtracting from the first time-delay snapshot,
$ \tilde{\mathbf{y}}\left({t}_0\right) $
, its mrDMD approximation
$ \hat{\mathbf{y}}\left({t}_0\right) $
, that is,
$ {\mathbf{r}}_0\hskip0.35em =\hskip0.35em \tilde{\mathbf{y}}\left({t}_0\right)-\hat{\mathbf{y}}\left({t}_0\right) $
. According to equation (12), if we set
$ {t}_0\hskip0.35em =\hskip0.35em 0 $
, the mrDMD approximations of the first time-delay snapshots can be explicitly written as
$$ \hat{\mathbf{y}}\left({t}_0\right)\hskip0.35em =\hskip0.35em \sum \limits_{l\hskip0.35em =\hskip0.35em 1}^L\sum \limits_{k\hskip0.35em =\hskip0.35em 1}^{m_L}{b}_k^{\left(l,1\right)}{\boldsymbol{\phi}}_k^{\left(l,1\right)}. $$
Specifically, the linear combination of the modes
$ {\phi}_k^{\left(l,1\right)} $
, computed at the first time-bin of each decomposition level, gives the mrDMD approximation of the first time-delay snapshot representing the evolution of the analyzed signal for a time-span as large as the delay. It is important to notice that to obtain the residuals
$ {\mathbf{r}}_0 $
only the modes computed at the first time-bin of each decomposition level are needed, and that each one of the modes
$ {\phi}_k^{\left(l,1\right)} $
is associated with a frequency
$ {\omega}_k^{\left(l,1\right)} $
.
In the scenario associated with wind speed of 13 m/s, we know that the wind turbulence mainly affects the signal at its highest frequencies, as it was observed in Antoniadou et al. (Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015) for the similar wind speed of 12 m/s. Figure 18 depicts three of those modes used to compute the reconstructed time-delay snapshots
$ {\hat{\mathbf{y}}}_0 $
required to obtain the residuals
$ {\mathbf{r}}_0 $
in Figure 14. The modes in Figure 18 have been computed applying the proposed analysis on vibration signal associated with wind speed of 13 m/s in the healthy and damaged scenario (Figure 7). Specifically, they are the first modes computed at the first time bins of the decomposition levels
$ l\hskip0.35em =\hskip0.35em 1,5,8 $
. Looking at Figure 18, we can make a number of observations that give us a better understanding of the modes’ properties and that generally hold for modes computed at different time bins and decomposition levels. Analyzing Figure 18, it is clear that the higher the decomposition level, the higher is the energy content
$ \mid \omega \mid $
associated with the modes computed in it. Here,
$ \omega $
is the frequency defined in Algorithm 1. This fact is empirical proof of what we already saw from a theoretical perspective in Section 3, where the mrDMD has been explained, and such a property of the modes was already introduced. Another consideration is that modes computed at the highest decomposition level are much noisier than those calculated at lower levels. This can be explained by the fact that high-energy modes incorporate the effects of the wind turbulence on the signal, which as we know from Antoniadou et al. (Reference Antoniadou, Manson, Staszewski, Barszcz and Worden2015) are more evident at high-frequencies. Thus, it further confirms that high-energy modes contain structures associated with high-frequencies, and that the residuals consist of those structures associated with the highest frequencies. Finally, we observe that in Figure 18 the modes computed in the healthy scenario and their associated energies (
$ \mid \omega \mid $
) are qualitatively the same as those calculated from signals that include the damage. This, can be motivated by the fact that the modes depicted are associated with frequencies not affected by the presence of the cracked tooth. More in general, we found that in this scenario, the effects of damage can be effectively identified and discerned from those related to wind turbulence only in those structures associated with the highest frequencies.
First modes computed at the first time bins of the decomposition levels
$ l=1,5,8 $
obtained applying the proposed analysis on simulation under varying load conditions with wind speed of 13 m/s and calculating
$ L=11 $
decomposition levels.
5.1.2.3. Automation
An advantage of the procedure we propose, compared to EMD, is that it may allow the automation of the damage detection process. The main reason for that is that it localizes the information related to damage.
Considering how EMD works, here, after a certain number of IMFs have been produced, a user made choice has to be performed to pick out all those IMFs that may contain damage features, and therefore need to be further analyzed. In general, we know that information related to the damage we are considering should be in the first IMFs, which are those associated with higher frequencies in the signal. However, this is only a rule of thumb, and because of the mode mixing problem, there is no principled way to perform the IMFs’ choice. Therefore, using EMD, there is no obvious method to automatize the damage identification process.
On the contrary, the procedure we propose localizes the information related to damage in the residual. Thus, there is no need for a modes’ selection, and the cracked tooth detection task reduces to identifying sudden increases of magnitude and anomalies within the residual. A broad spectrum of peak detection methodologies is available and should be applicable to automatize the damage identification process (Jordanov and Hall, Reference Jordanov and Hall2002; Du et al., Reference Du, Kibbe and Lin2006; Harmer et al., Reference Harmer, Howells, Sheng, Fairhurst and Deravi2008). Automatizing the damage identification process would be very advantageous in practice because it would allow to monitor the health condition of several gearboxes simultaneously and continuously, but it is out of the scope of this work.
5.1.2.4. Statistical analysis on simulation data
As a further analysis, we perform a statistical study by considering different stochastic winds and various excitation parameters. Specifically, we consider the two wind speeds of 5 and 13 m/s and for each of them 10 different shaft torque signals by using different random seeds. For each torque generated in a given wind scenario, we consider three distinct sets of internal excitation parameters (Table 4 of Supplementary Appendix A). Thus, for each wind speed, combining the 10 different torques with the three sets of internal excitation parameters, we simulated a total of 30 signals. In each of the wind scenarios we consider, we analyze the simulated acceleration signals employing the same setting we used to study the signals in Figures 6 and 7. Thus, the number of decomposition levels computed with the mrDMD is set to be
$ L\hskip0.35em =\hskip0.35em 11 $
and the length of the delay embedding is
$ d\hskip0.35em =\hskip0.35em 32,000 $
.
Figure 19 shows the results of a statistical analysis performed by applying the proposed strategy to the 30 signals generated in the damaged scenario, considering the wind speeds of 5 and 13 m/s. More precisely, Figure 19 shows the interquartile range (IQR) of the residuals’ amplitudes computed from the simulated acceleration signals and the median value of the residuals’ amplitudes at the angles where the cracked tooth interacts with the other gear. It is clearly illustrated that, independently of the wind speed, the spread of values in the residuals’ amplitudes is much higher at the angles where the cracked tooth affects the acceleration signals. Moreover, the median value of the peaks detecting the damage, that is, the red horizontal line in the figures, is clearly larger than the third quartile at the rotational angles where the damage does not occur. Thus, in both wind conditions, at least 50% of the residuals’ amplitudes show sudden increases in magnitude at the rotational degrees where the cracked tooth interacts with the other gear, hence detecting damage.
In blue, the interquartile range (IQR) of the residuals’ amplitudes computed from acceleration signals simulated under varying load conditions, with wind speed of 5 m/s (a) and 13 m/s (b). In red, the median value of the residuals’ amplitudes at the angles where the cracked tooth interacts the other gear. In gray, the a posteriori chosen threshold value to identify peaks detecting damage. For each wind speed condition, IQR and median have been calculated considering 30 acceleration signals simulated using different shaft torque signals and internal excitation parameters. The residuals have been computed calculating L = 11 decomposition levels.
This statistical analysis shows us that the results obtained on the single time series previously analyzed can be generalized to a more extensive set of simulated acceleration signals. Given sets of simulated acceleration signals associated with different wind speeds, wind fluctuations and internal excitation parameters, the proposed method can effectively detect the presence of damage in the simulated gearbox.
In applications, one would have access to acceleration measurements for several gears rotations. A practical approach based on mrDMD to detect the presence of a cracked tooth would be to consider a threshold of the residuals’ amplitudes as a value to identify damage. Suppose the degree-wise median over a certain time shows a peak at a particular rotational degree with a height clearly larger than the upper extremity of the IQR at the other rotational angles. In that case, we have that at the angle where the median’s peak arises, 50% of the analyzed residuals’ amplitudes show an increase of magnitude, thus indicating the presence of damage.
Notice that, for the setting we consider, independently of the wind speed, we can use a threshold for the peaks revealing the presence of damage smaller than the one provided by the median’s peak, but still significantly above the normal signals. For instance, in Figure 19, the gray line threshold of
$ 0.3\times {10}^{-3} $
to determine if a peak detects the existence of a cracked tooth would be effective for both wind speeds analyzed. In both Figure 19a,b, a repeatedly occurring peak height of
$ 0.3\times {10}^{-3} $
at a certain rotational degree would be a significant change of magnitude and thus indicate damage.
5.2. Results with experimental data
There is a lack of real or experimental data for the considered damage scenario under the presence of wind fluctuations. Instead, we now consider experimental data under controlled operating situations that offer other nontrivial challenges than the previously studied simulation data. Specifically, in the experimental data, the acceleration signals are mixed with other signals from different interacting gearbox components and are significantly attenuated by the dynamics of the case and other elements in the transmission path. Note that with the results on this experimental data, we want to show that the proposed mrDMD approach can overcome these additional challenging effects. We do not claim that for these controlled operating situations, our method can detect damages more effectively than other gear diagnostic techniques.
5.2.1. Experimental dataset
We consider accelerations signals from the double stage reduction gearbox described in PHMSociety (2009)Footnote 3 and shown in Figure 20. Data were sampled synchronously from accelerometers mounted on both the input and output shaft retaining plates, and were collected at 30, 35, 40, 45, and 50 Hz input shaft speed, under high and low loading. The sampling frequency is set to 66,666.67 Hz and the acquisition time is 4 s. Additionally, the runs were repeated twice for each load and speed. Both spur and helical gears are used separately in the setup of the gearbox to obtain different datasets.
Schematic of the gearbox (PHMSociety, 2009).
In this initial study, we analyze the data from the accelerometer mounted on the output shaft retaining plate of the spur gearbox. To demonstrate our method’s capability of detecting signs of damage, we compare how it operates under five different operating conditions given by different input shaft speeds and loading. The geometry of the gearbox and the operating conditions we consider are summarized in Table 1. As can be seen in Table 2, we consider two different health conditions. Either all the components of the gearbox are healthy, or there is a broken tooth together with an eccentric gear and a damaged bearing’s ball in the input shaft.
Geometry double stage reduction gearbox with spur gears set-up and its operating conditions.
Abbreviations: HL, high load condition; LL, low load condition.
Health conditions of the double stage reduction gearbox with spur gears set-up.
Abbreviations: ID, idler shaft; IS, input shaft; Is, input side; OS, output shaft; Os, output side; T, teeth.
Figure 21 shows representative acceleration signals, from the experimental dataset, generated with different load and speed conditions. Comparing the acceleration signals in the damaged and undamaged cases, the damage effects are plainly visible as sudden periodical magnitude increases in the signals. This is a characteristic shared by all the acceleration signals we analyze, independently of the operational condition considered. Notice that, the x-axis in Figure 21 are labeled with time and not degrees of rotation as previously done with simulation data. This is because in the damaged scenarios the eccentric gear on the idler shaft’s output side may modulate the speed of the gear on the output shaft (Mba et al., Reference Mba, Marchesiello, Fasana and Garibaldi2017), which is the one with the broken tooth. Consequently, in the damaged scenario, the actual angular speed of the gear on the output shaft may differ from the theoretical angular speed reported in Table 1. Since we do not know the entity of the damage, it is not possible to determine a priori degrees of rotation of the output gear for a fixed time interval.
Acceleration signals from the double stage reduction gearbox with spur gears set-up operating at different conditions. (a, b) 45 Hz input shaft speed with high load. (c, d) 50 Hz input shaft speed with low load.
For each operating condition we analyze, we consider a portion of the signal associated with three rotations of the output shaft in healthy condition plus the following 15,000 snapshots. This will ensure that when we apply our method with a time delay embedding associated with three rotations of the healthy output shaft we have 15,000 available time-delay snapshots in each scenario we consider.
5.2.2. mrDMD-based approach for damage detection on experimental data
We now apply the proposed procedure to acceleration signals from the double stage reduction gearbox with spur gears set-up operating at 30, 35, 40, 45, and 50 Hz input shaft speed, under high and low loading. Our goal is to show that the anomalies caused by the broken tooth and other damages, that affect high frequency structures of the signal, can be detected by analyzing the residuals even in a scenario in which the acceleration measurements include the dynamics of the case and other elements in the transmission path. In each operating condition, we compute
$ L\hskip0.35em =\hskip0.35em 8 $
decomposition levels, and consider a delay
$ d $
associated with three rotations of the healthy output shaft. Thus, looking at the residuals instantaneous amplitudes, we are analyzing a portion of the signal representing its temporal evolution for the time needed from the healthy output shaft for a rotation of
$ 1,{080}^{\circ } $
. Since each operating condition has a different output shaft speed, the delay
$ d $
will be given by different amounts of snapshots. The parameter choices have been summarized in Table 5 of Supplementary Appendix B.
5.2.2.1. mrDMD-residual analysis on experimental data
Figures 22–26 represent the results of the mrDMD-based approach on the experimental data. Specifically, they show the square instantaneous amplitudes of the residuals obtained by applying the proposed procedure to acceleration signals operating at 30, 35, 40, 45, and 50 Hz input shaft speed, under high and low loading.
Square instantaneous amplitudes of the residuals
$ {\mathbf{r}}_0 $
(L = 8) obtained from acceleration signals from the double stage reduction gearbox with spur gears set-up operating at 30 Hz input shaft speed with high load. (a) Undamaged case: healthy gearbox. (b) Damaged case: one gear of the gearbox has a broken tooth.
Square instantaneous amplitudes of the residuals
$ {\mathbf{r}}_0 $
(L = 8) obtained from acceleration signals from the double stage reduction gearbox with spur gears set-up operating at 35 Hz input shaft speed with high load. (a) Undamaged case: healthy gearbox. (b) Damaged case: one gear of the gearbox has a broken tooth.
Square instantaneous amplitudes of the residuals
$ {\mathbf{r}}_0 $
(L = 8) obtained from acceleration signals from the double stage reduction gearbox with spur gears set-up operating at 40 Hz input shaft speed with low load. (a) Undamaged case: healthy gearbox. (b) Damaged case: one gear of the gearbox has a broken tooth.
Square instantaneous amplitudes of the residuals
$ {\mathbf{r}}_0 $
(L = 8) obtained from acceleration signals from the double stage reduction gearbox with spur gears set-up operating at 45 Hz input shaft speed with high load. (a) Undamaged case: healthy gearbox. (b) Damaged case: one gear of the gearbox has a broken tooth.
Square instantaneous amplitudes of the residuals
$ {\mathbf{r}}_0 $
(L = 8) obtained from acceleration signals from the double stage reduction gearbox with spur gears set-up operating at 50 Hz input shaft speed with low load. (a) Undamaged case: healthy gearbox. (b) Damaged case: one gear of the gearbox has a broken tooth.
As we have previously seen, independently of the operational condition, in the acceleration signal associated with the damaged scenario, the damage effects are visible as sudden periodical magnitude increases (Figure 21). Interestingly, such behavior is reflected in the high-frequency components of the signals represented by the residuals, as it can be seen by looking at their instantaneous amplitudes (Figures 22–26). More precisely, comparing the residuals instantaneous amplitudes obtained in the damaged and undamaged scenarios in the different operating conditions the effects of damage can be clearly detected as periodical sudden increases of magnitude.
From the mrDMD analysis on simulation data, we would expect only three main increases of magnitude in the instantaneous amplitudes of the residual associated with the damaged cases. Specifically, we would expect the magnitude increases to arise when the broken tooth part of the output gear interacts with the gears mounted on the idlers’ shaft. Notice that, the residuals instantaneous amplitudes associated with the damaged scenarios show more than three prominent magnitude increases. This can be motivated by the presence of the other damages. In particular, the eccentric gear in the idler shaft’s output may modulate the speed of the gear on the output shaft, which is the one with the broken tooth. Consequently, the broken tooth part of the output gear interacts with the eccentric gear more than three times during the time-span we are considering, generating more magnitude increases.
Nevertheless, also in the more challenging scenario offered by the experimental data, the developed procedure can highlight the effects of damage on the signal’s high-frequency components, providing a tool that enables us to clearly distinguish the healthy scenarios from the damaged ones in various operating conditions.
5.2.2.2. Modal analysis with mrDMD on experimental data
For the simulated data, we performed a modal analysis where we in particular investigated the first mode for several decomposition levels. Here, for additional insight, we instead consider at each level “
$ l $
” the signal’s reconstruction,
$ {\mathbf{p}}_0^l $
provided by the modes computed at that level, that is,
$$ {\mathbf{p}}_0^l\hskip0.35em =\hskip0.35em \sum \limits_{k\hskip0.35em =\hskip0.35em 1}^{m_l}{b}_k^{\left(l,1\right)}{\boldsymbol{\phi}}_k^{\left(l,1\right)}. $$
The subscript “0” is due to the fact that we are assuming we want to analyze the first time-delay snapshot, the one computed at time
$ {t}_0 $
, as we did in the mrDMD-residual analysis. Consequently, we are interested in the modes computed in the first time-bin of each level.
Each signal component
$ {\mathbf{p}}_0^l $
summarizes information of the modes computed at level
$ l $
. Note that, each mode is associated with a specific frequency. Therefore, each signal component
$ {\mathbf{p}}_0^l $
represents information related to a specific frequency band. Once the components
$ {\mathbf{p}}_0^l $
that highlight damage features have been identified, and the frequency band or the multiple frequency bands affected by damage have been assessed, one can focus only on the signal’s structures of interest, either through the residual or by considering the components
$ {\mathbf{p}}_0^l $
associated with the frequency bands of interest.
Let us give an example. Figure 27 illustrates 10 components
$ {\mathbf{p}}_0^l $
,
$ l\hskip0.35em =\hskip0.35em 0,1,\dots, 10, $
obtained applying the mrDMD approach we propose to the acceleration signal in Figure 21b generated in high load conditions at 45 Hz input shaft speed. The delay embedding
$ d $
is the same we previously considered. In the scenario we are considering, the first four signal components can be assumed to be mainly associated with the core part of the signal, and the effects that damages have on them are negligible. Thus, the relevant signal components for our analysis are those computed at higher decomposition levels. It can be seen that the characteristics of the last two signal components are very different from the previous ones. In particular, both the components computed at the eighth and the ninth decomposition level show a much more prominent spiky behavior than the previous ones, and the magnitude of the spikes is higher as well. Moreover, after the eighth level of decomposition, there is a sudden increase in the magnitude of the computed signal components. These two factors lead us to conclude that the signal components computed at the eighth and ninth decomposition levels are associated with features related to damage. To provide additional evidence, we analyze the signal reconstruction,
$ {\mathbf{y}}_0^{1-8} $
(Figure 28), obtained using the components computed at the first eight levels, it can be seen that the sudden increases of magnitude that characterize the presence of damage in the original acceleration signal (Figure 21b) are not represented. This shows us that the features related to damage are mainly highlighted in higher frequency structures, among which there are the structures represented by the eighth and ninth signal components, that is,
$ {\mathbf{p}}_0^9 $
and
$ {\mathbf{p}}_0^{10} $
. Recall that we computed eight decomposition levels in our mrDMD-residual analysis on experimental data. Thus, the signal components
$ {\mathbf{p}}_0^9 $
, and
$ {\mathbf{p}}_0^{10} $
, were included in the residual.
Signal components
$ {\mathbf{p}}_0^l $
,
$ l=0,1,2,\dots, 10 $
, computed from acceleration signal of double stage reduction gearbox with spur gears set-up operating at 45 Hz input shaft speed with high load.
Partial reconstruction,
$ {\mathbf{y}}_0^{1-8} $
, of acceleration signal from double stage reduction gearbox with spur gears set-up operating at 45 Hz input shaft speed with high load. The partial reconstruction is obtained with the first eight signal components, that is,
$ {\mathbf{y}}_0^{1-8}={\sum}_{l=1}^8{\mathbf{p}}_0^l $
.
6. Conclusion
In this work, we proposed a numerical procedure for damage detection, which exploits DMD’s strengths of being equation-free and data-driven. From our results with simulation data, we observed that performing damage detection with Fourier analysis, TSA or EMD can be very challenging in a scenario that includes varying load conditions. We also saw that both, our mrDMD-based method and the EMD strategy, decompose the original signal into modes that can be used to extract characteristics of the signal associated with different frequencies. In particular, both methods can highlight those high-frequency structures affected by wind turbulence and change in stiffness caused by the cracked tooth and that can be effectively used to identify damages. On the one hand, we saw that the IMFs, produced by EMD, represent structures in the signal that are not associated with specific frequencies but rather with frequency bands. Moreover, the information contained in an IMF can be related to a too broad frequency band, or the same frequency information can be contained in different IMFs. These facts affect the interpretability of the IMFs and the effectiveness of EMD to perform fault detection because the information related to damage is not localized. On the other hand, mrDMD associates each mode with a specific frequency. Exploiting this fact, the proposed numerical procedure was able to identify those structures in the signal where features related to damage were highlighted. In particular, differences in the residuals’ magnitude (Figures 13 and 14) as well as in their instantaneous amplitude (Figures 15 and 16), allowed the identification of the cracked tooth visually, independently of the wind condition considered. We also saw, from the results with experimental data, that while the proposed mrDMD-based strategy is influenced by accelerations signals affected by additional vibration sources, it still is a feasible approach for damage detection.
The first strength of the proposed analysis is that, as EMD, it does not consider the signals’ spectra, avoiding all the issues related to the nondeterministic time variation of the signals spectral properties caused by the varying load condition. The other advantage is that the mrDMD algorithm is not affected by the mode mixing problem, as it is for EMD-based strategies. Moreover, with the strategy here proposed, information related to the effects of the damage in the signal is localized in the residual.
Note that, the strategy we propose can be further improved to perform an effective analysis in more challenging scenarios. Specifically, several DMD improvements could be applied in the second step of Algorithm 2 (Héas and Herzet, Reference Héas and Herzet2016; Hemati et al., Reference Hemati, Rowley, Deem and Cattafesta2017; Matsumoto and Indinger, Reference Matsumoto and Indinger2017). A particularly beneficial DMD improvement in this context was developed in Hemati et al. (Reference Hemati, Rowley, Deem and Cattafesta2017). It consists of building unbiased and noise-aware DMD by explicitly accounting for noise in the signal. Such modification could allow our method to distinguish between the noise introduced by wind turbulence, transmission path vibrations and damage effects more easily. Furthermore, additional improvements can be introduced working on the indicator function in equation (11). In general, it is possible to introduce various functional forms for the indicator function that can be used advantageously. For instance, one could consider the indicator function to take the form of different wavelet bases, that is, Coiflet, Haar, Fejér-Korovkin, and so forth.
In conclusion, this work shows the potentialities of a DMD-based strategy to perform anomaly detection, analyzing simulated signals representing the vibration response of a gearbox under varying load condition and experimental data generated with various speed and load operating conditions. The results we reported here tell us that this DMD-based strategy promises to provide high quality results by overcoming spectral and mode mixing related problems in a scenario that includes the effects of wind turbulence. Further research on mrDMD in the context of fault detection and condition monitoring of wind turbines, and related machinery, is warranted.
List of Abbreviations
- DMD
-
dynamic mode decomposition
- EMD
-
empirical mode decomposition
- FFT
-
fast Fourier transform
- HT
-
Hilbert transform
- IMF
-
intrinsic mode function
- IQR
-
interquartile range
- mrDMD
-
multi-resolution dynamic mode decomposition
- TSA
-
time synchronous averaging
Supplementary Materials
To view supplementary material for this article, please visit http://doi.org/10.1017/dce.2022.34.
Acknowledgments
We are grateful for the technical assistance of Ivan Lecei, from Fraunhofer SCAI, that supported this work in the context of data curation and providing the data produced with the FAST code. A preprint version of this article was published as Climaco et al. (Reference Climaco, Garcke and Iza-Teran2021).
Author Contributions
Conceptualization: P.C., R.I.-T.; Data curation: P.C.; Methodology: P.C.; Supervision: J.G.; Visualization: P.C., J.G.; Writing—original draft: P.C., J.G.; Writing—review and editing: all authors. All authors approved the final submitted draft.
Competing Interests
The authors declare no competing interests exist.
Data Availability Statement
The simulated signal data is available at https://doi.org/10.24406/fordatis/191.
Ethics Statement
The research meets all ethical guidelines, including adherence to the legal requirements of the study country.
Funding Statement
This work was supported by the German Federal Ministry of Education and Research (BMBF; project “MADESI” FKZ 01IS18043A).
Open access
Comments
No Comments have been published for this article.