1 Introduction
The coating process is a means of depositing viscous liquid films of precise thickness and uniformity on a solid surface. Imperfections in the layer, especially as the speed of coating is increased or its thickness is decreased, may arise due to flow instabilities and various wetting and edge effects. That is why stability analyses of steady states in coating flows attracts industrial interest. In this work, we consider the problem of coating flows on a rotating sphere, focusing on the case when the centrifugal forces due to rotation compete with the gravitational draining, in the presence of surface tension. We examine both the transient dynamics of the film, and the various steady states that the system may attain, while maintaining an axisymmetric state.
The flow of a thin fluid film on a flat surface such as an inclined plane in the presence of gravity has been the subject of numerous investigations over the years; e.g. Benjamin (Reference Benjamin1957), Bankoff (Reference Bankoff1971), Hocking (Reference Hocking1990), Tuck & Schwartz (Reference Tuck and Schwartz1991), Giacomelli (Reference Giacomelli1999). Recently a new level of complexity has been added to the flat surface problem by coupling the flow to a surfactant spreading model (Shearer & Levy Reference Shearer and Levy2006; Levy, Shearer & Witelski Reference Levy, Shearer and Witelski2007; Peterson & Shearer Reference Peterson and Shearer2011). Including the effects of rotation, in a combined experimental and theoretical investigation, Parmar, Tirumkudulu & Hinch (Reference Parmar, Tirumkudulu and Hinch2009) have studied the problem of a viscous liquid film that coats a rotating vertical disk, obtaining the two-dimensional thickness profile that captures the effects of gravity, disk rotation and surface tension. Dynamics of viscous coating flows on inner and outer surfaces of a horizontally rotating cylinder have also been studied in a series of publications (Pukhnachev Reference Pukhnachev1977; Johnson Reference Johnson1988; O’Brien Reference O’Brien1998; Benilov, O’Brien & Sazonov Reference Benilov, O’Brien and Sazonov2003; Pukhnachov Reference Pukhnachov2005; Karabut Reference Karabut2007; Benilov, Benilov & Kopteva Reference Benilov, Benilov and Kopteva2008; Kelmanson Reference Kelmanson2009; Sahu & Kumar Reference Sahu and Kumar2014), after the pioneering result of Moffatt (Reference Moffatt1977) which was partially motivated by the process of keeping honey from spilling from a horizontal knife by rotating it. In that spirit, the present work can be regarded as a model for coating a rotating cherry with liquid chocolate. The problem of coating of a sphere with a thin viscous film has previously been studied by Takagi & Huppert (Reference Takagi and Huppert2010), albeit without consideration of surface tension and rotational effects.
Different regimes resulting in different steady states emerge due to the competition between the gravitational field that causes the downward draining of the liquid, the surface tension that tries to preserve a spherical shape for the coating layer and the centrifugal force from rotation about the vertical axis that tries to concentrate the liquid mass around the equator of the rotating sphere. Modelling of thin films on curved substrates would also have application to tear film dynamics, e.g. Li et al. (Reference Li, Braun, Maki, Henshaw and King-Smith2014).
In our work, we shall derive the lubrication equations in axisymmetric spherical coordinates. A previous derivation of these equations has been described by Myers, Charpin & Chapman (Reference Myers, Charpin and Chapman2002), who aimed to obtain a very general equation for thin films on an arbitrary two-dimensional substrate by using two principal coordinates along the substrate surface, together with the locally normal direction, thereby yielding a locally orthogonal curvilinear coordinate system. When their general equation is specialized to a spherical substrate (see their § III.D), the result is an equation that differs from ours in the placement of the various metric factors. We believe our derivation to be the proper one since we begin with the exact equations in spherical coordinates before applying the lubrication approximation.
The paper is structured as follows: in § 2, we describe the physical model and derive the governing nonlinear partial differential equation for the evolution of the film profile under the lubrication approximation. The model includes the effects of gravity, surface tension and centrifugal forces, but ignores those due to Coriolis forces. In § 3, we examine the energies associated with the three main effects that are considered, deriving the overall potential energy function that is minimized at equilibrium. In § 4, we consider the existence and stability of steady states and describe the numerical approach for solving the transient problem as the film approaches its energy-minimizing equilibrium. Section 5 provides a summary of the results and provides the outlook for possible extensions of the model. Some technical aspects of the lubrication approximation, together with the conditions under which Coriolis forces can be neglected while centrifugal effects are important, are provided in appendix A.
2 Model formulation
Consider the flow of a thin liquid film on a sphere of radius
$R$
under gravity
$\boldsymbol{g}$
, acting downward along the negative
$z$
-axis. The sphere is also allowed to rotate at a constant angular velocity
${\it\bf\Omega}$
about the
$z$
-axis and the equations of motion are obtained in the rotating frame of reference. Figure 1 provides a schematic drawing of this system. Let
$r$
denote the radial distance from the origin,
${\it\theta}$
the polar angle in spherical coordinates (with
${\it\theta}=0$
along the positive
$z$
-axis) and let
$h({\it\theta},t)$
represent the thickness of the thin film, with
$t$
denoting time. We assume the system to be axisymmetric and use the standard lubrication approximation. We start with the kinematic boundary condition and the incompressibility condition. The interface
$r=R+h({\it\theta},t)$
is the zero-level set of the function
$$\begin{eqnarray}\mathscr{F}(r,{\it\theta},t)=r-R-h({\it\theta},t)=0.\end{eqnarray}$$
The kinematic boundary condition,
$\text{D}\mathscr{F}/\text{D}t=0$
, then yields
$$\begin{eqnarray}\frac{\partial h}{\partial t}=v_{r}-\frac{v_{{\it\theta}}}{r}\frac{\partial h}{\partial {\it\theta}},\end{eqnarray}$$
at
$r=R+h({\it\theta},t)$
.
A thin liquid film under gravity on a sphere of radius
$R$
rotating about the
$z$
-axis.
Next take the continuity equation in spherical coordinates:
$$\begin{eqnarray}\frac{1}{r^{2}}\frac{\partial }{\partial r}(r^{2}v_{r})+\frac{1}{r\sin ({\it\theta})}\frac{\partial }{\partial {\it\theta}}(\sin ({\it\theta})v_{{\it\theta}})=0,\end{eqnarray}$$
multiply it by
$r^{2}\sin ({\it\theta})$
and integrate it at a fixed
${\it\theta}$
from
$r=R$
to
$r=R+h({\it\theta},t)$
. Upon applying the Leibniz rule to bring the
${\it\theta}$
-derivative outside the integral sign and the no-penetration condition
$v_{r}=0$
at
$r=R$
, we obtain
$$\begin{eqnarray}\sin ({\it\theta})\left[\left.(R+h)^{2}v_{r}\right|_{R+h}-\left.(R+h)v_{{\it\theta}}\right|_{R+h}\frac{\partial h}{\partial {\it\theta}}\right]+\frac{\partial }{\partial {\it\theta}}\left(\sin ({\it\theta})\int _{R}^{R+h}rv_{{\it\theta}}\,\text{d}r\right)=0.\end{eqnarray}$$
Using the kinematic boundary condition (2.2) reduces this equation to an evolution equation for
$h({\it\theta},t)$
, namely:
$$\begin{eqnarray}\frac{\partial h}{\partial t}+\frac{1}{(R+h)^{2}}\frac{1}{\sin ({\it\theta})}\frac{\partial }{\partial {\it\theta}}\left(\sin ({\it\theta})\int _{R}^{R+h}rv_{{\it\theta}}\,\text{d}r\right)=0.\end{eqnarray}$$
If the velocity component
$v_{{\it\theta}}$
can be related to
$h({\it\theta},t)$
via the lubrication approximation, the requisite integral in (2.5) can also be expressed in terms of
$h({\it\theta},t)$
, allowing one to track the subsequent evolution of the film thickness from the resulting nonlinear partial differential equation.
To that end, consider the
$r$
- and
${\it\theta}$
-components of Navier–Stokes equation under the lubrication approximation. Since we wish to include the effects of both gravity and centrifugal forces arising from the rotation of the sphere, it is convenient to introduce the modified pressure
$P$
, which is related to the regular pressure
$p$
in the thin film by
$$\begin{eqnarray}P=p-{\it\rho}\boldsymbol{g}\boldsymbol{\cdot }\boldsymbol{x}-{\textstyle \frac{1}{2}}{\it\rho}|{\it\bf\Omega}\times \boldsymbol{x}|^{2}=p+{\it\rho}gr\cos {\it\theta}-{\textstyle \frac{1}{2}}{\it\rho}{\it\Omega}^{2}r^{2}\sin ^{2}({\it\theta}).\end{eqnarray}$$
Here,
$\boldsymbol{x}$
is the position vector, gravity points in the negative
$z$
direction, so that
$\boldsymbol{g}=-g\hat{\boldsymbol{k}}$
and the angular velocity is about the
$z$
-axis, so that
${\it\bf\Omega}=\hat{\boldsymbol{k}}\,{\it\Omega}$
. The gradient of the modified pressure replaces the contributions from the regular pressure gradient and the body forces resulting from gravity and centrifugal forces in the rotating frame of reference. In terms of the modified pressure, the
$r$
- and
${\it\theta}$
-components of the momentum equation read
$$\begin{eqnarray}\frac{\partial P}{\partial r}=0,\end{eqnarray}$$
$$\begin{eqnarray}\frac{1}{r}\frac{\partial P}{\partial {\it\theta}}=\frac{{\it\mu}}{r^{2}}\frac{\partial }{\partial r}\left(r^{2}\frac{\partial v_{{\it\theta}}}{\partial r}\right).\end{eqnarray}$$
Note that as will be shown in appendix A, under the lubrication approximation, the inertial terms in the Navier–Stokes equation are negligible and it is possible to find a scaling where centrifugal forces are of the same order of magnitude as gravity and surface tension, while Coriolis forces remain negligible. As expected, in the thin (radial) direction, the pressure remains uniform in the film, while along the long direction (
${\it\theta}$
), the pressure gradient balances the dominant viscous term that involves the second derivative of
$v_{{\it\theta}}$
in
$r$
.
Equation (2.7) implies that the modified pressure is independent of
$r$
and may only depend on
${\it\theta}$
and
$t$
. The normal stress balance at
$r=R+h$
requires the pressure in the film to be equal to the uniform pressure
$p_{o}$
in the external air plus the capillary contribution associated with the surface tension
${\it\sigma}$
of the interface and its curvature
$\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\hat{\boldsymbol{n}}$
, i.e.
$p|_{r=R+h}=p_{o}+{\it\sigma}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\hat{\boldsymbol{n}}$
. Here,
$\hat{\boldsymbol{n}}=\boldsymbol{{\rm\nabla}}\mathscr{F}/|\boldsymbol{{\rm\nabla}}\mathscr{F}|$
(with
$\mathscr{F}$
given by (2.1)) is the unit normal at the interface pointing toward the exterior. When the film thickness
$h({\it\theta},t)$
is much smaller than the sphere radius
$R$
, the curvature can be approximated in a series in
$h$
, with the first two terms given by:
$$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\hat{\boldsymbol{n}}\approx \frac{2}{R}-\frac{1}{R^{2}}\left(2h+\frac{1}{\sin ({\it\theta})}\frac{\partial }{\partial {\it\theta}}\left(\sin ({\it\theta})\frac{\partial h}{\partial {\it\theta}}\right)\right).\end{eqnarray}$$
In terms of the modified pressure, the normal stress balance becomes:
$$\begin{eqnarray}P({\it\theta},t)=p_{o}+{\it\sigma}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\hat{\boldsymbol{n}}+{\it\rho}g(R+h)\cos ({\it\theta})-{\textstyle \frac{1}{2}}{\it\rho}{\it\Omega}^{2}(R+h)^{2}\sin ^{2}({\it\theta}).\end{eqnarray}$$
This fully determines the modified pressure throughout the film. Its
${\it\theta}$
-derivative provides the forcing term in the remaining component of the momentum equation. Since
$\partial P/\partial {\it\theta}$
is independent of
$r$
, one can easily integrate the
${\it\theta}$
-momentum equation to obtain
$$\begin{eqnarray}v_{{\it\theta}}=\frac{1}{2{\it\mu}}\frac{\partial P}{\partial {\it\theta}}r-\frac{C_{1}}{r}+C_{2}.\end{eqnarray}$$
To find
$C_{1}$
and
$C_{2}$
, we apply the no-slip boundary condition on the surface of the solid sphere and take the shear stress at the free surface to be zero. Namely, at
$r=R$
,
$v_{{\it\theta}}=0$
and at
$r=R+h$
,
$(\partial /\partial r)(v_{{\it\theta}}/r)\approx 0$
. Once the integration constants are obtained, the flow field is found to be:
$$\begin{eqnarray}v_{{\it\theta}}=\frac{1}{2{\it\mu}(R-h)r}\frac{\partial P}{\partial {\it\theta}}[(R-h)r^{2}+(R+h)R^{2}-2R^{2}r].\end{eqnarray}$$
In the evolution equation (2.5) for
$h({\it\theta},t)$
we need the integral
$\int _{R}^{R+h}rv_{{\it\theta}}\,\text{d}r$
. We evaluate this integral to find
$$\begin{eqnarray}\int _{R}^{R+h}rv_{{\it\theta}}\,\text{d}r=-\frac{1}{{\it\mu}}\frac{\partial P}{\partial {\it\theta}}\frac{(2R+h)h^{3}}{6(R-h)}.\end{eqnarray}$$
Substituting this result in the evolution equation for
$h({\it\theta},t)$
yields
$$\begin{eqnarray}\frac{\partial h}{\partial t}+\frac{1}{(R+h)^{2}}\frac{1}{\sin ({\it\theta})}\frac{\partial }{\partial {\it\theta}}\left[-\frac{\sin ({\it\theta})}{{\it\mu}}\frac{\partial P}{\partial {\it\theta}}\frac{(2R+h)h^{3}}{6(R-h)}\right]=0.\end{eqnarray}$$
From the expression (2.10) for the modified pressure
$P$
, the combination
$-(1/{\it\mu})(\partial P/\partial {\it\theta})$
appearing in the above evolution equation can be written as:
$$\begin{eqnarray}\displaystyle -\frac{1}{{\it\mu}}\frac{\partial P}{\partial {\it\theta}} & = & \displaystyle -\frac{{\it\rho}g}{{\it\mu}}\frac{\partial }{\partial {\it\theta}}\left[(R+h)\cos ({\it\theta})\right]+\frac{{\it\rho}{\it\Omega}^{2}}{2{\it\mu}}\frac{\partial }{\partial {\it\theta}}\left[(R+h)^{2}\sin ^{2}({\it\theta})\right]\nonumber\\ \displaystyle & & \displaystyle +\,\frac{{\it\sigma}}{{\it\mu}R^{2}}\frac{\partial }{\partial {\it\theta}}\left[2h+\frac{1}{\sin ({\it\theta})}\frac{\partial }{\partial {\it\theta}}\left(\sin ({\it\theta})\frac{\partial h}{\partial {\it\theta}}\right)\right].\end{eqnarray}$$
The combination of equations (2.14) and (2.15) is a fourth-order nonlinear PDE which fully determines the evolution of the film thickness under the influence of gravity, rotation and surface tension. We can simplify it somewhat since in the lubrication limit,
$h\ll R$
, some of the terms can be replaced by their leading-order approximations. This results in the equation:
$$\begin{eqnarray}\frac{\partial h}{\partial t}+\frac{1}{\sin ({\it\theta})}\frac{\partial }{\partial {\it\theta}}\left[\sin ({\it\theta})h^{3}\mathscr{Q}\right]=0,\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle \mathscr{Q}=\frac{{\it\rho}g}{3{\it\mu}R}\sin ({\it\theta})+\frac{{\it\rho}{\it\Omega}^{2}}{3{\it\mu}}\sin ({\it\theta})\cos ({\it\theta})+\frac{{\it\sigma}}{3{\it\mu}R^{4}}\frac{\partial }{\partial {\it\theta}}\left[2h+\frac{1}{\sin ({\it\theta})}\frac{\partial }{\partial {\it\theta}}\left(\sin ({\it\theta})\frac{\partial h}{\partial {\it\theta}}\right)\right]. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
2.1 Scaling and non-dimensionalization
In addition to the parameters that have already appeared above, let us define
$H$
to be the characteristic thickness of the liquid film. The lubrication approximation requires
$H\ll R$
, so we can define the small parameter
${\it\epsilon}=H/R$
for later scaling purposes. To render equations (2.16) and (2.17) dimensionless, choose a time scale
${\it\tau}$
(to be specified later) and define a dimensionless time variable
$\hat{t}=t/{\it\tau}$
and a dimensionless film thickness
${\hat{h}}=h/H$
. Substituting these into (2.16) suggests the dimensionless counterpart of function
$\mathscr{Q}$
should be defined by
$\hat{\mathscr{Q}}={\it\tau}H^{2}Q$
. Rewriting both equations (2.16) and (2.17) in dimensionless form and dropping the hats from all the variables for notational clarity then yields:
$$\begin{eqnarray}\frac{\partial h}{\partial t}+\frac{1}{\sin ({\it\theta})}\frac{\partial }{\partial {\it\theta}}\left[\sin ({\it\theta})h^{3}\mathscr{Q}\right]=0,\end{eqnarray}$$
$$\begin{eqnarray}\mathscr{Q}=a\,\sin ({\it\theta})+b\,\sin ({\it\theta})\cos ({\it\theta})+c\,\frac{\partial }{\partial {\it\theta}}\left[2h+\frac{1}{\sin ({\it\theta})}\frac{\partial }{\partial {\it\theta}}\left(\sin ({\it\theta})\frac{\partial h}{\partial {\it\theta}}\right)\right],\end{eqnarray}$$
where
$$\begin{eqnarray}a=\frac{{\it\tau}{\it\rho}gH^{2}}{3{\it\mu}R},\quad b=\frac{{\it\tau}{\it\rho}{\it\Omega}^{2}H^{2}}{3{\it\mu}},\quad c=\frac{{\it\tau}H^{3}{\it\sigma}}{3{\it\mu}R^{4}}.\end{eqnarray}$$
The no-flux boundary conditions that
$\sin ({\it\theta})h^{3}\mathscr{Q}=0$
at
${\it\theta}=0$
and
${\it\theta}={\rm\pi}$
lead to
$$\begin{eqnarray}\sin ({\it\theta})\frac{\partial h}{\partial {\it\theta}}=0\quad \text{and}\quad \sin ({\it\theta})\frac{\partial }{\partial {\it\theta}}\left(\frac{1}{\sin ({\it\theta})}\frac{\partial }{\partial {\it\theta}}\left(\sin ({\it\theta})\frac{\partial h}{\partial {\it\theta}}\right)\right)=0,\end{eqnarray}$$
at the two end points. Note that while
$\sin ({\it\theta})$
vanishes at these points, the no-flux conditions would, in principle, allow certain
${\it\theta}$
-derivatives of
$h$
to be singular at the boundaries as long as their products with
$\sin ({\it\theta})$
approach zero at the two end points.
Since we have assumed that all three effects, namely gravity, centrifugal force and surface tension, are of comparable magnitude and equally influence the dynamics, the three dimensionless constants
$a,b,c$
should all be of order unity; in practice, we also explore more extreme regimes when some of these parameters are small or large. We can always choose the time scale
${\it\tau}$
so as to make one of these three parameters equal to 1, and interpret the remaining two parameters as ratios of these effects. For instance, if we choose the time scale
${\it\tau}$
based on the gravitational drainage time, i.e. take
${\it\tau}=3{\it\mu}R/({\it\rho}gH^{2})$
, the parameters become
$$\begin{eqnarray}a=1,\quad b=\frac{{\it\Omega}^{2}R}{g},\quad c=\frac{{\it\sigma}H}{{\it\rho}gR^{3}}=\frac{{\it\epsilon}\,{\it\sigma}}{{\it\rho}gR^{2}}.\end{eqnarray}$$
Parameter
$b$
is a measure of the ratio of centrifugal force to gravity while parameter
$c$
is a ratio of surface tension forces to gravity (whose reciprocal is typically called the Bond number). We note, however, that in order for
$c$
to remain of order unity, surface tension may need to be relatively large, i.e.
${\it\sigma}=\mathit{O}({\it\rho}gR^{2}/{\it\epsilon})$
. That is indeed the distinguished limit that we are considering in this work.
At this point it is also useful to point out the scaling for the (modified) pressure
$P$
and the velocity field
$v_{{\it\theta}}$
. From expression (2.10) for the modified pressure, the non-constant part of the pressure scales as
${\it\rho}gR$
or
${\it\rho}{\it\Omega}^{2}R^{2}$
or
${\it\sigma}H/R^{2}$
, all three of which are of the same order of magnitude. Call this pressure scale
$P_{o}$
, e.g.
$P_{o}={\it\rho}gR$
. Note that the constant part of the pressure in (2.10),
$p_{o}+2{\it\sigma}/R$
, does not contribute to the dynamics so its order of magnitude is irrelevant. Then, expression (2.12) for
$v_{{\it\theta}}$
suggests that the velocity component parallel to the sphere surface should scale as
$V=H^{2}P_{o}/({\it\mu}R)$
. The easiest way to see this is to evaluate
$v_{{\it\theta}}$
at the free surface
$r=R+h$
, note that
$h$
scales with
$H$
, and consider the limit
$H\ll R$
. The relation between velocity scale
$V$
and pressure scale
$P_{o}$
can be rearranged to give
$P_{o}={\it\mu}RV/H^{2}={\it\mu}V/({\it\epsilon}^{2}R)$
. These are the same scalings that arise in standard lubrication theory, as will be shown in appendix A.
3 The energy functional
From a physical point of view, one might expect that any equilibrium solution to the thin-film equation should minimize the total energy of the system. We use the terms ‘equilibrium’ and ‘steady state’ interchangeably in this work. While in general it is possible to reach steady-state shape profiles within which one still has a non-zero flow field, in this work the steady states correspond to cases in which the flow field and the resultant flux are zero. We shall see that this is the case when we obtain the solution in the next section. However, it is useful to identify all the contributions to the potential energy of the system in advance, which is what we undertake in this section. There are three separate contributions to the overall potential energy in this system: gravitational energy
$E_{G}$
, rotational energy
$E_{R}$
and surface energy
$E_{S}$
.
The gravitational potential energy can be written as a volume integral of
${\it\rho}gz$
over the volume occupied by the fluid:
$$\begin{eqnarray}E_{G}=\int _{V}{\it\rho}gz\,\text{d}V={\it\rho}g\int _{0}^{2{\rm\pi}}\int _{0}^{{\rm\pi}}\int _{R}^{R+h}(r\cos {\it\theta})\,r^{2}\sin {\it\theta}\,\text{d}r\,\text{d}{\it\theta}\,\text{d}{\it\phi}.\end{eqnarray}$$
Similarly, the rotational energy is the volume integral of
$-(1/2){\it\rho}{\it\Omega}^{2}(r\sin ({\it\theta}))^{2}$
over the fluid volume:
$$\begin{eqnarray}E_{R}=-\frac{1}{2}{\it\rho}{\it\Omega}^{2}\int _{0}^{2{\rm\pi}}\int _{0}^{{\rm\pi}}\int _{R}^{R+h}(r\sin {\it\theta})^{2}r^{2}\sin {\it\theta}\text{d}r\,\text{d}{\it\theta}\,\text{d}{\it\phi}.\end{eqnarray}$$
Note that this energy is minimized when the mass distribution is as far away from the rotation axis as possible; hence the negative sign in front of the integral. Finally, the surface energy is the product of surface tension
${\it\sigma}$
(which can also be interpreted as the surface energy per unit area) with the total area of the free surface:
$$\begin{eqnarray}E_{S}={\it\sigma}\int _{S}\,\text{d}S=2{\rm\pi}{\it\sigma}\int _{0}^{{\rm\pi}}\sin {\it\theta}(R+h)\sqrt{(R+h)^{2}+(h_{{\it\theta}})^{2}}\,\text{d}{\it\theta}.\end{eqnarray}$$
While the total potential energy, consisting of the sum of the above three contributions, will reach a minimum at equilibrium, the total mass of the fluid must remain conserved, which is equivalent to the statement that
$$\begin{eqnarray}M={\it\rho}\int _{0}^{2{\rm\pi}}\int _{0}^{{\rm\pi}}\int _{R}^{R+h}r^{2}\sin {\it\theta}\,\text{d}r\,\text{d}{\it\theta}\,\text{d}{\it\phi}=\text{const.}\end{eqnarray}$$
To be able to quantify the gravitational, rotational and surface energies in terms of the
$\mathit{O}(1)$
dimensionless parameters
$a$
,
$b$
and
$c$
introduced in (2.20), we non-dimensionalize all three energies by dividing them by the energy scale:
$(6{\rm\pi}{\it\mu}R^{5})/({\it\tau}H^{2})$
. At the same time, we use the scaled form of the film thickness
${\hat{h}}=h/H$
, so that
$R+h=R(1+{\it\epsilon}{\hat{h}})$
with
${\it\epsilon}=H/R\ll 1$
as before. The total potential energy in the system is thus given by
$$\begin{eqnarray}\displaystyle \hat{E}=\frac{E_{G}+E_{R}+E_{S}}{(6{\rm\pi}{\it\mu}R^{5})/({\it\tau}H^{2})} & = & \displaystyle a\int _{0}^{{\rm\pi}}\frac{(1+{\it\epsilon}{\hat{h}})^{3}-1}{3}\cos {\it\theta}\sin {\it\theta}\,\text{d}{\it\theta}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{b}{2}\int _{0}^{{\rm\pi}}\frac{(1+{\it\epsilon}{\hat{h}})^{5}-1}{5}\sin ^{3}{\it\theta}\,\text{d}{\it\theta}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{c}{{\it\epsilon}}\int _{0}^{{\rm\pi}}(1+{\it\epsilon}{\hat{h}})\sqrt{(1+{\it\epsilon}{\hat{h}})^{2}+{\it\epsilon}^{2}\left(\frac{\partial {\hat{h}}}{\partial {\it\theta}}\right)^{2}}\sin {\it\theta}\,\text{d}{\it\theta}.\end{eqnarray}$$
Expanding all the terms on the right-hand side to
$\mathit{O}({\it\epsilon})$
results in:
$$\begin{eqnarray}\displaystyle \hat{E} & = & \displaystyle a\int _{0}^{{\rm\pi}}{\it\epsilon}\,{\hat{h}}\cos {\it\theta}\sin {\it\theta}\,\text{d}{\it\theta}-\frac{b}{2}\int _{0}^{{\rm\pi}}{\it\epsilon}\,{\hat{h}}\sin ^{3}{\it\theta}\,\text{d}{\it\theta}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{c}{{\it\epsilon}}\int _{0}^{{\rm\pi}}\left(1+2\,{\it\epsilon}\,{\hat{h}}+{\it\epsilon}^{2}{\hat{h}}^{2}+\frac{1}{2}{\it\epsilon}^{2}\left(\frac{\partial {\hat{h}}}{\partial {\it\theta}}\right)^{2}\right)\sin {\it\theta}\,\text{d}{\it\theta}+\mathit{O}({\it\epsilon}^{2}).\end{eqnarray}$$
The potential energy is only determined to within an arbitrary additive constant, therefore we may omit any part in
$\hat{E}$
that is constant. The requirement of mass conservation (3.4) given above, can be expressed in powers of
${\it\epsilon}$
, the first two terms of which result in
$$\begin{eqnarray}\int _{0}^{{\rm\pi}}({\it\epsilon}{\hat{h}}+{\it\epsilon}^{2}{\hat{h}}^{2})\sin {\it\theta}\,\text{d}{\it\theta}=\text{const.}\end{eqnarray}$$
Making use of this result to cancel the constant parts in
$\hat{E}$
allows us to express the non-constant part of
$\hat{E}$
as:
$$\begin{eqnarray}\hat{E}={\it\epsilon}\int _{0}^{{\rm\pi}}\left[a\,{\hat{h}}\cos {\it\theta}+\frac{b}{2}{\hat{h}}\cos ^{2}{\it\theta}+c\left(\frac{1}{2}\left(\frac{\partial {\hat{h}}}{\partial {\it\theta}}\right)^{2}-{\hat{h}}^{2}\right)\right]\sin {\it\theta}\,\text{d}{\it\theta}+\mathit{O}({\it\epsilon}^{2}).\end{eqnarray}$$
Under the change of variable
$x=-\!\cos {\it\theta}$
, the leading term in the potential energy functional can therefore be expressed as:
$$\begin{eqnarray}\frac{\hat{E}}{{\it\epsilon}}=\int _{-1}^{1}\left[-ax{\hat{h}}+\frac{b}{2}x^{2}{\hat{h}}+c\left(\frac{1}{2}(1-x^{2})\left(\frac{\partial {\hat{h}}}{\partial x}\right)^{2}-{\hat{h}}^{2}\right)\right]\text{d}x.\end{eqnarray}$$
As we shall see in the next section, this energy functional decreases monotonically as the film thickness evolves in time until the energy achieves its minimum value at equilibrium.
4 Transient and steady-state computations
The computations can be simplified by removing the singularity in
$1/\sin ({\it\theta})$
when
${\it\theta}$
is 0 or
${\rm\pi}$
in (2.18)–(2.19). Upon introducing the change of variable
$x=-\!\cos {\it\theta}$
with
${\it\theta}\in (0,{\rm\pi})$
, so that
$x\in (-1,1)$
, the evolution equation transforms into
$$\begin{eqnarray}\frac{\partial h}{\partial t}+\frac{\partial }{\partial x}\left[h^{3}(1-x^{2})\left(a-bx+c\frac{\partial }{\partial x}\left[2h+\frac{\partial }{\partial x}\left((1-x^{2})\frac{\partial h}{\partial x}\right)\right]\right)\right]=0,\end{eqnarray}$$
with dimensionless parameters
$a$
,
$b$
and
$c$
given earlier in (2.20).
The no-flux boundary conditions at the end points
$x=\pm 1$
take the form (no boundary conditions would be required if the solution is classical as
$(1-x^{2})$
vanishes at both poles):
$$\begin{eqnarray}\displaystyle & \displaystyle (1-x^{2})\frac{\partial h}{\partial x}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle (1-x^{2})\frac{\partial ^{2}}{\partial x^{2}}\left((1-x^{2})\frac{\partial h}{\partial x}\right)=0. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\int _{-1}^{1}h(x,t)\,\text{d}x=0.\end{eqnarray}$$
4.1 Transient computations
To perform numerical computations, the fourth-order equation (4.1) was replaced by a system of two second-order equations as follows. Define the function
$v(x,t)$
by
$$\begin{eqnarray}v(x,t)=\frac{\partial }{\partial x}\left((1-x^{2})\frac{\partial h}{\partial x}\right)\end{eqnarray}$$
in terms of which (4.1) becomes:
$$\begin{eqnarray}\frac{\partial h}{\partial t}+\frac{\partial }{\partial x}\left[h^{3}(1-x^{2})\left(a-bx+c\frac{\partial }{\partial x}(2h+v)\right)\right]=0.\end{eqnarray}$$
The appropriate boundary conditions can then be expressed as:
$$\begin{eqnarray}\left.\begin{array}{@{}l@{}}\displaystyle (1-x^{2})\partial h/\partial x=0\quad \text{at }x=\pm 1,\\ \displaystyle (1-x^{2})\partial v/\partial x=0\quad \text{at }x=\pm 1.\end{array}\right\}\end{eqnarray}$$
This system was integrated numerically using the finite element solver COMSOL Multiphysics. We used the module for one-dimensional General Form PDE with Direct UMFPACK Solver and our discretized system contained approximately twenty-five thousand degrees of freedom. For the initial condition we mainly used symmetric perturbations of the uniform distribution
$h(x,0)=h_{0}$
and
$v(x,0)=0$
with different values of the initial layer thickness
$h_{0}$
. (Note, however, that it is possible to scale
$h_{0}$
out of the formulation, taking it to be 1, by appropriately modifying the parameters
$a$
,
$b$
and
$c$
.) Test simulations with decreasing mesh sizes and relative or absolute tolerances were performed to insure that the results converged to the same steady-state solutions and thus were robust. We will illustrate the convergence to a steady state by displaying the asymptotic approach of the energy to a constant.
4.2 Energy-minimizing steady states
Multiplying (4.1) by the expression
$$\begin{eqnarray}{\it\varphi}(h)=\left(ax-\frac{bx^{2}}{2}\right)+c\left[2h+\frac{\partial }{\partial x}\left((1-x^{2})\frac{\partial h}{\partial x}\right)\right]\end{eqnarray}$$
and integrating the result over the interval
$(-1,1)$
, we obtain the energy dissipation
$$\begin{eqnarray}\frac{\text{d}}{\text{d}t}E(h)=-\int _{-1}^{1}h^{3}(1-x^{2})\left\{a-bx+c\left[2\frac{\partial h}{\partial x}+\frac{\partial ^{2}}{\partial x^{2}}\left((1-x^{2})\frac{\partial h}{\partial x}\right)\right]\right\}^{2}\text{d}x\leqslant 0,\end{eqnarray}$$
where the energy functional
$E(h)$
is given by:
$$\begin{eqnarray}E(h)=\int _{-1}^{1}\left[c\left(\frac{1}{2}(1-x^{2})\left(\frac{\partial h}{\partial x}\right)^{2}-h^{2}\right)+\left(\frac{bx^{2}}{2}-ax\right)h\right]\text{d}x,\end{eqnarray}$$
which is the same as (3.9). It should be noted that in terms of
${\it\varphi}$
, the evolution equation has the ‘gradient flow’ structure
$$\begin{eqnarray}\frac{\partial h}{\partial t}+\frac{\partial }{\partial x}\left[h^{3}(1-x^{2})\frac{\partial {\it\varphi}}{\partial x}\right]=0,\end{eqnarray}$$
while for the same
${\it\varphi}$
:
$$\begin{eqnarray}\int _{-1}^{1}{\it\varphi}\frac{\partial h}{\partial t}\,\text{d}x=-\frac{\text{d}E}{\text{d}t},\end{eqnarray}$$
leading to the energy-dissipation equation
$$\begin{eqnarray}\frac{\text{d}E}{\text{d}t}=-\int _{-1}^{1}h^{3}(1-x^{2})\left(\frac{\partial {\it\varphi}}{\partial x}\right)^{2}\,\text{d}x\leqslant 0.\end{eqnarray}$$
We now establish that the energy functional
$E(h)$
has a finite lower bound under the constraints
$$\begin{eqnarray}h\geqslant 0\quad \text{and}\quad \int _{-1}^{1}h(x)\,\text{d}x=M.\end{eqnarray}$$
Note that the dimensionless ‘volume’
$M$
on the right-hand side of the last equation is not an independent parameter and its value is related to the characteristic thickness
$H$
of the film that was introduced earlier. In particular, if
$H$
is chosen to be the average of the initial thickness profile, the dimensionless volume turns out to be
$M=2$
. It proves convenient to split the energy functional into the sum
$E(h)=E_{1}(h)+E_{2}(h)$
. The first term
$E_{1}(h)$
is taken to be
$$\begin{eqnarray}E_{1}(h)=\int _{-1}^{1}c\,\left(\frac{1}{2}(1-x^{2})\left(\frac{\partial h}{\partial x}\right)^{2}-h^{2}\right)\text{d}x=c\,(\mathscr{L}h,h)\end{eqnarray}$$
in which
$({\it\phi},{\it\psi})=\int _{-1}^{1}{\it\phi}(x){\it\psi}(x)\,\text{d}x$
is the standard inner product and the operator
$\mathscr{L}$
is related to the second-order Legendre differential operator and is defined by
$$\begin{eqnarray}\mathscr{L}h=-{\textstyle \frac{1}{2}}((1-x^{2})h^{\prime })^{\prime }-h.\end{eqnarray}$$
The eigenfunctions of this operator are the Legendre polynomials
$P_{\ell }(x)$
with corresponding eigenvalues
$\ell (\ell +1)/2-1$
with
$\ell =0,1,2,\ldots$
. The first few eigenvalues of
$\mathscr{L}$
are given by
$-1,0,2,5,9,\ldots$
. It follows from the mass conservation constraint
$\int _{-1}^{1}h(x)\,\text{d}x=M$
that
$\bar{h}=M/2$
. Let
$h(x)=\bar{h}+h_{1}(x)$
so that
$h_{1}(x)$
has a zero mean, i.e.
$h_{1}$
is orthogonal to any constant. Any series expansion of
$h_{1}(x)$
would thus not include the term corresponding to
$\ell =0$
and will start from the
$\ell =1$
term. Therefore,
$(\mathscr{L}h_{1},h_{1})\geqslant 0$
. As such,
$$\begin{eqnarray}(\mathscr{L}h,h)=(\mathscr{L}(h_{1}+M/2),(h_{1}+M/2))=(\mathscr{L}h_{1},h_{1})-\frac{M^{2}}{2}\geqslant -\frac{M^{2}}{2}.\end{eqnarray}$$
At the same time
$$\begin{eqnarray}E_{2}(h)=\int _{-1}^{1}\left(\frac{1}{2}bx^{2}-ax\right)\,h(x)\,\text{d}x\geqslant -\frac{a^{2}}{2b}M,\end{eqnarray}$$
which derives from the fact that the parabola
$(bx^{2}/2-ax)$
has the minimum
$-a^{2}/2b$
together with the non-negativity constraint on
$h(x)$
. As a result
$$\begin{eqnarray}E(h)\geqslant -\frac{c}{2}M^{2}-\frac{a^{2}}{2b}M.\end{eqnarray}$$
Since the energy functional
$E(h)(t)$
is decaying, bounded from below and lower semi-continuous it must have a minimizer,
$h(x)$
, which is continuous on
$(-1,1)$
(Burchard, Chugunova & Stephens Reference Burchard, Chugunova and Stephens2012, Lemma 2.1).
Minimizing the energy functional
$E(h)$
under the constraints (4.14) and solving the corresponding Euler–Lagrange equation
$$\begin{eqnarray}-c\left[\frac{\text{d}}{\text{d}x}\left((1-x^{2})\frac{\text{d}h}{\text{d}x}\right)+2h\right]-\left(ax-\frac{bx^{2}}{2}\right)={\it\lambda}\end{eqnarray}$$
gives as an explicit
${\it\lambda}$
-parameter dependent expression, obtained by Maple (available from Maplesoft), for the energy-minimizing steady state
$$\begin{eqnarray}h(x)=c_{1}x+c_{2}\left(x\ln \frac{1-x}{1+x}+2\right)+\frac{a\,x\ln (1-x^{2})}{6c}+\frac{b}{8c}\left(1-x^{2}\right)-\frac{{\it\lambda}}{2c},\end{eqnarray}$$
which is valid only for uniformly positive parts of the energy minimizer.
We can also show that if the touchdown point
$x_{1}\in (-1,1)$
exists where
$h(x_{1})=0$
, the first derivative of
$h$
must also be zero at that point. Let
$h^{\ast }(x)\geqslant 0$
with
$\!\int _{-1}^{1}h^{\ast }\,\text{d}x=M$
be the energy-minimizing steady state. Define
$$\begin{eqnarray}h^{{\it\epsilon}}=\frac{M}{M^{{\it\epsilon}}}(h^{\ast }+{\it\epsilon}\,{\it\phi}),\quad M^{{\it\epsilon}}=\int _{-1}^{1}\,(h^{\ast }+{\it\epsilon}\,{\it\phi})\,\text{d}x,\end{eqnarray}$$
where
${\it\phi}(x)$
is chosen such that
$h^{{\it\epsilon}}$
remains non-negative. We compute the first variation of the energy functional as
$$\begin{eqnarray}\frac{\text{d}}{\text{d}{\it\epsilon}}E(h^{{\it\epsilon}})|_{{\it\epsilon}=0}=\int _{-1}^{1}\,\left[c(1-x^{2})\frac{\text{d}h^{\ast }}{\text{d}x}\frac{\text{d}{\it\phi}}{\text{d}x}-2ch^{\ast }{\it\phi}-\left(ax-\frac{bx^{2}}{2}\right)\,{\it\phi}-{\it\lambda}\,{\it\phi}\right]\text{d}x,\end{eqnarray}$$
where
$$\begin{eqnarray}{\it\lambda}=\frac{1}{M}\left(2E(h^{\ast })+\int _{-1}^{1}\,\left(ax-\frac{bx^{2}}{2}\right)\,h^{\ast }\,\text{d}x\right).\end{eqnarray}$$
The requirement that the first variation of the energy functional should be non-negative for any admissible test function means that
$h^{\ast }$
must satisfy the following variational inequality:
$$\begin{eqnarray}c\left[\frac{\text{d}}{\text{d}x}\left((1-x^{2})\frac{\text{d}h^{\ast }}{\text{d}x}\right)+2h^{\ast }\right]+\left(ax-\frac{bx^{2}}{2}\right)\leqslant ~-{\it\lambda}.\end{eqnarray}$$
This inequality implies that there exists a positive constant
$C$
such that in the neighbourhood of a touchdown point
$x_{1}$
, the inequality
$0\leqslant h^{\ast }(x)\leqslant C\,(x-x_{1})^{2}$
holds (Burchard et al.
Reference Burchard, Chugunova and Stephens2012, Lemma 2.2). Hence
$h^{\ast }$
has a zero first derivative (i.e. a zero contact angle) at the touchdown point and, as a result, it is continuously differentiable on
$(-1,1)$
. The zero contact angle for the energy-minimizing steady state corresponds to the full wetting regime. To obtain minimizers with non-zero contact angles, one needs to modify the energy functional by adding terms that account for the energies of the solid–liquid and solid–gas surfaces; see for example Benilov & Oron (Reference Benilov and Oron2010).
(a) Film shape at steady state in the absence of rotation (
$a=1,\,b=0,\,c=1$
). (b) Uniformly positive steady state in the absence of gravity (
$a=0,\,b=1,\,c=1$
). The radius of the sphere is chosen for visualization only.
The energy-minimizing steady states must always have at least one dry zone if the gravitational term is present, i.e. for
$a>0$
, because in this case the non-negativity constraint is active. Without this constraint, the energy functional is not bounded from below and will tend to negative infinity along a family of linear functions of the form:
$h_{{\it\alpha}}=M/2+{\it\alpha}x$
corresponding to which
$E(h_{{\it\alpha}})=Mb/6-Mc-(2a/3)\,{\it\alpha}\rightarrow -\infty$
as
${\it\alpha}\rightarrow +\infty$
.
4.2.1 Uniformly positive steady states
We consider the special case when the energy-minimizing steady state is strictly positive:
$h(x)>0$
(e.g. see figure 2
b). Hence the gravity constant
$a$
and the coefficient
$c_{2}$
in (4.21) must both be zero. Value of the constant
$c_{1}$
is undetermined in the absence of the gravitational field. With an additional assumption that
$h(x)$
is symmetric with respect to
$x=0$
, we will have
$c_{1}=0$
and the expression for the steady state simplifies to
$$\begin{eqnarray}h(x)=\frac{b}{8c}(1-x^{2})-\frac{{\it\lambda}}{2c}.\end{eqnarray}$$
From the volume constraint
$\int _{-1}^{1}h(x)\,\text{d}x=M$
, we obtain that
$$\begin{eqnarray}{\it\lambda}=\frac{b-6Mc}{6}.\end{eqnarray}$$
To make sure that
$h(x)>0$
for all
$x$
and in particular at
$x=\pm 1$
, it follows that the total volume
$M$
must satisfy the inequality
$$\begin{eqnarray}M>\frac{b}{6c}\quad \text{or}\quad H>\frac{b}{12c}.\end{eqnarray}$$
The faster the rotation, the thicker the original uniform coating needs to be in order not to develop a second dry zone on the bottom of the sphere. It also follows from (4.26) that for parameter values
$a=0,\,\,c=1,\,\,M=2$
, the second dry zone appears when
$b\geqslant 6M$
or
$b\geqslant 12$
.
4.2.2 Steady states with one dry zone
Assume that
$x_{1}$
is the only touchdown point of
$h(x)$
(see figure 2
a for an example of a steady state with one dry zone). The energy-minimizing steady state
$h(x)$
is a continuous function on the closed finite interval. Hence it must be bounded
$h(x)<\infty$
. The steady state
$h(x)$
also satisfies the conditions:
$$\begin{eqnarray}h(x_{1})=0,\quad h^{\prime }(x_{1})=0,\quad \int _{x_{1}}^{1}h(x)\,\text{d}x=M.\end{eqnarray}$$
It follows from
$h(x)<\infty$
that
$\lim _{x\rightarrow 1}h(x)<\infty$
; hence
$$\begin{eqnarray}c_{2}\ln (1-x^{2})+\frac{a}{6c}\ln (1-x^{2})=0.\end{eqnarray}$$
As such
$c_{2}=-a/6c$
and the expression for the positive part of
$h(x)$
simplifies to
$$\begin{eqnarray}h(x)=\frac{a}{3c}\left(x\ln (1+x)-1\right)+\frac{b}{8c}\left(1-x^{2}\right)+c_{1}x-\frac{{\it\lambda}}{2c}.\end{eqnarray}$$
To solve for
$x_{1}$
,
${\it\lambda}$
and
$c_{1}$
we apply the constraints (4.29) to this form of the function
$h(x)$
and obtain a system of three nonlinear equations that are not displayed explicitly here for brevity. The solution to this system can be obtained by Maple. For instance, for the set of the parameters
$a=b=c=1$
and
$M=0.2$
, we obtain the following result:
$c_{1}\approx 0.221,\,\,{\it\lambda}\approx -0.502,\,\,x_{1}\approx -0.336$
. The energy value at the steady state is approximately
$E(h)\approx -0.089$
. Convergence of the energy to this minimal value, as obtained by performing the transient computations in COMSOL for this set of parameters, is illustrated in figure 4.
4.2.3 Steady states with two dry zones
For a fixed volume
$M$
and with gravity present, i.e. for
$a>0$
, as the speed of rotation (i.e. parameter
$b$
) increases, the energy-minimizing steady state first loses convexity at the bottom, (i.e. the derivative
$h^{\prime }(1)$
changes sign, going from negative to positive, so that at that point
$h^{\prime }(1)=0$
) and, after that, if the rotation rate continues to increase, a second dry zone appears at the bottom of the sphere (with
$h(1)=0$
at and above the threshold). This transition in shape is illustrated in different coordinate systems in figures 3 and 6.
We can map out the boundaries in parameter space where these transitions take place. For instance, if we fix the strengths of gravity and surface tension by taking
$a=c=1$
, but allow ‘mass’
$M$
to vary, we can identify how fast the rotation (parameter
$b$
) needs to be for a dimple to develop at the bottom of the interface, and to obtain a second dry zone at the bottom of the sphere. Figure 7(a) shows these two threshold values for the parameter
$b$
as a function of
$M$
. Below the dashed curve, which is almost linear for large
$M$
having a slope
$3/2$
, the steady state is a convex-shaped drop hanging from the bottom of the sphere; between the dashed curve and the solid curve, which is almost linear for large
$M$
having a slope
$6$
, the steady state has a dimple, with the dimple depth growing as
$b$
increases toward the solid curve, and above the solid curve the steady-state shape has two dry zones, one above and one below the rotating sphere. Figure 7(b) shows the nonlinear behaviour of both curves for extremely small values of
$M$
. It should be born in mind that the first touchdown point
$x_{1}$
also varies, moving toward
$-1$
(the top of the sphere), as
$M$
increases.
Change in the steady-state film profile upon increasing the rate of rotation. The parameters are
$a=1,\,b=5,\,c=1$
on (a) and
$a=1,\,b=20,\,c=1$
on (b).
Convergence of the energy to a constant value for the parameters
$a=1$
,
$b=1$
,
$c=1$
and
$M=0.2$
.
Convergence of the energy to a constant value for the parameters
$a=0$
,
$b=100$
,
$c=1$
and
$M=0.2$
.
Alternatively, we can also fix the ‘mass’
$M$
and surface tension parameter
$c=1$
and examine how these same transitions take place as the ratio
$b/a$
of rotational to gravitational forces increases. Figure 8(b) illustrates the dependence of the ratio
$b/a$
on the value of the touchdown point
$x_{1}$
for a fixed mass. The dashed line corresponds the loss of convexity and the solid line to the appearance of the second dry zone. The zoomed picture on the right side illustrates that the steady state has a convex shape and only one dry zone on the top of the sphere if
$c=1$
and
$b<a$
for any mass
$M$
. Note that for a fixed
$c$
, once a steady solution
$h(x)$
is found for a given set of parameters
$\{a,b,M\}$
, a scaled solution
$kh(x)$
also exists for the scaled parameters
$\{ka,kb,kM\}$
. We can thus fix the mass (e.g. at
$M=2$
) and only the ratio
$b/a$
determines the transitions in the shapes.
Consider the case when two dry zones are present (e.g. see figure 3
b), one on the top and one on the bottom of the rotating sphere. In this case the steady state has two touchdown points
$x_{1}$
and
$x_{2}$
. The set of equations that determine their positions consists of:
$$\begin{eqnarray}h(x_{1})=0,\quad h(x_{2})=0,\quad h^{\prime }(x_{1})=0,\quad h^{\prime }(x_{2})=0,\quad \int _{x_{1}}^{x_{2}}h(x)\,\text{d}x=M,\end{eqnarray}$$
together with the requirement that
$h$
should be regular or
$h(x)<\infty$
. Given,
$a$
,
$b$
,
$c$
and
$M$
, the resulting equations for
$x_{1}$
,
$x_{2}$
,
$c_{1}$
,
$c_{2}$
and
${\it\lambda}$
are nonlinear and lengthy (they are not displayed here for brevity) and challenging to solve.
Shape evolution of the steady state which leads to the development of the second dry zone due to increase in the speed of the rotation in Cartesian coordinates. Parameters (a–c) right
$(a=1,b=3,c=1)\,$
,
$(a=1,b=5,c=1)\,$
and
$(a=1,b=10,c=1)$
.
For
$a=c=1$
, (a) the dashed curve corresponds to the value of
$b(M)$
when the steady state loses convexity and the solid curve corresponds to the value
$b(M)$
when a second dry zone appears at the bottom of the sphere; (b) the zoom-in illustrates the nonlinear behaviour of both curves when
$M$
approaches zero.
For
$c=1$
, (a) the ratio
$b/a$
versus the value of the first touchdown point
$x_{1}$
; the dashed curve corresponds to the value of the ratio
$b/a$
when the steady state loses convexity and the solid curve to the value
$b/a$
when the second dry zone starts to appear; (b) the zoom-in illustrates that for
$b<a$
the steady state has a convex shape and only one dry zone on the top of the sphere for any value of
$M$
.
However, for the symmetric case with
$x_{1}=-x_{2}$
which corresponds to the absence of gravity
$a=0$
, Maple can find the solution to the nonlinear equations numerically. For instance, for parameters
$a=0$
,
$b=100$
,
$c=1$
and
$M=0.2$
, the numerical solution is
$c_{1}=0$
,
$c_{2}\approx -5.051$
,
${\it\lambda}\approx 4.091$
and
$x_{2}=-x_{1}\approx 0.518$
. The energy value at this steady state is approximately
$E(h)\approx 0.608$
. The asymptotic approach of the energy to this minimal value, obtained via transient computations in COMSOL, is illustrated in figure 5.
More generally, to obtain non-symmetric steady states with 2 dry zones, one can proceed as follows. First we note that parameter
$c$
can always be set to unity (provided surface tension is present) by taking the time scale
${\it\tau}$
to be given by
$3{\it\mu}R^{4}/H^{3}{\it\sigma}$
in (2.20). With that choice for the time scale, parameters
$a$
and
$b$
both become proportional to
$1/H$
. Therefore, as the value of
$H$
changes, the ratio of
$a$
to
$b$
does not change. Also, adjustments to
$H$
are equivalent to changing the volume
$M$
in (4.32). We thus proceed by setting
$b=1$
as well (to be adjusted later). Now, rather than specifying the parameters and looking for the touchdown points, we start by specifying the touchdown points and see if we can find the parameters that correspond to them. So, we choose the two values
$x_{1}$
and
$x_{2}$
for the touchdown points in the range
$-1<x_{1}<x_{2}<1$
. Next we impose three of the constraints
$$\begin{eqnarray}h^{\prime }(x_{1})=0,\quad h^{\prime }(x_{2})=0,\quad h(x_{1})=h(x_{2}).\end{eqnarray}$$
The constant
${\it\lambda}$
disappears from these three equations. But with
$c=1$
and
$b=1$
, the remaining coefficients
$a$
,
$c_{1}$
and
$c_{2}$
can be found from the three equations, which are linear in the three unknowns. At this point, constant
${\it\lambda}$
which simply shifts the profile
$h(x)$
up and down can be found by enforcing that
$h(x_{1})=0$
or
$h(x_{2})=0$
. The portion of the resulting curve
$h(x)$
that lies between
$x_{1}$
and
$x_{2}$
represents the film profile in the wet zone. To adjust its area to any desired value, we simply calculate the present value of the integral
$M^{\prime }=\int _{x_{1}}^{x_{2}}h(x)\,\text{d}x$
and multiply the function
$h(x)$
by
$M/M^{\prime }$
. This has the effect of adjusting the values of parameters
$a$
and
$b$
(and
${\it\lambda}$
,
$c_{1}$
and
$c_{2}$
), while keeping
$c=1$
. For definiteness, we can choose
$M=2$
, which corresponds to the interpretation of
$H$
as the average of the initial height profile. Our scaled numerical simulations can thus start with the uniform initial condition
$h(x,0)=1$
.
As an example, for the arbitrarily chosen case with touchdown points at
$x_{1}=0.25$
and
$x_{2}=0.75$
, the above procedure leads to the parameter values
$a=16\,084.8$
,
$b=31\,785.5$
and
$c=1$
corresponding to
$M=2$
. Transient COMSOL simulations starting with
$h(x,0)=1$
for these parameter values are shown in figure 9, verifying that the solution approaches the predicted steady profile having those same touchdown points for long time periods.
Time evolution of the coating shape in
$x$
,
${\it\theta}$
and polar coordinates is presented by two time snap shots
$t=2\times 10^{-5}$
on (a,c,e)
$t=5\times 10^{3}$
on (b,d,f).
5 Discussion
In this work, we have derived the evolution equation for a thin liquid film on a rotating sphere in the presence of gravity, surface tension and centrifugal forces. We have also obtained the overall potential energy functional for this system and shown that the energy-minimizing steady states of the system are global attractors.
Upon analysing the steady states in various parameter regimes, we can summarize the results as follows. In the absence of rotation but with gravity and surface tension both present, the steady state has a convex drop-like shape with a dry zone on the top of the sphere. If we parameterize this state in terms of the total volume of the liquid film, we find that as the volume decreases, the dry zone will increase in size and the wet region occupied by the drop becomes more concentrated near the bottom (i.e. south pole) of the sphere. By contrast, as the amount of liquid in the system increases, the dry zone ultimately disappears and the state becomes one with a uniformly positive film thickness all around the sphere. Adding rotation (or increasing dimensionless parameter
$b$
away from zero) changes the shape of the steady state in the following way: starting with the convex-shaped hanging drop with a dry zone on the top of the sphere, as the rotation rate
$b$
increases, the drop shape loses convexity and a dimple appears in the interface at the bottom. As the rotation rate increases even more, that dimple grows in size and eventually touches down at the south pole, initiating a second dry zone that grows in size with further increases of the rotation rate.
The steady states on such a rotating sphere are qualitatively different from the ones that have been studied on a rotating cylinder or vertical disk. In most of these (an exception being the work of Evans, Schwartz & Roy (Reference Evans, Schwartz and Roy2004) who showed that the centrifugal force can lead to instabilities and drop formation on a rapidly rotating cylinder), centrifugal forces are not taken into account. Rather, the kinematic effect of rotation along with the no-slip boundary condition provide an equilibrium film thickness with the gravitational draining flow counteracted by the upward motion of the surface on one half of the region, and assisted by the downward motion on the other side. In our analysis, in the rotating frame, the sphere surface is completely stationary. But the centrifugal forces due to rotation drive the fluid toward the equator. The resulting steady states are shapes with no flow occurring within them once equilibrium is established.
There are multiple avenues for extending and generalizing the present work. It would be interesting to relax the assumption of axisymmetry and allow for dependence on the azimuthal angle. In particular, whether the axisymmetric global attractors obtained in the present work are stable to azimuthal perturbations is an important question to address. It is also useful to study cases where Marangoni effects associated with the dependence of surface tension on temperature or surfactant concentration are taken into account. One can also study thin films on the inside of rotating spheres, and consider cases where the effects of Coriolis forces are not negligible.
Acknowledgements
This work was partially supported by a grant from the Simons Foundation (no. 275088 to M.C.).
Appendix A. Justification of the lubrication approach
Whereas § 2 derives the thin-film equation starting with the simplified dimensional forms of the momentum and continuity equations, here we present the more rigorous derivation starting with the complete equations of motion and boundary conditions, with proper scaling and non-dimensionalization of the variables resulting in the correct leading-order evolution equation for the film, showing it to be the same as the one obtained in § 2.
The governing incompressible Navier–Stokes equations describing flow on the outer surface of the rotating sphere in the rotating frame are given by
$$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v}=0,\end{eqnarray}$$
$$\begin{eqnarray}{\it\rho}\frac{\text{D}\boldsymbol{v}}{\text{D}t}=-\boldsymbol{{\rm\nabla}}p+{\it\rho}\boldsymbol{g}+{\it\mu}{\rm\nabla}^{2}\boldsymbol{v}-2{\it\rho}\,{\it\bf\Omega}\times \boldsymbol{v}-{\it\rho}\,{\it\bf\Omega}\times ({\it\Omega}\times \boldsymbol{r}).\end{eqnarray}$$
The last two terms in the momentum equation represent the Coriolis and centrifugal forces, respectively. We assume that the sphere is rotating with a constant angular velocity
${\it\bf\Omega}$
. In this paper, we are interested in the limit where centrifugal and gravity forces are of comparable magnitudes but the Coriolis forces are negligible in comparison. The parameter regime where this is possible is described below and it is only when the Coriolis term is negligible that one can obtain a two-dimensional axisymmetric flow (with only
$r$
- and
${\it\theta}$
-components) in the rotating frame.
The no-slip boundary condition on the sphere surface in the rotating frame is given by
$$\begin{eqnarray}\boldsymbol{v}=\mathbf{0}\quad \text{at }r=R.\end{eqnarray}$$
We also need to introduce the dynamical boundary conditions (i.e. the normal and tangential stress balances) on the free surface of the film, namely
$$\begin{eqnarray}-\hat{\boldsymbol{n}}\boldsymbol{\cdot }\unicode[STIX]{x1D640}\boldsymbol{\cdot }\hat{\boldsymbol{n}}=p_{o}+{\it\sigma}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\hat{\boldsymbol{n}},\quad \boldsymbol{t}\boldsymbol{\cdot }\unicode[STIX]{x1D640}\boldsymbol{\cdot }\hat{\boldsymbol{n}}=0,\quad \text{at }r=R+h({\it\theta},t),\end{eqnarray}$$
where, assuming no dependence on the azimuthal angle
${\it\phi}$
,
$$\begin{eqnarray}\hat{\boldsymbol{n}}=\frac{\boldsymbol{{\rm\nabla}}\mathscr{F}}{|\boldsymbol{{\rm\nabla}}\mathscr{F}|}=\frac{1}{\sqrt{1+\left(h_{{\it\theta}}/(R+h)\right)^{2}}}\left[\begin{array}{@{}c@{}}1\\ h_{{\it\theta}}/(R+h)\end{array}\right]\end{eqnarray}$$
is the unit normal vector,
$$\begin{eqnarray}\boldsymbol{t}=\left[\begin{array}{@{}c@{}}h_{{\it\theta}}/(R+h)\\ -1\end{array}\right]\end{eqnarray}$$
is the tangent vector and
$$\begin{eqnarray}\unicode[STIX]{x1D640}=\left[\begin{array}{@{}cc@{}}\displaystyle 2{\it\mu}\frac{\partial v_{r}}{\partial r}-p & \displaystyle {\it\mu}\left(r\frac{\partial }{\partial r}\left(\frac{v_{{\it\theta}}}{r}\right)+\frac{1}{r}\frac{\partial v_{r}}{\partial {\it\theta}}\right)\\ \displaystyle {\it\mu}\left(r\frac{\partial }{\partial r}\left(\frac{v_{{\it\theta}}}{r}\right)+\frac{1}{r}\frac{\partial v_{r}}{\partial {\it\theta}}\right) & \displaystyle 2{\it\mu}\left(\frac{1}{r}\frac{\partial v_{{\it\theta}}}{\partial {\it\theta}}+\frac{v_{r}}{r}\right)-p\end{array}\right]\end{eqnarray}$$
is the stress tensor in the viscous film. Here,
$h_{{\it\theta}}\equiv \partial h/\partial {\it\theta}$
. Therefore we can rewrite the normal stress balance as
$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{-2{\it\mu}}{1+\left({\displaystyle \frac{h_{{\it\theta}}}{R+h}}\right)^{2}}\left[\frac{\partial v_{r}}{\partial r}+\frac{h_{{\it\theta}}r}{R+h}\frac{\partial }{\partial r}\left(\frac{v_{{\it\theta}}}{r}\right)+\frac{h_{{\it\theta}}}{(R+h)r}\frac{\partial v_{r}}{\partial {\it\theta}}+\left(\frac{h_{{\it\theta}}}{R+h}\right)^{2}\left(\frac{1}{r}\frac{\partial v_{{\it\theta}}}{\partial {\it\theta}}+\frac{v_{r}}{r}\right)\right]\nonumber\\ \displaystyle & & \displaystyle \quad =p_{o}-p|_{R+h}+{\it\sigma}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{n}\end{eqnarray}$$
and the tangential stress balance as
$$\begin{eqnarray}\displaystyle & & \displaystyle \left(\frac{h_{{\it\theta}}}{R+h}\right)^{2}\left(r\frac{\partial }{\partial r}\left(\frac{v_{{\it\theta}}}{r}\right)+\frac{1}{r}\frac{\partial v_{r}}{\partial {\it\theta}}\right)+\frac{2h_{{\it\theta}}}{R+h}\left(\frac{\partial v_{r}}{\partial r}-\frac{v_{r}}{r}-\frac{1}{r}\frac{\partial v_{{\it\theta}}}{\partial {\it\theta}}\right)\nonumber\\ \displaystyle & & \displaystyle \quad -\,r\frac{\partial }{\partial r}\left(\frac{v_{{\it\theta}}}{r}\right)-\frac{1}{r}\frac{\partial v_{r}}{\partial {\it\theta}}=0,\end{eqnarray}$$
at
$r=R+h$
.
We also need the continuity equation in axisymmetric spherical coordinates:
$$\begin{eqnarray}\frac{1}{r^{2}}\frac{\partial }{\partial r}(r^{2}v_{r})+\frac{1}{r\sin ({\it\theta})}\frac{\partial }{\partial {\it\theta}}\left(\sin ({\it\theta})v_{{\it\theta}}\right)=0.\end{eqnarray}$$
In spherical coordinates, with
${\it\bf\Omega}={\it\Omega}\hat{\boldsymbol{k}}$
and with
$\boldsymbol{g}=-g\hat{\boldsymbol{k}}$
, we can define a modified pressure
$P$
by
$$\begin{eqnarray}P=p-{\textstyle \frac{1}{2}}{\it\rho}{\it\Omega}^{2}r^{2}\sin ^{2}({\it\theta})+{\it\rho}gr\cos ({\it\theta}),\end{eqnarray}$$
using which the momentum equation can be written simply as
$$\begin{eqnarray}{\it\rho}\frac{\text{D}\boldsymbol{v}}{\text{D}t}=-\boldsymbol{{\rm\nabla}}P+{\it\mu}{\rm\nabla}^{2}\boldsymbol{v},\end{eqnarray}$$
having dropped the Coriolis term (the justification for which is presented below). The
$r$
-component of the momentum equation reads:
$$\begin{eqnarray}\displaystyle & & \displaystyle {\it\rho}\left(\frac{\partial v_{r}}{\partial t}+v_{r}\frac{\partial v_{r}}{\partial r}+\frac{v_{{\it\theta}}}{r}\frac{\partial v_{r}}{\partial {\it\theta}}-\frac{v_{{\it\theta}}^{2}}{r}\right)=-\frac{\partial P}{\partial r}\nonumber\\ \displaystyle & & \displaystyle \quad +\,{\it\mu}\left[\frac{1}{r^{2}}\frac{\partial ^{2}}{\partial r^{2}}(r^{2}v_{r})+\frac{1}{r^{2}\sin {\it\theta}}\frac{\partial }{\partial {\it\theta}}\left(\sin {\it\theta}\frac{\partial v_{r}}{\partial {\it\theta}}\right)\right],\end{eqnarray}$$
while its
${\it\theta}$
-component is given by
$$\begin{eqnarray}\displaystyle & & \displaystyle {\it\rho}\left(\frac{\partial v_{{\it\theta}}}{\partial t}+v_{r}\frac{\partial v_{{\it\theta}}}{\partial r}+\frac{v_{{\it\theta}}}{r}\frac{\partial v_{{\it\theta}}}{\partial {\it\theta}}+\frac{v_{r}v_{{\it\theta}}}{r}\right)=-\frac{1}{r}\frac{\partial P}{\partial {\it\theta}}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,{\it\mu}\left[\frac{1}{r^{2}}\frac{\partial }{\partial r}\left(r^{2}\frac{\partial v_{{\it\theta}}}{\partial r}\right)+\frac{1}{r^{2}}\frac{\partial }{\partial {\it\theta}}\left(\frac{1}{\sin {\it\theta}}\frac{\partial }{\partial {\it\theta}}(v_{{\it\theta}}\sin {\it\theta})\right)+\frac{2}{r^{2}}\frac{\partial v_{r}}{\partial {\it\theta}}\right].\end{eqnarray}$$
We define the scaled non-dimensional variables:
$$\begin{eqnarray}\hat{v}_{{\it\theta}}=\frac{v_{{\it\theta}}}{V},\quad \hat{v}_{r}=\frac{v_{r}}{{\it\epsilon}V},\quad \hat{{\it\theta}}={\it\theta},\quad {\hat{y}}=\frac{r-R}{{\it\epsilon}R},\quad \hat{t}=\frac{Vt}{R}\quad \text{and}\quad \hat{P}=\frac{P}{P_{0}},\end{eqnarray}$$
where
$$\begin{eqnarray}V=\frac{{\it\rho}gH^{2}}{{\it\mu}}=\frac{{\it\epsilon}^{2}P_{0}R}{{\it\mu}}\quad \text{and}\quad P_{0}={\it\rho}gR.\end{eqnarray}$$
In terms of these variables, the non-dimensional form of the continuity equation (A 10) becomes:
$$\begin{eqnarray}\frac{{\it\epsilon}VR^{2}}{{\it\epsilon}R^{3}}\frac{1}{(1+{\it\epsilon}{\hat{y}})^{2}}\frac{\partial \left((1+{\it\epsilon}{\hat{y}})^{2}\hat{v}_{r}\right)}{\partial {\hat{y}}}+\frac{V}{R\sin \hat{{\it\theta}}}\frac{1}{1+{\it\epsilon}{\hat{y}}}\frac{\partial }{\partial \hat{{\it\theta}}}\left(\hat{v}_{{\it\theta}}\sin \hat{{\it\theta}}\right)=0.\end{eqnarray}$$
The
$r$
-component of the momentum equation reads:
$$\begin{eqnarray}\displaystyle & & \displaystyle \hspace{-12.0pt}\frac{V^{2}}{gR}\left({\it\epsilon}\frac{\partial \hat{v}_{r}}{\partial \hat{t}}+{\it\epsilon}\hat{v}_{r}\frac{\partial \hat{v}_{r}}{\partial {\hat{y}}}+{\it\epsilon}\frac{\hat{v}_{{\it\theta}}}{1+{\it\epsilon}{\hat{y}}}\frac{\partial \hat{v}_{r}}{\partial {\it\theta}}-\frac{\hat{v}_{{\it\theta}}^{2}}{1+{\it\epsilon}{\hat{y}}}\right)=-\frac{1}{{\it\epsilon}}\frac{\partial \hat{P}}{\partial {\hat{y}}}\nonumber\\ \displaystyle & & \displaystyle \hspace{-12.0pt}\qquad +\,\frac{{\it\mu}V}{{\it\rho}gR^{2}}\left[\frac{1}{{\it\epsilon}(1+{\it\epsilon}{\hat{y}})^{2}}\frac{\partial ^{2}}{\partial {\hat{y}}^{2}}((1+{\it\epsilon}{\hat{y}})^{2}\hat{v}_{r})+\frac{{\it\epsilon}}{(1+{\it\epsilon}{\hat{y}})^{2}\sin \hat{{\it\theta}}}\frac{\partial }{\partial \hat{{\it\theta}}}\left(\sin \hat{{\it\theta}}\frac{\partial \hat{v}_{r}}{\partial \hat{{\it\theta}}}\right)\right]\qquad\end{eqnarray}$$
while its
${\it\theta}$
-component is given by
$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{V^{2}}{gR}\left(\frac{\partial \hat{v}_{{\it\theta}}}{\partial \hat{t}}+\hat{v}_{r}\frac{\partial \hat{v}_{{\it\theta}}}{\partial {\hat{y}}}+\frac{\hat{v}_{{\it\theta}}}{1+{\it\epsilon}{\hat{y}}}\frac{\partial \hat{v}_{{\it\theta}}}{\partial \hat{{\it\theta}}}+\frac{\hat{v}_{r}\hat{v}_{{\it\theta}}}{1+{\it\epsilon}{\hat{y}}}\right)\nonumber\\ \displaystyle & & \displaystyle \quad =-\frac{1}{1+{\it\epsilon}{\hat{y}}}\frac{\partial \hat{P}}{\partial \hat{{\it\theta}}}+\frac{{\it\mu}V}{{\it\rho}gR^{2}}\left[\frac{1}{{\it\epsilon}^{2}(1+{\it\epsilon}{\hat{y}})^{2}}\frac{\partial }{\partial {\hat{y}}}\left((1+{\it\epsilon}{\hat{y}})^{2}\frac{\partial \hat{v}_{{\it\theta}}}{\partial {\hat{y}}}\right)\right.\nonumber\\ \displaystyle & & \displaystyle \left.\qquad +\,\frac{1}{(1+{\it\epsilon}{\hat{y}})^{2}}\frac{\partial }{\partial \hat{{\it\theta}}}\left(\frac{1}{\sin \hat{{\it\theta}}}\frac{\partial }{\partial \hat{{\it\theta}}}(\hat{v}_{{\it\theta}}\sin \hat{{\it\theta}})\right)+\frac{2}{(1+{\it\epsilon}{\hat{y}})^{2}}\frac{\partial \hat{v}_{r}}{\partial \hat{{\it\theta}}}\right].\end{eqnarray}$$
We use the lubrication approximation parameter
${\it\epsilon}=H/R\ll 1$
. Therefore we have
$V={\it\epsilon}^{2}P_{0}R/{\it\mu}\sim \mathit{O}({\it\epsilon}^{2})$
and
$\mathit{Re}={\it\rho}VR/{\it\mu}\sim \mathit{O}({\it\epsilon}^{2}).$
The dimensionless form of the no-slip boundary condition is
$$\begin{eqnarray}\hat{v}_{{\it\theta}}=0\quad \text{at }{\hat{y}}=0.\end{eqnarray}$$
After scaling the tangential stress balance boundary condition (A 9), we see that the
$r(\partial /\partial r)(v_{{\it\theta}}/r)$
term is dominant; thus this boundary condition can be simplified to
$$\begin{eqnarray}\frac{\partial \hat{v}_{{\it\theta}}}{\partial {\hat{y}}}\approx 0\quad \text{at }{\hat{y}}={\hat{h}}.\end{eqnarray}$$
In the normal stress balance, the terms on the left-hand side of (A 8) have scale no larger than
${\it\mu}V/R$
while the right-hand side has the scale
$P_{0}$
. So the right-hand side dominates and this boundary condition reduces to
$$\begin{eqnarray}\hat{P}|_{{\hat{y}}={\hat{h}}}=\frac{p_{o}}{{\it\rho}gR}-\frac{{\it\Omega}^{2}R}{2g}\sin ^{2}({\it\theta})+\cos ({\it\theta})+\frac{{\it\sigma}}{{\it\rho}gR}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\hat{\boldsymbol{n}}.\end{eqnarray}$$
The leading-order form of the continuity equation is
$$\begin{eqnarray}\frac{\partial \hat{v}_{r}}{\partial {\hat{y}}}+\frac{1}{\sin \hat{{\it\theta}}}\frac{\partial }{\partial \hat{{\it\theta}}}\left(\hat{v}_{{\it\theta}}\sin \hat{{\it\theta}}\right)=0.\end{eqnarray}$$
As
$V^{2}/(gR)={\it\epsilon}^{2}\mathit{Re}=\mathit{O}({\it\epsilon}^{4})$
, the leading-order form of the
$r$
-momentum equation reads:
$$\begin{eqnarray}\frac{\partial \hat{P}}{\partial {\hat{y}}}=0.\end{eqnarray}$$
Finally, the leading-order form of the
${\it\theta}$
-momentum equation becomes:
$$\begin{eqnarray}\frac{\partial \hat{P}}{\partial \hat{{\it\theta}}}=\frac{\partial ^{2}\hat{v}_{{\it\theta}}}{\partial {\hat{y}}^{2}}.\end{eqnarray}$$
With regard to the neglect of the Coriolis term, we see that it scales with
${\it\rho}{\it\Omega}V$
, which is negligible compared to the leading-order term
$\partial P/\partial r$
(with scale
${\it\rho}g/{\it\epsilon}$
) in the
$r$
-momentum equation and the leading-order term
$\partial P/\partial {\it\theta}$
(with scale
${\it\rho}gR$
) in the
${\it\theta}$
-momentum equation. In particular, ensuring that the Coriolis term is negligible compared to gravity means
${\it\Omega}\ll g/V$
, while keeping the Coriolis term small compared to the centrifugal term requires
${\it\Omega}\gg V/R$
. The simultaneous inequalities
$$\begin{eqnarray}\frac{V}{R}\ll {\it\Omega}\ll \frac{g}{V}\end{eqnarray}$$
can only be achieved if
$V$
itself satisfies the inequality
$$\begin{eqnarray}V^{2}\ll gR.\end{eqnarray}$$
Since the drainage velocity scale
$V={\it\rho}gH^{2}/{\it\mu}$
can be quite small for a thin and highly viscous film, it can be much smaller than the velocity scale
$\sqrt{gR}$
associated with free fall through distance
$R$
. In terms of this drainage velocity scale, the above requirement becomes:
${\it\epsilon}^{2}{\it\rho}\sqrt{gR^{3}}/{\it\mu}\ll 1$
, and the angular velocity
${\it\Omega}$
needs to be large compared to
$V/R$
but small compared to
$g/V$
.
We now can solve for
$\hat{v}_{{\it\theta}}$
from the
${\it\theta}$
-momentum equation and (A 20), (A 21):
$$\begin{eqnarray}\hat{v}_{{\it\theta}}=\frac{1}{2}\frac{\partial \hat{P}}{\partial \hat{{\it\theta}}}{\hat{y}}({\hat{y}}-2{\hat{h}}).\end{eqnarray}$$
The evolution equation for the thin film, obtained by integrating the continuity equation at a fixed
${\it\theta}$
across the radial direction and applying the kinematic boundary condition at the surface, reads
$$\begin{eqnarray}\frac{\partial {\hat{h}}}{\partial \hat{t}}=-\frac{1}{\sin \hat{{\it\theta}}}\frac{\partial }{\partial \hat{{\it\theta}}}\left(\sin \hat{{\it\theta}}\int _{0}^{{\hat{h}}}\hat{v}_{{\it\theta}}\,\text{d}{\hat{y}}\right).\end{eqnarray}$$
By substituting the solution for
$\hat{v}_{{\it\theta}}$
we find
$$\begin{eqnarray}\frac{\partial {\hat{h}}}{\partial \hat{t}}+\frac{1}{\sin \hat{{\it\theta}}}\frac{\partial }{\partial \hat{{\it\theta}}}\left(-\frac{\sin \hat{{\it\theta}}}{3}\frac{\partial \hat{P}}{\partial \hat{{\it\theta}}}{\hat{h}}^{3}\right)=0.\end{eqnarray}$$
To obtain the expression for
$\partial \hat{P}/\partial \hat{{\it\theta}}$
, we use the normal stress balance condition at the boundary (A 22)
$$\begin{eqnarray}\hat{P}|_{{\hat{y}}={\hat{h}}}=\frac{p_{o}}{{\it\rho}gR}-\frac{{\it\Omega}^{2}R}{2g}\sin ^{2}({\it\theta})+\cos ({\it\theta})+\frac{{\it\sigma}}{{\it\rho}gR}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\hat{\boldsymbol{n}},\end{eqnarray}$$
where
$$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\hat{\boldsymbol{n}}\approx \frac{2}{R}-\frac{1}{R^{2}}\left(2h+\frac{1}{\sin ({\it\theta})}\frac{\partial }{\partial {\it\theta}}\left(\sin ({\it\theta})\frac{\partial h}{\partial {\it\theta}}\right)\right).\end{eqnarray}$$
We thus find the leading-order dimensionless pressure gradient to be (dropping the hats for clarity):
$$\begin{eqnarray}-\frac{\partial P}{\partial {\it\theta}}=\sin ({\it\theta})+\frac{{\it\Omega}^{2}R}{g}\sin ({\it\theta})\cos ({\it\theta})+\frac{{\it\sigma}H}{{\it\rho}gR^{3}}\frac{\partial }{\partial {\it\theta}}\left[2h+\frac{1}{\sin ({\it\theta})}\frac{\partial }{\partial {\it\theta}}\left(\sin ({\it\theta})\frac{\partial h}{\partial {\it\theta}}\right)\right],\end{eqnarray}$$
which leads to the leading-order dimensionless evolution equation:
$$\begin{eqnarray}\frac{\partial h}{\partial t}+\frac{1}{\sin ({\it\theta})}\frac{\partial }{\partial {\it\theta}}\left[\sin ({\it\theta})h^{3}\mathscr{Q}\right]=0,\end{eqnarray}$$
with
$$\begin{eqnarray}\mathscr{Q}=a\,\sin ({\it\theta})+b\,\sin ({\it\theta})\cos ({\it\theta})+c\,\frac{\partial }{\partial {\it\theta}}\left[2h+\frac{1}{\sin ({\it\theta})}\frac{\partial }{\partial {\it\theta}}\left(\sin ({\it\theta})\frac{\partial h}{\partial {\it\theta}}\right)\right].\end{eqnarray}$$
where
$$\begin{eqnarray}a=\frac{{\it\tau}{\it\rho}gH^{2}}{3{\it\mu}R},\quad b=\frac{{\it\tau}{\it\rho}{\it\Omega}^{2}H^{2}}{3{\it\mu}},\quad c=\frac{{\it\tau}H^{3}{\it\sigma}}{3{\it\mu}R^{4}},\end{eqnarray}$$
in agreement with the result of § 2. Note that in arriving at this form, we have allowed for a more arbitrary scaling of the time variable with time scale
${\it\tau}$
as opposed to the original
$R/V$
.
Appendix B. Navier slip boundary condition
In this paper we restricted our analysis and computations to the no-slip condition between the surface of the sphere and the coating liquid layer. It is relatively straightforward to model some degree of slip at that surface. We shall present the derivation here, working with the dimensional form of the equations. The boundary condition at
$r=R$
changes from no slip (
$v_{{\it\theta}}=0$
) to one allowing for slip based on the Navier slip model, namely
$$\begin{eqnarray}v_{{\it\theta}}=\frac{{\it\beta}}{{\it\mu}}{\it\tau}_{r{\it\theta}}={\it\beta}R\frac{\partial }{\partial r}\left(\frac{v_{{\it\theta}}}{r}\right)\quad \text{at }r=R.\end{eqnarray}$$
Here,
${\it\tau}_{r{\it\theta}}$
is the shear stress evaluated on the sphere surface. Note that
${\it\beta}$
has units of length and can be interpreted as a ‘slip length’. We recover the no-slip condition when
${\it\beta}=0$
.
Solving the equation for
$v_{{\it\theta}}$
with the new boundary condition given above results in
$$\begin{eqnarray}v_{{\it\theta}}=\left(\frac{R^{4}+R^{3}(h-2r)+R^{2}r^{2}-Rhr^{2}-2{\it\beta}hr^{2}}{2(R^{2}-Rh-2{\it\beta}h)r}\right)\,\frac{1}{{\it\mu}}\frac{\partial P}{\partial {\it\theta}}.\end{eqnarray}$$
Calculating
$q=\int _{R}^{R+h}rv_{{\it\theta}}\,\text{d}r$
then yields
$$\begin{eqnarray}q=\frac{h}{6}\left(\frac{3R^{4}}{2{\it\beta}h-R(R-h)}+3R^{2}+3Rh+h^{2}\right)\,\frac{1}{{\it\mu}}\frac{\partial P}{\partial {\it\theta}}.\end{eqnarray}$$
When
${\it\beta}=0$
, we recover the no-slip result that
$$\begin{eqnarray}q=-\frac{h^{3}(2R+h)}{6(R-h)}\,\frac{1}{{\it\mu}}\frac{\partial P}{\partial {\it\theta}}\approx -\frac{h^{3}}{3}\frac{1}{{\it\mu}}\frac{\partial P}{\partial {\it\theta}}\quad \text{for }h\ll R.\end{eqnarray}$$
If we consider the case where
${\it\beta}$
and
$h$
are of the same order of magnitude and both are much smaller than
$R$
(i.e. with
${\it\beta}/R=\mathscr{O}({\it\epsilon})$
and
$h/R=\mathscr{O}({\it\epsilon})$
), the leading term in flux
$q$
takes the form
$$\begin{eqnarray}q\approx -\left(\frac{h^{3}}{3}+{\it\beta}h^{2}\right)\frac{1}{{\it\mu}}\frac{\partial P}{\partial {\it\theta}}\quad \text{for }h\approx {\it\beta}\ll R,\end{eqnarray}$$
which is the usual result for incorporating the Navier slip model in the thin-film equations. However, taking the formal limit
${\it\beta}\rightarrow \infty$
, which corresponds to perfect slip on the sphere surface, would result in
$$\begin{eqnarray}q=\frac{h}{6}(3R^{2}+3Rh+h^{2})\,\frac{1}{{\it\mu}}\frac{\partial P}{\partial {\it\theta}}\approx \frac{R^{2}h}{2}\frac{1}{{\it\mu}}\frac{\partial P}{\partial {\it\theta}}\quad \text{for }h\ll R.\end{eqnarray}$$
We see that the ‘mobility’ factor that relates the flux
$q$
to the pressure gradient scales as
$h^{3}$
for no slip and as
$h$
for perfect slip, but it changes sign in the two limits; with perfect slip, the flux appears to move ‘uphill’ against the pressure! The limit of perfect slip on the surface does not appear to be physically well posed.
Obviously, one needs to be careful with the scaling of the slip length
${\it\beta}$
. Physically
${\it\beta}$
is thought to be extremely small, so in all likelihood
${\it\beta}\ll h\ll R$
except where
$h\rightarrow 0$
at the touchdown points of the dry zones. On the other hand, to assume perfect slip we must take
${\it\beta}\rightarrow \infty$
, which would imply that
$h\ll R\ll {\it\beta}$
(the case that does not appear to be well posed). One can also consider intermediate regimes when
$h\ll {\it\beta}\ll R$
, or as we saw above, cases when
${\it\beta}$
and
$h$
are of comparable magnitudes and both small compared to
$R$
.
