Simplified transformation matrices of journal bearings in vertical application

.


Introduction
Fluid-film journal bearings, in general, are employed in rotating machines to provide radial support to the rotor.Journal bearings in horizontal rotors often operate under static load conditions due to the weight of the rotor.In vertical rotors, however, the weight acts vertically downward; hence radial load is relatively lower, leading to susceptibility to fluid-film-induced instability.Vertical rotors are used in many different applications, such as in hydropower units and pumps.
In rotordynamic analyses, a reliable fluid-film bearing model is important for accurately estimating machine dynamics.A thin fluid-film lubricant under a bearing is mathematically modeled by the Reynolds equation or Navier-Stokes equations.The governing equations are partial differential equations and must be solved either analytically or numerically.In the former, the model must be simplified by employing certain approximations, like in the case of short (Ocvirk, 1952;Ruggiero and Senatore, 2007;Miraskari et al., 2018) and long (Myers, 1984;Amamou and Mnaouar, 2014) bearings.On the other hand, in advanced simulations, the lubrication model is solved numerically using finite element methods, finite difference methods, or the Pseudo-Spectral method.
The fluid-film-induced forces are commonly described using eight linearized stiffness and damping coefficients, which are calculated using either infinitesimal (Lund and Thomsen, 1978;Lund, 1987) or finite perturbation approaches (Qiu and Tieu, 1996) at the equilibrium location.Childs (1993) and Ishida and Yamamoto (2012) described a method to calculate linearized bearing coefficients based on short and long-bearing approximations.Although the linearized coefficients are an adequate approximation of the bearing dynamics for many cases, they are not always valid.Several papers (Hattori, 1993;Muller-Karger et al., 1997;Machado et al., 2018;Zhao et al., 2005;Meruane and Pascual, 2008;Choy et al., 1991) demonstrated non-linearity in journal bearings that can influence the dynamics of the rotor-bearing system.Zhao et al. (2005) indicated the presence of non-linearity in bearing forces due to larger excitation amplitude.The authors proposed non-linear models based on Taylor series expansion and identified the bearing coefficients from the experiment.Using a similar approach as Refs.Zhao et al. (2005), Meruane and Pascual (2008) showed the importance of non-linear bearing models for predicting transient solutions and instabilities.
In many advanced simulations, numerical methods are employed to calculate the non-linear bearing forces, and the lubrication model is solved simultaneously with the rotor-bearing system at each time step.As this process can be computationally expensive, researchers have attempted to simplify the simulation procedure and reduce the amount of time required for computation.For systems with many DOFs and complex systems that require Multiphysics, the traditional simulation approach, carried out by solving the fluid film lubrication model at each time step, is impractical.To the authors' best knowledge, most commercial software available today perform different types of simulations based on bearing characteristics calculated at a static location, mainly influenced by the dead weight of the rotor.However, vertically oriented rotor systems operate at no or minimal static load as the weight is carried axially.The radial load is often rotary, and the dynamic simulation requires calculation of bearing coefficients at each journal location.Nässelqvist et al. (2014) presented a bearing model by pre-calculating the bearing coefficients of a four-shoe TPJB and representing them by one-dimensional polynomial equations.The proposed model, however, ignored the cross-coupling coefficients, which were later considered by Synnegård et al. (2016).This approach was extended by Benti et al. (2022), and the bearing model was further represented by two-dimensional polynomial equations, which allowed transient simulations under variable rotor speed.analysis, require many iterations, which take several hours or days of computation.Such simulations can be heavy even for the simplified model presented in Benti et al. (2022).Improving the computational efficiency would significantly benefit and enable Rotordynamic simulations within a reasonable time.
The simulation of the rotor-bearing dynamic system often includes the computation of bearing force coefficients, followed by coordinate transformation.This paper aims to simplify the transformation matrices and improve computational efficiency.The unbalance response of a vertical rotor with two TPJBs was simulated and analyzed for two operational conditions, i.e., constant and variable rotor speeds.The proposed model was compared with the original model and validated by experiments.

Experimentation
The test rig consisted of a vertically oriented mid-span rigid rotor supported by two identical flexible supporting structures, as shown in Fig. 1.The supporting structures were mounted at the top and bottom end of the rotor to provide radial support.A series of measurements were carried out, and the performances of the bearings were evaluated for different unbalances and operating speeds.Both constant and variable rotor speeds were considered.In the former, tests were carried out under three rotor's operational speeds (500 rpm, 1500 rpm, and 2500 rpm), whereas, in the latter, the angular velocity of the rotor ramped up to 2500 rpm by a constant angular acceleration of 2π rad/s 2 .Sections 2.1 and 2.2 present technical details of the rotor and the eight-shoe TPJBs, respectively.

Rotor description
The mid-span rotor was driven by an electric motor with an angular velocity of up to 2500 rpm.A slender stinger integrated with two jaw couplings was used to connect the rotor with the electric motor, and it was designed so that it does not introduce a radial force to the rotorsupport system.Each supporting structure consists of an eight-shoe TPJB and a bracket with a stiffness of 500 MN/m and no damping (Nässelqvist et al., 2012).Note that TPJBs installed in Swedish hydropower units consist of eight pads or more (Gustavsson et al., 2019).The technical specification of the test rig is presented in Table 1.
An unbalance mass was mounted on a disk at a distance of d = 0.7 m from the center of the rotor.The rotor trajectory of the journal center, and forces acting on the bearings at the two bearing locations were measured.Each bearing was equipped with inductive proximity sensors to measure the shaft displacement, whereas full Wheatstone bridge strain gauges were mounted on each bracket to measure the bearing reaction forces.Besides, an optical sensor (≈ ±1 rpm accuracy) was used to measure rotor's angular velocity.All the sensors were calibrated onsite using certified calibration equipment.

Bearing description
The hydrodynamic TPJBs have eight identical pads with rocker pivots radially located 45 • apart.Each pad was made of a white-metal (babbitt) lining and backed with steel-1650.The bearings were supplied with 0.01 MPa lubricant (Q8 Handel oil), and the lubricant temperature was measured during the test and ranged from 23 • C to 33 • C. Fig. 2 and Table 2 show the schematic representation and the detailed specification of the eight-shoe TPJB, respectively.

Numerical model
The simulation procedure was performed in two steps, which are discussed in Sections 3.1 and 3.2.Section 3.1 describes how the bearing models are developed by pre-calculating the stiffness and damping coefficients from the fluid film lubrication model.Besides, the simplified transformation matrices of the bearing force coefficients, which is the main objective of this paper, were briefly discussed.Section 3.2 presents the second step of the simulation procedure, which is the numerical simulation model of the rotor-support system.

Bearing model
The bearing was modeled as two DOFs, and the equation of motion consisted of damping and stiffness coefficient matrices.Eq. ( 1) shows the total bearing forces, which are the summation of the stiffness and damping force components. [ where, f i is the bearing force, x and y are displacements, ẋ and ẏ are velocities, C ij and K ij are the damping and stiffness coefficients, respectively.Subscripts i and j denote the directions in the fixed coordinate.The bearing coefficient matrices in Eq. ( 1) are in the fixed frame of reference and are calculated by coordinate transforming the local coefficients matrices using Eqs.( 2) and ( 3). [ where, The coordinate transformation requires matrix multiplications, which must be employed at each time step in the numerical simulation procedure.The process is computationally expensive, and the main aim of this study is to simplify the procedure by avoiding matrix multiplications.The first simplification modifies the stiffness matrix, and it comes from the fact that the displacement in the local η-axis is always zero.This is because, at each time step, the bearing coefficients are calculated at a stationary journal center location.Thus, the eccentricity line is always fixed to the ξ-axis and normal to the η-axis.Furthermore,  G.B. Benti et al. in the second simplification, the damping matrix was modified.The trajectory of the journal center is approximated to be centered circular for TPJBs with many pads when installed in vertical rotors, resulting in low velocity in the radial direction ( ξ ≈ 0).The bearing forces in Eq. ( 1) can be re-written as follows: Thus, excluding η from Eq. ( 4) removes the second column of the stiffness matrix.Besides, assuming ξ = 0, the first column of the damping matrix vanishes.Therefore, for η = ξ = 0, Eq. ( 4) is reduced to Eq. ( 5). [ Rearranging Eq. ( 5), the bearing coefficient matrices are expressed by the local stiffness (K ξξ and K ηξ ) and damping coefficients (C ηη and C ξη ) as shown in Eq. ( 7).The bearing coefficient matrices require no matrix multiplications, which was the case for the original matrices in Eq. ( 1). [

Polynomial representation of the maxima and minima bearing coefficients
A similar approach was employed as in Benti et al. (2022) to calculate the local stiffness and damping coefficients.As shown in Fig. 3, the bearing coefficients of the eight-shoe TPJB vary periodically due to the number of pads; hence, they can be described by harmonic functions, Eqs. ( 8) and ( 9).This, in other words, means once the maxima and minima values are known, the coefficients can be calculated for any given eccentricity angle (α).
where,  Seal Research, 2023).Each pad consisted of three degrees of freedom.However, the radial stiffness of the rocker pivot was modeled as rigid, and only the rotating DOF was retained, while the other two DOFs were eliminated.The full dynamic solution of the eight-she TPJB was computed and the eight Rotordynamics coefficients were calculated.Note that the influence of excitation frequency was not investigated in this paper.The bearing coefficients in the local ξ and η directions were calculated for certain number of angular velocities and eccentricity ratios.Computations were performed at six different rotor angular velocities (250 rpm, 500 rpm, 1000 rpm, …, 2500 rpm) and nine eccentricity ratios (0.01, 0.1, 0.2, …, 0.8), which gives a total of 36 individual computations.The calculations were performed under prescribed eccentricity, and it is relatively faster as it is less sensitive to convergence issue, which was the case for the bearing simulations under prescribed load.For instance, for rotor speed of 2000 rpm and 0.4  For each calculation, the maxima and minima values of the bearing coefficients were extracted, and the dataset was fitted with two-dimensional polynomial regression.A MATLAB inbuilt function (poly) was used to perform the regression.The choice of the order of the regression was based on the goodness-of-fit (Table A1), and the regression coefficients of the best-fitted model (poly51) are illustrated in Table 3.

Equivalent stiffness coefficient
When solving the equation of motion of the rotor-support system (Eq.( 11)), which will be discussed later in Section 3.2, linearized bearing stiffness coefficients are required.However, the maxima and minima local stiffness coefficients presented in Table 3 are nonlinear, and they cannot be directly used in Eq. ( 8).Therefore, for a given eccentricity ratio (ϵ 1 ), the linearized equivalent maxima and minima stiffness coefficients were calculated by integrating the local coefficients over the eccentricity ratio, as shown in Eq. ( 10).
where, Kk ij is the maxima or minima equivalent stiffness coefficients, and K k ij is the maxima or minima stiffness coefficients presented in Table 3 and ϵ 1 is a given eccentricity ratio.

Rotor-support model
The test rig was modeled as a rigid rotor with two isotropic flexible supports.The support system consisted of the bearing and bracket connected in series and was represented by springs and a damper.Eq. ( 11) is the mathematical representation of the forced response of the rotor-support system.The model had a total of eight DOF; four DOF for the rotor and four DOF for the two supports.
where M and G are the mass and gyroscopic matrix, K B and C B are the bearing stiffness and damping matrices, K S is the stiffness matrix of the supporting bracket, φ is the angular acceleration, which is the rate of change of the angular velocity ( φ).For simulations with the constant angular velocity of the rotor, the angular acceleration ( φ) is zero.Eq. ( 11) was solved numerically using a MATLAB inbuilt function (ode15s).
The force (F u ) is due to the unbalance mass attached to the disk at Ddistance from the center of the rotor and acted only in the x and y DOF.

Results and discussions
The dynamic responses of the rotor-support system were simulated for different unbalance magnitudes and when the rotor was spinning at constant and variable speeds.For the constant speed case, the simulations were carried out for 200 shaft revolutions, and the transient part of the response at the beginning of the simulation was disregarded.On the other hand, simulations under variable rotor speeds were carried out when the speed ramped up linearly to 2500 rpm.
Sections 3.1 and 3.2 compare the results from the proposed model with those from the original model and experiment, respectively.The bearing matrices of the original model were obtained by adopting coordinate transformation, which entails matrix multiplication, as demonstrated in Eqs. ( 1) -(3).On the other hand, the bearing matrices of the proposed model required no matrix multiplications for coordinate transformation.As shown in Eq. ( 7), the modified bearing model is described by four local coefficients, whereas the original bearing model, Eq. ( 1), is represented by the full eight coefficients in the fixed coordinate.

Original vs proposed model
As the center of the journal moved from the bearing center, the two TPJBs generated the fluid-film induced reaction forces, which were formulated by the stiffness and damping coefficients in the rotor-bearing simulations.Figs. 4 and 5 depict the coefficients of the upper/lower TPJB (the upper and lower bearings produced similar results) when the rotor was running at 2500 rpm and subjected to 1.7 × 10 − 3 kg⋅m and 5.9 × 10 − 3 kg⋅m unbalance magnitudes, respectively.The corresponding elements of the matrices of the two models are plotted together for Table 3 Computed coefficients of the two-dimensional polynomial equation (poly51) as a function of centered and scaled eccentricity ratio and rotor speed.The eccentric ratio was normalized by mean 0.4011 and std of 0.2589, whereas the rotor speed was normalized by mean 1292 rpm and std of 803.5.    on the force they generate, which is the product of the coefficient matrices and the corresponding displacement or velocity vectors.Therefore, the stiffness (f K x and f K y ) and damping (f C x and f C y ) force components were calculated.Figs. 6 and 7 portray the simulation results obtained from the two models, and the stiffness and damping force components are presented as a function of the eccentricity angle.For lower unbalance magnitude, the two models display similar results.Besides, the stiffness force component was dominant and the ratio of the resultant stiffness and damping force components (f C /f K ) ranges between 0.79-0.83.When the orbit radius and velocities of the journal center increased due to large unbalance magnitude, the results from the two models deviated, particularly the damping force components.Although the deviations are small and insignificant in terms of force amplitude, it apparently affected the shape of the force profile, which will be discussed later in this section.Similarly, the stiffness force component appeared dominant, and the ratio f C /f K ranges between 0.56-0.72,slightly less than the result from the simulation with the lower unbalance.
Fig. 8 shows the simulated trajectory of the journal center and the force acting on the bearings at Ω = 2500 rpm, and three unbalance magnitudes, i.e., m⋅d = 1.7 × 10 − 3 kg⋅m, 3.8 × 10 − 3 kg⋅m and 5.9 × 10 − 3 kg⋅m.For the proposed model, the modified bearing stiffness matrix was first used individually in the simulation (Modified K & Original C) and then together with the modified damping matrix (Modified K&C).The results indicate that the proposed model with the modified bearing stiffness matrix was the same as the original model, regardless of the rotor speed or orbit amplitude.This condition is generally true since η is always zero in the rotating coordinate.Similarly, at a lower orbit radius, no significant differences were observed between the results from the models with the Original K&C and Modified K&C.As the orbit radius increased, however, discrepancies in the force predictions emerged mainly in terms of phase.At the orbit relative eccentricity of 0.74, the phase of the octagonal-shaped force profile predicted by the proposed model preceded those from the other two models by around 3 • .Nevertheless, there was no significant difference between the two models in force amplitudes, and the amplitude ratio (r) of the maximum resultant force was 0.977.Besides, for some applications, the unbalances used in this paper might be too large for the deviation to occur.The most significant disagreement between the models appeared at 5.9 × 10 − 3 kg⋅m, which was ten times larger than the permissible residual unbalance suggested by Iso standard requirement for rigid rotor (ISO 1940/1, and G6.3 grade balancing) (ISO (International Organization for Standardization), 2003).
Moreover, the unbalance responses were investigated for variable rotor speed.A series of simulations were conducted when the rotor speed increased linearly from 250 rpm to 2500 rpm by angular accelerations of 0.5π rad/s 2 , 2π rad/s 2 , and 4π rad/s 2 .Fig. 9 shows the average orbit and force amplitudes calculate at 12 speed ranges with a block size of 200 rpm.It is found that both the orbit and force responses increased with the rotor speed and unbalance.Besides, the results from all three models were similar, and no significant differences were observed, regardless of the unbalance magnitude and angular acceleration.
The computational efficiency of the proposed model was investigated.For simulation, MATLAB 2019a software installed on a standard consumer laptop (Intel Core i7-8850H and CPU at 2.60 GHz) was used, and the elapsed computation time was measured by an inbuilt function (tic-toc).For an unbalance response simulation at Ω = 2500 rpm, and for 200 shaft revolutions, the simulation model with the modified bearing matrices (modified K&C) was two times faster than the model with the original bearing matrices.The time saving would be much more significant for long-duration simulations.For instance, bifurcation simulation usually performed for long duration or many shaft revolutions, ranging up to 10 6 .Reducing the computation time by halve in such cases would make a huge difference.

Experiment vs proposed model
The numerical results from the proposed model were further validated by experiment.Both the bearing reaction forces and the trajectory of the journal center at the two bearing positions were investigated, and the experimental and simulation results were compared for operations under constant and variable rotor angular velocities.All the measured Fig. 6.The stiffness and damping force components of the upper /lower TPJB in the x and y directions: Ω = 2500 rpm and m⋅d = 1.7 × 10 − 3 kg⋅m.
responses were filtered by a band-pass filter with a cut-off frequency between 3 Hz and 500 Hz.
Fig. 10 shows the orbit and force responses when the rotor was subjected to 5.9 × 10 − 3 kg⋅m unbalance magnitude and operated at three constant rotor speeds, i.e., 500 rpm, 1500 rpm, and 2500 rpm.The response measurements at the two bearings show a similar trend, albeit some differences in their amplitudes.Considering the filtered mean amplitudes (Fm), the deviation of the maximum orbit amplitudes at the two bearing positions can range up to 15%, whereas 9% for the bearing reaction forces.These deviations are not unexpected since the performance of the two bearings may differ due to unavoidable differences in bearing parameters, like geometry.For simulation, however, the performance of the two bearings was assumed to be the same, resulting in no significant difference between the responses at the two bearing positions.
The experimental and simulation results were compared and found to have similar trends.At 2500 rpm rotor speed, the average amplitude of the measured orbit and force deviated from those in the simulation up to 7% and 15%, respectively.Besides, both the orbits and forces were circular-shaped at low orbit amplitudes.As the rotor speed increased, the orbit amplitude increased.This resulted in an octagonal-shaped force with peaks and valleys, which were influenced by the number of pads.The maximum bearing force appeared at the pegs, where larger stiffness and damping coefficients are located.
Furthermore, the results from the experiment and simulation were compared for operations under variable rotor speeds.In Fig. 11, the upper bearing forces obtained from measurement are plotted along with the corresponding simulation results.For both cases, the force amplitudes increased with the rotor speed.Unlike results from the simulation, the experimental result contains an amplitude peak close to 1000 rpm, which was due to the test rig's structural resonance.Besides, as shown in Fig. 12, the profile of the simulated and measured forces at different rotor speed regions are similar, and the amplitudes agree with some discrepancies.Inaccurate bearing parameter prediction could be ascribed to these deviations.The fluid-film lubrication simulation was carried out on the eight-shoe TPJB, taking the nominal bearing  parameters.However, these parameters may vary within some tolerance range.For instance, the pads were assumed to be identical with uniform geometric characteristics, which is not true.Besides, the inevitable parameter difference between the upper and lower TPJBs was not considered in the simulation.

Conclusions
The transformation matrices of the bearing force coefficients were successfully simplified.Compared to the simulation model with the original transformation matrices, the proposed model has considerably increased computational efficiency while preserving high accuracy.Discrepancies, which emerged due to the reduced damping matrix, affected the shape of the bearing force, albeit there were no significant differences in amplitudes.For bearings with many pads, these discrepancies are insignificant, as the bearing coefficients became less dependent on the angular position of the journal center.Furthermore, the proposed model was validated by experiment, and the predicted responses agreed with those from the measurement.The simplifications that were employed on bearing stiffness and damping matrices improved the efficiency of the simulation model.The proposed model is, therefore, effective for rotor-dynamic simulations and response analysis.

Funding
The research presented in this thesis was carried out as a part of "Swedish Hydropower centre -SVC".SVC has been established by the Swedish Energy Agency, Energiforsk and Svenska Kraftnät together with Luleå University of Technology, KTH Royal Institute of Technology, Chalmers University of Technology, Uppsala University and Lund University.
and C ij are the stiffness and damping coefficients, i and j in the local ξ and η directions, Kmax ij and Kmin ij are the maxima and minima equivalent stiffness coefficients, C max ij and C min ij are the maxima and minima damping coefficient, ϵ is the eccentricity ratio, α is the eccentricity angle and Ω is the rotor speed.A Navier-Stokes equations-based fluid film lubrication model was solved numerically using commercial software, RAPPID (Rotordynamic

Fig. 4 .
Fig. 4. The stiffness and damping coefficients of the upper/lower TPJB as a function of eccentricity angle at Ω = 2500 rpm and m⋅d = 1.7 × 10 − 3 kg⋅m.

Fig. 5 .
Fig. 5.The stiffness and damping coefficients of the upper/lower TPJB as a function of eccentricity angle at Ω = 2500 rpm and m⋅d = 5.9 × 10 − 3 kg⋅m.

Fig. 7 .
Fig. 7.The stiffness and damping force components of the upper /lower TPJB in the x and y directions: Ω = 2500 rpm and m⋅d = 5.9 × 10 − 3 kg⋅m.

Fig. 8 .
Fig. 8.The orbits and forces obtained from the three models at Ω = 2500 rpm and three unbalance magnitudes.
Fig. B.1.The maxima local stiffness coefficients fitted with poly51 (MATLAB).The dots represent the computational result from the fluid-film lubrication model, RAPPID.

Fig. B. 2 .
Fig. B.2.The minima local stiffness coefficients fitted with poly51 (MATLAB).The dots represent the computational result from the fluid-film lubrication model, RAPPID.

Fig. B. 3 .
Fig. B.3.The maxima local damping coefficients fitted with poly51 (MATLAB).The dots represent the computational result from the fluid-film lubrication model, RAPPID.
Some simulations, like bifurcation

Table 1
Technical specification of the test rig.
Description Value Ds Shaft diameter (mm) 49.84 L Shaft length (mm) 500 Dd Disk diameter (mm) 168 t Disk thickness (mm) 100 Direction of rotation Counter clockwise Ks The stiffness of the bracket (MN/m) 500 Ms Rotor mass (kg) 24.74