Experimental and theoretical analysis of a rigid rotor supported by air foil bearings

The popularity of compressors utilizing foil bearings is increasing. Their mechanical design is challenging, and an accurate prediction of the bearing coefficients is important. A mathematical model taking into account the foil structure, and the detailed geometry of a three pad foil bearing are presented. The steady state solution and dynamic coefficients are obtained through zeroth and first order perturbed equations respectively. Analysis of the foil structure reveals the importance of distinguishing between a static foil stiffness for the zeroth order equation and a dynamic stiffness for the first order equation. Calculated bearing coefficients are compared to experimental results obtained from a dedicated test rig. Generally, good agreement is observed and minor discrepancies for the damping coefficients are discussed.


Introduction
After five decades of research, compliant foil bearings are now gaining more popularity in the industry than ever before.The current trend is, that foil bearings are being implemented in larger and heavier machinery such as turbo compressors and turbo expanders.Turbo compressors with rotor weights up to 50 kg are available and now widely used.However, mathematical modelling of foil bearings are still associated with significant uncertainties which makes the design of rotor bearing systems difficult and sometimes very costly.Without more accurate mathematical models for prediction of the mechanical behaviour, the technology will continue to be associated with a limited amount of applications and its true potential will not be explored.
Heshmat [1,2] was among the pioneers within this research field.He originally included the flexibility of the foil structure in the Reynolds equation by introducing a linear elastic displacement as function of the fluid film pressure, h c = K(p − p a ).The foil flexibility was based on the analytical expressions given by Walowit and Anno [3].This model is commonly referred to as the simple elastic foundation model (SEFM).The model has later been extended by several authors.Iordanoff [4] developed a more detailed analytical model for the bump foil stiffness which took into account the fixation of the first bump.However, a Corresponding author: jon@stadsen.dkit did not account for the state (stick-slip) of the individual contact points and the resulting interaction between the bumps.Peng and Carpino [5] employed a perturbation method to the SEFM [6] to obtain a steady state solution and equations for the linearised bearing coefficients.Similar work was later performed by Kim [7] and extended with theoretical parametric stability studies.Both used the foil flexibility given by Walowit and Anno [3].
Several authors worked on coupled fluid and structure models as well as including the friction.San Andrés and Kim [8,9] integrated finite element models of the top foil structure into the steady state solution and compared this result against experimental values [10].The bump foils were modelled using the analytical mathematical expressions developed by Iordanoff [4].Carpino et al. [11] developed a structural finite element model of the bump and top foils.Simultaneously, Peng and Carpino [12] investigated the effects of Coulomb friction on the linearised bearing coefficients by means of an equivalent viscous damping coefficient.Their joint effort resulted in the first fully coupled mathematical model [13] with a detailed foil finite element formulation and an equivalent viscous damping for the friction.Their work was of purely theoretical nature.Other authors have later introduced similar models e.g.Heshmat [14] who coupled the structural results obtained by a commercial FE program with the solution of the Reynolds equation for a thrust bearing, and Lee [15] who solved a fully coupled model of a journal foil bearing in the time domain.Lee et al. [16]  The fully coupled models are computationally heavy.As a consequence, accurate equivalent structural models are desirable.Le Lez et al. [17,18] developed equivalent structural models taking into account the Coulomb friction in the contact zones.The theoretical models were compared to experimental results with good agreement and underlined the importance of taking into account the bump interactions and their individual state (stick-slip).Feng and Kaneko [19] developed a similar equivalent foil model coupled with the fluid film equations and compared calculated film heights with experimental data [10].
While mathematical modelling of the gas film itself, using Reynolds equation, has proven to be very accurate compared to experiments [20,21], accurate mathematical modelling of the foil structure is still a main research topic among tribologists, especially when coupling these models with the fluid film equations.A large amount of the work related to the foil structure modelling has been purely theoretical.However the structural models presented by Ku and Heshmat [22][23][24] were compared to experimental data [25,26] showing reasonably good agreement.Theoretical and experimental investigations entirely focused on the bump foil mechanical behaviour were performed by Le Lez et al. [17,18] and Larsen et al. [27].Their theoretical bump foil models were able to accurately predict the steady state hysteresis loops, related to the friction in the sliding points which was previously identified by Ku and Heshmat [28] through isolated bump foil experiments.The experiments of Larsen et al. [27] indicated that the Coulomb friction model is insufficient to properly model the dynamic behaviour of the bump foil structure when subjected to loads at higher frequencies.
Overall experimental investigations of foil journal bearings were performed by several authors e.g.San Andrés et al. [29] who compared predicted unbalance response and critical speeds to experimental results.Dellacorte and Valco [30] experimentally investigated the load carrying capacity of foil bearings and derived a rule of thumb, and Howard [31] investigated the effect of misalignment on the bearing performance.However, the number of authors dealing with the experimental identification of the linearised dynamic stiffness and damping coefficients is limited.Howard et al. [32,33] experimentally identified the dynamic coefficients of a foil journal bearing and investigated the temperature dependency.Matta et al. [34] identified the bearing coefficients of a rigid journal gas bearing.Ertas et al. [35] developed a floating bearing test rig for the identification of foil journal bearing coefficients, and presented experimental results, and San Andrés [36] identified the dynamic coefficients for a hybrid flexural tilting pad gas bearing.
In this paper a theoretical model for obtaining the steady state solution and the linearised bearing coefficients of a journal foil bearing is presented.It takes into account the slope related to the leading edge fixation arrangement of the individual pads and is based on the SEFM which is perturbed in order to obtain a set of first order equations for the linearised bearing coefficients.The zeroth order equation is then coupled with a detailed structural model to obtain a more accurate steady state solution.The structural model has previously been described and compared to experimental results in [27].Results of this model show a significant difference between the static and dynamic foil stiffness, and indicate the need to distinguish between them.The static stiffness should be used when solving for the steady state solution i.e. the zeroth order equation, and the dynamic stiffness should be used when solving the first order equation, which is obtained by assuming small harmonic pressure oscillations.Furthermore, a dedicated test rig has been designed and built to experimentally identify the linearised bearing coefficients of two identical three pad foil bearings supporting a nearly symmetrical rotor.The rotor mass is approximately 21 kg.Bearing coefficients obtained experimentally are presented and compared to the theoretical results of both the SEFM and the presented coupled model.The experimental identification of the bearing coefficients, the inclusion and investigation of the inlet slope and the calculation of the static and dynamic structural stiffness used in the coupled fluid structure model constitute the main original contributions of this work.

Theoretical model
The derivation of theoretical model is based on the SEFM originally presented by Heshmat [1,2] and later extended by several authors [5,13,37,38].Here the derivation of the mathematical model is kept brief but generally follows [39].
For a journal bearing with the nomenclature as given in Figures 1 and 2, the compressible Reynolds equation can be written in vector form [40] as For a simple journal bearing, with the nomenclature as illustrated in Figure 1, the film height can be written as where and with h c being the foil deflection [1,2] and i = 1, 2, . . ., N p the pad number.A perturbation method is employed [6] by assuming that the shaft exhibits small harmonic oscillations around its equilibrium position in the bearing (e x0 , e y0 ).The shaft motion is given by e x = e x0 + Δe x e iωst and e y = e y0 + Δe y e iωst .(5) Assuming the amplitudes to be small Δe x C and Δe y C, a first order Taylor expansion of the pressure can be written as Substituting ( 5), ( 6) into ( 1) and ( 2), discarding second and higher order terms yields, upon separation of variables, the zeroth and first order equations: Zeroth order First order where γ = x, y and f x = cos(θ) and f y = sin(θ) and the film height h 0 is given by where Solving the zeroth order Equation ( 7) for an eccentricity (e x0 , e y0 ) yields the static film pressure p 0 .Then by solving the first order Equation ( 8), with respect to this pressure, the dynamic pressures p x and p y are obtained.Integrating these pressures over the bearing pad areas yields the linearised dynamic stiffness and damping coefficients.The Equations ( 7)-( 9) constitute the SEFM.Solutions of these equations are obtained numerically by use of an FE approach [39].

Foil flexibility
In this paper, the foil flexibility is predicted using two methods.One is to treat the foils linearly having constant flexibility K based on the analytical expression given by Walowit and Anno [3].The other is to model the entire foil structure using a non-linear finite element model.The development of the latter is thoroughly described in [27] and its derivation is only briefly introduced in the following.

Simple elastic foundation model
The steady state foil deflection is given by h c0 = K(p 0 − p a ) and the foil flexibility K can be assumed to be constant, based on the analytical expressions given by Walowit and Anno [3] as or it can be a scalar field K(θ, z, p 0 ) based on a closed form expression e.g.[39] and a mechanical loss factor η can be included by making K complex [5,[37][38][39].Here, (11) is used in complex form to include a loss factor.

Coupled fluid structure model
A more accurate prediction of h c0 can be obtained numerically, using a non-linear FE model.Such a model taking into account the sliding friction in contact points between bump foils and top foils and bump foils and bearing housing is presented and compared to experiments in [27].In Figure 3, load displacement diagrams obtained by (11) and this FE model for a strip of four bumps is compared.
It is clear that the flexibility obtained by ( 11) is overestimated [27].It forms a straight line in the load displacement plot as opposed to the curve predicted by the FE model which forms closed hysteresis loops.From the FE result, it is seen that a monotonic loading will result in a near linear deflection line with a slope k s .At 5, 10, 15 and 20 μm, the monotonic loading is substituted by small load oscillations which are seen to result in significantly higher slopes k d .Due to the significant difference between the stiffness k s and k d it is necessary to distinguish between them.A foil flexibility based on the slope represented by k s should be used when solving for the steady state solution i.e. the zeroth order Equation (7).However, the first order Equation (8) which is obtained by assuming small harmonic pressure oscillations should be solved using a foil flexibility based on the slope represented by k d .
The mesh and boundary conditions of the finite element model are illustrated in Figure 4a.The bump foils are modelled using iso-parametric plain strain elements based on the Green-Lagrange strain measure for large displacements and the top foil is modelled using a simple 2D-plate model with a nodal distribution equal to the fluid film mesh.In this work, the structural non-linear FE model is coupled with the fluid film FE model and solved iteratively in order to obtain a more accurate estimate of the steady state foil deflection h c0 .The coupling of the structure and fluid film is straight forward as the solution of the zeroth order steady state Equation ( 7) follows an iterative Newton-Raphson scheme [39].Between each iteration one may simply solve the structural model for the deflected height h c .However this is computationally demanding so a faster approach, which is used here, is to linearise by a numerical perturbation, each bump stiffness and add these equivalent stiffness k eq to the linear top foil model as shown in Figure 4b.
The equivalent foil structure model is then solved between each Newton-Raphson step until pressure convergence is obtained and the bump equivalent stiffness k eq is updated between each shaft eccentricity step.This approach is valid only if the bump foil behaves linearly.How-ever, as seen in Figure 3 this is the case under monotonic loading.
The coupled fluid-structure model (CFSM) is compared to the deflection obtained by the SEFM using the flexibility given by Walowit and Anno [3] as well as experimental data obtained from a dedicated test rig.

Assumptions
Both the SEFM and the CFSM are based on the assumption that the foil structure deflects evenly over the the axial length of the bearing.For both methods, the pressure applied on the foil structure is taken to be the arithmetic mean pressure along the length of the bearing.According to [1,2] sub-ambient pressures will lift the top foil of the bumps.Here the Gumbel boundary condition is used by simply setting any sub ambient pressures equal to p a .
Solving for a steady state solution (e x0 , e y0 , p 0 , h c0 ) using the CFSM with a coefficient of friction μ f = 0, may lead to several solutions dependant on the initial conditions and the loading path.With the CFSM it is not possible to know the loading path.The loading path for the CFSM is determined by the Newton-Raphson routine applied to solve for the steady state shaft equilibrium position.To overcome the problem of multiple solutions, the coefficient of friction is set to zero when solving for the steady state solution, thus eliminating the load path dependency.

Experimental setup
With the goal of identifying the linearised stiffness and damping coefficients of a foil bearing experimentally, under realistic loading and operation conditions, a dedicated test rig was constructed.The test rig is illustrated in Figure 5.
The test rig consists of a near symmetrical hollow shaft supported by two identical segmented foil journal bearings having 3 pads.Detailed data of the bearings is listed in Table 1.
A set of permanent magnets is press fit mounted in the centre of the shaft, and together with the stator windings this arrangement forms the electrical motor used for spinning the shaft.There are no axial bearings.High accuracy proximity sensors are mounted in horizontal and vertical directions close to each bearing location.In each shaft end a disk is mounted which has a stationary (non-rotating) inner part connected through high precision preloaded ball bearings.This stationary inner part is used in connection with an electromagnetic shaker to excite the shaft.The excitation force is measured using a piezoelectric force transducer mounted directly at the stationary inner part.The assembled shaft mass is approximately 21 kg and the operating speed range is 15 to 30 kRPM.The first free-free natural frequency (bending) of the assembled shaft is calculated to approximately 1050 Hz.This  is more than twice the 1X excitation frequency at maximum speed, hence the shaft can be treated as rigid in the operational speed range.The entire rotor assembly is balanced to below ISO grade G2.5.

Identification procedure
The identification of the linearised bearing coefficients, of the two bearings A and B, is achieved using frequency domain techniques combined with the method of structural joint parameter identification procedure [41,42].Harmonic forcing excitation is applied in the individual excitation points AE, BE of the rotor, while simultaneously obtaining the vibrations in the position sensor locations AP, BP as illustrated in Figure 6.Values of all where [H EP ] is the FRF matrix which is determined experimentally (see Appendices A and B) and {f E } = {f AEx , f AEy , f BEx , f BEy } T and {q P } = {q AP x , q AP y , q BP x , q BP y } T are the complex force and deflection vectors respectively.The FRF matrix (12) can be transformed to relate the forcing and mechanical vibrations to the bearing locations A and B where The transformation matrices [T 1 ] and [T 2 ] are given in Appendix A. The equation of motion of the rotor-bearing system can be written as where the damping and stiffness contributions from the bearings.The dynamic stiffness of the rotor-bearing system can then be stated as 16) Since [H AB ] is obtained experimentally, the stiffness and damping of the bearings can be obtained as −1 (17) where [H r AB ] −1 is the dynamic stiffness matrix of the rotor without bearings which can be obtained theoretically as shown in Appendix A.

If the values of the matrices [H AB ] and [H r
AB ] are exactly known, then from ( 17), the bearing coefficients can be exactly identified too, since no approximations have been introduced.However, both matrices are associated with uncertainties.The matrix [H AB ] is associated with measurement uncertainties and the matrix [H r AB ] is associated with modelling uncertainties.The systematic uncertainties of both have been evaluated using standard statistical methods as described in [43,44].Specifically, a computerized uncertainty analysis is imposed where the uncertainty of all variables, i.e. measurement transducers and geometrical properties like mass, inertia and lengths associated with the identification procedure is taken into account.The identified uncertainties are relatively small, generally below 25%.It should be kept in mind though, that this is under the assumption that the dynamics of the test rig can be precisely described by the mathematical model given in Appendix A i.e. the higher order (nonmodelled) dynamics does not play a role in the frequency range of interest.Any mechanical effects not described by the mathematical model may result in significant uncertainties, higher than what has been estimated.

Results
The experimental results obtained from the test rig are compared to theoretical predictions based on the mathematical model represented by the zeroth and first order Equations ( 7) and (8), respectively.In this comparison, the compliant height h c is based on both the compliance (11) given by Walowit and Anno [3] as well as the structural finite element model illustrated in Figure 4 and derived in [27].

Experimental results
The first four eigenfrequencies of the rotor-bearing system are identified between 90 and 210 Hz for the entire speed range.As mentioned, the test rig is operated below its first bending mode so the corresponding eigenmodes are a barrel and a conical mode with respectively backward and forward precession.Validation of the mathematical model is of particular interest around these modes, as they can potentially become unstable during operation.Therefore, the experimental identification of the bearing coefficients is carried out in the frequency range 25-250 Hz.Validation at the synchronous frequency is desirable too, but identification close to and at the synchronous frequency turned out difficult and even impossible with the current test setup.The shaft vibration amplitudes are below 5 μm 0-p meaning that the mechanical and magnetic run-out becomes significant.The frequency domain based method, described here, cannot compensate for the run-out or mass imbalance.They both results in fundamental vibration components at the synchronous frequency.This is the main reason for not identifying coefficients at synchronous frequency.In order to obtain the synchronous bearing coefficients with the current test rig, a very accurate run-out compensation, requiring additional instrumentation is necessary, and the knowledge of the exact rotor mass imbalance and the vibrations caused by the excitation bearing arrangement should be known.Hence it is decided to focus on the sub-synchronous frequency range around the rigid modes only.In Figure 7, the identified coefficients, at 20 kRPM for bearing B, are illustrated.Results for bearing A are very similar.All cross coupled stiffness and damping coefficients are positive in this sub-synchronous frequency range which is likely a consequence of the segmented bearing geometry.The stiffness coefficients are flat over the frequency range, whereas the damping is found to decrease asymptotically towards low values as the excitation frequency approaches the rotational frequency.The systematic uncertainties are seen to be low for all coefficients except for the direct stiffness coefficients at higher speeds.As mentioned, these uncertainties exclude the effects of "un-modelled" dynamics in the identification procedure.
Coefficients are identified at rotor speeds from 16 kRPM up to 25 kRPM and only minor differences are observed.The focus in this paper is therefore limited to one operating speed only i.e. 20 kRPM.

Theoretical results
By using the SEFM, as given in (7) to (9), with a foil flexibility based on (11), a set of bearing coefficients is obtained.The coefficients are illustrated in Figure 8. Comparing these coefficients to those experimentally obtained, a qualitative similarity is seen.However, the quantitative discrepancies are large.The SEFM significantly under predicts the bearing coefficients.As mentioned in Section 2.2, solving for a steady state solution, using the CFSM with a coefficient of friction μ f = 0, may lead to several solutions.Therefore, the coefficient of friction is set to zero when solving for the steady state solution.The resulting shaft eccentricity ratios predicted by using the two different methods are listed in Table 3.
Though it is not possible to measure the shaft eccentricity ratios experimentally, it is assumed that the results of the CFSM are the most accurate since this model is more detailed.Furthermore, previous comparisons have proven that the bump foil stiffness is far under estimated by the SEFM [27], hence the large difference in the shaft eccentricity ratios as illustrated in Table 3.
Having a steady state solution of better accuracy, the first order equation is solved with a flexibility based on the slope represented by k d , as illustrated in Figure 3, to yield the linearised coefficients.The steady state load on all pads and hence all bumps are then known, so a separate foil calculation can be performed in which the load is applied gradually with μ f = 0 until the steady state load is achieved.At this point, the foil deflection is equal to the results of the CFSM.Then by performing simultaneous load oscillations with an amplitude equal to a fraction of the steady state loads and with μ f = 0.2 [27], the flexibility for each bump is found.The resulting load displacement diagram for each of the bumps under the second bearing pad (the loaded pad) is illustrated in Figure 9. Here, the amplitudes of the individual load oscillations are given as {f b0 }/2, where {f b0 } is the vector containing the radial reaction forces at the top of the bumps (in the contact between bump and top foil).The slope related to the oscillations is found to be approximately a factor of 10 higher than the steady state slope.With the flexibility K, in the first order Equation ( 8), based on this, the bearing coefficients are recalculated.The result is illustrated in Figure 10.Comparing the stiffness coefficients obtained using the CFSM, against the coefficients experimentally obtained, presented in Figure 7, shows good agreement.However, discrepancies are seen when comparing the damping coefficients.Even though a loss factor of η = 0.25 is used when solving the first order equation, the experimentally identified damping coefficients are generally higher than the theoretically predicted.This is especially evident for the direct damping coefficient in the y direction, and generally at lower frequencies.It is important to highlight, that the choice of the structural loss factor η = 0.25 is based partly on a guess and partly on previous investigations [37,45].However, the equivalent structural loss factor is strongly dependant on both frequency and displacement amplitude.The latter is clearly seen in Figure 11 where the load/deflection amplitude is increased, compared to Figure 9.
It is clear that the area covered by the hysteresis loops increases leading to more energy dissipation.This increase in energy dissipation does not follow a linear relationship with the size of deflection amplitude [27].In fact, a small increase in deflection amplitude can lead to a 10 times higher energy dissipation.Choosing a higher structural loss factor η = 0.8 is found to increase the overall damping coefficients without significantly altering the calculated stiffness coefficients.This is illustrated in Figure 12.However, the direct damping coefficient in the horizontal direction is still under predicted.
Parametric studies indicates that the parameter h s , which defines the slope height at the leading edge of each pad, is quite important (see Fig. 3).In the Siemens bearing design, this slope height is originally introduced to avoid rotor to stator contact.The leading edge of the pads is not flexible due to the foil fixation arrangement.Hence the leading edge needs to be detracted somewhat (h s ).In Figure 13, the stiffness and damping coefficients obtained theoretically using the nominal value h s = 50 μm as well as a smaller value h s = 5 μm are compared.
Both the stiffness and damping coefficients are clearly altered by the variations in h s .In fact, discarding the slope in the theoretical model leads to large deviations in the bearing coefficients when compared to experimentally obtained results.

Conclusion
Experimentally obtained bearing coefficients of an air foil bearing are presented and compared to theoretical predictions.The theoretical predictions are based on the simple elastic foundation model as well as an improved model implementing a detailed foil structure finite element formulation [39] coupled with the zeroth order steady state equation.The finite element formulation takes into account large deflections and friction in the sliding contact points.Results of the structural foil model indicate the need to distinguish between the static and dynamic stiffness of the foil structure when solving the zeroth and first order equations respectively.Ignoring this, by using the simple elastic foundation model with a constant foil stiffness, a poor agreement between theoretical predictions and experimental results is observed.In contrast, by distinguishing between the static and dynamic foil stiffness, when solving the zeroth and first order equations respectively, the theoretically obtained coefficients show generally good agreement with the experimental results.However, a certain discrepancy related to the direct damping coefficient perpendicular to the loading direction is observed.This discrepancy cannot be explained by varying the equivalent loss factor alone.The loss factor is regarded constant for the entire foil structure.However, a lightly loaded pad can offer a significantly higher amount of damping than a heavily loaded pad since the sticking phase is less dominant [27].This is also seen as the wider hysteresis loops occurring in the lower left of the load-displacement plot Figure 3.A more accurate method for predicting the damping should be established in future work.
Finally, the influence of the leading edge slope, related to the foil fixation arrangement, is investigated and found to affect the bearing coefficients significantly.Careful modelling of this part of the bearing geometry is therefore important.

A.2 Rotor model matrices
The dynamic stiffness matrix of a rigid rotor without bearings can be written as:

A.3 Experimental FRF matrix
The experimentally obtained FRF matrix [H EP ] are defined as: where each of the individual FRF's f Ej q Pi are obtained experimentally and illustrated in Appendix B.

Appendix B: Experimental FRF
The experimental frequency response functions (FRF's) are obtained by linear 5-300 Hz chirp excitation with a duration of 9.5 s.All signals were simultaneously sampled at a frequency of 1706 Hz.To cancel out measurement noise, each FRF consists of an average of 50 measurements (individual chirps).1024 points were used in the FFT and 512 overlaps when calculating the auto power and cross power spectrum for each FRF.

Fig. 3 .
Fig. 3. Load displacement diagrams obtained using non-linear finite element model and the analytical expression of Walowit and Anno [3].

Fig. 4 .
Fig. 4. (a) First bump segment (at leading edge) of the structural finite element model of top and bump foils with boundary conditions.(b) Equivalent linear model.

Fig. 6 .
Fig. 6.Schematics and position nomenclature of the test shaft.
[M ], [C] and [K] are the mass, damping and stiffness matrices of the rotor alone, and [C b ] and [K b ] are

Fig. 7 .
Fig. 7. Experimentally obtained linear bearing coefficients versus excitation frequency ωs, and the associated uncertainties, experimentally obtained at 20 000 RPM for bearing B.

Fig. 9 .
Fig. 9. Calculated load displacement diagram for each of the 10 bumps under the second pad.Low slope μ f = 0, high slope μ f = 0.2.
coupled Article published by EDP Sciences J.S. Larsen et al.: Mechanics & Industry 16, 106 (2015) a detailed three dimensional structural model with the steady state solution of the Reynolds equation but did not include friction.

Table 1 .
Geometry, material properties and operating conditions of the Siemens air foil bearing.

Table 2 .
Geometry and parameters for the identification.

Table 3 .
Predicted steady state eccentricity ratios using SEFM and CFSM, respectively.