Application of Artificial Neural Networks to Calculation of Oil Film Reaction Forces and Dynamics of Rotors on Journal Bearings

Increase of energy efficiency and level of information system development of rotor machines in general requires improvement of theoretical approaches to research. In the present paper the problem of high-precision and high-performance computing programs development has been considered to simulate rotor vibrations. Based on two-layer feed-forward neural networks, numericalmodels have been developed to calculate oil film reaction forces to solve the rotor dynamics problems. Comparison has been done of linear and nonlinear approaches to solution of rotor dynamics problems, and a qualitative evaluation has been presented of accuracy and performance of a neural network approach compared to conventional approaches to rotor dynamics.


Introduction
Dynamic behavior of a rotating shaft on journal bearings is defined by the resonance vibration at critical speeds, self-excited vibrations due to the internal damping of the shaft's material, self-excited vibrations due to the oil film of journal bearings, oil whip, and flow-induced vibrations [1][2][3][4].Since the authors are closely connected with hydrodynamic lubrication, this paper is dedicated to oil whirl and oil whip problems.Oil film forces are obtained by means of solving the hydrodynamics equations and depend on many design and operational parameters of a bearing as well as on various specific physical effects that may occur, such as vibration, cavitation, and turbulence [2,4,5].The number of input parameters can reach several dozen [6].At the same time, development of technology in the direction of informatization requires simultaneous increase of accuracy and performance of the numerical models.One of the possible ways to provide both accuracy and high performance is using deterministic and stochastic approaches to modeling.Nowadays the methods and resources of implementation of neural network algorithms are developing quite intensively, which, based on the known data about the subject of study, allows building multifactor models of input and output parameters relations.Three basic levels of developed models could be distinguished.
Models of the first level have no direct interaction with the rotating machine; their training is based on the given numerical or experimental data.In [7] the results of study of an elastic rotor's dynamics on rigid bearings have been presented.Neural networks that have been applied are feedforward with two hidden layers.High accuracy has been demonstrated for models of rotor systems with roller-element bearings based on neural networks: accuracy from 82% up to 100% has been presented.Dynamics has been studied in a certain range of speed and for every case of speed separately, so the neural network had to be retrained.Some works are known, where neural network models are used to calculate the distributed characteristics, for example, pressure distribution in the oil film of the journal bearings [8,9].It should be noted that application of stochastic methods of rotor dynamics modeling allows skipping the stage of distributed values calculation.Instead, one could instantly determine the components of a resulting oil film force.In [10] dynamics of an imbalanced rigid rotor on fluid-film bearings has been studied.Authors determined the relation between three kinematic characteristics and two components of the resulting oil film force using a model of a neural network.The numerical tests were carried out for various angular speed and allowed determination of rotor trajectories.The synchronous, periodic, quasiperiodic, and chaotic motions for rotation speed within the range from 1000 to 10000 rpm were presented as the results of calculation.The authors suggest that there is a limitation to using linearized models in solving the problems of dynamics; however, they do not present a quantitative comparison of the numerical results obtained using linear and nonlinear models.Moreover, the reduction of the number of input variables has been done according to Poincare transformation of Reynolds equation, but possible outcomes of such substitution were not discussed in the paper.It is also not quite clear how one introduced kinematic variable allowed accounting for three given variables (two projections of rotor's lateral vibration speed and its rotating speed).It could be that the neural network had to be retrained for every particular case of the rotating speed.So, the results of this indeed interesting and important work remain in need of additional verification.A number of quite perspective papers should also be highlighted where neural network algorithms are used to account for complex specific phenomena [5,11] and to optimize the operational conditions of fluid-film bearings [12].Models of higher level can learn in real time using the processed measurement data.One of the most popular applications of such models is using neural networks to flaw detection of bearings in rotor machines [13][14][15][16].In [17] authors investigate the inverse model of the squeeze-film damper bearing (force input/displacement output) from empirical data and its application to a nonlinear inverse rotor-bearing problem.Comparison of experimental and numerical data showed high accuracy of the model based on recurrent neural networks.Models of highest level of complexity are able to generate some control impact [18,19].Today such approach is used to solve the problems of control in magnetic bearings [20][21][22].It should be, however, noted that the values of clearance in such bearings are an order of magnitude more than the same value in fluid-film bearings.
The present paper considers development of neural network algorithms of first level, that is, models with no direct contact to the rotor machine.Matters of quantitative evaluation of accuracy of application of neural network in comparison to known methods of oil film forces calculation and applicability of neural networks during simulation of nonlinear self-excited oscillations-oil whip-and during simulation of transient regimes of rotor's operation have been considered.

Mathematical Model
A rigid symmetrical imbalanced rotor is considered.Both ends of the rotor are supported by journal bearings of the same specification (see Figure 1).Given the chosen direction of coordinate axes, the pressure distribution and the resulting oil film force can be determined by means of numerically solving the Reynolds equation for a case of an isothermal flow of a viscous incompressible fluid [6]: where ℎ = ℎ( 1 ) is the oil film thickness (radial gap) function,  is the dynamic viscosity coefficient (assumed to be constant),  = ( 1 ,  3 ) is the pressure function, and   =   ( 1 ) is accordingly the tangential and normal components of the journal's surface velocity,  = 1, 2.
Radial gap function ℎ( 1 ) could be determined geometrically (see Figure 1) (coordinates and velocities of the center of the rotor are indicated with the capital letters, while the coordinates and the velocities of the fluid medium are indicated with the lower case letters): where ℎ 0 is the mean oil film thickness, ( 1 ,  2 ) are coordinates of the center of the journal,  = 2 1 / is the angular coordinate, and  is the bearing's diameter.The values of   over the surface of the journal could be determined using the following kinematic relations: where  is angular speed of the rotor and ( 1 ,  2 ) are projections of the velocity vector of the center of the rotor.Calculation of pressure distribution ( 1 ,  3 ) at some specific point in time is a boundary problem of solution of Reynolds equation ( 1) considering ( 2) and (3) with the following boundary conditions: the pressure on each side of the bearing is known to be ( 1 , 0) =  0 and ( 1 ,  3 ) =  1 ; along the  1 direction the following condition has to be met: (0,  3 ) = (, 3 ); / 1 |  1 =0 = / 1 |  1 = .Additionally, Gumbel's boundary condition for film rupture is used [6,23].
Given the pressure distribution, the components of the resulting oil film reaction force can be obtained as projections on the   (see Figure 1): where  3 is the length of the bearing.
The equation for lateral vibrations of the rotor's center can be presented in simple form [1][2][3][4]6]: where  is the mass of the rotor per one bearing,  is time, is the bearing load (centrifugal and inertia forces), Δ is the rotor's imbalance, and  is free fall acceleration.
Oscillations' velocity ⃗  of the rotor is significantly less than tangential velocity of journal's surface /2, so the solution of ( 5) is linked to solution of the steady-state fluid's flow problem solution.So, time as variable could be neglected during calculation of oil film force (4).This significantly simplifies the procedure of numerical solution of rotor dynamics problems [1,6].Three basic approaches have been considered and are characterized by different methods of oil film forces determination.
The first approach, called trajectory method [1][2][3][4]6], is based on solution of Reynolds equation ( 1) followed by integration of pressure distribution and determination of   (4).Such approach is considered referential due to high precision and however has low performance in computational terms, since the Reynolds equation has to be solved at each time step.
The second approach, called the method of linearization, of the oil film force determination [1][2][3][4]6] is based on the decomposition of reaction force components   into the Taylor series in the proximity of the equilibrium position    .The matrices of stiffness  and damping , calculated in the equilibrium position, shall be introduced as follows: Then the motion equation of the rotor (5) with given operational conditions of the rotor system and given the linearization of the film's reaction forces will be as follows [1-4, 6, 24]: The equilibrium point    can be determined using the trajectory approach.In many practical applications the equilibrium locus curve is calculated to solve problems of rotor dynamics [1][2][3][4]6].The points on this curve correspond to the equilibrium positions at various rotational speeds.Necessity to determine the equilibrium position or the equilibrium locus curve complicates the method of linearization.
The third approach is based on a neural network program module to approximate the oil film force's components   and to solve (5).A two-layer feed-forward neural network with inputs   ( = 1, . . .,  inp ) and outputs   ( = 1, . . .,  out ) is used (see Figure 2).The functionality of the network is determined by specifying the strengths of the connections, called weights , and the neuron transfer function  [25].The mathematical model of the used neural network is where  out () =  is a linear transfer function and  hid () = −1 + 2/(1 + exp(−2)) is a sigmoid symmetric transfer function,  out , and  hid , are weights,  out and  hid are thresholds of output and hidden layers, respectively, and   are the intermediate outputs of hidden neurons ( = 1, . . .,  hid ).
In the present paper the components of the oil film force   =   are considered as outputs, and four or five kinematic values, depending on the problem, are considered inputs: where  is the rotor's speed.
Neural network learning is essentially the solution of optimization problem of a minimum of the mean squared error function search: where  is a number of samples of inputs and outputs and Ỹ extremum of such function is implemented numerically.One of the criteria of neural network's quality is the correlation coefficient between outputs and targets   Ỹ [26].Target data Ỹ in the present paper is based on the numerical results, obtained using the trajectory approach.Stability analysis of the obtained results could be implemented based on the rotor's locus analysis using Poincare maps [10,[27][28][29].Moreover, the Routh-Hurwitz criterion could be applied to the linearization approach [1,10].

Simulation Results
Three basic approaches have been considered and characterized by different methods of oil film forces determination: the trajectory method, the linearization method, and the neural network method.Numerical solution of Reynolds equation (1) was implemented using the finite difference method.The obtained set of linear equations was solved using the Gauss method [30].The Levenberg-Marquardt backpropagation algorithm was used in training neural networks [25,31,32], which, in its turn, was based on the numerical results, obtained using the trajectory method.The rotor dynamics calculation program was based on the numerical solution of (5) or (7) using the explicit Runge-Kutta method of 4th and 5th order [30].During the development of the programs the GNU-Octave software and Octave's Neural Network Package [31] have been used.The GNU-Octave products are free software and use the Matlab compatible language.During simulations a PC with following characteristics has been used: Intel5 Core6 i5-4690K CPU 3.5 GHz, RAM 12 GB, x64 OS.
The results presented below correspond to various types of problems solved.

Rotor Dynamics with Fixed Speed and Low Vibrations.
A rotor system with the following parameters has been subjected to the study: rotor with 8.21 kg mass rotates at a constant speed of 4000 rpm under the load  and centrifugal force 1.8 and is supported by two plain journal fluid-film bearings with 40 mm diameter and 20 mm length with a mean radial gap of 100 m; lubricant is mineral oil with dynamic viscosity of 13 mPa⋅s; the pressure on each side of the bearing is 0.1 MPa.
First, trajectory approach has been applied.A preliminary numerical test neglecting the centrifugal force successfully  Then, a neural network has been developed for training with a number of neurons in the hidden layer  hid = 15.
As the input data a matrix of 4 × 10 5 elements was used, which consisted of numerical values of the coordinates   and the velocities   , with  0 = 1s calculated time of rotor's motion.The target data matrix of 2 × 10 5 elements consisted of the numerical values of oil film forces   .The matrices of input and output data were divided into three groups for training, validation and testing in fractions 0.6, 0.3, and 0.1, respectively.The procedure of adjustment has been implemented using application software [31].Training was done in 100 iterations, so-called epochs [25], in approximately 10 minutes.The convergence graph of a mean squared error of the network during training is presented in Figure 3. Training quality of a neural network is characterized by the following values of the mean squared error's value and correlation coefficients: MSE = 7.5 ⋅ 10 −4 N,   Ỹ = 9.9999 ⋅ 10 −1 .
The numerical tests for comparative calculation of the rotor's trajectories have been carried out for various values of  As a result of numerical tests a calculation time of 1 sec of rotor motion was determined given the action of dimensionless centrifugal force Δ 2 / = 1.6 (see Figure 4).It could be seen that the calculation time of rotor dynamics calculation using the linearization approach is the least in the whole range.Application of the neural network module increased the machine time up to 1.5-2 times under given conditions.Application of trajectory approach increased the calculation time up to 60-90 times compared to the linearization approach.It is indicative that calculation time of rotor dynamics using linearization approach and neural network approach is faster than the flow of real time in the whole range.This fact means that the application of linearization or neural network approaches in a real rotational machine control system would make it possible to obtain results of prognostic modeling [20,33,34].
As a result of another series of numerical tests, the trajectories have been calculated for a dimensionless centrifugal force Δ 2 / = 0.8, . . ., 2.5 (see Figure 5).Calculated rotor motion's time was set to 1 sec.The simulation started at the equilibrium point    , the trajectory approach taken as referential.The relative error of calculation has been chosen as quantitative criteria of accuracy evaluation: where  Traj  are coordinates, obtained using the trajectory approach, and  Lin(NN)  are coordinates, obtained using the linearization or neural network approach,  = 1, 2.
The quality analysis of stability of calculated trajectories has been implemented based on the Poincare maps [27][28][29].The point were displayed at moments of time that were multiples of the period of rotation of the rotor  = 1/.It could be seen from the Poincare maps (see Figure 5) that in the given range of loads the rotor's trajectories are quasiperiodic.This results from the fact that after a short transient process the trajectory's mapping becomes a dense pointed line, which closes if the calculation time increases by 100 times or more.This shows the presence of incommensurable frequencies in the oscillations, while the motion is not chaotic [27].It could also be seen from the Figure 5 that all three approaches to calculation of reaction forces gave the same result.The relative error of the center of mass location calculation with Δ 2 / = 0.8 during steadystate regime was 0.4% and 0.3% for the linearization approach and the neural network approach, respectively.With Δ 2 / = 1.6 the relative error remains small but accuracy differs by 10 times: 0.3% and 3% for the linearization approach and the neural network approach, respectively.Preservation of accuracy of the neural network approach could be explained by the fact that the model of training was carried out with value of the dimensionless centrifugal force Δ 2 / = 1.8 close to 1.6.Finally, with Δ 2 / = 2.5 the error was no more than 2% for the neural network approach and no more than 6% for the linearization approach.So, the neural network approach compared to the linearization approach given that the amplitude of oscillations is more than 0.1ℎ 0 has been proven to be significantly more accurate.
Further increase of the imbalance was not studied as it usually results in unstable trajectories.Increase of imbalance leads to increase of vibration's amplitude.That is why the linearization method could give inadequate results due to actual nonlinearity of the reaction forces at large eccentricities.The neural network approach gives good interpolation results, that is, within the learning sample.During extrapolation, deviation from the learning conditions could lead to unpredictable results and rapid increase of the error.

Rotor Dynamics for Various Speed Regimes and Transition from Low to High
Vibrations.The same rotor system as in Section 3.1 of the present paper has been studied.Imbalance was constant Δ = 10 −4 kg⋅m, and tests were run for various rotation speeds: 1000 rpm, 3000 rpm, 5000 rpm, and 7000 rpm.
The first series of simulations of rotor's behavior at various rotation speeds was implemented using the trajectory approach.Calculated motion time was set to 1 sec.Apart from trajectories, the location of equilibrium points has been determined, marked "+" in Figure 6(a).The second series of simulations has been implemented using the linearization approach.For that, matrices of stiffness and damping have been determined at four points on the equilibrium locus curve.
Then, four neural networks have been developed for training with  hid = 15 neurons in the hidden layer.For every network, matrices ⃗  = [ 1 ,  2 ,  1 ,  2 ] and ⃗  = [ 1 ,  2 ] of 4 × 10 4 and 2 × 10 4 elements were used as input and output Trajectory approach Linearization approach Neural network approach Trajectory approach Linearization approach Neural network approach Trajectory approach Linearization approach Neural network approach Trajectory approach Linearization approach Neural network approach Relative load, Δ • 휔 2 /g = 0.8

Trajectory approach Linearization approach Neural network approach
Trajectory approach Linearization approach Neural network approach Relative load, Δ • 휔 2 /g = 2.5   The calculation results showed that in the range of speeds from 1000 to 5000 rpm the vibrations of the rotor have the frequency equal to the rotation frequency (see Figure 6(b)).In this range of speed the oil whirl does not appear under given conditions, and the main reason of the occurrence of vibrations is the imbalance of the rotor.Additional calculations have shown that the critical speed of the rotor is 3300 rpm.In the range of speed from 6000 to 7000 rpm some increased vibrations occur that correspond to approximately twice critical speed (see Figures 6(a) and 6(b)).The vibrations of the rotor have the frequency close to the critical frequency (see Figure 6(b)).This is typical of oil whip [6].Here, linearization approach becomes unstable, and the neural network approach, on the other hand, remains stable with sufficient accuracy (see Table 1).
So, the neural network approach allows modeling of significantly nonlinear vibrations, such as oil whip.The necessity of development of a separate neural network for every rotation speed; however, it is a drawback of the applied approach.Consequently, in the next paragraph a neural network is considered with an additional input value, rotation speed of the rotor.

Rotor Dynamics during Transient
Regimes.Modeling of rotor dynamics during transient start/stop regime is more complex and costly in terms of calculation resources.The main advantages of the linearization approach that are essentially simple implementation and good computational performance manifest themselves weakly in this case.This is due to necessity to implement a number of intermediate calculations of equilibrium position, of stiffness and damping matrices, and due to interpolation procedures.In case of the neural network approach, such additional interpolation procedures are not needed.There is also a possibility of training a neural network based on one experiment and not on a series of them.This is why application of the linearization approach has not been considered in this section.
The transient regime has been studied using the test rig, shown in the Figure 7.The main parameters of the test rig It has to be noted that one of the alternative ways to train a neural network during monotonous smooth change of speed was not successful.The results obtained based on such neural networks were unstable.This could be explained by fact that every value of speed occurs in a monotonous function only once.In a step-like function, many intermediate values are absent, but the neural network manages to successfully interpolate them.
In Figures 8(a So, the solution of the rotor dynamics problems during transition regimes could be successfully implemented using neural networks.

Conclusion
An approach to neural network program module formation has been presented aimed at calculation of the reaction force's components of the fluid film depending on the position and oscillation velocities of the center of a rotor.Three approaches used in solving the rotor dynamics problem have been compared: the trajectory approach, the linearization approach, and the neural network approach to approximation of the oil film reaction forces.As a result of quantitative analysis, it has been determined that the use of linear and neural approximation allows decreasing the calculation time up to two orders of magnitude.It has also been determined that if the amplitude of the rotor's oscillation exceeds 10% of the mean radial gap, the neural network approach is up to 3 times more accurate compared to the linearization approach.Also, neural networks were used to simulate significantly nonlinear regimes of rotor's vibration with high vibrodisplacement and demonstrated good results in simulation of transient regimes.
The considered models of a rotor and the fluid film were relatively simple, which significantly improved performance using conventional methods and aided sampling for neural network's training.Further complication of the models, as well as training under condition of interaction with a real rotor machine, shall not significantly influence the main characteristics of artificial neural networks.

Figure 1 :
Figure 1: Calculation diagram of a journal bearing.

Figure 2 :
Figure 2: Feed-forward neural network with  inp number of inputs, with  hid neurons in hidden layer and with  out outputs.

Figure 5 :Figure 6 :
Figure 5: Comparative calculation results in the form of trajectories, Poincare maps, and error plots with different values of dimensionless load.
data, respectively.The matrices of input and output data were divided into three groups for training, validation and testing in fractions 0.6, 0.3, and 0.1, respectively.Training of each network was done in 1000 epochs in approximately 1-3 min.Training quality is characterized by the following values of mean squared errors and correlation coefficients: MSE = [2.43.1 2.7 14] ⋅ 10 −4 N,   Ỹ = [9.999.99 9.99 9.99] ⋅ 10 −1 .

2 )Figure 8 :
Figure8: Cascade frequency response during deceleration of the rotor from 3900 rpm down to 1000 rpm during 32 sec., calculated using the trajectory approach (a), neural network approach (b), and comparative results of the housing's vibration measurement and the rotor's vibration calculation (c).