An optimized bearing parameter identification approach from vibration response spectra

1. Introduction

Identification of rotor bearing parameters is an essential task in simplifying stability analysis procedure. Especially, the bearing parameters drastically affect vibration modes and responses of a rotating system. In practice, the high-speed rotors are often supported on various types of bearings and have unbalance and coupling forces leading to complex overall dynamics. Nowadays, the fluid film bearings are widely used in such rotors in reducing critical vibration amplitudes considerably due to their high damping and stiffness forces. These time-varying supporting forces often result in highly nonlinear response signals. Moreover, such force systems employ considerable computational memory which leads to relatively slower output performances. Equivalent linear bearing parameters if identified on the other way would result in better outcomes with respect to computational requirements. Further, these parameters are speed dependent [1] and require careful identification approaches. The identification studies of bearing parameters have been presented widely for plain and aerostatic journal bearings [2-4] estimated linear and nonlinear bearing stiffness of journal bearing using perturbation technique with two-dimensional Newton-Raphson iteration method. The methodology of prediction of sixteen dynamic coefficients for journal bearings in a rotor system from experimental unbalance responses was presented in Ref. [5, 6]. A method with multi-frequency excitation for measurement of equivalent stiffness and damping of the active magnetic bearing rotor was presented [7]. From the unbalance response [8, 9] and frequency characteristics [10], similar kinds of magnetic bearing parameter identification approaches were found. For ball bearing systems, the parameters were identified with simulated and experimental data [11]. Linear and nonlinear bearing coefficients of oil-free bearings including gas-foil [12-16] and gas-film bearings [17, 18] were obtained. Response based identification methodologies for tilting pad journal bearings were also noticed [19-21].

Hydrodynamic bearing forces are highly nonlinear and parametric in nature. For hydrodynamic journal bearings, field identification method for stiffness and damping characteristics was illustrated [22] using measured responses at both shaft and housing locations. An experimental approach was proposed [23] to estimate the stiffness and damping parameters via the least-square minimization under different operating conditions. A procedure to evaluate the rotor dynamic force coefficients of series bearing-supports was presented for impact and unbalance from the field measurements [24]. For identification of bearing parameters and unbalance from the measured responses, an optimization based strategy was proposed [25]. Qu et al. [26] explained the influence of the support stiffness on the engine vibration characteristics. Kriging surrogate model together with differential evolution optimization scheme was implemented [27, 28] to predict the bearing parameters. A modal parameter genetic time domain identification approach has been proposed [29] to study the characteristics of bearings using a multi-frequency signal decomposition technic. Prediction of bearing parameter information chart is; therefore, a very important task and a generalized methodology is, therefore, necessary to obtain the parameters conveniently. Although many studies are available in the literature, the estimation approaches based on correlating the real-time data with model-based outputs are found in limited papers. In the present work, the speed-dependent stiffness and damping parameters of the floating-ring bearing system are obtained from frequency response measurements followed by minimizing the average error in amplitudes between the actual and model-based response signals. Initially, the reference signal is obtained by analyzing the rotor-bearing system using three-dimensional beam element model of a rotor supported over floating ring bearings. The speed-dependent bearing coefficients are considered as variables and MPSO scheme is implemented to minimize the mean square error between the reference and linear idealized signals. Robustness of methodology is tested by introducing the noise into the measured reference signal. The main focus of this work is obtaining the unknown bearing stiffness and damping coefficients from the known responses in a less computational time. Measurements of lateral displacements at the bearings are simultaneously taken to achieve more accuracy in identification of parameters. Objective function considered is the sum of the mean square errors of displacements in both directions. Cross-coupled stiffness and damping coefficients are also calculated at the two bearings. In addition, the effect of stiffness of the bearing casing is accounted for obtaining the dynamic responses.

Remaining part of the paper is organized as follows: Section 2 describes the rotor-bearing system model and the expression of nonlinear bearing forces as well as the dynamic formulation of equivalent lumped parameter model of the rotor-bearing system. Section 3 presents the formulation of the objective function and optimization technique employed in the present work. Finally, the model validation along with optimization outcomes of a test case is illustrated in the results and discussion part.

2. Dynamic model of rotor bearing system

The more effective advantage of supporting action can be obtained at high speeds from dual-film hydrodynamic bearing systems such as full and semi floating-ring bearings. The response studies in high-speed rotor dynamic systems with floating ring journal bearings were thoroughly analyzed by simplified mathematical models [30-33]. The dynamic model of the flexible rotor dynamic system is formulated using quasi-finite element analysis with lumped floating ring masses considered at the bearing locations. The shaft is treated as flexible member and disks are treated as rigid. Each node has four degrees of freedom (DOF) including two translations ($u$, $v$) and two bending slopes (${\theta }_x$, ${\theta }_y$). By consideration of the bending and shearing effects the kinetic and potential energy of the shaft element can be expressed as:

The kinetic energy of each disk can be expressed as:

The virtual work done by unbalance forces at the disks can be expressed as:

From Hamilton principle:

with $\delta$ denotes the variational symbol, the equations of motion of the rotor alone is written as:

where, $q$ represents the displacement vector of size $N\times 1$, $F_u$ represents both the unbalance force and gravity force vector at the disks. $F_i$ is the inner oil film force and $\mathbf{M}$, $\mathbf{C}$, $\mathbf{G}$ and $\mathbf{K}$ denote respectively the assembled system mass, system damping, gyroscopic matrix of shaft and stiffness matrices assembly of size $N\times N$. The motion of the floating ring is identified with the help of inner and outer oil film hydrodynamic fluid forces, feed pressure of the lubricant and floating ring dead weight. The final motion of equation for the floating ring is written as:

where $q_R$ denotes the 4×1 displacement vector of the floating rings at both bearings, ${\mathbf{M}}_{\mathbf{R}}=diag\left(m_R,m_R,m_R,m_R\right)$ is the mass matrix of the rings with size 4×4. The outer and inner oil film forces of floating ring bearing can be represented as $F_o$ and $F_i$. By accounting the housing flexibility in the form of a single degree of freedom model with fixed mass, stiffness and damping coefficients, the combined simplified system is shown in Fig. 1.

Fig. 1Combined floating ring bearing forces with housing

The equation of motion of bearing housing of floating-ring is represented as:

where:

are the corresponding 4×4 mass, damping and stiffness matrices of the bearing housings at both left and right bearings, $q_B={\left[x_{B1}{y}_{B1}{x}_{B2}{y}_{B2}\right]}^T$denotes the 4×1 displacement vector of the bearing housings. By combining the system of equations for rotor and bearing, the assembled equations are represented as:

where:

are effective $\left(N+8\right)\times \left(N+8\right)$ square matrices, while $\tilde{q}={\left[\begin{array}{lll}q& q_R& q_B\end{array}\right]}^T$ is $\left(N+8\right)\times 1$ vector of displacements and $\tilde{F}={\tilde{F}}_u+{\tilde{F}}_b+{\tilde{F}}_g$ is resultant $\left(N+8\right)\times 1$ force vector of unbalance, bearing and gravity forces.

2.1. Floating-ring bearing with elastic housing

A floating ring bearing (FRB) has an annular ring placed in-between the journal and sleeve and there is a thin oil film in-between them. The bearing mid-plane consists of the circumferential feed grooves and the lubricant is fed from the journal to the sleeve via a bunch of feed holes located in the ring. The coordinate system considered is shown in Fig. 2.

Fig. 2Geometry of the floating-ring bearing with housing elasticity

With the help of the film pressure distribution, the hydrodynamic fluid forces can be derived from 2-D Reynold’s equation. Reynolds equations for both inner and outer lubricant films can be expressed as follows:

10

$\frac{1}{R_j^2}\frac{\partial }{\partial {\mathrm{\theta }}_i}\left(\frac{h_i^3}{12{\mathrm{\mu }}_i}\frac{\partial p_i}{\partial {\mathrm{\theta }}_i}\right)+\frac{\partial }{\partial z_i}\left(\frac{h_i^3}{12{\mathrm{\mu }}_i}\frac{\partial p_i}{\partial z_i}\right)=\frac{{\mathrm{\Omega }}_j+{\mathrm{\Omega }}_r}{2}\frac{\partial h_i}{\partial {\mathrm{\theta }}_i}+\frac{\partial h_i}{\partial t},$

11

$\frac{1}{R_o^2}\frac{\partial }{\partial {\mathrm{\theta }}_o}\left(\frac{h_o^3}{12{\mathrm{\mu }}_o}\frac{\partial p_o}{\partial z_o}\right)+\frac{\partial }{\partial z_o}\left(\frac{h_o^3}{12{\mathrm{\mu }}_o}\frac{\partial p_o}{\partial z_o}\right)=\frac{{\mathrm{\Omega }}_r}{2}\frac{\partial h_o}{\partial {\mathrm{\theta }}_o}+\frac{\partial h_o}{\partial t},$

where $p$ is the pressure of the oil film, and $\mu$ represents the viscosity of the lubricating oil, the subscripts $i$ and $o$ denote the parameters of inner oil film and outer oil film, respectively. While the subscripts $j$ and $r$ indicate the parameters between the journal and floating ring. $R_j$ and $R_o$ correspond to the journal and floating ring outer radius, respectively. $\theta$ is the angular coordinate for the inner and outer oil films. The axial coordinates of the inner and outer films are denoted by $z_i$ and $z_o$ respectively. The simplified expressions for oil film thicknesses and the film pressure distributions are expressed by considering Octvick’s theory of short bearings as:

12

$h_i\left({\theta }_i,t\right)=C_1-x_j\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_i-y_j\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_i,$

13

$h_o\left({\theta }_o,t\right)=C_2-X_r\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_o-Y_r\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_o,$

14

$p_i=\frac{3{\mu }_i}{h_i^3}\left(Z_i^2-\frac{L_i^2}{4}\right)\left\{\left[\left({\mathrm{\Omega }}_j+{\mathrm{\Omega }}_r\right)x_j-2{\dot{y}}_j\right]\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_i-\left[\left({\mathrm{\Omega }}_j+{\mathrm{\Omega }}_r\right)y_j+2{\dot{x}}_j\right]\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_i\right\},$

15

$p_o=\frac{3{\mu }_o}{h_o^3}\left(Z_o^2-\frac{L_o^2}{4}\right)\left[\left({\mathrm{\Omega }}_r{X}_r-2{\dot{Y}}_r{\mathrm{\Omega }}_r\right)\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_o-\left({\mathrm{\Omega }}_r{Y}_r+2{\dot{X}}_r\right)\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_o\right],$

16

$\left\{\begin{array}{l}{x}_j=X_j-X_r,\\ y_j=Y_j-Y_r,\end{array}\right.$

17

$\left\{\begin{array}{l}{\dot{x}}_j={\dot{X}}_j-{\dot{X}}_r,\\ {\dot{y}}_j={\dot{Y}}_j-{\dot{Y}}_r,\end{array}\right.$

where ($X_j$, $Y_j$) is the displacement vector of the journal center $O_j$ in the fixed reference frame. Also, $C_1$ and $C_2$ represent the static clearances of inner and outer film regions. The absolute displacement and velocity components of the floating ring center $O_r$ are ($X_r$, $Y_r$) and (${\dot{X}}_r$, ${\dot{Y}}_r$) while the absolute velocities of centers $O_j$ are denoted as $({\dot{X}}_j,{\dot{Y}}_j)$. The lower case letters $\left(x_j{y}_j\right)$, (${\dot{x}}_j$, ${\dot{y}}_j$) denote the displacement and velocity components of $O_j$ relative to $O_r$. The final expressions for inner and outer oil film force components are written as [34]:

18

$\left\{\begin{array}{l}{F}_{ix}\\ F_{iy}\end{array}\right\}={\mu }_i\left({\mathrm{\Omega }}_j+{\mathrm{\Omega }}_r\right)R_j{L}_i{\left(\frac{R_j}{C_1}\right)}^2{\left(\frac{L_i}{2R_j}\right)}^2\left\{\begin{array}{l}{f}_{ix}\\ f_{iy}\end{array}\right\},$

19

$\left\{\begin{array}{l}{F}_{ox}\\ F_{oy}\end{array}\right\}={\mathrm{\mu }}_{\mathrm{o}}{\mathrm{\Omega }}_r{R}_{ro}{L_o}_i{\left(\frac{R_o}{C_2}\right)}^2{\left(\frac{L_o}{2R_{ro}}\right)}^2\left\{\begin{array}{l}{f}_{ox}\\ f_{oy}\end{array}\right\},$

where the detailed expression for $f_{ix}$,$f_{iy}$,$f_{ox}$, $f_{oy}$ are given in the Appendix.

3. Methodology and optimization technique

The bearing force components are expressed in terms of displacements and velocities in bearing coordinates as:

where the terms $c$ and $k$ represent the unknown damping and stiffness bearing force coefficients. The direct and cross-coupled terms are denoted by suffices $xx$, $yy$ and $xy$ respectively. The suffix $B$ denotes the bearing support location. By substituting these bearing forces into the Eq. (9) and converting to the frequency domain, the system of equations can be rewritten as:

where, $H_{ij}\left(\omega \right)=k_{ij}\left(\omega \right)+i\omega c_{ij}\left(\omega \right)$ is impedance function, ($F_x$, $F_y$) and ($X$, $Y$) are the discrete Fourier transforms of external forces and displacements respectively. With the knowledge of the component displacements in bending directions at any location on the rotor, it is possible to compute the twelve force coefficients corresponding to each of the two bearings. In order to obtain the correct set of parameters, an error function defined in terms of $X$ and $Y$ amplitudes at the bearing nodes is considered at every operating speed. Mathematical formulation of the optimization problem in the current context is stated as:

where ${X_{NBF}}_i$ and ${Y_{NBF}}_i$ are the reference amplitudes obtained from the nonlinear force model, while ${X_L}_i$ and ${Y_L}_i$ are the corresponding displacement amplitude of frequency response via linear bearing forces. Here, $n$ denotes the total number of sample points considered in frequency-domain. This error $E$ is an implicit function of bearing coefficients which are defined with upper and lower bounds.

3.1. Particle swarm optimization (MPSO) with mutation

Conventional particle swarm optimization scheme is one of the robust metaheuristic optimization methods works on the behavior of flocking birds/fish during the food search [35]. Initially, with the random set of solutions, the system starts and searches for optimum value by updating the generations. The particle is described as each candidate solution and set of particles is known as a swarm. In a cooperative manner, they move in n-dimensional search space. The variable velocity performs the swarm movement of each particle. This velocity is influenced by social and local factors. Each particle moves through the search space based on the best positions found so far by itself ($P_{best}$)and the best position found by the swarm ($G_{best}$). In search space $R^k$ (with $k$ initial sets) for each particle, the objective function value is calculated. If $S\left(n\right)$ and $Vel\left(n\right)$ are position and velocity of each particle at $n$thiteration, movements of each particle are influenced by three factors (i) own direction search of Particle (ii) Particle beast position itself (iii) whole swarms best position. The position and velocity of every particle after iteration number $n$ is updated using the following equation:

23

$Vel\left(n+1\right)=u\times Vel\left(n\right)+c_1{r}_1\left(P_{best}-S\left(n\right)\right)+c_2{r}_2\left(G_{best}-S\left(n\right)\right),$

24

$S\left(n+1\right)=S\left(n\right)+Vel\left(n+1\right),$

where $u$ is called inertia factor of the particle, which often reduced in every cycle. Acceleration coefficients $c_1$ and $c_2$ describe the private (cognitive) and global (social) behavior of the system_. Also, $r_1$ and $r_2$ are random numbers between 0 and 1. $P_{best}$ is the best position of particle till the current iteration while $G_{best}$ is the best position of the group until current iteration. The algorithm converges to the best swarm by selecting the correct values of $u$, $c_1$ and $c_2$. Further, the velocity of the particle is bounded between the minimum and maximum values. Premature convergence can take place under different situations such as (i) the population has converged to local optima, (ii) the population has lost its diversity resulting in the search algorithm to proceed slowly. In this regard, to attain faster convergence without loss of accuracy, several modifications were suggested. To improve the population diversity and PSO’s performance, mutation is a powerful tool [36]. Here, a correction to the updated vector $S\left(n+1\right)$ in every cycle is introduced. This approach evaluates a mutation vector created from randomly selected three swarms (vectors, $Y_{B1}$, $Y_{B2}$ and $Y_{B3}$) in that generation. Fig. 3 shows the flowchart of the MPSO approach. The termination criterion employed in the present work is to achieve the maximum number of generations or to attain the error tolerance in successive objective function values whichever reaches earlier.

Mathematically, mutation vector ($Y_{mut}$) is expressed as [37]:

25

$Y_{mut}=Y_{B1}+\gamma \mathrm{*}\left(Y_{B2}-Y_{B3}\right),$

where ${\gamma }^{\mathrm{*}}\in$[0.9, 1] is the mutation constant.

The resultant population is modified using this mutation vector according to the following rule:

26

$Y_{new,i}=\left\{\begin{array}{l}{Y}_{mut,i},\left(rand<p_{cross}\right),\\ Y_i,\forall i=\mathrm{1,2},...k.\end{array}\right.$

Here, $k$ is the number of points in the population (swarm size). The $p_{cross}$ is crossover probability selected in the range of 0.1 to 0.9.

Fig. 3Flowchart of MPSO

4. Results and discussion

In order to examine the potentiality of the objective function, the simulated frequency response is initially generated from the rotor system with a linear bearing model having certain known input bearing force coefficients. This reference signal with input bearing force coefficients is provided to the optimization program and the bearing parameters are retrieved back through the error minimization procedure. The rotor is analyzed by finite element model using Timoshenko beam elements having two bending deflections and slopes at each node.

There are eight elements and nine nodes and total degrees of freedom are equal to 36. The discs are mounted at nodes 1 and 9, while the bearing nodes are at 3 and 7 as shown in Fig. 4.

Fig. 4Finite element model of the rotor- bearing system

Dimensional data of the rotor system is represented in Table 1.

The dynamic equations are solved by using fourth order Runge-Kutta time integration method with zero initial conditions. Fig. 5 shows the frequency response obtained at the left bearing node at a rotor speed of 5000 rpm. As the peak modes are occurring over the range 0-100 Hz, within this span the amplitudes are accounted in the objective function.

Fig. 5Frequency domain response at left bearing (kxx= 0.25 Mn/m, kxy= 0.12 Mn/m, kyy= 0.275 Mn/m, cxx= 300 Ns/m, cxy= 20 Ns/m, cyy= 399 Ns/m)

a)$X$-direction

b)$Y$-direction

The PSO parameters are taken as: $c_1=c_2=$ 2.1 and $u=u_{max}-j\left(\left(u_{max}-u_{min}\right)/j_{max}\right)$, where $u_{max}$ is maximum weight, $u_{min}$ is minimum weight, $j$ is iteration number and $j_{max}$ is maximum iterations. In present study, $u_{max}=$ 0.9 and $u_{min}=$ 0.4. The variable bounds are taken as: $k_{mn}\in$[10 kN/m, 5 MN/m], $c_{mn}\in$[10 N-s/m, 5000 N-s/m]. Fig. 6 shows the final error achieved for different swarm sizes at a rotor speed 5000 rpm. The minimum error occurs at a swarm size of 30.

Fig. 7 shows the fitness function convergence using proposed MPSO with a swarm size of 30 along with standard PSO scheme. It is seen that MPSO converges at a faster rate.

Table 2 shows the obtained bearing parameters as optimized design variables. The identified direct and cross-coupled parameters are found close to the reference values. Furthermore, the noise is added to the response signal for predicting the accuracy of identification. In this regard, random signal is added as a fraction of original signal. The percentage deviation is relatively small with the added input noise the frequency response.

Fig. 6Effect of swarm size in MPSO on converged error

Fig. 7Objective function convergence in PSO and MPSO

4.1. Rotor supported on floating-ring bearings

With floating ring nonlinear bearing forces, the rotor response is obtained from the same finite element model of the rotor. This time, the speed-varying equivalent bearing force parameters are obtained. Table 3 shows the bearing parameters employed in the simulation.

Table 3Floating-ring bearing parameters considered [32]

Parameter	Value
Bearing outer clearance $C_2$(m)	8×10-5
Bearing inner clearance $C_1$(m)	2×10-5
Mass of the ring $m_r$ (kg)	0.02
Inner film viscosity ${\mu }_i$(Pa-s)	0.006
Outer film viscosity ${\mu }_o$(Pa-s)	0.012

The time responses and frequency spectra in two directions are obtained at different operating speeds. The stiffness, damping, and mass of the bearing casing are respectively considered as 100 kN/m, 100 Ns/m and 0.1 kg. Fig. 8 shows the frequency spectra at a speed of 5000 rpm with and without accounting the casing stiffness. It is seen that there are two critical frequencies (at 83 Hz and 120 Hz) in the first case without casing flexibility. Furthermore, the subharmonic resonances resulting from hydrodynamic bearing forces are relatively small at this speed of operation [38]. With casing flexibility taken into account, it is observed that an additional mode is each direction results. The amplitude of the main dominating peak became small in both the directions as an absorber effect.

In order to identify the bearing coefficients, the response data is given as input to the optimization program. The same error function is further minimized by other well-known meta-heuristic optimization methods namely genetic algorithms (GA) [39] and simulated annealing (SA) [40]. In GA, population sizes of 30 along with mutation and crossover functions selected as constraint dependent and scattered type respectively. In SA also the function tolerance is considered as 10^-6. The annealing function used is fast annealing type and the reannealing interval considered as 100. The error function value obtained, and computation time taken in all three approaches are depicted in Table 4.

Fig. 8Frequency spectra at left bearing

a)$X$-direction

b)$Y$-direction

Table 4Performance comparison of various optimization schemes

Method	Error	Time (s)
MPSO	5.49e-17	140
GA	1.083e-14	155
SA	1.082e-11	185

It is seen that the MPSO is relatively good in terms of both accuracy and computational time. Table 5 shows the identified bearing stiffness coefficients at different operating speeds obtained from MPSO.

The corresponding identified damping parameters are given in Table 6.

Table 7 shows changes in the identified parameters at a rotor speed of 10,000 rpm with added noise in the frequency response. It is observed that the average error is well below 4 %.

Using the identified parameter data, the response at the left bearing node is further obtained from the finite element model and is shown in Fig. 9. The corresponding response with floating ring bearing forces at 10,000 rpm is also illustrated for comparison. The response spectrum is well matching with marked resonant peaks.

Table 5Identified stiffness coefficients of the bearings at different speeds

Sl. No	Speed (rpm)	Bearing-1			Bearing-2
		$k_{xx1}$ (MN/m)	$k_{yy1}$ (MN/m)	$k_{xy1}$ (MN/m)	$k_{xx2}$ (MN/m)	$k_{yy2}$ (MN/m)	$k_{xy2}$ (MN/m)
1	1000	3.1501	3.43912	1.2102	2.916827	3.145741	1.4651
2	2000	3.0110	3.45786	1.2113	2.81801	3.31343	1.4865
3	3000	3.0011	3.81297	1.2236	2.803257	3.621957	1.5135
4	4000	3.0001	4.16761	1.2258	2.758222	3.681405	1.5364
5	5000	2.9140	4.37417	1.3451	2.637015	3.83408	1.5412
6	6000	2.8466	4.62592	1.8684	2.45915	4.086198	1.5945
7	7000	2.8141	4.88159	2.2687	2.38048	4.15605	1.6124
8	8000	2.7957	5.27948	2.5587	2.193331	4.234727	1.6354
9	9000	2.8511	5.67308	2.6869	2.043291	4.327847	1.6589
10	10000	2.9534	5.95391	2.7246	2.006007	4.543091	1.6758

Table 6Identified damping coefficients of the bearings at different speeds

Sl.no	Speed (rpm)	Bearing-1			Bearing-2
		$c_{xx1}$ (MNs/m)	$c_{yy1}$ (MNs/m)	$c_{xy1}$ (MNs/m)	$c_{xx2}$ (MNs/m)	$c_{yy2}$ (MNs/m)	$c_{xy2}$ (MNs/m)
1	1000	0.00434	0.0049	0.00028	0.004291	0.004923	0.00035
2	2000	0.004238	0.00479	0.00029	0.003912	0.004835	0.00036
3	3000	0.004138	0.00492	0.00031	0.003839	0.004723	0.00039
4	4000	0.003924	0.0050	0.00035	0.003792	0.004957	0.00041
5	5000	0.003748	0.00511	0.00039	0.004001	0.005085	0.00042
6	6000	0.00386	0.00549	0.00041	0.004182	0.00533	0.00043
7	7000	0.00398	0.00568	0.00043	0.004314	0.005711	0.00046
8	8000	0.004002	0.00583	0.00043	0.004408	0.005877	0.00049
9	9000	0.004149	0.00606	0.00043	0.004693	0.005916	0.00050
10	10000	0.004316	0.00627	0.00043	0.004942	0.006178	0.00512

Table 7Identified bearing parameters with added noise at 10,000 rpm

Parameter	Without noise	With 5 % noise	With 10 % noise
$k_{xx1}$ (MN/m)	2.9534	2.974	2.991
$k_{xy1}$(MN/m)	2.7246	2.754	2.785
$k_{yy1}$(MN/m)	5.95391	5.995	6.024
$k_{xx2}$ (MN/m)	2.00600	2.012	2.098
$k_{xy2}$ (MN/m)	1.6758	1.694	1.701
$k_{yy2}$ (MN/m)	4.543091	4.578	4.597
$c_{xx1}$(MNs/m)	0.004316	0.00451	0.0047
$c_{xy1}$(MNs/m)	0.00043	0.00047	0.0005
$c_{yy1}$(MNs/m)	0.00627	0.00635	0.0064
$c_{xx2}$(MNs/m)	0.004942	0.00501	0.0052
$c_{xy2}$(MNs/m)	0.00512	0.00534	0.0054
$c_{yy2}$(MNs/m)	0.006178	0.00628	0.0064

Fig. 9Comparison of frequency response at left bearing node (10,000 rpm)

5. Conclusions

In this work, bearing force parameter identification procedure was illustrated with available frequency spectra in a rotor dynamic system. Flexible rotor system was analyzed by using a finite element model and the frequency domain responses at different rotor speeds were considered. An error-based formulation in terms of response amplitudes at bearing nodes was employed and the speed dependent stiffness and damping parameters of the bearings were identified via modified particle swarm optimization with the mutation. The methodology was found to be reliable and predicts the coefficients with limited computation effort. Robustness of methodology was tested by introducing the noise into the measured reference signals. The average error was not exceeding four percent. The average time for function evaluation can be further minimized by employing an approximate solution technique for solving nonlinear dynamic equations for achieving directly the frequency spectrum without the need of time domain analysis. Alternatively, surrogate models via neural networks may also be employed to avoid the complex time-domain calculations.