Dynamics of a high speed coned thrust bearing with a Navier slip boundary condition

A modified Reynolds equation for flow dynamically represented as incompressible is used to model the dynamics of a thin film bearing with slip flow and a rapidly rotating coned rotor. Previous studies including a Navier slip length shear condition on the bearing faces are extended to investigate applications with a coned bearing gap. A modified Reynolds equation for the film flow is coupled, through the pressure exerted by the fluid film, to the dynamic motion of the stator. Introducing a new variable leads to explicit analytical expressions for the pressure field and force on the stator with the equation for the time-dependent face clearance transformed to a nonlinear second-order non-autonomous ordinary differential equation. The face clearance for periodic axial motion of the coned rotor is obtained using a stroboscopic map solver; a focus is investigating bearing behaviour under extreme conditions. The coupled fluid flow and unsteady bearing dynamics are examined for a range of configurations to evaluate potential face contact over a range of bearing surface conditions.

Typical bearing geometries include a journal bearing with a supporting cylindrical sleeve containing a rotating cylindrical shaft separated by a thin air film, as studied by Belforte et al. [1]; a slider bearing comprising two non-parallel moving plates separated by a thin lubricating air film, as considered by Witelski [2]; and a thrust bearing which has an important significance dynamically in turbomachinery applications operating with very high rotational speeds.
In advanced engineering applications, an increased lubrication force on the stator can be obtained by incorporating a coned rotor. In early stability studies of coned bearings, Etison [3] studied the stability and response of a coned-face seal in the presence of rotor axial runout and assembly misalignment. An optimal coning angle was identified, as the one at which the stability is maximised. While examining a similar coned seal under a range of pressurisations, Green and Etison [4] concluded that axial stability exists if the film converges in the direction of the pressure drop or if there is a sufficiently high support stiffness in a divergent film. A divergent air film can lead to seal failure in cases of low support stiffness, low pressure gradient or large coning angles. Green [5] extended the bearing study, examining positive and negative coning angles consistent with a convergent air gap. To avoid face contact and maintain positive stiffness, a critical angle exists depending on geometric restrictions. Green and Barnsby [6] showed negative coning is resilient against hydrostatic instabilities, although still prone to dynamics instabilities.
Additional centrifugal inertia effects, which can be important in turbomachinery applications requiring operation at higher rotational speeds, were studied by Garratt et al. [7] for a thrust bearing. The bearing dynamics were investigated when one face undergoes prescribed periodic axial oscillations of amplitude smaller than the equilibrium film thickness, and the other face moves axially in response to the film dynamics. Considering similar geometry with incompressible flow and incorporating a coned rotor, Bailey et al. [8] provided more extensive analytical investigations and examined the effect of prescribed axial oscillations with amplitude larger than the equilibrium film thickness. Results indicated that the fluid film thickness can become very small causing the classical no-slip velocity condition to become invalid. Therefore, a slip condition on the bearing faces is included in the analysis of a parallel bearing with a small face clearance by Bailey et al. [9], with results giving no face contact, although the bearing gap can become very thin. In this work, the case of a coned bearing with a small face clearance is analysed to evaluate if contact can occur.
Fluid flows confined to micro-and nano scales are subject to surface effects beginning to dominate over volumerelated phenomena, requiring accurate details of the fluid-solid interaction. In the case of a gas film, the Knudsen number Kn = l/ĥ 0 , with l as the mean free path andĥ 0 the characteristic fluid thickness, is usually used to classify mathematical models for thin gas-flow regimes. The fluid is regarded as continuum for small Knudsen number (Kn ≤ 10 −3 ) where the no-slip boundary condition is valid. For larger Knudsen number between 10 −3 and 10 −1 , a continuum model is still valid, but slip boundary conditions are usually implemented, which is the flow regime of interest in the present work. The flow is in a transition region for a Knudsen number between 10 −1 and 10, where a modified continuum model is required and in the case of a larger Knudsen number (Kn ≥ 10), the free molecular flow can be modelled by molecular dynamics [10]. When classifying very thin flow regimes for a liquid, the Knudsen number can be used, but with the mean free path replaced by the intermolecular distance between molecules, giving exceedingly small fluid thickness using the above classification. When considering liquid flow over a hydrophilic surface, even on a nano scale, the classical no-slip boundary condition appears to hold; however, when considering a hydrophobic surface an apparent slip velocity has been observed just above the solid surface.
Slip flow has been studied in an inclined plain slider bearing by Burgdorfer [11] who devised a first-order slip model where the boundary slip velocity was identified at a mean free path distance from the wall. A non-axisymmetric thrust bearing with foil pads on the rotor face was examined by Park et al. [12] where the gas flow is coupled to the bearing structure. A Reynolds equation with slip was derived, and the bearing dynamics were investigated for rotor displacement of small amplitude in comparison to the film thickness, using existing perturbation analysis with results given for the no-slip and slip conditions. More extensive details on previous bearing studies incorporating slip flow are given in Bailey et al. [9].
The increasing performance advantages for bearings to operate with reduced gap between the rotor and stator motivated the study of a thin film model valid for slip flow. This study extends the formulation, analysis and predictions of the previous study of Bailey et al. [8] in which a coned thrust bearing is considered within a classical no-slip condition which predicts that face contact will not occur, although a period of very close proximity may exist. New bearing technologies are requiring bearing gap thicknesses to become very small and maintain appropriate dynamic properties to sufficiently resist the possibility of face contact and maintain a minimum gap thickness. A focus of this work is to incorporate a more general slip condition more relevant to modelling the fluid film under extreme operating conditions and evaluating the corresponding film flow and bearing behaviour. A slip length model is developed which couples the thrust bearing structure to flow dynamics induced by the fluid force exerted. An extended bearing model incorporates high-speed rotational bearing operation through centrifugal fluid inertia. The formulation of the coupled governing equations is presented in Sect. 2, including the Reynolds equation in terms of a slip boundary condition characterised by a slip length parameter, and the stator displacement equation. A representative single nonlinear second-order non-autonomous ordinary differential equation for the face clearance is derived in Sect. 3. To solve for the periodic bearing face clearance iteratively, a stroboscopic map solver is formulated and numerical procedure identified with an extension to the solver, identifying the value of slip length corresponding to a prescribed minimum face clearance (MFC) in the period. Detailed evaluation of the steady coned bearing is given in Sect. 4 and for the dynamic bearing in Sect. 5 to explore the dynamical behaviour and influence of increasing effects of slip boundary conditions through selected parameter studies. The possibility of contact through a parameter study, including the slip length and coning angle is examined in Sect. 6. To verify the numerical solution in the limit of small face clearance asymptotic analysis is carried out, and two different symbolic computations were used to verify the correctness of the lengthy algebra.

Geometric configuration
A simplified mathematical model of a fluid lubricated bearing containing a coned rotor, such as developed in Bailey et al. [8], is extended by incorporating a slip condition on the bearing faces. The coaxial annular rotor and stator have inner and outer dimensional radii ofr I andr O , respectively, and the rotor has angular rotationˆ . The rotor has a fixed coning angle, assumed to be small giving cosβ 1 and sinβ β , with cases of a positively and negatively coned bearings considered, referred to as PCB and NCB, respectively. Pressuresp I andp O are imposed at the inner and outer radii of the bearing, respectively, allowing a pressure gradient to drive the radial flow. The axisymmetric rotor-stator clearance is defined with reference to the point of minimum film thickness; at the inside, radius for a PCB given bŷ and at the outer, radius for an NCB is given bŷ Hereĥ s is the stator height, and the rotor height is given byĥ r =ĥ r (r j ,t) − (r −r j ) tanβ with j = I, O for the PCB and NCB, respectively, with a temporal dependenceĥ r (r j ,t) at the reference height and a spatial dependence (r −r j ) tanβ. The prescribed axial motion of the rotor is modelled byĥ r (r j ,t) = ĥ 0 sin ωt , with amplitude ĥ 0 whereĥ 0 is the equilibrium film thickness at the inner and outer radii for a PCB, NCB, respectively.
Deriving the velocity boundary conditions for a slip bearing requires the normal and tangential velocities on the rotor and stator to be specified, due to the additional slip component on the bearing face. Both PCB and NCB are considered separately in the coordinate system (ẑ,r ,θ) with velocities (ŵ,û,v). Configurations of the PCB and NCB are shown in Fig. 1.
The PCB and NCB have the same normal and azimuthal tangent given bŷ n r = (cosβ, sinβ, 0),t 2 = (0, 0, 1), respectively, whereas the axial tangents are given bŷ Fig. 1 Axially symmetric geometry and notation of a bearing with a PCBβ > 0 and NCBβ < 0 with normal vectorn r and two orthogonal tangential components t 1 and t 2 on the face of the rotor with the origin of the coordinate system at the inner radius of the annulus for the PCB and NCB, respectively, requiring separate radial velocity boundary conditions to be derived. The stator lies parallel to the radial direction, giving the normal vectorn s in the negativeẑ direction and two orthogonal tangential components, t 3 and t 4 on the face of the stator in the radial and azimuthal directions, respectively.
A first-order Navier formulation slip model is implemented as used by Burgdorfer [11] for the case of an inclined, plain slider bearing, where the slip velocity is proportional to the tangential component of the face shear. The velocity boundary conditions comprise a jump in the tangential velocity across the fluid-solid interface, which is equal to the slip velocity induced due to wall shear, and continuity of the normal velocity with a no flux condition. Thus, the velocity boundary conditions on the rotor and stator are given bŷ u ·t 1 −û r ·t 1 = 2l sêi jn jt1,i ,û ·t 2 −û r ·t 2 = 2l sêi jn jt2,i ,û ·n r −û r ·n r = 0 atẑ =ĥ r , u ·t 3 −û s ·t 3 = 2l sêi jn jt3,i ,û ·t 4 −û s ·t 4 = 2l sêi jn jt4,i ,û ·n s −û s ·n s = 0 atẑ =ĥ s .
The rotor velocity components denoted by superscript r, are given byû r = (∂ĥ r /∂t, 0,ˆ r ) corresponding to prescribed axial motion and fixed azimuthal rotation. The stator is allowed to move axially due to the interaction with the fluid, giving the stator velocity, denoted by superscript s, asû s = (∂ĥ s /∂t, 0, 0). Velocity boundary conditions on the rotor in (5) are given bŷ using the normal and tangent vectors given in (3) and (4) and the rate of strain tensor components as given in Batchelor [13, p. 602]. In an axisymmetric configuration, the azimuthal derivatives do not appear. The velocity boundary conditions on the stator (5), givê A model for the incompressible fluid flow through the bearing is derived from the Navier-Stokes momentum and continuity equations in axisymmetrical coordinates. Variables are non-dimensionalized with respect to the typical bearing radiusr 0 , rotor velocityˆ r and time scaleT = 1/ω, with dimensionless time variable t =ωt whereω is the angular frequency. Dimensionless velocities are taken to beû/Û ,v/ˆ r 0 andŵ/ĥ 0T −1 with dimensionless radius and height given by r =r /r 0 and z =ẑ/ĥ 0 , respectively. The dimensionless slip length is given by l s =l s /ĥ 0 . The associated radial and azimuthal Reynolds numbers and relative ratio are given, respectively, by The aspect ratio δ 0 , squeeze numberσ and Froude number Fr are defined as respectively, whereĝ is the acceleration due to gravity,ρ the density andμ the dynamic viscosity. For thin film bearings δ 0 1 and to ensure that the effects of viscosity are retained at leading order, the pressure is scaled asP =μr 0Û /ĥ 2 0 . Classical lubrication theory neglects inertia due to the reduced Reynolds number Re U δ 2 0 1; however, as was done in Garratt et al. [7], the centrifugal inertia is retained to include cases of high rotational-speed bearing operations for which terms of the order Re U δ 2 0 (Re * ) 2 are considered to be of O(1), with (Re * ) 2 1. The squeeze numberσ characterises any time-dependent effects, whilst the Froude number Fr parameterises the importance of the gravitational effects relative to the radial flow although gravity can be neglected with Re U δ 2 0 Fr −2 1; this is consistent with lubrication theory provided that the Froude number is O(1). Importantly for more general application, it can be shown that under the present approximation, the flow field with a very thin film air bearing (nano scale gap) the usual incompressible solenoidal condition remains valid for compressible flow when the following conditions are satisfied: as found in Bailey et al. [9]. To leading order, where terms of O(δ 0 ) are neglected, the governing lubrication equations become where the speed parameter η = Re U δ 2 0 (Re * ) 2 =ρr 0ĥ 2 0ˆ 2 /μÛ characterises the importance of centrifugal inertia. For the derivation, see Bailey et al. [8]. Using (11) and (12) and (63)-(65) from Appendix 1, analytic expressions can be found for u, v and w in terms of the unknown stator position h s (t).
In the case of a small coning angle β with prescribed periodic oscillations, the dimensionless film thickness is given by where h s (t) is the height of the axially displaced stator and a =r I /r 0 is a measure of the bearing width, 0 < a < 1. The coning angle is assumed small, giving sinβ β = O(δ 0 ) and cosβ = 1 + O(δ 0 ) leading to the scalinĝ β = βδ 0 with β = O(1), giving consistency with the lubrication condition, such that ∂h r /∂r −β.
Under these conditions, the leading-order non-dimensional rotor and stator velocity boundary conditions given in (6) and (7), respectively, become The only dependence on the coning angle is given in the axial rotor velocity boundary condition. The governing equation for flow in the bearing is obtained by integrating the continuity equation (12) between the rotor and the stator, and applying the Leibniz integral rule gives To leading order, imposing the velocity boundary conditions from (14) and using the expression for the radial velocity (63) gives the modified Reynolds equation where λ = 3η/10 and σ = 12σ . The Reynolds equation (16) has no explicit dependence on the coning angle and expresses the relationship between the pressure p and film thickness h. Pressure boundary conditions at the inner and outer radii of the bearing are defined in dimensionless variables as Axial displacement of the stator is modelled using a spring-mass-damper model giving the axial position of the stator in dimensional variables as with h s (t) defined at the origin of reference as shown in Fig. 1 with the film thickness at the inner radius taken as h = 1 to give the stator height h s = 1 for PCB and h s = 1 + (1 − a)β for an NCB. In (18), F(t) is the resultant dimensionless axial force on the stator defined by where p a is the ambient pressure above the stator and the force coupling dimensionless parameter given by α = μÛ /(mω 2 δ 3 0 ), considered to be of O (1). Bearing quantities D a =D a /(mω) and K z =K z /(mω 2 ) are coefficients or dimensionless linear damping and effective restoring force, respectively, withD a andK z as their corresponding dimensional values andm the mass of the stator.
The corresponding radial flux through the bearing is obtained from the integration of the radial flow velocity over an azimuthal bearing cross section and is given in Appendix 1, Eq. (66). The steady-state streamfunction along a radial cross section is obtained from the integration of the velocity components (u, w) and is given by (67).

Bearing gap equation
Displacement of the stator is coupled by the periodic forcing of the rotor through the fluid film flow. To solve the Reynolds equation (16) and stator equation (18) simultaneously, it is mathematically convenient to introduce a measure of the bearing gap through the time-dependent magnitude of the MFC as occurring at the inner radius for a PCB and outer radius for an NCB. It is shown the system-dependent quantities can be characterised by g(t) which is given as the solution of an explicit integro-differential equation.
Integrating the modified Reynolds equation (16) and imposing the pressure boundary conditions (17) gives where the integrals G(g, r ), H (g, r ) and L(g, r ), for a PCB, are calculated from with exact evaluations given in Appendix 2 as (69), (71)  The force on the stator is obtained from the integration of the pressure field (21) using the force integral in (19), giving Expressions for A and B are defined as The integrals G I , H I and L I are given by with exact evaluations for a PCB given in Appendix 2, Eqs.  18) gives a nonlinear second-order non-autonomous ordinary differential equation for the periodic time-dependent minimum bearing gap as where for a PCB. An NCB only differs in the term; In this way, a simple ordinary differential equation is defined for the minimum bearing gap, as a function of the unknown stator height, h s (t). Dynamically Eq. (29) corresponds to a harmonically forced oscillator with nonlinear damping coefficient D(g, β, l s ) and effective restoring force S(g, λ, β, l s ).
The total stiffness of the system is defined as comprising a structural component K z and a fluid stiffness term using the expression for A(g, λ, β, l s ) in (26). Details on solving Eq. (29) are given in the next subsection where a stroboscopic map solver is formulated, using Newton's method to find periodic solutions. To compute solutions for an increased value of the slip length l s + l s , an Euler-Newton scheme (parameter continuation) is developed. The Euler procedure can be directly extended to find threshold values of the slip length corresponding to a minimum value of the MFC over the period, with the limiting case of contact given by g min = 0.

Numerical method
Solutions to Eq. (29) for a fixed bearing configuration are denoted by the vector g(g(t), z(t)), for a given initial conditions g(t 0 ) = g 0 , z(t 0 ) = z 0 and fixed value of slip length l s , and are sought from the equivalent system of two first-order differential equations, It is anticipated that for a fixed value of slip length l s and period of prescribed rotor oscillation T , the system of equations (33) has periodic solutions. Thus, a mapping is constructed which advances any initial condition g 0 at time t 0 by a time T , defining a stroboscopic map φ(T ; g 0 , t 0 ). This R 2 → R 2 map, integrates the system of equations (33) forward through one period of the forcing term. To identify periodic solutions, the fixed points of the stroboscopic map g(t) = g(t + T ) are found iteratively, corresponding to the condition giving periodic solutions g(g(t), z(t)).
Solutions are found numerically from an iterative Newton's method, given an initial guess valueg 0 . Successively improved iterates g 0 are given from the numerical iterative scheme where the Jacobian matrix is given by To find the elements of the Jacobian matrix J (T ), the following auxiliary system of first-order differential equations is defined; with initial conditions The values of the elements of the Jacobian matrix are given by the values of the auxiliary variables at the end of the forcing period, t = T , corresponding to the time at which periodicity is tested for. Thus, for any given initial condition, a solution of the system of equations (33) and (37) for t 0 ≤ t ≤ T can be found using MATLAB's ode15s routine. This procedure is repeated recursively, with the improved value of g 0 used in the system of equations (33). The scheme is successively repeated until a prescribed tolerance, tol, is achieved, i.e. | g(T ) − g 0 (t 0 ) |≤ tol and a periodic solution is obtained.
To find a periodic solution for increasing slip lengths l s + l s , a new initial condition is needed. This is done via parameter continuation, which means the solution is now defined with the slip length as a parameter; g(T ) = φ(T ; g 0 , t 0 , l s0 ).
To define a new initial condition g 0 for l s + l s , first, the derivative of G(φ(T ; g 0 , t t ), g 0 ) = φ(T ; g 0 , t t )−g 0 = 0 is taken with respect to the slip parameter l s , i.e.
Then an Euler predictor step is performed and a new initial condition g 0 for l s + l s is defined. The inverse of the Jacobian matrix is as found previously, at the value of l s for which a periodic solution was obtained.
To find the values of ∂g(T )/∂l s the following auxiliary system of first-order differential equations is defined; with initial conditions ∂g ∂l s = 0, ∂z ∂l s = 0.
This system of equations can be coupled with the previous augmented system of equations and solved using the same MATLAB routine.
The above Newton procedure is iterated until convergence is achieved and a periodic solution for l s + l s is found. To ensure Newton's method converges it may be necessary to reduce the value of l s ; typically the initial value of l s was successively divided by two until a convergent solution was found.
A major advantage of the above numerical (Euler) procedure is that it can be directly extended to find threshold values of the slip length l s corresponding to a minimum value of the MFC over the period, g min , with the limiting case of contact given by g min = 0. In this approach, the slip length l s is taken as a new dependent variable in the Newton scheme to determine the value of the unknown vector g = (g(T ), z(T ), l s ) given an initial guess value g 0 = (g 0 , z 0 , l s0 ). Correspondingly additional constraint equation g min − g * = 0 is added, with g * as the prescribed g min and the corresponding Jacobian matrix given by The additional terms in the Jacobian matrix (43) with respect to (36), i.e. the last row and column in (43), are obtained from the solution of the augmented system of first-order differential equation (37) and (41). The values in the last row are determined at the time when g min is achieved. The threshold slip length value, at the specified g min , can be stepped down successively to the point of contact, g * = 0. Continuation is used for a new value of the specified g * to get with the Jacobian matrix as defined in (43) and g * < 0. The value of ∂g/∂g * is evaluated numerically by a firstorder forward finite difference approximation, in terms of the obtained solution g(g * ) and a new auxiliary solution g(g * ), corresponding to a specified g min ,g * = g * + g * , with 1. This leads to a new initial condition being defined, allowing a periodic solution for a decreased g min value g * + g * to be found.

Effects of slip on steady-state bearing flow
In this section, the influences of the slip effects on the fluid flow in a steady-state bearing are examined. The stator is fixed axially at h s = 1 and rotor axially fixed at h r = 0 at the inner or outer radius for a PCB or NCB, respectively, with g = 1 and = 0. The rotor has constant azimuthal velocity. Two classes of bearing pressurisation with typical values are used: internal pressurisation with p I = 2, p O = 1 and external pressurisation p I = 1, p O = 2. Results reported here correspond to a = 0.2, i.e. a wide bearing of width 1 − a = 0.8. Similar results occur for other values of a, with smaller variations on the pressure profile and larger variations on the velocity profiles as the value of a increases.
The pressure is investigated for the presence of non-monotonic behaviour by examining the derivative of the pressure (21) with respect to r , A bearing with negligible inertia effects λ = 0, under ambient pressure ( p I = p O ) has pressure gradient ∂ p/∂r = 0, while under internal and external pressurisation the pressure field minima are dictated by ∂G/∂r = 0, with Thus, both PCB and NCB have a monotonic pressure field for finite slip lengths, as a minimum occurs only at r = ∞. The location of the minima when including inertia effects λ = 0 is given by ∂ p/∂r = 0, requiring the numerical evaluation of 1 G(1, 1) and G(1, 1) and L(1, 1) as given in Appendix 2. Figure 2 shows the minimum position r min with the increasing speed parameter for the case of a PCB under external and internal pressurisations. A monotonic pressure field exists when inertia effects are negligible, i.e. λ = 0, independent of the slip length. Increasing the speed parameter, λ = 0, causes a minimum to develop at the inner and outer radii, for external and internal pressurisation respectively, before moving into the bearing. The value of r min decreases or increases with the increasing slip length in the case of external or internal pressurisations, respectively. Substituting the location of the minimum r min into the pressure equation (21) gives the minimum pressure decreasing monotonically with the increasing speed parameter, which can become negative. Increasing the slip length causes the pressure to increase, tending to the large finite slip length limit. The radial flux in a steady bearing can be examined by integrating the steady-state Reynolds equation from (16) to give where the left-hand side of (49) is proportional to the flux (66). The constant C is obtained by substituting the steady pressure field (21) into (49), giving the steady flux as A bearing with negligible inertia effects λ = 0 and ambient pressure has zero flux, corresponding to azimuthal flow only, as can be identified from velocity field expressions (63)-(65). Imposing internal or external pressurisation across the bearing gives the radial flux as being unidirectional outward or inward flow, respectively. For non-zero speeds λ = 0, a critical value λ c can arise where zero flux through the bearing is achieved; the flow driven by the pressure gradient exactly balances the flow due to high speed rotation. The critical speed parameter λ c , found from (50), is given by with L(1, 1, β, l s ) in (73). Numerical evaluations give L(1, 1, β, l s ) > 0 for the range of parameters used, giving zero flux existing for an externally pressurised bearing only, with low speed parameters λ < λ c having negative flux and high speed parameters λ > λ c having positive flux. The streamlines in Fig. 3 give the radial path of the fluid flow through the bearing for an externally pressurised PCB and NCB. Choosing the critical speed parameter for a no-slip bearing, which is consistent with the estimated value, λ c (l s = 0), gives zero flux through the bearing. Increasing the slip length, but keeping the same speed parameter, results in a negative flux in the bearing, increasing in magnitude as the slip length increases. In Fig. 3, a comparison for slip length l s = 1, shows that an NCB has unidirectional flow but a PCB retains some recirculation.
Due to the algebraic complexity of the equations, the fluid flow is studied by examining the velocity field in the asymptotic limit of a large finite slip length. To leading order the radial, azimuthal and axial velocity components are given by with G l s given by where subscripts + and − denote a PCB and NCB, respectively. Expressions in (53) and (54) have dependence on β and l s and are non-negative. The azimuthal velocity has rigid body motion taking the value of the average between the azimuthal velocity of the rotor and the stationary stator. The radial and axial velocity are proportional to the slip length, with the radial velocity independent of the axial coordinates and the axial velocity the only component to change sign with a PCB and NCB, since G l s > 0 always. These limits match the numerical analysis of the full expressions. Examining the flux in the limit of large slip length gives which gives the flux proportional to slip length through G l s , as given in (53) and (54) for PCB and NCB, respectively. The corresponding critical speed parameter is given by confirming only an externally pressurised bearing as having zero flux. Equation (56) is consistent with (51) in the limit of large l s , with the limit of L as L l s = 5(1 − a 2 )/12, which has been verified by the use of two different symbolic computations.

Parametric study of the bearing motion dynamics
This section investigates a parameter study of the dynamical behaviour of a bearing with slip effects for large amplitudes of the prescribed rotor oscillations and with the stator free to move axially in response to the fluid film dynamics; dynamical results for a bearing with a no-slip condition l s = 0, are given in Bailey et al. [8]. Figure 4 shows the periodic fluid stiffness, total damping, force, stator height and MFC for varying coning angle under external pressurisation in the case of a wide annulus 1 − a = 0.8 with parameter values l s = 0.1, = 1.2, λ = 1, α = 1, K z = 10, σ = 1 and D a = 1. An NCB has smaller g min than a PCB, with all cases having a very small MFC which is approximately constant during the time interval 1 < t < 2.5, where the total damping and fluid stiffness have a local maximum and extrema, respectively. The fluid stiffness has a maximum of small magnitude for a PCB and minimum of large magnitude for an NCB. The damping has a maximum, with largest amplitude in the case of an NCB. We note that the stator has larger amplitudes of oscillation for a NCB than a PCB of the same value of coning angle and the force extrema occurring when the MFC is small, with an NCB having larger oscillations than a PCB. Increasing the bearing parameters, such as the amplitude of the rotor oscillations, slip length or bearing width give appropriate modified values but overall similar dynamics.
A main focus of this work is to examine the MFC g min with a selection of values for the main bearing characteristics of parameters coning angle, slip length, speed parameter. Table 1 shows the obtained values of g min for a wide a = 0.2 and narrow a = 0.8 annulus in the case of external pressurisation and large amplitude rotor oscillations, = 1.2. A corresponding table for internal pressurisation is given in Appendix 4, Table 4. Under all pressurisations considered, values of g min are larger in a bearing with a wide annulus than with a narrow annulus for a given bearing configuration. An exception is the case of a PCB under internal pressurisation with large coning angle and speed parameter, as is the case of β = 0.3, l s = 0 and λ = 5 in Table 4. In this last case, for sufficiently large coning angle and speed parameter values, and since inertial effects are more significant in a wide annulus than a narrow one, it follows that g min is smaller in the case of a wide annulus. In all cases, increasing the magnitude of coning angle and slip length decreases the value of g min ; face contact is possible for some parameter choices, denoted by * . Under external pressurisation a PCB generally has larger values of g min than a NCB for a given set of parameters, an exception is a wide annulus with large speed parameter is considered, as shown in Table 1 for λ = 5. In this case, the increased inertia effects dominates over the differential pressure, giving larger values of g min for an NCB. This Entries denoted * correspond to predicted contact effect is minimal in a narrow annulus. In contrast with the case of external pressurisation, a bearing under internal pressurisation gives an NCB having a larger value of g min than a PCB. Increasing the speed parameter for a narrow annulus causes an NCB to have increasing g min , while a PCB has the situation reversed. In the case of a wide annulus, a similar trend occurs. However, an NCB with small angle has decreasing g min for the increasing speed parameter; this is italicised in Tables 1 and 4. Results are consistent with the expectation that, for the decreasing coning angle, both PCB and NCB should have common dynamics. A bearing with narrow annulus follows a similar trend to the wide annulus, but inertia effects are less pronounced. The change in behaviour with the increasing speed parameter occurs for smaller coning angle than those shown in Table 1.

Influence of slip effects on g min
An important aspect in bearing operation is the impact of a slip effect on maintaining a face clearance when the rotor is subjected to axial disturbances. Bearing gap behaviour has been investigated by Bailey et al. [8] to identify no contact occurs with a no-slip condition for amplitudes of the axial rotor displacement similar to or larger than the equilibrium MFC. Additional consideration is given to the dynamics of a parallel bearing to investigate if contact occurs for a wide range of slip lengths. Figure 5 shows a log-log plot of g min against increasing slip length within the range 10 −4 ≤ l s ≤ 10 6 , relative to a typical bearing gap, for axial rotor oscillations of amplitude 1.0 ≤ ≤ 1.2. The bearing has two distinctly different asymptotic behaviours for increasing slip length depending on the magnitude of axial rotor oscillations . In region I, < 1.05325, g min approaches a constant value for increasing slip length values. In region II > 1.05325, the plot of g min decreases approximately linearly as the slip length increases with gradient −1, giving contact occurring only in the limit of l s → ∞. The values for g min are obtained from the stroboscopic map solver at each value of the rotor oscillation , with increasing discrete values of the slip length between 10 −4 ≤ l s ≤ 10 6 . Investigation is extended to evaluate if, in general, contact can occur with a slip effect. Figure 6 compares g min against slip length for varying coning angle −0.2 ≤ β ≤ 0.2 in the case of a wide and narrow bearing under external pressurisation. A narrow bearing has contact in all cases, with an NCB having contact at smaller slip lengths than a PCB of the same angle, and coning angles of larger magnitude give contact at  Tables 2 and 3 for external and internal pressurisations, respectively. For values of l s larger than these values, contact will always occur. It is interesting to observe that a bearing under external pressurisation can have no contact for a PCB in the case of a wide annulus, denoted by o in Table 2. In all the cases considered, a bearing with wide annulus has contact at larger slip lengths than the one with a narrow annulus. An exception is a wide PCB under internal pressurisation with large speed parameter, italicised in Table 3. In this case, the large inertia effect is more significant in a wide annulus than a narrow annulus, giving g min smaller in a wide annulus for a sufficiently large coning angle. Increasing the magnitude of the coning angle decreases the values of slip length at contact and increasing the speed parameter gives the increase in the slip length value at contact in an NCB, but a decreased one in the case of a PCB.
As expected in the case of internal pressurisation, an NCB has a larger slip length value at contact than a PCB, while for external pressurisation this situation is reversed, except for the case of a wide bearing with large values of λ where inertial effects dominate the pressure gradient effects, italicised in Table 2.

Small face clearance in limit g(t) 1
To further investigate the case when the MFC is small and also verify the numerical solution, an asymptotic analysis is undertaken. At near contact, it is assumed that g β in the limit g(t) 1, so a new variable f = g/β is introduced, giving the limit to be examined for f The damping coefficient D as defined in (30) in the limit of f 1, is given for a PCB, denoted by subscript + and NCB, denoted by subscript −, by which are both of O(1) to leading order. The effective restoring force S as defined in (30), in the limit of f 1 to leading order is given by In the above,L and L −I are given in Appendix 3, by Eqs. (81) and (82), respectively, which are of O(1), giving the effective restoring force of O(1).
In the limit of small face clearance g/β 1, the nonlinear second-order non-autonomous ordinary differential equation (29) becomes a linear ordinary differential equation with constant coefficients, given to leading order by Solving Eq. (59) in the limit g/β 1 with initial conditions g(t 0 ) = g 0 and dg(t 0 )/dt = −ξ , with ξ > 0 say, gives with constants of integration The dependence of l s on the MFC is given by the coefficients D ± and A ± only. An NCB under internal pressurisation is examined since a PCB under external pressurisation does not have contact for large amplitude and slip lengths l s ≤ 4.5, see Fig. 6b. Rotor amplitude = 1.2 and slip length l s = 0.13 are used, giving a close comparison of the numerical solution and the asymptotic composite solution of the MFC as shown in Fig. 7, in the region 0.8 < t < 2.2. The composite solution has increasing discrepancy with the numerical solution as the gap grows, consistent with the asymptotics approximations; region III was not included in the analysis.
On the other hand, in the limit of increasing l s an asymptotic evaluation for g(t) 1 gives the same form as the asymptotic equation in (59). In this case, the damping coefficient D and the effective restoring force S are given, respectively, as to leading order. In this limit, the solution takes the same form as Eq. (60) with constants given by (61), and using D ± = D a and modified values for A ± as given in (62), which are all independent of l s . In consequence in the limit, the value of the MFC cannot depend on the slip length. This behaviour is shown in Fig. 6 when g min asymptotes of to a constant value when contact does not occur.

Summary and conclusions
Previous work by Bailey et al. [8] for a coned thrust bearing showed for axial disturbances a fluid film gap was maintained when subject to a no-slip condition, although this gap may become extremely small. This work identifies for a more general case, where a slip flow is relevant, that contact may occur, and for more practical purposes, maintaining a minimum gap thickness will depend on the coning angle, slip length, amplitude of rotor oscillations disturbance amplitude and direction of pressurisation.
A modified Reynolds equation is derived for a bearing containing a rigid coned rotor and a slip boundary condition imposed on the bearing faces characterised by a slip length, l s . An axisymmetric lubrication approximation is used and leading-order effects of centrifugal inertia are retained, which is relevant for very high-speed bearing rotation and flows that are dynamically incompressible. Axial structural response of the rotor are formulated corresponding to the stator motion, induced by periodic axial displacement of the rotor, through the fluid flow, acting as a springmass-damper structure.
The fully coupled unsteady bearing motion is examined where the stator displacement is coupled to the film dynamics via the axial force exerted on the stator by the fluid film. Rewriting the modified Reynolds equation and the stator equation in terms of the time-dependent MFC g(t) allows for explicit analytic expressions for the pressure field and force on the stator. Bearing solutions are obtained from solving a nonlinear second-order non-autonomous ordinary differential equation using a stroboscopic map to find periodic solutions. Subsequent calculations give the force on the stator, stator height, total damping and fluid stiffness.
A steady bearing case with axially fixed rotor and stator was investigated to examine the existence of nonmonotonic behaviour in the pressure field for increasing slip length. Results tending to the limit of large slip length are given for representative bearing parameter choices, showing an externally pressurised bearing can exist with no overall flux through the bearing at a critical speed parameter, corresponding physically to pressure effects directly balancing inertial effects.
The dynamics of the bearing are investigated for rotor displacement, in particular with amplitude larger than the equilibrium gap between a coupled rotor and stator system, simulating possible destabilising behaviour. Results show that if the rotor is forced into close proximity of the stator and a very thin fluid film is maintained, as the rotor continues on its forced path, the stator follows the rotor closely for a significant part of the periodic motion, preserving the fluid film. Dynamically the total damping, fluid stiffness and force are effectively constant for most of the time period, unless the MFC becomes small. Here, an enhanced localised force and a maximum in the total damping resist plate contact. Numerical solution of bearing dynamics is supported by asymptotic analysis.
Evaluation of possible contact occurring in a bearing is carried out when the rotor is subject to large axial disturbances. Increasing the coning angle decreases the value of slip length when contact first occurs, whereas increasing the speed parameter gives the slip length values at contact increasing for NCB but decreasing for a PCB. In general, a wide annulus has contact at larger slip lengths than a narrow annulus, with a wide annulus under external pressurisation having no contact for certain bearing configurations. − 3(z − h r ) 2 h − 6(z − h r )h 2 ) + 10l 3 s (−6(z − h r )h)) respectively. The radial flux through the bearing is given by The steady-state streamfunction is found using Eqs. (63) and (65). The film thickness in terms of the variable g(t) is given by h(r, β, t) = g(t) + (r − a)β for β > 0, h(r, β, t) = g(t) + (r − 1)β for β < 0.