Eﬀect of Uncertainty in External Forcing on a Fluid Lubricated Bearing

. Developments in industrial applications motivate improvements in ﬂuid lubricated bearing technology, enabling smaller face clearances and increased rotation rates. Associated ﬁlm lubrication technology aims to improve eﬃciency and reliability. A bearing model is developed to evaluate the eﬀect of external, potentially destabilising, random forcing applied to a pair of highly rotating axisymmetric bearing faces, separated by a thin ﬂuid ﬁlm. Two cases of random external force are examined. A ﬁrst study considers an imposed random force disturbance constrained to a ﬁxed period, where the average minimum face clearance together with the probability it reaches a speciﬁed gap tolerance. More general uncertainties are associated with more complex external forcing and takes the form of a white or coloured noise. In this case the average time for the face clearance to reach a prescribed tolerance is examined. Results can inform bearing design, providing an indication of the eﬀect of disturbances on the average lifetime and identify constraints on operating conditions for safe and reliable behaviour. technique is derived and implemented based on a simple random walk combined with a Monte Carlo technique. A detailed approach is provided for the case of a parallel bearing and the external force being taken as a Gaussian white noise source. The average time for the face clearance to reach a prescribed tolerance was examined for both white and coloured noise. Results show the force coupling parameter can have a signiﬁcant eﬀect on the average time for the face clearance to exceed a prescribed tolerance and provide an indication of the average lifetime of a bearing.


Introduction
Improvements in fluid lubricated bearing technology, also described as non contacting, air/film riding or oil-free, are motivated by ongoing developments in industrial applications. Increasingly aero-engine technology is characterised by high operating speeds, requiring low frictional losses and smaller clearances. The associated film lubrication technology requires improvements in efficiency while maintaining reliability under more demanding operating conditions. An important scenario to evaluate is the effect of external loading and disturbances that could act to destabilise the bearing.
Fluid lubricated bearings typically comprise a pair coaxial axisymmetric rotating and stationary components separated by a thin fluid film. Improvements to the bearing efficiency is associated with maintaining a thin fluid film, even when experiencing a range of external forces, in order to minimise contact between the faces. It is desirable to have comprehensive understanding of the bearing response characteristics, quantifying the possibility of contact between the bearing faces and bearing life span to identify safe and reliable operation.
Detailed predictions require a bearing model in which the fluid flow is coupled to the bearing structure [1]. Hydrodynamic and hydrostatic components of the air film pressure have been identified in a three dimensional coupled fluid lubricated device [2]; it was shown that the hydrostatic component could potentially maintain the clearance alone. The effect of axial disturbances, with significant amplitude, have been investigated by Salbu [3] using both theoretical and experimental techniques.
For a coupled high speed air lubricated bearing with parallel faces, Garratt et al. [4] examined the bearing dynamics in the case of one face having prescribed axial oscillation of amplitude less than the equilibrium film thickness and the second face being free to move axially in response; no face contact was predicted. A similar configuration was studied by Bailey et al. [5] with the inclusion of a coned rotor, incompressible flow and prescribed axial oscillations of large amplitude. Results indicate the clearance can become very small. Sayma et al. [6] identified for typical gaps of less than 10 µm steady-state flow solutions could not be obtained, indicating more analysis is required at this smaller scale.
For devices on a micro-and nanofluidics scale, surface effects start to dominate the volume related phenomena and classification of thin gas flow is through the Knudsen number. For Kn ≤ 10 −3 the normal continuum flow with a no-slip boundary condition hold, whereas for the range 10 −3 ≤ Kn ≤ 10 −1 , representing micro-and nanoscale flow, usually a continuum flow model with a slip boundary condition is used.
A non-axisymmetric coupled thrust bearing, with no-slip and a slip condition on the bearing surface, was examined by Park et al. for small amplitude rotor displacement [7], showing a slip condition reduced the load carrying capacity. A coupled high speed fluid lubricated bearing in the slip regime and incompressible flow was examined under prescribed periodic rotor oscillations for parallel faces [8] and tapered faces [9]. Under extreme operating conditions parallel faces do not give formal face contact but the gap becomes very small, whereas tapered faces give predictions of possible face contact. Similar behaviour is identified for the case of compressible flow [10].
Greater predictability in the dynamics of industrial bearing designs is required, especially as bearings are placed within larger operating environments and may therefore be subject to a broad range of more extensive and uncertain external disturbances. Correspondingly, uncertainty quantification approaches are of increased interest to improve the predictability of the dynamic behaviour in industrial bearing designs. A simplified model of the effect of periodic excitations on a fluid lubricated bearing has been studied where a random variable with prescribed mean and standard deviation was used to describe the amplitude of external force oscillations, giving predictions of the probability that the gap reaches a given tolerances [11]. Another method to analyse dynamical systems experienceing external excitations is to use computational methods based on Monte Carlo techniques [12], which sample from the random input parameter directly, then evaluates the deterministic model at each sample and a corresponding output is given. It is robust and capable of evaluating complex systems consisting a large number of uncertain parameters, however computationally a large number of realisations may be required to achieve sufficiently accurate results.
In this work, the response of a bearing is studied when its operation is subject to external disturbances which could destabilise the bearing behaviour. These are described for two characteristic types of random external forces applied to the rotor. For both, a fluid flow model is coupled to the motion of the bearing faces as given in Section 2. The first external force, examined in Section 3, has a fixed period as the bearing is envisaged to be part of a larger complex dynamical system characterised by a periodic motion; such as in the case of an aero-turbine engine. A stroboscopic map solver [8,9] is used combined with a Monte Carlo technique to identify the average minimum face clearance between the rotating and stationary faces together with the probability it reaches a given tolerance. The second external force, presented in Section 4, takes the form of a white or coloured noise, giving a more physically realistic description of general external force disturbances.The resulting model takes the form of stochastic differential equations (SDEs) and a numerical technique is implemented based on a simple random walk combined with a Monte Carlo technique. The quantity of interest is the average time for the face clearance to reach a small prescribed tolerance, giving an indication of the average lifetime of a bearing. Conclusions are contained in Section 5. A simplified mathematical model is formulated for a gas lubricated bearing, in which the fluid flow is coupled to the structural dynamics through the axial force of the fluid on the bearing faces. The rotor and stator are coaxial, axisymmetric annuli which can have an axial displacement. The rotor has a given rotational speed and a given pressure is imposed at the inner and outer radius. A parallel bearing geometry is shown in Figure 1 in a dimensionless coordinate system (r, θ, z).

Bearing Geometry and Mathematical Model
Governing equations for the fluid flow are expressed in dimensionless variables with typical air densityρ 0 , bearing radiusr 0 and film thicknessĥ 0 ; the aspect ratio is defined as δ 0 =ĥ 0 /r 0 . The characteristic rotor velocity isΩr and time is scaled byT , defining the time scale of an imposed force. Expressionsû/Û ,v/(Ωr 0 ) and w/(ĥ 0T −1 ) denote the dimensionless velocities, with (û,v,ŵ) as the corresponding dimensional velocities. Dimensionless variables for the radius, height, density and slip length are given by r =r/r 0 , z =ẑ/ĥ 0 , ρ =ρ/ρ 0 and l s =l s /ĥ 0 , respectively. The annulus inner radius is scaled to a =r I /r O with corresponding dimensionless pressure boundary condition p = p I and outer radius to 1 with p = p O .
For the case of tapered faces, the rotor has a coned angle β which is assumed to be small; sinβ = O(δ 0 ) and cosβ = 1 + O(δ 0 ). To ensure consistency with the lubrication condition, scalingβ = βδ 0 is used where β = O(1) and ∂h R /∂r = −β. A negatively coned bearing (NCB) and positively coned bearing (PCB) cross section is given in Figure 2 where the dimensionless coned rotor height h R (r, β, t) defined by The corresponding axisymmetric rotor-stator clearance is given by h(r, β, t) = h s (t)− h R (r, β, t), where h r (t) is defined at the outer radius for a NCB and inner radius for a PCB. It is assumed the tapered faces have arisen from over pressurisation of the bearing, therefore a NCB has external pressurisation (p I < p O ) and a PCB has internal pressurisation (p I > p O ); pressure gradient always corresponds to a diverging channel. 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 bearing g(t) = h(t). If the MFC always remains positive, i.e. g > 0, there is no bearing face contact and a fluid film is maintained. The leading order radial and azimuthal Reynolds numbers are given by Re U = Uĥ 0 /ν and Re Ω =r 2 0Ω /ν, respectively, with the Reynolds number ratio Re * = Ωr 0 δ 0 −1 /Û ; ν = µ/ρ 0 is the kinematic viscosity. Classical lubrication theory neglects inertia when the reduced Reynolds number Re U δ 0 1, as will be the case in the incompressible flow model. However for high-speed bearing operation, the ratio of the Reynolds numbers Re * is not always negligible and as such will be considered in the compressible flow model.

Fluid flow model
For thin film bearings, the aspect ratio is small δ 0 1, which justifies the use of a lubrication approximation. To retain the pressure gradient at leading order in the radial direction, the pressure is scaled as P = µr 0Û /ĥ 2 0 . The importance of gravitational effects relative to radial flow is parametrised by the Froude number F r =Û (gĥ 0 ) −1/2 ; g denotes the acceleration due to gravity, and taking F r = O(1) gives Re U δ 0 F r −2 1 and gravitational effects to be neglected. A measure of the significance of the viscous time scale is given by To derive the compressible flow model, the above approximations and assumptions are applied to the compressible Navier-Stokes and momentum equations and conservation of mass equation, resulting in corresponding leading order equations. A first-order Navier formulation slip model is applied, where the dimensionless velocity boundary conditions to leading order are at z = h R and z = h s , respectively. The fluid velocity components tangential to the wall are proportional to the relevant wall shear stress; the proportionality constant l s denotes a dimensionless slip length. For a full derivation see Bailey et al. [9].
Integrating the leading order Navier Stokes equations, together with the conservation of mass equation between the rotor and stator and applying the Leibniz integral rule with the velocity boundary conditions (2), results in the modified nonlinear Reynolds equation for compressible flow [10]: The effects of inertia due to rotation are characterised by speed parameter 2 /(νÛ ) and dimensionless ideal gas constant K s = Rτĥ 2 0 /(µr 0ÛM ) which relates the pressure to the density;M denotes the molar mass, R the gas constant and τ the constant temperature. Dependence on the coning angle, β, is implicit through the rotor-stator clearance h. If λ 1, the centrifugal inertia effects are neglected but the rotor and stator still have relative rotational motion because of the velocity boundary conditions. Similar methodology is followed to derive the incompressible flow model. In this case no-slip velocity boundary conditions are imposed on the bearing faces and the effects of inertia are neglected, however the rotor and stator still have relative ro-tational motion due to the azimuthal boundary condition. The resulting modified Reynolds equation is given by In order to provide a quantitative study of the interaction between the fluid film and the bearing structure of fluid lubricated technology, a fully coupled bearing model is required. Therefore, the above modified Reynolds equations are coupled to a structural model of the bearing faces.

Structural model
The rotor and stator axial displacements are modelled using standard spring-massdamper models which incorporate the pressure variation and the external random time dependent force N (t) experienced by the rotor. The dimensionless stator equation is given by and the rotor equation is for β ≤ 0 and β > 0, respectively. The axial force of the fluid on the rotor and stator is where p a =p a /P is the dimensionless atmospheric pressure. The dimensionless force coupling parameter is given byζ = µÛT 2 /(mδ 0 3 ),m is the mass of a face, D a = D aT /m the damping parameter and K z =K zT 2 /m the effective restoring force parameters.
In a physical application, random fluctuations reflect the uncertainty in the precise state of a system. Two different models are used to describe the external force N (t).

Random force with discrete spectrum
The bearing is envisaged to be part of a larger complex dynamical system, characterised by a periodic motion with periodT ; this is the case of an aero-turbine engine which is composed of several rotating parts. The motion of the individual parts and the interaction with the surrounding environment induces an external force on the rotor N (t) and is taken to be continuous, stationary, having a bounded amplitude, and a finite number of frequencies in its spectrum. This results in a random external force of the form The amplitudes A i are truncated normal random variables, X i are standard normal random variables andĀ > 0 is a constant. The phases φ i are random variables which are uniformly distributed on [0, 2π] and the frequencies B i are random variables taken to be discretely uniformly distributed on the set of non-negative integers: {0, 1, ...,B};B is an upper bound motivated by the need for the forcing time scale to be comparable or longer than the viscous time scale. The random variables X i , φ i , B i , i = 1, ..., n are mutually independent.
A realisation of the external force N (t) in (8) for given A i , B i , φ i , results in a periodic function of time. The force N (t) is a stationary stochastic process [13,14]. As n is finite and the amplitudes A i are bounded random variables, the external force N (t) is also bounded which is often a necessity for modelling a physical force.

Numerical technique
To solve the coupled bearing problem numerically, the modified Reynolds equation (3) is solved simultaneously with the structural equations (5)- (6). Similiar methodology was used in Bailey et al. [10], the main steps are given below for a NCB; the equations for a parallel bearing are achieved when β is set equal to zero.
The modified Reynolds equation (3) is rewritten in terms of the MFC g(t) = h s (t) − h r (t), together with the structural equations (5)- (6). To solve these coupled equations, the Reynolds equation is discretised in the spatial variable r by a second order central finite difference scheme. The problem can then be reduced to a system of first order random ordinary differential equations. For every random force realisation, the rotor experiences a periodic force, N (t) which is of period T = 2π. A map is developed which advances an initial condition at time t 0 by T = 2π, defining a stroboscopic map. Full details of the stroboscopic map solver can be found in [10]. For a given set of initial conditions and pressure boundary conditions, the set of equations can be solved using a stiff ordinary differential equation solver. Results shown follow from using a variable-step, variable-order solver based on the numerical differentiation formulas of orders 1 to 5; due to the nonlinearity of the system of equations, an iterative Newton's method was also used to find solutions numerically.
The stroboscopic map solver was combined with a Monte Carlo technique; a large number of independent random realisations of the stroboscopic map are simulated, with each realisation having a randomly generated external force N (t). The stroboscopic map was initially used to find the approximationḡ min of g min defined as The approximate average of g min was computed as whereĝ min is approximately the average ofḡ min .ḡ where var(ḡ min ) is estimated by and E[ḡ min ] ∈ (ĝ min − M C err ,ĝ min + M C err ), with a probability of 95%. To obtain sufficiently accurate results, an appropriately large number of independent realisations M was chosen to achieve a sufficiently small Monte Carlo error. The estimatê g min has bias E[g min ] − E[ḡ min ], corresponding to the error of the numerical approximationḡ(t) of the stroboscopic map g(t)and can be considered negligible [18]. This methodology can identify the average minimum face clearance between the rotating and stationary faces together with the probability it reaches a given tolerance [18]. The generality of the proposed random force is potentially transferable to other models of physical phenomena beyond this particular application.

Results
The effect of the random force and the probability the MFC clearance reaches a prescribed tolerance was examined. The configuration of the bearing was for that of a narrow annulus a = 0.8, with speed parameter λ = 0.0029 and slip length l s = 0.1. The ideal gas constant was taken to be K s = 1 and the force coupling parameter, ζ = 1.22. The structural parameters are given by K z = 12.2 and D a = 1.5. Ambient pressure was imposed at the inner radius p I = p a = 1 and the pressure at the outer radius was taken to be p O = 0.5 or p O = 1.5 for a NCB or PCB, respectively; the outer pressurisation of a parallel bearing depends on the case being examined. The external force (8) parameters were n = 20,Ā = 20 andB = 10. The parameter values used refer to a bearing with face clearanceĥ 0 = 0.3x10 −4 m, characteristic radiusr 0 = 0.01 m and typical pressureP = 5x10 3 Pa. Structural stiffness is taken asK z = 5x10 5 Nm −1 , damping asD a = 300 Nsm −1 and rotor rotational speed as 6000rpm. Representative values for the bearing face mass is 1Kg and fluid properties areρ 0 = 1.2922Kgm −1 , µ = 2.026x10 −5 Kgm −1 s −1 . Figure 3 gives the mean and variance of g min with corresponding Monte Carlo errors for both a NCB and a PCB; a common case of a parallel bearing is examined twice with the corresponding internal pressure to compare with a NCB and PCB. The number of realisations was chosen such that the Monte Carlo error is ≤ 2.5% giving a compromise between the computational time and accuracy of the results [18]. For the case of a PCB, the average value of g min remains of a similar value as the magnitude of the coning angle increases, whereas a NCB has a decrease in the average value from g min = 0.396 ± 0.0102 for β = 0 to g min = 0.308 ± 0.00799 for β = −0.5. The average value of g min is larger for a PCB than a NCB for all cases considered. The variance is small and of O(10 −2 ) for all cases considered.  Figure 4 shows the corresponding probability of g min becoming as small as a prescribed tolerance g tol for a NCB and PCB. In both cases, decreasing the prescribed tolerance causes the probability of contact to decrease, although results are similar for g tol = 0.001 and g tol = 0; as expected the probability of contact for β = 0 is zero. For a NCB, decreasing the magnitude of the coning angle gives a general decreasing trend of the probability of g min reaching the prescribed tolerance. Whereas a PCB has a general increase in the probability of g min reaching g tol , except for g tol = 0.1 where there is a nonmonotonic trend.
Overall results indicate that the coning angle may have a significant effect of the likelihood of face contact and deformation on the bearing faces due to over pressurisation should be considered in bearing designs and for constraining operating conditions.

Random force with continuous spectrum
In this section, an external random force applied to the bearing is described by a noise term with continuous spectrum. An advantage is that the random force is described more effectively by a limited number of parameters, resulting in a more robust representation of a given disturbance for a specific application which is not restricted by the discretisation.
The case of incompressible flow with a no-slip condition is considered where the modified Reynolds equation (4) describes the fluid flow. A bearing with parallel faces is examined, with structural equations as in (5)-(6) with β = 0. An analytical solution for the pressure field can be computed from solving the Reynolds equation (4), which leads to an analytical expression for the fluid force using (7). Substituting the force expression into the structural equations (5)- (6) and rewriting in terms of the MFC results in a second order ordinary differential equation of the form It is noted that the coefficient 2ζB/g 3 has a singularity when g → 0.
Equation (12) is written as a system of first order differential equations withB = γB = 2ζB andĀ = 1 + ζA/K z and initial conditions of the form g(0) = g 0 , dg/dt(0) = z 0 . A noise term N (t) represents the cumulative effects of multiple weakly coupled environmental degrees of freedom; the bearing may be subjected to disturbances from various sources including vibration, neighbouring components or contaminants. Discussion of the time scales of the physical system and noise are given by Hänggi and Jung [15]. Motivation for choosing coloured noise to model N (t) is that the noise term should have a well-defined characteristic time (correlation time). The noise time scale is dependent on the characteristic response time of the bearing systemT and ideal correlation times are comparable to or larger than the characteristic time of the bearing system to relax. In other words the correlation times should be of O(δ 0 ) or larger. White noise is a limit of coloured noise (see e.g. [15]), when the correlation time tends to zero and it can be a suitable model for external disturbances when the correlation time is very small. This limiting case is simpler to analyse mathematically. The numerical procedure is described in detail for the white noise case and it is similar for the coloured noise. The coloured noise will give a more representative model of the random force due to it having non-zero correlation time.

Numerical technique
In this section the methodology is given for white noise; similar can be implemented for coloured noise. For the case of the external force N (t) represented as a Gaussian white noise, the expected value is E[N (t)] = 0 and covariance is E[N (t)N (s)] = σ 2 δ(t − s);σ is a noise intensity parameter. The set of equations (13) are formulated in the differential form with initial conditions g(s) = g 0 > 0 and z(s) = z 0 . Here w(t) is a standard Wiener process defined on a filtered probability space (Ω, F, F t , P ) and the scaled noise intensity parameterσ = γσ/2. The associated solution notation for (14) is g(t) = g s,g0,z0 (t), z(t) = z s,g0,z0 (t), t ≥ s; g(s) = g 0 , z(s) = z 0 are the initial condition at t = s. If s = 0 the we write g g0,z0 (t), z g0,z0 (t) for compactness. Existence and uniqueness of the solution to (14), as well as boundedness of its moments was shown in [19]; this analysis formally demonstrated that g g0,z0 (t) can never become zero, thus that a positive face clearance is always maintained, although it may become arbitrary small. The time at which the face clearance g(t) first reaches a small given tolerance, δ > 0, is of interest. It is useful to define a space domain G δ = (δ, ∞) × R, leading to time-space domain Q δ = [0, T ) × G δ . The first exit time of the trajectory (g s,g0,z0 (t), z s,g0,z0 (t)), t ≥ s from domain G δ is defined as θ δ , thus θ δ = θ δ s,g0,z0 = inf{t : g s,g0,z0 (t) = δ}, is the first time g s,g0,z0 (t) reaches δ; the time when the face clearance value first becomes as small as the tolerance. The first exit time of the space-time trajectory (t, g s,g0,z0 (t), z s,g0,z0 (t)) from the domain Q δ is defined by τ δ = θ δ ∧ T , which is the time when the face clearance first reaches the prescribed tolerance δ or value of T if the tolerance is not reached over the time examined.
It is noted that at the initial time s = 0, the face clearance g 0 and velocity z 0 should be treated as uncertain. As neither the face clearance or velocity are know at the start of the simulation, because the modelled behaviour follows an initial start-up phase where the dynamics are unknown, the initial conditions g 0 , z 0 are treated as independent random variables; g 0 is chosen to be a truncated normal random variable to ensure its positivity and z 0 to be a normal random variable. These random variables are independent of the Wiener process w(t), t ≥ 0. The expectation of a quantity of interest can be expressed, for some functions φ and f , as [19]: To compute the quantities of interest U , the Monte Carlo technique is implemented together with a simplest random walk for the Dirichlet problem to simulate the SDEs (14) in bounded domains; following methodology in Milstein & Tretyakov [16,17]. The SDEs are simulated on a time interval [0, T ], which is uniformly discretised with time step T > 0; correspondingly the maximum number of time steps in each discrete trajectory in the given algorithms below isN = T / T with t k+1 = t k + T , k = 0, . . . ,N − 1.
Applying the weak Euler approximation with the simplest simulation of noise [17] to the system (14) for given initial conditions (g 0 , z 0 ) ∈ G at t = 0 gives Variables ξ k are mutually independent random variables taking the value of ±1 with probability 1/2. A Markov chain (t k , g k , z k ) is assembled to approximate (t k , g(t k ), z(t k )), where (g(t), z(t)) is the solution of (14), as described by the subsequent algorithm. Algorithm to simulating trajectory (t k , g k , z k ).
Step 3: If g k+1 > δ and k <N − 1, set k = k + 1 and go to step 2; otherwise go to step 4.
The stopping time τ δ is approximated by t κ obtained in the algorithm, which stops when either the approximate trajectory exits the domain G or the integration time reaches the end of time interval T . The approximate discrete trajectories (t k , g k , z k ) and the approximate stopping times t κ are substituted in (15) instead of the corresponding exact values resulting in the approximationŪ instead of U . It can be proved [17,16] that the algorithm is of first weak order, i.e. |Ū − U |≤ C T where C > 0 is a constant independent of T .
The expectationŪ is computed by exploiting the Monte Carlo technique. A large finite number M of independent discrete approximate trajectories are simulated and the corresponding sample meanÛ is found, approximatingŪ , together with the respective sample varianceV . The Monte Carlo errorŪ −Û is quantified by R mc = V / √ M givingŪ belonging to the confidence intervalŪ ∈ Û − 2R mc ,Û + 2R mc with probability 0.95 [17].
Simulated quantitiesÛ incur two errors; the first is the bias, equal toŪ − U , the value of which is determined by choosing a sufficiently small time step T and the second is the statistical error, controlled by selecting a sufficiently large number M of independent realisations.
The same methodology can be applied when the external force N (t) in (13) is modelled as coloured noise. Using the Ornstein-Uhlenbeck process, the following system of SDEs is obtained (see, e.g. [15,17,19]): where α > 0 is constant. As in the white noise case, the bearing gap g g0,z0,y0 (t) always remains positive, though it may become arbitrary small.
Coloured noise y(t) is a Gaussian process and has the properties The Ornstein-Uhlenbeck process y(t) is ergodic with invariant measure y (x) = (πα) −1/2 e −x 2 /α . For the noise term N (t) in (13) to be a stationary process, which is the natural requirement for a model of external disturbances, initial condition y 0 in (17) is chosen as a random variable distributed according to the invariant measure In (17), the characteristic correlation time (time scale) of the noise is given by 1/α. As α → ∞, i.e. the characteristic correlation time 1/α → 0, the correlation function of y(t) tends to the Dirac delta function, with the power spectrum becoming flat. In this limit, coloured noise converges weakly to the white noise [15].
Analogous to the white noise case, we obtain an algorthim for approximating (15) based on simple random walks and the Monte Carlo technique (see further details in [19]).

Results
Results are given for a fixed bearing configuration with a narrow annulus a = 0.8 with ambient pressure imposed at the inner and outer radius p a = p I = p O = 1. The structural parameters are taken to be K z = 55.56 and D a = 1, respectively, and the unscaled noise intensity parameter isσ = 20. The prescribed tolerance of g is taken to be δ = 0.05. The time horizon for examining the bearing is T = 1000; it is assumed that this is the maximum time that a bearing will experience each individual disturbance. If the mean time E[τ δ ] equals T = 1000, it is noted that the minimum gap does not become as small as the prescribed tolerance for the examined conditions with probability one. The time step used is t = 1×10 −5 which is sufficiently small to give a solution of the required accuracy and ensuring the numerical integration error does not exceed the Monte Carlo error. The number of independent trajectories M is chosen such that the Monte Carlo error is ≤ 5.5%. The initial condition g 0 is sampled from a truncated normal distribution, obtained using a normal distribution with mean 1, standard deviation 0.5 together with lower and upper truncation bounds g a = 0.1 and g b = 2. The initial condition z 0 is sampled from normal distribution with mean −1 and standard deviation 0.5.
Results are given for an external white noise term on the rotor and coloured noise with short correlation time α = 100. Of particular interest is the average time for the face clearance to reach the tolerance over a range of values for coupling parameter ζ.
The average time for the minimum gap to become as small as the given tolerance δ = 0.05 against increasing values of force coupling parameter ζ is given in Figure  5. The error bars denote the Monte Carlo error and dashed lines give one standard deviation below and above the mean value. Examining the white noise case shows as the force coupling parameter ζ increases the average time at which the face clearance becomes as small as the prescribed tolerance E[τ δ ] reduces. A similar trend is seen for coloured noise, with α = 100. Investigations show increasing the coloured noise parameter α, results in the average time at which the face clearance becomes as small as the prescribed tolerance E[τ δ ] tending to the values obtained for white noise [19]; results for white noise and coloured noise in Figure 5 are similar. These results can help aid bearing design by giving an indication of average bearing lifetime for types of external disturbances and force coupling parameters.

Conclusions
A model fluid lubricated bearing, comprising a pair of parallel or tapered, coaxial axisymmetric rotating and stationary components separated by a thin fluid film is examined under external loading subject to additional random excitations. A film flow model is derived as a modified Reynolds equation using an axisymmetric lubrication approximation and is coupled to the axial motion of the bearing faces, modelled dynamically as spring-mass-damper systems. Evaluation of the bearing dynamics when experiencing two generic forms of random external forces is presented; an advantage of this is that, dependent on the operational environment of the bearing, one model may be more representative than the other.
The first study considers a random axial force imposed on the bearing, characterised by a fixed period and is taken to be continuous, stationary, and restricted to having bounded amplitude and a finite number of frequencies in its spectrum. The form of this random force was selected as representative of industrially relevant harmonic disturbances and is transferable beyond this particular application. A stroboscopic map solver is used in combination with a Monte Carlo technique to enable evaluation of an average minimum face clearance and notably provide a probability that a minimum clearance tolerance is maintained. Results identify the bearing coning angle may have a significant effect on the probability of contact; parallel bearing faces maintain a positive face clearance, although this bearing gap may become arbitrary small whereas tapered faces may have a non-negligible probability of contact.
A second study evaluates the dynamic effect through a noise disturbance term described by a continuous discrete spectrum represented by a white or coloured noise. A numerical technique is derived and implemented based on a simple random walk combined with a Monte Carlo technique. A detailed approach is provided for the case of a parallel bearing and the external force being taken as a Gaussian white noise source. The average time for the face clearance to reach a prescribed tolerance was examined for both white and coloured noise. Results show the force coupling parameter can have a significant effect on the average time for the face clearance to exceed a prescribed tolerance and provide an indication of the average lifetime of a bearing.
Within industrial applications bearing disturbances may emanate from spurious disturbances associated with the rotational system's linked components. Additionally, accumulated noise from external vibrations may become significant as bearings operate with decreasing tolerances. The methodology and results from the studies presented can inform bearing design and constraints on the operating conditions for safe and reliable behaviour.