Probability of Face Contact for a High-Speed Pressurised Liquid Film Bearing Including a Slip Boundary Condition

An initial deterministic mathematical model for the dynamic motion of a simple pressurised liquid film bearing is derived and utilised to evaluate the possibility of bearing contact for thin film operation. For a very thin film bearing the flow incorporates a Navier slip boundary condition as parametrised by a slip length that in general is subject to significant variability and is difficult to determine with precision. This work considers the formulation of a modified Reynolds equation for the pressurised liquid flow in a highly rotating coned bearing. Coupling of the axial motion of the stator is induced by prescribed axial oscillations of the rotor through the liquid film. The bearing gap is obtained from solving a nonlinear second-order non-autonomous ordinary differential equation, via a mapping solver. Variability in the value of the slip length parameter is addressed by considering it as a random variable with prescribed mean and standard deviation. The method of derived distributions is used to exactly quantify the impact of variability in the slip length with a parametric study investigating the effect of both the deterministic and distribution parameters on the probability of contact. Additionally, as the axial rotor oscillations also have a random aspect due to possible varying excitations of the system, the probability of contact is investigated for both random amplitude of the periodic rotor oscillations and random slip length, resulting in a two parameter random input problem. The probability of contact is examined to obtain exact solutions and evaluate a range of bearing configurations. Lubricants 2015, 3 494


Introduction
Pressurised liquid film bearings are considered in this work consisting of two structural components, a rotor and stator, separated by a thin liquid film which experiences relative rotational motion.The thin liquid film is used to maintain a separation between the rotating and stationary elements under external axial loads on the bearing through the generation of a local film pressure.In practice, these bearings operate as a combined hydrostatic and hydrodynamic bearing where the fluid gap is maintained by the pressurised fluid and the dynamics of the bearings rotational motion.In some practical applications such bearings operate at very high rotational speeds, requiring the inclusion of additional centrifugal inertia effects which are not considered in the classical lubrication formulation of this type of problems (Reynolds equation).Under these conditions, the bearing design requires precise detailed knowledge of variation in the bearing gap during operating condition in order to evaluate possible contact between the rotor and stator.
Analysis of a pressurised fluid film bearing with parallel surfaces is given for a Newtonian flow with a no-slip condition by Garratt et al. [1] where the bearing structure was coupled to the pressurised fluid flow and the dynamics investigated when the lower face undergoes prescribed periodic axial oscillations.A more general geometry with the inclusion of a coned rotor (tapered surfaces) was considered by Bailey et al. [2], allowing more extensive bearing gap investigations.The bearing axial dynamics were examined for prescribed periodic axial oscillations including those of amplitude larger than the equilibrium film thickness.Results indicated that the film thickness can become very small, even reaching film gaps at nano and micro scales.At these scales, surface effects start to dominates over volume related phenomena, requiring accurate details of the flow surface interaction.In a continuous approximation of fluid motions at very small scales, nano and micro, surface slip velocity conditions are typically considered to take into account the additional surface effects.Incorporation of a slip condition on the bearing faces in the analysis of highly rotating parallel bearings with a small face clearance by Bailey et al. [3] showed the bearing faces do not have contact, although the fluid film can become very small.However, for the case of a coned bearing geometry with a slip boundary condition, Bailey et al. [4] showed face contact is possible.
Thin film gas flow regimes over an atomically smooth surface are classified by the Knudsen number, K n = l/ ĥ0 , where l is the fluid mean free path corresponding to the collision distance between molecules (for air at atmospheric conditions l = 68 nm) and ĥ0 is the characteristic fluid thickness.For flow with a Knudsen number between 10 −3 and 10 −1 , denoting micro and nano scales, a continuum flow model with a slip boundary condition is usually employed.Navier [9] was the first to predict the existence of velocity slip, proposing a slip model based on a linear relationship between the tangential shear rate and the fluid-solid tangential velocity difference, and with the slip length as the proportionality constant.
The slip condition for liquid flows are however not as well defined as in gas flows.For liquid flowing over a hydrophobic solid surface, an apparent slip velocity has been observed just above the solid surface, with a slip length of the order of 1 µm and even of the order of 50 µm in the case of super hydrophobic surfaces [11,12].For a hydrophilic surface, a no-slip boundary condition is usually employed in the analysis of liquid flows at micro and nano scale.However experimental studies examining specific chemical and electrochemical conditions of both the solid and liquid showed slip velocities can exist, which are attributed to the effect of the extra reactions on the wettability of the liquid and the potential absorption into the solid surface [13,14].Typical experimental and theoretical studies on the existence of slip velocity at micro and nano scales, are given by [5][6][7][8].
The numerical simulation of a thin film, with thickness of the order of magnitude of the irregular surface roughness, becomes extremely difficult due to the requirement of imposing the no-slip velocity condition on the corresponding irregular boundary.In these cases, it is possible to replace the no-slip boundary condition over the irregular surface by an effective slip boundary condition onto an equivalent smooth surface, with the slip length of the order of the size of the surface roughness (see [15] and [16]).An example is the average surface roughness, R a , of stainless steel that can range from 0.1µm ≤ R a ≤ 1µm, according with surface finishes technology employed.An approach employed by several authors is to find an asymptotically equivalent slip length for small-scale variations of the surface roughness boundary using homogenization theory, (see [17]).
More generally the mathematical formulation and associated numerical simulation of fluid film flow motions, at the nano and micro scale, makes use of slip velocity conditions, where the slip length can range from the size of the surface roughness to the fluid thickness or higher.
This work considers the dynamic analysis of liquid film bearing under extreme conditions, where the fluid gap can attain a thickness of nano or micro scale size.Included are hydrophobic surface bearings, with slip length of the order of the liquid film thickness or larger, as well as hydrophilic surfaces, where the slip length is of the order of the bearings surface roughness.
Precise evaluation of the slip length is restricted due to limitations of the simulations [18] and incorrect interpretation of results [19].Experimental investigation of the slip length requires length scales comparable to the slip length to be examined, therefore extremely accurate techniques with high spatial resolution, capable of interfacial flow measurements are needed [20].Indirect experimental methods measure macroscopic flow quantities [21] but can suffer from limited spacial resolution, with the slip velocity inferred from using techniques which can rely on tentative theoretical models [22].Direct methods also suffer from limited spatial resolution [23], with the fluid flow velocity, which may measure only a few nanometres, are typically determined in an interrogation area which is several hundred nanometres in size [24] leading to potential errors in measurements.
Increased interest exists on accounting for greater predictability of the dynamic behaviour in industrial bearing designs.These are typically complex, with associated mathematical models involving large numbers of parameters that may only be known approximately.Sources of uncertainty in a bearing model include the slip flow at the fluid-solid interface and from the random axial motion of the rotor due to external excitations.
Uncertainty in a bearing model and dynamical system originating from external excitations has been examined using a variety of methods, typically based on Monte Carlo techniques [25].In this case, sampling is directly from the random input parameter with the output of the deterministic model evaluated at each input sample.This is a robust process able to deal with complex situations and a large number of uncertain input parameters, however an extensive number of code executions may be required to produce sufficiency accurate results.To minimise cost and efficiency additional theoretical methods for uncertainty quantification have been developed, such as polynomial chaos, Gaussian process emulator and Markov Chain Monte Carlo method.
In Zigang et al. [26], the dynamic response and vibration characteristics of two in-line rotor bearing systems are examined to quantify the effects of uncertainty.The equations of motion of the rotor system are derived and stochastic modelling for uncertainty in the damping and nonlinear support stiffness is developed, based on the polynomial chaos expansion technique in a stochastic framework.Monte-Carlo simulation is used as a comparison, with good agreement reported between the two methods.Response statistics and the dynamic behaviour of the stochastic system are examined through the mean and probability density function.Results show the uncertainties in the nonlinear support stiffness and damping have a significant effect on the predicted behaviour of the rotor system, and increasing the random intensities causes the realisation amplitudes to span a wider band of amplitudes.
In the special case where a comprehensive deterministic mathematical model is available, a subjective interpretation of probability, the method of derived distributions [27], can be used to exactly quantify the effects of random input parameters.Describing the random input parameters by a probability distribution function and using the deterministic relationship between the input and output parameters enables the probability density function and cumulative density function of the output to be determined exactly.Statistical analysis can be conducted on the resulting realisations allowing the relevant statistical properties of the process to be extracted.
A main focus of this work is to analyse the probability of contact in an important application requiring a mathematical and numerical model for a high speed pressurised liquid film coned bearing with a slip boundary condition and prescribed periodic axial oscillations of the rotor; the deterministic model is based on the work by Bailey et al. [4].The governing Reynolds equation for slip flow is derived by incorporating a constant slip model with a linear relationship between the tangential shear rate and the fluid-wall velocity differences with a slip length as the proportionality constant.The probability of bearing face contact is examined using the method of derived distributions to investigate the lack of precise knowledge of the slip length.To account for variability in the slip length and amplitude of rotor oscillations, these parameters are considered as random variables with prescribed mean and standard deviation to investigate the probability of face contact through a parametric study.
Although in liquid film bearings design, the use of a tapered surface is considered to increase the liquid loading force in order to maintain a desirable film thickness at operating conditions, in this work it has been found that this can lead to undesirable aspects under extreme conditions.If during the dynamic response of the bearing, under an external force, the liquid gap reaches a thickness of nano or micro scale, the tapered surface can induce possible stator-rotor collision (contact) dependant on the liquid slip velocity condition.
An approach is provided for determining the probability of contact at different operating conditions to inform bearing design criteria.It is noted that for controlling the possible displacements of the rotor due to external forces, using structural constraints, the dynamic response of the bearing may encounter very small fluid gaps and the use of a tapered surface needs to be carefully evaluated.This may be only used where necessary and the corresponding angle needs to be limited to a small value.

Deterministic Mathematical Model
A simplified mathematical model of a high speed pressurised fluid film bearing for a coned rotor and incorporating a slip boundary condition is developed, following Bailey et al. [4].A model for the incompressible fluid flow through the bearing is derived from the Navier-Stokes momentum and continuity equations in axisymmetrical coordinates.To express the mathematical formulation of the problem in dimensionless variables, a characteristic bearing radius r0 , film thickness ĥ0 , rotor velocity Ωr and differential pressure ∆ P, due to imposed pressure at inner radius rI and outer radius rO are defined with dimensionless time variable t = ω t where ω is the angular frequency.Dimensionless velocities are taken to be û/ ĥ2 0 ∆ P(µ (r O − rI )) −1 , v/ Ωr 0 and ŵ/ ĥ0 T −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 = ls / ĥ0 and for a thin film bearing the aspect ratio δ 0 = ĥ0 /r 0 is small δ 0 1.The rotor has a fixed coning angle, assumed to be small giving cos β 1 and sin β β , leading to the scaling β = β δ 0 with β = O(1), giving consistency with the lubrication condition, such that Cases of a positively and negatively coned bearing are considered separately as is shown in Figure 1, referred to as a PCB and NCB, respectively.For a bearing containing a coned rotor undergoing prescribed periodic oscillations, the dimensionless film thickness is given by where h s (t) is the height of the axially displaced stator, ε sint the axial forcing of the rotor with amplitude ε and a = rI /r 0 is a measure of the bearing width, 0 < a < 1.The associated radial Reynolds number is given by Re U = ρ ĥ0 Û/ μ, where ρ is the density and μ the dynamic viscosity.To ensure that the effects of viscosity are retained at leading order the pressure is scaled as P = μ r0 Û/ ĥ2 0 .Classical lubrication theory neglects inertia due to the reduced Reynolds number Re U δ 0 1, however, as in Garratt et al. [1], the centrifugal inertia is retained to include cases of high rotational speed bearing operations for which terms of the order Re U δ 0 (Re * ) 2 are considered to be of O(1), with (Re * ) 2 = ( Ωr 0 / Û) 2  1.The squeeze number σ = r0 /( Û T ) characterises any time dependent effects whilst the Froude number Fr = Û/( ĝĥ 0 ) 1/2 , where ĝ is the acceleration due to gravity, parametrises the importance of the gravitational effects relative to the radial flow.However, gravity can be neglected given Re U δ 0 Fr −2 1; this is consistent with lubrication theory provided the Froude number is O(1).
In the case of slip flow, as considered by Bailey et al. [4], the velocity boundary condition on a bearing face comprise of a jump in the tangential velocities across the fluid-solid interface corresponding to the slip velocity induced due to wall shear, i.e., a Navier slip condition.Continuity of the normal velocity across the fluid-solid interface is imposed along with a no flux condition.Thus the velocity boundary conditions on the rotor and stator, denoted by superscript r and s, respectively, are given by û where the right hand side of the above equations are due to the surface shear stress with proportionality constant ls and rate of strain tensor êi j .The rotor velocity components are given by ûr = (∂ ĥr /∂ t, 0, Ωr) corresponding to prescribed axial motion and fixed azimuthal rotation and the stator is allowed to move axially due to the interaction with the fluid, giving the stator velocity as ûs = (∂ ĥs /∂ t, 0, 0).Due to the additional slip component on the bearing face both the normal and tangential velocities on the rotor and stator are required.A PCB and NCB are considered separately, denoted by subscripts + and −, respectively, in the coordinate system (ẑ, r, θ ) with velocities û = ( ŵ, û, v). Figure 1 shows the normal rotor vector nr and two orthogonal tangential vectors t1 and t2 on the rotor face in the cases of a PCB and NCB.The normal and azimuthal tangent are given by nr = (cos β , sin β , 0), t2 = (0, 0, 1) respectively, and the axial tangents are Cases of a PCB and NCB, respectively, requiring separate radial velocity boundary conditions to be derived.
The stator lies parallel to the radial direction, giving the normal vector ns in the negative ẑ direction and two orthogonal tangential components, t3 and t4 on the face of the stator in the radial and azimuthal direction, respectively.
Velocity boundary conditions on the rotor in (2) are given by using the normal and tangent vectors given in ( 3) and ( 4) and the rate of strain tensor components as given in Batchelor [28] (p.602).In an axisymmetric configuration the azimuthal derivatives do not appear.
The velocity boundary conditions on the stator (2), give Under these conditions, the leading order non-dimensional rotor and stator velocity boundary conditions given in ( 5) and ( 6), respectively, become as δ 0 1.The only dependence on the coning angle is given in the axial rotor velocity boundary condition.
A governing equation for a bearing with narrow gap, satisfying the slip boundary conditions, is readily obtained from the leading order thin film approximation of the continuity equation, where formally terms of O(δ 0 ) are neglected, giving the modified Reynolds equation as The speed parameter λ = 3Re U δ 0 (Re * ) 2 /10 characterises the importance of centrifugal inertia and squeeze number σ = 12 σ is rescaled.The Reynolds Equation ( 8) has no explicit dependence on the coning angle and expresses the relationship between the pressure p and film thickness h.Similar Reynolds equations have been derived for slip flow in a bearing [29][30][31], but neglecting the effects of centrifugal inertia and in [3,4] including inertia effects.Pressure boundary conditions at the inner and outer radii of the bearing are defined in Appendix A, Equation (37).
Axial displacement of the stator is modelled using a spring-mass-damper model, where the stator equilibrium position at steady flow conditions is given by ĥ0 with m ĝ as the stator weight and Kz as the restoring force coefficient (stiffness), respectively.The corresponding dynamic equation of the stator axial position in dimensionless variables is given by where D a = Da / m ω and K z = Kz / m ω2 are dimensionless linear damping and restoring force coefficients, respectively.In (10) Da is the damping coefficient and F(t) is the resultant dimensionless axial force on the stator, defined in Appendix A by Equation (45).The force coupling dimensionless parameter is given by α = μ Û/ m ω2 δ 3 0 and considered to be of O(1).The dynamic Equation (10) for the axial position of the stator, provides the corresponding dimensionless stator equilibrium bearing height, corresponding to the geometric configuration in Figure 1 with film thickness h = 1 at the inner radius for PCB and at the outer radius for NCB; giving h s = 1 for PCB and h s = 1 + (1 − a)β for a NCB, according to Equation (9).
To solve the Reynolds Equation ( 8) and stator Equation ( 10) simultaneously, it is mathematically convenient to introduce the time dependent minimum face clearance (MFC) as Integration of the Reynolds equation gives explicit analytic expressions for the pressure field and force on the stator, given in Appendix A, Equations ( 38) and (46), respectively.Substituting the force expression (46) and MFC (11) into the stator Equation ( 10) yields a nonlinear second-order non-autonomous ordinary differential equation for the MFC given by where for a PCB.Expression for A and B are given in Appendix A, Equations ( 47) and (48) respectively.A NCB only differs in the term; S(g, λ , β , l s ) = K z (g − (1 − a)β − 1) − απA(g, λ , β , l s ).Dynamically Equation ( 12) 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 K zT is found by differentiating the effective restoring force with respect to the MFC, giving where K z is the contribution from the structure and the fluid contribution is given by using the definition for A in (47).
The functions G and G I /g are given in Appendix A, Equations ( 42) and (50), respectively, with derivatives ∂ G/∂ g and ∂ G I /∂ g taken as the derivative with respect to the MFC.
In the case of a non-pressurised parallel bearing the fluid stiffness is given by It is noted that in the case of a non-pressurised parallel bearing using classical lubrication theory, i.e., neglecting inertial effect and slip condition, predictions give zero contribution of the stiffness of the system from the fluid component.

Numerical Scheme
It is anticipated that for a fixed value of the slip length, l s , and rotor oscillation of period T , the nonlinear non-autonomous dynamic system (12) admits periodic solutions.Numerical schemes that are able to predict periodic solutions of nonlinear differential equations are not always reliable and require special considerations.In this work a numerical technique is implemented using a stroboscopic map to obtain the solution in a prescribed period T .
Periodic solutions to Equation ( 12) are denoted by the vector g(g(t), f (t)) with f = dg/dt and corresponding values of (g(t), f (t)) at the beginning of the period given by g(t 0 ) = g 0 , and f (t 0 ) = f0 , which are denoted here as the initial conditions of the periodic motion and considered as unknown in the numerical algorithm for increasing slip lengths.Periodicity is obtained when the values of the solution at the end of the period are identical to their initial conditions, i.e., g(T ) = g 0 , and f (T ) = f0 .
In the numerical algorithm, the second order differential Equation ( 12) is expressed as the following system of first order coupled differential equations: The numerical scheme consists of developing a mapping function that advances any initial condition g 0 at time t 0 by a time T , defining a stroboscopic map φ (T ; g 0 ,t 0 ).This ℜ 2 → ℜ 2 map integrates the system of Equations ( 17) forward through one period of the rotor oscillations.To identify periodic solutions, the fixed points of the stroboscopic map g(t) = g(t + T ) are found iteratively which corresponds to the condition giving periodic solutions g(g(t), f (t)).Solutions are found numerically from an iterative Newton's method, given an initial guess value g0 .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 period, t = T .
For any given initial condition, a solution of the system of Equations ( 17) and ( 21) for t 0 ≤ t ≤ T can be found using Matlab ode15s.This procedure is repeated recursively, with the improved value of g 0 used in the system of Equation ( 17).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.
Periodic solutions for increasing slip lengths l s + l s , are determined from successive new initial conditions.An Euler scheme (parameter continuation), includes the slip length as a variable 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., An Euler predictor step is then performed as to give an updated initial condition g 0 for l s + l s .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 This system of equations can be coupled with the previous augmented system of equations and solved using the same Matlab routine.
For this new value of the initial condition, the above Newton's method is repeated, with successive refinement of l s if necessary, until convergence is achieved and a periodic solution for l s + l s is found.With the use of this approach a consistent new initial condition for an incremented value of the parameter is defined.
A major advantage of the formulated numerical (Euler) procedure in the stroboscopic map solver is that it can be directly extended to find threshold values of the slip length l s corresponding to a specified value of 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 requiring the value of the unknown vector g = (g(T ), f (T ), l s ) to be determined given an initial guess value g 0 = (g 0 , f0 , l s0 ).Correspondingly an additional constraint equation g min − g * = 0 is added, with g * as the prescribed value of g min and the corresponding Jacobian matrix is given by The additional terms in the Jacobian matrix (27) with respect to (20), i.e., the last row and column in (27), are obtained from the solution of the augmented system of first-order differential Equations ( 21)- (25).The values in the last row are determined at the time when g min is achieved.The threshold slip length value is calculated at g * , which can be stepped down successively to the point of contact, g * = 0. Continuation is used for a new value of g * to get with the Jacobian matrix as defined in ( 27) and g * < 0. The value of ∂ g/∂ g * is evaluated numerically by a first-order 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 specified value of g min , g * + g * to be found.

Deterministic Results
Deterministic results obtained from the stroboscopic map solver are examined to investigate the value of the slip length parameter at decreasing values of a prescribed minimum value of the MFC over the period (g min ).A PCB and NCB are examined under external (p O = 2, p I = 1) and internal (p O = 1, p I = 2) pressurisation, respectively, as these configurations generally give the better operating conditions.Figure 2 shows the MFC over a period of a rotor oscillation in the case of a PCB with non-zero slip length.The MFC follows the path of a distorted negative sine curve, which is effectively constant and of very small magntiude between t = 1.32 and t = 2.05, with g min = 4.97 × 10 −4 .Physical dimensions of the the bearing corresponding to the result in Figure 2   Table 1 shows the value of the slip length parameter when g min reaches a prescribed tolerance of decreasing magnitude 1 × 10 −3 ≥ g tol ≥ 0 in the case of a PCB and NCB with increasing amplitude of rotor oscillations, the smallest of which is close to where contact first occurs.Increasing the amplitude of rotor oscillations gives lower values of slip length for a given prescribed tolerance.The slip length values increase for decreasing prescribed tolerances, with tolerances g tol = 1 × 10 −5 and g tol = 0 having effectively the same values; the former is significantly more straightforward to compute requiring typically only one run of the code opposed to multiple runs of the code with decreasing step size in g min for the later.A numerical approximation of face contact is adopted, when g min becomes less than a prescribed tolerance taken in this case as g tol = 1 × 10 −5 .
Table 1.Values of slip length for a given tolerance 0.001 ≥ g tol ≥ 0 in the case of a PCB β = 0.2 and NCB β = −0.2 for increasing values of amplitude of rotor oscillations 1.25 ≤ ε ≤ 1.60 and 1.12 ≤ ε ≤ 1.60, respectively; a = 0.2, σ = 1, λ = 1, α = 1, K z = 10 and D a = 1.To compute the probability of contact, the values of the slip length and amplitude of rotor oscillations where contact first occurs are evaluated as shown in Figure 3. Configurations of a PCB and NCB, in the case of decreasing annulus and angle, are examined where face contact initially occurs for values of the slip length and amplitude of rotor oscillations on the curve in Figure 3.For parameter choices of slip length and amplitude of rotor oscillations above the curve, face contact occurs but a fluid film is maintained between the bearing faces for parameter choices under the curve.In the cases considered, for decreasing bearing width 0.8 ≥ (1 − a) ≥ 0.2 face contact first occurs at lower values of amplitude of rotor oscillations for a given slip length, with the exception of a NCB with sufficiently large angle and slip length.Decreasing the coning angle of the rotor gives contact first occurring at larger values of slip length and amplitude of rotor oscillations.The contact curves exhibit non-smooth behaviour which is more distinct in a NCB than a PCB.The existence of this feature has been verified by using values of the slip length and amplitude of rotor oscillations close to the contact curve in the stroboscopic map solver, resulting in solutions of g min tending to the contact tolerance g tol = 1 × 10 −5 , and thus to the current contact curve.Analogous contact plots are produced for slip length against amplitude of rotor oscillations for a PCB and NCB with decreasing angle and increasing speed parameter in the case of a wide and narrow annulus.These results will be used to calculated the probability of contact during the parameter study.From a bearing design aspect, if the fluid gap becomes very small the tapered surfaces need to be limited to a small angle.

Uncertainty in Slip Length
Initially the slip length parameter is considered as a random variable reflecting the typical lack of knowledge of its precise value and the resulting uncertainty on the magnitude of g min is determined.Taking a constant slip length along the surface of the rotor and stator whose value is not known with certainty, a subjective interpretation of probability is used to model this lack of knowledge.Using the method of derived distributions, the corresponding effect of variability in the slip length parameter can be exactly determined, when it is expressed as a random variable.To incorporate that the slip length is always positive, a truncated log-normal probability density function (pdf) distribution is used providing a non-negative value with possible value ranging from the magnitude of bearing surface roughness to the order of the film thickness, or larger according to the nature of the fluid and solid material.
The truncated log-normal pdf and cumulative distribution function (CDF) used for the slip length are given by in terms of the untruncated pdf, fL s and CDF, FL s , with mean µ l s and standard deviation σ l s .The untruncated log-normal CDF is defined in terms of the untruncated pdf by For example when a carbon surface is used, which is known to have large value of slip length, the log-normal distribution is right truncation at l s = 3. Initially the average value and standard deviation will be taken as 0.1, with the effect of alternative average values and standard deviations examined later.
Truncating the distribution however may lead to the specified mean not being representative of the average value of the distribution.Figure 4 shows the pdf and CDF with standard deviation σ l s = 0.1 and effective median M at the desired average value, M l s = 0.1, requiring the prescribed mean to be larger than the average value; µ l s = 0.127.The discrepancy between the effective median and prescribed mean increases as the standard deviation increases.Therefore throughout this work the prescribed mean of the distribution is such that the effective median is representative of the desired average value of the slip length.Describing the slip length by a pdf allows the corresponding pdf of g min to be found using the change of variables, f Gmin (g min ) =| ∂ l s /∂ g min | f L s (l s ), where the derivative is found from the stroboscopic map solver.The CDF of g min can also be calculated, giving the probability g min is less than or equal to a given tolerance, thus the probability of contact is given when g min ≤ 1x10 −5 .To calculate the CDF of g min , the set of deterministic slip lengths giving g min less than or equal to each possible specified value is required; namely A l s g min = {L s : χ ε (l s ) ≤ g min }.In this set, χ ε denotes the deterministic relation between the slip length and g min which can be found for different amplitude of rotor oscillations ε using the mathematical model.The untruncated CDF of the minimum face clearance is calculated by Thus the truncated CDF for g min is given by using the untruncated functions given in (30).
Comparison of the CDF of g min obtained from the Monte Carlo method and by using the method of derived distributions is given in Figure 5.In this case the Monte Carlo method is run with an increasing number of steps and the slip length values are taken at random from the truncated distribution.Obtaining the values of g min for each chosen slip length via the stroboscopic map solver allows the CDF and 95% upper and lower confidence bounds to be computed.The CDF computed by the method of derived distributions uses Equations (32) and the deterministic relation between the slip length and g min to find the set of values A l s g min .Increasing the number of steps used in the Monte Carlo method gives the CDF converging to the result computed via the method of derived distributions.The Monte Carlo method was computationally more expensive, requiring 10, 000 steps to achieve the same level of accuracy as the method of derived distributions.For further calculations the method of derived distributions was used to produce subsequent results.The bearing configuration in this case gives g min ≥ 0.210, giving a zero probability of contact.Figure 6 shows the deterministic relationship between the slip length and g min and the corresponding pdf and CDF of g min for increasing amplitude of rotor oscillations in the case of a PCB and slip length distribution; median M l s = 0.1, standard deviation σ l s = 0.1.The deterministic results show the amplitude of rotor oscillations taking ε = 1.00 and 1.07 with increasing slip length have g min decreasing asymptotically to a constant value.However for larger amplitude of rotor oscillations ε = 1.10, g min decreases until face contact at l s = 0.186.The pdf of g min shows smaller values are likely to occur as the amplitude of rotor oscillations increase, with the CDF showing the smallest value of g min with probably 1 decreases as the amplitude of rotor oscillations increase.The probability of contact is given when the CDF intercepts the y-axis, giving P(contact) = 0 for ε = 1.00 and 1.07 but P(contact) = 0.185 for ε = 1.10.
For parameter studies of the probability of contact, it is convenient to compute the probability of contact directly using where A l s g tol = {L s : χ ε (l s ) ≤ g tol } is the set of discrete values of the slip length parameter which give face contact for different deterministic values of the amplitude of rotor oscillations.Therefore the probability of contact is given as a function of the amplitude of rotor oscillations.The set A l s g tol can be found using the deterministic contact curves in Figure 3 and analogous plots for increasing speed parameter in the case of wide and narrow annuli.The probability of contact against amplitude of rotor oscillations for increasing speed parameter in the case of a PCB and NCB is shown in Figure 7.In the case of a PCB increasing the amplitude of rotor oscillations increases the probability of contact, with the probability of contact increasing from 0 to 1 over a small range of amplitude of rotor oscillations.For a given value of amplitude of rotor oscillations the probability of contact increases for an increasing speed parameter due to a PCB being externally pressurised and the effects of inertia overcoming the differential pressure to decrease g min and therefore increase the probability of contact.On the other hand a NCB has P(contact) remaining less than 1, see Figure 7b, although it becomes large at the upper limit of amplitude of rotor oscillations ε = 2. Increasing the speed parameter reduces the probability of contact for a given amplitude of rotor oscillations as a NCB is internally pressurised giving the effects of inertia coinciding with the differential pressure to increase g min and thus decrease the probability of contact.
The probability of contact against amplitude of rotor oscillations for a PCB and NCB in the case of decreasing angle is given in Figure 8, where increasing the amplitude of rotor oscillations increases the probability of contact.Decreasing the coning angle in a PCB and NCB decreases the probability of contact for a given amplitude of rotor oscillations, tending to the limit of a parallel bearing which has zero probability of contact under the conditions considered in this work.It is important to observe that the non-contact condition found for a parallel bearing is due to the inclusion of inertial effect not considered in classical lubrication theory.Under non-pressurisation conditions and zero inertial effect the fluid stiffness ( 16) is zero and consequently contact is expected under extreme conditions.A PCB has zero probability of contact for larger amplitude of rotor oscillations than a NCB, but larger probabilities of contact for amplitude of rotor oscillations larger than a given value which is less than ε = 1.50 for all the cases considered.The non-smooth behaviour in these plots is due to the discontinuity in the deterministic solutions as shown in Figure 3.
The effect of decreasing bearing width on the probability of contact is shown in Figure 9 over increasing amplitude of rotor oscillations for a PCB and NCB, where increasing the amplitude of rotor oscillations increases the probability of contact.Deceasing the bearing width gives an increased probability of contact for a given amplitude of rotor oscillations, with a PCB having larger increases than a NCB.In the case of a narrow annulus a = 0.8, a NCB has a larger probability of contact than a PCB for a given amplitude of rotor oscillations, whereas other bearing widths considered have a PCB or NCB having a larger probability of contact depending upon the magnitude of amplitude of rotor oscillations.
The effect of the probabilistic parameters on the probability of contact is examined by considering the case of a fixed median with increasing standard deviation and increasing median with fixed standard deviation as shown in Figure 10 for a PCB.Increasing the standard deviation of the slip length causes the probability of contact to increase for a given amplitude of rotor oscillations in the case of ε < ε c , but to decrease otherwise, with critical value ε c = 1.52 where the probability of contact curves intersect in Figure 10a.This situation arises as increasing the standard deviation gives slip lengths over a greater range, with ε < ε c leading to an increase in the probability of contact.Whereas for ε > ε c the range of slip lengths which do not give contact increases and thus decreases the probability of contact.The effect of slip lengths with increased median is shown in Figure 10b where a larger probability of contact for a given amplitude of rotor oscillation, with the transition of probability of contact increasing over a shorter range of amplitude of rotor oscillations.Corresponding plots for a NCB give a similar trend to a PCB with a small magnitude change.For sufficiently large amplitude of rotor oscillations a NCB has larger probability of contact than a PCB, otherwise the reverse situation is true.

Uncertainty in Slip Length and Amplitude of Rotor Oscillations
An additional source of uncertainty in bearing operation is from the axial oscillations of the rotor which are subject to a random excitation.A simplified case where oscillations are considered periodic with random amplitude is given by Hou et al. [32].For the current study the oscillations are considered periodic with constant amplitude, but the value is not known with certainty, allowing an exact solution for a two parameter random input problem to be obtained.Expressing the slip length and amplitude of rotor oscillations as random variables, and using the deterministic results for the relation between the values of slip length and amplitude of rotor oscillations at contact (for example Figure 3), allows the probability of contact to be calculated.
The rotor oscillations are described by h r = ε sint, where the amplitude ε is defined as a random variable instead of a fixed parameter.Following the approach adapted in the previous section the slip length is described by a truncated log-normal pdf.The amplitude of the rotor oscillations is modelled by a double truncated normal function, to include the expectation of little skewing.Truncation of the normal distribution is required to ensure the amplitude of rotor oscillations is non-negative and, as in realistic applications, the amplitude of rotor oscillations is restricted; in the results shown the upper limit is set at ε = 2 which is twice the equilibrium gap g = 1.Initially the average value of the amplitude of rotor oscillations is taken as ε = 1.2 to reflect extreme conditions, with the amplitude of rotor oscillations larger then the equilibrium gap.The standard deviation is taken as σ ε = 0.1, with the effects of changes in the average value and standard deviation examined later.
The truncated normal pdf and CDF for the amplitude of the rotor oscillations are given by with mean µ ε and standard deviation σ ε .The untruncated pdf, fE and CDF, FE are defined by It is assumed the amplitude of rotor oscillations and slip length random variables are statistically independent, giving the joint pdf as the product of their marginal pdf f E,L s (ε, l s ) = f E (ε) f L s (l s ).To compute the probability of contact the deterministic contact curves are used, as given in Figure 3, to find the set of amplitude of rotor oscillations A εg tol = {E : G min ≤ g tol } and slip lengths A l s g tol = {L s : G min ≤ g tol } which give g min less than or equal to the contact tolerance g tol = 1 × 10 −5 .The probability of contact is computed numerically using whose value is represented by the area above the contact curve in the probability space.The probability of contact is computed numerically using Matlab.
The probability of contact for differing bearing geometries in a PCB and NCB are shown in Figure 11 with decreasing coning angle and bearing width.Decreasing the bearing width (1 − a) increases the probability of contact, with greater increases as the bearing width decreases.Decreasing the coning angle in a PCB and NCB reduces the probability of contact for all given bearing widths.A PCB and NCB with small coning angle | β |= 0.1 have similar probabilities of contact, whereas for increasing coning angles a PCB has consistently lower probability of contact than a NCB.
Table 2 shows the probability of contact for a PCB and NCB in the case of a wide (a = 0.2) and narrow (a = 0.8) annulus with increasing speed parameter.Results show increasing the coning angle consistently increases the probability of contact and the probability of contact is greater in a bearing with a narrow annulus compared to a wide annulus.Increasing the speed parameter has an opposing effect between the two bearing types, a PCB has the probability of contact increasing whereas a NCB has a decreased probability of contact.As a PCB is externally pressurised, the effects of inertia overcome the differential pressure giving contact at a smaller value of l s and therefore a larger probability of contact as the speed parameter increases.Whereas a NCB is internally pressurised giving the effects of inertia coinciding with the differential pressure to give contact at a larger value of l s as the speed parameter increases and thus a smaller probability of contact.A NCB has a larger probability of contact than a PCB, with the exception of a sufficiently large speed parameter as indicated by the shaded region in Table 2.In these cases, the contact curves for a PCB sit below those for a NCB, with contact at smaller values of slip lengths and amplitude of rotor oscillations.
Table 2. Probability of contact for a NCB (β < 0) and PCB (β > 0) in the case of a wide a = 0.2 and narrow a = 0.8 annulus with increasing speed parameter 0 To evaluate the effect of values associated with the distribution parameters, the probability of contact is examined for increasing standard deviations with fixed medians and increasing medians with fixed standard deviations.Figure 12 shows the probability of contact against the standard deviation of the amplitude of rotor oscillations with increasing standard deviation of the slip length and against median of amplitude of rotor oscillations with increasing median of the slip length in the case of a PCB.Increasing the standard deviation of the amplitude of rotor oscillations increases the probability of contact.The effect of increasing the standard deviation of the slip length however depends on the size of the standard deviation of the amplitude of rotor oscillations with respect to a critical value σ εc , given as the value where the probability of contact curves intersect in Figure 12a.Increasing the standard deviation of the slip length gives the probability of contact increasing (decreasing) for standard deviations of the amplitude of rotor oscillations below (above) the critical value σ εc = 0.424, due to more values of the slip length giving contact (no contact) being included.Results in Figure 12b show increasing the median value of the amplitude of rotor oscillations gives negligible probability of contact until ε ≤ 1. Increasing the median of the slip length increases the probability of contact.Corresponding plots for a NCB show the probability of contact for increasing standard deviation of the slip length and amplitude of rotor oscillations has a similar trend but with a decrease in magnitude for a given choice of probabilistic parameters.In the case of increasing median of the slip length and amplitude of rotor oscillations a NCB has a similar trend to a PCB but with a small magnitude shift.For the cases considered in Figure 12; small slip length medians, M l s = 0.05 and 0.1 give a larger probability of contact in a PCB than a NCB, effectively the same for median M l s = 0.5 and for larger slip length median M l s = 1.0 a NCB has larger probability of contact than a PCB.

Conclusions and Summary
A deterministic mathematical model is extended for a high speed pressurised thin film bearing containing a rigid coned rotor with a slip length boundary condition imposed on the bearing faces.The axial motion of the stator is modelled as a spring-mass-damper system to investigate a fully coupled unsteady bearing when subject to uncertainty in amplitude of rotor oscillations and slip length parameter.In an initial deterministic analysis of the film the Reynolds equation and stator equation are given in terms of the time dependent minimum face clearance, allowing explicit analytic expressions for the pressure field, force on the stator and a nonlinear second-order non-autonomous ordinary differential modified stator equation to be derived.A stroboscopic map solver is implemented to find periodic solutions.The performance of two coned bearing configurations, a PCB and NCB are considered under external and internal pressurisation, respectively and examined under extreme operating conditions.
Recognising the realistic implications of uncertainty in the model parameter values the slip length parameter is taken as a random variable with prescribed mean and standard deviation together with similar attributes for the amplitude of the axial rotor oscillations.The value of the prescribed mean is such that the effective median of the random parameter is representative of the required average value.The deterministic mathematical model is used to find the relation between the slip length and g min , allowing the pdf of g min to be computed.The CDF of g min is found via the method of derived distributions from which the probability of contact can be found; a numerical approximation of face contact is used to significantly reduce the computational requirement without loss of accuracy.The probability of contact is determined directly for the parameter studies.
A detailed study is provided to evaluate the effects of uncertainty in values of the slip length parameter on potential bearing geometries associated with coning angle, bearing width and rotational speed.Results show an increase in the probability of contact as the amplitude of rotor oscillations increase, a wide bearing has a lower probability of contact than a narrow bearing and a small angle giving a reduced probability of contact compared to a large angle.Increasing the speed parameter causes a decreased probability of contact for a NCB but increased probability in a PCB due to the inertial effects overcoming the differential pressure effects.Examining the effects of different distribution parameters shows increasing the median of the slip length increases the probability of contact and increasing the standard deviation for an amplitude of rotor oscillation below a critical value ε c increases the probability of contact whereas for amplitude of rotor oscillation ε > ε c decreases the probability of contact.
To reflect operational bearing conditions the slip length and amplitude of rotor oscillations parameters are considered as random variables with a prescribed mean and standard deviation and the method of derived distributions is used to calculate the probability of contact.Results show decreasing the bearing width and increasing the coning angle increases the probability of contact.Increasing the speed parameter gives the probability of contact increasing in a PCB but decreasing in a NCB.Increasing the slip length standard deviation causes the probability of contact to increase for standard deviations of the amplitude of rotor oscillations below a critical value σ ε < σ εc , but decrease for σ ε > σ εc .Increasing the median of the slip length and amplitude of rotor oscillations increases the probability of contact.respectively.Similar expressions are computed for a NCB.Equations for a no-slip bearing are given in [2].The force on the stator is defined as

1 Figure 4 .
Figure 4. Probability density function (pdf) and cumulative distribution function (CDF) of slip length with standard deviation σ l s = 0.1 and M l s = 0.1, with prescribed mean µ l s = 0.127.

Figure 6 .
Figure 6.Probability density function (pdf) of slip length, deterministic relationship between g min and slip length, pdf of g min and CDF of g min for increasing amplitude of rotor oscillations (a) ε = 1.00,(b) ε = 1.07 and (c) ε = 1.10 in the case of a PCB β = 0.2 with narrow annulus and slip length distribution with median M l s = 0.1, standard deviation σ l s = 0.1; a = 0.8,

Figure 7 .
Figure 7. Probability of contact against amplitude of rotor oscillations for increasing speed parameter 0 ≤ λ ≤ 5 in the case of (a) PCB; β = 0.2 and (b) NCB β = −0.2 with wide annulus and slip length distribution median M l s = 0.1, standard deviation σ l s = 0.1; a = 0.2,

Figure 8 .
Figure 8. Probability of contact for (a) PCB and (b) NCB with decreasing coning angle with wide annulus and slip length distribution median M l s = 0.1, standard deviation σ l s = 0.1,

σFigure 12 .
Figure 12.Probability of contact against (a) standard deviation of the amplitude of rotor oscillations with increasing standard deviation of the slip length 0 ≤ σ l s ≤ 1; M ε = 1.2, M l s = 0.1 and (b) median of the amplitude of rotor oscillations with increasing median of the slip length 0.1 ≤ M ≤ 1.0; σ ε = 0.1, σ l s = 0.1 in the case of a PCB with wide annulus; a = 0.2, σ = 1, λ = 1, α = 1, K z = 10 and D a = 1.

1 a
(p − p a )rdr(45)where p a is the ambient pressure above the stator.The force on the stator is obtained in the formF(t) = π A(g(t)) + dg dt B(g(t))(46)using the force integral in (45) and pressure in (38).Expressions for A and B are defined asA(g, λ , l s , β ) = (1 − a 2 )(p I − p a ) + 2(p O − p I )