# Kriging Metamodeling in Rotordynamics: Application for Predicting Critical Speeds and Vibrations of a Flexible Rotor

Authors: J.-J. Sinou; L. Nechak; S. Besset
Journal: Complexity
 (2018)

---

## Abstract

Rotating machinery produces vibrations depending upon the design of the rotor systems as well as any faults or uncertainties in the machine that can increase the vibrations of such systems. This study illustrates the effectiveness of using surrogate modeling based on kriging in order to predict the vibrational behavior (i.e., the critical speeds and the vibration amplitudes) of a complex flexible rotor in the presence of uncertainties. The basic idea of kriging is to predict unknown values of a function by using a small size set of known data. The kriging estimation is based on a weighted average of the known values of the function in the neighborhood of the point for which the value of the function has to be calculated. The crucial dependence of a kriging predictor versus the correlation functions and different orders will be illustrated. This paper also shows that reducing the number of samples required to have predictive models can be achieved by performing an initial understanding of the mechanical system of interest and by considering certain characteristics directly deriving from the physics of the problem studied.


---

## Body

## 1. Introduction

Rotating machinery, which generally consists of a mechanical assembly that has rotating structures (i.e., rotors) and one or more supporting structures (i.e., stators), is an important element in many structures of engineering applications. Rotor dynamics is a specialized branch of applied mechanics concerned with the behavior and diagnosis of rotating structures. It is commonly used to analyze the behavior of structures ranging from jet engines and steam turbines to autoengines and computer disk storage. If this aspect is ignored, it might result in loss of the equipment, excessive wear and tear on the machinery, catastrophic breakage beyond repair, or even human injury and loss of lives. In order to avoid catastrophic breakage beyond repair of rotating machinery or even human injury, numerical simulations and advanced strategies must be developed to predict all the possible dynamic behaviors in a predefined design of interest including potential uncertainties on physical parameters. The objective here should be to avoid operations that are close to the critical and pass safely through them when in acceleration or deceleration. In order to improve the safety and efficiency of such rotating systems, a specific and sometimes complex design process has to be performed [1–3].Vibration behavior of rotating structures due to imbalance is one of the main aspects of rotating machinery which must be studied in detail and considered while designing. For instance, as the speed of rotation increases, vibrations may appear and the amplitude of vibration often passes through a maximum that is called a critical speed. This may result in machine failure due to the generation of excessive vibrations. Other aspects such as instability or other faults [4] (such as misalignment, bow, and asymmetric shaft) can be important issues to consider as well in order to avoid worst design. In the case of rotating systems composed of several rotating discs, the observed dynamic behavior is not always easily intuitive. For example, forward and backward critical speeds decrease and increase, respectively, with the gyroscopic stiffening effect of each rotating disc. So it generally requires many calculations to achieve a good design. Moreover the real dynamics of rotating machinery is sometimes difficult to model theoretically. Indeed, during a phase of preliminary design, many physical parameters can have a significant influence on the vibration response. For example, some parameters may vary throughout the life of a rotating machinery or some uncertainties have to be taken into account during the manufacturing monitoring. All these facts can drastically affect the vibration responses of a rotor system (see, e.g., [5–13]). And therefore it complicates the approach to be put in place in order to propose a robust design. Also, due to the complexity of the systems studied and the presence of uncertainties, it is generally quite difficult to deal with all possible cases in the design space. So the chief concern of engineers who design large rotors is to be able to predict the dynamic behavior of rotating machinery with a low computational cost (i.e., a reduced number of calculations). Thus the crucial objective is to establish the best dynamic performance in a chosen design space and in the presence of uncertainties. This study is then dedicated to the prediction of dynamic properties of a rotor system under parameter uncertainty. In fact the related variability characterized by first- and second-order moments is required for the defining of suitable strategies for robust design of rotor systems. For instance, strategies based on Tagushi rules use variance-to-mean ratio whereas others use a mathematical programming approach; for example, this approach minimizes the mean while keeping the standard deviation below a prespecified threshold [14].In order to considerably reduce the number of calculations, various strategies based on surrogate models have been developed. More specifically kriging formalism and advanced and recent concepts [14–20] have been used successfully for different branches of applied mechanics. The basic concept is to replace any complex model by a suitable surrogate model which offers a compromise between the accuracy of its predictions and the cost related to its implementation. So this surrogate model approximates a parameter-dependent function by using a small number of parameter samples generated according to a predefined experimental design. Of course the accuracy of the surrogate depends on the number and location of samples in the design space of interest and these surrogate models are characterized by some tuning parameters (such as the order of the regression model and the spatial correlation function) that control their accuracy.The main objective of the present study is to illustrate the efficiency of using kriging surrogate model to accurately predict the dynamic behavior (i.e., the critical speeds as well as the maximum vibrational amplitudes) of a flexible rotor with multiple uncertainties. A second goal is to show the usefulness of considering certain physical properties of the problem in order to reduce the number of samples required to achieve relevant kriging models. The paper is organized as follows: first the finite-element model of the flexible rotor under study is presented. Secondly kriging based modeling is described and discussed. Finally the method is applied to the rotor example. Numerical results in the case of multiple uncertainties are presented and commented on. The issue of computation time and major drawbacks of the proposed strategy are also addressed.


## 2. Rotor System

The flexible rotor under study is shown in Figure1(a). This rotor system is composed of a shaft with four discs. It is supported at the ends by bearings.Figure 1
Finite-element model of the rotor system: (a) rotor system and (b) coordinates in fixed frame.
(a) (b)


### 2.1. Shaft Elements

The shaft is discretized into 10 Euler beam finite elements with circular cross-sections. Each beam finite element has four degrees of freedom at each node (two lateral displacements and two rotations, the axial and torsional degrees of freedom being not considered). The nodal displacement vectors for theith shaft element defined as(1)xs,i=v1w1θ1ψ1v2w2θ2ψ2Tin fixed frame are given in Figure 1(b). The equation of motion for the ith beam finite element can be written in the following form:(2)MTs,i+MRs,ix¨s,it+Cs,i+ωGs,ix˙s,it+Ks,ixs,it=Ws,i,where the superscript s refers to the shaft. MTs,i and MRs,i are the translational and rotary mass matrices of the shaft element. Cs,i, Gs,i, and Ks,i are the external damping and gyroscopic and stiffness matrices, respectively. Ws,i represents the gravity contribution. ω is the rotational speed of the shaft. The shaft is 1 m long and 0.04 m in diameter. The shaft is made of steel with Young’s modulus of elasticity Es=21011 N·m−2 and the density ρs=7800 kg·m−3. Table 1 lists the shaft properties.Table 1
Value of the physical parameters for the shaft.
Notation Description Value Rs Radius of the rotor shaft 0.02  m Ls Length of the rotor shaft 1  m Es Young’s modulus of elasticity 2.1   1 0 11 N⋅m−2 Gs Shear modulus 7.7   1 0 10 N⋅m−2 ρs Density 7800 kg⋅m−3 νs Poisson ratio 0.3


### 2.2. Disk Elements

Each disk is modeled as a rigid disk. The equation of motion for thejth disc finite element can be written in the following form:(3)MTd,j+MRd,jx¨d,jt+ωGd,jx˙d,jt=Fd,j+Wd,j,where the superscript d refers to the disk. xd,j defines the nodal displacement vector for the jth disk element that is defined as vwθψ in fixed frame. MTd,j and MRd,j define the translational and rotary mass matrices for each disk. Gd,j is the gyroscopic matrix. Fd,j and Wd,j correspond to the unbalance and gravity on the disk. Each disk is made of steel with Young’s modulus of elasticity Ed=21011 N·m−2 and the density ρd=7800 kg·m−3. Table 2 lists the disks properties.Table 2
Value of the physical parameters for the four disks.
Notation Description Value R1d Outer radius of the 1st disk 0.25  m h1d Thickness of the 1st disk 30  m m R2d Outer radius of 2nd the disk 0.1875  m h2d Thickness of the 2nd disk 15  m m R3d Outer radius of 3rd the disk 0.1875  m h3d Thickness of the 3rd disk 22.5  m m R4d Outer radius of 4th the disk 0.125  m h4d Thickness of 4th the disk 37.5  m m Ed Young’s modulus of elasticity 2.1   1 0 11 N⋅m−2 ρd Density 7800 kg⋅m−3 me Mass unbalance for 1st disk 0.01  k g de Eccentricity of the mass unbalance for the 1st disk 0.01  m


### 2.3. Bearing Supports

The two bearing supports are added at each end of the flexible rotor. For simplicity the combined support and bearing elements of the shaft are modeled as two-mode linear elastic spring elements. The elementary stiffness matrices are denoted byKb,1 and Kb,2 for the first and second bearing supports, respectively.


### 2.4. Out-of-Balance Forces

Generally the most significant lateral forces are caused by an imperfect distribution of mass in the rotor system which can be considered as one of the most common faults in a rotor problem. Unbalance appears when the centre of mass of a disc (or a shaft element) is out of alignment with the centre of rotation of the rotor system. This force that is generated when the rotor spins about its equilibrium is called out-of-balance force. It gives to the rotor a wobbling movement characteristic of vibration of rotating machinery. It can be usually modeled by an eccentric mass on a disc (or on the shaft). Considering the degrees of freedomvwθψT in fixed frame (see Figure 1(b)), the out-of-balance force for a given node of the finite-element rotor system can be given by(4)fe=medeω2cosωt+ϕmedeω2sinωt+ϕ00T,where me and de are, respectively, the mass unbalance and the eccentricity (i.e., the distance between unbalance and the centre of rotation of the system). ϕ defines the initial angular position of the mass unbalance in relation to z-axis. In the present study, a residual unbalance is assumed to be present on the first disk. Table 2 gives the value of the two unbalance parameters (i.e., mass unbalance and eccentricity).


### 2.5. General Equation of the Rotor System

The overall system can be formed by adding the previous separate contributions of the rotor (i.e., the shaft and the four disks) and the bearings. After assembling the different elements, the equations of the rotor system can be written as(5)Mx¨t+C+ωGx˙t+Kxt=ft+w,where x defines the vibrational response of the rotor system. M, K, and G are the mass, stiffness, and gyroscopic matrices of the complete system (i.e., the shaft, the four rigid discs, and the two bearing supports), respectively. For example, the mass matrix M includes the contributions of each beam element (i.e., MTs,i and MRs,i for i=1,…,10) and each disk element (i.e., MTd,j and MRd,j for j=1,…,4). The matrix C includes the effects of the shafts internal damping and damping of the supports. ω is the rotational speed of the rotor system. f corresponds to the previous linear forces (i.e., unbalance forces) applied to the rotor system. Gravity denoted by w is also included in this equation of the rotor system.


## 2.1. Shaft Elements

The shaft is discretized into 10 Euler beam finite elements with circular cross-sections. Each beam finite element has four degrees of freedom at each node (two lateral displacements and two rotations, the axial and torsional degrees of freedom being not considered). The nodal displacement vectors for theith shaft element defined as(1)xs,i=v1w1θ1ψ1v2w2θ2ψ2Tin fixed frame are given in Figure 1(b). The equation of motion for the ith beam finite element can be written in the following form:(2)MTs,i+MRs,ix¨s,it+Cs,i+ωGs,ix˙s,it+Ks,ixs,it=Ws,i,where the superscript s refers to the shaft. MTs,i and MRs,i are the translational and rotary mass matrices of the shaft element. Cs,i, Gs,i, and Ks,i are the external damping and gyroscopic and stiffness matrices, respectively. Ws,i represents the gravity contribution. ω is the rotational speed of the shaft. The shaft is 1 m long and 0.04 m in diameter. The shaft is made of steel with Young’s modulus of elasticity Es=21011 N·m−2 and the density ρs=7800 kg·m−3. Table 1 lists the shaft properties.Table 1
Value of the physical parameters for the shaft.
Notation Description Value Rs Radius of the rotor shaft 0.02  m Ls Length of the rotor shaft 1  m Es Young’s modulus of elasticity 2.1   1 0 11 N⋅m−2 Gs Shear modulus 7.7   1 0 10 N⋅m−2 ρs Density 7800 kg⋅m−3 νs Poisson ratio 0.3


## 2.2. Disk Elements

Each disk is modeled as a rigid disk. The equation of motion for thejth disc finite element can be written in the following form:(3)MTd,j+MRd,jx¨d,jt+ωGd,jx˙d,jt=Fd,j+Wd,j,where the superscript d refers to the disk. xd,j defines the nodal displacement vector for the jth disk element that is defined as vwθψ in fixed frame. MTd,j and MRd,j define the translational and rotary mass matrices for each disk. Gd,j is the gyroscopic matrix. Fd,j and Wd,j correspond to the unbalance and gravity on the disk. Each disk is made of steel with Young’s modulus of elasticity Ed=21011 N·m−2 and the density ρd=7800 kg·m−3. Table 2 lists the disks properties.Table 2
Value of the physical parameters for the four disks.
Notation Description Value R1d Outer radius of the 1st disk 0.25  m h1d Thickness of the 1st disk 30  m m R2d Outer radius of 2nd the disk 0.1875  m h2d Thickness of the 2nd disk 15  m m R3d Outer radius of 3rd the disk 0.1875  m h3d Thickness of the 3rd disk 22.5  m m R4d Outer radius of 4th the disk 0.125  m h4d Thickness of 4th the disk 37.5  m m Ed Young’s modulus of elasticity 2.1   1 0 11 N⋅m−2 ρd Density 7800 kg⋅m−3 me Mass unbalance for 1st disk 0.01  k g de Eccentricity of the mass unbalance for the 1st disk 0.01  m


## 2.3. Bearing Supports

The two bearing supports are added at each end of the flexible rotor. For simplicity the combined support and bearing elements of the shaft are modeled as two-mode linear elastic spring elements. The elementary stiffness matrices are denoted byKb,1 and Kb,2 for the first and second bearing supports, respectively.


## 2.4. Out-of-Balance Forces

Generally the most significant lateral forces are caused by an imperfect distribution of mass in the rotor system which can be considered as one of the most common faults in a rotor problem. Unbalance appears when the centre of mass of a disc (or a shaft element) is out of alignment with the centre of rotation of the rotor system. This force that is generated when the rotor spins about its equilibrium is called out-of-balance force. It gives to the rotor a wobbling movement characteristic of vibration of rotating machinery. It can be usually modeled by an eccentric mass on a disc (or on the shaft). Considering the degrees of freedomvwθψT in fixed frame (see Figure 1(b)), the out-of-balance force for a given node of the finite-element rotor system can be given by(4)fe=medeω2cosωt+ϕmedeω2sinωt+ϕ00T,where me and de are, respectively, the mass unbalance and the eccentricity (i.e., the distance between unbalance and the centre of rotation of the system). ϕ defines the initial angular position of the mass unbalance in relation to z-axis. In the present study, a residual unbalance is assumed to be present on the first disk. Table 2 gives the value of the two unbalance parameters (i.e., mass unbalance and eccentricity).


## 2.5. General Equation of the Rotor System

The overall system can be formed by adding the previous separate contributions of the rotor (i.e., the shaft and the four disks) and the bearings. After assembling the different elements, the equations of the rotor system can be written as(5)Mx¨t+C+ωGx˙t+Kxt=ft+w,where x defines the vibrational response of the rotor system. M, K, and G are the mass, stiffness, and gyroscopic matrices of the complete system (i.e., the shaft, the four rigid discs, and the two bearing supports), respectively. For example, the mass matrix M includes the contributions of each beam element (i.e., MTs,i and MRs,i for i=1,…,10) and each disk element (i.e., MTd,j and MRd,j for j=1,…,4). The matrix C includes the effects of the shafts internal damping and damping of the supports. ω is the rotational speed of the rotor system. f corresponds to the previous linear forces (i.e., unbalance forces) applied to the rotor system. Gravity denoted by w is also included in this equation of the rotor system.


## 3. Objectives of the Problem and Numerical Results

### 3.1. Some Details on the Rotor Problem

First of all, Figure2 highlights the Campbell diagram that represents the evolutions of the forcing frequencies as a function of the rotor speed for two deterministic cases (with kvb,2=khb,2=0.1106 N·m−1 for the first case and kvb,2=khb,2=2106 N·m−1 for the second case). The values of the physical parameters used correspond to those previously indicated in Section 2. It is recalled that the coincidence of an excitation frequency with a resonance frequency of the rotor system corresponds to the critical speeds. In other words, this indicates the rotational speeds of the rotor for which large synchronous responses occur. In the present study, four backward and forward critical speeds are present in the rotor speed range of interest. We recall that the terms “backward” and “forward” refer to the whirling motion of the rotor that can be in the same direction as the shafts rotation (“forward whirl”) or in the opposite direction (“backward whirl”). The ith backward and forward critical speeds (i.e., the backward whirl and forward whirl modes) diverge as the rotor speed increases (see Figure 2).Figure 2
Campbell diagram: continuous blue lines forkvb,2=khb,2=0.1106 N⋅m−1 and black dotted lines for kvb,2=khb,2=2106 N⋅m−1.In the following sections, we investigate the potential of kriging based models to be efficiently considered as representative descriptions of the dynamic behavior of the rotor system. For this, we will seek to predict the four backward and forward critical speeds as well as the associated vibration amplitudes according to two deterministic parameters (i.e., the vertical and horizontal bearing support of the second supportkvb,2 and khb,2) and in the presence of uncertainties for seven other selected parameters: the independent thicknesses of the four disks, Young’s modulus of the rotor shaft, and the vertical and horizontal bearing support of the first support. Variability intervals of these seven input parameters are indicated in Table 3. For the proposed study, the deterministic parameters kvb,2 and khb,2 will be set independently for fixed values ranging from 0.1106 N·m−1 to 2106 N·m−1. For each pair of these fixed parameters kvb,2;khb,2, uncertainties according to the seven other parameters will be taken into account.Table 3
Variability intervals of input parameters.
Parameters Variability intervals Thickness of the 1st disk (m) [ . 9 ; 1.1 ] h 1 d Thickness of the 2nd disk (m) [ . 9 ; 1.1 ] h 2 d Thickness of the 3rd disk (m) [ . 9 ; 1.1 ] h 3 d Thickness of the 4th disk (m) [ . 9 ; 1.1 ] h 4 d Young’s modulus of the rotor shaft (Pa or N⋅m−2) [ . 95 ; 1.05 ] E s Vertical bearing stiffness of the first support (N⋅m−1) [ . 95 ; 1.05 ] k v b , 1 Horizontal bearing stiffness of the first support (N⋅m−1) [ . 95 ; 1.05 ] k h b , 1For the reader’s understanding, the explicit functions of theith backward and forward critical speeds relative to the two deterministic parameters kvb,2 and khb,2 are defined as(6)CSbackwardifh1,…,4d,fEs,fkvb,1,fkhb,1kvb,2,khb,2CSforwardifh1,…,4d,fEs,fkvb,1,fkhb,1kvb,2,khb,2for i=1,…,4. fh1,…,4d, fEs, fkvb,1, and fkhb,1 correspond to the uniform probability density functions associated with the uncertainties of the seven selected parameters. If the users have information about the most likely values of the uncertain parameters, then these users may specify nonuniform probability density functions with the most likely values.Similarly, the explicit functions of the associated maximal amplitudes for the rotor system (i.e., the maximum amplitude of the rotor system at the backward and forward critical speeds) are defined as(7)Abackwardifh1,…,4d,fEs,fkvb,1,fkhb,1kvb,2,khb,2=maxAh,backwardifh1,…,4d,fEs,fkvb,1,fkhb,1kvb,2,khb,2,Av,backwardifh1,…,4d,fEs,fkvb,1,fkhb,1kvb,2,khb,2Aforwardifh1,…,4d,fEs,fkvb,1,fkhb,1kvb,2,khb,2=maxAh,forwardifh1,…,4d,fEs,fkvb,1,fkhb,1kvb,2,khb,2,Av,forwardifh1,…,4d,fEs,fkvb,1,fkhb,1kvb,2,khb,2where Ah,backwardi and Av,backwardi (Ah,forwardi and Av,forwardi, resp.) refer to the maximum amplitudes on the horizontal and vertical directions for the ith backward critical speeds (the ith forward critical speeds, resp.). In the following, these explicit functions will be noted CSbackwardi, CSforwardi, Abackwardi, and Aforwardi (for i=1,…,4) in order to simplify the notations in the paper.


### 3.2. Kriging Metamodeling

The main idea with kriging based metamodeling is related to the possibility of estimating unknown values of a function by using a small size set of input/output data previously generated from a predefined experiment design [15]. The key principle consists in the exploiting of the spatial correlations between the function values to correct the average behavior of the function described by a regression model. For a more general description of the kriging formalism, the associated metamodeling method will be presented in the sequel by considering a parameter-dependent output y(p), where y may represent one of the maximum amplitudes, the backward or forward critical speeds defined previously while p∈Rl is a parameter vector which will refer to the vertical and horizontal bearing support kvb,2 and khb,2 (hence l=2). Different kriging models dedicated to each one of the outputs are performed. It may be mentioned that y may also represent multiple outputs (e.g., multiple responses of the maximum amplitudes or/and the critical speeds). However, multivariate kriging is more complicated and may result in inferior predictions because additional kriging parameters must be estimated (for more details see [14]).Based on the kriging theory, the parameter-dependent outputy(p) can be metamodeled as follows:(8)yp=hpTβ+Zp.The latter expression is obtained by juxtaposing two mathematical objects. The first one consists in a linear combination of polynomial functions hi (for i=1,…,Q) defining a regression with a predetermined order. The second object designated by Z(·) is a realization of a Gaussian process which is assumed to be with a zero-mean value and a covariance function given by(9)EZx,Zp=σ2Rθ,x,p,where E· denotes the expectation operator.The covariance function (9) is characterized by both σ2 and R which define the process variance and the spatial correlation function, respectively. The latter which is a monotone function within [0,1] measures how any two points x and p in the design space are spatially close to each other. It is defined by(10)Rθ,x,p=∏j=1lRjθj,xj-pj,where θ∈Rl denotes an l-dimensional vector with scaling parameters.The following parameters have to be determined to completely define a kriging predictor: the parameter vectorsβ and θ plus the scalar σ2. The obtaining of those parameters requires some preliminary choices, namely, the order of the regression model, the spatial correlation function, and the experiment design. The optimal choices are not obvious to determine. However, the nature of the function to be approximated in terms of its smoothness and nonlinearity can guide the users. Concerning the spatial correlation functions, numerous types exist. Some of them are given in Table 4.Table 4
Correlation functions.
Correlation model Rj(θ,xj,pj) Exponential e x p ( - θ j | x j - p j | ) Gaussian exp(-θj(xj-pj)2) Linear max{0,1-θj|xj-pj|} Cubic 1 - 3 ξ j 2 + 2 ξ j 3 with ξj=min(1,θj|xj-pj|)As mentioned previously, defining a kriging predictor requiresβ, σ2, and θ to be determined. For this aim, an input/output of the simulation model is required. The latter is generated according to a predefined experiment design such as a Latin Hypercube sampling. The latter is used in this paper. So, by considering the input parameter vector p1,…,pNk and the associated output vector Y=y1,…,yNk, the kriging parameters can be calculated as follows.First, theθ parameter vector is given by the solution of the following optimization problem:(11)minθ=R1/Nkσ2,where R is a (Nk×Nk) correlation matrix with entries defined by Rij=R(θ,p(i),p(j)) and |·| denotes the determinant of a matrix.The optimal regression parameterβ∗ is given by(12)β∗=HTR-1H-1HTR-1Yas the optimal solution to the following classical least square problem:(13)Hβ≃Y,where H is a (Q×Nk) matrix whose entries are given by Hij=hip(j).The variance processσ2 can then be estimated through(14)σ2=1NkY-Hβ∗TR-1Y-Hβ∗.All parameters being determined, the kriging predictor is completely defined and then can be written as follows:(15)Yx=hTxβ∗+rxTR-1Y-Hβ∗,where rx=Rθ,p(1),x,…,Rθ,pNk,xT is the vector of the correlation coefficients between the new input vector x and the Nk old sampled input vectors (p(i))i∈1,…,Nk.


### 3.3. Reference Solutions

First of all, a Monte Carlo simulation was carried out by considering two deterministic inputskvb,2 and khb,2 in the design interval of interest (i.e., these inputs are deterministic because they are assumed to be controlled by the users) and seven random inputs with uniform distributions (i.e., these inputs are not controlled by the users but are given by the environment of the rotor system, as seen in Table 3). For the reference solutions we examine 100×100 combinations of the two deterministic inputs by sampling the two-dimensional space parameter (kvb,2,khb,2); for each combination, 10000 simulations with the seven noise factors are performed (i.e., N=10000 samples for each combination of the two decision variables (kvb,2,khb,2) are calculated). Then for each deterministic couple (kvb,2,khb,2) the mean and the variance of the four backward and forward critical speeds as well as the associated vibration amplitudes are estimated. These moments are determined by using the classical formula for the mean value and variance given by(16)MeanYx≈1N∑k=1NYxkVarianceYx≈1N∑k=1NYxk-MeanYx2,where N is the number of samples for each combination of the two decision variables kvb,2 and khb,2.All these results are given in Figures3, 4, 5, and 6, respectively. Several comments may be made. First of all the numerical results show that both the values of critical speeds and the associated maximum magnitudes of rotor response are drastically affected by the deterministic couple (kvb,2,khb,2). For example, the mean values of the first and second backward critical speeds changed from 16 Hz to 33 Hz and from 44 Hz to 65 Hz, respectively. Similarly, the mean values of the associated maximum magnitudes vary considerably for all critical speeds as illustrated in Figure 5. The map shapes concerning the mean values and variances of the critical speeds are quite simple. A monotonic increasing of the mean values of each critical speed with the increase of the physical parameters kvb,2 and khb,2 is observed (see Figures 3(a)–3(h)). Evolution of the variance of each critical speed can also be simple with an increase (see Figures 4(a), 4(b), and 4(f)) or decrease (see Figures 4(g) and 4(h)) according to the evolution of kvb,2 or khb,2 or more complex as illustrated in Figures 4(c), 4(d), and 4(e).Figure 3
Referential map of the four backward and forward critical speeds: mean. (a)CSbackward1; (b) CSforward1; (c) CSbackward2; (d) CSforward2; (e) CSbackward3; (f) CSforward3; (g) CSbackward4; and (h) CSforward4.
(a) (b) (c) (d) (e) (f) (g) (h)Figure 4
Referential map of the four backward and forward critical speeds: variance. (a)CSbackward1; (b) CSforward1; (c) CSbackward2; (d) CSforward2; (e) CSbackward3; (f) CSforward3; (g) CSbackward4; and (h) CSforward4.
(a) (b) (c) (d) (e) (f) (g) (h)Figure 5
Referential map of the vibration amplitudes associated with the four backward and forward critical speeds: mean. (a)Abackward1; (b) Aforward1; (c) Abackward2; (d) Aforward2; (e) Abackward3; (f) Aforward3; (g) Abackward4; and (h) Aforward4.
(a) (b) (c) (d) (e) (f) (g) (h)Figure 6
Referential map of the vibration amplitudes associated with the four backward and forward critical speeds: variance. (a)Abackward1; (b) Aforward1; (c) Abackward2; (d) Aforward2; (e) Abackward3; (f) Aforward3; (g) Abackward4; and (h) Aforward4.
(a) (b) (c) (d) (e) (f) (g) (h)Regarding the evolution of amplitudes, the observed behavior is generally more complex with an increase or decrease of the mean and variance for each critical speed with increasing of each parameterkvb,2 or khb,2 (see Figures 5 and 6). For certain specific cases, the asymmetry of the bearing support (i.e., kvb,2≠khb,2) introduces further modifications to the vibration of the rotor. The mean values of the magnitude responses of the asymmetric rotor (i.e., kvb,2≠khb,2) at each backward critical speed are more important than the mean value of the vibration of the symmetric rotor (i.e., kvb,2=khb,2). These effects of the asymmetric support are highly visible in Figures 5(a), 5(c), 5(d), and 5(f). Finally, it is also observed that the evolution of the mean amplitudes of the asymmetric rotor can be difficult to predict intuitively (see, e.g., the evolutions of the maximum amplitudes for the second forward critical speeds in Figure 5(d)). Furthermore, the uncertainties introduced induce a nonnegligible variability on the maximum magnitudes of rotor response at each critical speed, as indicated in Figure 6. More or less complex evolutions of the variability on the maximum magnitudes are still observed.In conclusion, these results will serve as reference solutions for the validation of the kriging surrogate model. The objective being to be able to reproduce all the previously described evolutions concerning the mean and the variance of the four backward and forward critical speeds as well as the associated vibration amplitudes.As mentioned previously, reference solutions are performed by considering100×100 combinations of these two deterministic inputs kvb,2 and khb,2) with 10000 samples for each combination. This leads to a total of 108 calculations for the reference solutions. It can be noted that this number of calculations is prohibitive from a practical point of view. So the basic idea is to be able to drastically reduce this number of samples of the two controllable deterministic inputs kvb,2 and khb,2 by creating continuous surfaces from a predefined set of points. In other words, we will seek to predict the value at the point where the value is unknown by using an interpolation method (and by considering the values at a predefined set of sample points). In the following we will propose using a kriging method that allows the prediction of a variable of interest at a location where it has not been calculated, by taking some form of weighted average of the values at surrounding points.In conclusion, these results will serve as reference solutions for the validation of the kriging surrogate model. The objective being to be able to reproduce all the previously described evolutions concerning the mean and the variance of the four backward and forward critical speeds as well as the associated vibration amplitudes.


### 3.4. Kriging Based Predictions

In this subsection kriging models are used to estimate the referential maps of the four backward and forward critical speeds (shown in Figures3 and 4 for the reference model). In this perspective, the DACE software presented in [21] is exploited. Otherwise, there exist numerous other software packages that may be used for the same aim. Results are not necessarily the same. A comparative study of several software packages defined for kriging metamodeling can be found in [22].As mentioned previously, the construction of a kriging predictor requires preliminary choices of the regression order, the spatial correlation function, and the size of the experiment design. Hence, a study of the influence of these parameters on the accuracy of kriging based predictions of critical speedsCSbackwardi and CSforwardi and the vibration amplitudes Abackwardi and Aforwardi is performed. In a first step, one random experiment design with Nk=600 is fixed and the influence of the regression order and the spatial correlation function is analyzed. In a second step, the regression order and the spatial correlation function are fixed while the size of the experiment design is varied to observe the influence of the size of the experiment design.


#### 3.4.1. Prediction of Critical SpeedsCSbackwardi and CSforwardi

With a random plan of sizeNk=600, the mean values and variances of the critical speeds CSbackwardi and CSforwardi and the vibration amplitudes Abackwardi and Aforwardi are estimated against (kvb,2,khb,2) through kriging metamodels considering zero-, first-, and second-order regressions and the linear, Gaussian, exponential, and cubic correlation functions.The maximum and mean of the relative errors which are defined byerrorRMax and errorRMean (for both the mean and variance of each backward or forward critical speed) between the kriging metamodels and the reference model are given by(17)errorRMax=maxXk,ref-Xk,krigXk,refk=1,…,NerrorRMean=1N∑k=1NXk,ref-Xk,krigXk,ref,where Xref corresponds to the previous reference results. Xkrig are the results of the metamodel prediction value. N defines the number of points in the design space (here N=100∗100 as previously indicated in Section 3.3). It may be worth mentioning that each previous criterion errorRMax or errorRMean returns one unique value each for all the design space.Finally the accuracy of the metamodels is compared by using relative generalization error [23], which is defined as (18)errorG=1/N∑k=1NXk,ref-Xk,krig2varXref,where var denotes the variance. As for the two previous criteria, the relative generalization error errorG returns one unique value in the design space according to the variable of interest (i.e., one is interested on the calculation of the mean or variance of a given critical speed in the present study). In several cases, the size Nk of the design experiment is determined by 10 times the number of the considered parameters [24]. In this study, Nk was fixed in such a way that the mean relative errors corresponding to the predicted critical speeds and amplitudes do not exceed a predefined error level (1% in our study).Results of the two criteriaerrorRMax and errorRMean for the mean and variance on the design space are given in Tables 5 and 6 for all the estimated critical speeds, respectively. Tables 7 and 8 give the associated relative generalization error. The values marked in bold (in underlines, resp.) correspond to the minimum (the maximum, resp.) of the error values for each criterion (for a given critical speed and a given regression order) in order to visualize the best (the worst, resp.) spatial correlation functions to predict the desired quantities with kriging models. The values marked in italic bold correspond to the best metamodel. So it indicates the best spatial correlation functions and regression order in order to predict the mean and variance of each critical speed. In this part of the paper illustrations of the maps, which are obtained by using the best kriging metamodels, are not given because they are very similar to the results displayed previously for the reference calculations.Table 5
% of the relative errors(errorRMax;errorRMean) of the mean of critical speeds for different regression orders and correlation functions. Nonsymmetric case with 600 samples.
Gaussian Linear Exponential Cubic Zero-order regression CSbackward1 (17.86; 0.29) (28.58; 0.65) (25.30; 0.36) (40.20; 0.81) CSforward1 (10.94;0.10) (10.33;0.14) (12.82; 0.11) (8.58; 0.10) CSbackward2 (2.55; 0.02) (4.54; 0.16) (3.02; 0.08) (10.40; 0.27) CSforward2 (3.37;0.02) (4.07; 0.08) (0.99; 0.03) (11.01; 0.14) CSbackward3 (0.00; 0.00) (2.04; 0.06) (0.66; 0.02) (4.65; 0.11) CSforward3 (0.00; 0.00) (2.22; 0.05) (0.74; 0.02) (5.15; 0.10) CSbackward4 (0.00; 0.00) (0.54; 0.01) (0.18; 0.00) (1.23; 0.03) CSforward4 (0.00; 0.00) (0.31; 0.01) (0.10; 0.00) (0.70; 0.02)  First-order regression CSbackward1 (11.98; 0.25) (23.81; 0.50) (23.13; 0.30) (27.46; 0.56) CSforward1 (6.51; 0.09) (12.64;0.15) (14.28; 0.12) (7.42;0.08) CSbackward2 (1.54; 0.02) (1.93; 0.07) (2.04; 0.04) (3.66; 0.10) CSforward2 (0.86; 0.01) (1.33;0.04) (1.59; 0.03) (1.69; 0.03) CSbackward3 (0.00; 0.00) (0.39; 0.01) (0.27; 0.00) (0.70; 0.01) CSforward3 (0.00; 0.00) (0.11; 0.00) (0.07; 0.00) (0.26; 0.01) CSbackward4 (0.00; 0.00) (0.05; 0.00) (0.04; 0.00) (0.10; 0.00) CSforward4 (0.00; 0.00) (0.01; 0.00) (0.00; 0.00) (0.01; 0.00)  Second-order regression CSbackward1 (8.98; 0.20) (20.97;0.38) (21.40; 0.24) (16.80; 0.35) CSforward1 (6.35; 0.08) (12.26;0.14) (13.91; 0.12) (8.45; 0.08) CSbackward2 (0.63; 0.01) (1.22;0.04) (1.57; 0.03) (1.11; 0.03) CSforward2 (0.60; 0.01) (1.20;0.03) (1.43; 0.02) (1.14; 0.02) CSbackward3 (0.00; 0.00) (0.07; 0.00) (0.06; 0.00) (0.11; 0.00) CSforward3 (0.00; 0.00) (0.05; 0.00) (0.04; 0.00) (0.10; 0.00) CSbackward4 (0.00; 0.00) (0.00; 0.00) (0.00; 0.00) (0.00; 0.00) CSforward4 (0.00; 0.00) (0.00; 0.00) (0.00; 0.00) (0.00; 0.00)Table 6
% of the relative errors(errorRMax,errorRMean) on the variance of the critical speeds for different regression orders and correlation functions. Nonsymmetric case with 600 samples.
Gaussian Linear Exponential Cubic Zero-order regression CSbackward1 (151.85; 1.18) (76.12; 1.66) (58.14; 0.92) (166.04; 2.84) CSforward1 (48.74;0.78) (21.74; 0.67) (12.19; 0.39) (77.41; 0.72) CSbackward2 (10.58; 0.19) (11.05;0.33) (11.89; 0.22) (7.44; 0.17) CSforward2 (12.40; 0.23) (7.52;0.41) (8.02; 0.27) (6.32; 0.21) CSbackward3 (0.00; 0.00) (0.67; 0.02) (0.29; 0.01) (1.29; 0.04) CSforward3 (0.00; 0.00) (2.29; 0.07) (0.78; 0.02) (5.00; 0.14) CSbackward4 (0.00; 0.00) (1.01; 0.03) (0.69; 0.01) (2.27; 0.05) CSforward4 (0.00; 0.00) (0.49; 0.01) (0.16; 0.00) (1.11; 0.03)  First-order regression CSbackward1 (72.82; 0.72) (46.36; 0.86) (46.19; 0.55) (61.98;1.10) CSforward1 (39.81; 5.32) (13.20; 0.52) (18.12; 0.41) (12.17; 0.30) CSbackward2 (8.74; 0.20) (11.41;0.33) (12.06; 0.23) (7.91; 0.17) CSforward2 (10.79; 0.18) (9.42;0.40) (9.06; 0.26) (9.33;0.17) CSbackward3 (0.00; 0.00) (0.72; 0.01) (0.60; 0.01) (1.23; 0.02) CSforward3 (0.00; 0.00) (0.78; 0.01) (0.68; 0.01) (1.31; 0.01) CSbackward4 (0.00; 0.00) (0.34; 0.01) (0.25; 0.00) (0.60; 0.01) CSforward4 (0.00; 0.00) (0.02; 0.00) (0.01; 0.00) (0.03; 0.00)  Second-order regression CSbackward1 (15.47; 0.29) (35.98;0.50) (39.75; 0.38) (27.97; 0.33) CSforward1 (148.91; 4.25) (12.84; 0.42) (16.13; 0.37) (16.10;0.25) CSbackward2 (6.95; 0.15) (11.16; 0.32) (11.05; 0.21) (8.06;0.13) CSforward2 (7.37; 0.15) (7.06;0.31) (7.70; 0.23) (6.55; 0.13) CSbackward3 (0.00; 0.00) (0.19; 0.01) (0.16; 0.00) (0.33; 0.01) CSforward3 (0.00; 0.00) (0.08; 0.00) (0.06; 0.00) (0.13; 0.00) CSbackward4 (0.00; 0.00) (0.03; 0.00) (0.03; 0.00) (0.05; 0.00) CSforward4 (0.00; 0.00) (0.00; 0.00) (0.00; 0.00) (0.01; 0.00)Table 7
Relative generalization errorerrorG on the mean of the critical speeds for different regression orders and correlation functions. Nonsymmetric case with 600 samples.
Gaussian Linear Exponential Cubic Zero-order regression CSbackward1 1.87 10 −3 6.3710-3 3.1810-3 1.49   10 - 2 _ CSforward1 5.4710-3 5.3610-3 7.30   10 - 3 _ 4.09 10 −3 CSbackward2 6.14 10 −5 8.2310-4 2.8010-4 4.44   10 - 3 _ CSforward2 2.0610-4 5.0910-4 6.69 10 −5 3.43   10 - 3 _ CSbackward3 5.36 10 −12 5.8710-4 1.0710-4 3.54   10 - 3 _ CSforward3 7.71 10 −12 6.0310-4 9.83 10−5 3.36   10 - 3 _ CSbackward4 1.36 10 −10 5.8310-4 9.3310-5 3.29   10 - 3 _ CSforward4 3.90 10 −10 5.7510-4 8.8610-5 3.37   10 - 3 _  First-order regression CSbackward1 1.38 10 −3 3.8510-3 2.4010-3 6.79   10 - 3 _ CSforward1 3.2410-3 7.4910-3 9.54   10 - 3 _ 2.08 10 −3 CSbackward2 2.82 10 −5 1.4610-4 8.7010-5 5.20   10 - 4 _ CSforward2 2.86 10 −5 1.1210-4 9.5110-5 2.06   10 - 4 _ CSbackward3 8.26 10 −13 1.1610-5 3.8910-6 4.89   10 - 5 _ CSforward3 2.09 10 −13 2.7810-6 8.7610-7 1.48   10 - 5 _ CSbackward4 1.03 10 −10 3.1310-6 1.2110-6 1.16   10 - 5 _ CSforward4 3.51 10 −10 9.4010-8 3.7710-8 3.26   10 - 7 _  Second-order regression CSbackward1 8.93 10 −4 2.2510-3 1.7010-3 2.45   10 - 3 _ CSforward1 2.19 10 −3 6.6010-3 9.14   10 - 3 _ 3.2610-3 CSbackward2 8.02 10 −6 4.86   10 - 5 _ 4.6410-5 3.35 10−5 CSforward2 1.15 10 −5 7.26   10 - 5 _ 6.2710-5 5.2010-5 CSbackward3 2.51 10 −13 3.0210-7 1.9310-7 6.46   10 - 7 _ CSforward3 1.74 10 −13 4.4910-7 2.1310-7 1.80   10 - 6 _ CSbackward4 1.23 10 −10 4.6910-10 3.7610-10 9.57   10 - 10 _ CSforward4 3.37 10 −10 6.9610-9 2.6310-9 2.52   10 - 8 _Table 8
Relative generalization errorerrorG on the variance of the critical speeds for different regression orders and correlation functions. Nonsymmetric case with 600 samples.
Gaussian Linear Exponential Cubic Zero-order regression CSbackward1 1.2110-3 1.0810-3 3.97 10 −4 4.61   10 - 3 _ CSforward1 4.07   10 - 3 _ 1.5110-3 7.07 10 −4 3.8110-3 CSbackward2 9.9010-3 1.29   10 - 2 _ 1.0810-2 4.26 10 −3 CSforward2 7.00   10 - 3 _ 5.9710-3 4.1410-3 2.99 10 −3 CSbackward3 7.91 10 −12 2.8310-4 5.4610-5 1.56   10 - 3 _ CSforward3 4.16 10 −12 4.5310-4 7.6310-5 2.76   10 - 3 _ CSbackward4 6.50 10 −10 7.0710-4 1.7110-4 3.34   10 - 3 _ CSforward4 3.08 10 −9 5.7010-4 8.7910-5 3.30   10 - 3 _  First-order regression CSbackward1 3.7510-4 2.9010-4 1.71 10 −4 5.92   10 - 4 _ CSforward1 6.49   10 - 2 _ 9.6310-4 7.6310-4 5.33 10 −4 CSbackward2 8.1610-3 1.33   10 - 2 _ 1.1010-2 4.93 10 −3 CSforward2 4.1710-3 6.19   10 - 3 _ 4.3310-3 2.41 10 −3 CSbackward3 1.81 10 −11 1.1510-4 6.7010-5 3.86   10 - 4 _ CSforward3 5.73 10 −12 1.6210-5 1.0110-5 4.84   10 - 5 _ CSbackward4 7.46 10 −10 3.3310-5 1.1910-5 1.33   10 - 4 _ CSforward4 3.28 10 −9 4.2510-7 1.4510-7 1.94   10 - 6 _  Second-order regression CSbackward1 6.40 10 −5 1.31   10 - 4 _ 1.1310-4 6.9210-5 CSforward1 3.46   10 - 2 _ 7.0610-4 6.7210-4 3.25 10 −4 CSbackward2 4.8910-3 1.29   10 - 2 _ 9.4310-3 3.54 10 −3 CSforward2 2.5010-3 4.31   10 - 3 _ 3.5410-3 1.49 10 −3 CSbackward3 7.46 10 −12 1.5510-5 8.4310-6 5.14   10 - 5 _ CSforward3 9.90 10 −13 6.3610-7 3.5010-7 1.74   10 - 6 _ CSbackward4 7.25 10 −10 1.3410-7 1.1210-7 2.63   10 - 7 _ CSforward4 3.15 10 −9 8.3010-8 3.5310-8 3.05   10 - 7 _First, if some marginal cases are excluded, it can also be noted that, whatever the selected kriging patterns, satisfactory results are generally observed for all the three criteriaerrorRMax, errorRMean, and errorG. It also appears that the best kriging model (in terms of the prediction of the mean and variance of the critical speeds) is not always associated with the same regression order or the same correlation functions. These results illustrate the difficulty of prior choices of both the regression order and the spatial correlation function for a mechanical system. However, it is shown that the best kriging model is generally obtained with the Gaussian correlation function and the second-order regression. In practice the users may try different kriging metamodels (if these metamodels require much less computer time than the simulation code) and select the metamodel that gives best results for a separate (small) test set of inputs.Similarly it appears that, according to the criterion chosen (errorRMax, errorRMean, or errorG), the best kriging model identified for a given output data to be predicted (i.e., a specific critical speed and according to the choice between its means or variances) is not necessarily the same. Moreover, it is also clearly shown that all the values of the maximum and mean of errorRMax, errorRMean, or errorG are not necessarily of the same order of magnitude according to the critical speeds considered. For example, a maximum relative error errorRMax of 8.63% is registered for the mean value of CSbackward1 whereas this value is inferior to 1% for CSbackward2 and CSforward2 and close to 1% for CSbackward3, CSforward3, CSbackward4, and CSforward4. Similarly the values of errorG for both the mean and variance (indicated in Tables 7 and 8) for the first four critical speeds CSbackward1, CSforward1, CSbackward2, and CSforward2 are very different (i.e., superior by a factor of at least 105) from those of the last four critical speeds CSbackward3, CSforward3, CSbackward4, and CSforward4. All these facts show that the definition of relevant criteria for a given mechanical system and the selection of the output data to be predicted play an important role in the construction of a kriging model.In conclusion the efficiency of the proposed kriging approach for the estimation of the four backward and forward critical speeds of the rotor system in presence of uncertainties is validated through the use of three different criteria.


#### 3.4.2. Discussion on Exploiting the Symmetry Property of the Problem Addressed

In the previous section, the number of samples used for the calculation of all the kriging predictors has been chosen as equal toNk=600. It may be noted that this number is not optimal due to the fact that Nk is not necessarily the same to ensure the same accuracy for the predictions of the mean values and variances of the critical speeds. In the previous section, this number has been chosen to have an error errorRMean of less than 1% for both the mean and the variance of all the critical speeds of interest.Even if a convergence study with respect toNk for each function needed to be approximated is not the objective of this study (which, for recall, aims to show the possibility of exploiting kriging metamodeling in the analysis of rotor dynamics), it can be said that the necessary number of samples will be smaller if the map to be approximated by the kriging model will have a simple evolution in terms of the parameters considered. This requires an a priori knowledge of the more or less complex evolutions of the quantities processed.So, in the following, a specific study on the influence of the number of samples will be conducted in conjunction with a specific property of the system to drastically reduce the number of samples to be considered. Indeed it can be noted that the size of the experiment design can be drastically reduced by exploiting the physics of the problem. More particularly, exploiting the property of symmetry of the rotor system (between the horizontal and vertical directionskhb,2 and kvb,2) may allow reducing the size of the experiment design, and the input combinations (design points) within that design, without decreasing the accuracy of the corresponding kriging predictions. First, the main idea is to divide the design space according to the identified symmetry property and then to generate an experiment design within the reduced design space. Kriging metamodels are then used to mimic the functions within the new design spaces. To illustrate this on the proposed rotor system, the previous design space (kvb,2,khb,2) with 0.1106 N·m−1≤kvb,2≤2106 N·m−1 and 0.1106 N·m−1≤khb,2≤2106 N·m−1 will be restricted to the following design space: 0.1106 N·m−1 ≤ kvb,2≤khb,2 and kvb,2≤khb,2≤2106 N·m−1. More precisely three cases will be investigated.Case 1.
A total of 350 points are distributed as follows: 50 samples are distributed uniformly (between0.1106 N·m−1 and 2106 N·m−1) at the boundary of the reduced design space (for khb,2=kvb,2) and 300 samples are chosen according to a LHS process restricted to the following design space: 0.1106 N·m−1≤kvb,2<khb,2 and kvb,2<khb,2≤2106 N·m−1.Case 2.
A total of 120 points are distributed as follows: 20 samples are distributed uniformly (between0.1106 N·m−1 and 2106 N·m−1) at the boundary of the reduced design space (for khb,2=kvb,2) and 100 samples are chosen according to a LHS process restricted to the following design space: 0.1106 N·m−1≤kvb,2<khb,2 and kvb,2<khb,2≤2106 N·m−1.Case 3.
A total of 70 points are distributed as follows: 20 samples are distributed uniformly (between0.1106 N·m−1 and 2106 N·m−1) at the boundary of the reduced design space (for khb,2=kvb,2) and 50 samples are chosen according to a LHS process restricted to the following design space: 0.1106 N·m−1≤kvb,2<khb,2 and kvb,2<khb,2≤2106 N·m−1.These three cases of study will enable us to test the performance associated with taking into account the property of symmetry and also to examine the convergence according to number of samples. As previously explained a convergence study with respect toNk for each of the critical speeds has to be performed in order to optimize the number of points required for accurate approximation. However the objective of the present part of the study is only to illustrate the potential of including thoughtful choices in the kriging simulation.It turns that exploiting the symmetry property enhances the accuracies of kriging predictions of the mean values and the variances of the critical speeds observed previously withNk=600. For example, the maximum errors of the three criteria errorRMax, errorRMean, and errorG by considering the best kriging models are given in Table 9 for the previous case (Section 3.4) and the three cases dealt with in this part of the paper. It is shown that the best kriging metamodel with only 120 samples and by using the symmetry property gives better results for most of the tested criteria than the best kriging metamodel without considering the symmetry property. Thus it is possible to considerably reduce the number of samples required to have predictive models if one is able to employ certain characteristics directly deriving from the physics of the mechanical system studied.Table 9
Maximum errors for all the critical speeds by considering the best kriging models.
e r r o r R M a x ; e r r o r R M e a n e r r o r G Mean values Case 1: 350 samples (1.20; 0.02) 4.8410-5 Case 2: 120 samples (1.99; 0.06) 2.1210-4 Case 3: 70 samples (18.61; 0.31) 3.6310-3  Variance values Case 1: 350 samples (4.80; 1.55) 6.5710-2 Case 2: 120 samples (12.27; 1.65) 6.6910-2 Case 3: 70 samples (30.99; 1.92) 6.6010-2Furthermore, by comparing the results of the three cases, it is clearly shown that increasing the number of samples significantly improves predictions of kriging models. To illustrate this, Table9 summarizes the maximum value of each criterion of errorRMax, errorRMean, and errorG for the prediction of the mean and variance of all the critical speeds.For the interested reader, AppendixA gives all the calculations of the three criteria errorRMax, errorRMean, and errorG for different regression orders and correlation functions. A detailed analysis of all these results leads to the following conclusions:(i)
Kriging metamodeling is a valuable tool to provide valuable approximate results in the field of rotordynamics. It is very interesting for practical mechanical problems where limited outputs data is available (due to expensive simulation codes).(ii)
The best kriging metamodel, for the prediction of the mean or the variance of a given critical speed, does not necessarily correspond to a specific correlation function or regression order. However, for the problem under study, it may be noted that the Gaussian correlation function and the second-order regression often correspond to the best compromise to obtain the most accurate kriging model versus all the criteria tested.(iii)
The necessary number of samples to reach a minimum according to a given criterion is not identical for each critical speed. For example, considering 70 samples is sufficient to have the errorerrorRMax less than 1% for the mean of CSbackward3, CSforward3, CSbackward4, and CSforward4 whereas 120 samples (350 samples, resp.) are needed for CSbackward2 and CSforward2 (CSforward1, resp.); the error errorRMax of the mean of CSbackward1 is still greater than 1% even with 350 samples.In conclusion, even if the choices of the correlation function and the regression order are crucial for an accurate kriging predictor, it is demonstrated that a drastic reduction in the number of samples necessary for a correct prediction of all the critical speeds of the rotor system is possible.


#### 3.4.3. Prediction of the Vibration AmplitudesAbackwardi and Aforwardi

The same study is performed for the estimation of the vibration amplitudesAbackwardi and Aforwardi by using kriging metamodels without considering the symmetry property. Indeed, the properties of symmetry on results of the vibratory responses are no longer valid due to the presence of gravity. Moreover, only the generalized error criterion errorG, which is defined in (18), can be used in a robust manner. Indeed, vibratory responses of very small amplitudes for certain configurations (i.e., the choice of parameters kvb,2 and khb,2) can generate very high error levels (the two criteria errorRMax and errorRMean are then maximized because of a division by a quantity close to zero).Results of the relatively generalized errors for the mean and variance on the design space are given in Tables10 and 11 for the prediction of all the vibration amplitudes Abackwardi and Aforwardi (for i=1,…,4). We recall that the values marked in bold (in underlines, resp.) correspond to the best (the worst) metamodel for a given regression order. As said previously, results show that the choice of spatial correlation functions and regression orders can affect efficiency of kriging metamodeling. Similarly, the best spatial correlation function and the best regression order are not necessarily identical for each of the vibration amplitudes (i.e., Abackwardi and Aforwardi for i=1,…,4). Moreover, for a given vibration amplitude, the best spatial correlation function and regression order for the mean value are not necessarily those which will also make it possible to minimize the error of the variance. This shows the difficulty of choosing a priori the best spatial correlation function and the best regression order for the prediction of the maximal amplitudes of a rotor system. However it can be noted that a cubic correlation with a second-order regression works well in most cases.Table 10
Relative generalization errorerrorG on the mean of the amplitude for different regression orders and correlation functions. Nonsymmetric case with 600 samples.
Gaussian Linear Exponential Cubic Zero-order regression Abackward1 9.43   10 - 1 _ 6.5210-2 6.07 10 - 2 8.4310-2 Aforward1 9.45   10 - 1 _ 2.1510-3 1.39 10 - 3 6.3710-3 Abackward2 1.1110-2 1.89   10 - 2 _ 1.1810-2 5.12 10 - 3 Aforward2 1.7810-2 1.79   10 - 2 _ 1.2510-2 7.05 10 - 3 Abackward3 1.0810-3 1.1010-3 5.58 10 - 4 3.63   10 - 3 _ Aforward3 3.45 10 −3 4.69   10 - 3 _ 4.2610-3 4.0510-3 Abackward4 5.3110-4 9.2710-4 4.43 10 - 4 3.29   10 - 3 _ Aforward4 2.4810-3 1.2710-3 1.09 10 - 3 3.60   10 - 3 _  First-order regression Abackward1 7.90   10 - 1 _ 7.3510-2 5.95 10 - 2 7.8610-2 Aforward1 5.29   10 - 2 _ 1.6810-3 1.42 10 - 3 2.8410-3 Abackward2 6.3410-3 1.0110-2 1.01   10 - 2 _ 3.41 10 - 3 Aforward2 3.92   10 - 1 _ 1.2410-2 1.1510-2 6.14 10 - 3 Abackward3 9.7810-4 6.4810-4 4.08 10 - 4 9.81   10 - 4 _ Aforward3 3.0810-3 4.1510-3 4.23   10 - 3 _ 1.50 10 - 3 Abackward4 4.1410-4 2.3410-4 1.12 10 - 4 7.45   10 - 4 _ Aforward4 1.0410-3 1.6010-3 1.69   10 - 3 _ 6.54 10 - 4  Second-order regression Abackward1 5.34   10 - 1 _ 5.9610-2 5.55 10 −2 6.9710-2 Aforward1 1.23   10 - 2 _ 1.4210-3 1.22 10 −3 1.6910-3 Abackward2 1.43   10 - 1 _ 9.7310-3 9.2210-3 3.28 10 −3 Aforward2 2.39   10 - 1 _ 9.3110-3 1.0910-2 5.66 10 −3 Abackward3 2.2510-4 4.02   10 - 4 _ 3.2410-4 1.52 10 −4 Aforward3 1.7810-3 2.6910-3 3.50   10 - 3 _ 1.35 10 −3 Abackward4 9.7610-5 1.36   10 - 4 _ 9.6810-5 7.95 10 −5 Aforward4 6.1010-4 1.0010-3 1.33   10 - 3 _ 5.22 10 −4Table 11
Relative generalization errorerrorG on the variance of the amplitude for different regression orders and correlation functions. Nonsymmetric case with 600 samples.
Gaussian Linear Exponential Cubic Zero-order regression Abackward1 9.42   1 0 - 1 _ 1.8610-2 1.3810-2 9.30 1 0 - 3 Aforward1 9.42   1 0 - 1 _ 4.0910-2 3.51 1 0 - 2 4.1710-2 Abackward2 9.35   1 0 - 1 _ 3.3910-1 1.7910-1 1.53 10 −1 Aforward2 9.44   1 0 - 1 _ 1.4510-1 8.8510-2 6.52 1 0 - 2 Abackward3 3.75 1 0 - 4 1.0110-3 4.8210-4 4.69   1 0 - 3 _ Aforward3 5.2810-3 5.36   1 0 - 3 _ 4.2810-3 3.39 1 0 - 3 Abackward4 1.0910-3 1.7010-3 9.23 1 0 - 4 4.39   1 0 - 3 _ Aforward4 4.58   1 0 - 3 _ 1.4110-3 1.07 1 0 - 3 3.9510-3  First-order regression Abackward1 3.42   1 0 - 1 _ 1.3110-2 1.2810-2 8.56 1 0 - 3 Aforward1 3.18   1 0 - 1 _ 3.2810-2 3.28 1 0 - 2 3.7110-2 Abackward2 7.53   1 0 - 1 _ 3.4510-1 1.8110-1 1.55 1 0 - 1 Aforward2 3.39   1 0 - 1 _ 1.4410-1 7.3910-2 5.28 10 −2 Abackward3 1.2710-4 1.7810-4 5.42 1 0 - 5 1.13   1 0 - 3 _ Aforward3 3.9710-3 5.28   1 0 - 3 _ 4.5010-3 2.02 1 0 - 3 Abackward4 5.0010-4 3.8910-4 1.89 1 0 - 4 1.21   1 0 - 3 _ Aforward4 1.96   1 0 - 3 _ 1.6010-3 1.6710-3 7.52 1 0 - 4  Second-order regression Abackward1 2.77   1 0 - 1 _ 1.2810-2 1.2610-2 8.44 10 −3 Aforward1 2.39   1 0 - 1 _ 3.10 10 −2 3.2110-2 3.5710-2 Abackward2 6.20   1 0 - 1 _ 3.4710-1 1.8810-1 1.67 1 0 - 1 Aforward2 3.08   1 0 - 1 _ 1.4410-1 7.4010-2 5.40 1 0 - 2 Abackward3 5.18   1 0 - 5 2.68 10 −5 2.7910-5 7.82   1 0 - 5 _ Aforward3 1.44   1 0 - 1 _ 3.3310-3 3.7210-3 1.64 10 −3 Abackward4 3.54   1 0 - 3 _ 3.0310-5 2.35 10 −5 4.6410-5 Aforward4 1.81   1 0 - 2 _ 1.4610-3 1.5010-3 7.14 10 −4For the interested reader and for a better visualization of the results that can be obtained via kriging metamodeling, AppendixB gives the best predicted maps of mean and variance of the maximal vibration amplitudes associated with the four backward and forward critical speeds.Finally, it can be concluded that the best kriging metamodeling (i.e., the best prediction based on the relative generalization errorerrorG; see the italic bold number for each amplitude in Tables 10 and 11) gives a very reliable approximation of the reference solution by performing a minimum number of calculations.


## 3.1. Some Details on the Rotor Problem

First of all, Figure2 highlights the Campbell diagram that represents the evolutions of the forcing frequencies as a function of the rotor speed for two deterministic cases (with kvb,2=khb,2=0.1106 N·m−1 for the first case and kvb,2=khb,2=2106 N·m−1 for the second case). The values of the physical parameters used correspond to those previously indicated in Section 2. It is recalled that the coincidence of an excitation frequency with a resonance frequency of the rotor system corresponds to the critical speeds. In other words, this indicates the rotational speeds of the rotor for which large synchronous responses occur. In the present study, four backward and forward critical speeds are present in the rotor speed range of interest. We recall that the terms “backward” and “forward” refer to the whirling motion of the rotor that can be in the same direction as the shafts rotation (“forward whirl”) or in the opposite direction (“backward whirl”). The ith backward and forward critical speeds (i.e., the backward whirl and forward whirl modes) diverge as the rotor speed increases (see Figure 2).Figure 2
Campbell diagram: continuous blue lines forkvb,2=khb,2=0.1106 N⋅m−1 and black dotted lines for kvb,2=khb,2=2106 N⋅m−1.In the following sections, we investigate the potential of kriging based models to be efficiently considered as representative descriptions of the dynamic behavior of the rotor system. For this, we will seek to predict the four backward and forward critical speeds as well as the associated vibration amplitudes according to two deterministic parameters (i.e., the vertical and horizontal bearing support of the second supportkvb,2 and khb,2) and in the presence of uncertainties for seven other selected parameters: the independent thicknesses of the four disks, Young’s modulus of the rotor shaft, and the vertical and horizontal bearing support of the first support. Variability intervals of these seven input parameters are indicated in Table 3. For the proposed study, the deterministic parameters kvb,2 and khb,2 will be set independently for fixed values ranging from 0.1106 N·m−1 to 2106 N·m−1. For each pair of these fixed parameters kvb,2;khb,2, uncertainties according to the seven other parameters will be taken into account.Table 3
Variability intervals of input parameters.
Parameters Variability intervals Thickness of the 1st disk (m) [ . 9 ; 1.1 ] h 1 d Thickness of the 2nd disk (m) [ . 9 ; 1.1 ] h 2 d Thickness of the 3rd disk (m) [ . 9 ; 1.1 ] h 3 d Thickness of the 4th disk (m) [ . 9 ; 1.1 ] h 4 d Young’s modulus of the rotor shaft (Pa or N⋅m−2) [ . 95 ; 1.05 ] E s Vertical bearing stiffness of the first support (N⋅m−1) [ . 95 ; 1.05 ] k v b , 1 Horizontal bearing stiffness of the first support (N⋅m−1) [ . 95 ; 1.05 ] k h b , 1For the reader’s understanding, the explicit functions of theith backward and forward critical speeds relative to the two deterministic parameters kvb,2 and khb,2 are defined as(6)CSbackwardifh1,…,4d,fEs,fkvb,1,fkhb,1kvb,2,khb,2CSforwardifh1,…,4d,fEs,fkvb,1,fkhb,1kvb,2,khb,2for i=1,…,4. fh1,…,4d, fEs, fkvb,1, and fkhb,1 correspond to the uniform probability density functions associated with the uncertainties of the seven selected parameters. If the users have information about the most likely values of the uncertain parameters, then these users may specify nonuniform probability density functions with the most likely values.Similarly, the explicit functions of the associated maximal amplitudes for the rotor system (i.e., the maximum amplitude of the rotor system at the backward and forward critical speeds) are defined as(7)Abackwardifh1,…,4d,fEs,fkvb,1,fkhb,1kvb,2,khb,2=maxAh,backwardifh1,…,4d,fEs,fkvb,1,fkhb,1kvb,2,khb,2,Av,backwardifh1,…,4d,fEs,fkvb,1,fkhb,1kvb,2,khb,2Aforwardifh1,…,4d,fEs,fkvb,1,fkhb,1kvb,2,khb,2=maxAh,forwardifh1,…,4d,fEs,fkvb,1,fkhb,1kvb,2,khb,2,Av,forwardifh1,…,4d,fEs,fkvb,1,fkhb,1kvb,2,khb,2where Ah,backwardi and Av,backwardi (Ah,forwardi and Av,forwardi, resp.) refer to the maximum amplitudes on the horizontal and vertical directions for the ith backward critical speeds (the ith forward critical speeds, resp.). In the following, these explicit functions will be noted CSbackwardi, CSforwardi, Abackwardi, and Aforwardi (for i=1,…,4) in order to simplify the notations in the paper.


## 3.2. Kriging Metamodeling

The main idea with kriging based metamodeling is related to the possibility of estimating unknown values of a function by using a small size set of input/output data previously generated from a predefined experiment design [15]. The key principle consists in the exploiting of the spatial correlations between the function values to correct the average behavior of the function described by a regression model. For a more general description of the kriging formalism, the associated metamodeling method will be presented in the sequel by considering a parameter-dependent output y(p), where y may represent one of the maximum amplitudes, the backward or forward critical speeds defined previously while p∈Rl is a parameter vector which will refer to the vertical and horizontal bearing support kvb,2 and khb,2 (hence l=2). Different kriging models dedicated to each one of the outputs are performed. It may be mentioned that y may also represent multiple outputs (e.g., multiple responses of the maximum amplitudes or/and the critical speeds). However, multivariate kriging is more complicated and may result in inferior predictions because additional kriging parameters must be estimated (for more details see [14]).Based on the kriging theory, the parameter-dependent outputy(p) can be metamodeled as follows:(8)yp=hpTβ+Zp.The latter expression is obtained by juxtaposing two mathematical objects. The first one consists in a linear combination of polynomial functions hi (for i=1,…,Q) defining a regression with a predetermined order. The second object designated by Z(·) is a realization of a Gaussian process which is assumed to be with a zero-mean value and a covariance function given by(9)EZx,Zp=σ2Rθ,x,p,where E· denotes the expectation operator.The covariance function (9) is characterized by both σ2 and R which define the process variance and the spatial correlation function, respectively. The latter which is a monotone function within [0,1] measures how any two points x and p in the design space are spatially close to each other. It is defined by(10)Rθ,x,p=∏j=1lRjθj,xj-pj,where θ∈Rl denotes an l-dimensional vector with scaling parameters.The following parameters have to be determined to completely define a kriging predictor: the parameter vectorsβ and θ plus the scalar σ2. The obtaining of those parameters requires some preliminary choices, namely, the order of the regression model, the spatial correlation function, and the experiment design. The optimal choices are not obvious to determine. However, the nature of the function to be approximated in terms of its smoothness and nonlinearity can guide the users. Concerning the spatial correlation functions, numerous types exist. Some of them are given in Table 4.Table 4
Correlation functions.
Correlation model Rj(θ,xj,pj) Exponential e x p ( - θ j | x j - p j | ) Gaussian exp(-θj(xj-pj)2) Linear max{0,1-θj|xj-pj|} Cubic 1 - 3 ξ j 2 + 2 ξ j 3 with ξj=min(1,θj|xj-pj|)As mentioned previously, defining a kriging predictor requiresβ, σ2, and θ to be determined. For this aim, an input/output of the simulation model is required. The latter is generated according to a predefined experiment design such as a Latin Hypercube sampling. The latter is used in this paper. So, by considering the input parameter vector p1,…,pNk and the associated output vector Y=y1,…,yNk, the kriging parameters can be calculated as follows.First, theθ parameter vector is given by the solution of the following optimization problem:(11)minθ=R1/Nkσ2,where R is a (Nk×Nk) correlation matrix with entries defined by Rij=R(θ,p(i),p(j)) and |·| denotes the determinant of a matrix.The optimal regression parameterβ∗ is given by(12)β∗=HTR-1H-1HTR-1Yas the optimal solution to the following classical least square problem:(13)Hβ≃Y,where H is a (Q×Nk) matrix whose entries are given by Hij=hip(j).The variance processσ2 can then be estimated through(14)σ2=1NkY-Hβ∗TR-1Y-Hβ∗.All parameters being determined, the kriging predictor is completely defined and then can be written as follows:(15)Yx=hTxβ∗+rxTR-1Y-Hβ∗,where rx=Rθ,p(1),x,…,Rθ,pNk,xT is the vector of the correlation coefficients between the new input vector x and the Nk old sampled input vectors (p(i))i∈1,…,Nk.


## 3.3. Reference Solutions

First of all, a Monte Carlo simulation was carried out by considering two deterministic inputskvb,2 and khb,2 in the design interval of interest (i.e., these inputs are deterministic because they are assumed to be controlled by the users) and seven random inputs with uniform distributions (i.e., these inputs are not controlled by the users but are given by the environment of the rotor system, as seen in Table 3). For the reference solutions we examine 100×100 combinations of the two deterministic inputs by sampling the two-dimensional space parameter (kvb,2,khb,2); for each combination, 10000 simulations with the seven noise factors are performed (i.e., N=10000 samples for each combination of the two decision variables (kvb,2,khb,2) are calculated). Then for each deterministic couple (kvb,2,khb,2) the mean and the variance of the four backward and forward critical speeds as well as the associated vibration amplitudes are estimated. These moments are determined by using the classical formula for the mean value and variance given by(16)MeanYx≈1N∑k=1NYxkVarianceYx≈1N∑k=1NYxk-MeanYx2,where N is the number of samples for each combination of the two decision variables kvb,2 and khb,2.All these results are given in Figures3, 4, 5, and 6, respectively. Several comments may be made. First of all the numerical results show that both the values of critical speeds and the associated maximum magnitudes of rotor response are drastically affected by the deterministic couple (kvb,2,khb,2). For example, the mean values of the first and second backward critical speeds changed from 16 Hz to 33 Hz and from 44 Hz to 65 Hz, respectively. Similarly, the mean values of the associated maximum magnitudes vary considerably for all critical speeds as illustrated in Figure 5. The map shapes concerning the mean values and variances of the critical speeds are quite simple. A monotonic increasing of the mean values of each critical speed with the increase of the physical parameters kvb,2 and khb,2 is observed (see Figures 3(a)–3(h)). Evolution of the variance of each critical speed can also be simple with an increase (see Figures 4(a), 4(b), and 4(f)) or decrease (see Figures 4(g) and 4(h)) according to the evolution of kvb,2 or khb,2 or more complex as illustrated in Figures 4(c), 4(d), and 4(e).Figure 3
Referential map of the four backward and forward critical speeds: mean. (a)CSbackward1; (b) CSforward1; (c) CSbackward2; (d) CSforward2; (e) CSbackward3; (f) CSforward3; (g) CSbackward4; and (h) CSforward4.
(a) (b) (c) (d) (e) (f) (g) (h)Figure 4
Referential map of the four backward and forward critical speeds: variance. (a)CSbackward1; (b) CSforward1; (c) CSbackward2; (d) CSforward2; (e) CSbackward3; (f) CSforward3; (g) CSbackward4; and (h) CSforward4.
(a) (b) (c) (d) (e) (f) (g) (h)Figure 5
Referential map of the vibration amplitudes associated with the four backward and forward critical speeds: mean. (a)Abackward1; (b) Aforward1; (c) Abackward2; (d) Aforward2; (e) Abackward3; (f) Aforward3; (g) Abackward4; and (h) Aforward4.
(a) (b) (c) (d) (e) (f) (g) (h)Figure 6
Referential map of the vibration amplitudes associated with the four backward and forward critical speeds: variance. (a)Abackward1; (b) Aforward1; (c) Abackward2; (d) Aforward2; (e) Abackward3; (f) Aforward3; (g) Abackward4; and (h) Aforward4.
(a) (b) (c) (d) (e) (f) (g) (h)Regarding the evolution of amplitudes, the observed behavior is generally more complex with an increase or decrease of the mean and variance for each critical speed with increasing of each parameterkvb,2 or khb,2 (see Figures 5 and 6). For certain specific cases, the asymmetry of the bearing support (i.e., kvb,2≠khb,2) introduces further modifications to the vibration of the rotor. The mean values of the magnitude responses of the asymmetric rotor (i.e., kvb,2≠khb,2) at each backward critical speed are more important than the mean value of the vibration of the symmetric rotor (i.e., kvb,2=khb,2). These effects of the asymmetric support are highly visible in Figures 5(a), 5(c), 5(d), and 5(f). Finally, it is also observed that the evolution of the mean amplitudes of the asymmetric rotor can be difficult to predict intuitively (see, e.g., the evolutions of the maximum amplitudes for the second forward critical speeds in Figure 5(d)). Furthermore, the uncertainties introduced induce a nonnegligible variability on the maximum magnitudes of rotor response at each critical speed, as indicated in Figure 6. More or less complex evolutions of the variability on the maximum magnitudes are still observed.In conclusion, these results will serve as reference solutions for the validation of the kriging surrogate model. The objective being to be able to reproduce all the previously described evolutions concerning the mean and the variance of the four backward and forward critical speeds as well as the associated vibration amplitudes.As mentioned previously, reference solutions are performed by considering100×100 combinations of these two deterministic inputs kvb,2 and khb,2) with 10000 samples for each combination. This leads to a total of 108 calculations for the reference solutions. It can be noted that this number of calculations is prohibitive from a practical point of view. So the basic idea is to be able to drastically reduce this number of samples of the two controllable deterministic inputs kvb,2 and khb,2 by creating continuous surfaces from a predefined set of points. In other words, we will seek to predict the value at the point where the value is unknown by using an interpolation method (and by considering the values at a predefined set of sample points). In the following we will propose using a kriging method that allows the prediction of a variable of interest at a location where it has not been calculated, by taking some form of weighted average of the values at surrounding points.In conclusion, these results will serve as reference solutions for the validation of the kriging surrogate model. The objective being to be able to reproduce all the previously described evolutions concerning the mean and the variance of the four backward and forward critical speeds as well as the associated vibration amplitudes.


## 3.4. Kriging Based Predictions

In this subsection kriging models are used to estimate the referential maps of the four backward and forward critical speeds (shown in Figures3 and 4 for the reference model). In this perspective, the DACE software presented in [21] is exploited. Otherwise, there exist numerous other software packages that may be used for the same aim. Results are not necessarily the same. A comparative study of several software packages defined for kriging metamodeling can be found in [22].As mentioned previously, the construction of a kriging predictor requires preliminary choices of the regression order, the spatial correlation function, and the size of the experiment design. Hence, a study of the influence of these parameters on the accuracy of kriging based predictions of critical speedsCSbackwardi and CSforwardi and the vibration amplitudes Abackwardi and Aforwardi is performed. In a first step, one random experiment design with Nk=600 is fixed and the influence of the regression order and the spatial correlation function is analyzed. In a second step, the regression order and the spatial correlation function are fixed while the size of the experiment design is varied to observe the influence of the size of the experiment design.


### 3.4.1. Prediction of Critical SpeedsCSbackwardi and CSforwardi

With a random plan of sizeNk=600, the mean values and variances of the critical speeds CSbackwardi and CSforwardi and the vibration amplitudes Abackwardi and Aforwardi are estimated against (kvb,2,khb,2) through kriging metamodels considering zero-, first-, and second-order regressions and the linear, Gaussian, exponential, and cubic correlation functions.The maximum and mean of the relative errors which are defined byerrorRMax and errorRMean (for both the mean and variance of each backward or forward critical speed) between the kriging metamodels and the reference model are given by(17)errorRMax=maxXk,ref-Xk,krigXk,refk=1,…,NerrorRMean=1N∑k=1NXk,ref-Xk,krigXk,ref,where Xref corresponds to the previous reference results. Xkrig are the results of the metamodel prediction value. N defines the number of points in the design space (here N=100∗100 as previously indicated in Section 3.3). It may be worth mentioning that each previous criterion errorRMax or errorRMean returns one unique value each for all the design space.Finally the accuracy of the metamodels is compared by using relative generalization error [23], which is defined as (18)errorG=1/N∑k=1NXk,ref-Xk,krig2varXref,where var denotes the variance. As for the two previous criteria, the relative generalization error errorG returns one unique value in the design space according to the variable of interest (i.e., one is interested on the calculation of the mean or variance of a given critical speed in the present study). In several cases, the size Nk of the design experiment is determined by 10 times the number of the considered parameters [24]. In this study, Nk was fixed in such a way that the mean relative errors corresponding to the predicted critical speeds and amplitudes do not exceed a predefined error level (1% in our study).Results of the two criteriaerrorRMax and errorRMean for the mean and variance on the design space are given in Tables 5 and 6 for all the estimated critical speeds, respectively. Tables 7 and 8 give the associated relative generalization error. The values marked in bold (in underlines, resp.) correspond to the minimum (the maximum, resp.) of the error values for each criterion (for a given critical speed and a given regression order) in order to visualize the best (the worst, resp.) spatial correlation functions to predict the desired quantities with kriging models. The values marked in italic bold correspond to the best metamodel. So it indicates the best spatial correlation functions and regression order in order to predict the mean and variance of each critical speed. In this part of the paper illustrations of the maps, which are obtained by using the best kriging metamodels, are not given because they are very similar to the results displayed previously for the reference calculations.Table 5
% of the relative errors(errorRMax;errorRMean) of the mean of critical speeds for different regression orders and correlation functions. Nonsymmetric case with 600 samples.
Gaussian Linear Exponential Cubic Zero-order regression CSbackward1 (17.86; 0.29) (28.58; 0.65) (25.30; 0.36) (40.20; 0.81) CSforward1 (10.94;0.10) (10.33;0.14) (12.82; 0.11) (8.58; 0.10) CSbackward2 (2.55; 0.02) (4.54; 0.16) (3.02; 0.08) (10.40; 0.27) CSforward2 (3.37;0.02) (4.07; 0.08) (0.99; 0.03) (11.01; 0.14) CSbackward3 (0.00; 0.00) (2.04; 0.06) (0.66; 0.02) (4.65; 0.11) CSforward3 (0.00; 0.00) (2.22; 0.05) (0.74; 0.02) (5.15; 0.10) CSbackward4 (0.00; 0.00) (0.54; 0.01) (0.18; 0.00) (1.23; 0.03) CSforward4 (0.00; 0.00) (0.31; 0.01) (0.10; 0.00) (0.70; 0.02)  First-order regression CSbackward1 (11.98; 0.25) (23.81; 0.50) (23.13; 0.30) (27.46; 0.56) CSforward1 (6.51; 0.09) (12.64;0.15) (14.28; 0.12) (7.42;0.08) CSbackward2 (1.54; 0.02) (1.93; 0.07) (2.04; 0.04) (3.66; 0.10) CSforward2 (0.86; 0.01) (1.33;0.04) (1.59; 0.03) (1.69; 0.03) CSbackward3 (0.00; 0.00) (0.39; 0.01) (0.27; 0.00) (0.70; 0.01) CSforward3 (0.00; 0.00) (0.11; 0.00) (0.07; 0.00) (0.26; 0.01) CSbackward4 (0.00; 0.00) (0.05; 0.00) (0.04; 0.00) (0.10; 0.00) CSforward4 (0.00; 0.00) (0.01; 0.00) (0.00; 0.00) (0.01; 0.00)  Second-order regression CSbackward1 (8.98; 0.20) (20.97;0.38) (21.40; 0.24) (16.80; 0.35) CSforward1 (6.35; 0.08) (12.26;0.14) (13.91; 0.12) (8.45; 0.08) CSbackward2 (0.63; 0.01) (1.22;0.04) (1.57; 0.03) (1.11; 0.03) CSforward2 (0.60; 0.01) (1.20;0.03) (1.43; 0.02) (1.14; 0.02) CSbackward3 (0.00; 0.00) (0.07; 0.00) (0.06; 0.00) (0.11; 0.00) CSforward3 (0.00; 0.00) (0.05; 0.00) (0.04; 0.00) (0.10; 0.00) CSbackward4 (0.00; 0.00) (0.00; 0.00) (0.00; 0.00) (0.00; 0.00) CSforward4 (0.00; 0.00) (0.00; 0.00) (0.00; 0.00) (0.00; 0.00)Table 6
% of the relative errors(errorRMax,errorRMean) on the variance of the critical speeds for different regression orders and correlation functions. Nonsymmetric case with 600 samples.
Gaussian Linear Exponential Cubic Zero-order regression CSbackward1 (151.85; 1.18) (76.12; 1.66) (58.14; 0.92) (166.04; 2.84) CSforward1 (48.74;0.78) (21.74; 0.67) (12.19; 0.39) (77.41; 0.72) CSbackward2 (10.58; 0.19) (11.05;0.33) (11.89; 0.22) (7.44; 0.17) CSforward2 (12.40; 0.23) (7.52;0.41) (8.02; 0.27) (6.32; 0.21) CSbackward3 (0.00; 0.00) (0.67; 0.02) (0.29; 0.01) (1.29; 0.04) CSforward3 (0.00; 0.00) (2.29; 0.07) (0.78; 0.02) (5.00; 0.14) CSbackward4 (0.00; 0.00) (1.01; 0.03) (0.69; 0.01) (2.27; 0.05) CSforward4 (0.00; 0.00) (0.49; 0.01) (0.16; 0.00) (1.11; 0.03)  First-order regression CSbackward1 (72.82; 0.72) (46.36; 0.86) (46.19; 0.55) (61.98;1.10) CSforward1 (39.81; 5.32) (13.20; 0.52) (18.12; 0.41) (12.17; 0.30) CSbackward2 (8.74; 0.20) (11.41;0.33) (12.06; 0.23) (7.91; 0.17) CSforward2 (10.79; 0.18) (9.42;0.40) (9.06; 0.26) (9.33;0.17) CSbackward3 (0.00; 0.00) (0.72; 0.01) (0.60; 0.01) (1.23; 0.02) CSforward3 (0.00; 0.00) (0.78; 0.01) (0.68; 0.01) (1.31; 0.01) CSbackward4 (0.00; 0.00) (0.34; 0.01) (0.25; 0.00) (0.60; 0.01) CSforward4 (0.00; 0.00) (0.02; 0.00) (0.01; 0.00) (0.03; 0.00)  Second-order regression CSbackward1 (15.47; 0.29) (35.98;0.50) (39.75; 0.38) (27.97; 0.33) CSforward1 (148.91; 4.25) (12.84; 0.42) (16.13; 0.37) (16.10;0.25) CSbackward2 (6.95; 0.15) (11.16; 0.32) (11.05; 0.21) (8.06;0.13) CSforward2 (7.37; 0.15) (7.06;0.31) (7.70; 0.23) (6.55; 0.13) CSbackward3 (0.00; 0.00) (0.19; 0.01) (0.16; 0.00) (0.33; 0.01) CSforward3 (0.00; 0.00) (0.08; 0.00) (0.06; 0.00) (0.13; 0.00) CSbackward4 (0.00; 0.00) (0.03; 0.00) (0.03; 0.00) (0.05; 0.00) CSforward4 (0.00; 0.00) (0.00; 0.00) (0.00; 0.00) (0.01; 0.00)Table 7
Relative generalization errorerrorG on the mean of the critical speeds for different regression orders and correlation functions. Nonsymmetric case with 600 samples.
Gaussian Linear Exponential Cubic Zero-order regression CSbackward1 1.87 10 −3 6.3710-3 3.1810-3 1.49   10 - 2 _ CSforward1 5.4710-3 5.3610-3 7.30   10 - 3 _ 4.09 10 −3 CSbackward2 6.14 10 −5 8.2310-4 2.8010-4 4.44   10 - 3 _ CSforward2 2.0610-4 5.0910-4 6.69 10 −5 3.43   10 - 3 _ CSbackward3 5.36 10 −12 5.8710-4 1.0710-4 3.54   10 - 3 _ CSforward3 7.71 10 −12 6.0310-4 9.83 10−5 3.36   10 - 3 _ CSbackward4 1.36 10 −10 5.8310-4 9.3310-5 3.29   10 - 3 _ CSforward4 3.90 10 −10 5.7510-4 8.8610-5 3.37   10 - 3 _  First-order regression CSbackward1 1.38 10 −3 3.8510-3 2.4010-3 6.79   10 - 3 _ CSforward1 3.2410-3 7.4910-3 9.54   10 - 3 _ 2.08 10 −3 CSbackward2 2.82 10 −5 1.4610-4 8.7010-5 5.20   10 - 4 _ CSforward2 2.86 10 −5 1.1210-4 9.5110-5 2.06   10 - 4 _ CSbackward3 8.26 10 −13 1.1610-5 3.8910-6 4.89   10 - 5 _ CSforward3 2.09 10 −13 2.7810-6 8.7610-7 1.48   10 - 5 _ CSbackward4 1.03 10 −10 3.1310-6 1.2110-6 1.16   10 - 5 _ CSforward4 3.51 10 −10 9.4010-8 3.7710-8 3.26   10 - 7 _  Second-order regression CSbackward1 8.93 10 −4 2.2510-3 1.7010-3 2.45   10 - 3 _ CSforward1 2.19 10 −3 6.6010-3 9.14   10 - 3 _ 3.2610-3 CSbackward2 8.02 10 −6 4.86   10 - 5 _ 4.6410-5 3.35 10−5 CSforward2 1.15 10 −5 7.26   10 - 5 _ 6.2710-5 5.2010-5 CSbackward3 2.51 10 −13 3.0210-7 1.9310-7 6.46   10 - 7 _ CSforward3 1.74 10 −13 4.4910-7 2.1310-7 1.80   10 - 6 _ CSbackward4 1.23 10 −10 4.6910-10 3.7610-10 9.57   10 - 10 _ CSforward4 3.37 10 −10 6.9610-9 2.6310-9 2.52   10 - 8 _Table 8
Relative generalization errorerrorG on the variance of the critical speeds for different regression orders and correlation functions. Nonsymmetric case with 600 samples.
Gaussian Linear Exponential Cubic Zero-order regression CSbackward1 1.2110-3 1.0810-3 3.97 10 −4 4.61   10 - 3 _ CSforward1 4.07   10 - 3 _ 1.5110-3 7.07 10 −4 3.8110-3 CSbackward2 9.9010-3 1.29   10 - 2 _ 1.0810-2 4.26 10 −3 CSforward2 7.00   10 - 3 _ 5.9710-3 4.1410-3 2.99 10 −3 CSbackward3 7.91 10 −12 2.8310-4 5.4610-5 1.56   10 - 3 _ CSforward3 4.16 10 −12 4.5310-4 7.6310-5 2.76   10 - 3 _ CSbackward4 6.50 10 −10 7.0710-4 1.7110-4 3.34   10 - 3 _ CSforward4 3.08 10 −9 5.7010-4 8.7910-5 3.30   10 - 3 _  First-order regression CSbackward1 3.7510-4 2.9010-4 1.71 10 −4 5.92   10 - 4 _ CSforward1 6.49   10 - 2 _ 9.6310-4 7.6310-4 5.33 10 −4 CSbackward2 8.1610-3 1.33   10 - 2 _ 1.1010-2 4.93 10 −3 CSforward2 4.1710-3 6.19   10 - 3 _ 4.3310-3 2.41 10 −3 CSbackward3 1.81 10 −11 1.1510-4 6.7010-5 3.86   10 - 4 _ CSforward3 5.73 10 −12 1.6210-5 1.0110-5 4.84   10 - 5 _ CSbackward4 7.46 10 −10 3.3310-5 1.1910-5 1.33   10 - 4 _ CSforward4 3.28 10 −9 4.2510-7 1.4510-7 1.94   10 - 6 _  Second-order regression CSbackward1 6.40 10 −5 1.31   10 - 4 _ 1.1310-4 6.9210-5 CSforward1 3.46   10 - 2 _ 7.0610-4 6.7210-4 3.25 10 −4 CSbackward2 4.8910-3 1.29   10 - 2 _ 9.4310-3 3.54 10 −3 CSforward2 2.5010-3 4.31   10 - 3 _ 3.5410-3 1.49 10 −3 CSbackward3 7.46 10 −12 1.5510-5 8.4310-6 5.14   10 - 5 _ CSforward3 9.90 10 −13 6.3610-7 3.5010-7 1.74   10 - 6 _ CSbackward4 7.25 10 −10 1.3410-7 1.1210-7 2.63   10 - 7 _ CSforward4 3.15 10 −9 8.3010-8 3.5310-8 3.05   10 - 7 _First, if some marginal cases are excluded, it can also be noted that, whatever the selected kriging patterns, satisfactory results are generally observed for all the three criteriaerrorRMax, errorRMean, and errorG. It also appears that the best kriging model (in terms of the prediction of the mean and variance of the critical speeds) is not always associated with the same regression order or the same correlation functions. These results illustrate the difficulty of prior choices of both the regression order and the spatial correlation function for a mechanical system. However, it is shown that the best kriging model is generally obtained with the Gaussian correlation function and the second-order regression. In practice the users may try different kriging metamodels (if these metamodels require much less computer time than the simulation code) and select the metamodel that gives best results for a separate (small) test set of inputs.Similarly it appears that, according to the criterion chosen (errorRMax, errorRMean, or errorG), the best kriging model identified for a given output data to be predicted (i.e., a specific critical speed and according to the choice between its means or variances) is not necessarily the same. Moreover, it is also clearly shown that all the values of the maximum and mean of errorRMax, errorRMean, or errorG are not necessarily of the same order of magnitude according to the critical speeds considered. For example, a maximum relative error errorRMax of 8.63% is registered for the mean value of CSbackward1 whereas this value is inferior to 1% for CSbackward2 and CSforward2 and close to 1% for CSbackward3, CSforward3, CSbackward4, and CSforward4. Similarly the values of errorG for both the mean and variance (indicated in Tables 7 and 8) for the first four critical speeds CSbackward1, CSforward1, CSbackward2, and CSforward2 are very different (i.e., superior by a factor of at least 105) from those of the last four critical speeds CSbackward3, CSforward3, CSbackward4, and CSforward4. All these facts show that the definition of relevant criteria for a given mechanical system and the selection of the output data to be predicted play an important role in the construction of a kriging model.In conclusion the efficiency of the proposed kriging approach for the estimation of the four backward and forward critical speeds of the rotor system in presence of uncertainties is validated through the use of three different criteria.


### 3.4.2. Discussion on Exploiting the Symmetry Property of the Problem Addressed

In the previous section, the number of samples used for the calculation of all the kriging predictors has been chosen as equal toNk=600. It may be noted that this number is not optimal due to the fact that Nk is not necessarily the same to ensure the same accuracy for the predictions of the mean values and variances of the critical speeds. In the previous section, this number has been chosen to have an error errorRMean of less than 1% for both the mean and the variance of all the critical speeds of interest.Even if a convergence study with respect toNk for each function needed to be approximated is not the objective of this study (which, for recall, aims to show the possibility of exploiting kriging metamodeling in the analysis of rotor dynamics), it can be said that the necessary number of samples will be smaller if the map to be approximated by the kriging model will have a simple evolution in terms of the parameters considered. This requires an a priori knowledge of the more or less complex evolutions of the quantities processed.So, in the following, a specific study on the influence of the number of samples will be conducted in conjunction with a specific property of the system to drastically reduce the number of samples to be considered. Indeed it can be noted that the size of the experiment design can be drastically reduced by exploiting the physics of the problem. More particularly, exploiting the property of symmetry of the rotor system (between the horizontal and vertical directionskhb,2 and kvb,2) may allow reducing the size of the experiment design, and the input combinations (design points) within that design, without decreasing the accuracy of the corresponding kriging predictions. First, the main idea is to divide the design space according to the identified symmetry property and then to generate an experiment design within the reduced design space. Kriging metamodels are then used to mimic the functions within the new design spaces. To illustrate this on the proposed rotor system, the previous design space (kvb,2,khb,2) with 0.1106 N·m−1≤kvb,2≤2106 N·m−1 and 0.1106 N·m−1≤khb,2≤2106 N·m−1 will be restricted to the following design space: 0.1106 N·m−1 ≤ kvb,2≤khb,2 and kvb,2≤khb,2≤2106 N·m−1. More precisely three cases will be investigated.Case 1.
A total of 350 points are distributed as follows: 50 samples are distributed uniformly (between0.1106 N·m−1 and 2106 N·m−1) at the boundary of the reduced design space (for khb,2=kvb,2) and 300 samples are chosen according to a LHS process restricted to the following design space: 0.1106 N·m−1≤kvb,2<khb,2 and kvb,2<khb,2≤2106 N·m−1.Case 2.
A total of 120 points are distributed as follows: 20 samples are distributed uniformly (between0.1106 N·m−1 and 2106 N·m−1) at the boundary of the reduced design space (for khb,2=kvb,2) and 100 samples are chosen according to a LHS process restricted to the following design space: 0.1106 N·m−1≤kvb,2<khb,2 and kvb,2<khb,2≤2106 N·m−1.Case 3.
A total of 70 points are distributed as follows: 20 samples are distributed uniformly (between0.1106 N·m−1 and 2106 N·m−1) at the boundary of the reduced design space (for khb,2=kvb,2) and 50 samples are chosen according to a LHS process restricted to the following design space: 0.1106 N·m−1≤kvb,2<khb,2 and kvb,2<khb,2≤2106 N·m−1.These three cases of study will enable us to test the performance associated with taking into account the property of symmetry and also to examine the convergence according to number of samples. As previously explained a convergence study with respect toNk for each of the critical speeds has to be performed in order to optimize the number of points required for accurate approximation. However the objective of the present part of the study is only to illustrate the potential of including thoughtful choices in the kriging simulation.It turns that exploiting the symmetry property enhances the accuracies of kriging predictions of the mean values and the variances of the critical speeds observed previously withNk=600. For example, the maximum errors of the three criteria errorRMax, errorRMean, and errorG by considering the best kriging models are given in Table 9 for the previous case (Section 3.4) and the three cases dealt with in this part of the paper. It is shown that the best kriging metamodel with only 120 samples and by using the symmetry property gives better results for most of the tested criteria than the best kriging metamodel without considering the symmetry property. Thus it is possible to considerably reduce the number of samples required to have predictive models if one is able to employ certain characteristics directly deriving from the physics of the mechanical system studied.Table 9
Maximum errors for all the critical speeds by considering the best kriging models.
e r r o r R M a x ; e r r o r R M e a n e r r o r G Mean values Case 1: 350 samples (1.20; 0.02) 4.8410-5 Case 2: 120 samples (1.99; 0.06) 2.1210-4 Case 3: 70 samples (18.61; 0.31) 3.6310-3  Variance values Case 1: 350 samples (4.80; 1.55) 6.5710-2 Case 2: 120 samples (12.27; 1.65) 6.6910-2 Case 3: 70 samples (30.99; 1.92) 6.6010-2Furthermore, by comparing the results of the three cases, it is clearly shown that increasing the number of samples significantly improves predictions of kriging models. To illustrate this, Table9 summarizes the maximum value of each criterion of errorRMax, errorRMean, and errorG for the prediction of the mean and variance of all the critical speeds.For the interested reader, AppendixA gives all the calculations of the three criteria errorRMax, errorRMean, and errorG for different regression orders and correlation functions. A detailed analysis of all these results leads to the following conclusions:(i)
Kriging metamodeling is a valuable tool to provide valuable approximate results in the field of rotordynamics. It is very interesting for practical mechanical problems where limited outputs data is available (due to expensive simulation codes).(ii)
The best kriging metamodel, for the prediction of the mean or the variance of a given critical speed, does not necessarily correspond to a specific correlation function or regression order. However, for the problem under study, it may be noted that the Gaussian correlation function and the second-order regression often correspond to the best compromise to obtain the most accurate kriging model versus all the criteria tested.(iii)
The necessary number of samples to reach a minimum according to a given criterion is not identical for each critical speed. For example, considering 70 samples is sufficient to have the errorerrorRMax less than 1% for the mean of CSbackward3, CSforward3, CSbackward4, and CSforward4 whereas 120 samples (350 samples, resp.) are needed for CSbackward2 and CSforward2 (CSforward1, resp.); the error errorRMax of the mean of CSbackward1 is still greater than 1% even with 350 samples.In conclusion, even if the choices of the correlation function and the regression order are crucial for an accurate kriging predictor, it is demonstrated that a drastic reduction in the number of samples necessary for a correct prediction of all the critical speeds of the rotor system is possible.


### 3.4.3. Prediction of the Vibration AmplitudesAbackwardi and Aforwardi

The same study is performed for the estimation of the vibration amplitudesAbackwardi and Aforwardi by using kriging metamodels without considering the symmetry property. Indeed, the properties of symmetry on results of the vibratory responses are no longer valid due to the presence of gravity. Moreover, only the generalized error criterion errorG, which is defined in (18), can be used in a robust manner. Indeed, vibratory responses of very small amplitudes for certain configurations (i.e., the choice of parameters kvb,2 and khb,2) can generate very high error levels (the two criteria errorRMax and errorRMean are then maximized because of a division by a quantity close to zero).Results of the relatively generalized errors for the mean and variance on the design space are given in Tables10 and 11 for the prediction of all the vibration amplitudes Abackwardi and Aforwardi (for i=1,…,4). We recall that the values marked in bold (in underlines, resp.) correspond to the best (the worst) metamodel for a given regression order. As said previously, results show that the choice of spatial correlation functions and regression orders can affect efficiency of kriging metamodeling. Similarly, the best spatial correlation function and the best regression order are not necessarily identical for each of the vibration amplitudes (i.e., Abackwardi and Aforwardi for i=1,…,4). Moreover, for a given vibration amplitude, the best spatial correlation function and regression order for the mean value are not necessarily those which will also make it possible to minimize the error of the variance. This shows the difficulty of choosing a priori the best spatial correlation function and the best regression order for the prediction of the maximal amplitudes of a rotor system. However it can be noted that a cubic correlation with a second-order regression works well in most cases.Table 10
Relative generalization errorerrorG on the mean of the amplitude for different regression orders and correlation functions. Nonsymmetric case with 600 samples.
Gaussian Linear Exponential Cubic Zero-order regression Abackward1 9.43   10 - 1 _ 6.5210-2 6.07 10 - 2 8.4310-2 Aforward1 9.45   10 - 1 _ 2.1510-3 1.39 10 - 3 6.3710-3 Abackward2 1.1110-2 1.89   10 - 2 _ 1.1810-2 5.12 10 - 3 Aforward2 1.7810-2 1.79   10 - 2 _ 1.2510-2 7.05 10 - 3 Abackward3 1.0810-3 1.1010-3 5.58 10 - 4 3.63   10 - 3 _ Aforward3 3.45 10 −3 4.69   10 - 3 _ 4.2610-3 4.0510-3 Abackward4 5.3110-4 9.2710-4 4.43 10 - 4 3.29   10 - 3 _ Aforward4 2.4810-3 1.2710-3 1.09 10 - 3 3.60   10 - 3 _  First-order regression Abackward1 7.90   10 - 1 _ 7.3510-2 5.95 10 - 2 7.8610-2 Aforward1 5.29   10 - 2 _ 1.6810-3 1.42 10 - 3 2.8410-3 Abackward2 6.3410-3 1.0110-2 1.01   10 - 2 _ 3.41 10 - 3 Aforward2 3.92   10 - 1 _ 1.2410-2 1.1510-2 6.14 10 - 3 Abackward3 9.7810-4 6.4810-4 4.08 10 - 4 9.81   10 - 4 _ Aforward3 3.0810-3 4.1510-3 4.23   10 - 3 _ 1.50 10 - 3 Abackward4 4.1410-4 2.3410-4 1.12 10 - 4 7.45   10 - 4 _ Aforward4 1.0410-3 1.6010-3 1.69   10 - 3 _ 6.54 10 - 4  Second-order regression Abackward1 5.34   10 - 1 _ 5.9610-2 5.55 10 −2 6.9710-2 Aforward1 1.23   10 - 2 _ 1.4210-3 1.22 10 −3 1.6910-3 Abackward2 1.43   10 - 1 _ 9.7310-3 9.2210-3 3.28 10 −3 Aforward2 2.39   10 - 1 _ 9.3110-3 1.0910-2 5.66 10 −3 Abackward3 2.2510-4 4.02   10 - 4 _ 3.2410-4 1.52 10 −4 Aforward3 1.7810-3 2.6910-3 3.50   10 - 3 _ 1.35 10 −3 Abackward4 9.7610-5 1.36   10 - 4 _ 9.6810-5 7.95 10 −5 Aforward4 6.1010-4 1.0010-3 1.33   10 - 3 _ 5.22 10 −4Table 11
Relative generalization errorerrorG on the variance of the amplitude for different regression orders and correlation functions. Nonsymmetric case with 600 samples.
Gaussian Linear Exponential Cubic Zero-order regression Abackward1 9.42   1 0 - 1 _ 1.8610-2 1.3810-2 9.30 1 0 - 3 Aforward1 9.42   1 0 - 1 _ 4.0910-2 3.51 1 0 - 2 4.1710-2 Abackward2 9.35   1 0 - 1 _ 3.3910-1 1.7910-1 1.53 10 −1 Aforward2 9.44   1 0 - 1 _ 1.4510-1 8.8510-2 6.52 1 0 - 2 Abackward3 3.75 1 0 - 4 1.0110-3 4.8210-4 4.69   1 0 - 3 _ Aforward3 5.2810-3 5.36   1 0 - 3 _ 4.2810-3 3.39 1 0 - 3 Abackward4 1.0910-3 1.7010-3 9.23 1 0 - 4 4.39   1 0 - 3 _ Aforward4 4.58   1 0 - 3 _ 1.4110-3 1.07 1 0 - 3 3.9510-3  First-order regression Abackward1 3.42   1 0 - 1 _ 1.3110-2 1.2810-2 8.56 1 0 - 3 Aforward1 3.18   1 0 - 1 _ 3.2810-2 3.28 1 0 - 2 3.7110-2 Abackward2 7.53   1 0 - 1 _ 3.4510-1 1.8110-1 1.55 1 0 - 1 Aforward2 3.39   1 0 - 1 _ 1.4410-1 7.3910-2 5.28 10 −2 Abackward3 1.2710-4 1.7810-4 5.42 1 0 - 5 1.13   1 0 - 3 _ Aforward3 3.9710-3 5.28   1 0 - 3 _ 4.5010-3 2.02 1 0 - 3 Abackward4 5.0010-4 3.8910-4 1.89 1 0 - 4 1.21   1 0 - 3 _ Aforward4 1.96   1 0 - 3 _ 1.6010-3 1.6710-3 7.52 1 0 - 4  Second-order regression Abackward1 2.77   1 0 - 1 _ 1.2810-2 1.2610-2 8.44 10 −3 Aforward1 2.39   1 0 - 1 _ 3.10 10 −2 3.2110-2 3.5710-2 Abackward2 6.20   1 0 - 1 _ 3.4710-1 1.8810-1 1.67 1 0 - 1 Aforward2 3.08   1 0 - 1 _ 1.4410-1 7.4010-2 5.40 1 0 - 2 Abackward3 5.18   1 0 - 5 2.68 10 −5 2.7910-5 7.82   1 0 - 5 _ Aforward3 1.44   1 0 - 1 _ 3.3310-3 3.7210-3 1.64 10 −3 Abackward4 3.54   1 0 - 3 _ 3.0310-5 2.35 10 −5 4.6410-5 Aforward4 1.81   1 0 - 2 _ 1.4610-3 1.5010-3 7.14 10 −4For the interested reader and for a better visualization of the results that can be obtained via kriging metamodeling, AppendixB gives the best predicted maps of mean and variance of the maximal vibration amplitudes associated with the four backward and forward critical speeds.Finally, it can be concluded that the best kriging metamodeling (i.e., the best prediction based on the relative generalization errorerrorG; see the italic bold number for each amplitude in Tables 10 and 11) gives a very reliable approximation of the reference solution by performing a minimum number of calculations.


## 3.4.1. Prediction of Critical SpeedsCSbackwardi and CSforwardi

With a random plan of sizeNk=600, the mean values and variances of the critical speeds CSbackwardi and CSforwardi and the vibration amplitudes Abackwardi and Aforwardi are estimated against (kvb,2,khb,2) through kriging metamodels considering zero-, first-, and second-order regressions and the linear, Gaussian, exponential, and cubic correlation functions.The maximum and mean of the relative errors which are defined byerrorRMax and errorRMean (for both the mean and variance of each backward or forward critical speed) between the kriging metamodels and the reference model are given by(17)errorRMax=maxXk,ref-Xk,krigXk,refk=1,…,NerrorRMean=1N∑k=1NXk,ref-Xk,krigXk,ref,where Xref corresponds to the previous reference results. Xkrig are the results of the metamodel prediction value. N defines the number of points in the design space (here N=100∗100 as previously indicated in Section 3.3). It may be worth mentioning that each previous criterion errorRMax or errorRMean returns one unique value each for all the design space.Finally the accuracy of the metamodels is compared by using relative generalization error [23], which is defined as (18)errorG=1/N∑k=1NXk,ref-Xk,krig2varXref,where var denotes the variance. As for the two previous criteria, the relative generalization error errorG returns one unique value in the design space according to the variable of interest (i.e., one is interested on the calculation of the mean or variance of a given critical speed in the present study). In several cases, the size Nk of the design experiment is determined by 10 times the number of the considered parameters [24]. In this study, Nk was fixed in such a way that the mean relative errors corresponding to the predicted critical speeds and amplitudes do not exceed a predefined error level (1% in our study).Results of the two criteriaerrorRMax and errorRMean for the mean and variance on the design space are given in Tables 5 and 6 for all the estimated critical speeds, respectively. Tables 7 and 8 give the associated relative generalization error. The values marked in bold (in underlines, resp.) correspond to the minimum (the maximum, resp.) of the error values for each criterion (for a given critical speed and a given regression order) in order to visualize the best (the worst, resp.) spatial correlation functions to predict the desired quantities with kriging models. The values marked in italic bold correspond to the best metamodel. So it indicates the best spatial correlation functions and regression order in order to predict the mean and variance of each critical speed. In this part of the paper illustrations of the maps, which are obtained by using the best kriging metamodels, are not given because they are very similar to the results displayed previously for the reference calculations.Table 5
% of the relative errors(errorRMax;errorRMean) of the mean of critical speeds for different regression orders and correlation functions. Nonsymmetric case with 600 samples.
Gaussian Linear Exponential Cubic Zero-order regression CSbackward1 (17.86; 0.29) (28.58; 0.65) (25.30; 0.36) (40.20; 0.81) CSforward1 (10.94;0.10) (10.33;0.14) (12.82; 0.11) (8.58; 0.10) CSbackward2 (2.55; 0.02) (4.54; 0.16) (3.02; 0.08) (10.40; 0.27) CSforward2 (3.37;0.02) (4.07; 0.08) (0.99; 0.03) (11.01; 0.14) CSbackward3 (0.00; 0.00) (2.04; 0.06) (0.66; 0.02) (4.65; 0.11) CSforward3 (0.00; 0.00) (2.22; 0.05) (0.74; 0.02) (5.15; 0.10) CSbackward4 (0.00; 0.00) (0.54; 0.01) (0.18; 0.00) (1.23; 0.03) CSforward4 (0.00; 0.00) (0.31; 0.01) (0.10; 0.00) (0.70; 0.02)  First-order regression CSbackward1 (11.98; 0.25) (23.81; 0.50) (23.13; 0.30) (27.46; 0.56) CSforward1 (6.51; 0.09) (12.64;0.15) (14.28; 0.12) (7.42;0.08) CSbackward2 (1.54; 0.02) (1.93; 0.07) (2.04; 0.04) (3.66; 0.10) CSforward2 (0.86; 0.01) (1.33;0.04) (1.59; 0.03) (1.69; 0.03) CSbackward3 (0.00; 0.00) (0.39; 0.01) (0.27; 0.00) (0.70; 0.01) CSforward3 (0.00; 0.00) (0.11; 0.00) (0.07; 0.00) (0.26; 0.01) CSbackward4 (0.00; 0.00) (0.05; 0.00) (0.04; 0.00) (0.10; 0.00) CSforward4 (0.00; 0.00) (0.01; 0.00) (0.00; 0.00) (0.01; 0.00)  Second-order regression CSbackward1 (8.98; 0.20) (20.97;0.38) (21.40; 0.24) (16.80; 0.35) CSforward1 (6.35; 0.08) (12.26;0.14) (13.91; 0.12) (8.45; 0.08) CSbackward2 (0.63; 0.01) (1.22;0.04) (1.57; 0.03) (1.11; 0.03) CSforward2 (0.60; 0.01) (1.20;0.03) (1.43; 0.02) (1.14; 0.02) CSbackward3 (0.00; 0.00) (0.07; 0.00) (0.06; 0.00) (0.11; 0.00) CSforward3 (0.00; 0.00) (0.05; 0.00) (0.04; 0.00) (0.10; 0.00) CSbackward4 (0.00; 0.00) (0.00; 0.00) (0.00; 0.00) (0.00; 0.00) CSforward4 (0.00; 0.00) (0.00; 0.00) (0.00; 0.00) (0.00; 0.00)Table 6
% of the relative errors(errorRMax,errorRMean) on the variance of the critical speeds for different regression orders and correlation functions. Nonsymmetric case with 600 samples.
Gaussian Linear Exponential Cubic Zero-order regression CSbackward1 (151.85; 1.18) (76.12; 1.66) (58.14; 0.92) (166.04; 2.84) CSforward1 (48.74;0.78) (21.74; 0.67) (12.19; 0.39) (77.41; 0.72) CSbackward2 (10.58; 0.19) (11.05;0.33) (11.89; 0.22) (7.44; 0.17) CSforward2 (12.40; 0.23) (7.52;0.41) (8.02; 0.27) (6.32; 0.21) CSbackward3 (0.00; 0.00) (0.67; 0.02) (0.29; 0.01) (1.29; 0.04) CSforward3 (0.00; 0.00) (2.29; 0.07) (0.78; 0.02) (5.00; 0.14) CSbackward4 (0.00; 0.00) (1.01; 0.03) (0.69; 0.01) (2.27; 0.05) CSforward4 (0.00; 0.00) (0.49; 0.01) (0.16; 0.00) (1.11; 0.03)  First-order regression CSbackward1 (72.82; 0.72) (46.36; 0.86) (46.19; 0.55) (61.98;1.10) CSforward1 (39.81; 5.32) (13.20; 0.52) (18.12; 0.41) (12.17; 0.30) CSbackward2 (8.74; 0.20) (11.41;0.33) (12.06; 0.23) (7.91; 0.17) CSforward2 (10.79; 0.18) (9.42;0.40) (9.06; 0.26) (9.33;0.17) CSbackward3 (0.00; 0.00) (0.72; 0.01) (0.60; 0.01) (1.23; 0.02) CSforward3 (0.00; 0.00) (0.78; 0.01) (0.68; 0.01) (1.31; 0.01) CSbackward4 (0.00; 0.00) (0.34; 0.01) (0.25; 0.00) (0.60; 0.01) CSforward4 (0.00; 0.00) (0.02; 0.00) (0.01; 0.00) (0.03; 0.00)  Second-order regression CSbackward1 (15.47; 0.29) (35.98;0.50) (39.75; 0.38) (27.97; 0.33) CSforward1 (148.91; 4.25) (12.84; 0.42) (16.13; 0.37) (16.10;0.25) CSbackward2 (6.95; 0.15) (11.16; 0.32) (11.05; 0.21) (8.06;0.13) CSforward2 (7.37; 0.15) (7.06;0.31) (7.70; 0.23) (6.55; 0.13) CSbackward3 (0.00; 0.00) (0.19; 0.01) (0.16; 0.00) (0.33; 0.01) CSforward3 (0.00; 0.00) (0.08; 0.00) (0.06; 0.00) (0.13; 0.00) CSbackward4 (0.00; 0.00) (0.03; 0.00) (0.03; 0.00) (0.05; 0.00) CSforward4 (0.00; 0.00) (0.00; 0.00) (0.00; 0.00) (0.01; 0.00)Table 7
Relative generalization errorerrorG on the mean of the critical speeds for different regression orders and correlation functions. Nonsymmetric case with 600 samples.
Gaussian Linear Exponential Cubic Zero-order regression CSbackward1 1.87 10 −3 6.3710-3 3.1810-3 1.49   10 - 2 _ CSforward1 5.4710-3 5.3610-3 7.30   10 - 3 _ 4.09 10 −3 CSbackward2 6.14 10 −5 8.2310-4 2.8010-4 4.44   10 - 3 _ CSforward2 2.0610-4 5.0910-4 6.69 10 −5 3.43   10 - 3 _ CSbackward3 5.36 10 −12 5.8710-4 1.0710-4 3.54   10 - 3 _ CSforward3 7.71 10 −12 6.0310-4 9.83 10−5 3.36   10 - 3 _ CSbackward4 1.36 10 −10 5.8310-4 9.3310-5 3.29   10 - 3 _ CSforward4 3.90 10 −10 5.7510-4 8.8610-5 3.37   10 - 3 _  First-order regression CSbackward1 1.38 10 −3 3.8510-3 2.4010-3 6.79   10 - 3 _ CSforward1 3.2410-3 7.4910-3 9.54   10 - 3 _ 2.08 10 −3 CSbackward2 2.82 10 −5 1.4610-4 8.7010-5 5.20   10 - 4 _ CSforward2 2.86 10 −5 1.1210-4 9.5110-5 2.06   10 - 4 _ CSbackward3 8.26 10 −13 1.1610-5 3.8910-6 4.89   10 - 5 _ CSforward3 2.09 10 −13 2.7810-6 8.7610-7 1.48   10 - 5 _ CSbackward4 1.03 10 −10 3.1310-6 1.2110-6 1.16   10 - 5 _ CSforward4 3.51 10 −10 9.4010-8 3.7710-8 3.26   10 - 7 _  Second-order regression CSbackward1 8.93 10 −4 2.2510-3 1.7010-3 2.45   10 - 3 _ CSforward1 2.19 10 −3 6.6010-3 9.14   10 - 3 _ 3.2610-3 CSbackward2 8.02 10 −6 4.86   10 - 5 _ 4.6410-5 3.35 10−5 CSforward2 1.15 10 −5 7.26   10 - 5 _ 6.2710-5 5.2010-5 CSbackward3 2.51 10 −13 3.0210-7 1.9310-7 6.46   10 - 7 _ CSforward3 1.74 10 −13 4.4910-7 2.1310-7 1.80   10 - 6 _ CSbackward4 1.23 10 −10 4.6910-10 3.7610-10 9.57   10 - 10 _ CSforward4 3.37 10 −10 6.9610-9 2.6310-9 2.52   10 - 8 _Table 8
Relative generalization errorerrorG on the variance of the critical speeds for different regression orders and correlation functions. Nonsymmetric case with 600 samples.
Gaussian Linear Exponential Cubic Zero-order regression CSbackward1 1.2110-3 1.0810-3 3.97 10 −4 4.61   10 - 3 _ CSforward1 4.07   10 - 3 _ 1.5110-3 7.07 10 −4 3.8110-3 CSbackward2 9.9010-3 1.29   10 - 2 _ 1.0810-2 4.26 10 −3 CSforward2 7.00   10 - 3 _ 5.9710-3 4.1410-3 2.99 10 −3 CSbackward3 7.91 10 −12 2.8310-4 5.4610-5 1.56   10 - 3 _ CSforward3 4.16 10 −12 4.5310-4 7.6310-5 2.76   10 - 3 _ CSbackward4 6.50 10 −10 7.0710-4 1.7110-4 3.34   10 - 3 _ CSforward4 3.08 10 −9 5.7010-4 8.7910-5 3.30   10 - 3 _  First-order regression CSbackward1 3.7510-4 2.9010-4 1.71 10 −4 5.92   10 - 4 _ CSforward1 6.49   10 - 2 _ 9.6310-4 7.6310-4 5.33 10 −4 CSbackward2 8.1610-3 1.33   10 - 2 _ 1.1010-2 4.93 10 −3 CSforward2 4.1710-3 6.19   10 - 3 _ 4.3310-3 2.41 10 −3 CSbackward3 1.81 10 −11 1.1510-4 6.7010-5 3.86   10 - 4 _ CSforward3 5.73 10 −12 1.6210-5 1.0110-5 4.84   10 - 5 _ CSbackward4 7.46 10 −10 3.3310-5 1.1910-5 1.33   10 - 4 _ CSforward4 3.28 10 −9 4.2510-7 1.4510-7 1.94   10 - 6 _  Second-order regression CSbackward1 6.40 10 −5 1.31   10 - 4 _ 1.1310-4 6.9210-5 CSforward1 3.46   10 - 2 _ 7.0610-4 6.7210-4 3.25 10 −4 CSbackward2 4.8910-3 1.29   10 - 2 _ 9.4310-3 3.54 10 −3 CSforward2 2.5010-3 4.31   10 - 3 _ 3.5410-3 1.49 10 −3 CSbackward3 7.46 10 −12 1.5510-5 8.4310-6 5.14   10 - 5 _ CSforward3 9.90 10 −13 6.3610-7 3.5010-7 1.74   10 - 6 _ CSbackward4 7.25 10 −10 1.3410-7 1.1210-7 2.63   10 - 7 _ CSforward4 3.15 10 −9 8.3010-8 3.5310-8 3.05   10 - 7 _First, if some marginal cases are excluded, it can also be noted that, whatever the selected kriging patterns, satisfactory results are generally observed for all the three criteriaerrorRMax, errorRMean, and errorG. It also appears that the best kriging model (in terms of the prediction of the mean and variance of the critical speeds) is not always associated with the same regression order or the same correlation functions. These results illustrate the difficulty of prior choices of both the regression order and the spatial correlation function for a mechanical system. However, it is shown that the best kriging model is generally obtained with the Gaussian correlation function and the second-order regression. In practice the users may try different kriging metamodels (if these metamodels require much less computer time than the simulation code) and select the metamodel that gives best results for a separate (small) test set of inputs.Similarly it appears that, according to the criterion chosen (errorRMax, errorRMean, or errorG), the best kriging model identified for a given output data to be predicted (i.e., a specific critical speed and according to the choice between its means or variances) is not necessarily the same. Moreover, it is also clearly shown that all the values of the maximum and mean of errorRMax, errorRMean, or errorG are not necessarily of the same order of magnitude according to the critical speeds considered. For example, a maximum relative error errorRMax of 8.63% is registered for the mean value of CSbackward1 whereas this value is inferior to 1% for CSbackward2 and CSforward2 and close to 1% for CSbackward3, CSforward3, CSbackward4, and CSforward4. Similarly the values of errorG for both the mean and variance (indicated in Tables 7 and 8) for the first four critical speeds CSbackward1, CSforward1, CSbackward2, and CSforward2 are very different (i.e., superior by a factor of at least 105) from those of the last four critical speeds CSbackward3, CSforward3, CSbackward4, and CSforward4. All these facts show that the definition of relevant criteria for a given mechanical system and the selection of the output data to be predicted play an important role in the construction of a kriging model.In conclusion the efficiency of the proposed kriging approach for the estimation of the four backward and forward critical speeds of the rotor system in presence of uncertainties is validated through the use of three different criteria.


## 3.4.2. Discussion on Exploiting the Symmetry Property of the Problem Addressed

In the previous section, the number of samples used for the calculation of all the kriging predictors has been chosen as equal toNk=600. It may be noted that this number is not optimal due to the fact that Nk is not necessarily the same to ensure the same accuracy for the predictions of the mean values and variances of the critical speeds. In the previous section, this number has been chosen to have an error errorRMean of less than 1% for both the mean and the variance of all the critical speeds of interest.Even if a convergence study with respect toNk for each function needed to be approximated is not the objective of this study (which, for recall, aims to show the possibility of exploiting kriging metamodeling in the analysis of rotor dynamics), it can be said that the necessary number of samples will be smaller if the map to be approximated by the kriging model will have a simple evolution in terms of the parameters considered. This requires an a priori knowledge of the more or less complex evolutions of the quantities processed.So, in the following, a specific study on the influence of the number of samples will be conducted in conjunction with a specific property of the system to drastically reduce the number of samples to be considered. Indeed it can be noted that the size of the experiment design can be drastically reduced by exploiting the physics of the problem. More particularly, exploiting the property of symmetry of the rotor system (between the horizontal and vertical directionskhb,2 and kvb,2) may allow reducing the size of the experiment design, and the input combinations (design points) within that design, without decreasing the accuracy of the corresponding kriging predictions. First, the main idea is to divide the design space according to the identified symmetry property and then to generate an experiment design within the reduced design space. Kriging metamodels are then used to mimic the functions within the new design spaces. To illustrate this on the proposed rotor system, the previous design space (kvb,2,khb,2) with 0.1106 N·m−1≤kvb,2≤2106 N·m−1 and 0.1106 N·m−1≤khb,2≤2106 N·m−1 will be restricted to the following design space: 0.1106 N·m−1 ≤ kvb,2≤khb,2 and kvb,2≤khb,2≤2106 N·m−1. More precisely three cases will be investigated.Case 1.
A total of 350 points are distributed as follows: 50 samples are distributed uniformly (between0.1106 N·m−1 and 2106 N·m−1) at the boundary of the reduced design space (for khb,2=kvb,2) and 300 samples are chosen according to a LHS process restricted to the following design space: 0.1106 N·m−1≤kvb,2<khb,2 and kvb,2<khb,2≤2106 N·m−1.Case 2.
A total of 120 points are distributed as follows: 20 samples are distributed uniformly (between0.1106 N·m−1 and 2106 N·m−1) at the boundary of the reduced design space (for khb,2=kvb,2) and 100 samples are chosen according to a LHS process restricted to the following design space: 0.1106 N·m−1≤kvb,2<khb,2 and kvb,2<khb,2≤2106 N·m−1.Case 3.
A total of 70 points are distributed as follows: 20 samples are distributed uniformly (between0.1106 N·m−1 and 2106 N·m−1) at the boundary of the reduced design space (for khb,2=kvb,2) and 50 samples are chosen according to a LHS process restricted to the following design space: 0.1106 N·m−1≤kvb,2<khb,2 and kvb,2<khb,2≤2106 N·m−1.These three cases of study will enable us to test the performance associated with taking into account the property of symmetry and also to examine the convergence according to number of samples. As previously explained a convergence study with respect toNk for each of the critical speeds has to be performed in order to optimize the number of points required for accurate approximation. However the objective of the present part of the study is only to illustrate the potential of including thoughtful choices in the kriging simulation.It turns that exploiting the symmetry property enhances the accuracies of kriging predictions of the mean values and the variances of the critical speeds observed previously withNk=600. For example, the maximum errors of the three criteria errorRMax, errorRMean, and errorG by considering the best kriging models are given in Table 9 for the previous case (Section 3.4) and the three cases dealt with in this part of the paper. It is shown that the best kriging metamodel with only 120 samples and by using the symmetry property gives better results for most of the tested criteria than the best kriging metamodel without considering the symmetry property. Thus it is possible to considerably reduce the number of samples required to have predictive models if one is able to employ certain characteristics directly deriving from the physics of the mechanical system studied.Table 9
Maximum errors for all the critical speeds by considering the best kriging models.
e r r o r R M a x ; e r r o r R M e a n e r r o r G Mean values Case 1: 350 samples (1.20; 0.02) 4.8410-5 Case 2: 120 samples (1.99; 0.06) 2.1210-4 Case 3: 70 samples (18.61; 0.31) 3.6310-3  Variance values Case 1: 350 samples (4.80; 1.55) 6.5710-2 Case 2: 120 samples (12.27; 1.65) 6.6910-2 Case 3: 70 samples (30.99; 1.92) 6.6010-2Furthermore, by comparing the results of the three cases, it is clearly shown that increasing the number of samples significantly improves predictions of kriging models. To illustrate this, Table9 summarizes the maximum value of each criterion of errorRMax, errorRMean, and errorG for the prediction of the mean and variance of all the critical speeds.For the interested reader, AppendixA gives all the calculations of the three criteria errorRMax, errorRMean, and errorG for different regression orders and correlation functions. A detailed analysis of all these results leads to the following conclusions:(i)
Kriging metamodeling is a valuable tool to provide valuable approximate results in the field of rotordynamics. It is very interesting for practical mechanical problems where limited outputs data is available (due to expensive simulation codes).(ii)
The best kriging metamodel, for the prediction of the mean or the variance of a given critical speed, does not necessarily correspond to a specific correlation function or regression order. However, for the problem under study, it may be noted that the Gaussian correlation function and the second-order regression often correspond to the best compromise to obtain the most accurate kriging model versus all the criteria tested.(iii)
The necessary number of samples to reach a minimum according to a given criterion is not identical for each critical speed. For example, considering 70 samples is sufficient to have the errorerrorRMax less than 1% for the mean of CSbackward3, CSforward3, CSbackward4, and CSforward4 whereas 120 samples (350 samples, resp.) are needed for CSbackward2 and CSforward2 (CSforward1, resp.); the error errorRMax of the mean of CSbackward1 is still greater than 1% even with 350 samples.In conclusion, even if the choices of the correlation function and the regression order are crucial for an accurate kriging predictor, it is demonstrated that a drastic reduction in the number of samples necessary for a correct prediction of all the critical speeds of the rotor system is possible.


## 3.4.3. Prediction of the Vibration AmplitudesAbackwardi and Aforwardi

The same study is performed for the estimation of the vibration amplitudesAbackwardi and Aforwardi by using kriging metamodels without considering the symmetry property. Indeed, the properties of symmetry on results of the vibratory responses are no longer valid due to the presence of gravity. Moreover, only the generalized error criterion errorG, which is defined in (18), can be used in a robust manner. Indeed, vibratory responses of very small amplitudes for certain configurations (i.e., the choice of parameters kvb,2 and khb,2) can generate very high error levels (the two criteria errorRMax and errorRMean are then maximized because of a division by a quantity close to zero).Results of the relatively generalized errors for the mean and variance on the design space are given in Tables10 and 11 for the prediction of all the vibration amplitudes Abackwardi and Aforwardi (for i=1,…,4). We recall that the values marked in bold (in underlines, resp.) correspond to the best (the worst) metamodel for a given regression order. As said previously, results show that the choice of spatial correlation functions and regression orders can affect efficiency of kriging metamodeling. Similarly, the best spatial correlation function and the best regression order are not necessarily identical for each of the vibration amplitudes (i.e., Abackwardi and Aforwardi for i=1,…,4). Moreover, for a given vibration amplitude, the best spatial correlation function and regression order for the mean value are not necessarily those which will also make it possible to minimize the error of the variance. This shows the difficulty of choosing a priori the best spatial correlation function and the best regression order for the prediction of the maximal amplitudes of a rotor system. However it can be noted that a cubic correlation with a second-order regression works well in most cases.Table 10
Relative generalization errorerrorG on the mean of the amplitude for different regression orders and correlation functions. Nonsymmetric case with 600 samples.
Gaussian Linear Exponential Cubic Zero-order regression Abackward1 9.43   10 - 1 _ 6.5210-2 6.07 10 - 2 8.4310-2 Aforward1 9.45   10 - 1 _ 2.1510-3 1.39 10 - 3 6.3710-3 Abackward2 1.1110-2 1.89   10 - 2 _ 1.1810-2 5.12 10 - 3 Aforward2 1.7810-2 1.79   10 - 2 _ 1.2510-2 7.05 10 - 3 Abackward3 1.0810-3 1.1010-3 5.58 10 - 4 3.63   10 - 3 _ Aforward3 3.45 10 −3 4.69   10 - 3 _ 4.2610-3 4.0510-3 Abackward4 5.3110-4 9.2710-4 4.43 10 - 4 3.29   10 - 3 _ Aforward4 2.4810-3 1.2710-3 1.09 10 - 3 3.60   10 - 3 _  First-order regression Abackward1 7.90   10 - 1 _ 7.3510-2 5.95 10 - 2 7.8610-2 Aforward1 5.29   10 - 2 _ 1.6810-3 1.42 10 - 3 2.8410-3 Abackward2 6.3410-3 1.0110-2 1.01   10 - 2 _ 3.41 10 - 3 Aforward2 3.92   10 - 1 _ 1.2410-2 1.1510-2 6.14 10 - 3 Abackward3 9.7810-4 6.4810-4 4.08 10 - 4 9.81   10 - 4 _ Aforward3 3.0810-3 4.1510-3 4.23   10 - 3 _ 1.50 10 - 3 Abackward4 4.1410-4 2.3410-4 1.12 10 - 4 7.45   10 - 4 _ Aforward4 1.0410-3 1.6010-3 1.69   10 - 3 _ 6.54 10 - 4  Second-order regression Abackward1 5.34   10 - 1 _ 5.9610-2 5.55 10 −2 6.9710-2 Aforward1 1.23   10 - 2 _ 1.4210-3 1.22 10 −3 1.6910-3 Abackward2 1.43   10 - 1 _ 9.7310-3 9.2210-3 3.28 10 −3 Aforward2 2.39   10 - 1 _ 9.3110-3 1.0910-2 5.66 10 −3 Abackward3 2.2510-4 4.02   10 - 4 _ 3.2410-4 1.52 10 −4 Aforward3 1.7810-3 2.6910-3 3.50   10 - 3 _ 1.35 10 −3 Abackward4 9.7610-5 1.36   10 - 4 _ 9.6810-5 7.95 10 −5 Aforward4 6.1010-4 1.0010-3 1.33   10 - 3 _ 5.22 10 −4Table 11
Relative generalization errorerrorG on the variance of the amplitude for different regression orders and correlation functions. Nonsymmetric case with 600 samples.
Gaussian Linear Exponential Cubic Zero-order regression Abackward1 9.42   1 0 - 1 _ 1.8610-2 1.3810-2 9.30 1 0 - 3 Aforward1 9.42   1 0 - 1 _ 4.0910-2 3.51 1 0 - 2 4.1710-2 Abackward2 9.35   1 0 - 1 _ 3.3910-1 1.7910-1 1.53 10 −1 Aforward2 9.44   1 0 - 1 _ 1.4510-1 8.8510-2 6.52 1 0 - 2 Abackward3 3.75 1 0 - 4 1.0110-3 4.8210-4 4.69   1 0 - 3 _ Aforward3 5.2810-3 5.36   1 0 - 3 _ 4.2810-3 3.39 1 0 - 3 Abackward4 1.0910-3 1.7010-3 9.23 1 0 - 4 4.39   1 0 - 3 _ Aforward4 4.58   1 0 - 3 _ 1.4110-3 1.07 1 0 - 3 3.9510-3  First-order regression Abackward1 3.42   1 0 - 1 _ 1.3110-2 1.2810-2 8.56 1 0 - 3 Aforward1 3.18   1 0 - 1 _ 3.2810-2 3.28 1 0 - 2 3.7110-2 Abackward2 7.53   1 0 - 1 _ 3.4510-1 1.8110-1 1.55 1 0 - 1 Aforward2 3.39   1 0 - 1 _ 1.4410-1 7.3910-2 5.28 10 −2 Abackward3 1.2710-4 1.7810-4 5.42 1 0 - 5 1.13   1 0 - 3 _ Aforward3 3.9710-3 5.28   1 0 - 3 _ 4.5010-3 2.02 1 0 - 3 Abackward4 5.0010-4 3.8910-4 1.89 1 0 - 4 1.21   1 0 - 3 _ Aforward4 1.96   1 0 - 3 _ 1.6010-3 1.6710-3 7.52 1 0 - 4  Second-order regression Abackward1 2.77   1 0 - 1 _ 1.2810-2 1.2610-2 8.44 10 −3 Aforward1 2.39   1 0 - 1 _ 3.10 10 −2 3.2110-2 3.5710-2 Abackward2 6.20   1 0 - 1 _ 3.4710-1 1.8810-1 1.67 1 0 - 1 Aforward2 3.08   1 0 - 1 _ 1.4410-1 7.4010-2 5.40 1 0 - 2 Abackward3 5.18   1 0 - 5 2.68 10 −5 2.7910-5 7.82   1 0 - 5 _ Aforward3 1.44   1 0 - 1 _ 3.3310-3 3.7210-3 1.64 10 −3 Abackward4 3.54   1 0 - 3 _ 3.0310-5 2.35 10 −5 4.6410-5 Aforward4 1.81   1 0 - 2 _ 1.4610-3 1.5010-3 7.14 10 −4For the interested reader and for a better visualization of the results that can be obtained via kriging metamodeling, AppendixB gives the best predicted maps of mean and variance of the maximal vibration amplitudes associated with the four backward and forward critical speeds.Finally, it can be concluded that the best kriging metamodeling (i.e., the best prediction based on the relative generalization errorerrorG; see the italic bold number for each amplitude in Tables 10 and 11) gives a very reliable approximation of the reference solution by performing a minimum number of calculations.


## 4. Conclusion

This study describes the use of a kriging based model in a rotordynamic problem. The efficiency of the proposed strategy has been tested successfully for a finite-element rotor model in the case of multiple uncertain parameters. Numerical results are validated by comparison with the results obtained by direct Monte Carlo simulations.It is shown that the methodology returns very accurate results for both the critical speeds and the maximum vibrational amplitudes even in the case of complex evolutions of the dynamical behavior of the rotor system. The use of kriging surrogate model has also a great efficiency in terms of the number of calculations needed. It was also illustrated that it is potentially possible to reduce the number of samples required to have predictive models by using some physical properties of the mechanical problem under study.Indeed, whatever the choice of the correlation function and the choice of the regression order, the use of kriging predictor allows a drastic reduction in the number of samples necessary for estimating all the possible dynamic behaviors in a predefined design of interest including potential uncertainties on physical parameters.From a practical point of view, the definition of relevant criteria and the nature of the output data to be predicted in connection with an initial understanding of the dynamic behavior of the mechanical system of interest have to be considered as the first step for an effective use of kriging approach. Even if it can be difficult to select the order and the correlation function for a specific simulation response, it is observed that in our problem under study the combination of a Gaussian correlation with a second-order regression (a cubic correlation with a second-order regression, respectively) works well in most cases for the prediction of all the critical speeds (all the maximum amplitudes, resp.).Finally, the numerical results on the prediction of critical speeds and vibration amplitudes also illustrate the fact that uncertainties have to be taken into account for proper robust designs of rotating machinery.



---

*Source: 1264619-2018-03-20.xml*