# Reliable crack detection in a rotor system with uncertainties via advanced simulation models based on kriging and Polynomial Chaos Expansion

Jean-Jacques Sinou, Enora Denimal

► To cite this version:

Jean-Jacques Sinou, Enora Denimal. Reliable crack detection in a rotor system with uncertainties via advanced simulation models based on kriging and Polynomial Chaos Expansion. European Journal of Mechanics - A/Solids, 2022, 92, pp.1-18. 10.1016/j.euromechsol.2021.104451. hal-03425672

HAL Id: hal-03425672

https://hal.science/hal-03425672

Submitted on 11 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Reliable crack detection in a rotor system with uncertainties via advanced simulation models based on kriging and Polynomial Chaos Expansion

J.-J. Sinou<sup>a,c,*</sup>, E. Denimal<sup>b</sup>

## ARTICLE INFO

### Keywords:

Cracked rotor

Nonlinear vibration

Kriging

Polynomial chaos expansion

Hybrid meta-modelling

Uncertainties

## ABSTRACT

In this study the authors propose to take into account the nonlinear effects induced by the presence of a transverse crack to carry out vibratory monitoring and detect transverse cracks in rotating systems subject to model uncertainties. More precisely, we focus more particularly on the global complexity of the nonlinear dynamic behaviour of cracked rotors and the evolution of their harmonic components as a function of the parameters of a transverse breathing crack (its position and depth) when numerous uncertainties are considered. These random uncertainties correspond to random geometric imperfections (two disc thicknesses), random material properties (Young modulus and material density) and boundary conditions uncertainty (two bearing stiffnesses). The objective of the present work is to identify robust indicators capable of determining the presence of a crack and its status even though numerous uncertainties are present.

To conduct such a study, an advanced surrogate modelling technique based on kriging and Polynomial Chaos Expansion (PCE) is proposed for the prediction of both the critical speeds and the harmonic components $n \times$ during passage through sub-critical resonances. An extensive study to ensure the validation of the surrogate models and a relevant choice of both the parametric and random Design of Experiments (i.e. kriging DoE and PCE DoE) is proposed. The proposed methodology is applied on a flexible rotor with a transverse breathing crack and subjected to random geometric imperfections and fluctuations in material properties of the rotor system.

## 1. Introduction

Numerous works on the dynamic behaviour of cracked rotors have been performed in the past and a number of reviews for crack detection from changes in the measurement of linear and nonlinear vibrations have been published on this topic. The interested reader can refer to the following state-of-the-art papers (Wauer, 1990; Dimarogonas, 1996; Sinou, 2009b; Kumar and Rastogi, 2009; Bachschmid et al., 2010; Fan and Qiao, 2011; Bovsunovsky and Surace, 2015) for more details and comprehensive overviews on this subject.

Some approaches emphasise methodologies based on linear measurements and coupling vibration measurements of a rotating cracked shaft. The use of coupled vibrations is based on the fact that when a crack is present in a structure, the response to an excitation is not only observed in the direction of the excitation but also in other directions. For example various studies have shown the effectiveness of such a technique based on the coupling phenomena of longitudinal and bending vibrations (Papadopoulos and Dimarogonas, 1987a; Collins

et al., 1991) or bending and torsional vibrations (Papadopoulos and Dimarogonas, 1987b; Ostachowicz and Krawczuk, 1992). Studies have also shown that using coupling of longitudinal, torsional, and bending vibrations as an identification tool for cracked rotors is effective even for small cracks (Papadopoulos and Dimarogonas, 1992; Giannopoulos et al., 2015). However, it is necessary to make sure that the detected vibration coupling is not due to other undesirable phenomena such as external excitation on the rotor for example.

In addition, many developments based on the appearance and tracking of the nonlinear behaviour of cracked rotors have been proposed to allow a more efficient and robust detection or identification of cracks in complex rotor systems. The first contributions were proposed by Gasch (1976, 1993), Henry and Okah (1976) and Mayes and Davies (1976). They demonstrated that the appearances of new resonance peaks when the rotational speeds of the rotor reach $\frac{1}{n}$ of the critical speeds are

<footnote><sup>a</sup> Laboratoire de Tribologie et Dynamique des Systèmes, UMR CNRS 5513, École Centrale de Lyon, 36 avenue Guy de Collongue, 69134 Écully Cedex, France</footnote>

<footnote><sup>b</sup> Univ. Gustave Eiffel, Inria, COSYS/SII, 14S, Campus de Beaulieu, 35042 Rennes, France</footnote>

<footnote><sup>c</sup> Institut Universitaire de France, 75005 Paris, France</footnote>

<footnote>* Corresponding author at: Laboratoire de Tribologie et Dynamique des Systèmes, UMR CNRS 5513, École Centrale de Lyon, 36 avenue Guy de Collongue, 69134 Écully Cedex, France.</footnote>

<footnote>E-mail addresses: jean-jacques.sinou@ec-lyon.fr (J.-J. Sinou), enora.denimal@inria.fr (E. Denimal).</footnote>
key indicators for the detection of crack in rotating systems. Subsequently, many studies have proposed to detect damage by considering not only the nonlinear response and appearances of new sub-critical resonances, but also the evolution of each harmonic component for $n \times$ orders (Schmied and Kramer, 1984; Liao and Gasch, 1992; Darpe et al., 2004b; Sinou and Lees, 2005; Sinou, 2008a; Al-Shudeifat et al., 2010; Jun, 2012). Some experimental and numerical studies have also shown that one of the direct consequences of the evolution of the nonlinear response during passage through sub-critical resonances leads to the appearances of the loops phenomena for rotating structures (Henry and Okah, 1976; Schmied and Kramer, 1984; Adewusi and Al-Bedoor, 2002; Sinou and Lees, 2007; Sinou, 2009a; Al-Shudeifat et al., 2010; Al-Shudeifat and Butcher, 2011; Guo et al., 2017). These phenomena result in a change of the rotor orbit shape from a single loop to a double loop (triple loop, respectively) when the rotational speed passes through half (one third, respectively) of a critical speed. The double and triple loops are related to the presence of the second and third harmonic components due to damage (Schmied and Kramer, 1984; Sinou and Lees, 2007). Moreover the distortion of the orbit as well as the sizes of the inside and outside loops are drastically affected by the crack depth (Sinou, 2009a). Other research has focused on the appearances of multi-harmonics frequencies by adding periodic or impulse external excitation to detect crack in rotors (Iwasubo et al., 1992; Darpe et al., 2002, 2004a; Ishida and Inoue, 2006; Darpe, 2007; Sawicki et al., 2011; Guo et al., 2017). On the other hand, some authors investigated the evolution of the transient vibrational response of a cracked rotor according to the increasing or decreasing of the rotor speed and the use of transient signals and wavelet transform to detect the presence of damage (Sekhar and Prabhu, 1998; Adewusi and Al-Bedoor, 2001; Sekhar, 2003a,b; Zou and Chen, 2004; Sinou, 2008b; Lin and Chu, 2009). These researches show in a very convincing way that the evolution of the transient response when the rotor is passing through one-half of critical speeds can be a very efficient indicator for crack detection. This change in the vibrational responses is directly related to the evolution of the second harmonic component as a function of the rotor speed.

So nowadays, a precise monitoring of the nonlinear dynamic behaviour of a rotor and the use of the appearances of the nonlinear harmonic components are considered as one of the most efficient and robust indicator for crack detection in rotating machinery. However, in practice, vibration monitoring in rotating systems is done classically on a rotor for which certain physical parameters are not known beforehand. As a result, the detection of the presence of cracks and their propagation can be more difficult due to a lack of knowledge of the healthy rotor system and the associated dynamic behaviour. The present study therefore proposes to answer this point and to examine the potential of some indicators based on the evolution of the linear and nonlinear vibration signatures to detect a crack in a rotor. Indeed very little research has been performed to analyse the efficiency of the approaches based on a tracking process of the nonlinear behaviour of cracked rotors to detect cracks in the presence of uncertainties or lack of knowledge of the physical parameters of the system in question. The only recent studies on the subject are limited to examining and discussing the impact of uncertainties on the nonlinear dynamics of the rotor with a breathing crack for given sizes and positions (Sinou and Faverjon, 2012; Fu et al., 2018, 2019; Yongfeng et al., 2019; Fu et al., 2020). The possibility of conducting a complete study to investigate the evolution of the nonlinear dynamics of a cracked rotor as a function not only of the size and position of a crack but also in the presence of uncertainties is indeed too complex to realise with classical techniques such as Monte Carlo Simulation (MCS) because of the prohibitive computational time associated with such study. Faced with these observations, one of the original contributions of the present study is to illustrate the capacity of a previously developed approach based on meta-modelling techniques to address the possibility of investigating in details the effects of crack parameters

(position and size) on the nonlinear dynamics of a cracked rotor in the presence of modelling uncertainties. One of the novelties of the proposed research work is to propose the combination of a nonlinear method, namely the Harmonic Balance Method (HBM), with an original advanced kriging-PCE surrogate modelling technique for predicting the nonlinear responses of a cracked rotor. One of the aims of this study is also to illustrate a comprehensive engineering approach to demonstrate the effectiveness of the proposed surrogate modelling technique by providing a full convergence study for the kriging and PCE Design of Experiments (DoE) and a final validation of the full hybrid meta-model. In addition, scientific contributions in the field of crack detection in rotors are proposed: one of the major issues is not only to highlight the impact of uncertainties on the evolution of critical speeds and harmonic components but also to define which indicators, based on the linear or nonlinear vibration signature, are the most robust for the detection of a defect on a rotor system with uncertainties. To be noted that some studies have already highlighted the use of a kriging surrogate model in a context of crack parameters identification for operating rotors (Wang et al., 2018; Lu et al., 2019) or crack identification of beams and cantilever plates (Gao et al., 2012, 2013). However, in the present study, additional uncertain parameters related to uncertain physical parameters are present and must be considered for different crack parameters. If the latter are parametric, the former are usually modelled with random parameters and a Probability Density Function (PDF). Thus, the authors propose to employ a strategy based on the combination of kriging surrogate model and the Polynomial Chaos Expansion (PCE) (Denimal et al., 2018; Denimal and Sinou, 2021b), which allows to model both random and parametric uncertainties in a single surrogate model and has already proved to give satisfactory results for the study of uncertain rotors (Denimal and Sinou, 2021b). To the authors' knowledge, no study has been proposed concerning the use of both kriging surrogate models and PCE (in combination with HBM) to determine the evolution of uncertainties in the dynamical rotor system for various crack parameters and when there is probabilistic uncertainties in the rotor physical parameters. The use of such surrogate model speeds up the simulations and numerous uncertain parameters can be considered and the impact of the uncertainties on the indicators for crack detection can finally be assessed.

The paper is organised as follows: firstly, a brief description of the cracked rotor as well as the resolution of the associated equation of motion and a brief reminder of the main nonlinear dynamic characteristics of a rotor with a breathing crack are proposed. Secondly, the uncertainty propagation methodology based on a hybrid formulation combining kriging and PCE is discussed. Special attention is paid to the validation of the surrogate models and the associated convergence study of both the kriging and PCE DoEs. Finally, simulations are presented and the possibility of crack detection in the presence of uncertainties is discussed in detail.

## 2. Model of the cracked rotor

The layout of the two-bearing flexible cracked rotor is shown in Fig. 1. The system consists of a steel shaft of 0.5 m long and 0.01 m diameter with two rigid circular discs. The first disc is located at the middle of the shaft while the second disc is located at a quarter of the left side. The rotor is excited by an out-of-balance force on each disc of amplitude $m_u e_u$. The bearings are elastically supported in both horizontal and vertical directions at the ends of the shaft (bearing 1 is at the left end, and bearing 2 at the right end of the shaft). The stiffness may be different in the vertical and horizontal directions and are given by $k_{x,i}$ and $k_{y,i}$ for the ith bearing. Values of the material properties and dimensions of the rotor are given in Table 1.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Model of the cracked rotor.</figcaption></figure>

## 2.1. Equation of motion of the uncracked rotor

The uncracked shaft is modelled by 20 beam elements with four degrees of freedom (dof) per node (i.e. the horizontal and vertical displacements and the two associated rotations). The axial and torsional dof are not considered here. After assembling the shaft elements, the two rigid discs and the bearing supports, the equations of the uncracked rotor system dynamics write:

$$M\ddot{x} + (C + \omega G)\dot{x} + Kx = f + q \quad (1)$$

where the matrices **M** and **G** combine the mass and gyroscopic contributions of the shaft and the two rigid discs, respectively, such as:

$$M = M_s + M_{d1} + M_{d2} \quad (2)$$

$$G = G_s + G_{d1} + G_{d2} \quad (3)$$

where the subscripts **s**, **d1**, **d2**, refer to the contributions of the rotor shaft and the first and second discs, respectively. Of course, the matrices **M**<sub>d1</sub>, **M**<sub>d2</sub>, **G**<sub>d1</sub> and **G**<sub>d2</sub> are non-null only at the dof associated with the first and second discs, respectively. The matrix **K** includes the stiffness matrices of the shaft and the supports, together with the circulatory matrix which accounts for the internal damping of the shaft, such as:

$$K = K_s + K_{b1} + K_{b2} \quad (4)$$

where the subscripts **s**, **b1** and **b2** refer to the contributions of the rotor shaft and the first and second bearing supports, respectively. Here again the matrices **K**<sub>b1</sub> and **K**<sub>b2</sub> are non-null only at the dof associated with the first and second bearing supports, respectively. The vector **q** contains the gravitational forces. The unbalance forces are located at the position of the two rigid discs and are denoted by the vector **f**. To be noted that we apply a classical Rayleigh damping for the shaft. This viscous damping model expresses the damping matrix **C** as a linear combination of the mass **M**<sub>s</sub> and stiffness **K**<sub>s</sub> matrices for the rotor shaft (i.e. **C** = α**M**<sub>s</sub> + β**K**<sub>s</sub> where α and β are constants of proportionality). The two reference vibration modes selected for the calculation of the two proportional damping coefficients are the first and fourth modes of the rotor system (i.e. the first backward mode and the second forward mode) with a damping ratio of 0.5%.

<table><caption>Table 1
Parametric values of the rotor system.</caption><thead><tr><th>Notation</th><th>Description</th><th>Value</th></tr></thead><tbody><tr><td><i>R</i></td><td>Radius of the shaft</td><td>0.005 m</td></tr><tr><td><i>L</i></td><td>Length of the shaft</td><td>0.5 m</td></tr><tr><td><i>ρ</i></td><td>Density</td><td>7800 kg/m³</td></tr><tr><td><i>E</i></td><td>Young's modulus of elasticity</td><td>2 × 10<sup>11</sup> N/m²</td></tr><tr><td><i>ν</i></td><td>Poisson ratio</td><td>0.3</td></tr><tr><td><i>α</i></td><td>First Rayleigh damping coefficient</td><td>0.52</td></tr><tr><td><i>β</i></td><td>Second Rayleigh damping coefficient</td><td>1 × 10<sup>-6</sup></td></tr><tr><td><i>R</i><sub>1</sub></td><td>Outer radius of disc 1</td><td>0.025 m</td></tr><tr><td><i>R</i><sub>2</sub></td><td>Outer radius of disc 2</td><td>0.01 m</td></tr><tr><td><i>t</i><sub>1</sub></td><td>Thickness of disc 1</td><td>0.015 m</td></tr><tr><td><i>t</i><sub>2</sub></td><td>Thickness of disc 2</td><td>0.015 m</td></tr><tr><td><i>ρ</i><sub>1</sub></td><td>Density of disc 1</td><td>7800 kg/m³</td></tr><tr><td><i>ρ</i><sub>2</sub></td><td>Density of disc 2</td><td>7800 kg/m³</td></tr><tr><td><i>m</i><sub>u</sub></td><td>Mass unbalance</td><td>0.001 kg</td></tr><tr><td><i>e</i><sub>u</sub></td><td>Eccentricity of the mass unbalance</td><td>0.01 m</td></tr><tr><td><i>k</i><sub>x,1</sub></td><td>Stiffness of the left support (vertical direction)</td><td>5 × 10<sup>5</sup> N/m</td></tr><tr><td><i>k</i><sub>y,1</sub></td><td>Stiffness of the left support (horizontal direction)</td><td>7.5 × 10<sup>5</sup> N/m</td></tr><tr><td><i>k</i><sub>x,2</sub></td><td>Stiffness of the right support (vertical direction)</td><td>5 × 10<sup>5</sup> N/m</td></tr><tr><td><i>k</i><sub>y,2</sub></td><td>Stiffness of the right support (horizontal direction)</td><td>5 × 10<sup>5</sup> N/m</td></tr></tbody></table>

## 2.2. Modelling of the cracked element and the breathing mechanism

The modelling of the transverse crack is based on the variation of the stiffness properties of the crack cross section, as proposed by Davies and Mayes (1984) and Mayes and Davies (1984). This model considers that a transverse crack introduces local flexibility due to strain energy concentration in the vicinity of the tip of the crack under load. The reduction of the second moment of area $\Delta I$ of the element that contains the transverse crack can be defined by

$$\Delta I = I_0 \left( \frac{\frac{R}{l} (1 - v^2) F(\mu)}{1 + \frac{R}{l} (1 - v^2) F(\mu)} \right) \quad (5)$$

where $I_0$, $R$, $l$, and $v$ are the second moments of area, beam radius, length of the section and Poisson's ratio, respectively. $\mu$ corresponds to the non-dimensional crack depth defined by $\mu = \frac{h}{R}$ where $h$ is the crack depth of the beam, as illustrated in Fig. 1. $F(\mu)$ defines the nonlinear compliance as a function of variations in non-dimensional crack depth $\mu$. Davies and Mayes (1984), Mayes and Davies (1984) proposed to obtain the evolution of $F(\mu)$ through a series of experiments with chordal cracks. The local flexibility due to the presence of one crack leads to
an additional contribution stiffness denoted by $\mathbf{K}_{crack}$ at the crack's location. For the interested reader, the whole process to obtain the complete expressions of the stiffness matrix $\mathbf{K}_{crack}$ are given in Sinou and Lees (2005).

Assuming that the breathing mechanism of a crack in rotor is mainly due to the shaft weight, two well-known models of cracks are classically used, namely the switching and breathing crack models. The switching model (also known as hinge model) considers that, based on the consideration of the static deflection of the rotor, the stiffness of the rotor switches directly from the fully closed crack (i.e. uncracked shaft) to the fully open crack during one periodic rotation of the shaft. One more realistic crack breathing model is defined by the fact that the stiffness variations are synchronous with the rotational angles of the rotor. In other words the breathing crack model closely imitates the breathing behaviour of a real crack. By assuming that the gravity force is much greater than the imbalance force, we use in this paper one of the classical form describing the breathing crack (Gasch, 1976) (i.e. the periodic opening and closing of the crack when a cracked rotor rotates). In the present work, the breathing crack function is approximated by a cosine function $g(t)$ (Mayes and Davies, 1984; Sinou and Lees, 2005)

$$g(t) = \frac{1 - \cos\omega t}{2} \qquad (6)$$

where $\omega$ defines the rotational speed of the rotor. For $g(t) = 1$, the crack is fully open and for $g(t) = 0$, the crack is totally closed and the cracked rotor stiffness is equal to the uncracked rotor stiffness (i.e. no effect due to the crack on the dynamic behaviour of the rotor system). Of course, more precise models of the breathing mechanism exist and their impacts on linear and nonlinear vibrations of cracked rotors have been proposed by researchers (Bachschmid et al., 2000, 2008; Al-Shudeifat et al., 2010; Al-Shudeifat and Butcher, 2011; Spagnol et al., 2020; Ganguly and Roy, 2021). In this study we assume that the well-known breathing model defined in Eq. (6) is sufficient to approximate the nonlinear phenomena of a cracked rotor as a function of the crack size and position. Furthermore, the proposed advanced surrogate modelling technique can be adapted to other crack breathing models without difficulties or loss of generality in the proposed future conclusions on the impact of uncertainties on the crack detection indicators.

## 2.3. General equation of motion of the cracked rotor and its resolution

The equations of the rotor system with a breathing crack writes

$$M\ddot{\mathbf{x}} + (\mathbf{C} + \omega\mathbf{G})\dot{\mathbf{x}} + (\mathbf{K} - g(t)\mathbf{K}_c)\mathbf{x} = \mathbf{f} + \mathbf{q} \qquad (7)$$

where the global stiffness matrix $\mathbf{K}_c$ due to the presence of the crack situated at the ith beam location is given by

$$\operatorname{diag}(\mathbf{K}_c) = \begin{array}{c} (\mathbf{O} \quad \cdots \quad \mathbf{O} \quad \mathbf{K}_{\text{crack}} \quad \mathbf{O} \quad \cdots \quad \mathbf{O}) \\ \uparrow \\ \text{i-th element} \end{array} \qquad (8)$$

where $\mathbf{O}$ defines the 8 × 8 null matrix.

Eq. (7) has time-dependent coefficients (i.e. parametric terms) due to the breathing behaviour of the crack when the system rotates. A very classical way to determine the periodic solutions of the cracked rotor system is to assume that the dynamical responses can be approximated by a truncated Fourier series of order $m$ with a fundamental frequency (Sinou and Faverjon, 2012)

$$\mathbf{x}(t) = \mathbf{A}_0 + \sum_{k=1}^{m} (\mathbf{A}_k \cos(k\omega t) + \mathbf{B}_k \sin(k\omega t)) \qquad (9)$$

where $\mathbf{A}_0$, $\mathbf{A}_k$ and $\mathbf{B}_k$ (with $k = 1, ..., m$) define the unknown coefficients of the finite Fourier series. $\omega$ is the rotational speed of the rotor system. The number of harmonic coefficients $m$ is selected on the basis of the number of significant harmonics expected in the dynamical response. In general, three or four harmonic components are sufficient to describe the dynamics of the cracked rotor (Sinou and Faverjon, 2012).

Then, the calculation of periodic solutions of Eq. (7) can be determined by solving the problem in the frequency domain instead of the time domain. The Harmonic Balance Method (HBM) is one of the most popular frequency domain method for computing periodic solutions of ordinary differential equations. The rewriting of Eq. (7) in the frequency domain gives a set of $(2m+1) \times n$ linear equations (where $n$ is the dof number of the rotor system)

$$(\mathbf{\Lambda} - \mathbf{\Lambda}_c) \mathbf{\Theta} = \mathbf{\Gamma} \qquad (10)$$

with $\mathbf{\Theta}$ the unknown coefficients of the finite Fourier series for the periodic solution such as

$$\mathbf{\Theta} = [\mathbf{A}_0 \quad \mathbf{A}_1 \quad \mathbf{B}_1 \quad \cdots \quad \mathbf{A}_k \quad \mathbf{B}_k \quad \cdots \quad \mathbf{A}_m \quad \mathbf{B}_m]^T \qquad (11)$$

To be noted that the unknown Fourier coefficients $\mathbf{\Theta}$ can be found directly by solving Eq. (10) due to the fact that the initial system in the time domain corresponds to a system with time-dependent coefficients (see Eq. (7)). $\mathbf{\Gamma}$ defines the contribution of the gravitational and unbalance forces. The vector $\mathbf{\Gamma}$ is given by $\mathbf{\Gamma} = [\mathbf{C}_0^f \quad \mathbf{C}_1^q \quad \mathbf{S}_1^q \quad \mathbf{0} \quad \cdots \quad \mathbf{0}]^T$ due to the fact that the gravitational and unbalance forces can be defined by finite Fourier series with constant components (i.e. $\mathbf{f} = \mathbf{C}_0^f$) and first-order periodic components (i.e. $\mathbf{q} = \mathbf{C}_1^q \cos(\omega t) + \mathbf{S}_1^q \sin(\omega t)$) in the frequency domain, respectively. $\mathbf{\Lambda}$ corresponds to the contributions of the mass, damping and stiffness matrices of the uncracked rotor system such as

$$\mathbf{\Lambda} = \operatorname{diag}(\mathbf{K} \quad \mathbf{A}_1 \quad \cdots \quad \mathbf{A}_k \quad \cdots \quad \mathbf{A}_m) \quad \text{with} \\ \mathbf{\Lambda}_k = \begin{bmatrix} \mathbf{K} - k^2 \omega^2 \mathbf{M} & k\omega \mathbf{D} \\ -k\omega \mathbf{D} & \mathbf{K} - k^2 \omega^2 \mathbf{M} \end{bmatrix} \qquad (12)$$

where $\mathbf{D}$ is the global damping matrix including both the damping matrix $\mathbf{C}$ and the gyroscopic matrix $\mathbf{G}$. The matrix $\mathbf{\Lambda}_c$ defines the contribution of the parametric terms $g(t) \mathbf{K}_c$ due to the presence of the breathing crack. $\mathbf{\Lambda}_c$ can be written in the previously predefined base given in Eq. (11) in the form:

$$\mathbf{\Lambda}_c = \frac{1}{4} \begin{bmatrix} 2\mathbf{K}_c & -\mathbf{K}_c & & & \\ -2\mathbf{K}_c & 2\mathbf{K}_c & \mathbf{0} & -\mathbf{K}_c & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & 2\mathbf{K}_c & \mathbf{0} & -\mathbf{K}_c \\ & & \ddots & & \ddots \\ & & & -\mathbf{K}_c & \mathbf{0} & 2\mathbf{K}_c & \mathbf{0} \\ & & & \mathbf{0} & -\mathbf{K}_c & \mathbf{0} & 2\mathbf{K}_c \\ & & & & \ddots & & \ddots \\ & & & & & -\mathbf{K}_c & \mathbf{0} & 2\mathbf{K}_c & \mathbf{0} \\ & & & & & \mathbf{0} & -\mathbf{K}_c & \mathbf{0} & 2\mathbf{K}_c \end{bmatrix} \qquad (13)$$

Considering this last expression, it can be noted that the crack leads to direct interactions between the static deflection (i.e. $\mathbf{A}_0$) and some contribution of the first order (i.e. $\mathbf{A}_1$) and the second order (i.e. $\mathbf{A}_2$) due to the presence of the gravitational forces.

## 2.4. Preamble — deterministic case

Before studying the problem of crack detection by investigating the effects of the two parameters of the crack (i.e. the crack depth $\mu$ and the crack location $L_{crack}$) and uncertainties in the rotor system parameters, a brief summary is given of the main nonlinear dynamic characteristics of a rotor with breathing crack.

Fig. 2 shows the vertical nonlinear response for the cracked rotor system according to the rotation speed of the cracked rotor, with a breathing crack of depth $\mu = 1$ and positioned at $L_{crack}=0.3$ m. Firstly, it is observed that the response is composed by not only the first harmonic component but also by the components of higher orders (in the present case only the first four orders are shown in Fig. 2). The interaction between a rotating and breathing crack, gravity and the presence of unbalances leads to the occurrence of super harmonics of
<figure><img src="image_3.png" /><figcaption>Fig. 2. Nonlinear vertical amplitudes at 0.1 m from the left side of the rotor system with a breathing crack ($\mu = 1$ and $L_{crack} = 0.3$ m from the left side) (black = order 1, red = order 2, magenta = order 3, blue = order 4).</figcaption></figure>

<table><caption>Table 2
Critical speeds of the uncracked rotor system and the rotor with a breathing crack ($\mu = 1$, $L_{crack}=0.3$ m).</caption><thead><tr><th>Case</th><th>f<sub>1</sub> (Hz)<br>(1st backward mode)</th><th>f<sub>2</sub> (Hz)<br>(1st forward mode)</th><th>f<sub>3</sub> (Hz)<br>(2nd backward mode)</th><th>f<sub>4</sub> (Hz)<br>(2nd forward mode)</th></tr></thead><tbody><tr><td>Uncracked rotor</td><td>49.0</td><td>49.2</td><td>252.3</td><td>266.1</td></tr><tr><td>Cracked rotor</td><td>48.1</td><td>49.0</td><td>250.4</td><td>265.3</td></tr></tbody></table>

order $j \times$ (with $j = 2, \dots, 4$) leading to amplitude peaks during rotation speeds equal to approximately $\frac{1}{j-n}$ (with $n = 0, \dots, j-1$) of the critical speeds of the cracked rotor. For example we can clearly observe two peaks corresponding to passages of sub-critical resonance peaks $\frac{1}{2}$ and $\frac{1}{3}$ for the first bending mode in the vertical direction (see the marks 3 and 2 in Fig. 2, respectively). Likewise, other peaks can also be seen (see the marks 5, 6 and 7), corresponding to sub-critical resonance peaks $\frac{1}{2}$ and $\frac{1}{3}$ for the second backward and forward bending modes. The presence of peaks at $\frac{1}{4}$ of the critical speeds are also visible even if the associated amplitudes remain low (see the marks 1 and 4). Amplitude peaks of orders $j \times$ (with $j = 2, \dots, 4$) are also present at the passages of the critical speeds of the cracked rotor (see the marks 8, 9 and 10). It is recalled that the response of the rotor without a breathing crack is only composed of the first harmonic component (i.e. the vibratory response of the healthy rotor is linear).

Secondly, the presence of the crack induces a small decrease in the critical speeds of the rotor system (denoted by $f_i$ for $i = 1, \dots, 4$) due to the reduction in system stiffness, as indicated in Table 2. For the interested reader, the associated natural frequencies of the rotor at rest are given by $f_1^{\omega=0} = 49$ Hz, $f_2^{\omega=0} = 49.2$ Hz, $f_3^{\omega=0} = 256.3$ Hz and $f_4^{\omega=0} = 262.3$ Hz. Thus, by comparing these values to those of the critical speeds given in Table 2, it is clearly shown that the gyroscopic effects are noticeable for the second backward and forward modes. This example allows to put forward some conclusions classically admitted for the detection of cracks in rotors. The critical speed differences between cracked and uncracked rotor configurations are generally too small to effectively detect the presence of a crack on a rotating system, even for a crack with a significant depth. A breathing crack generates prominent frequency components of $2 \times$ and $3 \times$ near to one-half and one-third of the critical speeds, which are well known indicators for detecting a transverse crack in a rotating system.

Nevertheless the crack depth and location can drastically affect both the critical speeds of the cracked rotor and the vibrational amplitudes of both the critical speeds and the sub-critical resonances. To briefly

illustrate this point, Fig. 3 shows the dependence of these outputs by drawing the evolution surface of each of these quantities as a function of the position and size of the crack (for an unbalance given in terms of intensity and dephasing in relation to the front of the crack). The figures in the left column in Fig. 3 illustrate this reduction in critical speeds associated with the first four modes. An increase in the size of the crack leads to a reduction of the critical speeds. The closer the crack is to a vibration node, the more the critical speed of the associated mode draws closer to the critical speed of a healthy rotor. This result is also valid for sub-critical resonances (not presented here for the sake of brevity). Other figures in Fig. 3 provide the results and trends for the evolutions of amplitude peaks of orders $j \times$ at $\frac{1}{j}$ (with $j = 1, \dots, 3$) for the first four critical speeds. Noted that the visualisation of the horizontal vibration responses (the vertical vibration responses, respectively) has been chosen for the backward modes $f_1$ and $f_3$ (the forward modes $f_2$ and $f_4$). This choice is guided by the fact that these configurations correspond to the largest amplitudes in the vertical and horizontal directions for a specific quantity of interest and will be kept for the rest of the study. The obtaining of the additional directions is straightforward and results presented here could be easily extended. It can be observed that the differences of evolution of amplitudes $j \times$ (with $j = 1, \dots, 3$) for the first two forward and backward modes are characteristic of interactions between the crack (position and depth), the unbalance and gravity. The presence of a crack introduces dissimilar contributions according to the modes due to the evolution of the centroid of the cracked element. The mechanism associated with the dynamic behaviour of the crack is in fact directly linked to the global nonlinear dynamic response (which depends on the unbalance and gravity), leading to dissymmetry in the breathing effects of the crack in the horizontal and vertical directions and therefore differences in amplitudes when passing through critical speeds and sub-critical resonances $\frac{1}{j}$. For the rest of the study, the amplitude at one critical speed will be denoted by $X_{1\times}^{f_i}$ where $i$ corresponds to the $i$th critical speed (with $i = 1, \dots, 4$) and the amplitude at one sub-critical resonance
will be denoted by $X_{j×}^{f_i}$ where i corresponds to the ith critical speed (with i = 1, ..., 4) and j corresponds to the jth order (with j = 2, 3).

However, caution is required regarding these initial observations, since the presence of uncertainties can make the detection of defects based on the evolution of the critical speeds and the nonlinear vibratory behaviour more delicate. Therefore, the major objective of the rest of the study is to be able to provide an answer to this open question which has never been addressed before.

Before going further in the study it is important to point out that the computation time via the HBM can be very costly as the critical speeds $f_i$ and the maximum amplitudes $X_{j×}^{f_i}$ (for i = 1, ..., 4 and j = 1, ..., 3) must be identified accurately for the rest of the study. Therefore, a classical HBM computation over the whole rotational speed interval of interest must be conducted with a very small rotational speed step, which represents a substantial number of simulations and so the major drawback for this approach. As an illustration, it requires more than 10 min for a reference calculation in order to have sufficient accuracy on the desired outputs. So, an iterative process with dichotomy over twelve pre-selected rotational speed intervals (i.e. one interval around each critical speed and also each sub-critical resonances $\frac{1}{j}$ with j = 2, 3) has been implemented to accurately identify the first four critical speeds $f_i$ (with i = 1, ..., 4), the associated maximum amplitudes $X_{j×}^{f_i}$ for each critical speed, and the maximum amplitudes $X_{j×}^{f_i}$ at sub-critical resonances (where j corresponds to the jth order, j = 2, 3). This dichotomy procedure reduces drastically the calculation time for each reference calculation as the rotational speed ranges of interest for each of the quantities sought are reduced, giving a final computational time of about 30 s.

# 3. Uncertainty propagation methodology

## 3.1. Uncertain parameters

In the present study, six parameters describing the rotor system are considered as uncertain: the Young modulus and density of the shaft (E and $\rho$), the thickness of the two discs ($t_1$ and $t_2$) and the vertical and horizontal stiffnesses of the right bearing support ($k_{x,2}$ and $k_{y,2}$). Their characteristics and their contributions on the different matrices of the rotor system are given in Table 3. They are all modelled with a uniform PDF and vary up to 5% around their mean value for E and $\rho$, and 10% around their mean value for $t_1$, $t_2$, $k_{x,2}$ and $k_{y,2}$. The choice of the uncertain parameters is guided by the fact that the matrices impacted by uncertainties are different according to the uncertain parameters, making a general framework to test the validity of our approach. To be noted that the damping matrix C that is defined by $C = \alpha M_s + \beta K_s$ is still estimated so that the calculation of the two proportional damping coefficients $\alpha$ and $\beta$ lead to a damping ratio of 0.5% for the first and fourth modes of the rotor system. So, even if taking uncertainties into account can modify the contributions of $M_s$ and $K_s$ and so the damping matrix C, we choose to keep a constant damping ratio whatever the origin or the value of the uncertain parameters.

In addition, a parametric variation of both the crack depth $\mu$ and crack location $L_{crack}$ are considered. The crack depth $\mu$ and the crack location $L_{crack}$ can take values in [0.05; 1] (as a reminder $\mu$ = 1 corresponds to the loss of half the shaft area at the crack location) and in [0; L], respectively. The final objective is to assess the evolutions of the first and second forward and backward critical speeds, denoted $f_1$, $f_2$, $f_3$ and $f_4$, and of the nonlinear amplitudes $X_{j×}^{f_i}$ (with i = 1, ..., 4 and j = 1, ..., 3) according to the position and depth of the crack and in the presence of uncertainties. The main interest is to be able to decide on the possibility of considering the evolution of nonlinear signatures for the detection of a breathing crack on a rotor with uncertainties. The different critical speeds $f_i$, i ∈ [1, 4] and the nonlinear amplitudes $X_{j×}^{f_i}$, i ∈ [1, 4], j ∈ [1, 3] depend on the crack parameters (location $L_{crack}$ and depth $\mu$), that are parametric and grouped in the vector $\mathbf{x} =$

$[L_{crack}, \mu]$, and they also depend on the random uncertain parameters grouped in the vector $\xi = [\xi_E, \xi_\rho, \xi_f, \xi_{t_2}, \xi_{k_{x,2}}, \xi_{k_{y,2}}]$. To study the impact of the uncertain parameters of the model on the signature of the cracked rotor, one could consider doing Monte Carlo Simulation to predict the evolutions of the critical speeds $f_i(\mathbf{x}, \xi)$ and the vibration amplitudes $X_{j×}^{f_i}(\mathbf{x}, \xi)$. This would however represents an unaffordable computation cost in practice, and so surrogate modelling techniques are preferred. It builds a meta-model (also called surrogate model) from a limited number of runs of the reference model. Thus the objective is to approximate each critical speed $f_i(\mathbf{x}, \xi)$ and each amplitude $X_{j×}^{f_i}(\mathbf{x}, \xi)$ by the surrogate models $\hat{f}_i(\mathbf{x}, \xi)$ and $\hat{X}_{j×}^{f_i}(\mathbf{x}, \xi)$, respectively. These surrogate models will be built from a hybrid formulation coupling kriging and PCE. Their general mathematical background is briefly reminded to the reader for the sake of concision. For an extended description of each method, the reader could refer to Wiener (1938), Jones et al. (1998), Lophaven et al. (2002), Rasmussen (2003) and Blatman and Sudret (2011). The hybrid formulation has been introduced first in Denimal et al. (2018).

## 3.2. Kriging mathematical background

For parametric uncertainties, surface response methods proved to be efficient in terms of numerical cost. The best unbiased predictor is the kriging formulation (Jones et al., 1998), which has proven its efficiency for many nonlinear problems. Let $\alpha$ denotes the function to be approximated of argument $\mathbf{x} \in \mathbb{R}^n$ (here $n = 2$ and $\mathbf{x} = [L_{crack}, \mu]$). A kriging approximation $\hat{\alpha}$ of $\alpha$ writes:

$$ \hat{\alpha}(\mathbf{x}) = \sum_{k=1}^{q} \beta_k g_k(\mathbf{x}) + Z(\mathbf{x}) \quad (14) $$

where the ($g_k$), $k \in [1, q]$ are q real known regression functions and $\beta$ are the regression coefficients to be determined. The functions are chosen a priori and often taken as polynomial functions of low order (up to order 2). Z is a zero-mean Gaussian process of unknown variance $\sigma^2$. Its covariance function is $C(\mathbf{x}, \mathbf{x}') = \sigma^2 R(\theta, \mathbf{x}, \mathbf{x}')$, where $\mathbf{x}$ and $\mathbf{x}'$ are two points in $\mathbb{R}^n$, and R the correlation function. The latter is usually unknown and constructed as the tensor product of a 1D-kernel function k:

$$ R(\theta, \mathbf{x}, \mathbf{x}') = \prod_{j=1}^{n} k(\theta_j, d_j) \quad (15) $$

with $\theta_j$ the value of $\theta$ in dimension j, and $d_j = x_j - x'_j$ the distance between $\mathbf{x}$ and $\mathbf{x}'$ in dimension j. For anisotropic kriging, $\theta_j$ is different in each dimension and must be computed. Classical kernel functions are the Gaussian kernel, the exponential kernel and the different Matèrn 5/2 kernels (Rasmussen, 2003).

When constructing a kriging surrogate model, N evaluations of the model are required, i.e. a set of inputs $\{\mathbf{x}^{(1)}, \dots, \mathbf{x}^{(N)}\}$ and their evaluations $y = \{\alpha(\mathbf{x}^{(1)}, \dots, \alpha(\mathbf{x}^{(N)}))\}$. This set of inputs and outputs is usually called *Design of Experiments* (DoE). If G denotes the regression matrix of coefficients $G_{ij} = g_j(\mathbf{x}^{(i)})$ and R the correlation matrix of coefficients $R_{ij} = R(\theta, \mathbf{x}^{(i)}, \mathbf{x}^{(j)})$, then $\theta$ is determined by solving a maximum likelihood problem and $\beta$ and $\sigma^2$ are directly deduced (Jones et al., 1998; Lophaven et al., 2002). The kriging predictor at a new point $\mathbf{x}_0$ writes:

$$ \hat{\alpha}(\mathbf{x}_0) = g_0^T \beta + r_0^T \mathbf{R}^{-1} (\mathbf{y} - \mathbf{G}\beta) \quad (16) $$

with $g_0$ the vector of $g_j(\mathbf{x}_0)$ and $r_0$ the vector of $R(\theta, \mathbf{x}^{(i)}, \mathbf{x}_0)$. For an extensive description of the kriging formulation and creation, the interested reader may refer to Jones et al. (1998) and Lophaven et al. (2002).

## 3.3. Polynomial Chaos mathematical background

Let $Y(\xi)$ be a second-order random parameter function, which depends on $\xi$ a vector of m independent random variables $\xi_j$, $j \in [1, m]$ (here
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Evolution of the deterministic critical speeds $f_i$ (with $i = 1, ..., 4$), amplitudes $X_{1\times}^{f_i}$ at the critical speeds (with $i = 1, ..., 4$) and amplitudes $X_{j\times}^{f_i}$ at sub-critical resonances $n\times$ (with $i = 1, ..., 4$ and $j = 2, ..., 3$) for the first two forward and backward modes for different crack locations and depths.</figcaption></figure>

<table><caption>**Table 3**
Random parameters description.</caption><thead><tr><th>Parameter</th><th>Notation</th><th>% variation</th><th>Law</th><th>Contribution</th></tr></thead><tbody><tr><td>Young modulus - shaft</td><td><i>E</i></td><td>±5%</td><td>Uniform</td><td><b>K</b><sub>s</sub></td></tr><tr><td>Density - shaft</td><td><i>ρ</i></td><td>±5%</td><td>Uniform</td><td><b>M</b><sub>s</sub> and <b>q</b></td></tr><tr><td>Thickness - disc 1</td><td><i>t</i><sub>1</sub></td><td>±10%</td><td>Uniform</td><td><b>G</b><sub>d1</sub>, <b>M</b><sub>d1</sub> and <b>q</b></td></tr><tr><td>Thickness - disc 2</td><td><i>t</i><sub>2</sub></td><td>±10%</td><td>Uniform</td><td><b>G</b><sub>d2</sub>, <b>M</b><sub>d2</sub> and <b>q</b></td></tr><tr><td>Stiffness - right support (vertical direction)</td><td><i>k</i><sub>x,2</sub></td><td>±10%</td><td>Uniform</td><td><b>K</b><sub>b2</sub></td></tr><tr><td>Stiffness - right support (horizontal direction)</td><td><i>k</i><sub>y,2</sub></td><td>±10%</td><td>Uniform</td><td><b>K</b><sub>b2</sub></td></tr></tbody></table>

$m = 6$ and the random variables correspond to those given in Table 3). According to the PCE theory (Wiener, 1938; Xiu and Karniadakis, 2002, 2003), it can be approximated by a convergent in mean expansion. For numerical reasons, this expansion is truncated and writes:

$$Y(\xi) \approx \sum_{J=0}^{p-1} \alpha_J \Phi_J(\xi) \quad (17)$$

where $P$ is the number of terms in the expansion and corresponds to the size of the PCE basis. $\Phi_J$, $J \in [0, P-1]$ is the multivariate orthogonal polynomial basis and $(\alpha_J)$, $J \in [0, P-1]$ are the weighting coefficients to be determined.

Each random variable $(\xi_i)$, $i \in [1, m]$ is described by its PDF $p_i$ and the joined PDF of $\xi$ is $p = \prod_{i=1}^{m} p_i$. The multivariate orthogonal polynomial $\Phi_J$ are constructed by tensorisation of the $m$ mono-variate

polynomial families $(P_d^{(i)})$, $i \in [1, m]$ orthogonal to the PDF $p_i$, where $d$ denotes the degree of the polynomial. It writes:

$$\Phi_J(\xi) = \prod_{i=1}^{m} P_{J_i}^{(i)}(\xi_i) \quad (18)$$

where $J = \sum_{i=1}^{m} J_i$ is the degree of $\Phi_J$. The mono-variate orthogonal polynomials $(P_d^{(i)})$, $i \in [1, m]$ are chosen based on the Askey-scheme, which gives the correspondence between classical PDF and their associated orthogonal polynomials family. Here, only the uniform law is considered and so only sets of Legendre polynomials are used.

To truncate the PCE, one chooses the chaos order $\ell$ which corresponds to the maximal polynomial degree to keep in the expansion. It means that all polynomials that satisfy $J = \sum_{i=1}^{m} J_i \le \ell$ are kept in the expansion (17). The total number of terms $P$ in the expansion
is equal to $P = \binom{m+\ell}{\ell}$. One can clearly see here how the number of terms explodes when either the dimension $m$ increases or when the chaos order $\ell$ increases. Strategies have been developed to limit the number of terms in the expansion. In the present work the hyperbolic truncation scheme is adopted as it has proven to be efficient for rotordynamics (Denimal and Sinou, 2021b). This truncation scheme is based on the idea that models and physics are often driven by low order interaction effects. This norm is defined based on the scalar $\gamma \in [0, 1]$ fixed, and polynomials kept in Eq. (17) are those for which the degree satisfies (Blatman and Sudret, 2011):

$$ \|J\|_{\gamma} = \left( \sum_{i=1}^{m} J_{i}^{\gamma} \right)^{1/\gamma} \le \ell \qquad (19) $$

The coefficients $(\alpha_J)$, $J \in [0, P-1]$ are computed in a non-intrusive way here, based on the regression method (Sudret, 2008; Blatman and Sudret, 2011). To do so, a set of $M$ inputs $\Xi = \{\xi^{(1)}, \dots, \xi^{(M)}\}$ and their evaluations $\{Y(\xi^{(1)}), \dots, Y(\xi^{(M)})\}$ is generated. The $(\alpha_J)$, $J \in [0, P-1]$ are determined by minimising in the least-square sense the distance between the random function $Y(\xi)$ and the PCE evaluation at the $M$ points.

### 3.4. Hybrid formulation

Each function $f_i(\mathbf{x}, \xi)$ and $X_{j \times}^{f_i}(\mathbf{x}, \xi)$ is approximated by a surrogate model that couples kriging and PCE. First, the random part of each function is decomposed on a PCE approximation:

$$ f_i(\mathbf{x}, \xi) \approx \hat{f}_i(\mathbf{x}, \xi) = \sum_{J=0}^{P-1} \alpha_{J,f_i}(\mathbf{x}) \Phi_J(\xi) \qquad (20) $$

$$ X_{j \times}^{f_i}(\mathbf{x}, \xi) \approx \widehat{X_{j \times}^{f_i}(\mathbf{x}, \xi)} = \sum_{J=0}^{P-1} \alpha_{J,X_{j \times}^{f_i}}(\mathbf{x}) \Phi_J(\xi) \qquad (21) $$

where the chaos coefficients $(\alpha_{J,f_i})$, $J \in [0, P-1]$ and $(\alpha_{J,X_{j \times}^{f_i}})$, $J \in [0, P-1]$ depend on the parametric variables. As each PCE coefficient depends on the parametric variable $\mathbf{x}$, it can be approximated by kriging:

$$ \alpha_{J,f_i}(\mathbf{x}) = \mathbf{g}_{f_i}^{(J)}(\mathbf{x})^T \boldsymbol{\beta}_{f_i}^{(J)} + \mathbf{Z}_{f_i}^{(J)}(\mathbf{x}) \qquad (22) $$

and similarly for the $\alpha_{J,X_{j \times}^{f_i}}(\mathbf{x})$. By introducing these in Eqs. (20) and (21), it writes:

$$ f_i(\mathbf{x}, \xi) \approx \hat{f}_i(\mathbf{x}, \xi) = \sum_{J=0}^{P-1} \left( \mathbf{g}_{f_i}^{(J)}(\mathbf{x})^T \boldsymbol{\beta}_{f_i}^{(J)} + \mathbf{Z}_{f_i}^{(J)}(\mathbf{x}) \right) \Phi_J(\xi) \qquad (23) $$

and a similar expression is obtained for $\widehat{X_{j \times}^{f_i}(\mathbf{x}, \xi)}$. The construction of the hybrid surrogate model is based on a few evaluations of the expensive function (here $f_i(\mathbf{x}, \xi)$ or $X_{j \times}^{f_i}(\mathbf{x}, \xi)$). Two input sets are generated: one of $N$ values $\mathbf{x}^{(k)}$ based on a Latin Hypercube Sampling (LHS) with low discrepancy, and one of $M$ values of $\xi^{(k)}$ based on an LHS. The final input set is obtained by the tensorisation of the two sets and is composed of $N \times M$ points. For each of those points, the critical speeds $f_i$ and the amplitudes $X_{j \times}^{f_i}$ are computed. For each $\xi^{(k)}$, the $M$ evaluations of $f_i$ and $X_{j \times}^{f_i}$ for different $\xi$ are used for the PCEs construction. Thus, for each frequency and each amplitude, $N$ PCEs are constructed and so $N$ values of each $(\alpha_J)$, $J \in [0, P-1]$ are available and used for the kriging construction. The full workflow is given in Fig. 4. The method has been initially developed in Denimal et al. (2018), and extended in Denimal and Sinou (2021a,b). The interested reader could refer to these three references for more details.

### 3.5. PCEs coefficients exploitation

From the PCE coefficients, the stochastic moments of the considered output can be obtained directly. With the current modelling,

they directly depend on the parametric variable $\mathbf{x}$ and are directly obtained from the kriging evaluation without any expensive Monte Carlo Simulations (MCS). More precisely, the average writes:

$$ \mathbb{E}[f_i(\mathbf{x})] = \alpha_{0,f_i}(\mathbf{x}) = \mathbf{g}_{f_i}^{(0)}(\mathbf{x})^T \boldsymbol{\beta}_{f_i}^{(0)} + \mathbf{Z}_{f_i}^{(0)}(\mathbf{x}) \qquad (24) $$

and similarly for $\mathbb{E}[X_{j \times}^{f_i}(\mathbf{x})]$. The variance writes:

$$ \sigma_{f_i}^2(\mathbf{x}) = \sum_{J=0}^{P-1} \alpha_{J,f_i}(\mathbf{x})^2 \|\Phi_J(\xi)\|^2 = \sum_{J=0}^{P-1} \left( \mathbf{g}_{f_i}^{(J)}(\mathbf{x})^T \boldsymbol{\beta}_{f_i}^{(J)} + \mathbf{Z}_{f_i}^{(J)}(\mathbf{x}) \right)^2 \|\Phi_J(\xi)\|^2 \qquad (25) $$

and similarly for $\sigma_{X_{j \times}^{f_i}}^2(\mathbf{x})$.

### 3.6. Creation and validation of the surrogate models

This next subsection is devoted to the calibration of the different surrogate models and more specifically to the construction of the DoE. Because of the important computational time, doing a convergence study of the DoE for the hybrid model directly is not affordable. For this reason, the convergence of the PCE DoE and of the kriging DoE are done separately. The choice of the final size of the DoE for each one is based on these separate convergence studies. To ensure that the final hybrid meta-model with the tensorized DoE is accurate enough, a final validation by comparison with some exact reference values is done. The first part focuses on the parametric DoE (i.e. for kriging), and the second part on the random DoE (i.e. for the PCE).

#### 3.6.1. Convergence study of the kriging DoE

The size of the parametric DoE is chosen based on a convergence study on the deterministic case, i.e. only the location and the depth of the crack vary and all random parameters are set to their mean value. This is based on the idea that if the kriging is accurate enough for the deterministic case, it should also be accurate enough for each PCE coefficient. This is especially true for the first one as it follows the mean value, which should be close to the deterministic case.

For the convergence study, different DoE sizes are considered. It is worth noting here that the parameter $L_{crack}$ is discrete as it can only take value at the centre of the elements. All DoE are created as the union between two sets: one generated by LHS whose size can vary between 5 and 200 points, and a second set of constant points. For the latter, two strategies are considered and compared in the following. The first one consists in adding 16 points homogeneously distributed on the border of the design space: one at each angle and 3 equally spread on each segment. The second strategy corresponds to this first 16 points set completed by 20 points added in staggered rows on the lines $\mu = 0.9$ and $\mu = 1$. The addition of the first 16 points is based on the fact that kriging experiences convergence issues on the domain limits and by adding a few points, its performances are substantially improved (Forrester et al., 2008). The addition of the 20 other points comes from the fact that we know the crack has a strong influence on the nonlinear signature when it is deep, and so large variations of the amplitudes at sub-critical resonances and critical speeds are expected for high values of $\mu$ (Sinou and Lees, 2007). Of course, the main idea of this second strategy is to propose a better compromise between the number of points retained and the quality of the approximate solution. By densifying the DoE in this area, a higher efficiency of the kriging is expected. The two strategies are illustrated for one case with a 20 points LHS in Fig. 5 where red points correspond to the constant set and blue points to points generated by LHS.

From a practical point of view, the convergence of the kriging w.r.t the DoE size is studied. The LHS can vary between 5 and 200 points, and for each size 100 LHS are generated and merged with one of the two constant sets. To limit the computational time, reference simulations are done on a 20 x 20 grid, and DoE are projected on it so nonlinear simulations are done only once. For each physical output, a kriging surrogate model is computed using a constant trend and a
<figure><img src="image_5.png" /><figcaption>Fig. 4. General workflow of the method.</figcaption></figure>

Matérn 5/2 correlation kernel. The latter is used to predict the physical value on the 400 points of the 20 × 20 grid and the prediction error is computed using the relative Mean Squared Error (MSE):

$$e_{\text{rMSE}} = \frac{\sum_{n=1}^{N_{\text{val}}} (\hat{y}_n - y_n)^2}{\sum_{n=1}^{N_{\text{val}}} (y_n - \bar{y})^2} \quad (26)$$

where $\hat{y}_n$ is the kriging prediction at point n of the grid, $y_n$ is the reference value at point n and $\bar{y}$ is the average over the $N_{\text{val}} = 400$

points. This error is defined for each critical speed $i$, and denoted $e_{\text{rMSE}}(f_i)$, and for each amplitude $X_{j\times}^{f_i}$, and denoted $e_{\text{rMSE}}(X_{j\times}^{f_i})$. For each case, the average and the variance of this error over the 100 LHS for each DoE size are computed. The average operator is denoted $\mathbf{E}$ and the variance $\mathbf{Var}$. Their evolution is given in Fig. 6, where the first sampling strategy is in solid lines and the second in dashed lines, and where each colour corresponds to a mode. The shift between the two sampling strategies comes from the different number of points added to the LHS.

One can see that in all cases, when the DoE size increases then the relative MSE decreases in average and variance. More precisely,
<figure><img src="image_6.png" /><figcaption>Fig. 5. Illustration of the DoE construction for kriging where (•) are generated by LHS and (•) are fixed - (a) strategy 1: 16 points added - (b) strategy 2: 31 points added.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Evolution of the average (top figures) and variance (bottom figures) of the relative MSE error w.r.t the DoE size for the critical speeds (a) and the nonlinear amplitudes of the first order (b), second order (c) and third order (d). Mode 1 (—), Mode 2 (—), Mode 3 (●) and Mode 4 (■). First sampling strategy (—●—) and second sampling strategy (■□-■□-■□).</figcaption></figure>

the precision of the kriging (average values) and its sensitivity to the DoE (variance values) are improved of several orders of magnitude (up to 8 orders of magnitude for the average and up to 14 orders of magnitude for the variance). If for the critical speeds a convergence is

reached at about 160 points for all modes, no convergence is reached for the vibration amplitudes (i.e. a continuous decrease is observed with only a slight decrease in the slope of the curve) and so one must make a compromise between the numerical cost and the accuracy.
Comparing the two sampling strategies, it is clear that the second one (dashed lines) improves substantially the kriging performances. Indeed, in almost all cases the level of error is several order of magnitude lower for the second sampling strategy than for the first sampling strategy (see results for the first critical speed for all cases for example). As the error is lower for both the average and variance it means that the kriging performs better and is also less sensitive to the points distribution in the input space. This means that when the amount of information is limited (only a limited number of training points), it is better to use a priori knowledge of the problem to place the training points intelligently.

Based on these results, in the following the second sampling strategy is employed and the final DoE size is 92. The LHS is generated to maximise the minimum distance between points by optimising the ($\phi_p$) criterion (Pronzato and Müller, 2012). A specific attention has been put to the fact that $L_{crack}$ is discrete.

### 3.6.2. Convergence study of the PCE DoE

In a second phase, the size of the DoE for the PCE has been determined with a convergence study. To do so, an LHS of 150 random points has been generated. Then, for each of the 92 points of the parametric DoE, the 150 random points have been evaluated, which represents a total of 13,800 simulations. For each point of the parametric space (i.e. each of the 92 points) up to 120 points are used to create a PCE and the last 30 points are used for validation by comparison between reference values and predictions. The relative error $e_r$ for all the PCE is computed (i.e. based on $N = 30 \times 92 = 2760$ validation points) for different sizes of the DoE. For a given DoE size, 10 different LHS are generated so the average and the variance of the relative error are computed and displayed in Fig. 7. Results are presented here for a PCE of order 5 with a $\gamma$-norm of 0.5 (40 terms in the PCE). Different PCE properties have been compared and this one present the best compromise between the PCE size and the accuracy and so is kept for the rest of the study.

Looking at Fig. 7, one can see that the average and the variance of the relative error decreases when the DoE size increases for each critical speed $f_i$ and each amplitude $X_{i\times}^{f_i}$. The level of error for the critical speeds and the $1\times$ amplitudes are lower than for higher orders. This is explained by the fact that amplitudes of higher orders take extremely small values, which tends to increase the relative error. The convergence is reached after 100 points in the DoE, and so 100 points are taken in the following for the hybrid surrogate model creation.

### 3.6.3. Final validation

To ensure that the hybrid meta-model is accurate enough with the chosen tensorized DoE, a final validation is done by comparison with reference values. To do so, 200 input values are generated randomly, i.e. 200 vectors [$E, \rho, t_1, t_2, k_{x,2}, k_{y,2}, \mu, L_{crack}$]. On the one hand, the exact dynamic response is computed with the HBM, giving the output $y_n$ for the nth input point (y can be a critical frequency or an amplitude). On the other hand, the hybrid surrogate models are evaluated for the 200 points, giving the predicted outputs $\hat{y}_n$. They are finally compared, and the average absolute error over the 200 points is computed:

$$e_{final} = \frac{1}{N_{val}} \sum_{k=1}^{N_{val}} |\hat{y}_n - y_n| \quad (27)$$

the errors are given in Table 4. One can see they are all very low. As a visual validation, the comparisons between the references and the predictions are given in Fig. 8 where the predictions are in red stars and the reference in black circles. One can observe the predictions are accurate and that the outputs are correctly predicted and so the different hybrid surrogate models are trustful and can be used to simulate the rotordynamics for a negligible numerical cost.

### 3.6.4. Computational time of the method

A few words must be said about the computational efficiency of the approach. Simulation times are given for a basic laptop and for in-house

<table><caption>**Table 4**
Average absolute error of the hybrid surrogate model for the different critical speeds $f_i$ and different amplitudes $X_{i\times}^{f_i}$.</caption><thead><tr><th>Mode number</th><th>$f_i$</th><th>$X_{i\times}^{f_i}$</th><th>$X_{2\times}^{f_i}$</th><th>$X_{3\times}^{f_i}$</th></tr></thead><tbody><tr><td>1</td><td>1.41e-2</td><td>8.26e-06</td><td>4.26e-07</td><td>7.61e-08</td></tr><tr><td>2</td><td>1.45e-2</td><td>2.74e-07</td><td>1.29e-08</td><td>3.08e-10</td></tr><tr><td>3</td><td>6.18e-2</td><td>7.68e-06</td><td>1.76e-07</td><td>3.16e-08</td></tr><tr><td>4</td><td>7.24e-2</td><td>2.33e-06</td><td>9.39e-08</td><td>7.87e-09</td></tr></tbody></table>

algorithms that might not be the most efficient. As previously explained in Section 2.4, the evaluation of the full model requires about 30 s for one simulation performed via the HBM with a dichotomy process to identify the sixteen quantities of interest. If a MC strategy with 500 points for each crack characteristic were adopted, the obtaining of the mean and average would have required about $500 \times 50 \times 50 \times 30 = 37500000$ s $\approx 434$ days for a $50 \times 50$ grid. The construction of the hybrid surrogate models requires 9200 evaluations of the full model, i.e. 276000 s $\approx 3.2$ days. The construction of the hybrid surrogate models takes about 10 s. Once they are constructed:

* the evaluation of one hybrid surrogate model for one point takes about 7e-4 s;

* the construction of the map for the average and variance requires only the kriging prediction part and takes about 0.375s for one variable for a $50 \times 50$ grid;

* the construction of the 95% confidence interval requires an MCS at each crack characteristic, here 500 points, which represents about 3 min at the end for a $25 \times 25$ grid (12 min for a $50 \times 50$ grid).

One can see here the interest of the formulation when only stochastic information is needed as no MCS is required on the surrogate model (see the 0.375 s versus 12 min for same grid sizes). Finally, the interest of using such approach and the efficiency of the approach is obvious: the global simulation time is reduced from 434 days to about 3.2 days, and so is divided by 135. Comparing the simulation time of the full model to the prediction time of the hybrid surrogate model, it is divided by 4.2e4 (30 s versus 7e-4 s).

## 4. Results and discussions

The objective of this section is to discuss the possibility to perform a robust and reliable crack detection in a rotor system in an uncertain context based on the analysis of the linear and nonlinear signatures versus the parametric variation of both the crack depth $\mu$ and crack location $L_{crack}$ (with $\mu$ in $[0.05; 1]$ and $L_{crack}$ in $[0; L]$) and more specifically the evolution of the critical speeds $f_i$, the associated amplitudes $X_{i\times}^{f_i}$ and the maximum amplitudes $X_{i\times}^{f_i}$ at sub-critical resonances. To achieve this, the stochastic nonlinear response of the rotor is first analysed for different crack depth $\mu$ and crack location $L_{crack}$ (with $\mu$ in $[0.05; 1]$ and $L_{crack}$ in $[0; L]$). To get an affordable computational time, the hybrid surrogate models constructed previously are employed to simulate the unbalanced response of the rotor, where parametric parameters are the crack depth and location (influence through kriging) and random parameters are the model uncertainties (influence through PCE). The mean and variance of the critical speeds $f_i$ and of the amplitudes $X_{i\times}^{f_i}$ over the geometrical imperfections, random material properties and uncertain boundary conditions are determined for a parametric variation of the crack depth and location. From this, the impact of the crack, in terms of depth and location, on the stochastic response of the rotor are deeply analysed. Finally, the evolution of the 95% interval of these quantities (critical speeds $f_i$ and amplitudes $X_{i\times}^{f_i}$) are computed in order to see if these critical speeds or these amplitudes could be a reliable and robust indicator for crack detection when many uncertainties are present.
<figure><img src="image_8.png" /><figcaption>Fig. 7. Evolution of the average (a) and of the variance (b) of the average relative error for the PCE for different learning set sizes.</figcaption></figure>

## 4.1. Evolution of the critical speeds and vibration amplitudes

The average and variance of the critical speeds $f_i$ and amplitudes $X_{j×}^{f_i}$ are directly computed from Eq. (24) and Eq. (25), respectively. The mean values over the parametric space are given in Fig. 9 and the variance values in Fig. 10. Looking at Fig. 9, the trends of the evolution of the mean value for each quantity $f_i$ and $X_{j×}^{f_i}$ as a function of the size and position of the crack are of course in agreement with the evolution of the previous results presented in the deterministic case. For example, an increase in the size of the crack for a given crack position, increases the variation of the quantity considered compared to the healthy rotor. Moreover, the most important variations of the critical speeds $f_i$ and vibration amplitudes $X_{j×}^{f_i}$ are observed when the crack is located at an anti-node of the considered mode. For the vibration responses of orders 2 and 3 (i.e. $X_{2×}^{f_i}$ and $X_{3×}^{f_i}$), one can see that they are more sensitive to the existence of a crack than the critical speeds $f_i$ and the amplitude of first order $X_{1×}^{f_i}$ as they experience larger variation in the presence of a crack.

Considering the evolution of the variances given in Fig. 10, their evolution is more complex than the averages and depends on each critical speed and amplitude. Globally, variances are relatively large for the critical speeds $f_i$ and represent a non-negligible influence of the random parameters on the critical speeds, much larger than the variations observed on average. For $f_1$ and $f_2$, the variance remains almost constant over the parametric space with a small evolution that follows approximately the same evolution as the mean (i.e. a small increase of the variance when the crack size increases and the maximum of the variance is observed for a crack in the middle of the rotor). For $f_3$ and $f_4$, we can distinguish a non negligible evolution but which is not continuous depending on the size or position of the crack. Thus, these results indicate that the modelling uncertainties have a non-negligible impact on the critical speeds.

To go into details and illustrate the complexity brought by the addition of modelling uncertainties, one can observe that the variance evolution of the vibration amplitudes $X_{j×}^{f_i}$ is completely different from one case to another and that the impact of the modelling uncertainties depends on the variable under consideration. For the vibration amplitudes $X_{1×}^{f_2}$, $X_{1×}^{f_3}$ and $X_{1×}^{f_4}$, variances are relatively large in comparison

to the variations observed on average, with an almost constant value over the parametric space for $X_{1×}^{f_2}$ and variations for $X_{1×}^{f_3}$ and $X_{1×}^{f_4}$. More specifically, for $X_{1×}^{f_3}$ and for $X_{1×}^{f_4}$, the largest variations of the variance when the crack depth increases are observed when the crack is located at the anti-node of the modeshape, and the minimum of variation of the variance are observed when the crack is located at a node of the modeshape (for e.g. see the line $L_{crack} = 0.25$ m in Fig. 10(j) and the constant value of the variance). For the second backward mode, when the crack is located at an anti-node, the variance of $X_{1×}^{f_3}$ tends to decrease when the crack becomes deeper, whereas for the forward mode the variance of $X_{1×}^{f_4}$ tends to increase when the crack becomes deeper. So depending on the mode, the system is more or less sensitive to the modelling uncertainties in the presence of a crack. For the other cases, the evolution of the variance are rather low compared to the values of the associated mean. So for these cases, the modelling uncertainties have a lower impact on the cracked rotor. The variances experience strong variations with maxima reached in the vicinity of the anti-node of the corresponding modeshape. For these cases, the variance is also clearly affected by the increase in the crack size and increases when the crack is deeper and located around an anti-node.

In conclusion, the presence of uncertainties in the system leads to very different fluctuations of the mean and variance depending on the quantities considered. This can make quite difficult at first sight to find the quantities $f_i$ or $X_{j×}^{f_i}$ that represent robust indicators for fault detection in a system with uncertainties. In fact it is of course sufficient to combine the mean and variance and their evolution for a given quantity $f_i$ or $X_{j×}^{f_i}$ in order to be able to examine the possibility of detecting a crack in the presence of uncertainties for this selected quantity. In the following paragraph we propose to develop this concept through the construction of a Confidence Interval (CI).

## 4.2. Discussion on the possibility of crack detection in the presence of uncertainties

Here, we focus on monitoring the nonlinear vibrations of a rotating crack system and the detection of its associated crack by carrying out a detailed analysis of the nonlinear dynamics involved in the presence of uncertainties. The objective of this part is to show the characteristic
<figure><img src="image_9.png" /><figcaption>Fig. 8. Comparison between the reference points (black) and the hybrid surrogate models predictions (red) for the different critical speeds $f_i$ and amplitudes $X_{j×}^{f_i}$.</figcaption></figure>

differences that make it possible to determine the presence of a crack in a rotor in the presence of uncertainties, as well as the limitations and precautions to be taken to avoid an erroneous diagnostic.

In order to achieve such an objective, we propose to build the commonly used 95% Confidence Interval (CI) of the outputs. A 95% confidence interval is a range of values for which there is 95% certainty that it contains the true population mean. For each $x = [\mu, L_{crack}]$, the 95% confidence interval is deduced from the simulation of an LHS of 1000 points. The 95% confidence intervals of the physical output $a$ is denoted $CI_a(\mu, L_{crack})$. The following interpretations of the confidence interval can be given:

* if $CI_a(\mu, L_{crack})$ has little variation in terms of the two parameters $\mu$ and $L_{crack}$, then the quantity $a$ is not a robust indicator for the detection of a crack. Indeed, in this case, the presence of uncertainties generates a confidence interval that includes the proper variations of the quantity $a$ according to the two parameters $\mu$ and $L_{crack}$. In other words, the variation of the quantity $a$ due to uncertainties is greater than the variation of the quantity $a$ depending on the two parameters related to the presence of the crack (i.e. $\mu$ and $L_{crack}$), which makes the parameter $a$ ineffective as a robust indicator of the crack detection in the presence of uncertainties.

* if the variation of $CI_a(\mu, L_{crack})$ versus the two parameters $\mu$ and $L_{crack}$ is visible, then $CI_a$ evolves by a sufficiently large amount, even in the presence of uncertainties, which makes it possible to corroborate the evolution of the parameter $a$ with the presence of a crack.

Fig. 11 gives the 95% confidence interval as a function of the position and size of the crack for the critical speeds $f_i$ (with $i = 1, ..., 4$) and the amplitudes $X_{j×}^{f_i}$ (with $i = 1, ..., 4$ and $j = 1, ..., 3$). The red surface represents the lower bound and the blue surface the upper bound. The black lines correspond to the projection of the minimum and maximum values. First of all, it is observed that the quantities $C_{f_i}$ are almost unchanged in terms of the two parameters $\mu$ and $L_{crack}$. The evolution of frequencies $f_i$ (for $\mu = [0.05; 1]$ and $L_{crack} = [0; 0.5]$ m) previously shown in Fig. 3 is totally included in the associated confidence interval $C_{f_i}$. This has a substantial impact as it means that the decrease in frequencies with increasing $\mu$ is too low to be used as a robust criterion during the detection of a crack in a rotor system in the presence of uncertainties.

Then regarding the validity of the amplitudes $X_{j×}^{f_i}$ for robust crack detection, different conclusions can be drawn which can be classified in three categories:
<figure><img src="image_10.png" /><figcaption>Fig. 9. Evolution of the average value over the parametric space for the different critical speeds $f_i$ and amplitudes $X_{jX}^{f_i}$.</figcaption></figure>

• Category 1: some quantities lead to the same conclusions as for the evolution of the frequencies $f_i$ with uncertainties (see more specifically the evolutions of $C_{X_{1X}f_2}$, $C_{X_{1X}f_3}$ and $C_{X_{1X}f_4}$). The variations of the minimum and maximum bounds of confidence interval are too small compared to the width of the latter. This means that the quantities $X_{1X}^{f_2}$, $X_{1X}^{f_3}$ and $X_{1X}^{f_4}$ cannot be used for the crack detection in the presence of uncertainties.

• Category 2: some variables such as $X_{1X}^{f_1}$, $X_{2X}^{f_3}$, $X_{2X}^{f_4}$, $X_{3X}^{f_3}$ and $X_{3X}^{f_4}$ can be used to detect a crack in an approximate way (i.e. these indicators will only detect a crack without leading to a possibility to identify the size of the crack for example). Indeed, the evolutions of these quantities $X_{jX}^{f_i}$ versus the two crack parameters are predominant in comparison with the width of the confidence interval for non-negligible variations of $\mu$ and $L_{crack}$. However, for a small variation of one of these two parameters, the presence of uncertainties implies that the evolution of the signature $X_{jX}^{f_i}$ remains within the previous confidence interval.

• Category 3: some variables appear to be reliable indicators for ensuring the detection of a crack (even of small size). The effect of a crack on the monitoring of these super harmonic orders at the passages of the associated sub-critical resonances can be seen very clearly even in the presence of uncertainties. See for example the evolutions of $C_{X_{2X}f_1}$, $C_{X_{3X}f_1}$, $C_{X_{2X}f_2}$ and $C_{X_{3X}f_2}$ as a function of $\mu$ and

$L_{crack}$ in Figs. 11. In addition two of them (i.e. $X_{2X}^{f_2}$ and $X_{3X}^{f_2}$) can be used in the presence of uncertainties even for a crack of small size.

If we consider more specifically the quantities $X_{jX}^{f_i}$ falling into the last two categories, notwithstanding the addition of uncertainties, we also find some of the classical results of the dynamics of cracked systems:

• increasing the crack size $\mu$ induces an increase of the amplitudes which can be different from one quantity to another. For example, a rapid evolution of the amplitudes $X_{1X}^{f_1}$ and $X_{2X}^{f_1}$ is observed when the crack depth is close to $\mu = 1$ while a more progressive and continuous evolution as a function of the crack size is observed in the other cases (see for example $X_{2X}^{f_2}$, $X_{3X}^{f_2}$ and $X_{2X}^{f_4}$).

• the closer the crack is to a vibration node of the mode $i$th, the more the value of the associated amplitude $X_{jX}^{f_i}$ tends towards zero. In other words, we find here the fact that when the crack is located close to a node (for a given mode), the associated super harmonic orders are affected only slightly by the presence of the crack.

• the effect of a crack on the amplitude $X_{jX}^{f_i}$ is maximal when the crack is located on one of the anti-nodes of the mode $i$th considered.
<figure><img src="image_11.png" /><figcaption>Fig. 10. Evolution of the variance over the parametric space for the different critical speeds $f_i$ and amplitudes $X_{j×}^{f_i}$.</figcaption></figure>

So the appearance of new peaks $X_{j×}^{f_i}$ (with $j \ge 2$) at passages through the sub-critical resonances $\frac{1}{j}$ of a cracked rotor can be used as a robust indicator of the presence of a crack, even if uncertainties or lack of knowledge of the physical parameters of the rotor are accepted. It is also interesting to note that even if some peaks $X_{j×}^{f_i}$ are not predominant in relation to the complete nonlinear dynamic response of the rotor system, their evolutions with increasing the crack depth can be non negligible which constitutes an interesting and robust indicator for detecting a crack on a rotor in the presence of uncertainties. However, caution is required regarding these observations, since the intensity of each of the super harmonic orders $j \times$ (with $j \ge 2$) depends on the nonlinear efforts of the crack which correspond to a parametric nonlinearity directly linked to the unbalances of the rotor, gravity and thus implicitly to the dephasing of the latter regarding the front of the breathing crack.

The global objective was to investigate the possibility to get reliable and robust indicators for crack detection in rotors with numerous uncertainties due to lack of knowledge or environmental variations. The previous discussions have demonstrated that the tracking of the critical speeds $f_i$ or the amplitudes $X_{1×}^{f_i}$ for the health monitoring of rotor systems with uncertainties is not an appropriate and efficient

strategy. However, the nonlinear signature could be used as an indicator robust to these uncertainties. More specifically, the increase of vibration amplitudes $X_{j×}^{f_i}$ (with $j \ge 2$) when passing through sub-critical resonances $\frac{1}{j}$ can be used as robust indicators for detecting a transverse crack in a rotating system and its evolution with time (due to the increase of the crack depth with time for example). The monitoring of super-harmonic orders $j \times$ (with $j \ge 2$) therefore appears the most reliable criterion for ensuring the rapid detection of the occurrence of a crack (even a small one). However, caution is necessary regarding the amplitudes $X_{j×}^{f_i}$ to be selected as robust indicators in order to achieve this.

## 5. Conclusion

The research and results presented in this paper show the essential contribution of taking into account the nonlinear dynamic behaviours of rotors to diagnose the presence of a transverse breathing crack in rotating systems and permit better vibration monitoring of rotors in the presence of uncertainties.

The monitoring of orders $j \times$ with $j \ge 2$ and more particularly the evolutions of the associated amplitudes when passing through
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Envelop of the 95% confidence interval (blue surface for upper bound and red surface for lower bound) over the parametric space for the different critical speeds $f_i$ and amplitudes $X_{j×}^{f_i}$ and projection of the minima and maxima values (black lines).</figcaption></figure>

sub-critical resonances provides indicators that can be used to detect the occurrence of cracks in a rotating system in the presence of uncertainties. The evolution of the resonance frequencies does not provide interesting information for the detection of a breathing crack because these evolutions are generally hidden by the variations of the frequencies according to the uncertainties. Nonetheless, it appears that the detection of a crack of small size and a precise estimation of the characteristics associated with it (localisation and size) can be more difficult to perform due to the presence of uncertainties in the system. In the latter case, a fine selection of the amplitude indicators $X_{j×}^{f_i}$ to be monitored is necessary. These results taken together also show the considerable difficulty and complexity of interpreting the physical phenomena observed and the identification of crack parameters (depth and position) due to the influence of a large number of physical parameters that can lead to variations of the different amplitudes $j×$. Also, although it appears clearly possible to detect the presence of a crack for a rotating system comprising uncertainties, the phase of identifying the position and depth of a crack remains an open and difficult problem when considering the presence of uncertainties or when gaps in knowledge of the healthy system remain.

All these results and the resulting conclusions were made possible by the use of an advanced computational technique based on kriging and Polynomial Chaos Expansion for the prediction of the critical speeds and the nonlinear signatures of the vibration responses of the

cracked rotor with uncertainties. Indeed, due to the high cost of computational models with time-dependent coefficients and the resolution based on the HBM for computing periodic solutions of ordinary differential equations, classical techniques such as Monte Carlo simulation are not applicable.

Some non-exhaustive interesting outlooks for future research are as follows:

* the proposed study uses a simplified breathing behaviour for the crack based on the research work of Davies and Mayes (1984) and Mayes and Davies (1984). It should be interesting to apply the proposed advanced hybrid approach for more complex efficient breathing crack model and possibly crack propagation techniques. A potential drawback of the proposed study is that if the basic model is assumed to be approximate, the associated results could also tend to differ from what might have been accurate.

* it would be interesting to apply the proposed hybrid approach on a real industrial structure and to compare the results with full experimental tests. One of the most important challenges may be to identify the different sources of uncertainty but also to be able to characterise and quantify the uncertainty related to each parameter for a complex system. It means being able to determine the type of PDF that follows each random parameter as well as the range of variation in a real environment and in real conditions.
* considering the methodology, it would be interesting to consider more efficient sampling strategies to reduce the size of the training set, which is currently large. Indeed, its construction is based on a tensorisation, which makes explode the size of the final training set. One challenge would be to find a more efficient way to train this hybrid surrogate model and avoids this expensive tensorisation.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgement

J.-J. Sinou acknowledges the support of the Institut Universitaire de France.

## References

Adewusi, S., Al-Bedoor, B., 2001. Wavelet analysis of vibration signals of an overhang rotor with a propagating transverse crack. J. Sound Vib. 246 (5), 777-793.

Adewusi, S., Al-Bedoor, B., 2002. Experimental study on the vibration of an overhanged rotor with a propagating transverse crack. Shock Vib. 9, 91-104.

Al-Shudeifat, M.-A., Butcher, E.-A., 2011. New breathing functions for the transverse breathing crack of the cracked rotor system: Approach for critical and subcritical harmonic analysis. J. Sound Vib. 330 (3), 526-544.

Al-Shudeifat, M.-A., Butcher, E.-A., Stern, C.-R., 2010. General harmonic balance solution of a cracked rotor-bearing-disk system for harmonic and sub-harmonic analysis: Analytical and experimental approach. Internat. J. Engrg. Sci. 48 (10), 921-935, Structural Health Monitoring in the Light of Inverse Problems of Mechanics.

Bachschmid, N., Pennacchi, P., Tanzi, E., 2008. Some remarks on breathing mechanism, on non-linear effects and on slant and helicoidal cracks. Mech. Syst. Signal Process. 22, 879-904.

Bachschmid, N., Pennacchi, P., Tanzi, E., 2010. Cracked rotors: A survey on static and dynamic behaviour including modelling and diagnosis. In: Springer Science and Business Media. Berlin, Germany, pp. 129-148.

Bachschmid, N., Pennacchi, P., Tanzi, E., Vania, A., 2000. Identification of transverse crack position and depth in rotor systems. Meccanica 35, 563-582.

Blatman, G., Sudret, B., 2011. Adaptive sparse polynomial chaos expansion based on least angle regression. J. Comput. Phys. 230 (6), 2345-2367.

Bovsunovsky, A., Surace, C., 2015. Non-linearities in the vibrations of elastic structures with a closing crack: A state of the art review. Mech. Syst. Signal Process. 62-63, 129-148.

Collins, K., Plaut, R., Wauer, J., 1991. Detection of cracks in rotating Timoshenko shafts using axial impulses. J. Vib. Acoust. 113, 74-78.

Darpe, A., 2007. Coupled vibrations of a rotor with slant crack. J. Sound Vib. 305, 172-193.

Darpe, A., Gupta, K., Chawla, A., 2002. Analysis of the response of a cracked jeffcott rotor to axial excitation. J. Sound Vib. 249, 429-445.

Darpe, A., Gupta, K., Chawla, A., 2004a. Coupled bending, longitudinal and torsional vibrations of a cracked rotor. J. Sound Vib. 269, 33-60.

Darpe, A., Gupta, K., Chawla, A., 2004b. Transient response and breathing behaviour of a cracked jeffcott rotor. J. Sound Vib. 272, 207-243.

Davies, W., Mayes, I., 1984. The vibrational behaviour of a multi-shaft, multi-bearing system in the presence of a propagating transverse crack. Transactions of the ASME Journal of Vibration, Acoustics, Stress, and Reliability in Design 106, 146-153.

Denimal, E., Nechak, L., Sinou, J.-J., Nacivet, S., 2018. A novel hybrid surrogate model and its application on a mechanical system subjected to friction-induced vibration. J. Sound Vib. 434, 456-474.

Denimal, E., Sinou, J.-J., 2021a. Advanced computational technique based on kriging and polynomial chaos expansion for structural stability of mechanical systems with uncertainties. J. Eng. Math. 130 (1), 1-19.

Denimal, E., Sinou, J.-J., 2021b. Advanced kriging-based surrogate modelling and sensitivity analysis for rotordynamics with uncertainties. Eur. J. Mech. A Solids 104331.

Dimarogonas, A., 1996. Vibration of cracked structures: A state of the art review. Eng. Fract. Mech. 55, 831-857.

Fan, W., Qiao, P., 2011. Vibration-based damage identification methods: A review and comparative study. Struct. Health Monit. 10, 83-111.

Forrester, A., Sobester, A., Keane, A., 2008. Engineering Design Via Surrogate Modelling: A Practical Guide. John Wiley & Sons.

Fu, C., Ren, X., Yang, Y., Lu, K., Qin, W., 2019. Steady-state response analysis of cracked rotors with uncertain-but-bounded parameters using a polynomial surrogate method. Commun. Nonlinear Sci. Numer. Simul. 68, 240-256.

Fu, C., Ren, X., Yang, Y., Lu, K., Wang, Y., 2018. Nonlinear response analysis of a rotor system with a transverse breathing crack under interval uncertainties. Int. J. Non-Linear Mech. 105, 77-87.

Fu, C., Xu, Y., Yang, Y., Lu, K., Gu, F., Ball, A., 2020. Dynamics analysis of a hollow-shaft rotor system with an open crack under model uncertainties. Commun. Nonlinear Sci. Numer. Simul. 83, 105102.

Ganguly, K., Roy, H., 2021. A novel geometric model of breathing crack and its influence on rotor dynamics. J. Vib. Control 1-10.

Gao, H., Guo, X., Hu, X., 2012. Crack identification based on kriging surrogate model. Struct. Eng. Mech. 41 (1), 191-206.

Gao, H., Guo, X., Ouyang, H., Han, F., 2013. Crack identification of cantilever plates based on a kriging surrogate model. J. Vib. Acoust. 135, 1-10.

Gasch, R., 1976. Dynamic behaviour of a simple rotor with a cross sectional crack. In: International Conference on Vibrations in Rotating Machinery, Paper C178/76, IMechE, Vol. 1, pp. 123-128.

Gasch, R., 1993. A survey of the dynamic behaviour of a simple rotating shaft with a transverse crack. J. Sound Vib. 160 (2), 313-332.

Giannopoulos, G., Georgantzinos, S., Anifantis, N., 2015. Coupled vibration response of a shaft with a breathing crack. J. Sound Vib. 336, 191-206.

Guo, C., Yan, J., Yang, W., 2017. Crack detection for a jeffcott rotor with a transverse crack: An experimental investigation. Mech. Syst. Signal Process. 83, 260-271.

Henry, T., Okah, B., 1976. Vibration in cracked shafts. In: IMechE Conference on Vibrations in Rotating Machinery, C162/76, Vol. 1, pp. 15-19.

Ishida, Y., Inoue, T., 2006. Detection of a rotor crack using a harmonic excitation and nonlinear vibration analysis. J. Vib. Acoust. 128, 741-749.

Iwasubo, T., Arii, S., Oks, A., 1992. Detection of a transverse crack in rotor shaft by adding external force. In: International Conference on Vibrations in Rotating Machinery, Paper C432/093, IMechE, Vol. 1, pp. 275-282.

Jones, D., Schonlau, M., Welch, W., 1998. Efficient global optimization of expensive black-box functions. J. Global Optim. 13, 455-492.

Jun, O.-S., 2012. Dynamic behavior analysis of cracked rotor based on harmonic motion. Mech. Syst. Signal Process. 30, 186-203.

Kumar, C., Rastogi, V., 2009. A brief review on dynamics of a cracked rotor. Int. J. Rotating Mach. 1-6.

Liao, M., Gasch, R., 1992. Crack detection in rotating shafts - An experimental study. In: International Conference on Vibrations in Rotating Machinery, Paper C432/106, IMechE, Vol. 1, pp. 289-295.

Lin, Y., Chu, F., 2009. Numerical and experimental investigations of flexural vibrations of a rotor system with transverse or slant crack. J. Sound Vib. 324, 107-125.

Lophaven, S., Nielsen, H., Søndergaard, J., 2002. DACE: A Matlab Kriging Toolbox, Vol. 2. Citeseer.

Lu, Z., Lv, Y., Ouyang, H., 2019. A super-harmonic feature based updating method for crack identification in rotors using a kriging surrogate model. Appl. Sci. (ISSN: 2076-3417) 9 (12).

Mayes, I., Davies, W., 1976. The vibrational behaviour of a rotating system containing a transverse crack. In: IMechE Conference on Vibrations in Rotating Machinery, C/168/76, Vol. 1, pp. 53-64.

Mayes, I., Davies, W., 1984. Analysis of the response of a multi-rotor-bearing system containing a transverse crack in a rotor. Trans. ASME J. Vib. Acoust. Stress Reliab. Des. 106, 139-145.

Ostachowicz, W., Krawczuk, M., 1992. Coupled torsional and bending vibrations of a rotor with an open crack. Arch. Appl. Mech. 62, 191-201.

Papadopoulos, C., Dimarogonas, A., 1987a. Coupled longitudinal and bending vibrations of a rotating shaft with an open crack. J. Sound Vib. 117, 81-93.

Papadopoulos, C., Dimarogonas, A., 1987b. Coupling of bending and torsional vibration of a cracked Timoshenko shaft. Arch. Appl. Mech. 57 (4), 257-266.

Papadopoulos, C., Dimarogonas, A., 1992. Coupled vibrations of cracked shafts. J. Vib. Acoust. 114, 461-467.

Pronzato, L., Müller, W.G., 2012. Design of computer experiments: space filling and beyond. Stat. Comput. 22 (3), 681-701.

Rasmussen, C., 2003. Gaussian processes in machine learning. In: Summer School on Machine Learning. Springer, pp. 63-71.

Sawicki, J., Friswell, M.-I., Kulesza, Z., Wroblewski, A., Lekki, J.-D., 2011. Detecting cracked rotors using auxiliary harmonic excitation. J. Sound Vib. 330 (7), 1365-1381.

Schmied, J., Kramer, E., 1984. Vibrational behaviour of a rotor with a cross-sectional crack. In: IMechE Conference on Vibrations in Rotating Machinery, C279/84, Vol. 1, pp. 15-19.

Sekhar, A.S., 2003a. Crack detection and monitoring in a rotor supported on flube film bearing: start-up vs run-down. Mech. Syst. Signal Process. 17 (4), 897-901.

Sekhar, A., 2003b. Crack detection through wavelet transform for a run-up rotor. J. Sound Vib. 259, 461-472.

Sekhar, A., Prabhu, B., 1998. Condition monitoring of cracked rotors through transient response. Mech. Mach. Theory 33 (8), 1167-1175.

Sinou, J.-J., 2008a. Detection of cracks in rotor based on the 2X and 3X super-harmonic frequency components and the crack-unbalance interactions. Commun. Nonlinear Sci. Numer. Simul. 13 (9), 2024-2040.

Sinou, J.-J., 2008b. An experimental investigation of condition monitoring for notched rotors through transient signals and wavelet transform. Exp. Mech. 49, 683-695.
Sinou, J.-J., 2009a. Experimental response and vibrational characteristics of a slotted rotor. Commun. Nonlinear Sci. Numer. Simul. 14, 3179–3194.

Sinou, J.-J., 2009b. A review of damage detection and health monitoring of mechanical systems from changes in the measurement of linear and non-linear vibrations. In: Sapri, R.R. (Ed.), Mechanical Vibrations: Measurement, Effects and Control, Chapter 13. Nova Science Publishers, ISBN: 978-1-60692-037-4, pp. 643–702.

Sinou, J.-J., Faverjon, B., 2012. The vibration signature of chordal cracks in a rotor system including uncertainties. J. Sound Vib. 331 (1), 138–154.

Sinou, J.-J., Lees, A., 2005. Influence of cracks in rotating shafts. J. Sound Vib. 285 (4-5), 1015–1037.

Sinou, J.-J., Lees, A., 2007. A non-linear study of a cracked rotor. J. Eur. Mech. - A/Solids 26 (1), 152–170.

Spagnol, J., Wu, H., Yang, C., 2020. Application of non-symmetric bending principles on modelling fatigue crack behaviour and vibration of a cracked rotor. Appl. Sci. 10 (2), 717.

Sudret, B., 2008. Global sensitivity analysis using polynomial chaos expansions. Reliab. Eng. Syst. Saf. 93 (7), 964–979.

Wang, D., Hua, C., Dong, D., He, B., Lu, Z., 2018. Crack parameters identification based on a kriging surrogate model for operating rotors. Shock Vib. ID 9274526, 1–12.

Wauer, J., 1990. Dynamics of cracked rotors: Literature survey. Appl. Mech. Rev. 43, 13–17.

Wiener, N., 1938. The homogeneous chaos. Amer. J. Math. 60 (4), 897–936.

Xiu, D., Karniadakis, G., 2002. The Wiener–Askey polynomial chaos for stochastic differential equations. SIAM J. Sci. Comput. 24 (2), 619–644.

Xiu, D., Karniadakis, G., 2003. Modeling uncertainty in flow simulations via generalized polynomial chaos. J. Comput. Phys. 187 (1), 137–167.

Yongfeng, Y., Qinyu, W., Yanlin, W., Weiyang, Q., Kuan, L., 2019. Dynamic characteristics of cracked uncertain hollow-shaft. Mech. Syst. Signal Process. 124, 36–48.

Zou, J., Chen, J., 2004. A comparative study on time-frequency feature of cracked rotor by wigner-Ville distribution and wavelet transform. J. Sound Vib. 276, 1–11.