Fluid dynamics of the slip boundary condition for isothermal rimming flow with moderate inertial effects

Motivated by evaluating coating oil ﬁlms within bearing chambers in an aero-engine application, an analysis is presented for the ﬂuid dynamics relevant in their dual capacity as both the coolant and lubricant in highly sheared ﬂows that may approach microscale thickness. An extended model is developed for isothermal rimming ﬂow driven by substantial surface shear within a stationary cylinder. In particular, a partial slip condition replaces the no-slip condition at the wall whilst retaining inertial effects relevant to an intrinsic high speed operation. A depth-averaged formulation is presented that includes appropriate inertial effects at leading-order within a thin ﬁlm approximation that encompasses a more general model of assessing the impact of surface slip. Non-dimensional mass and momentum equations are integrated across the ﬁlm depth yielding a one dimensional problem with the a priori assumption of local velocity proﬁles. The ﬁlm ﬂow solutions for rimming ﬂow with wall slip are modeled to a higher order than classical lubrication theory. We investigate the impact of wall slip on the transition from pooling to uniform ﬁlms. Numerical solutions of ﬁlm proﬁles are provided for the progressively increased Reynolds number, within a moderate inertia regime, offering evaluation into the effect of ﬁlm slippage on the dynamics of rimming ﬂow. We ﬁnd that slip allows non-unique solution regions and existence of multiple possible steady state solutions evaluated in transforming from smooth to pooling ﬁlm solutions. Additionally, boundary slip is shown to enhance the development of recirculation regions within the ﬁlm which are detrimental to bearing chamber ﬂows.


I. INTRODUCTION
Thin film flows are widely used in engineering applications, notably as a lubricant between moving surfaces and in providing a thermal conduit for heat exchange. The specific application we have in mind is in the modeling of oil rimming flows within a bearing chamber by the thin film approximation. When subjected to a high surface shear, complex flow dynamics develop that negate a classical lubrication approach and require high order modeling techniques. A standard element of conventional analyses is in the no-slip assumption as a boundary condition; however, through surface treatments, there is growing attention to situations such as progressively thinner films that may invalidate the no-slip condition, as described by Lauga, Brenner, and Stone. 1 This paper investigates the effects of the slip boundary condition for thin fluid films driven by surface shear.
A two dimensional stationary bearing chamber, see Fig. 1, is considered whilst driven by an imposed shear stress at the free surface causing thin films of oil to develop on the chamber surface. In maintaining a steady state, typically, fluid volume is replenished by droplet impacts from above the film, with an equal mass subtracted via an extraction sump, but is outside of the current scope. These conditions describe a rimming flow, seen in Ruschak and Scriven's work on cream separation, 2 or an equivalent ceramic coating of pipes by Menekse, Wood, and Riley. 3 The addition of a slip boundary condition can be thought of as a first step in modeling film flow within a porous cylinder. Beavers and Joseph 4 show the shear stress in the film at the wall is proportional to the difference in tangential velocity between the film and the porous flow at that point. The slip condition used in this model has the film shear stress at the wall proportional to the speed of the film at the wall (the Navier slip condition) which is often applied to micro-or nanoscale fluids.

Physics of Fluids
For rimming flows, the no-slip condition has been questioned only recently. On the topic of bearing chambers, Kay 5 exhibited film heights of 100 µm, where surface effects of the chamber wall may influence the developing film. We are not aware of the slip condition being used for large scale fluid rimming flows, though a comparative fluid flow can be found in the work of Barber et al. 6 or Bowles and Ducker 7 for air lubricated bearings. The surface to volume ratio leads to surface effects dominating the film dynamics. An air lubricated thrust bearing with the Navier slip condition is presented by Bailey et al., 8 where the key parameters are comparable for extremely small fluid thickness and the intermolecular distance.
Boundary conditions for flow in a highly sheared bearing chamber, as depicted in Fig. 1, typically adopt the widely accepted no-slip boundary condition, used from as early as 1738 by Bernoulli for the fluid-solid interface. 1 A partial slip boundary condition was an early notion by Navier according to Neto et al.,9 in which a liquid may slip across a solid surface whilst encountering an opposing force relative to the velocity of the solid arising from friction. This partial slip is generally defined by a slip length, widely accepted in quantifying the slip of liquid at the solid surface, analogous to defining the depth of penetration for the velocity profile into the solid and is the type of slippage used in this paper. Furthermore, there is also the stagnant layer boundary condition and full slip boundary conditions for completeness. 9,10 Research on the slip condition is heavily linked to surface treatments on the solid surface, characterised by their hydrophobicity as either hydrophobic or hydrophilic surfaces, and is explored for a partial slip boundary condition by Choi et al. 11 They found that a hydrophobic surface produced longer slip lengths for water than a hydrophilic surface, up to several µm. Further studies on hydrophilic and hydrophobic surfaces on Pyrex were made by Cottin-Bizonne et al. 12 While hydrophilic surfaces were bounded by the no-slip condition, hydrophobic surfaces had superior slip lengths to several µm, conjecturing nano-bubbles may coat hydrophobic surfaces and improve slip.
An interesting proposal for the slip boundary condition is in hydraulic fracturing of shale, where Javadpour, McClure, and Naraghi 13 documented slip lengths in micro to nano flows for organic pores of a shale matrix. Their work highlighted the importance of including and understanding the effect of slip in shale nanopores and hydraulic fracturing. Hu and Granick's review briefly mentioned that liquid slip may occur at very low shear rates and the importance of identifying surface effects as tribology approaches nanometer film thicknesses. 14 The Navier slip model is defined as a constant slip model, with a linear relationship between the tangential shear rate and the interface velocity differences, in which the slip length acts as the constant of proportionality. The true condition has yet to be fully verified for all applications due to the complexity involved with individual effects, for instance, surface roughness of the interface or the wettability of the surface.
The work by Ruschak and Scriven 2 on fluid films driven by high rotation speeds for a horizontal cylinder equated to a solidbody approximation for a suitably low fill volume and negligible gravity. However, under conditions where the fluid is not driven sufficiently, pooling and recirculation of the fluid is possible and leads to the breakdown of the solid-body approximation. The film is then no longer at uniform height as regions of film acceleration and deceleration due to gravity lead to thinning and thickening, known as pooling. The original approach by Moffatt 15 used a lubrication approximation for a two dimensional coating flow with a thin film approximation, and for a low fill fraction with surface tension and inertia negligible, they derived a flux balance equation. This equation predicts smooth film profiles up to a critical mass, beyond which fluid pooling invalidates the lubrication approximation. For rimming flow with a lubrication approximation, Johnson 16 arrived at the same cubic equation of Moffatt 15 whilst including a balance between gravitational and viscous forces to leading order.
On the premise of lubrication theory for rimming flows, the effects of neglecting surface tension and inertia were investigated by Tirumkudulu and Acrivos. 17 They showed discontinuous regions between the steady-state film height profiles, leading to an infinite pressure gradient. Upon incorporation of a hydrostatic pressure term at high order, they found reverse flow regions (recirculation zones), 15,16 which were not possible to obtain using classical lubrication theory.
Ashmore, Hosoi, and Stone 18 detailed the effects of surface tension on rimming flow in terms of a ratio of gravitational to viscous forces, λ. They found that the flow can be divided into 3 distinct regimes related to gravitational effects, defined by λ. In the first regime (0 < λ ≤ 2), viscous forces dominate the flow, resulting in smooth profiles. In the next regime (2 < λ ≤ 5), where gravity terms are of the same order as viscous terms, discontinuous profiles are found. In the last regime for λ > 5, surface tension becomes important and the fluid pools at the cylinder base due to high gravitational effects.
The modified lubrication equation of Acrivos and Jin 19 involves first order gravitational terms to model the hydrostatic pressure gradients that become considerable in the sharp gradient regions of fluid pooling and recirculation. The inclusion of higher order terms is also explored by Wilson, Hunt, and Duffy 20 by analytically calculating the smoothing from first order gravity on the Stokes equations to eliminate discontinuities at critical rotation rates, and also show the significant potential of surface tension.
The fundamental balancing of surface tension, gravity, and rotation is at the core of coating flows, with a good review of the subsequent experimental work on coating flows can be found in the work of Evans, Schwartz, and Roy, 21 whose later lubrication model 22 modeled the axial features of drops, undulations, and rings presumably due to surface tension.
The combined effects of incorporating surface tension and pressure are presented by Villegas-Díaz, Power, and Riley 23 for their stabilising effects on the presence of discontinuous profiles. Earlier work by Villegas-Díaz et al. 24 utilised a combination of surface shear and rotational driving forces, studying their effect on the stability both numerically and analytically, where the viscous and gravitational effects were of the same order.
In classical lubrication theory applied to rimming flows, such as in Refs. 15 and 16, the typical assumption of negligible inertia leads to a zero Reynolds number. This is contrary to later studies 25,26 which cover a large range of Reynolds numbers, analyzing the instabilities caused by inertia and highlighting complex pattern formations such as shark-teeth. The work of Hosoi and Mahadevan 26 extended 15,16 for surface tension, pressure gradients, and inertial effects in response to earlier work in lubrication theory by Thoroddsen and Mahadevan. 25 The addition of inertia in a higher order lubrication model demonstrates significant instability in the profile of the film, leading to values of the velocity profile that accounts for regions of recirculation in rimming flow.
The role of inertia on the Moffatt-Puknachov coating flow is described by Kelmanson 27 for a range of surface tension limits; a critical Reynolds number, above which no steady-state solution exists, is calculated due to the destabilizing role inertia plays.
The first order reduced Reynolds number model of Noakes, King, and Riley 28 focused on inertial effects for both rimming and coating flows in a three dimensional study, where complexity of nonlinear inertia at higher reduced Reynolds numbers can be posed similar to those of a weaker inertial form. The loading of fluid both external and internal to a horizontal cylinder in the analytical approach by Leslie, Wilson, and Duffy 29 calculates the existence of a critical loading threshold due to cylinder rotation before dewetting will occur. The constraints of fluid loading or rotational speed lead to a limitation of the sustainable speed or fluid volume, respectively.
A model to the first order of the viscous, gravitational, inertial, and capillarity effects was proposed by Pougatch and Frigaard, 30 noting minor changes to the location and amplitude of the film thickness profile in two and three dimensional flow. They suggest that gravity provides a stabilising force, whereas surface tension and inertia tend to destabilise the film for lower levels of inertia. On this instability caused by inertia, the work of Benilov and O'Brien 31 argued that surface tension terms provide a smoothing effect. Upon balancing the inertia and surface tension, instabilities arising from inertia are suppressed by surface tension leading to a smoothed solution for short wave disturbances, whereas the opposite is true for long scale disturbances which are counteracted by viscosity, see the work of Benilov 32 or Benilov and O'Brien 31 Recently, Kay 5 evaluated a thin film approximation for the effects of inertia, where low Reynolds numbers exist within the realm of lubrication theory and high Reynolds numbers corresponded to uniform film profiles. The instability of inertia is seen to be dampened by the presence of gravitational and Marangoni effects in a first order lubrication model. 33 The addition of inertia or viscous drag shifts the fluid pool in the direction of rotation within the lubrication approximation. Impregnating the fluid with insoluble surfactants, Kumawat and Tiwari 34 model rimming flow with the effects of gravity, viscosity, surface tension, and inertia. The inertial instability that appears in the form of oscillations is dampened by higher Marangoni stresses that arise from inhomogeneous surfactant gradients at the surface.
Very small cylinders may suppress this instability by surface tension using the full Navier-Stokes equations; however, the resulting unstable disturbances can be absorbed within a shock according to Benilov and Lapin 35 while the smoother films are unstable. This stabilising role of surface tension is in the regularization of rimming flow, where regularizating terms are typically intended to eliminate physically meaningless solutions and are also beneficial for stability analysis. 36 Solutions for the rimming flow problem are feasible with methods such as direct numerical simulation of the governing equations, as seen in the work of Benilov, Lapin, and O'Brien. 37 However, an alternative technique was proposed by Kay, Hibberd, and Power; 38 by the use of hydraulics theory applied to the lubrication concept via a depth-averaged formulation, retention of inertial terms to leading order is viable. On using the Karman-Polhausen depthaveraging of the governing equations between the flow boundary layers, the system is reduced in dimensionality obtaining the commonly known Shkadov model when using a parabolic velocity profile.
The technique of depth averaging the governing equations has been performed for thin film approximations on both planar flows and rimming flows. The work of Nguyen and Balakotaiah 39 formed two coupled depth-averaged equations from the mass and momentum equations by solving for the mean film properties, in the a priori selection of a local velocity profile to be substituted into the resulting thin film model. The selection of a parabolic film profile was adapted by Nguyen and Balakotaiah 39 and Chang, Demekhin, and Saprikin 40 on the basis of Stokes flow where Re ≪ 1 is having an exact quadratic solution.
A good formation of a depth-averaged model can be found in Benilov's work for a single depth-averaged governing equation, 41 with a quadratic velocity profile on an incline planar geometry. Kay et al. 38 covered the depth-averaging method in a rimming flow context, with a priori use of quadratic profile or cubic velocity profile when incorporating cylinder roughness for reducing the dimensionality of the Navier-Stokes equations.

II. THE ROLE OF THE BEARING CHAMBER
The operating condition within bearing chambers varies due to the nature of progressive technological deployment of aeroengine components. The chaotic and turbulent flows are taken as time averaged to develop steady profiles, with the imposed surface shearing as an average due to the transience and turbulent nature of the two phase flow. The coating flows of Evans, Schwartz, and Roy 21 report very large thicknesses of 0.11 m in a heavily loaded case; the typical properties for the films of interest within the chamber are similar to those of Kakimpa, Morvan, and Hibberd, 42 shown in Table I. Relevant studies on bearing chambers that are experimental or numerical provide a breadth of the operating conditions expected, of which report the thinnest films of 10 µm. 33,[43][44][45][46][47] Unlike the bearing chamber geometries and flow characteristics, the research into slip length predictions is less consistent due to the intrinsic difficulty of measuring nano-and microscale effects. The flow of water across hydrophobic and specially engineered super hydrophobic surfaces poses significant difficulties to qualitatively predict interface velocities, in part due to wettability and surface roughness. 1 In the presence of an air only interface, an infinite slip length is expected, 9 but interfaces of the substrate to support the 0.024 µm-0.063 µm Tan et al. 57 1.1 at ε = 0.5 3.8 at ε = 0.9 fluid reduce this with a variety of surfaces and fluids detailed in Table II for observed slip. Alongside these coated substrates, the free flow of fluid over a porous medium details a correlation between the internal and external velocities. This correlation typically assumes a constant slip length or fraction of the bulk velocity, see Table III. We assume a surface coating applied on the cylinder will act to provide fluid slip typically l * s ∼ 50 µm on the free fluid for this study, though a variety of slip lengths are used to observe their influence on films.
We present a new effort in rimming flows by applying a slip boundary condition to the inner surface of a cylinder, using a depthaveraged approach to account for inertial features. The mathematical limits of cylinder surfacing for drag reduction are detailed alongside the physical slip potential for a bearing chamber. Section III details the formulation of the depth-averaged model with slip. Film profiles are calculated by numerical differentiation using a pseudospectral method, see Refs. 58 and 59, that is beneficial to the periodic rimming flow problem for fast convergence and as a global approach. Section IV presents films consistent with the previous literature with further detailing on the slip boundary condition repercussions studied. In Sec. V, our main conclusions and summaries are conferred.

A. Film flow in cylindrical coordinates
We assume a spatially two dimensional configuration consistent with a simplified bearing chamber geometry, shown in Fig. 1. A liquid film wholly coats the chamber interior of constant distance r 0 from the centre while being subjected to a constant surface shear stress, τ * a . The superscript * indicates a dimensional variable, while t * represents time. The liquid film has a constant density ρ and viscosity µ, being treated as an incompressible fluid with a pressure p * a at the film surface. In cylindrical polar coordinates, the angle θ * is measured positive from the vertical downwards direction anticlockwise, with r * measuring the radial distance from the cylinder centre. In this coordinate system, gravity acts in the downwards direction defined as g = (g θ , gr), with the components defined as gr = g cos θ * and g θ = −g sin θ * for a per unit mass gravitational force.
The film free surface is described by r * = r 0 − h * t * , θ * , with a surface tension coefficient σ * . The two dimensional film has a velocity field u * = u * θ , u * r and a pressure field p * . The unit tangent and unit normal vectors of the free surface are defined as l = l * rr +l * θθ and n = n * rr + n * θθ with in which N * and H * are defined as Therefore, the governing equations for the conservation of mass and momentum for the film under the outlined conditions are given as With the dimensional gradient ∇ * and D * /D * t * , the dimensional convected derivative terms are defined as At the film interface, at r * = r 0 − h * θ * , t * , there acts a stress Tensor T * ensuring continuity along with a kinematic equation described as The film curvature κ * at the free surface interface is thus defined by Typically, the wall surface boundary condition is described by the no-slip condition; we define a slip boundary condition for a stationary wall surface, given in (12) along with the no-penetration condition of (13) at r * = r 0 , In order to ensure that the quantity of fluid volume inside the cylinder remains consistent, a variable A is introduced, known as the filling fraction, presented as In the studies of Ashmore et al., 18 the filling fraction was used to reduce the governing equations by acting as a small parameter; we use the film aspect ratio but retain A for comparisons to other studies and ensure film mass continuity.

B. Non-dimensional equations
The following scalings are introduced to non-dimensionalise the model for the further study. The pressure field is scaled to retain azimuthal pressure gradients by the dynamic pressure, ρU 2 0 , leading to a potentially stabilising mechanism. 17,18,38 We also introduce the scaling ε = h 0 /r 0 , known as the film aspect ratio, in which h 0 is the characteristic film height. The subsequent scalings for the model are The slip length l * s is scaled by the film height h 0 due to the nature of slip length as a representation of the film velocity penetration into the boundary in a partial slip condition. The typical film height, h 0 = 500 µm, leads to l * s being microscale. 42 Some key dimensionless factors to the film flow are the capillary number, Ca, the Reynolds number, Re, and a ratio between the gravitational to viscous forces, λ, seen as Re/Fr In moving to the non-dimensional coordinate system, the governing Navier-Stokes equations (5) and (6) separated into radial and azimuthal components are The tangential stress component along with the normal stress component of the interface boundary condition (9) at y = h becomes The kinematic condition (10) at the fluid free surface is 22) and the curvature at the film surface (11) is where ′ , ″ , and ‴ represent first, second, and third order derivatives with respect to s. The terms (3) and (4) non-dimensionalise as The filling fraction of the fluid within the cylinder (14) equates to The boundary conditions of (12) and (13) for the stationary wall at y = 0 reduce to

C. Thin film equations
Within the aero-engine bearing chambers of interest, oil films may comprise inertial effects of modest intensity which can be of the leading order, corresponding to εRe ∼ O(1) within a thin film approximation of ε ≪ 1. The model in this paper is appropriate to order ε when the non-dimensional parameters with moderate inertial effects are in the following range: Imposing these limits allows the use of the classical lubrication theory in this formulation. When taken to O(1), (18) will yield a constant pressure throughout the film. 38 This is inadequate in regions of large film gradients due to hydrostatic effects in pooling solutions leading to sharp pressure fluctuations. To include these effects, gravitational terms of O(ε) must be kept, resulting in an O(ε) accurate model. For consistency, the inertial terms (centrifugal) to O(ε) must also be included; we note that terms of O(ε) are factors of h ′ , allowing for their elevation to leading order in areas of large film gradients. In satisfying this approximation, the continuity equation (17) The radial momentum equation (18) equates to and azimuthal momentum equation (19) can be written as The boundary conditions (27) and (28) at the wall surface y = 0 to O(ε) are given as and the filling fraction (26) is At the film free surface, the kinematic condition (22), normal stress (21), and tangential stress (20) components at y = h become The curvature of the free surface (23) becomes

ARTICLE scitation.org/journal/phf
where O ε 2 have been retained in the normal stress condition due to h ″ acting as a stabilising mechanism that will be explored later upon depth averaging.

D. Depth averaging
Due to the retention of inertial terms approaching leading order, an appropriate resolution of the inertial terms is vital. The chosen method reduces the dependency on the velocity by depthaveraging, 38 giving s and t as the only independents. By the integration of the radial momentum equation across the film depth, the distribution of pressure across the film is obtained for the azimuthal momentum equation. Upon integration across the film height, the continuity equation (30) becomes with a dimensionless variable per unit length, q, defined as the volume flux The pressure distribution to leading order across the film is obtained upon integration of the radial momentum equation (31). the integration constant is evaluated at the interface from the normal stress component (37). This leads to a pressure distribution through the film that is primarily hydrostatic, with additional effects from surface tension and centrifugal forces. The depth-averaged azimuthal momentum equation (32) becomes Use of the pressure from (42) evaluates the pressure gradient term. Surface tension terms are retained up to O ε 3 due to the importance of h ‴ terms in pooling or shock solutions, where they provide a stabilising effect according to Refs. 18 and 23. An additional effect of surface tension is that short-wave perturbations are dampened when obtaining numerical solutions, conferred in the work of Kay et al. 38

E. Resolving velocity profiles
In a highly sheared flow within a bearing chamber, film thicknesses can be micro-scale while inertial effects can still be substantial. In the case of negligible inertia, i.e., the Stokes flow for Re ≪ 1, a quadratic velocity profile provides an analytical solution to the Navier-Stokes equations. For a constant flux q, if the pressure gradients, surface tension, and inertia are negligible, there is a balance between the gravity and the net shear at the free surface and at the wall. Following the depth-averaging of the governing equations, the integrals of the azimuthal momentum equation (43) must be solved. Based on this, a quadratic form of the velocity profile is assumed to be O(1) for the velocity profile, having also been used under similar conditions, 38 u = u 0 + εu 1 , where components u 0 and u 1 in (44) are given as Using the definition of the film flux q (41), the surface shear condition (38) and the cylinder surface slip condition (33) to first and leading order of the velocity coefficients become This set of simultaneous Eqs. (47)-(49) leads to a determination of the velocity coefficients of (45) and (46) defined in terms of h and q as outlined in Appendix A.
Having defined the azimuthal velocity profile u, in terms of the first and leading order components, upon integration of the continuity equation prior to depth-averaging (40) gives the radial velocity in terms of the azimuthal velocity (44) Using the chain rule to differentiate the velocity components, we obtain the derivatives of each coefficient with respect to h ′ and q ′ .

F. Film model equations
Following the substitution of (43) with the velocity polynomials (44) and (50), each integral is evaluated in Appendix B. A modified Reynolds number, Re * = εRe, is introduced, known as the reduced Reynolds number with Re * ∼ O(1) according to (29).

ARTICLE scitation.org/journal/phf
The momentum equation can be represented by contributions of the forces included where the full coefficients of α, β, γ, ω, ψ, N 0 , and N 1 are given in Appendix C. The key features are distinctly inertia, gravity, surface tension, pressure gradients, and the net viscous shear. 38 Here, we have retained the εh ′ terms because of the contribution at local levels for sharp film gradients. Further consolidation of the terms reveals the final model defined by the momentum equation (51), continuity equation (40), and volume conservation equation (35). These coupled temporal equations can be solved for the film height, h(s, t), and film flux, q(s, t).

G. Numerical method
For obtaining solutions to the system of equations proposed in (51), (40), and (35), a numerical approach is required to solve the system of non-linear differential equations as described in this section.
We solve the above system using a pseudospectral method based on a transformed uniform discretisation and Fourier differentiation matrices, i.e., a global approach for resolving spatial differentials. For resolving temporal terms, a first order fully implicit forward difference approximation is used. A steady-state solver was used to explore the parameter space with a transient solver determining the nature of the solutions. An initial Fourier differentiation matrix as proposed by Trefethen, 59 based on a uniform distribution of N nodes in the domain of η = [−1, 1], is created. A stretching is applied using a hyperbolic sine coordinate transformation, which clustered nodes near the formation of capillary ripples at the centre of the domain, Here, parameter c is used to control the nodal distribution from uniform, for c = 1, to densely clustered in the region where capillary ripples are formed at c > 1. The definition of L is half the domain length such that L = π. In transforming the Fourier differentiation matrices, Boyd 58 details a general derivative transformation such that d ds whereḟ,f , and ... f are the 1st, 2nd, and 3rd derivatives of (52) with respect to η. In order to maintain a consistent volume of fluid within the cylinder, the filling fraction is calculated using the trapezium rule approximation within the grid where hj is the film height at the j node of the domain and ∆sj,j +1 is the spacing between the jth and jth + 1 nodes with hN = h 0 due to the cylindrical domain. Having established the coupled azimuthal momentum, continuity, and volume equations that govern the film listed above, a steady state equation can be formulated, with the depth-averaged continuity equation establishing the film flux or filling fraction as a constant value. In solving the system of non-linear equations, the package from MATLAB for handling systems of equations fsolve is used for solving the N + 1 or 2N coupled equations for determining the film height and film flux at the N nodes.

H. Film profiles within classical lubrication theory
This section focuses on the classical lubrication theory limit, which within ε ≪ 1 leads to a gravity-viscosity balance in the film with all other forces (inertia, surface tension, and pressure gradients) neglected. In this case, the governing equation with a partial slip boundary condition reduces to This equation is important as the additional forces which are neglected can be considered as perturbations so that a lot about the film behavior can be determined from it. It is also algebraic rather than differential so can be analyzed simply. Analysis of (57) was first carried out by Moffatt 15 and later by Ref. 16 for the case of a rotating cylinder with no-slip conditions. By exploiting the boundedness of the sine function, they determined that the film flow was separated into three regimes corresponding to flows below, at, and above a critical value of height and flow rate (flux) which correspond to the values of h and q at s = π/2. Upon treating λ and τ as constant values, the depressed cubic equation (57) is readily solvable with the slip condition included, and the critical values are Typically, qc is used as the parameter to distinguish between the regimes though hc can be equivalently used. Previous studies with a no-slip boundary condition yield the well-known critical film flux of qc = τ 3 (6λ) 2 and the critical film height hc = τ/λ. An asymptotic analysis for an increasing slip length demonstrates that the critical film height tends toward hc = τ/2λ and the critical flux approaches qc = ∞ on an idealized drag free model. The three distinct classifications of film profiles for a specified film flux are termed sub-critical, critical, and super-critical, defined as 1. Sub-critical case, q < qc. Two branches of real and positive solutions exist. One branch represents a completely wetting solution across the entire domain, in which the cylinder is completely coated in a fluid. The second branch represents a non-physical solution in the region 0 ≤ s ≤ π, with an unbounded film height at the limits of this region. 2. Critical case, q = qc. The previously mentioned two branches converge at the point s = π/2, forming a completely  wetting solution branch connected to an unbounded solution that exists within 0 ≤ s ≤ π. At this value of the film flux, the film height forms a peak known as the critical film height, hc, and provides the maximum fluid loading for a continuous profile under the lubrication approximation, defined as the critical filling fraction, Ac. 3. Super-critical case, q > qc. No solutions exist in this region that completely coat the cylinder's internal surface for lubrication theory. Two solution branches exist with a disconnected region focused at s = π/2 in which no real and positive film solution can exist. Inclusion of higher order effects has been shown to bridge this discontinuous region outside of the typical lubrication theory assumptions.

Physics of Fluids
Within these classifications, a single physical solution exists in the region s ∈ [−π, 0], while for s ∈ [0, π], the solution is dependent upon the ratio of q/qc with physical, unbounded, or discontinuous solutions. These are explored in Fig. 3 and evaluated with moderate inertial effects to form fully coated rimming flows. The variation in critical values with increasing slip lengths are given in Table IV. As the slip length increases, the effects on the critical film values asymptotically approach their limiting values. The inclusion of a wall slip is found to limit the maximum sustainable filling fraction to half of the value in a no-slip case. It should be noted that for a fixed filling fraction A or film flux q, an increasing slip length ls alters the transition from sub-critical to super-critical cases, as discussed in Sec. IV. It can be seen that the diminishing film critical height hc becomes apparent when ls ∼ O(ε), as shown in Fig. 2(a), closely matching the trend for the critical filling fraction Ac, as shown in Fig. 2(c). By contrast, the film flux is seen to exponentially increase with ls in Fig. 2(b).

A. Smoothed solutions subjected to inertia and slippage
In accordance with classical lubrication theory, the existence of fully coating flows beyond q > qc is introduced by higher order effects when taking surface tension to be of first order such that ε 3 Ca ∼ O(ε). This provides a smoothing effect, bridging across the non-physical solution region in the super-critical profile as detailed by Benilov, Benilov, and Kopteva, 60 and Kay et al. 38 This smoothing allows for completely wetting film profiles when within the defined super-critical realm, as shown in Fig. 3(c), unattainable within lubrication theory. The presence of solutions within the super-critical class allows for the presence of shock-type solutions to be explored. This is categorised by a film profile that occurs on the smooth solution with a peak, typically located on the cylinder's rising side that connects to the unbounded solution found in sub-critical cases as in Fig. 3(a). The addition of inertial effects has smoothed the peak film height, seen at the critical film profile of Fig. 3(b), with the solution branches intersecting.

B. Increasing gravitational effects with negligible inertia
In treating the inertial effects as negligible, films with filling fractions equating to A > Ac are explored. In an equivalent sheardriven rimming flow by Ashmore et al., 18 they studied the effect of gravity-dominated solutions in the formation of pooling films. They varied the equivalent parameter of λ for the ratio of gravitational to viscous forces. Next, we will study the effect of gravity on film profiles with slip. At the first instance, the slip values are chosen to be small to obtain film heights comparable to Ref. 18. As noted previously, increasing λ transitions the film profile from a smooth rimming flow toward a deepening pool solution at the cylinder base. In accounting for the differences of rotational and shear driven flow, the film profiles in Fig. 4 are comparable to Ref. 18 for both high and low gravitational effects consistent with those reported. When the gravitational effects are weak (λ = 1 − 2), the film is seen to spread relatively uniformly across the cylinder, as in Fig. 4. When λ > 5 however, we observe a solution pooling deeply at the bottom, where much of the film coating the cylinder becomes exceedingly thin. Therefore, we expect the effect of slip to become substantial in the latter case. In Sec. IV C, we will investigate how an increasing slip length influences the transition from a smooth profile to a pooling solution.

Physics of Fluids
The pooling solutions, typically where higher film heights exist around the cylinder base, are synonymous to a recirculation region, highlighted in Fig. 5 for the flow streamlines of Fig. 4 with λ = 1 and λ = 3. The bulk of moving fluid from the thin region is driven over the recirculation by driving the air shear.

C. On the balance of slip and gravity
We will now look at how the transition between smooth and deeply pooled solutions (presented in Fig. 4) is affected by the slip the film outside of 0 < s < π, with the fluid mass accumulating into a raised peak observed at s = π/2. When the slip length increases toward a near drag free setting of ls = 1, a substantial thinning of the film is observed, with a reduction of the thickness up to half of that reported with ls = 0.01, and a significant deepening of the film for 0 < s < π/2, highlighting the importance of surface effects at ls ∼ O(1). Figure 6(b) shows the same collection of film profiles with different ls, but with a relatively stronger gravitational force (λ = 2). Here, we obtain a pooled solution for all slip lengths; however, increasing ls results in a significant thinning of the film in areas of uniform film height (for −π/4 > s and s > π/2). The pooling of fluid within the domain is attributed to the fixed filling fraction, thereby resembling a pooling solution even at low gravitational forces. In Sec. III H, the largest changes to critical film values are observed around O(ε), with Figs. 6(a) and 6(b) demonstrating a significant divergence from those of Ashmore, Hosoi, and Stone 18 as ls transitions beyond O(ε). A further increase in the slip length encourages further thinning of the flow profile, leading to excessively further pooled flows within the low gravitational force range. It is within the limit of Re * ≪ 1 and λ ≤ 5 that one observes the most significant effects of wall slip, where even small slip lengths, l * s ∼ O h * , result in observable pooling and thinning.
For exceedingly larger λ, the subsequent effect of ls is seen to diminish. At λ = 10, even the longer slip lengths of ls ∼ O(1) lead to negligible changes in the film profile compared to the low slip conditions. This is due to the overwhelmingly strong gravitational forces affecting the flow, which results in exceedingly low film heights across the domain outside of the pool at −1 < s < 1.
The controlled fixed filling fraction, such that A > Ac for all cases, allows for super-critical film states due to stability from the addition of surface tension, allowing for deeply pooled and shock formations to exist.

D. Inertial effects for transitioning to smooth solutions
It is within the realm of A > Ac that the most interesting features can be explored for solutions subjected to inertia undefined by lubrication theory. Previous studies have highlighted the effects of increased inertia with the no-slip condition, presenting rich features on film solutions. 31,38 Upon exploration of the parameter space for the reduced Reynolds number, Re * , and film flux ratio, q/qc, the flow with slip at the wall is compared to its equivalent no-slip flow. The influence of slip length on the critical filling fraction, Ac, leads to an adjusted filling fraction ratio A/Ac being maintained for comparable results on the flux ratio q/qc for no-slip and partial slip profiles.
Summarising a rimming flow of no-slip film subjected to low inertial effects, Re * = 0.01 for a fixed filling fraction of A > Ac; a typical pool film develops as seen in Fig. 7(a). Greater inertial effects lead to the eventual smoothing of the flow, as shown in Fig. 7(f) Re * = 100 . Moreover, for the pooled solutions at q/qc < 1, the fluid accumulates on the rising cylinder wall exhibiting a sharp front; see Fig. 7(a). Within the regime of the lubrication approximation (29), the increase in Re * leads to the formation of small wavelength disturbances at the sharp front, observable first in Fig. 7(b), which are eventually smoothed due to the presence of surface tension, inertial, and pressure gradients. These small wavelength disturbances, identified as capillary waves, are seen to grow for greater inertial effects, encapsulated in Figs. 7(c) and 7(d), where the sharp front seen in low Re * is replaced by these capillary waves. Upon moving toward more dominant inertial flows, of Re * > 15, the small capillary ripples are seen to be greatly diminished in Fig. 7(e) and have become non-existent in Fig. 7(f). The limitation of (29) could be breached at higher Re * unless the film remains smooth, where h ′ ≪ 1 which leads to limited inertial effects.
The possibility of smoothed solutions for super-critical filling fractions arises from the leveling of the supplementary inertial terms, therefore satisfying (29) as the film mass is almost uniformly distributed across the domain. Addition of a slip length, ls = 0.01, is examined for the behavior in obtaining pool solutions for q/qc < 1, shock solutions at q/qc = 1, and solutions for q/qc > 1 that exist outside of lubrication theory.
In establishing the initial solutions for slip flow, displayed in Fig. 8(a), the sharp front at the start of the solution pool equivalent to that of no-slip is found on the rising cylinder side, existing within the pool solution regime for q/qc < 1. While containing steep fronts, akin to a shock solution, regions of recirculation on the rising cylinder side within a pool are possible, seen in Fig. 5. Upon advancing Re * , the pooled solution depth is diminished as the fluid is spread across a larger portion of the rising cylinder wall. As Re * → 5, a number of small waves begin to form in front of the previously sharp front, see Fig. 8(c), which has been reduced in depth further due to inertial effects causing a smoothing across the film. These small waves reach a peak disturbance upon the film around Re * ∼ 10, where a non-insignificant portion of the film is under the influence of these waves, known to be an artefact of the additional surface tension terms, described as capillary ripples, which can pose challenges in resolving numerically.
At Re * ∼ 10, it is possible that (29) is infracted due to h ′ not being negligible. Increasing inertial effects of Re * > 10, a transition toward a smoothed solution occurs. The initially pooled solution becomes further smoothed due to inertial effects, with fluid mass rising further along the cylinder wall along with a stretching of the capillary ripples, shown in Fig. 8(e). This smoothing is culminated in an almost uniform rimming flow of Fig. 8(f), where h ′ ≪ 1; therefore, (29) is satisfied as inertia becomes minimal.
For Re * > 20, the inertial effects allow for smooth solutions to exist for q/qc > 1, contrary to classical lubrication theory's detached solution branches. Equally, the effects of inertia allow for film flux ratios of q/qc < 1, seen to be within the sub-critical regime for films that maintain a super-critical filling fraction A > Ac.
A large slip is shown to increase the rate at which the transitional region between pooled and smooth solutions appears, with the onset of capillary ripples occurring at lower Re * . This presents an interesting case where the increased inertial effects are restricted in their ability to smooth the pool solution, requiring larger effects to raise the fluid along the cylinder wall and begin the eventual removal of the capillary ripples. The overall trends on the pooled, shock, and smooth solutions are typically unchanged in appearance, notably as inertial smoothing becomes the dominant mechanism for Re * > 20; only minor ripples are apparent, which dissipate shortly thereafter. With larger slip lengths, capillary ripple formations are exceedingly harder to resolve numerically.
This highlights the benefit of the slip boundary condition in forming a levelled profile, which is not subjected to recirculation regions seen at low inertia, and hence avoids problems in the maintenance of fluid properties.

E. On the presence of non-unique solution branches
By the implementation of high order surface terms into (51), a region within a parameter space has been shown to contain nonunique solutions. 38 Their added wall friction factor highlighted novel characteristics of the flow dynamics, having employed a cubic polynomial for improvements in the depth-averaging method. In a similar fashion, the implementation of a slip boundary condition was explored across a parameter space, shown in Fig. 9.
Next, we will investigate the solutions at varying filling fraction A/Ac but at a fixed Reynolds number. We will show that non-unique solutions exist for a flow with wall slip, while the no-slip flow solution remains unique for the same parameters. An initial profile of a fixed low filling fraction passed through the steady-state solver achieves a smooth rimming flow profile seen in Fig. 9(a). Let us first focus on the no-slip case (solid line). When the filling fraction is increased from zero, the smooth solution develops toward the critical flux ratio q = qc in Fig. 9(b) with the fluid pooling on the lower rising side of the cylinder. By raising the filling fraction, transitioning above q > qc and A > Ac, seen in Fig. 9, shows further fluid accumulating on the rising cylinder wall in Fig. 9(c). This development can be understood by the smooth solution profile connecting the bounded and unbounded solutions from lubrication theory of (57) for the sub-critical classification of q < qc. Subsequent increments of the filling fraction lead to a minor reduction in the fluid flux for film profiles which remain consistently within the sub-critical classification. The increased fluid within the domain accumulates within a deepening pool at the lower half of the rising cylinder wall, captured in Fig. 9(d); further addition to the filling fraction leads to capillary ripples becoming more prevalent. This parameter space for no-slip contains only one unique solution for any A/Ac. The hump for the no-slip parameter space is enlarged with the inclusion of low reduced Reynolds. Next, we will retrace the filling fraction-space with a moderate slip length, ls = 0.5. The profile at low filling fraction is show in Fig. 9(e) and is very similar to the noslip profile for the same A/Ac. The profiles approaching A = Ac bear similarities, beyond which the parameter space and profiles differ considerably. When increasing the filling fraction above the critical (A Ac > 1), the flux ratio q/qc increases toward a maximum attainable value, for which the profile is shown in Fig. 9(g). Careful adjustment of the filling fraction then leads to a reduction in the fluid flux to Fig. 9(h), beyond which increases in the filling fraction fail to resolve solutions.

Physics of Fluids
Enforced reduction of the flux allows for a parameter space to be explored between Figs. 9(h) and 9(j) with the steady state solver for a defined flux rather than filling fraction. Beyond Fig. 9(j), increasing the filling fraction leads to accumulation of the extra fluid pooling on the rising cylinder wall exhibited in Fig. 9(k) comparable to the no-slip scenario with larger filling fractions [ Fig. 9(d)]; the observed shock structure experiences substantial short wavelength formations ahead of the shock front located on the rising side of the cylinder. This parameter space shows a region of defined filling fractions which lead to multiple solution profiles attainable, defined by the non-unique region between 1 < A/Ac < 1.2. The rimming flow profile presented in Fig. 9(i) was found to be stable; indeed, between solutions in Figs. 9(h) and 9(j), all solutions appeared to be stable, following analysis by the transient solver. This is in agreement with the theory on the location of shock solution stability, seen in the work of Villegas-Díaz et al., 24 where the fluid on the rising cylinder side is stable when on the lower cylinder section. By advent of slip, the maximum attainable flux ratio is increased, with a larger fluid filling fraction ratio required before reaching the transition region between the smooth and shock solutions. The onset of capillary ripples is more profound across the entire parameter range for slip, exhibited between Figs. 9(g) and 9(k), unlike the no-slip results. Both slip and no-slip cases experience a flux reduction, located in the transitional region from a smooth profile into a shock. The accumulation of fluid in the lower cylinder as a deepening shock solution provides the reduction in fluid flux required for maintaining of a smooth rimming flow. The non-unique region reported by Kay et al. 38 when implementing a wall friction factor affirms the use of surface effects at developing non-unique solution regions.

V. CONCLUSIONS
For the use of highly sheared thin liquid films in aero-engine components, specifically bearing chambers in which film thicknesses can be exceptionally thin, we have investigated the film flow dynamics where inertial and surface effects can be substantial. A model for rimming flow that retains leading-order inertial effects with a slip boundary condition imposed upon the chamber surface has been developed. The model consistency is discussed and compared to cases with negligible inertial effects or no-slip condition at the wall from previous studies of rimming flows.
Based upon simplification of the governing equations from hydraulics theory, a depth-averaged model is presented. An assumed quadratic velocity profile, recognised from the Stokes flow with a constant slip boundary condition, yields a physically relevant model for a range of slip conditions. Inertial terms are retained by depthaveraging the Navier-Stokes equations under a thin film assumption. The resulting system of equations is solved numerically, with a coordinate transformation to cluster points near the center. The computed film height and flux profiles are compared for no-slip and slip conditions at the wall by a compact numerical approach, with a global pseudospectral method in obtaining film height and flux profiles. The method to include high order effects such as leading order surface shear and inertial effects by depth-averaging is consistent with previous studies. 38 The addition of slip results in non-unique solution regions, observed through analysis of a parameter space for comparison between the no-slip condition of purely unique solutions. Existence of three simultaneously possible steady solutions is verified by a steady state solver, and by transient solver, all of them are determined to be wholly stable. Solution profiles with wall slip exhibit deeper pooling of films, arising from the thinning at the outer bounds of the domain, significant when the slip length is comparable to the film height.

ARTICLE scitation.org/journal/phf
Further analysis within a lubrication theory framework suggests modification of the well-known lubrication equation by wall slip. The asymptotic analysis of the resulting equation with an increasing slip length demonstrated reductions in allowable film volume alongside an escalation of the possible film flux, as a result of the modification of critical film values where hc asymptotically approaches τ/2λ and qc continuously increases. Upon fixation of the filling fraction, the addition of moderate slip lengths leads to a transition from a sub-critical to super-critical flow.
When inertial and gravitational effects become sufficiently large to dominate the flow dynamics, slip length has reduced impact on film profiles with results comparable to previous studies at low to no slip. The addition of slip leads to a swifter onset of capillary ripples whilst having a deferred transition toward smooth solutions for growing inertial effects. An introduction of boundary slip demonstrated a thinning of film profiles toward domain boundaries with a deeper pooling solution forming and introducing recirculation regions which are detrimental to bearing chamber flows.

ACKNOWLEDGMENTS
This work was supported by the Engineering and Physical Sciences Research Council (Grant No. EP/M507519/1) of the United Kingdom through the scholarship provided by the Industrial CASE Account for the University of Nottingham. The authors also thank Rolls-Royce PLC for its financial and technical support as part of the Industrial CASE Award.
We would like to acknowledge the commitment of Professor Henry Power, who was a great mentor and friend, who sadly passed away before this article was completed, and who will be dearly missed by all who had the pleasure of knowing him. The views expressed in this paper are those of the authors and not necessarily those of Rolls-Royce PLC.