Evaluation of the minimum face clearance of a high-speed gas-lubricated bearing with Navier slip boundary conditions under random excitations

Motivated by ongoing developments in aero-engine technology, a model for a coupled gas-lubricated bearing is developed in terms of an extended dynamical system. A slip boundary condition, characterised by a slip length, is incorporated on the bearing faces which can be relevant for operation in non-ideal extreme conditions, notably where external vibrations or disturbances could destabilise the bearing. A modiﬁed Reynolds equation is formulated to model the gas ﬂow, retaining the effects of centrifugal inertia which is increasingly important for high-speed operation, and is coupled to the structural equations; spring-mass-damper systems model the axial stator and rotor displacements. A novel model is developed corresponding to a bearing experiencing an external random force to evaluate the resulting induced displacements of the bearing components. The minimum face clearance is obtained from a mapping solver for the modiﬁed Reynolds equation and structural equations simultaneously. In the case of random excitations, the solver is combined with a Monte Carlo technique. Evaluation of the average value of the minimum gap and the probability of the gap reaching a prescribed tolerance are provided. Extensive insight is given on the effect of key bearing parameters on the corresponding bearing dynamics.


Introduction
Fluid-lubricated bearings typically utilise a thin fluid film to maintain a clearance between rotating and stationary components when an axial load is imposed. Increasingly industrial applications of bearings are characterised by high operating speeds, requiring low frictional losses and smaller clearances. Associated film lubrication technology is designed to enable further improvements in bearing efficiency. In this high-speed and small-clearance scenario, comprehensive and accurate predictions of variation in the face clearance are needed, especially when the bearing experiences a range of external loading and disturbances. Quantifying the possibility of contact between the bearing faces is essential.
The dynamics of an air bearing slider, which employs a thin lubricating air film to separate two non-parallel moving plates, has been studied in [1] when the bearing has relative normal motion between the two plates. The interaction of a compressible gas flow with a movable rigid surface was presented, together with linear stability analysis, showing separation of the plates can be maintained. Analytical methods were used in [2] to study the dynamics and stability of tapered air bearing sliders, providing insight into the slider behaviour.
Taylor and Saffman [3] performed an early theoretical study examining a coaxial parallel rotor and stator separated by a thin air film, experiencing normal motion. A model was derived from the compressible Navier-Stokes equations for negligible rotational inertia effects and small amplitude disturbances. For detailed predictions of the bearing dynamics, a coupled bearing model in which the fluid flow is coupled to the bearing structure is required, such that the air film is an integrated component of the bearing [4]. Etison [5] examined a coupled fluid-lubricated device in three dimensions, identifying hydrodynamic and hydrostatic components of the air film pressure; the hydrostatic component includes the squeeze film pressure. Results indicate that the squeeze film behaviour has a significant contribution to the dynamics and can potentially maintain the air film alone. For highly vibrating operational environments, Salbu [6] examined the effect of significant disturbances in the axial direction through theoretical analysis and experimental investigations. Describing the rotor-stator clearance by an oscillatory motion gives predictions of increasing pressure and force in the air film as the amplitude of axial oscillations increases.
The bearing gap dynamics in a coupled high-speed air-lubricated bearing with parallel faces were studied by Garratt et al. [7]. In this case, the axial displacement of one face has periodic axial oscillations of smaller amplitude than the equilibrium film thickness, and the other face moves axially in response to the film dynamics. Bailey et al. [8] considered a similar bearing configuration with the inclusion of a coned rotor and incompressible flow. In this case, an analytical solution of the corresponding Reynolds equation can be found, leading to an explicit analytical expression for the fluid force. The coupled bearing problem then reduces to a second-order non-autonomous ordinary differential equation for the bearing gap. For compressible flow, this direct approach is not feasible due to the nonlinearity of the resulting Reynolds equation. The effect of the rotor undergoing prescribed axial oscillations with amplitude larger than the equilibrium film thickness was examined. Results indicate that the fluid film can become very small, possibly invalidating the classical no-slip velocity condition. Comparison of a CFD model of an air-riding seal with experimental studies by Sayma et al. [9] identified that steady-state flow solutions could not be obtained for typical gaps of under 10 μm indicating more analysis is required at this smaller scale to comprehend the corresponding seal dynamics.
The reduced scale of micro-and nanofluidics devices results in the fluid flow being characterised by the confined fluid at micro-and nanoscales, where surface effects begin to dominate over volume-related phenomena. Mathematical models of thin gas flow regimes are usually classified by the Knudsen number K n = l/ĥ 0 , where l is the mean free path (for air at atmospheric conditions l = 68 nm) andĥ 0 is the characteristic film thickness. Knudsen numbers in the range 10 −3 ≤ K n ≤ 10 −1 , represent flow at micro-and nano-scales, where a continuum flow model with slip boundary condition is usually employed. The existence of velocity slip was first predicted by Navier [10], who proposed a constant slip model with a linear relationship between the fluid-wall velocity difference and the tangential shear rate, i.e. the velocity slip is proportional to the derivatives of the surface fluid velocity, with the slip length as the proportionality constant (a first-order model). Using Maxwell's [11] classical theory, the slip coefficient is taken as proportional to (2− f )/ f and is usually of the order of the mean free path; f is the momentum accommodation coefficient.
Use of a slip flow condition has been previously incorporated into models of various bearing geometries. A gaslubricated inclined plane slider bearing has been studied in the slip flow regime using a first-order slip model [12], second-order slip model [13] and modified second-order slip model incorporating additional physical considerations, referred to as a 1.5 slip model [14]. For each case, the corresponding Reynolds equation for compressible flow was developed and analytical results were obtained.
A non-axisymmetric thrust bearing with slip flow and foil pads on the rotor face is examined by Park et al. [15], where the gas flow modelled by a classical Reynolds equation is coupled to the bearing structure. Bearing dynamic predictions are studied for small amplitude rotor displacement using perturbation analysis in the case of a no-slip and slip condition imposed on the bearing faces. A slip condition results in a reduced hydrodynamic pressure and thus load carrying capacity when compared with a no-slip condition, as well as decreased stiffness and damping coefficients in the case of axial perturbations. Incorporating a slip boundary condition in the analysis of a coupled high-speed fluid-lubricated parallel bearing with small face clearance in the limit of incompressible flow [16] shows that face contact does not formally occur when the rotor undergoes prescribed periodic oscillations but the bearing gap can become very small. An increased lubrication force on the stator can be attained by the use of a coned rotor with stability studies identifying an optimal coning angle [5]. However, detailed analysis of a coned bearing in the slip regime under extreme operating conditions gives predictions of possible face contact [17]. A similar trend is also identified for the case of a coupled high-speed fluid-lubricated bearing in the case of compressible flow with slip boundary conditions under extreme operating conditions [18]. Possible external disturbances imposed on the rotor can cause the rotor to be displaced by magnitudes larger than the initial equilibrium film thickness. Detailed results show a parallel rotor and stator maintain a face clearance, although it can become very small, where as a coned bearing configuration permits contact between the rotor and stator.
For increasingly more compact designs at higher rotational speeds, there exists a requirement for greater predictability of the dynamic behaviour in industrial bearing designs. Sources of uncertainty in a bearing model typically arise from the random axial motion of the rotor due to external excitations. In [16] and [17], the effect of external excitations on the bearing is studied by prescribing the axial motion of the rotor as periodic oscillations and utilising a comprehensive deterministic mathematical model for the bearing dynamics. The probability of bearing face contact is correspondingly examined using the method of derived distributions to incorporate the lack of precise knowledge in the amplitude of rotor oscillations [19]. The probability density function (transforming functions) and cumulative distribution function of the minimum bearing gap are readily determined when the amplitude of rotor oscillations is described by a specified probability distribution function. Utilising the deterministic relationship between the amplitude of rotor oscillations and minimum value of the bearing gap subsequently leads to the evaluation of the probability of contact.
Analysis of dynamical systems and bearing models with uncertainty arising from external excitations use computational methods which are typically based on Monte Carlo techniques [20]. The Monte Carlo method samples directly from the random input parameter and the deterministic model is evaluated for each sample, giving the corresponding output. This robust procedure is able to evaluate complex cases and a large number of uncertain parameters; however, an extensive number of code executions may be needed to generate sufficiently accurate results. To minimise the computational cost and increase efficiency, other computational approaches for uncertainty quantification have been developed; examples include polynomial chaos, Gaussian process emulator and a Quasi-Monte Carlo method, see, e.g. [21][22][23].
The dynamic response and vibration characteristics of two in-line rotor bearing systems are examined by Li et al. [24] to quantify the effects of uncertainty. A model for the rotor system is derived using stochastic modelling based on the polynomial chaos expansion technique for uncertainty in the damping and nonlinear support stiffness. The results obtained are reported to have good agreement with Monte Carlo simulation. Response statistics demonstrate that uncertainty in the nonlinear support stiffness and damping has a significant effect on the predicted dynamic response of the rotor system.
Increasingly environmental and efficiency considerations are leading to the next generation of bearing and sealing technologies having an extremely small gap, resulting in an almost contact design. In this study, we outline a model for a gas-lubricated bearing incorporating a slip boundary condition. The compressible flow model is coupled to the structural dynamics through the axial force of the fluid on the bearing faces. Previous analysis is extended to incorporate operation in a non-ideal environment; external vibrations or disturbances could act to destabilise the bearing which are modelled as a random external force applied to the rotor. In Sect. 2, the coupled governing equations are formulated where a Reynolds equation for compressible flow is derived that retains the effects of centrifugal inertia for high-speed operation and incorporates slip boundary conditions, characterised by a slip length parameter, on the bearing faces. The axial stator and rotor displacements are modelled by spring-mass-damper systems. Of significant importance is the analysis of a random external force imposed on the rotor, modelled as a continuous, stationary, bounded process with finite modes of frequencies in its spectrum. The generality of the proposed random force is potentially transferable to other models of physical phenomena beyond this particular application. To numerically solve the coupled bearing model, the modified Reynolds equation is solved simultaneously with the structural equations via a mapping solver, details of which are given in Sect. 3. The Monte Carlo technique is implemented to determine the cumulative distribution density function of the minimum gap (face clearance) due to the imposed random destabilising external force, and also to find the average of the minimum gap and probability that it reaches a prescribed tolerance. Results are presented in Sect. 4 examining the effect of key model parameters.

Mathematical model
A simplified mathematical model is derived for a gas-lubricated bearing with a thin fluid film modelled as compressible flow and a slip boundary condition imposed on the bearing faces. The coaxial, axisymmetric annular rotor and stator can move axially, giving the respective displacement heights asĥ r andĥ s . The rotor also has a prescribed rotational speedΩ and the bearing has imposed pressures at the inner and outer radius. This configuration is motivated by aerospace applications and similar technology in the literature is referred to as air-riding [7], film-riding [25], oil-free [26] or noncontacting [4].
The governing Navier-Stokes equations and boundary conditions for the gas flow in the bearing are expressed in dimensionless variables by taking a typical bearing radiusr 0 and film thicknessĥ 0 . The characteristic rotor velocity is taken asΩr , typical air density is taken to beρ 0 and time is scaled byT defined by the time scale of an imposed force. The dimensionless velocities are defined asû/Û ,v/(Ωr 0 ) andŵ/(ĥ 0T −1 ), where (û,v,ŵ) are the dimensional radial, azimuthal and axial velocities, respectively. The dimensionless radius and height are given by r =r /r 0 and z =ẑ/ĥ 0 , respectively. The dimensionless density is ρ =ρ/ρ 0 and dimensionless slip length is l s =l s /ĥ 0 . The bearing geometry for parallel faces is shown in Fig. 1 in the dimensionless coordinate system (r, θ, z); corresponding geometry for a coned rotor is given in Fig. 2. The coning angle β is assumed small, such that sinβ = O(δ 0 ) and cosβ = 1 + O(δ 0 ). Correspondingly, a scalingβ = βδ 0 with β = O(1) gives consistency with the lubrication condition such that ∂h r /∂r = −β. Figure 3 shows a cross sectional view of a negatively coned bearing (NCB) and positively coned bearing (PCB). The dimensionless coned rotor and plane stator heights are given by h R (r, t; β) and h s (t), respectively, and the axisymmetric rotor-stator clearance is defined by h(r, as shown in Fig. 3; h r (t) is defined at the outer radius for a NCB and inner radius for a PCB. The inner radius is rescaled to a =r I /r O and the outer radius to 1, with corresponding dimensionless pressure boundary conditions This study considers explicitly the coning angle arising from deformation of the rotor due to over pressurisation of the bearing, resulting in a NCB having external pressurisation ( p I < p O ) and a PCB having internal pressurisation ( p I > p O ). In both cases, this configuration results in a pressure gradient corresponding to a diverging channel.

123
Evaluation of the minimum face clearance of a high-speed… The rotor experiences an external axial force N (t), which is assumed to be axisymmetric and imposed at r = a. It is taken to be a random process, which represents external excitations on the bearing that could act to destabilise the bearing operation. For the case of non-axisymmetric forcing, for example, point loading, the bearing faces may become misaligned and experience tilt and/or swash. In turn, this may lead to coning angle instabilties over time. A dynamic analysis for a noncontacting face seal with coned faces has been examined in three dimensions under the assumption of assembly face misalignment and rotor axial run-out. Optimal conditions to minimise face contact and maximise the stability are addressed in [27][28][29].
A key quantity in the bearing dynamics is the time-dependent minimum face clearance (MFC), defined by the variable g(t) = h s (t) − h r (t), evaluated at the outer radius for a NCB and inner radius for a PCB; for a parallel

Air flow model
The air flow model is derived from the compressible Navier-Stokes momentum equations and conservation of mass equation for a compressible flow in axisymmetric coordinates; the flow is assumed to be an isothermal ideal gas. Thermal effects in similar configurations have been studied in the absence of external forcing by numerous authors, for example see [30,31], suggesting thermal deformation of the faces is more significant than viscous heating. A study by Minikes and Bucher [32] concluded that temperature variations were less than five per-cent of the ambient temperature and consequently the flow could be an isothermal process. Correspondingly, the thermal effects are neglected in the paper as they are secondary to the effects of the dynamic motion. Under the above approximation the radial and azimuthal Reynolds numbers and the Reynolds number ratio Re * are given respectively as where ν = μ/ρ 0 is the kinematic viscosity. The aspect ratio δ 0 =ĥ 0 /r 0 is small for thin film bearings δ 0 1, justifying the use of the lubrication approximation. The pressure is scaled as P = μr 0Û /ĥ 2 0 to ensure that the pressure gradient in the radial direction is retained at leading order and the Froude number Fr is defined as Fr =Û (gĥ 0 ) −1/2 , where g denotes the acceleration due to gravity and parametrises the importance of gravitational effects relative to the radial flow. Gravity is neglected because taking the Froude number of O(1), ensuring consistency with lubrication theory, leads to Re U δ 0 Fr −2 1. Classical lubrication theory neglects inertia when the reduced Reynolds number Re U δ 0 1; however, in the case of high-speed bearing operation, an additional term corresponding to the ratio of the Reynolds numbers Re * must be considered, which is not always negligible. A measure of the importance of the viscous time scale is given by gives the time scale of the order of the time it takes viscosity effects to diffuse over the characteristic film thickness. Thus giving the acceleration term of the same order as the viscous term. In the current study, the time scale is taken to be of O(ĥ 0r0 /ν), resulting in κ = O(δ 0 ), justifying a quasi-static approximation of the momentum equations.
Applying a lubrication condition and the given assumptions to the compressible Navier-Stokes momentum equations, results in the leading-order momentum equations A speed parameter λ = Re U δ 0 Re * 2 =r 0ĥ 2 0Ω 2 /(νÛ ) characterises effects of inertia due to rotation. If λ = 0, the standard lubrication equations for compressible flow in axisymmetric cylindrical coordinates are recovered.
Similarly, the conservation of mass equation and equation of state become where σ = κ/Re U δ 0 =r 0 /ÛT . The dimensionless ideal gas constant in (5) is given by K s = Rτĥ 2 0 /μr 0Û M, which relates the pressure field to the density field, where R is the gas constant and M is the molar mass.
A first-order Navier formulation slip model is exploited, as used in Bailey et al. [16], where a slip velocity is introduced that is proportional to the tangential component of the surface traction. The velocity boundary conditions consist of components of the fluid which are continuous across the fluid-solid boundary with a modification due to slip flow induced by the wall shear, together with a no flux condition. Thus the dimensionless velocity boundary conditions, to leading order, are taken as The fluid velocity components tangential to the wall in (6) are proportional to the relevant wall shear stress where the proportionality constant l s , denotes a dimensionless slip length. The limit l s = 0 corresponds to no-slip conditions. A derivation of the slip velocity conditions is given in Bailey et al. [17]. The air flow velocities are readily determined from the leading order Navier-Stokes equations (4). Integrating the conservation of mass equation (5) between the rotor and stator, applying the Leibniz integral rule and the velocity boundary conditions (6), gives the modified nonlinear Reynolds equation for compressible flow as The modified Reynolds equation expresses the relationship between the film pressure and rotor-stator clearance h, subsequently yielding the flow characteristics when solved. Dependence on the coning angle is given implicitly in the rotor-stator clearance h. Taking λ = 0 results in the centrifugal inertia effects being neglected;however, the rotor and stator still experience relative rotational motion due to the velocity boundary conditions.

Structural dynamics
The axial displacements of the stator and rotor are modelled using a standard spring-mass-damper model incorporating the bearing pressure variation and an external random, time-dependent force imposed on the rotor N (t). In dimensionless form, the stator equation is and the rotor equation is either The rotor equations (9) are defined at r = a, resulting in separate equations for cases of a NCB and PCB, to ensure that they have a similar volume of gas in the bearing for a given magnitude of coning angle. In this configuration, a PCB and a NCB have different equilibrium rotor heights at r = a.
The corresponding fluid axial force of the fluid on the bearing faces is given by Here p a =p a /P is the dimensionless atmospheric pressure, α = νÛ /mσ 2 δ 3 0 is the dimensionless force coupling parameter,m is the mass of a face, D a =D a /mσ is the dimensionless damping parameter and K z =K z /mσ 2 is the effective restoring force parameters.
In a physical application, random fluctuations reflect the uncertainty in the precise state of a system. In the present work, the bearing is envisaged as part of a complex dynamical system characterised by a periodic motion of periodT ; this will be the case for an aero-turbine engine. The system is composed of several rotating parts where the motion of each individual part and its interaction with the surrounding environment induce an external force on the rotor surface of the bearing. As the external force is not known with any precision, it is considered to be a random force. Due to restrictions on the mathematical formulation of the problem as well as physical and operating constraints, the external force is taken to be continuous, stationary, having a bounded amplitude, and a finite number of frequencies in its spectrum. Correspondingly, the random external force is chosen of the form and amplitudes A i are truncated normal random variables given by where X i are standard normal random variables; the phases φ i are uniformly distributed random variables on [0, 2π ]; and the frequencies B i are taken as discrete random variables uniformly distributed on the finite set of non-negative integers: {0, 1, . . . ,B}. The choice ofB is motivated by the need for the time scale of the forcing to be comparable or longer than the viscous time scale. The random variables X i , B i , φ i , i = 1, . . . , n are mutually independent. The form of the random external force N (t) contrasts with Langevin-type equations associated with white noise models that are typically unbounded, discontinuous and aperiodic [33,34]. Boundedness of the noise and frequency in the current work, together with continuity, is important for a realistic model of an external force in aero-applications. Correspondingly, a realisation of N (t) from (11), determined from a realisation of A i , B i , φ i , provides a periodic function of time. Noting that φ i is a uniform random variable on [0, 2π ], then N (t) is a stationary stochastic process [35,36]. Since A i are bounded random variables and n is finite, N (t) is bounded which is a natural requirement for modelling a physical force. The properties of the proposed random force (11), namely that it is continuous, stationary, boundedness in physical space and has a finite number of frequencies in its spectrum, is potentially transferable to other models; the proposed random force may be useful when modelling other physical phenomena which experience random fluctuations.

Numerical methods
Consideration of a fully coupled bearing model provides a quantitative study of the interaction between the air film and the rotor-stator structure of fluid-lubricated technology. Investigations into the operational requirements for dynamic bearing behaviour are undertaken in the presence of an external force subject to uncertainties on the rotor. To numerically solve the coupled bearing model requires the modified Reynolds equation (7) to be simulated simultaneously with the structural equations (8)- (9). Equations are derived for a NCB, and similar methodolgy can be followed for those for a PCB, see Bailey et al. [18]. Substituting the MFC g(t) = h s (t) − h r (t) into the modified Reynolds equation (7) gives for a NCB. The equation for a parallel bearing is obtained from (13) by taking β = 0. Re-writing the stator equation (8) in terms of the MFC gives the system of structural equations for a NCB or parallel bearing as To couple numerically the structural and fluid dynamics, the Reynolds equation is discretised in the spatial variable and approximated by a second-order central finite difference scheme. Solutions to the coupled system, Reynolds equation (13) and structural equations (14), are denoted by the vector g(g(t), h r (t), Z (t), Y (t), p i (t)), where p i denotes the pressure at radial points r i ; a ≤ r i ≤ 1, with given initial conditions: In this way, the problem reduces to the following system of five first-order random ordinary differential equations. For a NCB, the set of equations is given by The Reynolds equation is discretised in the spatial variable and approximated by a second-order central finite difference scheme for i = 2 : ξ −1, where ξ is the number of radial points used in the finite difference approximation of equation (13). The boundary conditions for the pressure are given by p 1 = p I and p ξ = p O . The force of the fluid on the faces is denoted by where w i is a weighting function. The set of equations (13) are solved using a stiff ordinary differential equation solver, namely a variable-step, variable-order solver based on the numerical differentiation formulas of orders 1 to 5.
For every realisation of the random force (11), the rotor experiences a periodic force which has a period of T = 2π . Periodic solutions of the system of equations (16), with fixed A i , B i , φ i in (11), are sought; the solver used only converges to periodic solutions. Therefore, a Re 2 → Re 2 map is formulated, advancing an initial condition g 0 at time t 0 by the time T = 2π , defining a stroboscopic map (T ; g 0 , t 0 ) which integrates the system of equations (16) forward through one period of the forcing. A variable-order solver based on the numerical differentiation formulas is implemented due to the equations becoming locally stiff when the face clearance becomes very small. Periodic solutions are identified via the fixed points of the stroboscopic map g(t) = g(t + T ) and found iteratively, corresponding to the condition giving periodic solutions g(g(t), h r (t), Z (t), Y (t), p i (t)). Due to the nonlinear character of the system of equations, an iterative Newton's method is used to find solutions numerically, given an initial guess valueg 0 . Successively improved iterates of the initial guess g 0 are computed by the numerical iterative scheme with the Jacobian matrix The elements of the Jacobian matrix J (T ) are found by solving similar systems to the auxiliary system of equations (16), but which are the derivative with respect to each of the initial guess valueg 0 ; corresponding initial conditions are used. An alternative method, as the modified Reynolds equation (13) does not depend explicitly on the rotor height h r , is to solve directly for the pressure and face clearance only and find the rotor height in post processing using F(t) and N (t). This will eliminate the second and fourth equations in (16) and reduce the Jacobian (19) by removing the second and fourth rows and columns. However, if the number of radial points ξ is large, the reduction will be marginal.
The stroboscopic map is combined with the Monte Carlo technique, where a large number of independent random realisations of the stroboscopic map are computed, each with a randomly generated external force imposed on the rotor.
Initially, the stroboscopic map is used to obtain the approximationḡ min of the minimum face clearance g min as Then the approximate average of g min can be computed as whereĝ min is the approximate average ofḡ min andḡ (m) min is the value ofḡ min for independent realisation m; m = 1, . . . , M. The Monte Carlo error is where var(ḡ min ) is estimated as and E[ḡ min ] ∈ ĝ min − MC err ,ĝ min + MC err (23) with a probability of 95% (confidence interval). To achieve a sufficiently accurate result, the Monte Carlo error was made sufficiently small through choosing an appropriately large number M of independent realisations. We note that the estimateĝ min has bias E[g min ] − E[ḡ min ], corresponding to the error of the numerical approximationḡ(t) of the stroboscopic map g(t). The numerical solution was compared to an analytical asymptotic solution in the case of small face clearance in Bailey et al. [8]. Results showed the proposed numerical solution has negligible error of the minimum face clearance compared to the asymptotic solution, indicating that the bias of the estimateĝ min can be considered to be negligible. An estimate for the variance of g min is computed as in (22), and its Monte Carlo error is evaluated as

Results and discussion
A coupled gas bearing is examined to investigate the effects of key parameters of the bearing dynamics and to identify possible face contact corresponding to a generic aero-engine industrial application. The bearing configuration is taken to be a narrow bearing a = 0.8 with σ = 0.821 and ambient pressure imposed at the outer radius p O = p a = 1. Initially ambient pressure is considered at the inner radius p I = 1, the bearing faces taken to be parallel β = 0 and the external force (11) has parameters n = 20,B = 10, andĀ = 20; the effect of these parameters on potential face contact is examined. The speed parameter is considered to be λ = 0.0029, dimensionless slip length to be l s = 0.1 and the ideal gas constant as K s = 1. The structural parameters are given by K z = 12.2, D a = 1.5 and the force coupling parameter by α = 1.22.
As examples, two different realisations of the bearing dynamics are given in Fig. 4 showing the pressure in the fluid film, force of the fluid on the bearing faces, the stator and rotor heights, face clearance and external force imposed on the rotor. It is noted that the realisation in Fig. 4b is an extreme realisation where the face clearance becomes very small, and a more typical representation of the evolution of the system is given in Fig. 4a. The external axial force on the rotor causes axial displacement of the rotor: however, the stator only has axial displacement if the rotor displacements are sufficiently large such that the face clearance becomes sufficiently small. When the rotor is in close proximity of the stator, resulting in a very small face clearance, as shown in Fig. 4b, there is a peak in the fluid force and pressure of the fluid ensuring the fluid film is maintained through axial displacement of the stator. The stator oscillations decrease over time unless the rotor comes into close proximity again (as in Fig. 4b), where there is another peak in the fluid force and further axial displacement of the stator. The smaller the face clearance becomes, the larger the peaks in the fluid force and pressure in the fluid film become. The mean value and variance of the minimum face clearance are computed for an increasing number of Monte Carlo simulations and are given in Fig. 5 for a parallel bearing together with the corresponding Monte Carlo errors. Increasing the number of Monte Carlo realisations gives converging values for the mean and variance of the minimum face clearance of μ = 0.447 ± 0.00340 and σ 2 = 0.0432 ± 0.000885. Increasing the number of realisation from M = 1000 to M = 4000 effectively halves the Monte Carlo error for mean value from 3.00% to 1.48%. Increasing the number of realisations further to M = 15, 000 reduces the Monte Carlo error to 0.760%, signifying a balance between the accuracy of the results and the excessive computation requirement that is desirable. For all subsequent cases considered, the number of realisations will be chosen such that the Monte Carlo error is ≤ 2.5% for the mean value of g min to limit computation time without a significant loss in accuracy. A histogram for g min is presented in Fig. 6 for the case of M = 15, 000 in Fig. 5 showing the frequency distribution. The realisations give a right-skewed distribution due to the dynamics of the bearing when the face clearance becomes very small, i.e. in the limit of contact. The average value of g min is 0.447 ± 0.00340; however, there is a non-negligible possibility of the face clearance reaching a small prescribed tolerance, with the limiting case of contact g min = 0.
Evaluation of the parallel bearing is examined for increasing internal pressure imposed at the inner radius; for p I = 1 an ambient pressure is imposed across the bearing, p I = 0.5 corresponds to external pressurisation and p I = 1.5 to internal pressurisation. The mean and variance of g min are given in Fig. 7, together with the corresponding Monte Carlo error, for a parallel bearing with increasing pressure at the inner radius p I . In general, as the pressure at the inner radius is increased, the average value of g min increases. It is noted that the case of p I = 1.1 has a value less than p I = 1, however with the Monte Carlo error this case fits with the overall trend of the mean value of g min increasing as the value of p I increases. In all cases, the variance of g min is small; σ 2 is of the order of 10 −2 .
The corresponding probability of g min reaching a prescribed tolerance g tol is given in Fig. 8, including the Monte Carlo errors. Of note is that the probability decreases as the prescribed tolerance decreases, with negligible probability of face contact, i.e. g tol ∼ 0. For each prescribed tolerance, increasing the pressure at the inner radius leads to a general trend of the probability of g min reaching the value of g tol decreasing. There are larger decreases in the probability of g min reaching a prescribed tolerance for large values of g tol . Figure 9 gives the mean and variance of g min with the corresponding Monte Carlo errors in the case of increasing values for the truncated amplitude of the components of the random force. Corresponding values of the probability of g min reaching a prescribed tolerance of decreasing value g tol for increasingĀ are given in Fig. 10. The boundĀ of the truncated random variable A i models the amplitude of the individual components of the random force N (t); the bound is finite due to the structural limitation of the bearing geometry and its local environment. To test our numerical solution, values 20 ≤Ā ≤ 100 are examined; the mean value of g min is similiar for all cases examined, and the variance has a small value of the order of 10 −2 in all cases. The probability of g min reaching a given tolerance g tol remains similiar when the values ofĀ increase, with negligible probability of face contact, i.e. g tol ∼ 0. The mean and variance of g min with the corresponding Monte Carlo errors are shown in Table 1 for the case of a parallel bearing with increasing number of frequencies (modes) in the forcing term (11), together with the probability of g min reaching a prescribed tolerance g tol . For an increasing number of frequencies, the number of Monte Carlo realisations required increases, and the mean value of g min decreases, while the variance increases. Decreasing the prescribed tolerance causes the probability of g min becoming that small to decrease, with negligible  probability of g min reaching zero. Increasing the number of frequencies in the forcing term leads to an increase in the probability of g min reaching the prescribed tolerance. This can be observed as the increasing the number of random models causes the mean value of the minimum gap to significantly decrease with a corresponding increase in the variance. Analysis of a coned bearing is undertaken to examine the effects of possible rotor deformation arising due to the possible over pressurisation of the inner radius of the bearing. Correspondingly, a NCB is examined under external  pressurisation; p I = 0.5 and a PCB under internal pressurisation; p I = 1.5. Table 2 gives the mean value and variance of g min (μ and σ 2 , respectively) with the corresponding Monte Carlo errors in the case of a NCB and PCB; for comparison, a parallel bearing is examined under the respective pressure boundary conditions. In general, a PCB requires less realisations than a NCB to obtain the required level of accuracy and has larger values for the mean value of g min and variance. A NCB has a monotonic decrease in the mean value of g min as the magnitude of the coning angle increases, whereas a PCB has a non-monotonic trend. The non-monotonic trend in the case of a PCB may be due to the pressure gradient across the bearing resulting in an internally pressurised bearing as opposed to a NCB being externally pressurised.
The cumulative distribution function (CDF) is examined which gives the probability that g min is less than or equal to a given value of g min , i.e. F g min (x) = Prob(g min ≤ x). The cases of β = 0 and | β |= 0.3 for a NCB and PCB are shown in Fig. 11 with 95% upper and lower confidence bounds; the probability of contact is given by the CDF when g min = 0. The probability that g min becomes equal to or less than a prescribed tolerance g tol is shown in Table  3; note that the parallel bearing is examined under different pressurisations for comparison. Both scenarios have predictions for a parallel bearing having negligible probability of contact. In contrast, a coned bearing has a non-zero probability of contact for all the cases considered. A NCB has an increase in the probability of g min reaching a given tolerance as the magnitude of the coning angle increases. However, in the case of a PCB this situation occurs only for g tol = 0.001 and g tol = 0; there is a non-monotonic trend otherwise. A NCB has a larger probability of g min reaching the prescribed tolerance than the corresponding PCB, resulting in a higher probability of contact at the outer radius of the bearing than at the inner radius (see geometry in Fig. 3). The asymmetry between the behaviour of a NCB and PCB is noted; a NCB has a monotonic trend always, whereas a PCB can have a non-monotonic trend as seen in Tables 2 and 3. This type of analysis can be used to identify safe operating conditions or limitations on bearing geometries under operating conditions.

Conclusions
A dynamical model for a gas-lubricated bearing with a very small face clearance is derived, which is capable of predicting the bearing dynamics under destabilising operating conditions. A modified compressible Reynolds equation is formulated to model a gas film, which is then coupled to the bearing structure through an axial force of the gas on the bearing faces. The governing compressible Navier-Stokes equations are reduced through using an axisymmetric lubrication approximation, but maintain leading-order effects of centrifugal inertia, relevant for highspeed flows. Modified surface boundary conditions corresponding to a slip length formulation are incorporated, relevant to the case of very small face clearance. The axial displacement of the rotor and stator are modelled as spring-mass-damper systems. Typically the rotor experiences an external random force due to fluctuations in the bearings surrounding environment. The generality of the proposed random force, which is continuous, stationary, bounded in physical space and has a finite number of frequencies in its spectrum, is potentially transferable to other models of physical phenomena beyond this particular application.
To numerically solve the coupled bearing problem, the modified Reynolds equation and structural equations are simulated simultaneously within a period. Computed results directly give the face clearance and rotor height, whilst post-processing calculations are used to determine the force of the gas on the bearing faces and axial displacement of the stator. Monte Carlo simulations are combined with a derived stroboscopic map solver to evaluate approximations of the minimum face clearance g min . The probability associated with a specified minimum face clearance g min is examined, with face contact corresponding to a zero value of g min .
Utilising the bearing dynamics allows an investigation into the effect of the external force on the face clearance. Axial displacement of the rotor is due to the imposed external force, however axial displacement of the stator arises from the rotor being forced into close proximity of it. In this case, a peak in the force of the fluid causes the faces to move apart in order to maintain a fluid film. Results are given for an increasing number of realisations in a Monte Carlo simulation, yielding the mean and variance of the minimum face clearance.
The effect of key parameters have been examined, including the pressure gradient across the bearing, value of the bounds of the random variables describing the amplitude of the random force components and the coning angle of the rotor. Outcomes of this type of analysis provide insight to a given bearing geometry to identify safe operating conditions and inform ideal bearing geometries such as in the case of establishing a maximum coning angle for safe operation under given working conditions.
give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.