Research Article  Open Access
Kriging Metamodeling in Rotordynamics: Application for Predicting Critical Speeds and Vibrations of a Flexible Rotor
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.
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 secondorder moments is required for the defining of suitable strategies for robust design of rotor systems. For instance, strategies based on Tagushi rules use variancetomean 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 parameterdependent 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 finiteelement 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 Figure 1(a). This rotor system is composed of a shaft with four discs. It is supported at the ends by bearings.
(a)
(b)
2.1. Shaft Elements
The shaft is discretized into 10 Euler beam finite elements with circular crosssections. 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 the th shaft element defined asin fixed frame are given in Figure 1(b). The equation of motion for the th beam finite element can be written in the following form:where the superscript refers to the shaft. and are the translational and rotary mass matrices of the shaft element. , , and are the external damping and gyroscopic and stiffness matrices, respectively. 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 Nm^{−2} and the density kgm^{−3}. Table 1 lists the shaft properties.

2.2. Disk Elements
Each disk is modeled as a rigid disk. The equation of motion for the th disc finite element can be written in the following form:where the superscript refers to the disk. defines the nodal displacement vector for the th disk element that is defined as in fixed frame. and define the translational and rotary mass matrices for each disk. is the gyroscopic matrix. and correspond to the unbalance and gravity on the disk. Each disk is made of steel with Young’s modulus of elasticity Nm^{−2} and the density kgm^{−3}. Table 2 lists the disks properties.

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 twomode linear elastic spring elements. The elementary stiffness matrices are denoted by and for the first and second bearing supports, respectively.
2.4. OutofBalance 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 outofbalance 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 freedom in fixed frame (see Figure 1(b)), the outofbalance force for a given node of the finiteelement rotor system can be given bywhere and 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 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 aswhere defines the vibrational response of the rotor system. , , and 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 includes the contributions of each beam element (i.e., and for ) and each disk element (i.e., and for ). The matrix includes the effects of the shafts internal damping and damping of the supports. is the rotational speed of the rotor system. corresponds to the previous linear forces (i.e., unbalance forces) applied to the rotor system. Gravity denoted by 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, Figure 2 highlights the Campbell diagram that represents the evolutions of the forcing frequencies as a function of the rotor speed for two deterministic cases (with N·m^{−1} for the first case and 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 th backward and forward critical speeds (i.e., the backward whirl and forward whirl modes) diverge as the rotor speed increases (see Figure 2).
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 support and ) 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 and will be set independently for fixed values ranging from N·m^{−1} to N·m^{−1}. For each pair of these fixed parameters , uncertainties according to the seven other parameters will be taken into account.

For the reader’s understanding, the explicit functions of the th backward and forward critical speeds relative to the two deterministic parameters and are defined asfor . , , , and 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 aswhere and ( and , resp.) refer to the maximum amplitudes on the horizontal and vertical directions for the th backward critical speeds (the th forward critical speeds, resp.). In the following, these explicit functions will be noted , , , and (for ) 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 parameterdependent output , where may represent one of the maximum amplitudes, the backward or forward critical speeds defined previously while is a parameter vector which will refer to the vertical and horizontal bearing support and (hence ). Different kriging models dedicated to each one of the outputs are performed. It may be mentioned that 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 parameterdependent output can be metamodeled as follows:The latter expression is obtained by juxtaposing two mathematical objects. The first one consists in a linear combination of polynomial functions (for ) defining a regression with a predetermined order. The second object designated by is a realization of a Gaussian process which is assumed to be with a zeromean value and a covariance function given bywhere denotes the expectation operator.
The covariance function (9) is characterized by both and which define the process variance and the spatial correlation function, respectively. The latter which is a monotone function within measures how any two points and in the design space are spatially close to each other. It is defined bywhere denotes an 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 . 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.

As mentioned previously, defining a kriging predictor requires , , 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 and the associated output vector , the kriging parameters can be calculated as follows.
First, the parameter vector is given by the solution of the following optimization problem:where is a correlation matrix with entries defined by and denotes the determinant of a matrix.
The optimal regression parameter is given byas the optimal solution to the following classical least square problem:where is a matrix whose entries are given by .
The variance process can then be estimated throughAll parameters being determined, the kriging predictor is completely defined and then can be written as follows:where is the vector of the correlation coefficients between the new input vector and the old sampled input vectors .
3.3. Reference Solutions
First of all, a Monte Carlo simulation was carried out by considering two deterministic inputs and 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 combinations of the two deterministic inputs by sampling the twodimensional space parameter ; for each combination, simulations with the seven noise factors are performed (i.e., samples for each combination of the two decision variables are calculated). Then for each deterministic couple 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 bywhere is the number of samples for each combination of the two decision variables and .
All these results are given in Figures 3, 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 . For example, the mean values of the first and second backward critical speeds changed from Hz to Hz and from Hz to 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 and 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 or or more complex as illustrated in Figures 4(c), 4(d), and 4(e).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(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 parameter or (see Figures 5 and 6). For certain specific cases, the asymmetry of the bearing support (i.e., ) introduces further modifications to the vibration of the rotor. The mean values of the magnitude responses of the asymmetric rotor (i.e., ) at each backward critical speed are more important than the mean value of the vibration of the symmetric rotor (i.e., ). 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 considering combinations of these two deterministic inputs and with samples for each combination. This leads to a total of 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 and 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 Figures 3 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 speeds and and the vibration amplitudes and is performed. In a first step, one random experiment design with 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 Speeds and
With a random plan of size , the mean values and variances of the critical speeds and and the vibration amplitudes and are estimated against through kriging metamodels considering zero, first, and secondorder regressions and the linear, Gaussian, exponential, and cubic correlation functions.
The maximum and mean of the relative errors which are defined by and (for both the mean and variance of each backward or forward critical speed) between the kriging metamodels and the reference model are given bywhere corresponds to the previous reference results. are the results of the metamodel prediction value. defines the number of points in the design space (here as previously indicated in Section 3.3). It may be worth mentioning that each previous criterion or 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 where denotes the variance. As for the two previous criteria, the relative generalization error 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 of the design experiment is determined by 10 times the number of the considered parameters [24]. In this study, 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 ( in our study).
Results of the two criteria and 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.


