- Email: [email protected]
Improved identification of squeeze-film damper models for aeroengine vibration analysis
ARTICLE IN PRESS Tribology International 43 (2010) 1639–1649
Contents lists available at ScienceDirect
Tribology International journal homepage: www.elsevier.com/locate/triboint
Improved identification of squeeze-film damper models for aeroengine vibration analysis Keir Harvey Groves a,, Philip Bonello b a b
School of Mechanical, Aerospace and Civil Engineering, Desk E5, F-Floor, Pariser Building, University of Manchester, Manchester, UK School of Mechanical, Aerospace and Civil Engineering, C14, Pariser Building, University of Manchester, Manchester, UK
a r t i c l e in fo
abstract
Article history: Received 3 February 2010 Received in revised form 22 March 2010 Accepted 23 March 2010 Available online 30 March 2010
Numerical solution of the Reynolds equation imposes a prohibitive computational cost on the dynamic analysis of practical squeeze film damped turbomachinery. To surmount this problem, the present paper develops the use of Chebyshev polynomial fits to identify finite difference (FD) solution of the incompressible Reynolds equation. The proposed method manipulates the Reynolds equation to allow efficient and accurate identification in the presence of cavitation, the feed-groove, feed-ports, end-plate seals and supply pressure. The ability of Chebyshev polynomials to rapidly reproduce FD routines is demonstrated. The bearing models developed are experimentally proven to give more accurate results than alternative analytical bearing models. & 2010 Elsevier Ltd. All rights reserved.
Keywords: Squeeze film damper Reynolds equation Identification
1. Introduction Squeeze film damper (SFD) bearings are an effective solution to the problem of attenuating vibration and transmitted forces caused by rotor unbalance in gas turbine engines. Their simplicity of construction and consequent robustness mean they are commonplace in modern aircraft gas turbine engines. For example, the engines of a leading aero-engine manufacturer typically have 5–6 SFDs. However, their inclusion requires careful unbalance response calculations that take into account of the SFD’s nonlinearity, to ascertain smooth running [1]. Methods used to analyse the nonlinear response of rotordynamic systems may be classified as either time domain or frequency domain techniques. A time domain method marches the equations of motion forward in time until a steady-state solution is reached. This solution may not necessarily be periodic. This process requires the numerical solution of the coupled nonlinear equations of motion at each of a very large number of time steps. A frequency domain approach like the harmonic balance (HB) method solves nonlinear algebraic equations to extract steady-state solutions that are assumed to be periodic of given fundamental frequency. The relative merits of these two complementary approaches are discussed in [1]. A typical engine model requires consideration of many hundreds of modes, posing
Corresponding author. Tel.: + 44 161 3062605.
E-mail addresses: [email protected], [email protected] (K.H. Groves), [email protected] (P. Bonello). 0301-679X/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.triboint.2010.03.010
prohibitive demands on conventional time/frequency domain methods. The issue of the large number of modes was recently resolved through the development of novel fast time/frequency domain methods [1]. However, like conventional methods, these new methods still require a number of SFD force computations per time step (time domain) or iteration (HB). In the case of HB, one iteration requires the calculation of the SFD forces at each of an adequate number of time points over the period of vibration in order to obtain their Fourier coefficients. Moreover, the iterative process requires a Jacobian matrix that is obtained through repeated calculation of these Fourier coefficients. The SFD model can cripple a time/frequency domain solver unless the SFD forces are rapidly computed. This has lead to a trade-off between the capability/reliability of the SFD model and its computational cost. Della Pietra and Adiletta [2] provide a comprehensive review of SFD modelling techniques, where the above-mentioned tradeoff is evident among most of the research works cited. Advanced SFD models based on the Navier–Stokes equation to account for fluid inertia were proposed in [3,4]. However, in order to enable efficient solution for just one SFD, the system dynamics had to be restricted to circular centred journal orbits or small-amplitude orbits about the bearing centre. These restrictions render such an approach unsuitable for real engines like that considered in [1]. The typical and more practical approach in SFD modelling is based on the incompressible Reynolds lubrication equation to describe the relationship between the pressure in a fluid film and the journal motion. Even then, analytical solution is only achievable if the equation is simplified. Most simplifications are based upon neglecting pressure gradients in one direction, resulting in solutions termed long and short bearing approximations and their
ARTICLE IN PRESS 1640
K.H. Groves, P. Bonello / Tribology International 43 (2010) 1639–1649
Nomenclature aj ,bj ,cj ,dj Ah c Ci d e F h I lp ls L Lg p p0 ps pi,j
grouped constants—see Section 2.1 area of feed-port, m2 clearance of a centred SFD, m Chebyshev polynomial coefficients gap between sealing plate and journal, m journal eccentricity, m unbalance force amplitude, N squeeze film thickness, m moment of inertia about pivot, kg m2 length feed-port, m effective length of end-plate seal, m land length, m groove length, m pressure, Pa outlet pressure, Pa supply pressure, Pa pressure at node i,j in the mesh, Pa
variants (e.g. [5]). The most prominent of these variants involves the combination of the short and long pressure solutions through an ‘end-leakage factor’ l that accounts for the degree of end-sealing [5,6]. This model is popular with industry and was used recently for whole-engine analysis with 5 SFDs [1]. The factor l can be described as an empirical ‘fudge factor’ that in reality can only be related to the parameters of a given bearing when it is hosted in a given experimental setup. Attempts to theoretically relate it to the SFD parameters have proved unsuccessful [5]. To solve the Reynolds equation in its full two-dimensional form, with arbitrary boundary conditions, a numerical scheme like finite difference [4], or finite element (FE) and its variants [7], is necessary. The direct use of such a numerical scheme within a time/frequency domain solver is computationally prohibitive for practical applications. Hence, the recently added SFD routine in NastranTM is based on a one-dimensional approximation of the FD model that is restricted to long and short unsealed SFDs [8]. The full 2-D model can be rapidly deployed within a solver via an identification scheme, which usually involves an interpolation procedure [9,10]. Rodrigues et al. [9] use Chebyshev polynomial fits of a numerical bearing model to relate the three variables e, e_ _ to the two orthogonal output force vectors Q and Q and c R T (Fig. 1). The assumption of a uniform boundary condition at the groove in [9] meant that, in the absence of fluid inertia effects, QR,T were only functions of the aforementioned three variables. The inclusion of feed-ports into a bearing model necessitates the inclusion of c in the input space. Identification of functions with
Fig. 1. Polar coordinate system.
Q R t Ti u U vf w Dx, Dz Cj, Ej, Dj,
e y
m c
o ðÞ^ ðÞ_ ðÞ0
squeeze film force, N radius of SFD journal, m time, s the Chebyshev polynomials fluid velocity in the x direction, m/s2 unbalance, kg m average fluid velocity in feed-port, m/s2 fluid velocity in the z direction, m/s2 node spacing in the x and z directions, m Rj FD matrices non-dimensional journal eccentricity angular location, rad, see Fig. 2 dynamic viscosity, Pa s attitude angle, rad rotational speed, rad/s indicates a non-dimensional property differentiation with respect to time t differentiation with respect to tð ¼ otÞ
more than three variables is possible, yet the computation becomes numerically cumbersome [9]. It follows that the inclusion of feed-ports would not be practical using the technique of Rodrigues et al. Before the identification procedure may take place, the range of each of the input variables must be prescribed. This is simple in the case of e (bounded by radial clearance c) and c (if the feedports are equi-spaced by angle a around the circumference then QR,T are periodic in c, period a). However, selecting a suitable _ is more difficult, since these variables have no range for e_ and c natural bounds. For example in [9] the upper limits of the magnitudes of these two variables were set arbitrarily to 0.7 m/s and 20,000 rad/s, respectively. The judicious prescription of the range of these two input variables is critical to maintaining accuracy of the fit. Widening their range lowers the quality of fit, all other things kept equal. A further consideration is that part of the ‘training’ data required for the calculation of the Chebyshev coefficients involves unrealistic combinations of the input vari_ for which QR-N. This condition ables, namely high e and high e, is unrealistic since QR increasingly acts to reduce e_ as e grows. The inclusion of such unrealistic input combinations skews the identification. The work of Chen et al. [10], although using a crude identification technique based on linear interpolation, surmounts the problems associated with defining the input range by _ for a applying a reduction technique that exchanges e_ and c single variable with limits of 71. Moreover, the ability to reduce the number of input variables by one makes it practical to use Chebyshev polynomials to identify SFD models that necessitate inclusion of c. The major shortcoming of the reduction technique, as applied by Chen et al., is that it cannot accommodate fluid film rupture at a non-atmospheric pressure (e.g. absolute zero) and cannot handle nonzero boundary conditions (e.g. oil supply pressure). The research of the present paper addresses the need for reliable identification of numerical models. It combines the relative strengths of the work in [9,10] and uses two novel techniques to overcome their limitations. Firstly, Chebyshev interpolation, while using the reduction technique, is performed on the pressure. This allows post-interpolation truncation of the pressure to account for cavitation. Secondly, static and dynamic elements of the pressure function are separated. This treatment allows the identification of bearing models that have
ARTICLE IN PRESS K.H. Groves, P. Bonello / Tribology International 43 (2010) 1639–1649
non-axisymmetric and nonzero boundary conditions. This enables their rapid implementation in a full engine analysis, with no question of operation outside the range of identification. The identification method is equally applicable to both the FD and FE methods of solution of the Reynolds equation. The FD scheme is used here. The FD solution is obtained using boundary conditions that are developed to model end-plate seals, feed-grooves and feed-ports. To show the accuracy and speed of the identification, time/frequency domain analysis is performed on a simple rotor system and FD solutions are compared to their equivalent Chebyshev solutions. Comparisons are also made with experimental results and predictions using simple analytical models.
2. Theory The theory will refer to the end-plate-sealed SFD bearing, shown schematically in Fig. 2, which is typical of aeroengine applications. Fluid inertia effects are neglected, as in, [1,5,6,8,10,11] and most of [9]. The instantaneous pressure distribution pðz, yÞ is given by the solution to the incompressible Reynolds equation [12]: 1 @ @p @2 p _ sin y þ e_ cos yÞ þ h3 2 ¼ 12mðec h3 ð1Þ @y R2 @y @z y where h ¼ c þ e cos y. The forces QR,T in Fig. 1 are then given by: " # Z L þ Lg Z 2p QR cos y dz dy ¼ 2RL pt ðz, yÞ ð2Þ QT sin y 0 0 ywhere pt ðz, yÞ is a truncated pressure distribution that accounts for oil film rupture due to cavitation: ( pðz, yÞ if pðz, yÞ 4 pcav ð3Þ pt ðz, yÞ ¼ pcav if pðz, yÞ r pcav The cavitation pressure pcav in the present research is set to 101.325 kPa (absolute zero) [1,13,14]. In this paper, the untruncated pressure distribution p(z,y) is obtained by three different approaches: FD; Chebyshev fitting of the FD-computed p(z,y) and approximate analytical solution of Eq. (1). 2.1. Finite difference By converting partial derivatives into central difference formulas, the Reynolds equation is expressed in FD format: ð2bj 2cj Þpi,j þ cj pi þ 1,j þ cj pi1,j þ ðaj þ bj Þpi,j þ 1 ðaj þ bj Þpi,j1 ¼ dj
1641
To prepare the problem for solution, Eq. (4) must be converted into the following format: Cj pj þ Ej pj1 þ Dj pj þ 1 ¼ Rj ,
j ¼ 1 . . . Nj
ð5Þ
ywhere pj is a column vector of the pressures p1,j , . . . pNi ,j , Cj, Dj and Ej are Ni Ni matrices and Ri is an Ni 1 vector. Eq. (5) may then be solved for the Ni Nj pressure mesh by Castelli’s method [11]. To establish a feed-port boundary condition, a bearing with no feed-groove but a number of feed-ports is first considered. Assuming symmetry about the central cross-section, one can analyse half the axial length of the bearing. Fig. 3 shows a pressure mesh of such an arrangement with three feed-ports. The pressure gradient in the z direction at the centerline points are obtained by assuming a quadratic relationship between the three pressures in the z direction that are closest to the centreline: 3pNi ,j 4pNi 1,j þ pNi 2,j dp ð6Þ ¼ 2ðDzÞ dz Ni ,j Since the bearing is mirrored along the centreline, it is reasonable to assume that this gradient is zero at nodes that are not feedports. Hence, for j not coinciding with a feed-port: 3pNi ,j 4pNi 1,j þ pNi 2,j ¼ 0
ð7Þ
To model the feed-groove an extra land is added onto the existing model. For clarity, FD quantities within the feed-groove have subscript ‘b’ while their equivalents within the main land have subscript ‘a’. The pressure mesh is now ½ðNi Þa þ ðNi Þb Nj , where ðNi Þb is the number of points along half the axial length of the groove. The relationship between the land and groove is established by continuity of flow in the z direction at the boundary between them. Equating flow rates [12] gives: ! ! 1 dp 1 dp 3 ð8Þ ðhj Þa ¼ ðhj Þ3b 12m dz Ni ,j 12m dz 1,j a
b
Eq. (8) is then converted to FD form by substituting three-point approximations of the type of Eq. (6) for the two boundary pressure gradients in Eq. (8) and noting that ðpNi ,j Þa ¼ ðp1,j Þb . The addition of an end-plate seal boundary condition is achieved by assuming a constant pressure gradient within the seal and constant outlet pressure as per [5]. By performing a volume flow balance at the seal entrance it may be shown that: ! pout ðp1,j Þa 3 1 dp 1 ð9Þ d 1,j ðhj Þ3a ¼ ls 12m dz 12m a
ð4Þ
ywhere subscripts i and j refer to discrete locations in the z and y directions, respectively, pi,j ¼ pðzi , yj Þ, aj ¼ ð3=R2 Þh2j ð@h=@yÞj ð1=2DyÞ, _ sin y þ bj ¼ ð1=R2 Þh3j ð1=ðDyÞ2 Þ, cj ¼ h3j =ðDzÞ2 and dj ¼ 12mðec j _e cos yj Þ.
Fig. 2. Longitudinal cross-section of sealed SFD.
Fig. 3. Centreline boundary conditions.
ARTICLE IN PRESS 1642
K.H. Groves, P. Bonello / Tribology International 43 (2010) 1639–1649
Eq. (9) is then converted to FD form by substituting a three-point approximation of the type of Eq. (6) for the boundary pressure gradient on the left hand of Eq. (9). The simplest feed-port boundary condition states that, for the grid-point coinciding with (or closest to) the feed-port, the pressure is fixed and equal to the supply pressure: ðpNi ,j Þb ¼ ps
ð10Þ
This will be referred to as a ‘fixed pressure’ feed-port. A more accurate model that allows for a variable feed-port pressure can be derived by considering the following flow balance (Fig. 4): vf Ah ¼ 2w dy Dxu1 dy Dz þ u2 dy Dz þ
dh Dx Dz dt
2.2. Chebyshev polynomial fitting The Chebyshev polynomials are a family of orthogonal polynomials that are particularly suited to use in function approximation problems. Chebyshev identification of a function with two independent variables is illustrated below. However, it is simple to expand the input space to three variables. Identification in this paper is performed on systems with both two and three variables. A Chebyshev polynomial of degree r is defined as Tr ðxÞ ¼ cosðr arccosðxÞÞ [16], for 1 r x r1 and r ¼ 0,1,2, . . . ,m. Letting ffit ðe, cÞ denote the Chebyshev polynomial fit of the
ð11Þ
ywhere vf is average flow velocity in the feed-port of crosssectional area Ah. Using a pipe flow relationship [15], it is possible to express vf as a function of the variable pressure at the feed-port outlet and the constant supply pressure ps in the pipe leading to the feed-port:
Table 2 Highest degree of polynomials used to perform fit. Pressure function (see Eq. 20)
Variable
Degree of polynomial
for pdyni,j
e
15
q1/2
1 3 4 4
c
1 Ah fps ðpNi ,j Þb g vf ¼ pmlp 8
ð12Þ
for pstati,j
e
c
Substituting Eq. (12) into Eq. (11), integrating the velocity profiles of u1, u2 and w and converting derivatives into numerical format gives the pipe flow boundary condition when position j coincides with a feed-port: ð2bj þ 2cj þ lc ÞðpNi ,j Þb 2cj ðpNi 1,j Þb bj ðpNi ,j þ 1 Þb bj ðpNi ,j1 Þb ¼ lc ps dj
ð13Þ
ywhere lc ¼ 3A2h =ð2pDyDzRlp Þ.
Fig. 6. High eccentricity orbit plot.
0.6 0.4 Fig. 4. Feed-port flow balance boundary condition.
0.2 0
0
0.5
1
1.5
2
2.5
3
0
0.5
1
1.5
2
2.5
3
0.6 Fig. 5. Test rig schematic.
0.4
Table 1 SFD parameters. Journal radius Radial clearance Land length Groove width Groove depth Seal gap Seal length Oil viscosity
80.120 mm 0.254 mm 11.25 mm 2.0 mm 2.0 mm 0.8–0.5 mm 3.25 mm 0.006 Pa s
0.2 0
Fig. 7. Harmonic balance computation speed comparison: (a) finite difference solution—103,111 s; and (b) Chebyshev polynomial solution—438 s.
ARTICLE IN PRESS K.H. Groves, P. Bonello / Tribology International 43 (2010) 1639–1649
function f ðe, cÞ, then: ffit ðe, cÞ ¼
m X n X
~Þ ~ s ðc Cr,s Tr ðeÞT
ð14Þ
r ¼0s¼0
~ are e and where Cr,s are the coefficients to be determined, e~ and c c normalised over the range 71: e~ ¼ 2ðe1=2ðemax þ emin ÞÞ=ðemax emin Þ
ð15Þ
~ . It has been shown [16] that by invoking the ysimilarly for c orthogonality of the polynomials: ( ) m X n .n o X ~ ~ f ðer , cs ÞTr ðe r ÞTs ðc s Þ JTr J2 JTs J2 ð16Þ Cr,s ¼ r ¼0s¼0
ywhere JðÞJ represents the Euclidean norm: 2r þ1 p ~ ¼ cos 2s þ 1 p , c e~ r ¼ cos s mþ1 2 nþ1 2
ð17a; bÞ
Hence, calculation of the coefficients Cr,s requires ‘training’ data f ðer , cs Þ at input sets ðer , cs Þ. In the bearing models developed in Section 2.1, pðz, yÞ is a _ and (if feed-ports are included) c. The variables _ c function of e, e, _ _ c have no clear bounds. A technique is described that reduces e, the input space from four variables to three and allocates these variables finite bounds. Following the work of Chen et al. [10], Eq. (1) is rewritten in _ j: _ is greater than jec two alternative forms. If jej 2 1 @ @p 1 @ p 1 þ h3 h3 ¼ 12mðq1 sin y þ cos yÞ ð18Þ @y @z2 R2 @y
1643
_ =e. _ j is greater than jej: _ If jec _ where p 1 ¼ p=e_ and q1 ¼ ec 2 1 @ @p 2 @ p2 þ h3 h3 ¼ 12mðsin y þ q2 cos yÞ @y @z2 R2 @y
_ and q ¼ e=e _ . Given instantaneous values e, q _ c where p 2 ¼ p=ec 2 1/2 and c, either of Eqs. 18 and 19 can be solved for instantaneous distribution of p 1=2 (the influence of c is felt through the application of a feed-port boundary condition). Hence, p 1=2 are functions of only three variables e, c and q1 or q2. Moreover, q1/2 are bounded between 71. In this research, unlike [9,10], it is the (pre-truncated) gridpoint pressures pi,j that are identified, rather than the forces QR,T. This allows truncation of the grid of interpolated pressures at any arbitrary rupture pressure (as per Eq. (3)) and its subsequent quadrature (as per Eq. (2)) to obtain the forces QR,T. Moreover, the exclusion of cavitation from the polynomial approximation process reduces the complexity of the pressure relationship, thereby reducing the number of coefficients required to model it. The identified grid point pressures are given by: 8 _j _ Z jec < e_ p 1i,j ðe,q1 , cÞ, jej _ _ c , cÞ ¼ ð20Þ pi,j ðe, e, _j _ p ðe,q , cÞ, jej : ec _ o jec 2 2i,j In Eq. (20) it is the p 1=2 functions that are actually fitted as Chebyshev polynomials. It is clear that for each grid point two Chebyshev polynomial fits will be necessary. The coefficients of these fits are computed using, training, data from FD solutions of Eqs. (18) and (19). Although the coefficients of both fits have to be held in memory, only one of these is used for any single force
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.5
1
1.5
2
2.5
3
0
0.5
1
1.5
2
2.5
3
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
ð19Þ
Fig. 8. Short bearing two land model (Eq. (22)): (a) ps ¼239.2 kPa and U^ ¼ 0:128; and (b) ps ¼115.1 kPa and d^ ¼ 0:45.
ARTICLE IN PRESS 1644
K.H. Groves, P. Bonello / Tribology International 43 (2010) 1639–1649
computation. Hence, the force computation time is not significantly increased. When solving Eq. (18) or (19) for p 1=2 in the reduced ðe,q1=2 , cÞ space, it is not possible to accommodate an arbitrary boundary pressure pB since p B would be unknown as a consequence of _ from the input space. Hence, the having removed e_ and c reduction technique is only valid for zero (gauge) supply and outlet pressure. To allow use of arbitrary boundary conditions, the solution pi,j to the Reynolds equation is divided into its ‘dynamic’ and ‘static’ parts: _ Þ þpstat ðe, cÞ _ c, c pi,j ¼ pdyni,j ðe, e, i,j
ð21Þ
The first term of Eq. (21) is given by Eq. (20). It is the solution obtained under the actual vibratory conditions but with supply and outlet pressures set to zero. The second term of Eq. (21) is the solution obtained under the given boundary pressure conditions _ set to zero). This _ c for a static journal orientation of e, c (i.e. e, introduces an extra Chebyshev fit. However, since it is a function of only two variables, it adds only a small time penalty. The polynomial coefficients of the identified FD pressure function of a particular SFD are then deployed to compute the pressure instead of FD in any rotor-dynamic assembly where that bearing is installed. 2.3. Short sealed bearing approximations By dropping the pressure gradient in the y direction in Eq. (1) and applying the end-plate seal boundary condition of Eq. (9), it is
possible to derive analytical solutions to the Reynolds equation for a sealed bearing. Two such solutions are considered. The first model assumes a sealed bearing with two lands of length L each separated by a deep groove of fixed pressure ps. The pressure is then given as: p^ ¼
3 6E L 2 2 1 1 1 p^ s h^ 6EðL=RÞ2 ^þ ^ ^ ^s z z þ þ p z 3 R 3 3 4 2 2 h^ h^ þ d^ ðL=ls Þ
0 where E ¼ e0 cos y þ ec sin y, h^ ¼ 1 þ e cos y, e ¼ e=c, z^ ¼ z=L, d^ ¼ d=c 2 and p^ ¼ p=fmoðR=cÞ g. The second model assumes that the two lands of the bearing effectively form one large land of length 2L and that the bearing is sealed at either end of this land. The justification for such a model is that under tight sealing the feed-groove is somewhat bridged by the higher pressures caused by the sealing. The pressure is then given as: 3
2 ^ dÞ ^ 3 ðls =LÞ p^ ¼ ð6E=h^ ÞðL=RÞ2 ½ðz^ 1=4Þðh=
3.1. Description of the test rig The test rig used to produce the experimental readings presented in this paper was designed to be representative of a small business jet aeroengine [17]. Experimental results presented in this section are taken from reference [17], where
1 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0.5
1
1.5
2
2.5
3
0
0.5
1
1.5
2
2.5
3
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
ð23Þ
3. Testing
0.9
0
ð22Þ
Fig. 9. Experimental results (Ref. [17]): (a) ps ¼ 239.2 kPa and U^ ¼ 0:128; and (b) ps ¼ 115.1 kPa and d^ ¼ 0:45.
ARTICLE IN PRESS K.H. Groves, P. Bonello / Tribology International 43 (2010) 1639–1649
construction details of the rig may also be found. A schematic of the test rig is shown in Fig. 5. It consists of a heavy rotor, driven by an electric motor via a flexible coupling at its left hand end. At the driven end, the rotor runs in a rigidly housed self-aligning bearing. At the right hand end the rotor runs in a rollingelement bearing with a spring-supported housing. An SFD is located approximately at mid-section. Table 1 presents the SFD parameters used in the testing. The lowest undamped critical speed of the rig is 32.4 rev/s [17]. In this mode the rotor performs a conical motion with practically no flexure. Since the operating speed range is well below the first flexure critical speed (218 rev/s [17]) the equations of motion can be written as: kx ðl2 =l1 Þ2 xSFD þ Qx þðl3 =l1 ÞFx ðtÞ ¼ ðI=l21 Þx€ SFD ky ðl2 =l1 Þ2 ySFD þQy þðl3 =l1 ÞFy ðtÞ ¼ ðI=l21 Þy€ SFD
ð24Þ
where I ¼ 23:144 kg m2 . In this rig, the polar moment of inertia of the rotor is sufficiently small for gyroscopic effects to be negligible over the operational range of the rotational speed. Such an assumption was also made in the predictions of [17]. The above equations were solved using time/frequency domain methods to obtain the response for a given unbalance U, quantified here by the non-dimensional parameter U^ ¼ Ul1 l3 =ðIcÞ. All calculations in the following section were performed in Matlabs on a standard desktop pc with a 2.1 GHz dual core processor.
1645
3.2. Chebyshev vs FD To test the effectiveness of the proposed identification technique, results achieved using FD bearing models were tested against equivalent Chebyshev models. To efficiently identify the SFD over the input space, differing degrees of Chebyshev polynomial T (Eq. (14)) were required for each of the input variables. The highest degree of the Chebyshev polynomial used for each variable in the function is displayed in Table 2. It was discovered that the polynomials Tðq1=2 Þ were never above a degree of one; if a higher degree was forced into the fit, the identification process invariably returned zero values for the corresponding coefficients. This unexpected benefit of using the reduction technique with Chebyshev polynomials significantly enhances computation speeds both when executing the polynomials and achieving the coefficients. The one-off computation time for obtaining the coefficients corresponding to degrees in Table 2 was 85.7 s for identifying an FD pressure grid of ðNi Þa ¼ 31, ðNi Þb ¼ 7 and Nj ¼ 97. Fig. 6 shows an orbit from a time domain solution obtained using Matlab’s ode23s&. An FD bearing model was used for the SFD that incorporated end-plate seals, feed-groove and three equi-spaced pipe-flow feed-ports (Eq. (13)). In this example the spring support has been removed from the test rig model, thereby forcing the identified bearing model to operate in a highly nonlinear regime. Eqs. (24) were modified by removing the spring force term and including gravity in the second equation. Two polynomial fits were used to produce Fig. 6, as including
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.5
1
1.5
2
2.5
3
0
0.5
1
1.5
2
2.5
3
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
Fig. 10. FD–Chebyshev, two land deep groove end-plate sealed model: (a) ps ¼ 239.2 kPa and U^ ¼ 0:128; and (b) ps ¼ 115.1 kPa and d^ ¼ 0:45.
ARTICLE IN PRESS 1646
K.H. Groves, P. Bonello / Tribology International 43 (2010) 1639–1649
extremely high eccentricity degrades the fit at lower eccentricity and requires more polynomial coefficients. The fits were performed over 0 r e o 0:95 and 0:95 r e o 0:99. Splitting the fit into two meant the system needed to hold twice as many coefficients in memory; however, this caused a negligible time penalty. The plot is of a transient 0.1 s run from default initial conditions at 60 rev/s with an applied non-dimensional unbalance U^ ¼ 0:256, supply pressure of 239:2 kPa and seal gap of 0.45c. The respective orbit, computed using Chebyshev polynomial fits, is not displayed since the two plots are indistinguishable. The key difference between the two orbits is the execution time. The FD model took 3706 s to run, while the Chebyshev equivalent required only 137 s i.e. 3.7% of the time. Fig. 7 shows plots of amplitude against speed for the springsupported rig with a centralized SFD. The response was solved by harmonic balance (HB) using the arc-length continuation technique developed in [13]. The bearing tested had end-plate seals, feed-groove and three equi-spaced fixed pressure feedports (Eq. (10)). The supply pressure was set at 239.2 kPa with a seal gap of 0.45c and unbalance of U^ ¼ 0:128. It is clear that the use of Chebyshev polynomials is extremely well suited to the approximate solution method of HB. The smoothing effect of the identification technique clearly allows rapid, stable operation of the HB code. The FD solution points can be seen to cluster, reducing arc length to maintain accuracy and therefore the process progresses extremely slowly. The time saving is excellent; the FD solution took 103,111 s to run, while the Chebyshev model only required 438 s, i.e. 0.42% of the computation time.
3.3. FD models vs analytical models The previous section verified that Chebyshev identification provides a reliable means of implementing FD SFD models. Chebyshev will be henceforth be used to implement the bearing models of Section 2.1. The unbalance response of the test rig obtained with these ‘FD-Chebyshev’ models are compared in this section to the results obtained using the analytical bearing models of Section 2.3. Fig. 8 shows the response obtained with the analytical two land end-plate sealed model (Eq. (22)). These graphs exhibit a ‘spring-softening’ effect (i.e. a lowering of the critical speed). This effect was observed as being due to a combination of tight sealing and high supply pressure. The cause of this effect is the second term of Eq. (22), which contains the effect of the supply pressure ^ This and sealing and varies with y due to the presence of E and h. variation means that it does not tend to cancel out when integrated to give the forces, unlike the third term of Eq. (22). When ps is low or d^ is large, the second term of Eq. (22) becomes small and the softening virtually disappears, as evident from Fig. 8(a). The softening effect was not observed in the measurements of Fig. 9. Nor was it observed in predictions obtained for this test rig in [17] using the ‘l theory’, that assumes a pressure distribution of the form [5,6]: pðz, yÞ ¼ pshort,unsealed ðz, yÞ þ lplong ðyÞð0:5z^ Þ
ywhere l is the end-leakage factor. The first term of Eq. (25) is ^ Eq. (22) but with the seal gap ratio d-1 (i.e. unsealed SFD). The
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.5
1
1.5
2
2.5
3
0.5
1
1.5
2
2.5
3
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
ð25Þ
Fig. 11. Short bearing one land model (Eq. (23)): (a) ps ¼239.2 kPa and U^ ¼ 0:128; and (b) ps ¼115.1 kPa and d^ ¼ 0:45.
ARTICLE IN PRESS K.H. Groves, P. Bonello / Tribology International 43 (2010) 1639–1649
second term of Eq. (25) contains the ‘long bearing’ solution to the Reynolds equation, obtained by dropping its axial pressure gradient term (second term on LHS of Eq. (1)). It is noted that the model of Eq. (22) was obtained by dropping the circumferential pressure gradient term from the Reynolds equation. However, it was ascertained that this had nothing to do with the softening behaviour of Fig. 8. Indeed, Eq. (22) was verified to be computationally sound through FD. Fig. 10 shows the amplitude-speed plot for an end-plate sealed SFD model with a central deep groove of fixed pressure ps computed by FD-Chebyshev. This plot is almost indistinguishable from Fig. 8, despite the fact that the results of Fig. 10 are based on the full (two-dimensional) Reynolds equation. This means that the omission of the circumferential pressure gradient to solve the Reynolds equation was verified to have negligible effect for an end-plate sealed damper whose L/R is similar to that in [5,6]. This fact appears to contradict the l-factor theory of Eq. (25). However, comparing Eq. (25) with the (FD-verified) Eq. (22), it is observed that both expressions involve the addition of a circumferentially varying pressure distribution to the short, unsealed SFD expression. In the case of Eq. (22), this additional distribution is its second term. Since this term has no direct equivalence to the second term of Eq. (25), it is proved that it is not possible to find a theoretical link between l and the SFD parameters. This explains the unsuccessful attempt in [5]. Fig. 11 shows the results obtained using the analytical one land end-plate sealed model (Eq. (23)). The softening effect is absent since there is no lubricant supply included in the model. Comparing with the experimental results of Fig. 9, the results of
1647
Fig. 11 are clearly over-damped. This shows that the two lands do not quite act as one under tight sealing. On the other hand, the previous results of Figs. 8 and 10, based on two independent lands and neglecting any damping contributions from the groove, are under-damped. Fig. 12 presents the results obtained using an FD sealed bearing model that incorporates both axial and circumferential pressure gradients, the effect of the feed-groove but without any feed-ports. The results show marked improvement upon the previous models. Damping in the system is increased when compared to Fig. 8 or 10, while being significantly less than in the one-land model of Fig. 11, therefore proving the importance of feed-groove effects. The results of Fig. 12 show no softening since no feed-ports were included. The effect of adding three equi-spaced feed-ports of fixed pressure of ps (Eq. (10)) is shown in Fig. 13. It is clear that the artificial softening effect at tight sealing has reappeared and at loose sealing the critical speed shifts above the undamped critical. These results confirm that the assumption that a feed-port pressure is always fixed, irrespective of the journal dynamics, is unrealistic and results in artificial shifting of the critical speeds. The results presented in Fig. 14 use the pipe flow boundary condition (Eq. (13)), with a feed-port of length of lp ¼15 mm and area of Ah ¼0.7 mm2. It is apparent that modifying the lubricant supply boundary condition has almost removed the critical speed shifting effect on the predictions. It is observed that the results of Fig. 14 show little improvement on the results of Fig. 12. Moreover, comparing the curve for
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.5
1
1.5
2
2.5
3
0
0.5
1
1.5
2
2.5
3
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
Fig. 12. FD-Chebyshev model with seal and groove, no feed-ports: (a) ps ¼ 239.2 kPa and U^ ¼ 0:128; and (b) ps ¼115.1 kPa and d^ ¼ 0:45.
ARTICLE IN PRESS 1648
K.H. Groves, P. Bonello / Tribology International 43 (2010) 1639–1649
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.5
1
1.5
2
2.5
3
0
0.5
1
1.5
2
2.5
3
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
Fig. 13. FD-Chebyshev model with seal, feed-groove and fixed pressure feed-ports: (a) ps ¼239.2 kPa and U^ ¼ 0:128; and (b) ps ¼115.1 kPa and d^ ¼ 0:45.
d^ ¼ 0:45, U^ ¼ 0:128 in Fig. 14(a) with the corresponding one in Fig. 14(b), it is clear that increasing the supply pressure did not improve the predicted damping capacity of the SFD. The reason for this is that, for all the simulations, the SFD pressure was invariably above the chosen cut-off pressure (absolute zero), unlike the experimental observations (Fig. 9). The presence of supply pressure only influences the damping effect if there is cavitation present in the system. If a more demanding test case were used (i.e. one that forced the bearing models to cavitate), the positive effect of supply pressure upon the damping characteristics would be observed in the predictions. A key benefit of carrying out pressure truncation post-identification is that the simulations could be redone at a variable cut-off pressure, as in [5,6] where experimental cut-off pressures were used.
4. Conclusions This paper considered the modelling of end-plate sealed SFD bearings with groove and feed-ports. Finite difference models based on the incompressible Reynolds equation were first developed. A novel identification technique was then presented that provided a means for the FD models to be rapidly deployed. No assumption was made about the range of the input variables, the cavitation pressure or the axial symmetry of the pressure boundary conditions. Rather, arbitrary boundary pressure conditions may be applied at any location in the film while still restricting the number of input variables to a manageable number
(three). The Chebyshev technique has been proven to be particularly suited to use in harmonic balance solution schemes, where staggering time savings were observed. This work opens the possibility of using bearing models of increased accuracy within multi-SFD full aeroengine models. A highly instructive study was made on a simple rig in which the unbalance response was calculated using identified FD SFD models and simple analytical SFD models. It has been demonstrated that the effect of omitting the circumferential pressure gradient to solve the Reynolds equation is negligible for a short SFD, even if it has substantial sealing. This proves that it is not possible to theoretically relate the end-leakage factor (l) to the SFD parameters. The common usage of a fixed pressure condition at a groove or feed-port was shown to be inadequate since it resulted in artificial shifting of the critical speed as the sealing was tightened. The use of an improved feed-port boundary condition virtually eliminated the problem and gave superior agreement with the experimental observations. Despite the improved accuracy of the SFD model, the level of damping was somewhat underestimated relative to the experimental results. This could be partly attributed to the inherent degree of uncertainty in the experimental parameters, e.g. the seal gap dimension, which is only a fraction of the radial clearance for tight sealing. Another possible reason for the discrepancy is the simplicity of the groove model, which neglected inertia effects. Current efforts are being directed towards seeking a correction for fluid inertia while preserving the advantages of the presented identification scheme.
ARTICLE IN PRESS K.H. Groves, P. Bonello / Tribology International 43 (2010) 1639–1649
1649
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.5
1
1.5
2
2.5
3
0
0.5
1
1.5
2
2.5
3
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
Fig. 14. FD-Chebyshev model with seal, feed-groove and pipe flow feed-ports: (a) ps ¼239.2 kPa and U^ ¼ 0:128; and (b) ps ¼ 115.1 kPa and d^ ¼ 0:45. tx
References [1] Bonello P, Hai PM. Computational studies of the unbalance response of a whole aero-engine model with squeeze-film bearings. Trans ASME J Eng Gas Turbines Power 2010;132(3):032504 [7pp.]. [2] Pietra LD, Adiletta G. The squeeze film damper over four decades of investigations. Part I: characteristics and operating features. Shock Vibr Dig 2002;34:3–26. [3] San Andres LA, Vance JM. Effect of fluid inertia on squeeze-film damper forces for small-amplitude circular-centered motions. ASLE Trans 1987;30:63–8. [4] Zhang JX, Roberts JB. Solutions for the combined motion of finite length squeeze film dampers around the bearing center. ASME J Tribol 1996;118:617. [5] Dede MM, Dogan M, Holmes R. The damping capacity of a sealed squeeze film bearing. ASME J Tribol 1985;107:411–8. [6] Holmes R, Dogan M. The performance of a sealed squeeze-film bearing in a flexible support structure. Proc IMechE Part C 1985;199:1–9. [7] Jiangang Y, Rui G, Yongwei T. Hybrid radial basis function/finite element modelling of journal bearing. Tribol Int 2008;41:1169–75. [8] MSC. Nastran release guide. MSC. Software corporation; 2005. p. 150–6.
[9] Rodrigues FA, Thouverez F, Gibert C, Jezequel L. Chebyshev polynomials fits for efficient analysis of finite length squeeze film damped rotors. ASME J Eng Gas Turbines Power 2003;125:175–83. [10] Chen Z, Jiao Y, Xia S, Huang W, Zhang Z. An efficient calculation method of nonlinear fluid film forces in journal bearing. STLE Tribol Trans 2002;45:324–9. [11] Castelli V, Shapiro W. Improved method for numerical solutions of the general incompressible fluid film lubrication problem. ASME J Lubr Technol 1967;89:211–8. [12] Cameron A. Basic lubrication theory. London: Ellis Horwood; 1971. [13] Bonello P, Brennan MJ, Holmes R. Non-linear modelling of rotor dynamic systems with squeeze film dampers—an efficient integrated approach. J Sound Vibr 2002;249:743–73. [14] Feng NS, Hahn EJ. Effects of gas entrainment on squeeze film damper performance. ASME J Tribol 1987;109:149–54. [15] Munson BR, Young DF, Okiishi TH. Fundamentals of fluid mechanics. New York: John Wiley and Sons; 1998. [16] Dahlquist G, Bjorck A. Numerical methods. New York: Prentice-Hall; 1974. [17] Levesley MC, Holmes R. Experimental investigation into the vibration response of an aero-engine rotor-damper assembly. Proc IMechE Part G 1994;208:52–66.
Contents lists available at ScienceDirect
Tribology International journal homepage: www.elsevier.com/locate/triboint
Improved identification of squeeze-film damper models for aeroengine vibration analysis Keir Harvey Groves a,, Philip Bonello b a b
School of Mechanical, Aerospace and Civil Engineering, Desk E5, F-Floor, Pariser Building, University of Manchester, Manchester, UK School of Mechanical, Aerospace and Civil Engineering, C14, Pariser Building, University of Manchester, Manchester, UK
a r t i c l e in fo
abstract
Article history: Received 3 February 2010 Received in revised form 22 March 2010 Accepted 23 March 2010 Available online 30 March 2010
Numerical solution of the Reynolds equation imposes a prohibitive computational cost on the dynamic analysis of practical squeeze film damped turbomachinery. To surmount this problem, the present paper develops the use of Chebyshev polynomial fits to identify finite difference (FD) solution of the incompressible Reynolds equation. The proposed method manipulates the Reynolds equation to allow efficient and accurate identification in the presence of cavitation, the feed-groove, feed-ports, end-plate seals and supply pressure. The ability of Chebyshev polynomials to rapidly reproduce FD routines is demonstrated. The bearing models developed are experimentally proven to give more accurate results than alternative analytical bearing models. & 2010 Elsevier Ltd. All rights reserved.
Keywords: Squeeze film damper Reynolds equation Identification
1. Introduction Squeeze film damper (SFD) bearings are an effective solution to the problem of attenuating vibration and transmitted forces caused by rotor unbalance in gas turbine engines. Their simplicity of construction and consequent robustness mean they are commonplace in modern aircraft gas turbine engines. For example, the engines of a leading aero-engine manufacturer typically have 5–6 SFDs. However, their inclusion requires careful unbalance response calculations that take into account of the SFD’s nonlinearity, to ascertain smooth running [1]. Methods used to analyse the nonlinear response of rotordynamic systems may be classified as either time domain or frequency domain techniques. A time domain method marches the equations of motion forward in time until a steady-state solution is reached. This solution may not necessarily be periodic. This process requires the numerical solution of the coupled nonlinear equations of motion at each of a very large number of time steps. A frequency domain approach like the harmonic balance (HB) method solves nonlinear algebraic equations to extract steady-state solutions that are assumed to be periodic of given fundamental frequency. The relative merits of these two complementary approaches are discussed in [1]. A typical engine model requires consideration of many hundreds of modes, posing
Corresponding author. Tel.: + 44 161 3062605.
E-mail addresses: [email protected], [email protected] (K.H. Groves), [email protected] (P. Bonello). 0301-679X/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.triboint.2010.03.010
prohibitive demands on conventional time/frequency domain methods. The issue of the large number of modes was recently resolved through the development of novel fast time/frequency domain methods [1]. However, like conventional methods, these new methods still require a number of SFD force computations per time step (time domain) or iteration (HB). In the case of HB, one iteration requires the calculation of the SFD forces at each of an adequate number of time points over the period of vibration in order to obtain their Fourier coefficients. Moreover, the iterative process requires a Jacobian matrix that is obtained through repeated calculation of these Fourier coefficients. The SFD model can cripple a time/frequency domain solver unless the SFD forces are rapidly computed. This has lead to a trade-off between the capability/reliability of the SFD model and its computational cost. Della Pietra and Adiletta [2] provide a comprehensive review of SFD modelling techniques, where the above-mentioned tradeoff is evident among most of the research works cited. Advanced SFD models based on the Navier–Stokes equation to account for fluid inertia were proposed in [3,4]. However, in order to enable efficient solution for just one SFD, the system dynamics had to be restricted to circular centred journal orbits or small-amplitude orbits about the bearing centre. These restrictions render such an approach unsuitable for real engines like that considered in [1]. The typical and more practical approach in SFD modelling is based on the incompressible Reynolds lubrication equation to describe the relationship between the pressure in a fluid film and the journal motion. Even then, analytical solution is only achievable if the equation is simplified. Most simplifications are based upon neglecting pressure gradients in one direction, resulting in solutions termed long and short bearing approximations and their
ARTICLE IN PRESS 1640
K.H. Groves, P. Bonello / Tribology International 43 (2010) 1639–1649
Nomenclature aj ,bj ,cj ,dj Ah c Ci d e F h I lp ls L Lg p p0 ps pi,j
grouped constants—see Section 2.1 area of feed-port, m2 clearance of a centred SFD, m Chebyshev polynomial coefficients gap between sealing plate and journal, m journal eccentricity, m unbalance force amplitude, N squeeze film thickness, m moment of inertia about pivot, kg m2 length feed-port, m effective length of end-plate seal, m land length, m groove length, m pressure, Pa outlet pressure, Pa supply pressure, Pa pressure at node i,j in the mesh, Pa
variants (e.g. [5]). The most prominent of these variants involves the combination of the short and long pressure solutions through an ‘end-leakage factor’ l that accounts for the degree of end-sealing [5,6]. This model is popular with industry and was used recently for whole-engine analysis with 5 SFDs [1]. The factor l can be described as an empirical ‘fudge factor’ that in reality can only be related to the parameters of a given bearing when it is hosted in a given experimental setup. Attempts to theoretically relate it to the SFD parameters have proved unsuccessful [5]. To solve the Reynolds equation in its full two-dimensional form, with arbitrary boundary conditions, a numerical scheme like finite difference [4], or finite element (FE) and its variants [7], is necessary. The direct use of such a numerical scheme within a time/frequency domain solver is computationally prohibitive for practical applications. Hence, the recently added SFD routine in NastranTM is based on a one-dimensional approximation of the FD model that is restricted to long and short unsealed SFDs [8]. The full 2-D model can be rapidly deployed within a solver via an identification scheme, which usually involves an interpolation procedure [9,10]. Rodrigues et al. [9] use Chebyshev polynomial fits of a numerical bearing model to relate the three variables e, e_ _ to the two orthogonal output force vectors Q and Q and c R T (Fig. 1). The assumption of a uniform boundary condition at the groove in [9] meant that, in the absence of fluid inertia effects, QR,T were only functions of the aforementioned three variables. The inclusion of feed-ports into a bearing model necessitates the inclusion of c in the input space. Identification of functions with
Fig. 1. Polar coordinate system.
Q R t Ti u U vf w Dx, Dz Cj, Ej, Dj,
e y
m c
o ðÞ^ ðÞ_ ðÞ0
squeeze film force, N radius of SFD journal, m time, s the Chebyshev polynomials fluid velocity in the x direction, m/s2 unbalance, kg m average fluid velocity in feed-port, m/s2 fluid velocity in the z direction, m/s2 node spacing in the x and z directions, m Rj FD matrices non-dimensional journal eccentricity angular location, rad, see Fig. 2 dynamic viscosity, Pa s attitude angle, rad rotational speed, rad/s indicates a non-dimensional property differentiation with respect to time t differentiation with respect to tð ¼ otÞ
more than three variables is possible, yet the computation becomes numerically cumbersome [9]. It follows that the inclusion of feed-ports would not be practical using the technique of Rodrigues et al. Before the identification procedure may take place, the range of each of the input variables must be prescribed. This is simple in the case of e (bounded by radial clearance c) and c (if the feedports are equi-spaced by angle a around the circumference then QR,T are periodic in c, period a). However, selecting a suitable _ is more difficult, since these variables have no range for e_ and c natural bounds. For example in [9] the upper limits of the magnitudes of these two variables were set arbitrarily to 0.7 m/s and 20,000 rad/s, respectively. The judicious prescription of the range of these two input variables is critical to maintaining accuracy of the fit. Widening their range lowers the quality of fit, all other things kept equal. A further consideration is that part of the ‘training’ data required for the calculation of the Chebyshev coefficients involves unrealistic combinations of the input vari_ for which QR-N. This condition ables, namely high e and high e, is unrealistic since QR increasingly acts to reduce e_ as e grows. The inclusion of such unrealistic input combinations skews the identification. The work of Chen et al. [10], although using a crude identification technique based on linear interpolation, surmounts the problems associated with defining the input range by _ for a applying a reduction technique that exchanges e_ and c single variable with limits of 71. Moreover, the ability to reduce the number of input variables by one makes it practical to use Chebyshev polynomials to identify SFD models that necessitate inclusion of c. The major shortcoming of the reduction technique, as applied by Chen et al., is that it cannot accommodate fluid film rupture at a non-atmospheric pressure (e.g. absolute zero) and cannot handle nonzero boundary conditions (e.g. oil supply pressure). The research of the present paper addresses the need for reliable identification of numerical models. It combines the relative strengths of the work in [9,10] and uses two novel techniques to overcome their limitations. Firstly, Chebyshev interpolation, while using the reduction technique, is performed on the pressure. This allows post-interpolation truncation of the pressure to account for cavitation. Secondly, static and dynamic elements of the pressure function are separated. This treatment allows the identification of bearing models that have
ARTICLE IN PRESS K.H. Groves, P. Bonello / Tribology International 43 (2010) 1639–1649
non-axisymmetric and nonzero boundary conditions. This enables their rapid implementation in a full engine analysis, with no question of operation outside the range of identification. The identification method is equally applicable to both the FD and FE methods of solution of the Reynolds equation. The FD scheme is used here. The FD solution is obtained using boundary conditions that are developed to model end-plate seals, feed-grooves and feed-ports. To show the accuracy and speed of the identification, time/frequency domain analysis is performed on a simple rotor system and FD solutions are compared to their equivalent Chebyshev solutions. Comparisons are also made with experimental results and predictions using simple analytical models.
2. Theory The theory will refer to the end-plate-sealed SFD bearing, shown schematically in Fig. 2, which is typical of aeroengine applications. Fluid inertia effects are neglected, as in, [1,5,6,8,10,11] and most of [9]. The instantaneous pressure distribution pðz, yÞ is given by the solution to the incompressible Reynolds equation [12]: 1 @ @p @2 p _ sin y þ e_ cos yÞ þ h3 2 ¼ 12mðec h3 ð1Þ @y R2 @y @z y where h ¼ c þ e cos y. The forces QR,T in Fig. 1 are then given by: " # Z L þ Lg Z 2p QR cos y dz dy ¼ 2RL pt ðz, yÞ ð2Þ QT sin y 0 0 ywhere pt ðz, yÞ is a truncated pressure distribution that accounts for oil film rupture due to cavitation: ( pðz, yÞ if pðz, yÞ 4 pcav ð3Þ pt ðz, yÞ ¼ pcav if pðz, yÞ r pcav The cavitation pressure pcav in the present research is set to 101.325 kPa (absolute zero) [1,13,14]. In this paper, the untruncated pressure distribution p(z,y) is obtained by three different approaches: FD; Chebyshev fitting of the FD-computed p(z,y) and approximate analytical solution of Eq. (1). 2.1. Finite difference By converting partial derivatives into central difference formulas, the Reynolds equation is expressed in FD format: ð2bj 2cj Þpi,j þ cj pi þ 1,j þ cj pi1,j þ ðaj þ bj Þpi,j þ 1 ðaj þ bj Þpi,j1 ¼ dj
1641
To prepare the problem for solution, Eq. (4) must be converted into the following format: Cj pj þ Ej pj1 þ Dj pj þ 1 ¼ Rj ,
j ¼ 1 . . . Nj
ð5Þ
ywhere pj is a column vector of the pressures p1,j , . . . pNi ,j , Cj, Dj and Ej are Ni Ni matrices and Ri is an Ni 1 vector. Eq. (5) may then be solved for the Ni Nj pressure mesh by Castelli’s method [11]. To establish a feed-port boundary condition, a bearing with no feed-groove but a number of feed-ports is first considered. Assuming symmetry about the central cross-section, one can analyse half the axial length of the bearing. Fig. 3 shows a pressure mesh of such an arrangement with three feed-ports. The pressure gradient in the z direction at the centerline points are obtained by assuming a quadratic relationship between the three pressures in the z direction that are closest to the centreline: 3pNi ,j 4pNi 1,j þ pNi 2,j dp ð6Þ ¼ 2ðDzÞ dz Ni ,j Since the bearing is mirrored along the centreline, it is reasonable to assume that this gradient is zero at nodes that are not feedports. Hence, for j not coinciding with a feed-port: 3pNi ,j 4pNi 1,j þ pNi 2,j ¼ 0
ð7Þ
To model the feed-groove an extra land is added onto the existing model. For clarity, FD quantities within the feed-groove have subscript ‘b’ while their equivalents within the main land have subscript ‘a’. The pressure mesh is now ½ðNi Þa þ ðNi Þb Nj , where ðNi Þb is the number of points along half the axial length of the groove. The relationship between the land and groove is established by continuity of flow in the z direction at the boundary between them. Equating flow rates [12] gives: ! ! 1 dp 1 dp 3 ð8Þ ðhj Þa ¼ ðhj Þ3b 12m dz Ni ,j 12m dz 1,j a
b
Eq. (8) is then converted to FD form by substituting three-point approximations of the type of Eq. (6) for the two boundary pressure gradients in Eq. (8) and noting that ðpNi ,j Þa ¼ ðp1,j Þb . The addition of an end-plate seal boundary condition is achieved by assuming a constant pressure gradient within the seal and constant outlet pressure as per [5]. By performing a volume flow balance at the seal entrance it may be shown that: ! pout ðp1,j Þa 3 1 dp 1 ð9Þ d 1,j ðhj Þ3a ¼ ls 12m dz 12m a
ð4Þ
ywhere subscripts i and j refer to discrete locations in the z and y directions, respectively, pi,j ¼ pðzi , yj Þ, aj ¼ ð3=R2 Þh2j ð@h=@yÞj ð1=2DyÞ, _ sin y þ bj ¼ ð1=R2 Þh3j ð1=ðDyÞ2 Þ, cj ¼ h3j =ðDzÞ2 and dj ¼ 12mðec j _e cos yj Þ.
Fig. 2. Longitudinal cross-section of sealed SFD.
Fig. 3. Centreline boundary conditions.
ARTICLE IN PRESS 1642
K.H. Groves, P. Bonello / Tribology International 43 (2010) 1639–1649
Eq. (9) is then converted to FD form by substituting a three-point approximation of the type of Eq. (6) for the boundary pressure gradient on the left hand of Eq. (9). The simplest feed-port boundary condition states that, for the grid-point coinciding with (or closest to) the feed-port, the pressure is fixed and equal to the supply pressure: ðpNi ,j Þb ¼ ps
ð10Þ
This will be referred to as a ‘fixed pressure’ feed-port. A more accurate model that allows for a variable feed-port pressure can be derived by considering the following flow balance (Fig. 4): vf Ah ¼ 2w dy Dxu1 dy Dz þ u2 dy Dz þ
dh Dx Dz dt
2.2. Chebyshev polynomial fitting The Chebyshev polynomials are a family of orthogonal polynomials that are particularly suited to use in function approximation problems. Chebyshev identification of a function with two independent variables is illustrated below. However, it is simple to expand the input space to three variables. Identification in this paper is performed on systems with both two and three variables. A Chebyshev polynomial of degree r is defined as Tr ðxÞ ¼ cosðr arccosðxÞÞ [16], for 1 r x r1 and r ¼ 0,1,2, . . . ,m. Letting ffit ðe, cÞ denote the Chebyshev polynomial fit of the
ð11Þ
ywhere vf is average flow velocity in the feed-port of crosssectional area Ah. Using a pipe flow relationship [15], it is possible to express vf as a function of the variable pressure at the feed-port outlet and the constant supply pressure ps in the pipe leading to the feed-port:
Table 2 Highest degree of polynomials used to perform fit. Pressure function (see Eq. 20)
Variable
Degree of polynomial
for pdyni,j
e
15
q1/2
1 3 4 4
c
1 Ah fps ðpNi ,j Þb g vf ¼ pmlp 8
ð12Þ
for pstati,j
e
c
Substituting Eq. (12) into Eq. (11), integrating the velocity profiles of u1, u2 and w and converting derivatives into numerical format gives the pipe flow boundary condition when position j coincides with a feed-port: ð2bj þ 2cj þ lc ÞðpNi ,j Þb 2cj ðpNi 1,j Þb bj ðpNi ,j þ 1 Þb bj ðpNi ,j1 Þb ¼ lc ps dj
ð13Þ
ywhere lc ¼ 3A2h =ð2pDyDzRlp Þ.
Fig. 6. High eccentricity orbit plot.
0.6 0.4 Fig. 4. Feed-port flow balance boundary condition.
0.2 0
0
0.5
1
1.5
2
2.5
3
0
0.5
1
1.5
2
2.5
3
0.6 Fig. 5. Test rig schematic.
0.4
Table 1 SFD parameters. Journal radius Radial clearance Land length Groove width Groove depth Seal gap Seal length Oil viscosity
80.120 mm 0.254 mm 11.25 mm 2.0 mm 2.0 mm 0.8–0.5 mm 3.25 mm 0.006 Pa s
0.2 0
Fig. 7. Harmonic balance computation speed comparison: (a) finite difference solution—103,111 s; and (b) Chebyshev polynomial solution—438 s.
ARTICLE IN PRESS K.H. Groves, P. Bonello / Tribology International 43 (2010) 1639–1649
function f ðe, cÞ, then: ffit ðe, cÞ ¼
m X n X
~Þ ~ s ðc Cr,s Tr ðeÞT
ð14Þ
r ¼0s¼0
~ are e and where Cr,s are the coefficients to be determined, e~ and c c normalised over the range 71: e~ ¼ 2ðe1=2ðemax þ emin ÞÞ=ðemax emin Þ
ð15Þ
~ . It has been shown [16] that by invoking the ysimilarly for c orthogonality of the polynomials: ( ) m X n .n o X ~ ~ f ðer , cs ÞTr ðe r ÞTs ðc s Þ JTr J2 JTs J2 ð16Þ Cr,s ¼ r ¼0s¼0
ywhere JðÞJ represents the Euclidean norm: 2r þ1 p ~ ¼ cos 2s þ 1 p , c e~ r ¼ cos s mþ1 2 nþ1 2
ð17a; bÞ
Hence, calculation of the coefficients Cr,s requires ‘training’ data f ðer , cs Þ at input sets ðer , cs Þ. In the bearing models developed in Section 2.1, pðz, yÞ is a _ and (if feed-ports are included) c. The variables _ c function of e, e, _ _ c have no clear bounds. A technique is described that reduces e, the input space from four variables to three and allocates these variables finite bounds. Following the work of Chen et al. [10], Eq. (1) is rewritten in _ j: _ is greater than jec two alternative forms. If jej 2 1 @ @p 1 @ p 1 þ h3 h3 ¼ 12mðq1 sin y þ cos yÞ ð18Þ @y @z2 R2 @y
1643
_ =e. _ j is greater than jej: _ If jec _ where p 1 ¼ p=e_ and q1 ¼ ec 2 1 @ @p 2 @ p2 þ h3 h3 ¼ 12mðsin y þ q2 cos yÞ @y @z2 R2 @y
_ and q ¼ e=e _ . Given instantaneous values e, q _ c where p 2 ¼ p=ec 2 1/2 and c, either of Eqs. 18 and 19 can be solved for instantaneous distribution of p 1=2 (the influence of c is felt through the application of a feed-port boundary condition). Hence, p 1=2 are functions of only three variables e, c and q1 or q2. Moreover, q1/2 are bounded between 71. In this research, unlike [9,10], it is the (pre-truncated) gridpoint pressures pi,j that are identified, rather than the forces QR,T. This allows truncation of the grid of interpolated pressures at any arbitrary rupture pressure (as per Eq. (3)) and its subsequent quadrature (as per Eq. (2)) to obtain the forces QR,T. Moreover, the exclusion of cavitation from the polynomial approximation process reduces the complexity of the pressure relationship, thereby reducing the number of coefficients required to model it. The identified grid point pressures are given by: 8 _j _ Z jec < e_ p 1i,j ðe,q1 , cÞ, jej _ _ c , cÞ ¼ ð20Þ pi,j ðe, e, _j _ p ðe,q , cÞ, jej : ec _ o jec 2 2i,j In Eq. (20) it is the p 1=2 functions that are actually fitted as Chebyshev polynomials. It is clear that for each grid point two Chebyshev polynomial fits will be necessary. The coefficients of these fits are computed using, training, data from FD solutions of Eqs. (18) and (19). Although the coefficients of both fits have to be held in memory, only one of these is used for any single force
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.5
1
1.5
2
2.5
3
0
0.5
1
1.5
2
2.5
3
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
ð19Þ
Fig. 8. Short bearing two land model (Eq. (22)): (a) ps ¼239.2 kPa and U^ ¼ 0:128; and (b) ps ¼115.1 kPa and d^ ¼ 0:45.
ARTICLE IN PRESS 1644
K.H. Groves, P. Bonello / Tribology International 43 (2010) 1639–1649
computation. Hence, the force computation time is not significantly increased. When solving Eq. (18) or (19) for p 1=2 in the reduced ðe,q1=2 , cÞ space, it is not possible to accommodate an arbitrary boundary pressure pB since p B would be unknown as a consequence of _ from the input space. Hence, the having removed e_ and c reduction technique is only valid for zero (gauge) supply and outlet pressure. To allow use of arbitrary boundary conditions, the solution pi,j to the Reynolds equation is divided into its ‘dynamic’ and ‘static’ parts: _ Þ þpstat ðe, cÞ _ c, c pi,j ¼ pdyni,j ðe, e, i,j
ð21Þ
The first term of Eq. (21) is given by Eq. (20). It is the solution obtained under the actual vibratory conditions but with supply and outlet pressures set to zero. The second term of Eq. (21) is the solution obtained under the given boundary pressure conditions _ set to zero). This _ c for a static journal orientation of e, c (i.e. e, introduces an extra Chebyshev fit. However, since it is a function of only two variables, it adds only a small time penalty. The polynomial coefficients of the identified FD pressure function of a particular SFD are then deployed to compute the pressure instead of FD in any rotor-dynamic assembly where that bearing is installed. 2.3. Short sealed bearing approximations By dropping the pressure gradient in the y direction in Eq. (1) and applying the end-plate seal boundary condition of Eq. (9), it is
possible to derive analytical solutions to the Reynolds equation for a sealed bearing. Two such solutions are considered. The first model assumes a sealed bearing with two lands of length L each separated by a deep groove of fixed pressure ps. The pressure is then given as: p^ ¼
3 6E L 2 2 1 1 1 p^ s h^ 6EðL=RÞ2 ^þ ^ ^ ^s z z þ þ p z 3 R 3 3 4 2 2 h^ h^ þ d^ ðL=ls Þ
0 where E ¼ e0 cos y þ ec sin y, h^ ¼ 1 þ e cos y, e ¼ e=c, z^ ¼ z=L, d^ ¼ d=c 2 and p^ ¼ p=fmoðR=cÞ g. The second model assumes that the two lands of the bearing effectively form one large land of length 2L and that the bearing is sealed at either end of this land. The justification for such a model is that under tight sealing the feed-groove is somewhat bridged by the higher pressures caused by the sealing. The pressure is then given as: 3
2 ^ dÞ ^ 3 ðls =LÞ p^ ¼ ð6E=h^ ÞðL=RÞ2 ½ðz^ 1=4Þðh=
3.1. Description of the test rig The test rig used to produce the experimental readings presented in this paper was designed to be representative of a small business jet aeroengine [17]. Experimental results presented in this section are taken from reference [17], where
1 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0.5
1
1.5
2
2.5
3
0
0.5
1
1.5
2
2.5
3
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
ð23Þ
3. Testing
0.9
0
ð22Þ
Fig. 9. Experimental results (Ref. [17]): (a) ps ¼ 239.2 kPa and U^ ¼ 0:128; and (b) ps ¼ 115.1 kPa and d^ ¼ 0:45.
ARTICLE IN PRESS K.H. Groves, P. Bonello / Tribology International 43 (2010) 1639–1649
construction details of the rig may also be found. A schematic of the test rig is shown in Fig. 5. It consists of a heavy rotor, driven by an electric motor via a flexible coupling at its left hand end. At the driven end, the rotor runs in a rigidly housed self-aligning bearing. At the right hand end the rotor runs in a rollingelement bearing with a spring-supported housing. An SFD is located approximately at mid-section. Table 1 presents the SFD parameters used in the testing. The lowest undamped critical speed of the rig is 32.4 rev/s [17]. In this mode the rotor performs a conical motion with practically no flexure. Since the operating speed range is well below the first flexure critical speed (218 rev/s [17]) the equations of motion can be written as: kx ðl2 =l1 Þ2 xSFD þ Qx þðl3 =l1 ÞFx ðtÞ ¼ ðI=l21 Þx€ SFD ky ðl2 =l1 Þ2 ySFD þQy þðl3 =l1 ÞFy ðtÞ ¼ ðI=l21 Þy€ SFD
ð24Þ
where I ¼ 23:144 kg m2 . In this rig, the polar moment of inertia of the rotor is sufficiently small for gyroscopic effects to be negligible over the operational range of the rotational speed. Such an assumption was also made in the predictions of [17]. The above equations were solved using time/frequency domain methods to obtain the response for a given unbalance U, quantified here by the non-dimensional parameter U^ ¼ Ul1 l3 =ðIcÞ. All calculations in the following section were performed in Matlabs on a standard desktop pc with a 2.1 GHz dual core processor.
1645
3.2. Chebyshev vs FD To test the effectiveness of the proposed identification technique, results achieved using FD bearing models were tested against equivalent Chebyshev models. To efficiently identify the SFD over the input space, differing degrees of Chebyshev polynomial T (Eq. (14)) were required for each of the input variables. The highest degree of the Chebyshev polynomial used for each variable in the function is displayed in Table 2. It was discovered that the polynomials Tðq1=2 Þ were never above a degree of one; if a higher degree was forced into the fit, the identification process invariably returned zero values for the corresponding coefficients. This unexpected benefit of using the reduction technique with Chebyshev polynomials significantly enhances computation speeds both when executing the polynomials and achieving the coefficients. The one-off computation time for obtaining the coefficients corresponding to degrees in Table 2 was 85.7 s for identifying an FD pressure grid of ðNi Þa ¼ 31, ðNi Þb ¼ 7 and Nj ¼ 97. Fig. 6 shows an orbit from a time domain solution obtained using Matlab’s ode23s&. An FD bearing model was used for the SFD that incorporated end-plate seals, feed-groove and three equi-spaced pipe-flow feed-ports (Eq. (13)). In this example the spring support has been removed from the test rig model, thereby forcing the identified bearing model to operate in a highly nonlinear regime. Eqs. (24) were modified by removing the spring force term and including gravity in the second equation. Two polynomial fits were used to produce Fig. 6, as including
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.5
1
1.5
2
2.5
3
0
0.5
1
1.5
2
2.5
3
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
Fig. 10. FD–Chebyshev, two land deep groove end-plate sealed model: (a) ps ¼ 239.2 kPa and U^ ¼ 0:128; and (b) ps ¼ 115.1 kPa and d^ ¼ 0:45.
ARTICLE IN PRESS 1646
K.H. Groves, P. Bonello / Tribology International 43 (2010) 1639–1649
extremely high eccentricity degrades the fit at lower eccentricity and requires more polynomial coefficients. The fits were performed over 0 r e o 0:95 and 0:95 r e o 0:99. Splitting the fit into two meant the system needed to hold twice as many coefficients in memory; however, this caused a negligible time penalty. The plot is of a transient 0.1 s run from default initial conditions at 60 rev/s with an applied non-dimensional unbalance U^ ¼ 0:256, supply pressure of 239:2 kPa and seal gap of 0.45c. The respective orbit, computed using Chebyshev polynomial fits, is not displayed since the two plots are indistinguishable. The key difference between the two orbits is the execution time. The FD model took 3706 s to run, while the Chebyshev equivalent required only 137 s i.e. 3.7% of the time. Fig. 7 shows plots of amplitude against speed for the springsupported rig with a centralized SFD. The response was solved by harmonic balance (HB) using the arc-length continuation technique developed in [13]. The bearing tested had end-plate seals, feed-groove and three equi-spaced fixed pressure feedports (Eq. (10)). The supply pressure was set at 239.2 kPa with a seal gap of 0.45c and unbalance of U^ ¼ 0:128. It is clear that the use of Chebyshev polynomials is extremely well suited to the approximate solution method of HB. The smoothing effect of the identification technique clearly allows rapid, stable operation of the HB code. The FD solution points can be seen to cluster, reducing arc length to maintain accuracy and therefore the process progresses extremely slowly. The time saving is excellent; the FD solution took 103,111 s to run, while the Chebyshev model only required 438 s, i.e. 0.42% of the computation time.
3.3. FD models vs analytical models The previous section verified that Chebyshev identification provides a reliable means of implementing FD SFD models. Chebyshev will be henceforth be used to implement the bearing models of Section 2.1. The unbalance response of the test rig obtained with these ‘FD-Chebyshev’ models are compared in this section to the results obtained using the analytical bearing models of Section 2.3. Fig. 8 shows the response obtained with the analytical two land end-plate sealed model (Eq. (22)). These graphs exhibit a ‘spring-softening’ effect (i.e. a lowering of the critical speed). This effect was observed as being due to a combination of tight sealing and high supply pressure. The cause of this effect is the second term of Eq. (22), which contains the effect of the supply pressure ^ This and sealing and varies with y due to the presence of E and h. variation means that it does not tend to cancel out when integrated to give the forces, unlike the third term of Eq. (22). When ps is low or d^ is large, the second term of Eq. (22) becomes small and the softening virtually disappears, as evident from Fig. 8(a). The softening effect was not observed in the measurements of Fig. 9. Nor was it observed in predictions obtained for this test rig in [17] using the ‘l theory’, that assumes a pressure distribution of the form [5,6]: pðz, yÞ ¼ pshort,unsealed ðz, yÞ þ lplong ðyÞð0:5z^ Þ
ywhere l is the end-leakage factor. The first term of Eq. (25) is ^ Eq. (22) but with the seal gap ratio d-1 (i.e. unsealed SFD). The
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.5
1
1.5
2
2.5
3
0.5
1
1.5
2
2.5
3
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
ð25Þ
Fig. 11. Short bearing one land model (Eq. (23)): (a) ps ¼239.2 kPa and U^ ¼ 0:128; and (b) ps ¼115.1 kPa and d^ ¼ 0:45.
ARTICLE IN PRESS K.H. Groves, P. Bonello / Tribology International 43 (2010) 1639–1649
second term of Eq. (25) contains the ‘long bearing’ solution to the Reynolds equation, obtained by dropping its axial pressure gradient term (second term on LHS of Eq. (1)). It is noted that the model of Eq. (22) was obtained by dropping the circumferential pressure gradient term from the Reynolds equation. However, it was ascertained that this had nothing to do with the softening behaviour of Fig. 8. Indeed, Eq. (22) was verified to be computationally sound through FD. Fig. 10 shows the amplitude-speed plot for an end-plate sealed SFD model with a central deep groove of fixed pressure ps computed by FD-Chebyshev. This plot is almost indistinguishable from Fig. 8, despite the fact that the results of Fig. 10 are based on the full (two-dimensional) Reynolds equation. This means that the omission of the circumferential pressure gradient to solve the Reynolds equation was verified to have negligible effect for an end-plate sealed damper whose L/R is similar to that in [5,6]. This fact appears to contradict the l-factor theory of Eq. (25). However, comparing Eq. (25) with the (FD-verified) Eq. (22), it is observed that both expressions involve the addition of a circumferentially varying pressure distribution to the short, unsealed SFD expression. In the case of Eq. (22), this additional distribution is its second term. Since this term has no direct equivalence to the second term of Eq. (25), it is proved that it is not possible to find a theoretical link between l and the SFD parameters. This explains the unsuccessful attempt in [5]. Fig. 11 shows the results obtained using the analytical one land end-plate sealed model (Eq. (23)). The softening effect is absent since there is no lubricant supply included in the model. Comparing with the experimental results of Fig. 9, the results of
1647
Fig. 11 are clearly over-damped. This shows that the two lands do not quite act as one under tight sealing. On the other hand, the previous results of Figs. 8 and 10, based on two independent lands and neglecting any damping contributions from the groove, are under-damped. Fig. 12 presents the results obtained using an FD sealed bearing model that incorporates both axial and circumferential pressure gradients, the effect of the feed-groove but without any feed-ports. The results show marked improvement upon the previous models. Damping in the system is increased when compared to Fig. 8 or 10, while being significantly less than in the one-land model of Fig. 11, therefore proving the importance of feed-groove effects. The results of Fig. 12 show no softening since no feed-ports were included. The effect of adding three equi-spaced feed-ports of fixed pressure of ps (Eq. (10)) is shown in Fig. 13. It is clear that the artificial softening effect at tight sealing has reappeared and at loose sealing the critical speed shifts above the undamped critical. These results confirm that the assumption that a feed-port pressure is always fixed, irrespective of the journal dynamics, is unrealistic and results in artificial shifting of the critical speeds. The results presented in Fig. 14 use the pipe flow boundary condition (Eq. (13)), with a feed-port of length of lp ¼15 mm and area of Ah ¼0.7 mm2. It is apparent that modifying the lubricant supply boundary condition has almost removed the critical speed shifting effect on the predictions. It is observed that the results of Fig. 14 show little improvement on the results of Fig. 12. Moreover, comparing the curve for
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.5
1
1.5
2
2.5
3
0
0.5
1
1.5
2
2.5
3
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
Fig. 12. FD-Chebyshev model with seal and groove, no feed-ports: (a) ps ¼ 239.2 kPa and U^ ¼ 0:128; and (b) ps ¼115.1 kPa and d^ ¼ 0:45.
ARTICLE IN PRESS 1648
K.H. Groves, P. Bonello / Tribology International 43 (2010) 1639–1649
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.5
1
1.5
2
2.5
3
0
0.5
1
1.5
2
2.5
3
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
Fig. 13. FD-Chebyshev model with seal, feed-groove and fixed pressure feed-ports: (a) ps ¼239.2 kPa and U^ ¼ 0:128; and (b) ps ¼115.1 kPa and d^ ¼ 0:45.
d^ ¼ 0:45, U^ ¼ 0:128 in Fig. 14(a) with the corresponding one in Fig. 14(b), it is clear that increasing the supply pressure did not improve the predicted damping capacity of the SFD. The reason for this is that, for all the simulations, the SFD pressure was invariably above the chosen cut-off pressure (absolute zero), unlike the experimental observations (Fig. 9). The presence of supply pressure only influences the damping effect if there is cavitation present in the system. If a more demanding test case were used (i.e. one that forced the bearing models to cavitate), the positive effect of supply pressure upon the damping characteristics would be observed in the predictions. A key benefit of carrying out pressure truncation post-identification is that the simulations could be redone at a variable cut-off pressure, as in [5,6] where experimental cut-off pressures were used.
4. Conclusions This paper considered the modelling of end-plate sealed SFD bearings with groove and feed-ports. Finite difference models based on the incompressible Reynolds equation were first developed. A novel identification technique was then presented that provided a means for the FD models to be rapidly deployed. No assumption was made about the range of the input variables, the cavitation pressure or the axial symmetry of the pressure boundary conditions. Rather, arbitrary boundary pressure conditions may be applied at any location in the film while still restricting the number of input variables to a manageable number
(three). The Chebyshev technique has been proven to be particularly suited to use in harmonic balance solution schemes, where staggering time savings were observed. This work opens the possibility of using bearing models of increased accuracy within multi-SFD full aeroengine models. A highly instructive study was made on a simple rig in which the unbalance response was calculated using identified FD SFD models and simple analytical SFD models. It has been demonstrated that the effect of omitting the circumferential pressure gradient to solve the Reynolds equation is negligible for a short SFD, even if it has substantial sealing. This proves that it is not possible to theoretically relate the end-leakage factor (l) to the SFD parameters. The common usage of a fixed pressure condition at a groove or feed-port was shown to be inadequate since it resulted in artificial shifting of the critical speed as the sealing was tightened. The use of an improved feed-port boundary condition virtually eliminated the problem and gave superior agreement with the experimental observations. Despite the improved accuracy of the SFD model, the level of damping was somewhat underestimated relative to the experimental results. This could be partly attributed to the inherent degree of uncertainty in the experimental parameters, e.g. the seal gap dimension, which is only a fraction of the radial clearance for tight sealing. Another possible reason for the discrepancy is the simplicity of the groove model, which neglected inertia effects. Current efforts are being directed towards seeking a correction for fluid inertia while preserving the advantages of the presented identification scheme.
ARTICLE IN PRESS K.H. Groves, P. Bonello / Tribology International 43 (2010) 1639–1649
1649
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.5
1
1.5
2
2.5
3
0
0.5
1
1.5
2
2.5
3
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
Fig. 14. FD-Chebyshev model with seal, feed-groove and pipe flow feed-ports: (a) ps ¼239.2 kPa and U^ ¼ 0:128; and (b) ps ¼ 115.1 kPa and d^ ¼ 0:45. tx
References [1] Bonello P, Hai PM. Computational studies of the unbalance response of a whole aero-engine model with squeeze-film bearings. Trans ASME J Eng Gas Turbines Power 2010;132(3):032504 [7pp.]. [2] Pietra LD, Adiletta G. The squeeze film damper over four decades of investigations. Part I: characteristics and operating features. Shock Vibr Dig 2002;34:3–26. [3] San Andres LA, Vance JM. Effect of fluid inertia on squeeze-film damper forces for small-amplitude circular-centered motions. ASLE Trans 1987;30:63–8. [4] Zhang JX, Roberts JB. Solutions for the combined motion of finite length squeeze film dampers around the bearing center. ASME J Tribol 1996;118:617. [5] Dede MM, Dogan M, Holmes R. The damping capacity of a sealed squeeze film bearing. ASME J Tribol 1985;107:411–8. [6] Holmes R, Dogan M. The performance of a sealed squeeze-film bearing in a flexible support structure. Proc IMechE Part C 1985;199:1–9. [7] Jiangang Y, Rui G, Yongwei T. Hybrid radial basis function/finite element modelling of journal bearing. Tribol Int 2008;41:1169–75. [8] MSC. Nastran release guide. MSC. Software corporation; 2005. p. 150–6.
[9] Rodrigues FA, Thouverez F, Gibert C, Jezequel L. Chebyshev polynomials fits for efficient analysis of finite length squeeze film damped rotors. ASME J Eng Gas Turbines Power 2003;125:175–83. [10] Chen Z, Jiao Y, Xia S, Huang W, Zhang Z. An efficient calculation method of nonlinear fluid film forces in journal bearing. STLE Tribol Trans 2002;45:324–9. [11] Castelli V, Shapiro W. Improved method for numerical solutions of the general incompressible fluid film lubrication problem. ASME J Lubr Technol 1967;89:211–8. [12] Cameron A. Basic lubrication theory. London: Ellis Horwood; 1971. [13] Bonello P, Brennan MJ, Holmes R. Non-linear modelling of rotor dynamic systems with squeeze film dampers—an efficient integrated approach. J Sound Vibr 2002;249:743–73. [14] Feng NS, Hahn EJ. Effects of gas entrainment on squeeze film damper performance. ASME J Tribol 1987;109:149–54. [15] Munson BR, Young DF, Okiishi TH. Fundamentals of fluid mechanics. New York: John Wiley and Sons; 1998. [16] Dahlquist G, Bjorck A. Numerical methods. New York: Prentice-Hall; 1974. [17] Levesley MC, Holmes R. Experimental investigation into the vibration response of an aero-engine rotor-damper assembly. Proc IMechE Part G 1994;208:52–66.