Motivated by applications in aero-engines, steady two-dimensional thin-film flow on the inside of a circular cylinder is studied when the film surface is subject to mass and momentum transfer from impacting droplets. Asymptotic analysis is used systematically to identify distinguished limits that incorporate these transfer effects at leading order and to provide a new mathematical model. Applying both analytical and numerical approaches to the model, a set of stable steady, two-dimensional solutions that fit within the rational framework is determined. A number of these solutions feature steep fronts and associated recirculating pools, which are undesirable in an aeroengine since oil may be stripped away from the steep fronts when there is a core flow external to the film, and recirculation may lead to oil degradation. The model, however, provides a means of investigating whether the formation of the steep fronts on the film surface and of internal recirculation pools can be delayed, or inhibited altogether, by designing jets to deliver prescribed distributions of oil droplets or by the judicious siting of oil sinks. Moreover, by studying pathlines, oil-residence times can be predicted and systems optimized.

This study concerns rimming flow in a thin liquid film lining the inside of a horizontal circular cylinder forming the outer part of a cylindrical annulus. Rimming flow models are relevant to applications as diverse as cream separation, Rushak and Scriven,1 and the ceramic coating of pipes, Menekse et al.2 The particular motivation for the current study is an application to high-speed gas turbines, where rimming flows feature in the processes of the cooling and lubrication of aero-engine bearing chambers. In this application the motion of thin oil films is affected or effected, depending on whether the substrate is moving or not, by the impact of oil droplets. The cylindrical walls of the chamber comprise two moving shafts – the HP (high pressure) and IP (intermediate pressure) shafts, which may co-rotate or contra-rotate – and stationary sections. The fluid in the body of the chamber is a mixture of air and oil, driven round by the rotation of the central IP shaft: oil enters the chamber from an injector, forms a mist in the core air, and eventually impacts as small droplets on the film lining the outer wall. The oil then flows around the film and leaves the chamber through sinks located on the outer cylindrical wall and collects in external sumps. The function of the oil is to act as a lubricant for the bearings, and also as a coolant of the hot metal surfaces. It is important that a constant flow of oil is maintained around the wall, since stagnant or dry patches may lead to overheating and, in extreme cases, to coking.

It is the behaviour of the oil film that forms on the inside of the outer shaft of an aero-engine bearing chamber that is the primary concern in this study, which aims to develop a self-consistent model of droplet impact in rimming flow, and hence to derive and analyze expressions for film thickness. Thus the problem posed is to investigate the flow of a thin film of oil on the inside of a smooth rigid cylinder, which may be rotating or stationary, under the influences of gravity and a core flow consisting of a rotating mixture of air and oil droplets; thermal effects are not considered in the current study.

In an aero-engine, there will be a start-up phase where droplets impact on an initially dry surface and spread; as more droplets impact the surface will become coated with a continuous film and ultimately will reach a steady-state operating condition. The problem of the formation of films via droplet impact has received recent attention by Driscoll et al.3 In the “steady state,” droplets interact with the film in a complex manner, rebounding, splashing, and breaking up, as well as being absorbed, Farrall,4 Farrall et al.,5 Farrall et al.6 Such impact problems are extremely challenging, but theoretical progress has been made in the case of an isolated droplet, see, for example, Hicks and Purvis7 and references cited therein. Therefore, it is both reasonable and pragmatic to consider a continuous distribution obtained by averaging the droplet and air core flows spatially, rather than a discrete one. A similar approach was adopted previously in two limiting models of droplet impact at the surface of rimming film flow. In earlier investigations of draining and filling flows, Noakes,8 Noakes et al.9 ignored momentum, and regarded mass as being added continuously at the film surface and removed continuously at a sink. In contrast, Villegas-Díaz,10 Villegas-Díaz et al.,11 Villegas-Díaz et al.12 model the two-phase core flow as a fluid of low density which imparts momentum to the film via a constant surface shear along the film interface and droplet mass transfer is neglected. An adhoc modelling approach would simply be to combine the independent contributions of mass and momentum. This approach, however, is not adopted, instead continuity of mass and momentum is sought across the interface between the film and a prescribed core flow. The equation governing the film height h is determined by identifying distinguished asymptotic limits and the resulting equations’ relationship to previous models established.

It is well established that a complex variety of fluid-dynamical behaviour arises in rimming flow, see, for example, Thoroddsen and Mahadevan.13 There is correspondingly an ever-growing literature on the mathematical analysis of this problem: in particular, Moffatt,14 Pukhnachev,15 Johnson16 and O’Brien and Gath17 focus on a simple lubrication model; Ashmore et al.,18 Hosoi and Mahadevan,19 Noakes,8 Tirumkudulu and Acrivos,20 Wilson et al.,21 and Acrivos and Jin22 consider various higher order theories; some in addition investigate surface-tension effects. More recent studies, Benilov and O’Brien,23 Noakes et al.,24 and Kelmanson25 have focussed on inertial effects, while others, Hinch and Kelmanson,26 Hinch et al.,27 and Noakes et al.28 have used intricate multiple-timescale analysis to describe the evolution of instabilities in rimming flow. (Some of the papers cited here in turn cite many more papers, many of which are concerned with coating flow, the flow on the exterior of a cylinder. For our purposes this brief list of references is sufficient.)

In the simplest model of rimming flow a balance between viscosity and gravity is assumed and the effects of inertia and surface tension are neglected. Applying the ideas of lubrication theory to this situation, Moffatt14 studied the maximum load of liquid that can be supported on the outside of a cylinder (his analysis applies to both coating and rimming flows with appropriate modifications in interpretation). The problem of rimming flow was revisited by Johnson16 and by Benjamin, Pritchard and Tavener in unpublished work, with both studies addressing the question of the existence of weak solutions containing shocks (where the free-surface height changes rapidly in a very narrow region of the domain). It was Johnson16 who first recognized the relevance of discontinuous solutions of the equation determining the steady-state film profile as a function of the azimuthal angle. Benjamin et al. made a deeper study and, inter alia, considered the stability of the solutions that Johnson16 had identified and the influence of surface tension, the latter effect being investigated in more detail by Ashmore et al.18 A more recent, comprehensive study of the effects of surface tension and higher order gravitational effects has been undertaken by Villegas-Díaz et al.12 in a study of rimming flow subject to constant surface shear with the aim of resolving issues regarding existence, uniqueness, and stability of such weak solutions. (For a review of experimental work, see Evans et al.29)

There is a long-standing issue of correctly embedding shock structures in thin films within rational asymptotic schemes. Some progress has been made, see, for example, Noakes,8 Benilov et al.,30 and Benilov et al.,31 but the problem remains unresolved and is outside the scope of the current study. Finally, mention should be made of the use of CFD packages to study the current problem. In brief, while promising progress is being made on developing sophisticated packages to compute flow in aero-engine bearing chambers (see Lu et al.,32 for example), the authors know of no code that has satisfactorily resolved the interface between the liquid film and core two-phase flow and details of the film flow itself. Thus detailed analysis, such as that presented in this study, is currently the best means of gaining physical insight into this complex flow.

The structure of the paper is as follows. In Sec. II (and Appendix  A), the governing system of equations is described. In Sec. III (and Appendix  B), a systematic asymptotic analysis of the equations is developed and the key generic film evolution equation, (21), together with its relationship to previous studies, Table II, presented. In Sec. IV, the primary validation of the numerical code is undertaken. This numerical code together with theoretical analysis is applied in Sec. V to the case of a film on a stationary cylinder with significant droplet momentum but negligible mass transfer at the film surface, while in Sec. VI numerical solutions and their stability when droplet momentum is dominant are discussed. In Sec. VII, the case of significant mass and momentum transfer and the effects of varying mass flux, the film mass and sink position are considered. Finally in Sec. VIII, conclusions are summarized.

Table II.

Coefficients appearing in

$\frac{\partial h}{\partial t} + \frac{\partial }{\partial \theta }\left(c_1\left( \frac{\partial h}{\partial \theta }+\frac{\partial ^3 h}{\partial \theta ^3} \right)\frac{h^3}{3} + c_2\sin \theta \frac{h^3}{3} + \left(c_3\frac{\partial h}{\partial \theta }+c_4\right)\frac{h^2}{2} \right) + \left(c_5\frac{\partial h}{\partial \theta } +c_6\right)\break + \frac{\partial }{\partial \theta }\left(u_{cyl}h\right) - v_{cyl}= 0$
ht+θc1hθ+3hθ3h33+c2sinθh33+c3hθ+c4h22+c5hθ+c6+θucylhvcyl=0 and the corresponding coefficients featured in other studies; × not considered ✓ considered.

 c1c2c3c4c5/ucylc6/vcylHigher order terms
Leading-order cubic × ✓ × ✓ ✓ ×   
Momentum only (see Sec. IVζ −1 αu   
Mass and momentum RefWef−1 −RefFr−1 
$\mathrm{Re_f}a u_E^2$
Re fauE2
 
RefauEvE auE avE   
Noakes et al. (Ref. 9) (rotation and mass) ε −1 λI   
Villegas-Díaz (Ref. 10) (rotation and shear) B−1 −Γ γ   
Moffatt (Ref. 14× ✓ × ✓ ✓ × × 
Black (Ref. 34) Smoothing by               
surface tension ✓ ✓ × ✓ × × × 
Black (Ref. 34) Smoothing by gravity × ✓ × ✓ × × ✓ 
Black (Ref. 34) Rotating × ✓ × ✓ ✓ × × 
Benjamin et al. (unpublished) ✓ ✓ × ✓ ✓ × ✓ 
Acrivos and Jin (Ref. 22),               
Tirumkudulu and Acrivos (Ref. 20✓ ✓ × × ✓ × ✓ 
Hinch and Kelmanson (Ref. 26),               
Hinch et al. (Ref. 27✓ ✓ × × ✓ × × 
 c1c2c3c4c5/ucylc6/vcylHigher order terms
Leading-order cubic × ✓ × ✓ ✓ ×   
Momentum only (see Sec. IVζ −1 αu   
Mass and momentum RefWef−1 −RefFr−1 
$\mathrm{Re_f}a u_E^2$
Re fauE2
 
RefauEvE auE avE   
Noakes et al. (Ref. 9) (rotation and mass) ε −1 λI   
Villegas-Díaz (Ref. 10) (rotation and shear) B−1 −Γ γ   
Moffatt (Ref. 14× ✓ × ✓ ✓ × × 
Black (Ref. 34) Smoothing by               
surface tension ✓ ✓ × ✓ × × × 
Black (Ref. 34) Smoothing by gravity × ✓ × ✓ × × ✓ 
Black (Ref. 34) Rotating × ✓ × ✓ ✓ × × 
Benjamin et al. (unpublished) ✓ ✓ × ✓ ✓ × ✓ 
Acrivos and Jin (Ref. 22),               
Tirumkudulu and Acrivos (Ref. 20✓ ✓ × × ✓ × ✓ 
Hinch and Kelmanson (Ref. 26),               
Hinch et al. (Ref. 27✓ ✓ × × ✓ × × 

In the idealized model, the outer shaft of the aero-engine bearing chamber is represented as a cylinder of radius r0 with a thin film of liquid on its inner surface. Define cylindrical polar coordinates

$\hat{r}$
r̂⁠, with origin at the centre of the horizontal cylinder, and θ, measured anti-clockwise from the downward vertical, see Figure 1(a). The cylinder wall is described by
$\hat{r}={r}_0$
r̂=r0
, and the film surface by
$\hat{r}=r_0-\hat{h}(\theta ,\hat{t})$
r̂=r0ĥ(θ,t̂)
where
$\hat{h} \ge 0$
ĥ0
.

The effects of liquid droplets on the film are modelled as contributions to the surface forces and to the mass flux across the film surface. The droplets are assumed to be of the same constant density ρ and viscosity μ as the liquid film. The flow is assumed two-dimensional and the fluid incompressible and Newtonian.

In this study, we are not concerned with the detail of the two-phase core flow, so the behaviour “far” (on the scale of the film thickness) from the cylinder wall is not included in the model. For simplicity, the two-phase core fluid is modelled as an inviscid mixture (in an average sense) and the flow is assumed to have a constant local pressure

$\hat{p_E}$
pÊ outside the film and an average droplet velocity near the film surface
$\mathbf {\hat{u}_E}$
ûE
at a constant angle relative to the cylinder wall and at a constant speed such that the film may be driven against gravity; no rebounding, splashing, or breaking up of droplets is considered. Consistent with the core-flow averaging, mass, and momentum transfer are considered across a smoothed film surface. Subscripts E and I are used to refer to the two-phase gas-droplet mix external to the film and the liquid in the film, respectively; when there is no ambiguity, the subscripts will be dropped in what follows.

In the film, the fluid has kinematic viscosity ν = μ/ρ and the flow is subject to a gravitational field g = −gj, where j is a unit vector pointing vertically upwards and g is the acceleration due to gravity, thus the velocity,

${\mathbf {\hat{u}_I}}=(\hat{u}_{Ir},\hat{u}_{I\theta })$
ûI=(ûIr,ûIθ)⁠, and pressure,
$\hat{p}_I$
p̂I
, satisfy the continuity and Navier–Stokes equations

\begin{eqnarray}\hspace*{10.8pc} {\mathbf \nabla } \cdot {\mathbf {\hat{u}_I}} = 0,\end{eqnarray}
·ûI=0,
(1)
\begin{eqnarray}\frac{\partial {\mathbf {\hat{u}_I}}}{\partial t} + {\mathbf {\hat{u}_I} } \cdot {\mathbf {\nabla }} {\mathbf {\hat{u}_I} } = -\rho ^{-1} {\mathbf {\nabla } }{\hat{p}_I} + \nu \nabla ^2 {\mathbf {\hat{u}_I}} + {\mathbf g};\end{eqnarray}
ûIt+ûI·ûI=ρ1p̂I+ν2ûI+g;
(2)

2π-periodic solutions will be sought to (1) and (2).

The core two-phase flow comprises air with density ρg and a volume fraction α of liquid droplets with density ρ (and ρI ≡ ρ). ρE = αρ + (1 − α)ρg is the average density of the mixture; given that (1 − α)ρg ≪ αρ, the two-phase fluid is modelled as continuous and of constant density ρE = αρ.

Droplets impacting on the film surface, which is modelled as a single interface S(t), described by χ(x, t) = 0, transfer mass as well as momentum into the film. This interface is a surface at which density and velocity are discontinuous (see Appendix  A for further details).

The velocity of the interface

$\hat{U}_j$
Ûj follows from conservation of mass (cf. Sec. 1 in Appendix  A):

\begin{equation}\rho (\alpha -1) \hat{U}_j n_j=\rho (\alpha \hat{u}_{Ej}-\hat{u}_{Ij}) n_j.\end{equation}
ρ(α1)Ûjnj=ρ(αûEjûIj)nj.
(3)

The normal and tangential components of the stress boundary conditions on the interface are obtained from the conservation of momentum (see Sec. 2 in Appendix  A) as

\begin{eqnarray}\hspace*{-.8pc} -\sigma \kappa =\displaystyle { \frac{\alpha \rho }{1-\alpha }[\hat{u}_{i}][\hat{u}_j]n_j n_i - [\hat{T}_{ij}]n_jn_i},\end{eqnarray}
σκ=αρ1α[ûi][ûj]njni[T̂ij]njni,
(4)
\begin{eqnarray}0=\displaystyle {\frac{\alpha \rho }{1-\alpha }[\hat{u}_{i}][\hat{u}_j]n_jt_i - [\hat{T}_{ij}]n_jt_i};\end{eqnarray}
0=αρ1α[ûi][ûj]njti[T̂ij]njti;
(5)

[.] denotes the jump of quantities across the boundary between the two-phase flow and the film, so, for example, [u] = uEuI; σ is an effective surface tension coefficient, ti is the unit tangent vector at the interface such that niti = 0, and

$\hat{T}_{ij}$
T̂ij are the fluids’ stress tensors.

In the current case, the film surface is

$\chi \equiv -\hat{r} + r_0 -\hat{h}(\theta ,\hat{t})=0$
χr̂+r0ĥ(θ,t̂)=0⁠. Correspondingly, n = (nr, nθ), the unit normal pointing towards the centre of the cylinder has components

\begin{eqnarray}n_r=-N ,\quad n_\theta = -\frac{N}{\hat{r}}\frac{\partial \hat{h}}{\partial \theta },\end{eqnarray}
nr=N,nθ=Nr̂ĥθ,
(6)

where

\begin{eqnarray}N=\left(1+\left(\frac{1}{\hat{r}}\frac{\partial \hat{h}}{\partial \theta }\right)^2\right)^{-1/2},\end{eqnarray}
N=1+1r̂ĥθ21/2,
(7)

and the surface curvature is given by

$\kappa =\nabla \cdot \mathbf {n}|_{\hat{r}=r_0-\hat{h}}$
κ=·n|r̂=r0ĥ⁠, i.e.,

\begin{eqnarray}\kappa =\nabla \cdot \mathbf {n} = -\frac{N^3}{\hat{r}} \left(1+\frac{2}{\hat{r}^2}\left(\frac{\partial \hat{h}}{\partial \theta }\right)^2 +\frac{1}{\hat{r}}\frac{\partial ^2\hat{h}}{\partial \theta ^2} \right).\end{eqnarray}
κ=·n=N3r̂1+2r̂2ĥθ2+1r̂2ĥθ2.
(8)

The final expressions of the stress matching conditions (A19) and (A20) are given in Sec. 2 in Appendix  A.

The interface position χ = 0 is a material surface, hence

\begin{eqnarray}\frac{\partial \chi }{\partial t}+\hat{U}_j\frac{\partial \chi }{\partial x_j}=0,\end{eqnarray}
χt+Ûjχxj=0,
(9)

which, using (3), yields the kinematic condition

\begin{eqnarray}\left(1-\alpha \right)\frac{\partial \hat{h}}{\partial \hat{t}} +\left(\hat{u}_{I\theta }\frac{1}{\hat{r}}\frac{\partial \hat{h}}{\partial \theta }-\hat{u}_{Ir}\right) =\alpha \left(\hat{u}_{ E\theta }\frac{1}{\hat{r}}\frac{\partial \hat{h}}{\partial \theta }-\hat{u}_{Er}\right),\quad \mathrm{on}\quad \hat{r}=r_0-\hat{h}.\end{eqnarray}
1αĥt̂+ûIθ1r̂ĥθûIr=αûEθ1r̂ĥθûEr, on r̂=r0ĥ.
(10)

Finally, mass leaves a bearing chamber through one or more outlets, which are modelled as sinks distributed on the chamber wall. The sinks are described mathematically by specifying the normal velocity at the wall

\begin{eqnarray}\hat{u}_{Ir}=\hat{v}_{cyl}.\end{eqnarray}
ûIr=v̂cyl.
(11)

Here,

$\hat{v}_{cyl}$
v̂cyl is uniform across the sinks such that at steady state the total influx at the film surface is balanced by the total outflux at the sinks, while
$\hat{v}_{cyl}=0$
v̂cyl=0
elsewhere on the cylinder. For simplicity, the azimuthal velocity is taken to match the motion of the cylinder wall everywhere around the boundary (i.e., sinks are ignored), so that
$\hat{u}_{I\theta } = \hat{u}_{cyl}$
ûIθ=ûcyl
where
$\hat{u}_{cyl} = \Omega r_0$
ûcyl=Ωr0
for a cylinder rotating anticlockwise with angular velocity Ω; for a film driven purely by the impacting droplets’ momentum, the cylinder is stationary Ω = 0.

The study concentrates on thin film flows for which a typical film thickness such that

$\hat{h}_0 \ll {r}_0$
ĥ0r0 can conveniently be defined by

\begin{eqnarray}\hat{h}_0=\frac{1}{2 \pi } \int _0^{2\pi } \hat{h}(\theta ,0) \mbox{d}\theta .\end{eqnarray}
ĥ0=12π02πĥ(θ,0)dθ.
(12)

It is natural to employ a local coordinate system using a non-dimensional radial coordinate, y measuring distance inwards from the cylinder's surface, see Figure 1(ii). The scalings for the film variables are quite standard in rimming-flow analysis (U is a characteristic film speed); the scalings of the impact velocity have been chosen to ensure a balance between the radial flux of azimuthal droplet momentum

$\sim \hat{u}_{E\theta }(\alpha {\rho } \hat{u}_{Er})$
ûEθ(αρûEr) and the surface viscous shear
$\sim \mu \partial \hat{u}_{I\theta }/\partial \hat{r}$
μûIθ/r̂
,

\begin{eqnarray}\hat{t} =r_0 t/U,\quad \hat{r} = r_0(1-\epsilon y), \quad \hat{h} = \hat{h}_0 h, \quad \alpha = \epsilon ^k a,\quad \epsilon =\hat{h}_0/r_0, \nonumber \\[3pt](\hat{p}_{I},\hat{p}_{E}) = \rho U^2 \epsilon ^{-1}(p,p_E-\epsilon ^{-1}\mathrm{We_f}^{-1}),\quad (\hat{u}_{Ir},\hat{u}_{I\theta }) = U (-\epsilon v,u), \nonumber \\[3pt](\hat{u}_{Er},\hat{u}_{E\theta }) = U \epsilon ^{-(k+\beta )/2}(-\epsilon ^\beta v_E,u_E), \quad \quad (\hat{v}_{cyl},\hat{u}_{cyl}) = U(v_{cyl},u_{cyl}); \\[3pt]\mathrm{We_f}= U^2\rho r_0^3 / \sigma \hat{h}_0^2,\quad \mathrm{Re_f}= U\hat{h}_0/\nu , \quad \mathrm{Fr}= U^2/g\hat{h}_0, \nonumber\end{eqnarray}
t̂=r0t/U,r̂=r0(1εy),ĥ=ĥ0h,α=εka,ε=ĥ0/r0,(p̂I,p̂E)=ρU2ε1(p,pEε1 We f1),(ûIr,ûIθ)=U(εv,u),(ûEr,ûEθ)=Uε(k+β)/2(εβvE,uE),(v̂cyl,ûcyl)=U(vcyl,ucyl); We f=U2ρr03/σĥ02, Re f=Uĥ0/ν, Fr =U2/gĥ0,
(13)

where it is assumed that a = O(1), Ref = O(1), Wef = O(1), and Fr = O(1) as ε → 0. The timescale for motion in the film is that of a particle of fluid moving around the cylinder.

FIG. 1.

Film variables in cylindrical and scaled coordinates.

FIG. 1.

Film variables in cylindrical and scaled coordinates.

Close modal

Note that k controls the strength of α, the volume fraction of droplets in the two-phase flow; β controls the impact angle of the droplets since

$\hat{u}_{Er}/\hat{u}_{E\theta } = O(\epsilon ^\beta )$
ûEr/ûEθ=O(εβ)⁠; the pressure in the two-phase flow is measured relative to (εWef)−1 for later convenience; Wef, Ref, Fr, and ε denote reduced Weber number, reduced Reynolds number, Froude number, and aspect ratio, respectively. It should be emphasized that Wef may in practice be quite large, but surface tension provides an important smoothing mechanism and is thus retained in the analysis, cf. Benilov et al.33 

The dimensional film area is

\begin{eqnarray}\hat{A} = \int ^{\pi }_{-\pi }\int ^{r_0}_{r_0 - \hat{h}} \hat{r} \mathrm{d}\hat{r} \mathrm{d}\theta = \hat{h}_0 r_0 \int ^{\pi }_{-\pi } h + \frac{\epsilon h^2}{2} \mathrm{d}\theta ,\end{eqnarray}
Â=ππr0ĥr0r̂dr̂dθ=ĥ0r0ππh+εh22dθ,
(14)

which naturally leads to a leading-order non-dimensional measure

\begin{eqnarray}A= \int _{-\pi }^{\pi } h \mathrm{d}\theta .\end{eqnarray}
A=ππhdθ.
(15)

Note that the mass of fluid in the film per unit axial length is given simply by

$M = \rho \hat{A}$
M=ρÂ⁠.

The details of the derivation of the asymptotic form of the film-evolution equation are relegated to Appendix  B. To leading order, the dimensionless profile h is

\begin{eqnarray}\frac{\partial h_0}{\partial t} -\frac{\partial }{\partial \theta }\left(\mathrm{Re_f}\left(\frac{\partial \mathcal {N}_{0}}{\partial \theta }+ {\mathrm{Fr^{-1}}}\sin \theta \right)\frac{h_0^3}{3}\right) +\frac{\partial }{\partial \theta }\left(\mathcal {T}_0\frac{h_0^2}{2}\right)+\frac{\partial }{\partial \theta }\left(u_{cyl}h_0\right)-v_{cyl}-\mathcal {M}_0=0, \nonumber \\\end{eqnarray}
h0tθ Re fN0θ+ Fr 1sinθh033+θT0h022+θucylh0vcylM0=0,
(16)

where the effects of droplet impact are embedded in

$\mathcal {N}_{0}$
N0⁠,
$\mathcal {T}_0$
T0
, and
$\mathcal {M}_0$
M0
, the leading-order terms in the small-ε expansions of
$\mathcal {N}$
N
,
$\mathcal {T}$
T
, and
$\mathcal {M}$
M
. Analysis of (B7)–(B9) reveals that contributions arise from these terms at O(1) as ε → 0 provided that

\begin{eqnarray}|\beta| \le 1 \quad k \ge 2-\beta\end{eqnarray}
|β|1k2β
(17)

in which case

\begin{eqnarray}\hspace*{-4.5pc} \mathcal {N} \sim -\mathrm{We_f}^{-1}\left(h_0+\frac{\partial ^2 h_0}{\partial \theta ^2}\right)+a\epsilon ^{1+\beta } v_E^2 + \ldots ,\end{eqnarray}
N We f1h0+2h0θ2+aε1+βvE2+...,
(18)
\begin{eqnarray}\mathcal {T} \sim &-&\mathrm{Re_f}a\left(u_Ev_E - \epsilon ^{1-\beta } u_E^2\frac{\partial h_0}{\partial \theta } + \epsilon ^{1+\beta } v_E^2\frac{\partial h_0}{\partial \theta } - \epsilon ^\frac{k+\beta }{2}uv_E \right. \nonumber \\&&\quad \qquad + \left.\epsilon ^\frac{k-\beta +2}{2}\left(2uu_E \frac{\partial h_0}{\partial \theta }-vu_E\right)+ \ldots \right),\end{eqnarray}
T Re fauEvEε1βuE2h0θ+ε1+βvE2h0θεk+β2uvE+εkβ+222uuEh0θvuE+...,
(19)
\begin{eqnarray}\hspace*{-4pc} \mathcal {M} \sim a\left(\epsilon ^k \frac{\partial {h_0}}{\partial {t}}+ \epsilon ^\frac{k-\beta }{2}u_E\frac{\partial h_0}{\partial \theta } -\epsilon ^\frac{k+\beta -2}{2}v_E + \ldots \right).\end{eqnarray}
Maεkh0t+εkβ2uEh0θεk+β22vE+....
(20)

Note that the first three terms on the right-hand side of (19) arise from the flux across the surface of the droplets’ momentum in the direction tangential to that surface, namely, a(vEuEh/∂θ)(vEh/∂θ + uE). The exact form of the evolution equation (16) depends on the droplet fraction in the core flow (k) and droplets’ inclination angle (β). In view of (17), it follows that k ⩾ 1 and hence the model is restricted to dilute droplet volume fractions (consistent with the assumption that the two-phase flow may be treated as inviscid).

Droplets approach the film at a prescribed angle ϕ to the local tangent to the cylindrical wall, where tan ϕ = εβvE/uE. Physically, as illustrated in Figure 2, β = 1, 0, − 1 correspond to impacts at shallow, moderate, and steep angles, respectively.

FIG. 2.

Droplet impact angle is characterized by β.

FIG. 2.

Droplet impact angle is characterized by β.

Close modal

Substituting (18)–(20) into (16), and omitting straightforward details, it is possible to write the generic film evolution equation corresponding to all possible cases that can exist on the parameter space defined in (17) as

\begin{eqnarray}\frac{\partial h}{\partial t} +\frac{\partial }{\partial \theta }\left( \overbrace{c_1\left( \frac{\partial h}{\partial \theta }+\frac{\partial ^3 h}{\partial \theta ^3} \right)\frac{h^3}{3} }^{\mathit {\text{Surface} \,\text{Tension}}} + \overbrace{c_2\sin \theta \frac{h^3}{3}}^{\mathit {\text{Gravity}}} + \overbrace{\left(c_3(\theta )\frac{\partial h}{\partial \theta }+c_4(\theta )\right)\frac{h^2}{2}}^{\mathit {\text{Droplet}\,\text{Momentum}}} \right)\nonumber \\+\underbrace{\left(c_5(\theta )\frac{\partial h}{\partial \theta } +c_6(\theta )\right)}_{\mathit {\text{Droplet} \,\text{Mass}}} + \underbrace{\frac{\partial }{\partial \theta }\left(u_{cyl}h\right)}_{\mathit {\text{Wall} \, \text{Rotation}}} - \underbrace{v_{cyl}}_{\mathit {\text{Sink}}} =0.\end{eqnarray}
ht+θc1hθ+3hθ3h33Surface Tension+c2sinθh33Gravity+c3(θ)hθ+c4(θ)h22Droplet Momentum+c5(θ)hθ+c6(θ)Droplet Mass+θucylhWall RotationvcylSink=0.
(21)

The coefficients appearing in (21) are given in Table I for a selection of admissible values of β and k satisfying (17); note that β and k need not be integer, but that the selected integer values illustrate the full range of interesting qualitative behaviours. At a shallow inclination (β = 1), momentum is dominant for lower droplet fraction (k = 3), as c5 and c6 are negligible, while for higher droplet fraction (k = 1) both momentum and mass transfer are significant, and so on.

Table I.

Coefficients appearing in (21) for typical admissible values of β and k.

kβc1c2c3c4c5c6
−1 RefWef−1 −RefFr−1 
$-\mathrm{Re_f}a v_E^2$
Re favE2
 
−RefauEvE avE 
RefWef−1 −RefFr−1 −RefauEvE 
RefWef−1 −RefFr−1 
$\mathrm{Re_f}au_E^2$
Re fauE2
 
−RefauEvE 
RefWef−1 −RefFr−1 −RefauEvE avE 
RefWef−1 −RefFr−1 
$\mathrm{Re_f}a u_E^2$
Re fauE2
 
−RefauEvE 
RefWef−1 −RefFr−1 
$\mathrm{Re_f}a u_E^2$
Re fauE2
 
−RefauEvE auE avE 
kβc1c2c3c4c5c6
−1 RefWef−1 −RefFr−1 
$-\mathrm{Re_f}a v_E^2$
Re favE2
 
−RefauEvE avE 
RefWef−1 −RefFr−1 −RefauEvE 
RefWef−1 −RefFr−1 
$\mathrm{Re_f}au_E^2$
Re fauE2
 
−RefauEvE 
RefWef−1 −RefFr−1 −RefauEvE avE 
RefWef−1 −RefFr−1 
$\mathrm{Re_f}a u_E^2$
Re fauE2
 
−RefauEvE 
RefWef−1 −RefFr−1 
$\mathrm{Re_f}a u_E^2$
Re fauE2
 
−RefauEvE auE avE 

In Eq. (21), the term proportional to c5 corresponds to the azimuthal component of the mass flux and c6 to the radial component. In contrast, the whole term involving c3 and c4 is a contribution from droplet momentum, but the individual terms proportional to c3 and c4 do not lend themselves to a simple physical interpretation.

The nature of the constant-shear model of Villegas-Díaz et al.12 (hereafter referenced simply as VPR) in the context of the current generalized model is clear from (17)–(20): it is for ∣β∣ < 1 and k > 2 − β that droplets impart momentum equivalent to the effect of a constant surface shear force (i.e., the c4 term alone). If ∣β∣ = 1 and k > 2 − β, the c3 term in the droplet momentum is significant but is absent in the model of VPR. Finally, we note that if −1 ⩽ β < 1 and k = 2 − β, a constant mass flux term (i.e., the c6 term) of the form modelled in Noakes et al.9 features, but, in contrast to Noakes's model,8 momentum is not negligible. Table II summarizes the terms that arise in previous studies that lead to an equation of similar form to (21).

In Secs. IV–VI, we consider cases where the flow field is driven solely by the transfer of droplet momentum, with a stationary cylinder wall, i.e., ucyl = 0, and the mass flux into the film is negligibly small (c5 = c6 = 0). This implies that ∣β∣ = 1 and k > 2 − β; specifically we consider β = 1 and k = 2 and 3. Later in Sec. VII, we study the case of significant droplet momentum transfer and mass flux at the film surface, i.e., ∣β∣ ⩽ 1 and k = 2 − β, in order to complete the exploration of the whole of the parameter space defined in (17).

To simplify the governing equation in the case of negligible mass transfer

\begin{eqnarray}\frac{\partial h}{\partial t} + \mathrm{Re_f}\frac{\partial }{\partial \theta }\left( {\mathrm{We_f}^{-1}\left( \frac{\partial h}{\partial \theta }+\frac{\partial ^3 h}{\partial \theta ^3} \right)\frac{h^3}{3} } +{\mathrm{Fr}^{-1}} \sin \theta \frac{h^3}{3} + a u_E {\left(u_E\frac{\partial h}{\partial \theta }-v_E\right)\frac{h^2}{2}} \right) =0,\end{eqnarray}
ht+ Re fθ We f1hθ+3hθ3h33+ Fr 1sinθh33+auEuEhθvEh22=0,
(22)

in cases in parameter space where ∣β∣ = 1 and k > 2 − β, i.e., dominant droplet momentum with c5 = c6 = 0, we define the characteristic film thickness as

$\hat{h}_{0}=(U \nu / g)^\frac{1}{2}$
ĥ0=(Uν/g)12 so that RefFr−1 = 1, and the characteristic velocity as
$U=(\nu g/a^{2} u_{E}^{2} v_{E}^{2})^{1/3}$
U=(νg/a2uE2vE2)1/3
. The final term in (22) which represents azimuthal momentum becomes

\begin{eqnarray}\mathrm{Re_f}a u_E \left(u_E\frac{\partial h}{\partial \theta } -v_E\right)=\mathrm{Re_f}a u_E &#x007c;v_E&#x007c; \left(\frac{u_E}{&#x007c;v_E&#x007c;}\frac{\partial h}{\partial \theta }+1 \right)= \mathrm{sgn}(u_E) \left(\alpha _u\frac{\partial h}{\partial \theta }+1 \right)\end{eqnarray}
Re fauEuEhθvE= Re fauE|vE|uE|vE|hθ+1= sgn (uE)αuhθ+1
(23)

with sgn(uE) = uE/|uE| and αu = uE/|vE|; note that vE is always negative. On introducing ζ = RefWef−1 = σε3U, (22) becomes

\begin{eqnarray}\frac{\partial h }{\partial t }+\frac{\partial q }{\partial \theta }=0, \quad q =\zeta \frac{h ^3}{3}\left( \frac{\partial h }{\partial \theta }+\frac{\partial ^3 h }{\partial \theta ^3}\right)-\sin \theta \frac{h ^3}{3} +\mathrm{sgn}(u_E) \left(\alpha _u\frac{\partial h }{\partial \theta }+1\right)\frac{h ^2}{2}.\end{eqnarray}
ht+qθ=0,q=ζh33hθ+3hθ3sinθh33+ sgn (uE)αuhθ+1h22.
(24)

To summarize, on choosing

$\hat{h}_0$
ĥ0 and U as above, the parameter values in (21) are c1 ≡ ζ, c2 ≡ −1, c3 ≡ αu, and c4 ≡ 1.

Steady-state solutions are first sought for a dominant balance between surface shear and gravity (ζ, αu = 0) as the problem is purely algebraic and an explicit solution exists. The solutions for small non-zero ζ, αu are expected to be close to these, except where derivatives of h are large, causing the terms involving ζ, αu to have a leading-order effect.

For definiteness, we take the azimuthal velocity to be positive, i.e., sgn(uE) = +1, and so

\begin{equation}q=-\sin \theta \frac{h^3}{3}+\frac{h^2}{2}.\end{equation}
q=sinθh33+h22.
(25)

Black34 was the first to present in-depth analysis of this particular flow. When 0 < q < 1/6, there is a single, bounded, neutrally stable film profile which is smooth and completely wets the cylinder's surface. This profile has vertical symmetry, its thickness being a maximum at θ = π/2, i.e., at the equator on the “rising side” of the flow, and a minimum at the point diagonally opposite. There are two other profiles: a C-profile which is totally negative and physically infeasible, and an unbounded profile consisting of two disconnected positive and negative branches. When A = Ac = 3.9509 and q = 1/6, however, the profile remains symmetric and continuous, but loses smoothness, developing a corner with a discontinuous slope at θ = π/2. At this point the positive unbounded branch effectively coalesces with the physical branch. For q = 1/6 and Ac < A, no continuous profile exists, but asymmetric (weak) solutions can be constructed from the two solution branches by introducing a shock on the side of the cylinder where the fluid flows upwards (i.e., by jumping between the two branches).

To leading order, with q = 1/6 and a shock at θ = θK, we have

\begin{equation}A=\frac{3}{2}\ln \left(\frac{2\cos \left(\theta _K/3-\pi /6\right)+\sqrt{3}}{2\cos \left(\theta _K/3-\pi /6\right)-\sqrt{3}} \right);\end{equation}
A=32ln2cosθK/3π/6+32cosθK/3π/63;
(26)

when θK → 0, A → ∞ according to

\begin{equation}A=\frac{3}{2}\ln \left(\frac{6\sqrt{3}}{\theta _K}\right)+\frac{\theta _K}{\sqrt{3}}+O(\theta _K^2).\end{equation}
A=32ln63θK+θK3+O(θK2).
(27)

As highlighted by Black,34 the result that A → ∞ as the shock position approaches the lowest point of the cylinder (θ = 0) contrasts with that for the problem of flow around a rotating cylinder with zero surface stress when A asymptotes to a finite value.

When q > 1/6, there is no completely wetting solution, but so-called curtain solutions which become unbounded at the north and south poles of the cylinder may exist – these solutions are outside the scope of the present study. (Such solutions are briefly discussed in Black,34 who also gives explicit expressions for the real solutions of the cubic equation.)

In the case of negative azimuthal velocity, sgn(uE) = −1, the governing equation is

\begin{equation}q=-\sin \theta \frac{h^3}{3}-\frac{h^2}{2},\end{equation}
q=sinθh33h22,
(28)

solutions to which are obtained by taking the mirror image of the solutions (i.e., θ → −θ) to (25) with q < 0. Without loss of generality, the analysis hereafter is restricted to uE > 0.

Of key physical interest is the study of bounded, completely-wetting, steady-state film profiles. There are many such profiles, including, in principle, infinitely many featuring multiple shocks. If shocks are present, then for q < 1/6 there must be an even number of them. For q = 1/6 there must be an odd number. Whether such profiles are physically realisable, however, will depend on their stability.

In a physical system, surface tension and the droplet momentum locally smooth the shock profiles. Even when 0 < ζ, αu ≪ 1, the terms αuh/∂θ and ζ∂3h/∂θ3 in (24) can be significant in the region of a shock and can determine whether such shock solutions are feasible in the context of a mathematically rational asymptotic structure. By looking at inner steady-state solutions of (24) in a region around a shock at θ = θK, and considering the possibility of matching them asymptotically with the outer solutions defined by the cubic polynomial (25), which in 0 < θ < π has two positive real roots h = hj say, where j = 1, 2 on the upper and lower branches, respectively, and a negative root h3, it is possible to analyze which shock structures are feasible (for more details about this type of analysis see Black,34 Villegas-Díaz et al.12 and Benilov et al.33).

Omitting all the detail, the results on smoothing can be summarized as follows. When A exceeds Ac = 3.9509, then for

  • αu = O1/3) as ζ → 0 (surface tension and droplet momentum both significant): smoothed backward- (forward-) facing shocks located in 0 < θ < π are feasible – the inner solution has an oscillatory (monotonic) match with the outer. This also holds if surface tension is dominant (i.e.,

    $\alpha _u^3 \ll \zeta \ll 1$
    αu3ζ1 as ζ → 0);

  • $\zeta \ll \alpha _u^3 \ll 1$
    ζαu31 as αu → 0 (droplet momentum dominant): a smoothed forward-facing shock located in π/2 < θ < π is feasible – the inner solution has a monotonic match with the outer.

Some possible shock combinations for q = 1/6 are sketched in Figure 3. Not all the structures highlighted in Figure 3 as being feasible according to the local asymptotic matching analysis are indeed realisable.

FIG. 3.

Profiles of h against θ illustrating various possible shock combinations when surface tension and droplet momentum are both insignificant. (a) stable; (b)–(d) unstable; (e) and (f) stable according to kinematic wave theory, but not seen in numerical solutions (see later). The arrows indicate the direction of travel of infinitesimal perturbations.

FIG. 3.

Profiles of h against θ illustrating various possible shock combinations when surface tension and droplet momentum are both insignificant. (a) stable; (b)–(d) unstable; (e) and (f) stable according to kinematic wave theory, but not seen in numerical solutions (see later). The arrows indicate the direction of travel of infinitesimal perturbations.

Close modal

By looking at the evolution of an infinitesimally small perturbation ξ(θ, t) to a steady state solution hs of (24), cf. (32) with ζ, αu = 0, as suggested by Moffatt14 and by Benjamin et al. (see also Villegas-Díaz et al.12 and Benilov et al.33), it is possible to show that three of the solutions reported in Figure 3 are unstable. The motion of a small perturbation is described by the kinematic wave equation with wave velocity w given by

\begin{eqnarray}w= (\mathrm{d}q/\mathrm{d}h)_{h=h_s(\theta )} = - h_s^2 \sin \theta + h_s = h_s^3 \cos \theta /3 (\mathrm{d}h_s/\mathrm{d}\theta ),\end{eqnarray}
w=(dq/dh)h=hs(θ)=hs2sinθ+hs=hs3cosθ/3(dhs/dθ),
(29)

provided the derivative is non-zero in the final identity. Now w < 0 on the upper solution branch of cubic equation (25) and w > 0 on the lower, hence disturbances travel towards a backward-facing shock and away from a forward-facing shock.

Significant progress can be made directly by noting that waves are unable to travel through a point where w = 0 and, if waves travel towards that point, energy is focussed and the amplitude of the disturbance grows unboundedly. Thus, in terms of linear theory, such unbounded growth might lead to interface rupture or blow up and the interface is unstable. Furthermore, if waves travel towards a shock, then that shock can adjust slightly and absorb the disturbance, leaving a stable structure.

When q = 1/6, w = 0 at θ = π/2, and waves may focus at this critical point. Thus profiles with a single backward-facing shock in the lower quadrant (0 < θ < π/2) are stable, profiles with a single forward-facing shock in the upper quadrant (π/2 < θ < π) unstable, and profiles featuring multiple shock solutions may be stable (unstable) if the profile height is decreasing (increasing) as θ increases through π/2. The directions of travel of infinitesimally small disturbances are superimposed on the profiles in Figure 3. Cases (b)–(d) are unstable because of the focusing effect at θ = π/2, and will not be seen in practice; (a) is stable; while kinematic wave theory does not rule out the other two cases, (e) and (f), later numerical studies reveal that they are unstable and evolve to structures similar to that in case (a).

To study the shock structures described above in more depth, steady-state solutions to (24) were found using a transient solver with values of A, αu, and ζ for which shock structures are expected to exist. In the present study, the results were determined by solving the film-evolution equation numerically using a fully implicit backward-time finite-difference scheme, with the resulting fourth-order differential equation in space, subject to a periodicity condition at θ = π and −π, being solved by the MATLAB boundary-value routine bvp4c.

In principle, the initial film profile may be arbitrarily specified. An initial profile having uniform height is the easiest way to impose a specific value of A, but a “Moffatt” solution (i.e., one arising from the algebraic equation) with a shock imposed near the expected location was occasionally used. Continuation from a smooth solution found for a neighbouring parameter set was also employed to reduce the time to compute the steady-state solution from that initiated from a uniform or from a shock solution. The condition that A is conserved through time was checked using the trapezoidal rule.

A timestep of 0.1 provided sufficient accuracy in most cases. The simulation was terminated and a steady solution assumed when the maximum change in h between one timestep and the next was less than 0.001.

To validate the numerical scheme, the cases of constant surface traction in a rotating cylinder reported by VPR, i.e., when c3, c5, c6, and vcyl in (21) are set to zero, were calculated and excellent agreement with VPR obtained. Having validated the scheme, the effect of the parameter c3 in the droplet momentum neglected by VPR was investigated for the particular parameter set c1 = 10−3, c4 = −2.149, and A = 3.739. The results obtained show very little change when c3 increases from 0 to 10−3|c4|, but when c3 is increased to 0.215 = 10−2|c4| the film profile peak sharpens, and for larger c3 a steady-state film profile cannot be found. These results suggest that the effect of the c3 contribution in the droplet momentum on the evolution of the film profile is destabilizing. The analysis of this destabilizing effect is one of the significant contributions of the current study and will be considered in more detail in Sec. VI.

To study further the effects of the droplet momentum transfer and surface tension, i.e., the effect of the terms proportional to αu and ζ in (24), we take the cylinder wall to be stationary, ucyl = 0, and look for film profiles driven only by interaction with the external multiphase flow. Steady smooth shock solutions have been found for parameter values indicated in Figure 4 when A = 4.382, the film volume per unit axial length equivalent to that for a profile with a single (forward-) backward-facing shock at (3π/4) π/4 when ζ, αu = 0, here c1 = ζ and c3 = αu. The computed profiles agree well with the forms indicated by the asymptotic smoothing analysis, except for the cases corresponding to droplet momentum being dominant, i.e.,

$\zeta <\alpha _u^3$
ζ<αu3⁠. Solutions close to the
$\zeta =\alpha _u^3$
ζ=αu3
boundary are difficult to compute and require continuation from nearby solutions. For
$\zeta <\alpha _u^3$
ζ<αu3
the asymptotic analysis indicates forward-facing shocks are feasible in π/2 < θ < π, but kinematic wave theory suggested these will be unstable due to focusing of disturbance energy at θ = π/2; it is unsurprising therefore that the numerical code failed to find them.

FIG. 4.

ζ, αu parameter space where steady solutions have been sought for A = 4.382 using the transient solver. ○ converged steady solutions, × no steady solution found. Solid line — is

$\zeta =\alpha _u^3$
ζ=αu3⁠, the boundary suggested by asymptotic analysis.

FIG. 4.

ζ, αu parameter space where steady solutions have been sought for A = 4.382 using the transient solver. ○ converged steady solutions, × no steady solution found. Solid line — is

$\zeta =\alpha _u^3$
ζ=αu3⁠, the boundary suggested by asymptotic analysis.

Close modal

Similar results to those reported in Figure 4 were obtained for other film volumes, corresponding to different shock profiles.

As expected, the effect of the surface tension parameter ζ is to smooth the film profiles. For small surface tension, the resulting profile is close to that of the cubic-shock solution, with a small capillary wave at the foot of the shock, as predicted by the asymptotic analysis of the smoothed shock. As ζ increases, the peak decreases and flattens, with the excess mass being spread around the cylinder, mostly in 0 < θ < π.

On the other hand, on increasing A beyond the critical value A = 3.9509, and keeping all the other parameters fixed, the shock amplitude increases and its shape tends to the shock solutions of (25), apart from a slight overshoot at the top of the peak and a capillary wave at the foot.

Finally, the effect of varying αu in the droplet momentum on the interfacial profiles calculated for A = 4.382 is reported in Figure 5. This is equivalent to exploring along a horizontal cut through the parameter space of Figure 4. Increasing αu has little quantitative effect on the film profile; qualitatively, there is a slight steepening of the shock, until the parameter passes through the boundary

$\alpha _u^3=\zeta$
αu3=ζ⁠, after which no steady solutions have been obtained. Although the αu component of the droplet momentum, i.e., the term proportional to c3 in (21), does not have a major effect on the film profile, it has a significant effect on its stability, as previously commented and is considered further in Sec. VI.

FIG. 5.

Profiles of h against θ illustrating the effect of increasing αu; ζ = 10−3 and A = 4.382.

FIG. 5.

Profiles of h against θ illustrating the effect of increasing αu; ζ = 10−3 and A = 4.382.

Close modal

In Sec. V various steady solutions to (24) governing the evolution of a film driven by droplet momentum were found using a transient solver. It may be anticipated that such steady-state solutions are stable, as the method inevitably involves numerical errors that effectively simulate disturbances.

In order to confirm this stability, further analysis, which highlights some interesting features, was undertaken. The transient solver is used to follow the evolution of a prescribed perturbation to one of the calculated steady solutions. In the case of perturbations that do not add mass to, or subtract from, the original total mass of the film, it was found that the perturbation decays and the film profile returns to its initial condition, i.e., the solution returns to its original stable form. If a perturbation to the film adds or subtracts mass, however, the film profile generally evolves to a new steady solution that is a close neighbour to the original. Figure 6 shows the effect of adding a substantial mass to a steady-state solution, using a perturbation of the form

\begin{eqnarray}\xi (\theta ,0) = \left\lbrace \begin{array}{ll}\cal {A} \sin \left(2\pi \frac{\theta +\pi }{\cal {W}}\right), &\quad -\pi <\theta <-\pi +\frac{\cal {W}}{2} \\[6pt]0, &\quad \mathrm{otherwise}. \end{array} \right.\end{eqnarray}
ξ(θ,0)=Asin2πθ+πW,π<θ<π+W20, otherwise .
(30)

The original profile has A = 4.382, and therefore features a smoothed shock centred on θ = π/4. The perturbation is in the form of half-sine wave with

${\cal {A}}=1.45$
A=1.45 and
${\cal {W}}=2$
W=2
. As the profile evolves with time, the extra mass flows down the falling film and joins the existing shock. In this case, the shock accommodates the additional mass by growing in size and moving down the cylinder towards θ = 0, creating a new profile with A = 5.302.

FIG. 6.

Evolution of a perturbation to a steady-state solution (A = 4.382) which adds mass, moving it to a new steady state (A = 5.302). Parameters ζ = 10−4, αu = 10−3.

FIG. 6.

Evolution of a perturbation to a steady-state solution (A = 4.382) which adds mass, moving it to a new steady state (A = 5.302). Parameters ζ = 10−4, αu = 10−3.

Close modal

To prove that a steady solution hs(θ) of (24) is linearly stable, all small perturbations must be shown to decay. Substituting

\begin{eqnarray}h=h_s(\theta )+\xi (\theta ,t),\end{eqnarray}
h=hs(θ)+ξ(θ,t),
(31)

in (24) and neglecting nonlinear terms in ξ gives

\begin{equation}\frac{\partial \xi }{\partial t}+\underbrace{\frac{\partial }{\partial \theta }\left(w\xi \right)}_{{\rm convection}} =\underbrace{\frac{\partial }{\partial \theta }\left(D\frac{\partial \xi }{\partial \theta }\right)}_{{\rm diffusion}} \ \ -\underbrace{\frac{\partial }{\partial \theta }\left(\zeta \frac{h_s^3}{3}\frac{\partial ^3\xi }{\partial \theta ^3}\right)}_{{\rm higher\, order\, diffusion}},\end{equation}
ξt+θwξ convection =θDξθ diffusion θζhs333ξθ3 higher order diffusion ,
(32)

a convection-diffusion equation with

\begin{eqnarray}w=\frac{\partial }{\partial h_s}\left(q_s\left(h_s,\frac{\mathrm{d} h_s}{\mathrm{d} \theta },\frac{\mathrm{d} ^3h_s}{\mathrm{d} \theta ^3}\right)\right)=\zeta h_s^2\left(\frac{\mathrm{d} ^3h_s}{\mathrm{d} \theta ^3}+\frac{\mathrm{d} h_s}{\mathrm{d} \theta }\right)-h_s^2\sin \theta + \left(1+\alpha _u \frac{\mathrm{d} h_s}{\mathrm{d} \theta }\right)h_s,\end{eqnarray}
w=hsqshs,dhsdθ,d3hsdθ3=ζhs2d3hsdθ3+dhsdθhs2sinθ+1+αudhsdθhs,
(33)
\begin{eqnarray}D=-\left(\alpha _u \frac{h_s^2}{2}+\zeta \frac{h_s^3}{3}\right)<0,\end{eqnarray}
D=αuhs22+ζhs33<0,
(34)

and

\begin{eqnarray}q_s=\zeta \frac{h_s^3}{3}\left(\frac{\mathrm{d} ^3h_s}{\mathrm{d} \theta ^3}+\frac{\mathrm{d} h_s}{\mathrm{d} \theta }\right) -\frac{h_s^3}{3}\sin \theta + \left(1+\alpha _u \frac{\mathrm{d} h_s}{\mathrm{d} \theta }\right)\frac{h_s^2}{2}.\end{eqnarray}
qs=ζhs33d3hsdθ3+dhsdθhs33sinθ+1+αudhsdθhs22.
(35)

In (32), D is a negative diffusion and always has a destabilizing effect: here both the αu contribution from the droplet momentum and the surface tension ζ act as destabilizing agents. The surface tension, however, can also be stabilizing, acting through the higher order diffusivity term

$\zeta {h_s^3}/{3}$
ζhs3/3⁠. Without higher order diffusivity, (32) is an ill-posed problem and is unstable for all possible wavenumbers, but with it, (32) is an almost ill-posed problem, which may be stable for some wavelengths and unstable for some others, as discussed by Hsieh and Ho.35 

Taking a perturbation of the form ξ = η(θ)eσt, where the initial amplitude η is 2π-periodic and the exponent σ is complex, the stability problem was reduced to a discrete eigenvalue problem of the form

$\sigma \vec{\eta } = B \vec{\eta }$
ση=Bη⁠, by using a spectral collocation method, Trefethen.36 

All the eigenvalues of B were found by using MATLAB routine eig.m, or a few eigenvalues by using eigs.m; the routine also returns the corresponding eigenmodes. If the eigenvalues σ have positive real part, then the perturbations grow in time and the corresponding eigenmodes are unstable. If there are no eigenvalues with positive real part, then the solution hs is stable. Note that if there is only one eigenvalue with positive real part and the integral of the corresponding eigenmode is non-zero, then the solution is linearly stable to perturbations that conserve the mass of the film. This is the situation in all the cases that we have investigated.

For example, for the case with αu = 10−3, ζ = 10−3, and A = 4.382, there is a real eigenvalue, 4.543; the corresponding unstable eigenmode has non-zero total mass, in strong contrast to all the stable eigenmodes which have zero total mass. This result is typical for the stable solutions from the parameter range

$\alpha _u^3<\zeta$
αu3<ζ with A > Acrit – the film volume is supercritical and, as
$\alpha _u^3<\zeta$
αu3<ζ
, this case lies in the parameter range where a surface-tension-smoothed-shock solution is feasible by the matching argument. We conjecture that all steady solutions found using the transient code are linearly stable to perturbations with zero total mass. To explore further, a series of tests was conducted, including increasing the number of grid points, comparing results from transient and steady-state numerical solvers, and comparing results across a range of values of αu and ζ. All eigenmodes with zero total mass were found to be stable: the region of ζ, αu parameter space for which solutions are stable with A = 4.382 is shown in the stability map in Table III. The symbol × indicates that no steady solution was found, 1 indicates that a single positive real eigenvalue with an associated non-zero-total-mass unstable eigenmode was found. There is an exact correspondence between the present results and those in Figure 4, in the sense that steady-state solutions that are unstable to perturbations adding or subtracting mass to the system, evolve to new steady solutions with corresponding new mass.

Table III.

Parameter map of ζ and αu indicating number of eigenvalues with positive real part, for A = 4.382.

ζ\αu10−510−410−310−210−11
× 
10−1 × 
10−2 × × 
10−3 × × 
10−4 × × 
10−5 × × × 
10−6 × × × 
ζ\αu10−510−410−310−210−11
× 
10−1 × 
10−2 × × 
10−3 × × 
10−4 × × 
10−5 × × × 
10−6 × × × 

With the c5 and c6 terms of (21) different from zero, the droplets transfer mass as well as momentum to the film flow. As before, the focus is on cases where the film is driven solely by the impacting droplets, i.e., the wall is stationary.

The parameter values c1 to c4, and the film volume A, that have been shown previously (in Sec. V) to produce a stable steady solution when c5, c6 = 0, are all fixed. Solutions with values of c5 ⩽ 0 and c6 ⩽ 0, corresponding to mass addition to the film from droplets travelling in the azimuthal and radial directions, respectively, are explored. Recall that, according to the distinguished-limit analysis, mass transfer can become as significant as momentum transfer by increasing the average density or increasing the inclination angle of the incoming droplets.

Whereas, in the steady model studied, there is no net mass flux into the film and the emphasis has been on the flow into the film through the surface and out of the film through the wall sink, there is nothing in principle to prevent mass leaving the film through its surface as well. On a backslope of the interface, the term auEh0/∂θ is negative which might suggest that mass leaves the surface of the film. This flux should, however, be considered in conjunction with the positive term −avE, which, in the model, always adds mass to the film. Even when the gradient of h is negative, mass will not leave if

$\partial {h_0}/\partial {\theta }>v_E/u_E$
h0/θ>vE/uE⁠, although the net flux into the film will be lower than when the gradient of h is positive. For all of the cases considered herein the net local flux into the film is everywhere positive.

As previously, the evolution equation is solved numerically by a fully implicit backward-time-difference scheme, using the MATLAB function bvp4c. The code incorporates a single sink modelled numerically as a small finite section of the cylinder wall, across which there is a constant outflow of fluid. The outflow function used is

\begin{eqnarray}v_{cyl}= \left\lbrace \begin{array}{ll}0 &\theta <\theta _s, \\[6pt]v_{out} & \theta _s<\theta <\theta _e,\\[6pt]0 &\theta _e<\theta , \\\end{array} \right.\end{eqnarray}
vcyl=0θ<θs,voutθs<θ<θe,0θe<θ,
(36)

where for the cases of a single sink located at θ = 0, θs = −π/100, and θe = π/100 are used; in all cases ucyl = 0.

To study steady solutions, the total mass outflow is required to match the total mass inflow. The total mass entering the system is given by

\begin{eqnarray}\int _0^{2\pi }(c_5\frac{\partial h}{\partial \theta }+c_6 )\,\mathrm{d}\theta = 2 \pi c_6 ,\end{eqnarray}
02π(c5hθ+c6)dθ=2πc6,
(37)

since h is 2π-periodic and c5 and c6 are constants, and the outflow velocity at the sink is

\begin{equation}v_{out} = \frac{2 \pi }{\theta _e-\theta _s}c_6.\end{equation}
vout=2πθeθsc6.
(38)

As before, the numerical code was validated against the results reported by Noakes et al.9 for added mass but negligible droplet momentum.

In order to calculate residence times and to visualize the flow field, to check, for example, whether there is any recirculation, information on the streamlines is required. From (B21) and (B23), the leading-order streamfunction ψ is

\begin{eqnarray}\psi = & - & \left(c_1\left(\frac{\partial h_0}{\partial \theta }+\frac{\partial ^3 h_0}{\partial \theta ^3}\right) +c_2\sin \theta \right)\frac{y^3}{6} \nonumber \\& + & \left(\left(c_1\left(\frac{\partial h_0}{\partial \theta }+\frac{\partial ^3 h_0}{\partial \theta ^3}\right) +c_2\sin \theta \right)h_0+c_3\frac{\partial h_0}{\partial \theta }+c_4\right)\frac{y^2}{2}-\int v_{cyl}(\theta ) \mathrm{d}\theta .\end{eqnarray}
ψ=c1h0θ+3h0θ3+c2sinθy36+c1h0θ+3h0θ3+c2sinθh0+c3h0θ+c4y22vcyl(θ)dθ.
(39)

Figure 7 illustrates the steady-state streamlines when mass is extracted at a sink at θ = 0, a = 10−6, uE = 102, vE = −104, A = 4.382, ζ = 10−3, αu = 10−2, k = 1, and β = 1. The surface profile for the corresponding dominant-droplet-momentum case (without a sink) with k = 2 is shown in Figure 5. By setting k = 1, the droplet volume fraction in the core flow is increased and with all the other parameters fixed, the mass transfer is made as significant as momentum (see Table I) and results in the flow field shown in Figure 7. The dip in the surface elevation as the film passes over the sink is similar to that found by Momoniat et al.37 in the radial spreading under gravity and surface tension of a thin film under the influence of slot suction. By following a streamline, the path of a given fluid particle from its entry into the film surface until it exits at the sink can be traced – the path is shown by the thick black line in Figure 7. The fluid entering the surface near K goes around the cylinder between 2 and 3 times before leaving at the sink at N. There is recirculation, in a small region between G and H at the cylinder wall around θ = 1.

FIG. 7.

Streamlines: sink at θ = 0, a = 10−6, uE = 102, vE = −104, A = 4.382, ζ = 10−3, αu = 10−2, k = 1, and β = 1.

FIG. 7.

Streamlines: sink at θ = 0, a = 10−6, uE = 102, vE = −104, A = 4.382, ζ = 10−3, αu = 10−2, k = 1, and β = 1.

Close modal

To be consistent with the formulation, values of c5 and c6 are considered for which the ratios c5 : c6 and c3 : c4 are the same (see Table I). For a given value of A and with c2 = −1, c3 ≡ αu = 10−3, and c4 = 1, the values of c1 ≡ ζ are chosen so that stable momentum-dominant solutions exist; also c5 = 10−3 × c6 so that the azimuthal mass contribution is small, but physically consistent. Again fluid is withdrawn through the cylinder wall around θ = 0 at a rate to match the input. From the results in Sec. V for q = 1/6, we expect to see solutions with a shock in 0 < θ < π/2 if the volume of the film exceeds the critical value A = 3.951.

For fixed A = 4.382 and increasing mass flow into and out of the film, the mass of fluid ahead of the sink increases to provide the required pressure head for the given outflow, and correspondingly the mass behind the sink decreases, see Figure 8. For |c6| ⩾ 4 × 10−2 a backward-facing shock no longer features.

FIG. 8.

Film profiles h when A = 4.382, various c5 = 10−3 × c6. Other parameters are c1 = 10−4, c2 = −1, c3 = 10−3, c4 = 1, ucyl = 0.

FIG. 8.

Film profiles h when A = 4.382, various c5 = 10−3 × c6. Other parameters are c1 = 10−4, c2 = −1, c3 = 10−3, c4 = 1, ucyl = 0.

Close modal

Figure 9(a) shows the streamlines for a case with small droplet mass transfer when A = 4.382 and c6 = −0.01. As in the case of dominant droplet momentum, there is a broad recirculation region, and stagnation points on the boundary. The injected fluid is resident for 3.5–4.5 cycles. Figure 9(b) illustrates streamlines for a film with the same film volume but much greater droplet mass input – c6 = −0.07, the minimum for which a solution can be found – with a maximum residence of just over one cycle. In this case, so much fluid is built up on the upstream of the sink (θ < 0) that there is no shock in θ > 0.

FIG. 9.

Streamlines when (a) c5 = −10−5, c6 = −10−2, (b) c5 = −7 × 10−5, c6 = −7 × 10−2, and A = 4.382. Other parameters are c1 = 10−4, c2 = −1, c3 = 10−3, c4 = 1, ucyl = 0.

FIG. 9.

Streamlines when (a) c5 = −10−5, c6 = −10−2, (b) c5 = −7 × 10−5, c6 = −7 × 10−2, and A = 4.382. Other parameters are c1 = 10−4, c2 = −1, c3 = 10−3, c4 = 1, ucyl = 0.

Close modal

Figure 10 shows the film profiles for a range of A, the film volume per unit axial length, when droplet mass parameters are c5 = −10−5 and c6 = −10−2. Steady-state solutions were found for A = 2.0 to A = 5.991. In all cases, there is a small build up of fluid mass before the sink. Smooth profiles are obtained for lower values of A, and, as A increases, the profile develops a hump around θ = π/2, until a value of A ≈ 4.0, when it is close to the critical profile. A further increase in mass is incorporated into a backward-facing shock in 0 < θ < π/2. The loss of mass through the sink makes the film thin immediately after θ = 0; note that very thin films were difficult to compute. Continuation from nearby profiles was implemented, but no films were found with A < 2.0. For Amax = 5.991, the largest film volume for which a steady profile with these parameters could be obtained numerically, the shock is so close to the sink that the capillary wave at its foot is close to interacting with the film hump before the sink. This result demonstrates the sharp change that results from the inclusion of mass transfer, since A → ∞ when it is absent (cf. (27)).

FIG. 10.

Film profile h when c5 = −10−5, c6 = −10−2, various A. Other parameters are c1 = 10−4, c2 = −1, c3 = 10−3, c4 = 1, ucyl = 0.

FIG. 10.

Film profile h when c5 = −10−5, c6 = −10−2, various A. Other parameters are c1 = 10−4, c2 = −1, c3 = 10−3, c4 = 1, ucyl = 0.

Close modal

In designing an engine bearing chamber, it may be possible to re-position the sink through which oil leaves the chamber. Thus far we have assumed that the sink is at the bottom of the cylinder, at θ = 0; other sink positions are now examined to establish how varying the sink position affects the film. For all the simulations in this section, the parameters are A = 4.382, c1 = 10−3, c2 = −1, c3 = 10−2c4 = 1, c5 = −10−4, c6 = −10−2, and vcyl is defined by (38) for different values of sink position.

Figure 11 shows the film profiles and streamlines for sink locations ranging from θ = −3π/4 to θ = 3π/4. For sink positions −π < θ < 0, the sink is well away from the ascending shock and the solution behaves in a very similar way to that already discussed; there is no qualitative effect on the film profile in moving the sink. When the sink is in the quadrant 0 < θ < π/2 then it affects the shock. For shock positions in π/4 < θ < 1.25 the sink causes the height of the shock to diminish, but there is still recirculation downstream of the sink position. For θ = 1.3, the peak begins to emerge upstream of the sink, and there is recirculation on both sides of the sink. For θ ⩾ 1.35, all the recirculation is upstream of the sink. So for these parameters there is no sink position that completely removes the recirculation, although it diminishes greatly when the sink is under the shock.

FIG. 11.

Streamlines for various sink positions. Parameters are c1 = 10−3, c2 = −1, c3 = 10−2, c4 = 1, c5 = −10−4, c6 − 10−2, A = 4.382.

FIG. 11.

Streamlines for various sink positions. Parameters are c1 = 10−3, c2 = −1, c3 = 10−2, c4 = 1, c5 = −10−4, c6 − 10−2, A = 4.382.

Close modal

In this study, a model of a two-dimensional liquid film flowing around the inside boundary of a moving or stationary smooth rigid circular cylinder has been developed in the case when there is a dilute two-phase flow in the core of the cylinder and a strong interaction with droplets of the same liquid impacting on the free surface. The bulk flow of the film is taken as Newtonian and laminar and the central core of two-phase flow to be inviscid. Interface conditions based on conservation of mass and momentum have been derived.

Estimates of appropriate physical parameters to determine scalings and the development of simplified film models through distinguished limits are presented. In particular, the fact that the film has small aspect ratio (film thickness to azimuthal length scale) is exploited to simplify greatly the governing equations and boundary conditions. Various parameter ranges for the droplets are identified, and in each case the corresponding governing equation for the film evolution is derived. A single evolution equation for the film profile, (21),

\begin{eqnarray*}\frac{\partial h}{\partial t} +\frac{\partial }{\partial \theta }\left( c_1\left( \frac{\partial h}{\partial \theta }+\frac{\partial ^3 h}{\partial \theta ^3} \right)\frac{h^3}{3} + c_2\sin \theta \frac{h^3}{3} + \left(c_3(\theta )\frac{\partial h}{\partial \theta }+c_4(\theta )\right)\frac{h^2}{2}\right)\nonumber \\+ \left(c_5(\theta )\frac{\partial h}{\partial \theta } +c_6(\theta )\right) + \frac{\partial }{\partial \theta }\left(u_{cyl}h\right) - v_{cyl}=0, \nonumber\end{eqnarray*}
ht+θc1hθ+3hθ3h33+c2sinθh33+c3(θ)hθ+c4(θ)h22+c5(θ)hθ+c6(θ)+θucylhvcyl=0,

is obtained from which specific cases can be obtained by appropriate choice of coefficients.

The c1 coefficient may be small but, being a function of Wef, is independent of other coefficients and it is retained in order to study its effect, since it acts on the highest derivative of h. The term involving c2 corresponds to the contribution of gravity. The coefficient c4 in the droplet momentum is identical in all cases and is analogous to an imposed shear on the surface, such as that studied by Villegas-Díaz et al.12 

The term with coefficient c3 is analogous to the temperature-dependent surface-tension term that arises in the problem of flow down a heated plane, see Thiele and Knobloch.38 The terms c5 and c6 arise from the contribution of droplet mass flux. The derivative term arises from the azimuthal component of droplet mass, and the constant term from the radial component.

In this study the focus is on cases where the impact velocity of droplets into a film is sufficiently high for the droplet momentum to drive the motion of the film in the absence of rotation of the cylinder wall.

To fix on the key terms of surface tension and droplet momentum, the equation was re-scaled and a new set of coefficients, c1 = ζ, c2 = −1, c3 = αu, c4 = 1, c5 = 0, c6 = 0, in (21). When ζ, αu = 0, the flux equation reduces to a cubic, which has at most two positive real solutions, depending on A, the volume of the film. When AAc = 3.9509, there is a single completely wetting solution for the film. When A > Ac, shock solutions may occur, with the mass in excess of critical accumulating in a shock somewhere in 0 < θ < π.

A matching argument suggests that, for ζ, αu > 0, only certain combinations of shocks are mathematically feasible. If

$\alpha _u^3 \le O(\zeta )$
αu3O(ζ) as ζ → 0, smoothed backward- and forward-facing shocks located in 0 < θ < π are feasible; if
$\zeta \ll \alpha _u^3 \ll 1$
ζαu31
as αu → 0, a smoothed forward-facing shock located in π/2 < θ < π is feasible. Kinematic wave theory suggests that profiles with a single backward-facing shock in the lower quadrant (0 < θ < π/2) are stable, profiles with a single forward-facing shock in the upper quadrant (π/2 < θ < π) unstable, and profiles featuring multiple shock solutions may be stable (unstable) if the profile height is decreasing (increasing) as θ increases through π/2.

These results are strongly supported by the numerical results obtained from a transient solver. Linear stability analysis confirmed the theoretical and transient-solver results and highlighted the significance of the

$\zeta = \alpha _u^3$
ζ=αu3 boundary in parameter space. Interestingly, the stability analysis revealed positive real eigenvalues, and hence the existence of potentially unstable perturbation modes. On closer examination, however, the modes were found not to be mass-conserving in the sense that a system of fixed mass fraction cannot admit such modes.

In addition to the effect of droplet momentum transfer, the effect of mass transfer has also been considered, i.e., the effect of the terms c5 and c6 (azimuthal and radial components of droplet mass flux) in the film-height evolution equation (21). The values of c5 and c6 were varied independently but in such a way that they were consistent, in terms of the physics, with c3 and c4. Having chosen a set of coefficients c1 to c4 and A for which momentum-dominant steady-state solutions are known to be stable, mass transfer effects were considered by including c5 and c6. In order to achieve a steady-state condition, the film volume, A, is kept constant by introducing a sink at the cylinder wall. When c5 and c6 are constant, as in the cases considered in this study, the mass entering the system is linearly proportional to c6 and the strength of the sink is therefore trivial to determine. Plots of the calculated streamlines reveal the underlying steady-flow behaviour and enable the residence time of fluid that enters through the film surface to be determined. Importantly from the perspective of oil degradation, closed recirculation zones that never receive fresh fluid feature behind the main shock structures. The extent of this closed recirculation can be diminished, but not eliminated, by re-positioning the sink to lie beneath the shock.

On increasing the mass flux through the film, i.e., increasing the value of c6, the mass of fluid ahead of the sink increases to provide the required pressure head for the required balancing outflow; correspondingly, the mass behind the sink decreases. It was observed that an increase in droplet mass flux reduces the footprint of the backward-facing shock and recirculation zone.

The effect of c5 is more subtle. By increasing c5 substantially, steady-state solutions may cease to exist, as has been reported previously in cases of a film on a rotating wall without droplet impact. It was found that Acrit, the film volume above which a shock solution is required, and Amax, the maximum film volume sustained by a stable steady shock solution, are dependent on c5.

In the film evolution equation (21) the terms c5h/∂θ and ucylh/∂θ have the same form in h and are therefore interchangeable if the film height profile h is of interest. However, the first form arises from the azimuthal mass contribution at the film surface, and the second from the rotating wall boundary condition. Therefore, they are distinct in the calculation of the underlying film velocities u and v, and the two flows with c5 ≠ 0, ucyl = 0, and c5 = 0, ucyl ≠ 0 have distinct behaviours.

The authors wish to acknowledge the financial support provided by Rolls-Royce plc, Aerospace Group as part of University Technology Centre in Gas Turbine Transmission Systems at the University of Nottingham. The views expressed in this paper are those of the authors and not necessarily those of Rolls-Royce plc, Aerospace Group. The authors would like to thank one of the referees for very careful reading of the article and helpful feedback.

Consider conservation of mass and momentum near the surface separating the two-phase flow of air and droplets and the film flow (in reality the two flows are separated by a complex, narrow transitional layer, which is modelled as a simple surface, see Edwards et al.39). This interface is represented as a moving surface, a surface at which density and velocity are discontinuous and mass and momentum transfer across it (see Slattery40). Consider a material pillbox straddling the moving and deforming interface between the two-phase flow and the fluid film; the pillbox has endcaps of small arbitrary length parallel to the interface, and thickness δ across the interface.

In contrast to interface conditions between immiscible fluids, where, as the interface develops, fluid does not cross the interface, in the present case the droplets in the two-phase flow may cross the interface to become part of the liquid film. In order to apply the mass and momentum conservation laws in the present case, the pillbox comprises two elements, one with volume VE(t) containing an inviscid mixture of air and droplets and the other with volume VI(t) containing the viscous fluid, which includes any contribution from droplets impacting on the interface. (Note that for simplicity the hats on dimensional variables are dropped in this section.)

1. Conservation of mass

Mass conservation in the control volume VEVI gives

\begin{equation}\frac{D}{D t}\int _{V_E(t)}\rho _E \, {\rm d}V \; + \; \frac{D}{D t}\int _{V_I(t)}\rho _I \, {\rm d}V=0.\end{equation}
DDtVE(t)ρEdV+DDtVI(t)ρIdV=0.
(A1)

By using Reynolds's transport theorem and taking ρE = αρ and ρI = ρ, the above two terms on the left-hand side can be rewritten as

\begin{equation}\int _{V_E(t)}\frac{\partial (\alpha \rho )}{\partial t}+\frac{\partial (\alpha \rho u_{j})}{\partial x_j} \mathrm{d}V \;\; {\mbox{and}} \;\; \int _{V_I(t)}\frac{\partial \rho }{\partial t}+\frac{\partial (\rho u_{j})}{\partial x_j} \, {\rm d}V,\end{equation}
VE(t)(αρ)t+(αρuj)xjdVandVI(t)ρt+(ρuj)xjdV,
(A2)

or, on using Gauss's theorem,

\begin{equation}\int _{V_E(t)}\frac{\partial (\alpha \rho )}{\partial t} \, {\rm d}V + \int _{S_E(t)}\alpha \rho u_{j} \tilde{n}_j \, {\rm d}S \;\; {\mbox{and}} \;\; \int _{V_I(t)}\frac{\partial \rho }{\partial t} \, {\rm d}V + \int _{S_I(t)}\rho u_{j} \tilde{n}_j \, {\rm d}S,\end{equation}
VE(t)(αρ)tdV+SE(t)αρujñjdSandVI(t)ρtdV+SI(t)ρujñjdS,
(A3)

where SE(t)and SI(t) are the surface areas enclosing VE(t) and VI(t), respectively, and

$\tilde{n}_j$
ñj denotes the normal pointing out of the volume. In the limit δ → 0, the volumetric integrals vanish and these two expressions reduce to

\begin{equation}\int _{S_{int}(t)}\alpha \rho (u_{Ej} {n}_j-U_j {n}_j) \, {\rm d}S\end{equation}
Sint(t)αρ(uEjnjUjnj)dS
(A4)

and

\begin{equation}\int _{S_{int}(t)}\rho (U_j {n}_j-u_{Ij} {n}_j) \, {\rm d}S,\end{equation}
Sint(t)ρ(UjnjuIjnj)dS,
(A5)

where Sint is the common surface at the interface between the two volumes, Uj is the interface velocity, and nj denotes the normal to the interface pointing into the two-phase flow.

The velocity fields may not be continuous across the interface: let uEj/uIj denote the respective limiting forms of the two-phase/film velocity fields at a point approaching Sint from within the respective volumes.

Re-combining (A4) and (A5), we obtain

\begin{equation}\int _{S_{int}(t)}[\rho (1-\alpha ) U_j n_j+\rho (\alpha u_{Ej}-u_{Ij}) n_j] \, {\rm d}S =0.\end{equation}
Sint(t)[ρ(1α)Ujnj+ρ(αuEjuIj)nj]dS=0.
(A6)

Since the interface segment Sint is arbitrary,

\begin{equation}\rho (\alpha -1) U_j n_j=\rho (\alpha u_{Ej}-u_{Ij}) n_j.\end{equation}
ρ(α1)Ujnj=ρ(αuEjuIj)nj.
(A7)
2. Conservation of momentum

Applying the principle of conservation of momentum to the two volumes defined above, we obtain

\begin{equation}\frac{\mathrm{D}}{\mathrm{D}t}\int _{V_E(t)} \alpha \rho u_{i} \, {\rm d}V= \int _{V_E(t)}\alpha \rho g_i \, {\rm d}V + \int _{S_E(t) \setminus S_{int}(t)} T_{ij}\tilde{n}_j \, {\rm d}S + \int _{S_{int}(t)} \widetilde{T}_{Eij}\tilde{n}_j \, {\rm d}S\end{equation}
DDtVE(t)αρuidV=VE(t)αρgidV+SE(t)Sint(t)TijñjdS+Sint(t)T̃EijñjdS
(A8)

and

\begin{equation}\frac{\mathrm{D}}{\mathrm{D}t}\int _{V_I(t)} \rho u_{i} \, {\rm d}V= \int _{V_I(t)} \rho g_i \, {\rm d}V + \int _{S_I(t) \setminus S_{int}(t)} T_{ij}\tilde{n}_j \, {\rm d}S + \int _{S_{int}(t)} \widetilde{T}_{Iij}\tilde{n}_j \, {\rm d}S.\end{equation}
DDtVI(t)ρuidV=VI(t)ρgidV+SI(t)Sint(t)TijñjdS+Sint(t)T̃IijñjdS.
(A9)

Here, the volume integrals on the right-hand-sides correspond to the gravitational body force, and the first surface integrals on the left-hand-sides to the surface traction: Tij = −pIδij + 2μeij is the Newtonian stress tensor in the viscous film, Tij = −pEδij in the inviscid mixture of air and droplets. In the above,

$\widetilde{T}_{Eij}$
T̃Eij and
$\widetilde{T}_{Iij}$
T̃Iij
are representative interfacial stresses such that
$(\widetilde{T}_{Eij}-\widetilde{T}_{Iij})n_j = -\sigma \kappa n_i$
(T̃EijT̃Iij)nj=σκni
, where σ is an effective surface tension coefficient.

As before, on using Reynolds's transport and Gauss's theorems, these two identities reduce to

\begin{eqnarray}\int _{V_E(t)}\frac{\partial (\alpha \rho u_{i})}{\partial t} \, {\rm d}V+ \int _{S_E(t)}(\alpha \rho u_{i})u_{j} \tilde{n}_j \, {\rm d}S = \int _{V_E(t)}\alpha \rho g_i \, {\rm d}V + \int _{S_E(t)\setminus S_{int}(t)} T_{ij}\tilde{n}_j \, {\rm d}S \nonumber \\+ \int _{S_{int}(t)} \widetilde{T}_{Eij}\tilde{n}_j \, {\rm d}S\end{eqnarray}
VE(t)(αρui)tdV+SE(t)(αρui)ujñjdS=VE(t)αρgidV+SE(t)Sint(t)TijñjdS+Sint(t)T̃EijñjdS
(A10)

and

\begin{eqnarray}\int _{V_I(t)}\frac{\partial (\rho u_{i})}{\partial t} \, {\rm d}V+ \int _{S_I(t)}(\rho u_{i})u_{j}\tilde{n}_j \, {\rm d}S = \int _{V_I(t)}\rho g_i \, {\rm d}V + \int _{S_I(t)\setminus S_{int}(t)} T_{ij}\tilde{n}_j \, {\rm d}S \nonumber \\+ \int _{S_{int}(t)} \widetilde{T}_{Iij}\tilde{n}_j \, {\rm d}S.\end{eqnarray}
VI(t)(ρui)tdV+SI(t)(ρui)ujñjdS=VI(t)ρgidV+SI(t)Sint(t)TijñjdS+Sint(t)T̃IijñjdS.
(A11)

Taking δ → 0, the above two expressions reduce to

\begin{equation}\int _{S_{int}(t)}(\alpha \rho u_{Ei})(u_{Ej}-U_j) n_j \, {\rm d}S +\int _{S_{int}(t)}(\widetilde{T}_{Eij}-T_{Eij})n_j \, {\rm d}S=0\end{equation}
Sint(t)(αρuEi)(uEjUj)njdS+Sint(t)(T̃EijTEij)njdS=0
(A12)

and

\begin{equation}\int _{S_{int}(t)}(\rho u_{Ii})(U_j-u_{Ij}) n_j \, {\rm d}S +\int _{S_{int}(t)}(T_{Iij}-\widetilde{T}_{Iij})n_j \, {\rm d}S=0;\end{equation}
Sint(t)(ρuIi)(UjuIj)njdS+Sint(t)(TIijT̃Iij)njdS=0;
(A13)

the corresponding volume integrals tend to zero as δ → 0. In the above limit we take into account the jump property of the velocity field across the interface and a similar jump condition for the surface traction as the endcap surfaces of the pillboxes tend to the interface Sint, with limiting values TEij and TIij. In the limit we allow for a discontinuity in the traction across the interface, i.e.,

$\widetilde{T}_{Iij}n_j$
T̃Iijnj and
$\widetilde{T}_{Eij}n_j$
T̃Eijnj
, respectively.

By combining the above two expressions, we obtain

\begin{eqnarray}\int _{S_{int}(t)}[-\rho (\alpha u_{Ei}-u_{Ii})U_j n_j+ & \rho (\alpha u_{Ei}u_{Ej}-u_{Ii}u_{Ij})n_j \nonumber \\+ (\widetilde{T}_{Eij}-\widetilde{T}_{Iij})n_j- & (T_{Eij}-T_{Iij})n_j] \, {\rm d}S=0.\end{eqnarray}
Sint(t)[ρ(αuEiuIi)Ujnj+ρ(αuEiuEjuIiuIj)nj+(T̃EijT̃Iij)nj(TEijTIij)nj]dS=0.
(A14)

Thus, for arbitrary Sint,

\begin{equation}-\rho (\alpha u_{Ei}-u_{Ii})U_j n_j+\rho (\alpha u_{Ei}u_{Ej}-u_{Ii}u_{Ij})n_j - (T_{Eij}-T_{Iij})n_j=-(\widetilde{T}_{Eij}-\widetilde{T}_{Iij})n_j\end{equation}
ρ(αuEiuIi)Ujnj+ρ(αuEiuEjuIiuIj)nj(TEijTIij)nj=(T̃EijT̃Iij)nj
(A15)

or

\begin{equation}-\rho (\alpha u_{Ei}-u_{Ii})U_j n_j+\rho (\alpha u_{Ei}u_{Ej}-u_{Ii}u_{Ij})n_j - (T_{Eij}-T_{Iij})n_j=\sigma \kappa n_i.\end{equation}
ρ(αuEiuIi)Ujnj+ρ(αuEiuEjuIiuIj)nj(TEijTIij)nj=σκni.
(A16)

Using the mass interface condition (A7) and grouping terms to further simplify the resulting expression leads to the condition for momentum conservation across the interface

\begin{equation}\frac{\alpha \rho }{1-\alpha } (u_{Ei}-u_{Ii})(u_{Ej}-u_{Ij})n_j-(T_{Eij}-T_{Iij})n_j=\sigma \kappa n_i\end{equation}
αρ1α(uEiuIi)(uEjuIj)nj(TEijTIij)nj=σκni
(A17)

or, in jump notation,

\begin{equation}[T_{ij}]n_j=-\sigma \kappa n_i +\frac{\alpha \rho }{1-\alpha }[u_{i}][u_j]n_j.\end{equation}
[Tij]nj=σκni+αρ1α[ui][uj]nj.
(A18)

In (A18), the final term arises from conservation of momentum, and the left-hand-side from pressure and viscous forces.

Since the droplet fraction is low, the two-phase flow is taken as inviscid, and consequently the above surface traction jump condition can be written in terms of its normal and tangential components as

\begin{eqnarray}-\sigma \kappa = &\displaystyle {\frac{\alpha \rho }{1-\alpha }\left(u_{E\theta }\frac{1}{r_0-h}\frac{\partial h}{\partial \theta } +u_{Er}-u_{I\theta }\frac{1}{r_0-h}\frac{\partial h}{\partial \theta }-u_{Ir}\right)^2N^2 +p_E-p_I} \nonumber \\& + \displaystyle {\quad 2\mu \left( \left(\frac{1}{r_0-h}\frac{\partial h}{\partial \theta }\right)^2\left(\frac{1}{r_0-h}\frac{\partial u_{I\theta }}{\partial \theta } +\frac{u_{Ir}}{r_0-h}\right) \right.} \nonumber \\& \quad \quad \quad + \displaystyle {\left. \left(\frac{1}{r_0-h}\frac{\partial h}{\partial \theta }\right)\left(\frac{\partial u_{I\theta }}{\partial r} -\frac{u_{I\theta }}{r_0-h}+\frac{1}{r_0-h}\frac{\partial u_{Ir}}{\partial \theta }\right) +\frac{\partial u_{Ir}}{\partial r} \right)N^2}\end{eqnarray}
σκ=αρ1αuEθ1r0hhθ+uEruIθ1r0hhθuIr2N2+pEpI+2μ1r0hhθ21r0huIθθ+uIrr0h+1r0hhθuIθruIθr0h+1r0huIrθ+uIrrN2
(A19)

and

\begin{eqnarray}0&=& \displaystyle {\frac{\alpha \rho }{1-\alpha } \left(-u_{E\theta } +u_{Er}\frac{1}{r_0-h}\frac{\partial h}{\partial \theta } +u_{I\theta }-u_{Ir}\frac{1}{r_0-h}\frac{\partial h}{\partial \theta }\right)} \nonumber \\&&\quad \displaystyle {\times \left(u_{E\theta }\frac{1}{r_0-h}\frac{\partial h}{\partial \theta } +u_{Er} -u_{I\theta }\frac{1}{r_0-h}\frac{\partial h}{\partial \theta }-u_{Ir}\right)N^2} \nonumber \\&& \displaystyle {+ \mu \left(\left( \left(\frac{1}{r_0-h}\frac{\partial h}{\partial \theta }\right)^2 - 1\right) \left(\frac{\partial u_{I\theta }}{\partial r}-\frac{u_{I\theta }}{r_0-h} +\frac{1}{r_0-h}\frac{\partial u_{Ir}}{\partial \theta }\right) \right.} \nonumber \\&& \quad \displaystyle {- \left. 2\left(\frac{1}{r_0-h}\frac{\partial h}{\partial \theta }\right) \left(\frac{1}{r_0-h}\frac{\partial u_{I\theta }}{\partial \theta } +\frac{u_{Ir}}{r_0-h}-\frac{\partial u_{Ir}}{\partial r}\right)\right)N^2}.\end{eqnarray}
0=αρ1αuEθ+uEr1r0hhθ+uIθuIr1r0hhθ×uEθ1r0hhθ+uEruIθ1r0hhθuIrN2+μ1r0hhθ21uIθruIθr0h+1r0huIrθ21r0hhθ1r0huIθθ+uIrr0huIrrN2.
(A20)

In dimensionless form, the Navier-Stokes and continuity equations are

\begin{eqnarray}\epsilon \frac{\partial u}{\partial t}(1-\epsilon y)^2 + \epsilon (1-\epsilon y) u\frac{\partial u}{\partial \theta } + \epsilon v(1-\epsilon y)^2\frac{\partial u}{\partial y} \nonumber \\-\epsilon ^2(1-\epsilon y) u v + (1-\epsilon y)\frac{\partial p}{\partial \theta } \nonumber \\= {\mathrm{Re_f}^{-1}}\left(\epsilon ^2\frac{\partial ^2 u}{\partial \theta ^2} -2\epsilon ^3\frac{\partial v}{\partial \theta } +(1-\epsilon y)^2\frac{\partial ^2 u}{\partial y^2} -\epsilon (1-\epsilon y)\frac{\partial u}{\partial y} -\epsilon ^2 u \right) \nonumber \\-(1-\epsilon y)^2{\mathrm{Fr^{-1}}}\sin \theta ,\end{eqnarray}
εut(1εy)2+ε(1εy)uuθ+εv(1εy)2uyε2(1εy)uv+(1εy)pθ= Re f1ε22uθ22ε3vθ+(1εy)22uy2ε(1εy)uyε2u(1εy)2 Fr 1sinθ,
(B1)
\begin{eqnarray}\epsilon ^3(1-\epsilon y)^2\frac{\partial v}{\partial t} +\epsilon ^3(1-\epsilon y) u\frac{\partial v}{\partial \theta } +\epsilon ^3(1-\epsilon y)^2 v\frac{\partial v}{\partial y} \nonumber \\+\epsilon ^2(1-\epsilon y) u^2 +(1-\epsilon y)^2\frac{\partial p}{\partial y} \nonumber \\= {\epsilon ^2}{\mathrm{Re_f}^{-1}} \left(\epsilon ^2\frac{\partial ^2 v}{\partial \theta ^2} +2\epsilon \frac{\partial u}{\partial \theta } +(1-\epsilon y)^2\frac{\partial ^2 v}{\partial y^2} -\epsilon (1-\epsilon y)\frac{\partial v}{\partial y} -\epsilon ^2 v \right) \nonumber \\-\epsilon {\mathrm{Fr^{-1}}}\cos \theta (1-\epsilon y)^2,\end{eqnarray}
ε3(1εy)2vt+ε3(1εy)uvθ+ε3(1εy)2vvy+ε2(1εy)u2+(1εy)2py=ε2 Re f1ε22vθ2+2εuθ+(1εy)22vy2ε(1εy)vyε2vε Fr 1cosθ(1εy)2,
(B2)

and

\begin{eqnarray}\frac{\partial u}{\partial \theta } -\epsilon v +\frac{\partial v}{\partial y} -\epsilon y \frac{\partial v}{\partial y}=0.\end{eqnarray}
uθεv+vyεyvy=0.
(B3)

In terms of the above scaling, interface conditions (A19), (A20), and (10) at y = h are

\begin{eqnarray}(1-\epsilon h)^3(p-p_E) -{2\epsilon ^2}{\mathrm{Re_f}^{-1}}\left( \epsilon ^2\left(\frac{\partial h}{\partial \theta }\right)^2\left(\frac{\partial u}{\partial \theta }-\epsilon v\right) \right. \nonumber \\- \left.(1-\epsilon h)\left(\frac{\partial h}{\partial \theta }\right)\left((1-\epsilon h)\frac{\partial u}{\partial y}+ \epsilon u+\epsilon ^2\frac{\partial v}{\partial \theta }\right) +(1-\epsilon h)^3\frac{\partial v}{\partial y} \right)N^2 = \mathcal {N}(\theta ,t),\end{eqnarray}
(1εh)3(ppE)2ε2 Re f1ε2hθ2uθεv(1εh)hθ(1εh)uy+εu+ε2vθ+(1εh)3vyN2=N(θ,t),
(B4)
\begin{eqnarray}\left( -\epsilon ^2\left(\frac{\partial h}{\partial \theta }\right)^2\left((1-\epsilon h)\frac{\partial u}{\partial y}+\epsilon u+\epsilon ^2\frac{\partial v}{\partial \theta }\right) \right. \nonumber \\\left. + 2\epsilon ^2(1-\epsilon h)\left(\frac{\partial h}{\partial \theta }\right)\left(-\frac{\partial u}{\partial \theta }+\epsilon v+(1-\epsilon h)\frac{\partial v}{\partial y}\right)\right. \nonumber \\\left. + \left((1-\epsilon h)^3\frac{\partial u}{\partial y}+\epsilon u(1-\epsilon h)^2+\epsilon ^2(1-\epsilon h)^2\frac{\partial v}{\partial \theta }\right) \right)N^2 = \mathcal {T}(\theta ,t),\end{eqnarray}
ε2hθ2(1εh)uy+εu+ε2vθ+2ε2(1εh)hθuθ+εv+(1εh)vy+(1εh)3uy+εu(1εh)2+ε2(1εh)2vθN2=T(θ,t),
(B5)
\begin{eqnarray}{(1-\epsilon h)\frac{\partial h}{\partial t}+\left( u\frac{\partial h}{\partial \theta }-(1-\epsilon h) v\right) = \mathcal {M}(\theta ,t)},\end{eqnarray}
(1εh)ht+uhθ(1εh)v=M(θ,t),
(B6)

where

\[N=\left(1+\left(\frac{\epsilon }{1-\epsilon h}\frac{\partial h}{\partial \theta }\right)^2\right)^{-1/2}\]
N=1+ε1εhhθ21/2

and

\begin{eqnarray}\mathcal {N}=& \displaystyle { -(\epsilon \mathrm{We_f})^{-1} \left( (1-\epsilon h)^2 +2\epsilon ^2\left(\frac{\partial h}{\partial \theta }\right)^2 +\epsilon (1-\epsilon h)\frac{\partial ^2 h}{\partial \theta ^2} \right) N^3 +(1-\epsilon h)^3 (\epsilon \mathrm{We_f})^{-1} } \nonumber \\& \displaystyle {+ \frac{ a(1-\epsilon h)}{1-\epsilon ^k a}\left(-\epsilon ^{\frac{3-\beta }{2}} u_E\frac{\partial h}{\partial \theta }+(1-\epsilon h)\epsilon ^{\frac{1+\beta }{2}} v_E\right.} \nonumber \\& \displaystyle {\left.+\epsilon ^{\frac{3+k}{2}} \left(u\frac{\partial h}{\partial \theta }- (1-\epsilon h) v\right)\right)^2N^2,}\end{eqnarray}
N=(ε We f)1(1εh)2+2ε2hθ2+ε(1εh)2hθ2N3+(1εh)3(ε We f)1+a(1εh)1εkaε3β2uEhθ+(1εh)ε1+β2vE+ε3+k2uhθ(1εh)v2N2,
(B7)
\begin{eqnarray}\mathcal {T}&=& \displaystyle {-\mathrm{Re_f}\frac{a(1-\epsilon h)}{1-\epsilon ^k a}\left((1-\epsilon h)\epsilon ^{-\frac{\beta }{2}} u_E+\epsilon ^{\frac{\beta +2}{2}} v_E\frac{\partial h}{\partial \theta }- \epsilon ^{\frac{k}{2}}(1-\epsilon h)u- \epsilon ^{\frac{k+4}{2}} v\frac{\partial h}{\partial \theta }\right) } \nonumber \\[4pt]&&\qquad \displaystyle {\times \left(-\epsilon ^{\frac{2-\beta }{2}} u_E\frac{\partial h}{\partial \theta }+ (1-\epsilon h)\epsilon ^{\frac{\beta }{2}} v_E+ \epsilon ^{\frac{k+2}{2}}\left(u\frac{\partial h}{\partial \theta }- (1-\epsilon h) v\right)\right)},\end{eqnarray}
T= Re fa(1εh)1εka(1εh)εβ2uE+εβ+22vEhθεk2(1εh)uεk+42vhθ×ε2β2uEhθ+(1εh)εβ2vE+εk+22uhθ(1εh)v,
(B8)
\begin{eqnarray}\quad \mathcal {M} =\epsilon ^k a(1-\epsilon h)\frac{\partial h}{\partial t}+ a\left(\epsilon ^{\frac{k-\beta }{2}} u_E\frac{\partial h}{\partial \theta }-\epsilon ^{\frac{k+\beta -2}{2}} (1-\epsilon h)v_E\right).\end{eqnarray}
M=εka(1εh)ht+aεkβ2uEhθεk+β22(1εh)vE.
(B9)
1. Thin-film leading-order approximation

As ε → 0, asymptotic expansions for the film height h(θ, t), fluid velocity components u(y, θ, t), v(y, θ, t) and pressure p(y, θ, t) are sought

\begin{eqnarray}(u, \ v, \ p, \ h) = (u_0, \ v_0, \ p_0, \ h_0)+\epsilon (u_1, \ v_1, \ p_1, \ h_1) + O(\epsilon ^2).\end{eqnarray}
(u,v,p,h)=(u0,v0,p0,h0)+ε(u1,v1,p1,h1)+O(ε2).
(B10)

Substituting (B10) into (B1)–(B3) yields at leading order

\begin{eqnarray}{\mathrm{Re_f}^{-1}}\left(\frac{\partial ^2u_0}{\partial y^2}\right)-\frac{\partial p_0}{\partial \theta } ={\mathrm{Fr}^{-1}}\sin \theta ,\end{eqnarray}
Re f12u0y2p0θ= Fr 1sinθ,
(B11)
\begin{eqnarray}\hspace*{3.3pc} \frac{\partial p_0}{\partial y}=0,\end{eqnarray}
p0y=0,
(B12)
\begin{eqnarray}\quad \frac{\partial u_0}{\partial \theta }+\frac{\partial v_0}{\partial y}=0.\end{eqnarray}
u0θ+v0y=0.
(B13)

Expanding the non-homogeneous terms in (B4)–(B6),

\begin{eqnarray}\mathcal {N}(\theta ,t)=\mathcal {N}_{0}+\epsilon \mathcal {N}_{1}+O(\epsilon ^2),\end{eqnarray}
N(θ,t)=N0+εN1+O(ε2),
(B14)
\begin{eqnarray}\mathcal {T}(\theta ,t)=\mathcal {T}_{0}+\epsilon \mathcal {T}_{1}+O(\epsilon ^2),\end{eqnarray}
T(θ,t)=T0+εT1+O(ε2),
(B15)
\begin{eqnarray}\mathcal {M}(\theta ,t)=\mathcal {M}_{0}+\epsilon \mathcal {M}_{1}+O(\epsilon ^2),\end{eqnarray}
M(θ,t)=M0+εM1+O(ε2),
(B16)

where

$\mathcal {N}_{0}$
N0⁠, etc., are dependent on scalings yet to be fully determined, and substituting the expansions (B10) into the boundary conditions (B4)–(B6) and (11) yields the following leading-order conditions, on y = h0:

\begin{eqnarray}\hspace*{4.5pc} p_0-p_E= \mathcal {N}_{0},\end{eqnarray}
p0pE=N0,
(B17)
\begin{eqnarray}\hspace*{5.5pc} \frac{\partial u_0}{\partial y}=\mathcal {T}_0,\end{eqnarray}
u0y=T0,
(B18)
\begin{eqnarray}\frac{\partial h_0}{\partial t}+\left( u_0\frac{\partial h_0}{\partial \theta }- v_0\right)=\mathcal {M}_0,\end{eqnarray}
h0t+u0h0θv0=M0,
(B19)

and on y = 0,

\begin{eqnarray}u_0=u_{cyl},\quad v_0=v_{cyl}.\end{eqnarray}
u0=ucyl,v0=vcyl.
(B20)
2. Leading-order film-height equation

Direct integration of (B11)–(B13), subject to (B17), (B18), and (B20) leads to

\begin{equation}u_0=F_0\left(\frac{y^2}{2}-h_0y\right)+\mathcal {T}_0y+u_{cyl},\end{equation}
u0=F0y22h0y+T0y+ucyl,
(B21)

with

\begin{equation}F_0(\theta ,t)=\mathrm{Re_f}\left(\frac{\partial \mathcal {N}_{0}}{\partial \theta }+{\mathrm{Fr^{-1}}}\sin \theta \right),\end{equation}
F0(θ,t)= Re fN0θ+ Fr 1sinθ,
(B22)

and

\begin{equation}v_0 = -\left( \frac{\partial F_0}{\partial \theta }\frac{y^3}{6}+\frac{\partial }{\partial \theta }\left(-F_0h_0+\mathcal {T}_0\right)\frac{y^2}{2}\right)+v_{cyl}.\end{equation}
v0=F0θy36+θF0h0+T0y22+vcyl.
(B23)

On evaluating u0 and v0 at y = h0, the kinematic condition (B19) leads to the following film evolution equation:

\begin{eqnarray}\frac{\partial h_0}{\partial t} -\frac{\partial }{\partial \theta }\left(\mathrm{Re_f}\left(\frac{\partial \mathcal {N}_{0}}{\partial \theta }+{\mathrm{Fr^{-1}}}\sin \theta \right)\frac{h_0^3}{3}\right) +\frac{\partial }{\partial \theta }\left(\mathcal {T}_0\frac{h_0^2}{2}\right)+\frac{\partial }{\partial \theta }\left(u_{cyl}h_0\right)-v_{cyl}-\mathcal {M}_0=0, \nonumber \\\end{eqnarray}
h0tθ Re fN0θ+ Fr 1sinθh033+θT0h022+θucylh0vcylM0=0,
(B24)

where

$\mathcal {N}_{0}$
N0⁠,
$\mathcal {T}_0$
T0
, and
$\mathcal {M}_0$
M0
, the leading-order terms of (B7)–(B9) are dependent on specific choices of droplet parameters.

1.
K. J.
Ruschak
and
L. E.
Scriven
, “
Rimming flow of liquid in a rotating horizontal cylinder
,”
J. Fluid Mech.
76
,
113
125
(
1976
).
2.
O.
Menekse
,
J. V.
Wood
, and
D. S.
Riley
, “
Investigation of Fe2O3 − Al and Cr2O3 − Al reactions using high speed video camera
,”
Mat. Sci. Technol.
22
,
199
205
(
2006
).
3.
M. M.
Driscoll
,
C. S.
Stevens
, and
S. R.
Nagel
, “
Thin film formation during splashing of viscous liquids
,”
Phys. Rev. E
82
,
036303
1
(
2010
).
4.
M.
Farrall
, “
Numerical modelling of two-phase flow in a simplified bearing chamber
,” Ph.D. dissertation (
University of Nottingham
, Nottingham,
2000
).
5.
M.
Farrall
,
S.
Hibberd
,
K.
Simmons
, and
D.
Giddings
, “
Prediction of air/oil exit flows in a commercial aero-engine bearing chamber
,”
J. Aerosp. Eng.
220
,
197
202
(
2006
).
6.
M.
Farrall
,
S.
Hibberd
,
K.
Simmons
, and
P.
Gorse
, “
A numerical model for oil film flow in an aeroengine bearing chamber and comparison to experimental data
,”
Trans. ASME: J. Eng. Gas Turbines Power
128
,
111
117
(
2006
).
7.
P. D.
Hicks
and
R.
Purvis
, “
Air cushioning and bubble entrapment in three-dimensional droplet impacts
,”
J. Fluid Mech.
649
,
135
163
(
2010
).
8.
C. J.
Noakes
, “
The dynamics of liquid films on rotating surfaces
,” Ph.D. dissertation (
University of Nottingham
, Nottingham,
2001
).
9.
C. J.
Noakes
,
J. R.
King
, and
D. S.
Riley
, “
The effect of mass transfer on steady two-dimensional rimming flow
,”
J. Eng. Math.
71
,
223
236
(
2011
).
10.
M.
Villegas-Díaz
, “
Analytical and numerical studies of thin-film rimming flow subject to surface shear
,” Ph.D. dissertation (
University of Nottingham
, Nottingham,
2005
).
11.
M.
Villegas-Díaz
,
H.
Power
, and
D. S.
Riley
, “
On the stability of rimming flows to two-dimensional disturbances
,”
Fluid Dyn. Res.
33
,
141
172
(
2003
).
12.
M.
Villegas-Díaz
,
H.
Power
, and
D. S.
Riley
, “
Analytical and numerical studies of the stability of thin-film rimming flow subject to surface shear
,”
J. Fluid Mech.
541
,
317
344
(
2005
).
13.
S. T.
Thoroddsen
and
L.
Mahadevan
, “
Experimental study of coating flows in a partially-filled horizontally rotating cylinder
,”
Exp. Fluids
23
,
1
13
(
1997
).
14.
H. K.
Moffatt
, “
Behaviour of a viscous film on the outer surface of a rotating cylinder
,”
J. Méch.
16
,
651
673
(
1977
).
15.
V. V.
Pukhnachev
, “
Motion of a liquid film on a surface of a rotating cylinder in a gravitational field
,”
J. Appl. Mech. Tech. Phys.
18
,
344
351
(
1977
).
16.
R. E.
Johnson
, “
Steady-state coating flows inside a rotating horizontal cylinder
,”
J. Fluid Mech.
190
,
321
342
(
1988
).
17.
S. B. G.
O’Brien
and
E. G.
Gath
, “
The location of a shock in rimming flow
,”
Phys. Fluids
10
,
1040
1042
(
1998
).
18.
J. J.
Ashmore
,
A. E.
Hosoi
, and
H. A.
Stone
, “
The effect of surface tension on rimming flows in a partially filled rotating cylinder
,”
J. Fluid Mech.
479
,
65
98
(
2003
).
19.
A. E.
Hosoi
and
L.
Mahadevan
, “
Axial instability of a free-surface front in a partially filled horizontal rotating cylinder
,”
Phys. Fluids
11
,
97
106
(
1999
).
20.
M.
Tirumkudulu
and
A.
Acrivos
, “
Coating flows within a rotating horizontal cylinder: Lubrication analysis, numerical computations, and experimental measurements
,”
Phys. Fluids
13
,
14
19
(
2001
).
21.
S. K.
Wilson
,
R.
Hunt
, and
B. R.
Duffy
, “
On the critical solutions in coating and rimming flows on a uniformly rotating horizontal cylinder
,”
Q. J. Mech. Appl. Math.
55
,
357
387
(
2002
).
22.
A.
Acrivos
and
B.
Jin
, “
Rimming flows within a rotating horizontal cylinder: asymptotic analysis of the thin-film lubrication equations and stability of their solutions
,”
J. Eng. Math.
50
,
99
120
(
2004
).
23.
E. S.
Benilov
and
S. B. G.
O’Brien
, “
Inertial instability of a liquid film inside a rotating horizontal cylinder
,”
Phys. Fluids
17
,
1
16
(
2005
).
24.
C. J.
Noakes
,
J. R.
King
, and
D. S.
Riley
, “
On the development of rational approximations incorporating inertial effects in coating and rimming flows
,”
Q. J. Mech. Appl. Math.
59
,
163
190
(
2006
).
25.
M. A.
Kelmanson
, “
On inertial effects in the Moffatt-Pukhnachov coating-flow problem
,”
J. Fluid Mech.
633
,
327
353
(
2009
).
26.
E. J.
Hinch
and
M.
Kelmanson
, “
On the decay and drift of free-surface perturbations in viscous thin-film flow exterior to a rotating cylinder
,”
Proc. R. Soc. London
459
,
1193
1213
(
2003
).
27.
E. J.
Hinch
,
M.
Kelmanson
, and
P.
Metcalfe
, “
Shock-like free-surface perturbations in low-surface-tension, viscous, thin-film flow exterior to a rotating cylinder
,”
Proc. R. Soc. London
460
,
2975
2991
(
2004
).
28.
C. J.
Noakes
,
J. R.
King
, and
D. S.
Riley
, “
On three-dimensional stability of a uniform rigidly rotating film on a rotating cylinder
,”
Q. J. Mech. Appl. Math.
58
,
229
256
(
2005
).
29.
P.
Evans
,
L.
Schwartz
, and
R.
Roy
, “
Three-dimensional solutions for coating flow on a rotating horizontal cylinder: Theory and experiment
,”
Phys. Fluids
17
,
072102
(
2005
).
30.
E. S.
Benilov
,
S. J.
Chapman
,
J. B.
McLeod
,
J. R.
Ockendon
, and
V. S.
Zubkov
, “
On liquid films on an inclined plate
,”
J. Fluid Mech.
663
,
53
69
(
2010
).
31.
E. S.
Benilov
,
V. N.
Lapin
, and
S. B. G.
O’Brien
, “
On rimming flows with shocks
,”
J. Eng. Math.
(to be published).
32.
Y.
Lu
,
Z.
Liu
, and
T.
Xu
, “
Numerical simulation of aero-engine lubrication system
,”
J. Eng. Gas Turbines Power
131
,
034503
(
2009
).
33.
E. S.
Benilov
,
M. S.
Benilov
, and
S. B. G.
O’Brien
, “
Existence and stability of regularized shocks, with applications to rimming flows
,”
J. Eng. Math.
63
,
197
212
(
2009
).
34.
G. J. B.
Black
, “
Theoretical studies of thin-film flows
,” M.Phil. dissertation (
University of Strathclyde
, Glasgow,
2002
).
35.
D. Y.
Hsieh
and
S. P.
Ho
,
Wave and Stability in Fluids
(
World Scientific
,
Singapore
,
1994
).
36.
L. N.
Trefethen
, Spectral methods inMATLAB, (SIAM, Philadelphia,
2000
).
37.
E.
Momoniat
,
R.
Ravindran
, and
S.
Roy
, “
The influence of slot injection/suction on the spreading of a thin film flow under gravity and surface tension
,”
Acta Mech.
211
,
61
71
(
2010
).
38.
U.
Thiele
and
E.
Knobloch
, “
Thin liquid films on a slightly inclined heated plate
,”
Physica D
190
,
213
248
(
2004
).
39.
D. A.
Edwards
,
H.
Brenner
, and
D. T.
Wasan
,
Interfacial Transport Processes and Rheology
,
Butterworth-Heinemann Series in Chemical Engineering
(
Butterworth-Heinemann
,
Boston
,
1991
).
40.
J. C.
Slattery
,
Advanced Transport Phenomena
,
Cambridge Series in Chemical Engineering
(
Cambridge University Press
,
Cambridge, England
,
1999
).
41.
T. B.
Benjamin
,
W. G.
Pritchard
, and
S. J.
Tavener
, “
Steady and unsteady flows of a highly viscous liquid in a rotating horizontal cylinder
” (
1993
), preprint.