A.1 The paper by Benilov et al. (Reference Benilov, Benilov and Kopteva2008)

In the work by Benilov et al. (Reference Benilov, Benilov and Kopteva2008), the steady flow of a thin liquid film on the inner surface of a rotating cylinder was investigated using the following flux equation:

(A 1) $$\begin{eqnarray}q=h-\frac{h^{3}}{3}\cos \unicode[STIX]{x1D703}+\unicode[STIX]{x1D716}\frac{h^{3}}{3}\left(\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}^{3}h+\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}h\right),\end{eqnarray}$$

where $h$ is the film thickness made dimensionless with $\unicode[STIX]{x1D6FC}R$ with $\unicode[STIX]{x1D6FC}$ given by

(A 2) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}=\sqrt{\frac{\unicode[STIX]{x1D707}\unicode[STIX]{x1D714}}{\unicode[STIX]{x1D70C}gR}}.\end{eqnarray}$$

and

(A 3) $$\begin{eqnarray}\unicode[STIX]{x1D716}=\frac{\unicode[STIX]{x1D6FE}}{\unicode[STIX]{x1D70C}gR^{2}}\sqrt{\frac{\unicode[STIX]{x1D707}\unicode[STIX]{x1D714}}{\unicode[STIX]{x1D70C}gR}}\end{eqnarray}$$

represents a dimensionless surface tension parameter. The flux equation (A 1) must be solved together with the integral condition

(A 4) $$\begin{eqnarray}M=\int _{0}^{2\unicode[STIX]{x03C0}}h\,\text{d}\unicode[STIX]{x1D703},\end{eqnarray}$$

where $M$ is a (dimensionless) parameter characterising the net mass of fluid inside the cylinder, and the periodicity condition

(A 5) $$\begin{eqnarray}h(\unicode[STIX]{x1D703})=h(\unicode[STIX]{x1D703}+2\unicode[STIX]{x03C0}).\end{eqnarray}$$

Equations (A 1) and (A 4) are, essentially, the SLM equations (2.16) and (2.17) written in a different parametrisation and we note that they also can be derived using standard lubrication theory by assuming that $\unicode[STIX]{x1D6FC}\ll 1$ . Using a combination of numerics and asymptotics, Benilov et al. (Reference Benilov, Benilov and Kopteva2008) found that for some small surface tension values, there are certain values of mass $M$ that present up to three steady-state solutions. Figure 16 shows an example of such a situation, as obtained from the numerical solution of (A 1) for $\unicode[STIX]{x1D716}=10^{-3}$ . Figure 17 shows the three solutions found at $M=10.5$ , each with a different value of $q$ . Although the solutions look rather similar, there are clear distinctions in their peak and tail regions.

Figure 16.

Flux $q$ for rimming flow as a function of the net mass $M$ for $\unicode[STIX]{x1D716}=10^{-3}$ and $\unicode[STIX]{x1D6FC}=0.1$ . The solid line shows the numerical solution of the Stokes equations (A 6) and (A 7); the dashed line shows the numerical solution of the standard lubrication model (SLM), equation (A 1); the dotted line shows the numerical solution of the variational lubrication model (VLM), equation (A 13). The lower panel is an expanded view of the region marked by the box in the upper panel. The solutions at $M=10.5$ (marked by circles on the vertical dotted line) are shown in figure 17.

Figure 17.

Film profiles $h(\unicode[STIX]{x1D703})$ obtained from numerical solution of the Stokes equations (solid line), the standard lubrication model (dashed line) and the variational lubrication model (dashed line) for the rimming-flow problem studied by Benilov et al. (Reference Benilov, Benilov and Kopteva2008) for $\unicode[STIX]{x1D6FC}=0.1$ , $\unicode[STIX]{x1D716}=0.001$ and $M=10.5$ .

We shall revisit this result by solving the two-dimensional steady Stokes equations

(A 6) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}p=\unicode[STIX]{x1D6FC}^{3}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}-\unicode[STIX]{x1D6FC}\boldsymbol{j}, & \displaystyle\end{eqnarray}$$

(A 7) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

where $R$ , $\unicode[STIX]{x1D714}R$ and $\unicode[STIX]{x1D70C}gR/\unicode[STIX]{x1D6FC}$ were taken as length, velocity and pressure scales, respectively. We require the solutions of the above equations to satisfy the following boundary conditions

(A 8) $$\begin{eqnarray}\boldsymbol{u}=\boldsymbol{t}\quad \text{(no-slip condition)}\end{eqnarray}$$

at the surface of the cylinder,

(A 9) $$\begin{eqnarray}\displaystyle & \displaystyle \left\{-p\boldsymbol{I}+\frac{1}{\unicode[STIX]{x1D6FC}^{3}}\left[\unicode[STIX]{x1D735}\boldsymbol{u}+(\unicode[STIX]{x1D735}\boldsymbol{u})^{\text{T}}\right]\right\}\boldsymbol{\cdot }\boldsymbol{n}=\unicode[STIX]{x1D716}\unicode[STIX]{x1D705}\boldsymbol{n}\quad \text{(dynamic condition)} & \displaystyle\end{eqnarray}$$

(A 10) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}=0\quad \text{(kinematic condition)} & \displaystyle\end{eqnarray}$$

at the free surface and we also require that

(A 11) $$\begin{eqnarray}\oint \boldsymbol{x}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}s=\unicode[STIX]{x1D6FC}M\quad \text{(integral condition)},\end{eqnarray}$$

where the line integral is taken along the free-surface boundary. In addition, we also solve the corresponding variational lubrication model (VLM)

(A 12) $$\begin{eqnarray}\displaystyle & \displaystyle q=h-\frac{h^{3}}{3}\left[-\unicode[STIX]{x1D716}\unicode[STIX]{x1D705}+(1-\unicode[STIX]{x1D6FC}h)\sin \unicode[STIX]{x1D703}\right], & \displaystyle\end{eqnarray}$$

(A 13) $$\begin{eqnarray}\displaystyle & \displaystyle M=\int _{0}^{2\unicode[STIX]{x03C0}}\left(h-\unicode[STIX]{x1D6FC}\frac{h^{2}}{2}\right)\,\text{d}\unicode[STIX]{x1D703}. & \displaystyle\end{eqnarray}$$

A comparison of the fluxes obtained from the numerical solution of the Stokes equations as well as from the SLM and the VLM as a function of the net mass $M$ is shown in figure 16 for $\unicode[STIX]{x1D6FC}=0.1$ and $\unicode[STIX]{x1D716}=10^{-3}$ . The multiple solutions loop found by Benilov et al. (Reference Benilov, Benilov and Kopteva2008) is an artefact of the SLM, because it is neither found in the full model nor in the VLM at these parameter values. Moreover, we note from figure 17 that the film profiles $h(\unicode[STIX]{x1D703})$ obtained with the VLM and with the Stokes equations agree really well (shown for $M=10.5$ ). The solution loop is also absent in the ELM. If one continuously transforms SLM into ELM by means of a homotopy parameter that multiplies the additional terms, one finds that the loop of solutions vanishes in a cusp bifurcation where the limit points merge. We attribute this significant qualitative difference between the solutions obtained from the SLM and those from the ELM, VLM and Stokes equations to the fact that for these solutions a large ‘pool’ of fluid forms. In the ‘pool’ zone the film height peaks and large slopes in the free-surface height can be observed, which means that locally the hydrostatic pressure terms are no longer $O(\unicode[STIX]{x1D6FC})$ (a small perturbation), but actually $O(1)$ and therefore should be promoted to leading order. As the film thickness reduces the magnitude of the hydrostatic terms will decrease unless the slopes in the ‘pool’ region increase further to compensate. Hence, the multiple solutions found in the SLM would be expected to occur in the Stokes equations for sufficiently thin films. We have found no evidence that multiple solutions occur at these parameter values in the Stokes equations for film thicknesses down to $\unicode[STIX]{x1D6FC}=0.05$ , although the distinction between the solutions shown in figure 17 becomes very difficult to detect for such small absolute values of film thickness.

A.2 The paper by Benilov et al. (Reference Benilov, Lapin and O’Brien2012)

This paper examined two distinct regimes of rimming flows on a rotating cylinder: one where surface tension effects are taken into account and one where they are not. In the present work, we shall focus on the former where the authors showed, by solving the Stokes equations numerically, that if a dimensionless capillary coefficient, $\unicode[STIX]{x1D6FE}/(\unicode[STIX]{x1D707}\unicode[STIX]{x1D714}R)$ , exceeds a threshold value of 0.23, then a smoothed-shock steady solution always exists. However, the authors claim that if the surface tension is not sufficiently strong (i.e. the capillary coefficient does not exceed the given threshold value) there is a certain ‘forbidden’ region in the parameter space for which the problem does not have stable steady solutions.

In order to verify whether or not the solutions found in the present work are in agreement with the results found in Benilov et al. (Reference Benilov, Lapin and O’Brien2012) it is necessary first to rescale our variables (they chose $R$ , $\unicode[STIX]{x1D714}R$ , $\unicode[STIX]{x1D707}\unicode[STIX]{x1D714}$ and $\unicode[STIX]{x1D714}^{-1}$ as length, velocity, pressure and time scales) so that we can interpret our results in terms of the same dimensionless parameters. This can be achieved by introducing the following mappings which transform our dimensionless variables to those used by Benilov et al. (Reference Benilov, Lapin and O’Brien2012):

(A 14) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\boldsymbol{x}\rightarrow \boldsymbol{x},\quad \boldsymbol{u}\rightarrow \unicode[STIX]{x1D6FA}\boldsymbol{u},\quad p\rightarrow \unicode[STIX]{x1D6FA}p,\quad t\rightarrow \unicode[STIX]{x1D6FA}^{-1}t,\\ \unicode[STIX]{x1D748}\rightarrow \unicode[STIX]{x1D6FA}\unicode[STIX]{x1D748},\quad \unicode[STIX]{x1D705}\rightarrow \unicode[STIX]{x1D705}.\end{array}\right\}\end{eqnarray}$$

In terms of (A 14), the Stokes equations (2.1) take the form

(A 15) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}p=\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}-\frac{1}{\unicode[STIX]{x1D6FC}^{2}}\,\boldsymbol{j}, & \displaystyle\end{eqnarray}$$

(A 16) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

where the dimensionless parameter $\unicode[STIX]{x1D6FC}$ is given by

(A 17) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}=\sqrt{\unicode[STIX]{x1D6FA}}=\sqrt{\frac{\unicode[STIX]{x1D707}\unicode[STIX]{x1D714}}{\unicode[STIX]{x1D70C}gR}}.\end{eqnarray}$$

The no-slip condition (2.3) now corresponds to

(A 18) $$\begin{eqnarray}\boldsymbol{u}=\boldsymbol{t}\end{eqnarray}$$

at the surface of the rotating cylinder. Moreover, the kinematic condition (2.5) retains its original form

(A 19) $$\begin{eqnarray}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}=\unicode[STIX]{x2202}_{t}\boldsymbol{x}\boldsymbol{\cdot }\boldsymbol{n}\end{eqnarray}$$

at the free surface. On the other hand, the dynamic condition (2.6) becomes

(A 20) $$\begin{eqnarray}\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{n}=Ca^{-1}\unicode[STIX]{x1D705}\boldsymbol{n},\end{eqnarray}$$

where

(A 21) $$\begin{eqnarray}Ca^{-1}=\unicode[STIX]{x1D6FA}^{-1}Bo^{-1}=\frac{\unicode[STIX]{x1D6FE}}{\unicode[STIX]{x1D707}\unicode[STIX]{x1D714}R}\end{eqnarray}$$

is the inverse capillary number (note that $Ca^{-1}$ is the aforementioned ‘capillary coefficient’). To introduce the last parameter, it is convenient to note that if we scale the dimensionless film thickness $h$ with $\unicode[STIX]{x1D6FC}$ , then the mass conservation condition (2.20) requires that

(A 22) $$\begin{eqnarray}\int _{0}^{2\unicode[STIX]{x03C0}}\left(h-\unicode[STIX]{x1D6FC}\frac{h^{2}}{2}\right)\,\text{d}\unicode[STIX]{x1D703}=\frac{2\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D6FC}}\left(\unicode[STIX]{x1D706}-\frac{\unicode[STIX]{x1D706}^{2}}{2}\right)\equiv V,\end{eqnarray}$$

where $V$ is the volume of liquid within the cylinder. Moreover, it is possible to show (see, e.g. Benjamin et al. (Reference Benjamin, Pritchard and Tavener1993) and O’Brien & Gath (Reference O’Brien and Gath1998)), in the lubrication limit that a shock solution can only exist if

(A 23) $$\begin{eqnarray}V_{1}\leqslant V<V_{2},\end{eqnarray}$$

where $V_{1}\approx 4.44$ and $V_{2}\approx 6.93$ . We note that if $V<V_{1}$ Moffatt (Reference Moffatt1977) demonstrated that the film coats the inside of the rotating cylinder in a nearly uniform fashion. On the other hand, if $V>V_{2}$ then a pool starts to form at the bottom of the cylinder. For $\unicode[STIX]{x1D706}=0.1$ , equations (A 22)–(A 23) give

(A 24) $$\begin{eqnarray}0.086\lesssim \unicode[STIX]{x1D6FC}\lesssim 0.134\end{eqnarray}$$

and therefore, on noting that $\unicode[STIX]{x1D6FA}=\unicode[STIX]{x1D706}^{2}\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D706}}$ , and using (A 17), the relevant range for $\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D706}}$ is

(A 25) $$\begin{eqnarray}0.740\lesssim \unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D706}}\lesssim 1.796.\end{eqnarray}$$

In the present work, we investigated the following Bond numbers:

1. (i) $Bo_{\unicode[STIX]{x1D706}}=20$ ,

2. (ii) $Bo_{\unicode[STIX]{x1D706}}=120$ ,

3. (iii) $Bo_{\unicode[STIX]{x1D706}}=500$ ,

4. (iv) $Bo_{\unicode[STIX]{x1D706}}=2000$ ,

which corresponds, in view of (A 25) and given that $Bo_{\unicode[STIX]{x1D706}}=\unicode[STIX]{x1D706}^{-1}Bo$ , to

1. (i) $27.84\lesssim Ca^{-1}\lesssim 67.57$ ,

2. (ii) $4.64\lesssim Ca^{-1}\lesssim 11.26$ ,

3. (iii) $1.11\lesssim Ca^{-1}\lesssim 2.70$ ,

4. (iv) $0.28\lesssim Ca^{-1}\lesssim 0.68$ .

Note that in all cases, $Ca^{-1}$ is larger than 0.23 (the threshold value) and therefore the existence of steady solutions in these intervals is in accordance with the findings in Benilov et al. (Reference Benilov, Lapin and O’Brien2012).

Figure 18.

Example of a rimming flow with a shock-like region as obtained from the numerical solution of the Stokes equations for $\unicode[STIX]{x1D706}=0.1$ , $\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D706}}=1$ and $Bo_{\unicode[STIX]{x1D706}}=7800$ . The panel (b) is an expanded view of the shock-like region marked by the box in the panel (a).

Furthermore, we decided to investigate the so-called ‘forbidden’ region by solving the Stokes equations for $\unicode[STIX]{x1D706}=0.1$ , $\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D706}}=1$ and $Bo_{\unicode[STIX]{x1D706}}=7800$ , which yields $Ca^{-1}\approx 0.13$ and $\unicode[STIX]{x1D6FC}=0.1$ . As can be seen in figure 18, although $Ca^{-1}$ is below the threshold value, a stable shock-like (corner) steady-state solution is found for this set of parameters. We note that the numerical scheme used by Benilov et al. (Reference Benilov, Lapin and O’Brien2012) assumes that the interface position is given by a function $h(\unicode[STIX]{x1D703})$ , which cannot represent shock-like regions where $h(\unicode[STIX]{x1D703})$ becomes multi-valued. Moreover, Benilov et al. (Reference Benilov, Lapin and O’Brien2012) used a time-dependent solver to approach steady solutions which does not allow unstable solutions to be found. In addition, the restriction on interface shapes prevents computation of stable solutions that are reached via transient regions of interface overturning. We conclude that the results in Benilov et al. (Reference Benilov, Lapin and O’Brien2012) when surface tension is taken into account should be reviewed using a formulation that allows the computation of overturning free surfaces.

B.1 Onsager’s variational principle

Let $a_{1},a_{2},\ldots ,a_{n}$ denote the fluctuations from equilibrium of a set of macroscopic thermodynamic quantities of a closed system, in the following called ‘state variables’. Also, let $F=F(a_{1},a_{2},\ldots ,a_{n})$ be the free energy, and let $F_{e}=F(0,0,\ldots ,0)$ be its minimum at equilibrium. Then, $\unicode[STIX]{x0394}F\equiv F-F_{e}\geqslant 0$ can be expressed as a quadratic form

(B 1) $$\begin{eqnarray}\unicode[STIX]{x0394}F=\frac{1}{2}\mathop{\sum }_{i=1}^{n}\mathop{\sum }_{j=1}^{n}\unicode[STIX]{x1D6FD}_{ij}a_{i}a_{j},\end{eqnarray}$$

where $\unicode[STIX]{x1D737}$ is a real positive-definite symmetric matrix.

The thermodynamic conjugate to $a_{i}$ are the thermodynamic forces defined as

(B 2) $$\begin{eqnarray}X_{i}=-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x0394}F}{\unicode[STIX]{x2202}a_{i}}.\end{eqnarray}$$

For small fluctuations, they can be expressed as linear combinations of the state variables, that is

(B 3) $$\begin{eqnarray}X_{i}=-\mathop{\sum }_{j=1}^{n}\unicode[STIX]{x1D6FD}_{ij}a_{j}.\end{eqnarray}$$

Assuming that the system is not far from equilibrium, the state variables $a_{1},a_{2},\ldots ,a_{n}$ evolve according to the linear relation to the forces

(B 4) $$\begin{eqnarray}{\dot{a}}_{i}\equiv \frac{\text{d}}{\text{d}t}a_{i}(t)=\mathop{\sum }_{j=1}^{n}\unicode[STIX]{x1D613}_{ij}X_{j}(t),\end{eqnarray}$$

where the $\unicode[STIX]{x1D613}_{ij}$ are called kinetic coefficients. The Onsager principle (Landau & Lifshitz Reference Landau and Lifshitz1980), also known as principle of symmetry of kinetic coefficients, states that $\unicode[STIX]{x1D613}_{ij}=L_{ji}$ , that is, $\unicode[STIX]{x1D613}$ is a symmetric matrix. Based on this principle, Onsager formulated a variational principle applicable to some irreversible transport processes such as heat conduction, electrical conduction and diffusion.

For isothermal systems, Onsager’s variational principle states that the state evolution equation can be obtained by minimising the so-called Rayleighian (Doi Reference Doi2013)

(B 5) $$\begin{eqnarray}{\mathcal{R}}=\mathop{\sum }_{i=1}^{n}\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}a_{i}}{\dot{a}}_{i}+\frac{1}{2}\mathop{\sum }_{i=1}^{n}\mathop{\sum }_{j=1}^{n}\unicode[STIX]{x1D701}_{ij}{\dot{a}}_{i}{\dot{a}}_{j}\end{eqnarray}$$

with respect to $\dot{a_{1}},\dot{a_{2}},\ldots ,\dot{a_{n}}$ , that is the rates of change of the macroscopic variables $a_{1},a_{2},\ldots ,a_{n}$ . The two terms on the right-hand side are physically distinct. The first term is ${\dot{F}}\equiv \text{d}F/\text{d}t$ , which is the rate of change of the free energy of the system $F=F(a_{1},\ldots ,a_{n})$ . The second term is the energy dissipation function $\unicode[STIX]{x1D6F7}$ , which is a quadratic form with the friction coefficient matrix $\unicode[STIX]{x1D701}$ forming a positive definite and symmetric matrix.

Minimisation of ${\mathcal{R}}$ with respect to the rates ${\dot{a}}_{i}$ leads to the kinetic equations

(B 6) $$\begin{eqnarray}-\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}a_{i}}=\mathop{\sum }_{j=1}^{n}\unicode[STIX]{x1D701}_{ij}{\dot{a}}_{j},\end{eqnarray}$$

which expresses a balance between two forces: the reversible force $-\unicode[STIX]{x2202}F/\unicode[STIX]{x2202}a_{i}$ and the friction force $\sum _{j=1}^{n}\unicode[STIX]{x1D701}_{ij}{\dot{a}}_{j}$ . Inverting (B 6), one obtains (B 4), i.e. $\unicode[STIX]{x1D701}_{ij}^{-1}=\unicode[STIX]{x1D613}_{ij}$ . Also note, that substitution of (B 6) into (B 5) shows that

(B 7) $$\begin{eqnarray}{\dot{F}}=-2\unicode[STIX]{x1D6F7}.\end{eqnarray}$$

The Onsager variational principle can be then summarised as follows.

1. (i) Choose appropriate state variables $a_{1},a_{2},\ldots ,a_{n}$ .

2. (ii) Construct the free energy function $F=F(a_{1},\ldots ,a_{n})$ .

3. (iii) Construct the energy dissipation function $\unicode[STIX]{x1D6F7}$ .

4. (iv) Construct the Rayleighian ${\mathcal{R}}={\dot{F}}+\unicode[STIX]{x1D6F7}$ .

5. (v) Obtain the time evolution of the state (kinetic equations) by minimising ${\mathcal{R}}$ with respect to the rates $\dot{a_{i}}$ .

B.2 Application of Onsager’s variational principle to our problem

We note that all the following physical variables and quantities are, unless otherwise stated, dimensional.

The free energy (per unit length) function $F=F[h]$ is defined as

(B 8) $$\begin{eqnarray}F=\underbrace{\unicode[STIX]{x1D6FE}\oint \text{d}s}_{\mathit{surface}\,\mathit{energy}\,(\mathit{per}\,\mathit{unit}\,\mathit{length})}+\underbrace{\unicode[STIX]{x1D70C}g\iint r\sin \unicode[STIX]{x1D703}\,\text{d}A}_{\mathit{potential}\,\mathit{energy}\,(\mathit{per}\,\mathit{unit}\,\mathit{length})},\end{eqnarray}$$

where the line integral is taken over the fluid perimeter and the area integral is taken over the fluid domain. Alternatively, we can write the previous equation as

(B 9) $$\begin{eqnarray}F=\int _{0}^{2\unicode[STIX]{x03C0}}(\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D709}+f)\,\text{d}\unicode[STIX]{x1D703},\end{eqnarray}$$

where $\unicode[STIX]{x1D709}\,\text{d}\unicode[STIX]{x1D703}\equiv \sqrt{(R\pm h)^{2}+(\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}h)^{2}}\,\text{d}\unicode[STIX]{x1D703}$ corresponds to the line element $\text{d}s$ in polar coordinates and $f=f(h)$ is given by

(B 10) $$\begin{eqnarray}f=\pm \unicode[STIX]{x1D70C}g\int _{R}^{R\pm h}r^{2}\sin \unicode[STIX]{x1D703}\,\text{d}r=\pm \frac{\unicode[STIX]{x1D70C}g}{3}[(R\pm h)^{3}-R^{3}]\sin \unicode[STIX]{x1D703}.\end{eqnarray}$$

The rate of change of $F$ is given by

(B 11) $$\begin{eqnarray}{\dot{F}}=\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D703}\unicode[STIX]{x2202}_{t}h\frac{\unicode[STIX]{x1D6FF}F}{\unicode[STIX]{x1D6FF}h}=\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D703}\unicode[STIX]{x2202}_{t}h\frac{\unicode[STIX]{x1D6FF}}{\unicode[STIX]{x1D6FF}h}\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D703}^{\prime }(\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D709}+f),\end{eqnarray}$$

and by using the linearity property of the functional derivative, this can be simplified to

(B 12) $$\begin{eqnarray}{\dot{F}}=\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D703}\unicode[STIX]{x2202}_{t}h\frac{\unicode[STIX]{x1D6FF}}{\unicode[STIX]{x1D6FF}h}\left(\unicode[STIX]{x1D6FE}\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D703}^{\prime }\unicode[STIX]{x1D709}+\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D703}^{\prime }f\right).\end{eqnarray}$$

We note that

(B 13) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D6FF}}{\unicode[STIX]{x1D6FF}h}\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D703}^{\prime }\unicode[STIX]{x1D709} & = & \displaystyle \frac{\unicode[STIX]{x1D6FF}}{\unicode[STIX]{x1D6FF}h}\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D703}^{\prime }\left[\sqrt{(R\pm h)^{2}+\left(\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}h\right)^{2}}\right]\nonumber\\ \displaystyle & = & \displaystyle \pm \frac{(R\pm h)}{\sqrt{(R\pm h)^{2}+(\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}h)^{2}}}-\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}\left[\frac{\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}h}{\sqrt{(R\pm h)^{2}+(\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}h)^{2}}}\right]\nonumber\\ \displaystyle & = & \displaystyle \frac{\pm (R\pm h)-\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}^{2}h}{\sqrt{(R\pm h)^{2}+(\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}h)^{2}}}+\frac{(\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}h)^{2}\left[\pm (R\pm h)+\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}^{2}h\right]}{\left[(R\pm h)^{2}+\left(\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}h\right)^{2}\right]^{3/2}}\nonumber\\ \displaystyle & = & \displaystyle (R\pm h)\left\{\frac{\pm (R\pm h)^{2}\pm 2(\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}h)^{2}-(R\pm h)\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}^{2}h}{\left[(R\pm h)^{2}+(\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}h)^{2}\right]^{3/2}}\right\}\nonumber\\ \displaystyle & = & \displaystyle -(R\pm h)\unicode[STIX]{x1D705},\end{eqnarray}$$

where the last equality follows from (2.23), and that

(B 14) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D6FF}}{\unicode[STIX]{x1D6FF}h}\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D703}^{\prime }f={\displaystyle \frac{\text{d}f}{\text{d}h}}=\unicode[STIX]{x1D70C}g(R\pm h)^{2}\sin (\unicode[STIX]{x1D703}).\end{eqnarray}$$

The variation of the free energy with respect to the rate $\unicode[STIX]{x2202}_{t}h$ is then given by

(B 15) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}{\dot{F}}=\int _{0}^{2\unicode[STIX]{x03C0}}[-\unicode[STIX]{x1D6FE}(R\pm h)\unicode[STIX]{x1D705}+f^{\prime }]\unicode[STIX]{x1D6FF}(\unicode[STIX]{x2202}_{t}h)\,\text{d}\unicode[STIX]{x1D703},\end{eqnarray}$$

where $f^{\prime }=f^{\prime }(h)=\text{d}f/\text{d}h$ . The (viscous) dissipation function $\unicode[STIX]{x1D6F7}$ , which is the rate at which (mechanical) work is converted into internal thermal energy, is given by

(B 16) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}=\frac{1}{2}\iint \left(\boldsymbol{T}\boldsymbol{ : }\unicode[STIX]{x1D735}\boldsymbol{u}\right)\,\text{d}A.\end{eqnarray}$$

Note that, in polar coordinates, the double dot product $\boldsymbol{T}\boldsymbol{ : }\unicode[STIX]{x1D735}\boldsymbol{u}$ for an incompressible Newtonian liquid can be written as (Bird, Stewart & Lightfoot Reference Bird, Stewart and Lightfoot2002):

(B 17) $$\begin{eqnarray}\boldsymbol{T}\boldsymbol{ : }\unicode[STIX]{x1D735}\boldsymbol{u}=2\unicode[STIX]{x1D707}\left[(\unicode[STIX]{x2202}_{r}u)^{2}+\left(\frac{\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}v}{r}+\frac{v}{r}\right)^{2}\right]+\unicode[STIX]{x1D707}\left[\unicode[STIX]{x2202}_{r}v-\frac{v}{r}+\frac{\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}u}{r}\right]^{2}.\end{eqnarray}$$

In the thin-film limit, $\boldsymbol{T}\boldsymbol{ : }\unicode[STIX]{x1D735}\boldsymbol{u}$ is represented simply by $\unicode[STIX]{x1D707}(\unicode[STIX]{x2202}_{r}v)^{2}$ and, consequently,

(B 18) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}=\frac{\unicode[STIX]{x1D707}}{2}\iint (\unicode[STIX]{x2202}_{r}v)^{2}\,\text{d}A.\end{eqnarray}$$

The variation of $\unicode[STIX]{x1D6F7}$ with respect to $v$ is then given from the definition of variational derivative,

(B 19) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6F7} & = & \displaystyle {\displaystyle \frac{\unicode[STIX]{x1D707}}{2}}\iint \lim _{\unicode[STIX]{x1D716}\rightarrow 0}\left\{{\displaystyle \frac{\left[\unicode[STIX]{x2202}_{r}(v+\unicode[STIX]{x1D716}\unicode[STIX]{x1D6FF}v)\right]^{2}-(\unicode[STIX]{x2202}_{r}v)^{2}}{\unicode[STIX]{x1D716}}}\right\}\,\text{d}A\nonumber\\ \displaystyle & = & \displaystyle {\displaystyle \frac{\unicode[STIX]{x1D707}}{2}}\iint \lim _{\unicode[STIX]{x1D716}\rightarrow 0}\left\{2\unicode[STIX]{x2202}_{r}v\unicode[STIX]{x2202}_{r}(\unicode[STIX]{x1D6FF}v)+\unicode[STIX]{x1D716}\left[\unicode[STIX]{x2202}_{r}(\unicode[STIX]{x1D6FF}v)\right]^{2}\right\}\,\text{d}A\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D707}\iint \unicode[STIX]{x2202}_{r}v\unicode[STIX]{x2202}_{r}\left(\unicode[STIX]{x1D6FF}v\right)\,\text{d}A\nonumber\\ \displaystyle & = & \displaystyle \pm \unicode[STIX]{x1D707}\int _{0}^{2\unicode[STIX]{x03C0}}\int _{R}^{R\pm h}(r\unicode[STIX]{x2202}_{r}v)\unicode[STIX]{x2202}_{r}(\unicode[STIX]{x1D6FF}v)\,\text{d}r\,\text{d}\unicode[STIX]{x1D703}.\end{eqnarray}$$

A single integration by parts in the previous integral yields

(B 20) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6F7} & = & \displaystyle \left.\mp \unicode[STIX]{x1D707}\int _{0}^{2\unicode[STIX]{x03C0}}\int _{R}^{R\pm h}\unicode[STIX]{x2202}_{r}(r\unicode[STIX]{x2202}_{r}v)\unicode[STIX]{x1D6FF}v\,\text{d}r\,\text{d}\unicode[STIX]{x1D703}\pm \unicode[STIX]{x1D707}\int _{0}^{2\unicode[STIX]{x03C0}}r\unicode[STIX]{x2202}_{r}v\unicode[STIX]{x1D6FF}v\right|_{r=R}^{r=R\pm h}\,\text{d}\unicode[STIX]{x1D703}\nonumber\\ \displaystyle & = & \displaystyle \left.\mp \unicode[STIX]{x1D707}\int _{0}^{2\unicode[STIX]{x03C0}}\int _{R}^{R\pm h}r\left(\unicode[STIX]{x2202}_{r}^{2}v+\frac{\unicode[STIX]{x2202}_{r}v}{r}\right)\unicode[STIX]{x1D6FF}v\,\text{d}r\,\text{d}\unicode[STIX]{x1D703}\pm \unicode[STIX]{x1D707}\int _{0}^{2\unicode[STIX]{x03C0}}r\unicode[STIX]{x2202}_{r}v\unicode[STIX]{x1D6FF}v\right|_{r=R}^{r=R\pm h}\,\text{d}\unicode[STIX]{x1D703}.\end{eqnarray}$$

In the thin-film limit, the previous equation reduces to

(B 21) $$\begin{eqnarray}\left.\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6F7}=\mp \unicode[STIX]{x1D707}R\int _{0}^{2\unicode[STIX]{x03C0}}\int _{R}^{R\pm h}\unicode[STIX]{x2202}_{r}^{2}v\unicode[STIX]{x1D6FF}v\,\text{d}r\,\text{d}\unicode[STIX]{x1D703}\pm \unicode[STIX]{x1D707}\int _{0}^{2\unicode[STIX]{x03C0}}r\unicode[STIX]{x2202}_{r}v\unicode[STIX]{x1D6FF}v\right|_{r=R}^{r=R\pm h}\,\text{d}\unicode[STIX]{x1D703}.\end{eqnarray}$$

Next, we use the pressure $p(r,\unicode[STIX]{x1D703})$ as a Lagrange multiplier to impose the incompressibility constraint, $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0$ , by adding the following term $C$ to the Rayleighian ${\mathcal{R}}$ :

(B 22) $$\begin{eqnarray}C=\iint [-p(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u})]\,\text{d}A,\end{eqnarray}$$

which can be integrated by parts using the well-known Gauss’s divergence theorem to give

(B 23) $$\begin{eqnarray}C=\iint \unicode[STIX]{x1D735}p\boldsymbol{\cdot }\boldsymbol{u}\,\text{d}A-\oint pv_{n}\,\text{d}s\end{eqnarray}$$

(B 24) $$\begin{eqnarray}=\pm \int _{0}^{2\unicode[STIX]{x03C0}}\int _{R}^{R\pm h}\left(u\unicode[STIX]{x2202}_{r}p+\frac{v}{r}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}p\right)r\,\text{d}r\,\text{d}\unicode[STIX]{x1D703}\mp \int _{0}^{2\unicode[STIX]{x03C0}}pv_{n}\unicode[STIX]{x1D709}\,\text{d}\unicode[STIX]{x1D703},\end{eqnarray}$$

where $v_{n}=\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}$ is the outward normal velocity at the free surface. The variation of $C$ with respect to the rates $u$ and $v$ is then found to be

(B 25) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}C=\pm \int _{0}^{2\unicode[STIX]{x03C0}}\int _{R}^{R\pm h}\left(\unicode[STIX]{x1D6FF}u\unicode[STIX]{x2202}_{r}p+\frac{\unicode[STIX]{x1D6FF}v}{r}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}p\right)r\,\text{d}r\,\text{d}\unicode[STIX]{x1D703}\mp \int _{0}^{2\unicode[STIX]{x03C0}}p\unicode[STIX]{x1D6FF}v_{n}\unicode[STIX]{x1D709}\,\text{d}\unicode[STIX]{x1D703}.\end{eqnarray}$$

Using the kinematic condition

(B 26) $$\begin{eqnarray}v_{n}=(R\pm h)\unicode[STIX]{x1D709}^{-1}\unicode[STIX]{x2202}_{t}h\quad \text{at}~r=R\pm h,\end{eqnarray}$$

we can write $\unicode[STIX]{x1D6FF}C$ as

(B 27) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}C=\pm \int _{0}^{2\unicode[STIX]{x03C0}}\int _{R}^{R\pm h}\left(\unicode[STIX]{x1D6FF}u\unicode[STIX]{x2202}_{r}p+\frac{\unicode[STIX]{x1D6FF}v}{r}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}p\right)r\,\text{d}r\,\text{d}\unicode[STIX]{x1D703}\mp \int _{0}^{2\unicode[STIX]{x03C0}}p(R\pm h)\unicode[STIX]{x1D6FF}\left(\unicode[STIX]{x2202}_{t}h\right)\,\text{d}\unicode[STIX]{x1D703}.\end{eqnarray}$$

Hence, from the expressions obtained above we can write $\unicode[STIX]{x1D6FF}{\mathcal{R}}=\unicode[STIX]{x1D6FF}{\dot{F}}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6F7}+\unicode[STIX]{x1D6FF}C$ as a linear combination of $\unicode[STIX]{x1D6FF}u$ , $\unicode[STIX]{x1D6FF}v$ and $\unicode[STIX]{x1D6FF}(\unicode[STIX]{x2202}_{t}h)$ . Minimising ${\mathcal{R}}$ with respect to $\unicode[STIX]{x2202}_{t}h$ gives the pressure at the free surface:

(B 28) $$\begin{eqnarray}p=\mp \unicode[STIX]{x1D6FE}\unicode[STIX]{x1D705}+\unicode[STIX]{x1D70C}g(R\pm h)\sin \unicode[STIX]{x1D703}\quad \text{at}~r=R\pm h.\end{eqnarray}$$

Minimising ${\mathcal{R}}$ with respect to $u$ gives

(B 29) $$\begin{eqnarray}\unicode[STIX]{x2202}_{r}p=0,\end{eqnarray}$$

i.e. the pressure varies uniformly in the radial direction. Minimising ${\mathcal{R}}$ with respect to $v$ leads to

(B 30) $$\begin{eqnarray}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}p-\unicode[STIX]{x1D707}R\unicode[STIX]{x2202}_{r}^{2}v=0.\end{eqnarray}$$

Lastly, minimising ${\mathcal{R}}$ with respect to $v$ evaluated at the free surface gives

(B 31) $$\begin{eqnarray}\unicode[STIX]{x2202}_{r}v=0\quad \text{at}~r=R\pm h,\end{eqnarray}$$

which is the tangential stress boundary condition at the free surface. In the dimensionless scaled form (see § 2), the governing equations are given by

(B 32) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}p=0, & \displaystyle\end{eqnarray}$$

(B 33) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}p-\unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}^{2}v=0, & \displaystyle\end{eqnarray}$$

where the asterisks have been omitted, and they are to be solved together with the boundary conditions

(B 34) $$\begin{eqnarray}\displaystyle & \displaystyle p=-Bo_{\unicode[STIX]{x1D706}}^{-1}\unicode[STIX]{x1D705}+(1+\unicode[STIX]{x1D6FC}h)\sin \unicode[STIX]{x1D703}\quad \text{at}~\unicode[STIX]{x1D702}=h, & \displaystyle\end{eqnarray}$$

(B 35) $$\begin{eqnarray}\displaystyle & \displaystyle v=\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D706}}\quad \text{at}~\unicode[STIX]{x1D702}=0, & \displaystyle\end{eqnarray}$$

(B 36) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{\unicode[STIX]{x1D702}}v=0\quad \text{at}~\unicode[STIX]{x1D702}=h. & \displaystyle\end{eqnarray}$$

Solving (B 32) subject to the boundary condition (B 34), we obtain the following expression for the fluid pressure

(B 37) $$\begin{eqnarray}p=-Bo_{\unicode[STIX]{x1D706}}^{-1}\unicode[STIX]{x1D705}+(1+\unicode[STIX]{x1D6FC}h)\sin \unicode[STIX]{x1D703},\end{eqnarray}$$

which contains the contributions of both Laplace pressure $-Bo^{-1}\unicode[STIX]{x1D705}$ and hydrostatic pressures $(1+\unicode[STIX]{x1D6FC}h)\sin \unicode[STIX]{x1D703}$ . Substituting it in (B 33), and applying the boundary conditions (B 35) and (B 36), we find that

(B 38) $$\begin{eqnarray}v=\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D706}}+\left(\frac{\unicode[STIX]{x1D702}^{2}}{2}-\unicode[STIX]{x1D702}h\right)\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}[-Bo^{-1}\unicode[STIX]{x1D705}+(1+\unicode[STIX]{x1D6FC}h)\sin \unicode[STIX]{x1D703}].\end{eqnarray}$$

Hence, the expression for the dimensionless flux $q=\int _{0}^{h}v\,\text{d}\unicode[STIX]{x1D702}$ is of the form

(B 39) $$\begin{eqnarray}q=\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D706}}h+\frac{h^{3}}{3}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}\left[Bo_{\unicode[STIX]{x1D706}}^{-1}\unicode[STIX]{x1D705}-\left(1+\unicode[STIX]{x1D6FC}h\right)\sin \unicode[STIX]{x1D703}\right].\end{eqnarray}$$

Then, using the conservation law for $h$ given by $(1+\unicode[STIX]{x1D6FC}h)\unicode[STIX]{x2202}_{t}h=-\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}q$ (for derivation, see Benilov & O’Brien Reference Benilov and O’Brien2005), we arrive at the following fourth-order thin-film evolution equation

(B 40) $$\begin{eqnarray}(1+\unicode[STIX]{x1D6FC}h)\unicode[STIX]{x2202}_{t}h+\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}\left\{\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D706}}h+\frac{h^{3}}{3}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}\left[Bo_{\unicode[STIX]{x1D706}}^{-1}\unicode[STIX]{x1D705}-(1+\unicode[STIX]{x1D6FC}h)\sin \unicode[STIX]{x1D703}\right]\right\}=0.\end{eqnarray}$$

Observe that to the leading order the previous equation yields, after some algebra,

(B 41) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}h+\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}\left[\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D706}}h-\frac{h^{3}}{3}\cos \unicode[STIX]{x1D703}+Bo^{-1}\frac{h^{3}}{3}(\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}^{3}h+\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}h)\right]=0,\end{eqnarray}$$

where the $O(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FA},\unicode[STIX]{x1D6FC}Bo^{-1})$ terms have been omitted, which is in agreement with the SLM (2.13) presented in the § 2.

Finally, we highlight that the (B 40) can be presented in the gradient-dynamics form for the conserved field $h$ as

(B 42) $$\begin{eqnarray}(1+\unicode[STIX]{x1D6FC}h)\unicode[STIX]{x2202}_{t}h=-\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}\left\{M-Q\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}\left[(1+\unicode[STIX]{x1D6FC}h)^{-1}\frac{\unicode[STIX]{x1D6FF}F}{\unicode[STIX]{x1D6FF}h}\right]\right\},\end{eqnarray}$$

where the no-slip boundary condition at the surface of the rotating cylinder introduces the term $M=M(h)=\unicode[STIX]{x1D6FA}h$ and $Q=Q(h)=h^{3}/3$ is the mobility factor.