# On the liquid film instability of an internally coated horizontal tube

Shahab Eghbali<EMAIL_ADDRESS>Yves-Marie Ducimetière  Edouard Boujo
François Gallaire  Laboratory of Fluid Mechanics and Instabilities, École
Polytechnique Fédérale de Lausanne, Lausanne CH-1015, Switzerland

###### Abstract

We study numerically and theoretically the gravity-driven flow of a viscous
liquid film coating the inner side of a horizontal cylindrical tube and
surrounding a shear-free dynamically inert gaseous core. The liquid-gas
interface is prone to the Rayleigh-Plateau and Rayleigh-Taylor instabilities.
Here, we focus on the limit of low and intermediate Bond numbers, $Bo$, where
the capillary and gravitational forces are comparable and the Rayleigh-Taylor
instability is known to be suppressed. We first study the evolution of the
axially invariant draining flow, initiating from a uniform film thickness
until reaching a quasi-static regime as the bubble approaches the upper tube
wall. We then investigate the flow’s linear stability within two frameworks:
frozen time-frame (quasi-steady) stability analysis and transient growth
analysis. We explore the effect of the surface tension ($Bo$) and inertia
(measured by the Ohnesorge number, $Oh$) on the flow and its stability. The
linear stability analysis suggests that the interface deformation at large
$Bo$ results in the suppression of the Rayleigh-Plateau instability in the
asymptotic long-time limit. Furthermore, the transient growth analysis
suggests that the initial flow evolution does not lead to any considerable
additional amplification of initial interface perturbations, a posteriori
rationalising the quasi-steady assumption. The present study yields a
satisfactory prediction of the stabilisation threshold found experimentally by
Duclaux et al. [16].

††preprint: APS/123-QED

##  I Introduction

The gravity-driven flow of a viscous liquid film coating inside a solid
cylinder has received attention by virtue of its rich dynamics and numerous
applications. One major industrial application of such a flow is in two-phase
heat exchangers, like in evaporators [32] and vertical tube condensers [35],
where the dynamics of the gas plugs and capillary blockage affect the heat
transfer efficiency and pressure loss [15, 41]. It is also relevant in the
airways of the human lung where the coupling between the moving wall of the
airways and the liquid film coating the wall can cause airway closure [24, 9,
28].

The liquid-gas interface is prone to different instabilities due to surface
tension, shear, gravity, and inertia, which eventually lead to the emergence
of various patterns [17, 20]. For instance, the long column may break apart
into distinct plugs separated axially by collars, thus minimising the surface
energy. This phenomenon is classified in a large family of hydrodynamic
instabilities known as the Rayleigh-Plateau instability [34, 37]. Flow
characteristics of the viscous film coating the inner side of a tube depend
also on its orientation with respect to gravity. In the case of a vertical
tube, where gravity drives an axial flow, the liquid interface is mainly
destabilised by the Rayleigh-Plateau instability [21]. Different stages of
such instability have been investigated through a large number of numerical
and experimental studies. Goren [22] showed analytically that in the absence
of inertia, the length of plugs is set by the interface radius. Frenkel et al.
[19] demonstrated that despite the growth of small disturbances, nonlinear
saturation of instabilities can avoid the rupture of the liquid film in a
certain range of axial flow parameters. Using the long-wavelength
approximation, Camassa et al. [10, 11] showed the existence of non-trivial
axially traveling waves along the interface when the film thickness exceeds a
critical value. Moreover, these studies describe how plugs form via Hopf
bifurcation as the waves grow, and explore experimentally the
absolute/convective properties of the traveling waves. With this general
picture, further studies address the effects of flow parameters on the
dynamics of the flow instability, e.g. air driven flow [12] and multiple
liquid layers [30], wall porosity [29], Marangoni effect [14], and the
presence of surfactant [31].

Alternatively, in the case of liquid film coating the inner side of long
horizontal or inclined tubes, where gravity is not orthogonal to the cross-
section of the tube, a second instability, the Rayleigh-Taylor instability
[38, 40], may arise. As a result, the heavy liquid film accelerates into the
light gaseous core in the direction of the gravitational field, thus forming
suspended droplets, moving lenses, or rivulets [42, 3, 4], extending the case
of a flat overhanging thin liquid film [18]. Therefore, in the case of a
horizontal tube, both Rayleigh-Plateau and Rayleigh-Taylor instabilities may
potentially coexist depending on the film thickness [8, 7]. Trinh et al. [42]
investigated both experimentally and theoretically, in the absence of any
axial flow, the stability of a thin viscous film that coats the underside of a
tube. Focusing on the situations where gravitational forces dominate surface
tension forces, they did not evidence any manifestation of the Rayleigh-
Plateau instability, highlighting that finite wall curvature and sufficiently
high surface tension can suppress the Rayleigh-Taylor instability as the
liquid drains. Using a similar flow configuration, and also employing the
lubrication approximation, Balestra et al. [3] showed that when gravity
overcomes capillarity, the linear transient growth of the disturbances can
lead to the formation of tiny spanwise homogeneous traveling waves at the top
of the interface, which may pinch or decay as they travel downward. Balestra
et al. [4] however showed that spanwise periodic (but streamwise homogeneous)
structures called rivulets display significantly larger transient gains,
rationalising their experimental observations.

By increasing the film thickness and/or lowering the tube curvature, liquid
drainage is enhanced and the Rayleigh-Taylor instability is dampened, thus
allowing the Rayleigh-Plateau instability to occur. Through an experimental
investigation of a wider range of film thicknesses, Duclaux et al. [16]
evidenced that when capillary forces dominate over gravity forces, the
Rayleigh-Plateau instability sets in, but when they stop to dominate, the
Rayleigh-Plateau instability can be suppressed by gravity. In this case,
liquid drainage is faster than capillarity-induced interface undulations and
plug formation. Hence, gravity and surface tension play dual roles in
interface destabilisation. While gravity promotes the Rayleigh-Taylor
instability at the top of the interface, it opposes the formation of hanging
droplets by enhancing drainage progressively while draining down the curved
substrate. In contradistinction, surface tension promotes the Rayleigh-Plateau
instability, whereas it opposes the Rayleigh-Taylor instability.

Duclaux et al. [16] complemented their experimental measurements with a
theoretical analysis based on the lubrication theory and assuming a circular
interface. This analysis correctly predicted the scaling of the stability
threshold in terms of film thickness and Bond number and qualitatively
captured the increasing trend of the instability wavelength with the initial
radius of the liquid-air interface. It did not, however, achieve a fully
satisfactory quantitative agreement. In this study, we set to investigate
accurately how gravity-induced drainage suppresses the Rayleigh-Plateau
instability. Crucially, we consider the full Navier-Stokes equations and take
into account the temporal evolution of the interface. We study the linear
instability in two ways: (i) with a linear stability analysis performed at
each instant taken in isolation, assuming that the base flow is ”frozen”; (ii)
with a transient growth analysis rigorously accounting for the temporal
evolution of the base flow.

This paper is structured as follows. The methodology is presented in §II. The
problem formulation and the dimensionless governing equations are presented in
§II.1, from which the base flow is deduced and discussed in §II.2. In §II.3,
the stability analysis formulation and the linearised governing equations are
elaborated. Then, in §II.4 the formulation of the linear transient growth of
perturbations is detailed. Corresponding numerical methods are detailed in
§II.5. In §III, the results of the stability and transient growth analyses are
presented and discussed. Section III.1 summarises the effect of different
dimensionless parameters on the stability of the flow: §III.1.2-III.1.3
present the effect of surface tension and inertia, respectively, and §III.1.4
presents the linear stability diagram. The results of the transient growth
analysis are presented in §III.2 and finally, conclusions are drawn in §IV by
comparing the linear stability and transient growth analyses.

##  II Governing equations and methods

###  II.1 Problem formulation

Figure 1: Schematic of the liquid film coating the inner side of a horizontal
tube and the geometrical parameters. The thick solid black line shows the tube
wall of inner radius $R_{\text{w}}$, which is concentric with the coordinate
reference. Initially surrounded by a constant film thickness, the liquid-gas
interface is a cylinder of radius $R$, concentric with the tube. The dashed
black line represents the perturbed liquid-gas interface of local radius
$r_{\text{int}}$ and axial wavelength $\lambda R$. The inset shows the zoomed
cross-section of the initial perturbed interface. The gravitational field acts
vertically, perpendicular to the tube axis.

The inner side of a cylindrical tube of radius $R_{\text{w}}$ is coated with a
viscous liquid film. The schematic of the flow is shown in figure 1. The
standard Cartesian coordinates $(x,y,z)$ are considered with the origin
located at the center of the tube cross-section. In-plane coordinates are
$(x,y)$ and the gravity acceleration, ${\mathbf{g}}$, points in $-y$
direction. The liquid is Newtonian, of constant dynamic viscosity $\mu$,
surface tension $\gamma$, and density $\rho$, and surrounds a core bubble of
inviscid gas of density much smaller than that of the liquid. The bubble is
initially concentric with the tube and the liquid film thickness is constant
on the wall, $h_{0}=R_{\text{w}}-R$, where $R$ denotes the initial bubble
radius. The bubble interface can be parametrised without loss of generality in
cylindrical coordinates $(r,\theta,z)$ as $r_{\text{int}}(t,\theta,z)$ using
the tube center as the origin. It will be shown in § III.1 that it remains
radially representable at later times. The dimensionless state vector
${\mathbf{q}}=({\mathbf{u}},p,\mathcal{R}_{\text{int}})^{T}$ defines the flow
at time $t$, where ${\mathbf{u}}(t,x,y,z)=(u_{x},u_{y},u_{z})^{T}$ denotes the
three-dimensional velocity field, $p(t,x,y,z)$ denotes the pressure, and
$\mathcal{R}_{\text{int}}=r_{\text{int}}/R$ denotes the dimensionless
interface radius. The state vector and the governing equations are rendered
dimensionless by the intrinsic velocity scale associated with a viscous film
of thickness $h_{0}$ falling under its weight, inspired from Duclaux et al.
[16]. However, we choose differently the length, time, and pressure scales as
follows:

$\displaystyle\mathcal{L}$ $\displaystyle=R,$ $\displaystyle\mathcal{U}$
$\displaystyle=\frac{\rho gh_{0}^{2}}{\mu}=\frac{\rho
gR^{2}}{\mu}(\beta-1)^{2},$ (1) $\displaystyle\mathcal{P}$
$\displaystyle=\frac{\gamma}{R},$ $\displaystyle\mathcal{T}$
$\displaystyle=\frac{\mathcal{L}}{\mathcal{U}}=\frac{\mu}{\rho
gR}(\beta-1)^{-2},$

where $\beta={R_{\text{w}}}/{R}$ denotes the dimensionless tube radius. As a
result, the dimensionless value of the initial film thickness can be expressed
as $\delta=h_{0}/R=\beta-1$. The flow is governed by the incompressible
Navier-Stokes equations which in dimensionless form read

$\nabla\cdot\,{\mathbf{u}}=0,$ (2)
$\left(\frac{Bo}{Oh}\right)^{2}\delta^{4}\left(\partial_{t}+{{\mathbf{u}}}\cdot\nabla\right){\mathbf{u}}=\nabla\cdot\,\underline{\underline{\tau}}-Bo\text{
e}_{y},$ (3)

where $\partial_{j}$ denotes the partial derivative with respect to quantity
$j$, and the stress tensor $\underline{\underline{\tau}}$ reads

$\underline{\underline{\tau}}=-\,{p}\mathsfbi{I}+Bo\
\delta^{2}\left(\nabla{{\mathbf{u}}}+\nabla{{\mathbf{u}}}^{T}\right).$ (4)

The two other dimensionless numbers which appear in the governing equations
are the Ohnesorge number, $Oh={\mu}/{\sqrt{\rho\gamma R}}$, and the Bond
number, $Bo={\rho gR^{2}}/{\gamma}$. While $Oh$ compares the viscous forces to
the inertial and surface tension forces, $Bo$ compares the gravitational and
surface tension forces.

The no-slip boundary condition, $\mathbf{u}=0$, is applied at the tube wall,
$r=\beta$. At the shear-free liquid-gas interface, the dimensionless kinematic
and dynamic boundary conditions write

$\partial_{t}\mathcal{R}_{\text{int}}+\left({\mathbf{u}}\cdot\nabla\right)\mathcal{R}_{\text{int}}={\mathbf{u}\cdot
e}_{r}\quad\text{ at }\ r=\mathcal{R}_{\text{int}},$ (5)
$\underline{\underline{\tau}}\ {{\mathbf{n}}}=\kappa{{\mathbf{n}}}\quad\text{
at }\ r=\mathcal{R}_{\text{int}},$ (6)

respectively, where ${\mathbf{e}}_{r}$ denotes the unit radial vector,
${\mathbf{n}}=\nabla\left(r-\mathcal{R}_{\text{int}}\right)/\left\|\nabla\left(r-\mathcal{R}_{\text{int}}\right)\right\|$
denotes the interface unit normal vector pointing from the gas to the liquid,
$\left\|\cdot\right\|$ denotes the Euclidean norm, and
$\kappa=\nabla\cdot{{\mathbf{n}}}$ denotes the interface mean curvature.

###  II.2 Base flow

Figure 2: Base flow evolution: (a) Snapshots of the flow field: colour map
shows the velocity magnitude, arrows show the liquid velocity field, and the
point $N$ marks the north pole of the bubble, $\theta=\pi/2$; (b) Temporal
variation of the relative liquid film thickness at $N$; (c) Film thickness
variation as a function of the vertical velocity at $N$;
$Oh\rightarrow\infty,Bo=0.05,\delta/\beta=0.3$.

The base flow ${\mathbf{q}}^{0}(t,x,y)$ is the two-dimensional time-dependent
solution of the Navier-Stokes equations (2)-(6) initiated from rest at $t=0$,
with pressure $p^{0}(t=0)=-1-Bo\ y$, which consists of a constant contribution
from the interface Laplace pressure and a linear hydrostatic term. To begin
with, before even presenting the numerical tools in § II.5, we build up the
intuition about the flow by presenting an illustrative example of the base
flow evolution in figure 2 for a case where surface tension dominates gravity,
$Bo=0.05$, and a thick initial film, $\delta/\beta=0.3$, considering the
quasi-inertialess limit when $Oh\rightarrow\infty$. The liquid drains around
the core bubble under gravity, and consequently, the bubble moves upwards
under the buoyancy effect. The dynamics of the base flow, presented in figure
2(b-c), can be characterised by quantifying the liquid film thickness and the
vertical velocity at the north pole of the bubble ($N$, hereafter referred to
as the north pole), where the strongest gravitational effects are expected
[42]. In the quasi-inertialess limit, $Oh\rightarrow\infty$, drainage begins
with an immediate upward drift of the bubble (denoted as $t=0^{+}$ in figure
2(a)). Drainage and hence the rising velocity of the bubble decrease with
time, as the bubble approaches the top side of the tube.

###  II.3 Linear stability analysis

Since the base flow presented in §II.2 is time-dependent, in order to perform
the linear stability analysis, we look at the frozen frames of the unperturbed
flow, from which we investigate the evolution of small perturbations applied
to each time frame $t^{0}$, sequentially. Here, we assume that the
perturbations evolve much faster than the base flow itself, such that within
their evolution, the base flow can be approximated as quasi-steady, an
assumption to be verified a posteriori in §III.1.4. To this aim, the state
vector ${\mathbf{q}}=({\mathbf{u}},p,\mathcal{R}_{\text{int}})^{T}$ is
decomposed into the sum of the frozen base flow solution ${\mathbf{q}}^{0}$,
and the infinitesimal time-dependent perturbation
${\mathbf{q}}^{1}=\left({\mathbf{u}}^{1},p^{1},\eta^{1}\right)^{T}$ which
writes

${\mathbf{q}}={\mathbf{q}}^{0}|_{t^{0}}+\epsilon{\mathbf{q}}^{1}+\mathcal{O}(\epsilon^{2}),\quad\epsilon\ll
1,$ (7)

where the amplitude $\epsilon$ is assumed to be small. The ansatz of the
perturbation ${\mathbf{q}^{1}}$ with the longitudinal wavenumber $k$
(associated with the wavelength $\lambda={2\pi}/{k}$) reads

${\mathbf{q}}^{1}=\tilde{{\mathbf{q}}}(x,y)\ \mathrm{exp}\left[\sigma
t+\mathrm{i}kz\right]+\mathrm{c.c.},$ (8)

where c.c. denotes the complex conjugate. Any other function in terms of the
state vector can be decomposed in the same fashion, e.g.
$\underline{\underline{\tau}}=\underline{\underline{\tau}}^{0}+\epsilon\underline{\underline{\tau}}^{1}$,
${\mathbf{n}}={\mathbf{n}}^{0}+\epsilon{\mathbf{n}}^{1}$ and
$\kappa=\kappa^{0}+\epsilon\kappa^{1}$. For further details about the
formulation of ${\mathbf{n}}^{1}$ and $\kappa^{1}$, see appendix C. In the
asymptotic limit of large times, a normal eigenmode perturbation with complex
pulsation $\sigma=\sigma_{r}+\mathrm{i}\sigma_{i}$, is considered as unstable
if $\sigma_{r}>0$, i.e. if $\sigma$ lies in the unstable complex half-plane.
An unstable eigenmode grows exponentially in time with the growth rate
$\sigma_{r}$ (unless otherwise noted, the indices $r$ and $i$ denote the real
and imaginary parts of a complex number, respectively). By casting the
perturbed state (7) into the governing equations (2)-(3), and keeping the
first-order terms, the linearised equations are obtained as

$\nabla\cdot\,{\mathbf{u}}^{1}=0,$ (9)
$\left(\frac{Bo}{Oh}\right)^{2}\delta^{4}\left(\partial_{t}{\mathbf{u}}^{1}+\left({\mathbf{u}}^{0}\cdot\nabla\right){\mathbf{u}}^{1}+\left({\mathbf{u}}^{1}\cdot\nabla\right){\mathbf{u}}^{0}\right)=\nabla\cdot\,\underline{\underline{\tau}}^{1}.$
(10)

The no-slip condition implies ${\mathbf{\tilde{}}{u}=0}$ on the solid wall.
The interface conditions (5)-(6), applied on the perturbed liquid interface,
can be projected radially onto the base interface and ultimately linearised, a
process called flattening (see (25) in appendix A). The linearised kinematic
condition can be expressed as

$\sigma\tilde{\eta}+\underbrace{\left(-\partial_{r}u^{0}_{r}+\frac{\partial_{r}u^{0}_{\theta}\
\partial_{\theta}\mathcal{R}_{\text{int}}^{0}}{\mathcal{R}_{\text{int}}^{0}}-\frac{u^{0}_{\theta}\
\partial_{\theta}\mathcal{R}_{\text{int}}^{0}}{\left(\mathcal{R}_{\text{int}}^{0}\right)^{2}}\right)\tilde{\eta}+\frac{u^{0}_{\theta}}{\mathcal{R}_{\text{int}}^{0}}\partial_{\theta}\tilde{\eta}}_{-\mathsfbi{G}^{0}\tilde{\eta}}+\
\frac{\partial_{\theta}\mathcal{R}_{\text{int}}^{0}}{\mathcal{R}_{\text{int}}^{0}}\tilde{u}_{\theta}={\tilde{u}_{r}}\quad\text{
at}\ r=\mathcal{R}_{\text{int}}^{0},$ (11)

where $\left(u^{0}_{r},u^{0}_{\theta}\right)^{T}$ and
$\left(\tilde{u}_{r},\tilde{u}_{\theta}\right)^{T}$ denote the velocity
vectors of the base state and perturbations, respectively, represented in the
cylindrical coordinates. Introducing an eigenstate vector of the form (8) into
(9)-(10), combined with (11), leads to a generalised eigenvalue problem for
$\sigma$ and $\tilde{{\mathbf{q}}}$ that writes

$\mathsfbi{L}\tilde{{\mathbf{q}}}+\text{c.c.}=\sigma\mathsfbi{B}\tilde{{\mathbf{q}}}+\text{c.c.},$
(12)

where the linear operators $\mathsfbi{L}$ and $\mathsfbi{B}$ are defined as

$\displaystyle\mathsfbi{L}$
$\displaystyle=\begin{bmatrix}\left(\frac{Bo}{Oh}\right)^{2}\delta^{4}\mathsfbi{F}^{0}+Bo\delta^{2}\left(\tilde{\nabla}\cdot(\tilde{\nabla}+\tilde{\nabla}^{T})\right)&-\tilde{\nabla}&{\mathbf{0}}\\\
\tilde{\nabla}\cdot&0&0\\\
\left({\mathbf{e}}_{r}-\frac{\partial_{\theta}\mathcal{R}_{\text{int}}^{0}}{\mathcal{R}_{\text{int}}^{0}}{\mathbf{e}}_{\theta}\right)\cdot&0&\mathsfbi{G}^{0}\end{bmatrix},$
$\displaystyle\mathsfbi{B}$
$\displaystyle=\begin{bmatrix}\left(\frac{Bo}{Oh}\right)^{2}\delta^{4}\mathsfbi{I}&{\mathbf{0}}&{\mathbf{0}}\\\
{\mathbf{0}}&0&0\\\ {\mathbf{0}}&0&1\,\end{bmatrix},$ (13)

where
$\mathsfbi{F}^{0}\tilde{\mathbf{u}}=-\bigl{(}\left({\mathbf{u}}^{0}\cdot\tilde{\nabla}\right)\tilde{\mathbf{u}}+\left(\tilde{\mathbf{u}}\cdot\nabla\right){\mathbf{u}}^{0}\bigr{)}$,
$({\mathbf{e}}_{r},{\mathbf{e}}_{\theta},{\mathbf{e}}_{z})$ denote the vectors
of unit directions in the cylindrical coordinates $(r,\theta,z)$ used for
parameterising the interface, and the gradient operators and the velocity
gradient tensors in the Cartesian coordinates read

$\displaystyle{\nabla}$
$\displaystyle=(\partial_{x},\partial_{y},\partial_{z})^{T},$
$\displaystyle{\nabla}{\mathbf{u}}^{0}$
$\displaystyle=\begin{bmatrix}\partial_{x}{u}^{0}_{x}&\partial_{y}{u}^{0}_{x}&0\\\
\partial_{x}{u}^{0}_{y}&\partial_{y}{u}^{0}_{y}&0\\\ 0&0&0\,\end{bmatrix},$
$\displaystyle\tilde{\nabla}$
$\displaystyle=(\partial_{x},\partial_{y},\mathrm{i}k)^{T},$
$\displaystyle\tilde{\nabla}\tilde{{\mathbf{u}}}$
$\displaystyle=\begin{bmatrix}\partial_{x}\tilde{u}_{x}&\partial_{y}\tilde{u}_{x}&\mathrm{i}k\tilde{u}_{x}\\\
\partial_{x}\tilde{u}_{y}&\partial_{y}\tilde{u}_{y}&\mathrm{i}k\tilde{u}_{y}\\\
\partial_{x}\tilde{u}_{z}&\partial_{y}\tilde{u}_{z}&\mathrm{i}k\tilde{u}_{z}\,\end{bmatrix}.$
(14)

The interface dynamic condition (6), once linearised, can be expressed as

$\underline{\underline{\tau}}^{0}\tilde{\mathbf{n}}+\tilde{\eta}\,\partial_{r}\underline{\underline{\tau}}^{0}{\mathbf{n}}^{0}+\underline{\underline{{\tilde{\tau}}}}\
{\textbf{n}}^{0}=\kappa^{0}\tilde{\mathbf{n}}+\tilde{\kappa}{\textbf{n}}^{0}\quad\text{
at}\ r=\mathcal{R}_{\text{int}}^{0}.$ (16)

For further details on the derivation of the interface conditions and their
implementation, see appendices A and B, respectively.

###  II.4 Transient growth analysis

In contrast to the linear stability analysis, the transient growth analysis
tolerates any kind of temporal dependency for both the base flow and the
perturbation and does not rely on a separation of time scales between their
respective evolutions. Moreover, it accounts for the so-called ”nonmodal”
mechanisms arising from the nonnormality (non-commutativity with the adjoint)
of the linear operator $\mathsfbi{L}$. Owing to these mechanisms, small-
amplitude initial perturbations may experience a large transient amplification
due to an intricate cooperation between a possibly large number of eigenmodes;
therefore, reducing the dynamics to the leading (least stable or most
unstable) eigenmode might be irrelevant at finite time. If nonnormality is
often inherited from the linearization of the advective term, the operator
$\mathsfbi{L}$ expressed in (II.3) can be nonnormal even in the quasi-
inertialess limit $Oh\rightarrow\infty$, due to the presence of the interface.

The ansatz (8) is generalised as

${\mathbf{q}}^{1}=\bar{\mathbf{q}}(t,x,y)\
\mathrm{exp}\left[\mathrm{i}kz\right]+\mathrm{c.c.},$ (17)

where an exponential temporal behavior is not enforced. The temporal evolution
of $\bar{\mathbf{q}}(t,x,y)$ is an extension of (12)

$\mathsfbi{L}\bar{\mathbf{q}}+\text{c.c.}=\mathsfbi{B}\frac{\partial\bar{\mathbf{q}}}{\partial
t}+\text{c.c.},$ (18)

We recall that $\mathsfbi{L}$ is time-dependent and parameterised by $k$. We
are interested in the initial perturbation of the interface
$\bar{\mathbf{q}}(0,x,y)=[{\mathbf{0}},0,\bar{\eta}(0)]^{T}$ that is the most
amplified by (18) after a temporal horizon $t=T$, where $T$ is named a
temporal horizon. A given initial condition might project on the optimal one,
and can be amplified strongly enough to lead to a nonlinear regime. We follow
the methodology proposed in Del Guercio et al. [13], and take advantage of the
fact that the base flow interface is axisymmetric at $t=0$ (see figure 2) to
expand $\bar{\eta}(0)$ as a finite series of Fourier modes in $\theta$

$\bar{\eta}(0)=\sum_{m=-N}^{N}\alpha_{m}e^{\mathrm{i}m\theta}=a_{0}+\sum_{m=1}^{N}[a_{m}\cos(m\theta)+b_{m}\sin(m\theta)]$
(19)

with

$\alpha_{-m}=\alpha_{m}^{\star}$

where the superscript $\star$ denotes the complex conjugate,
$a_{0}=\alpha_{0}$, and for $m\geq 1$

$a_{m}=\alpha_{m}+\alpha_{m}^{\star},\quad\text{and}\quad
b_{m}=i(\alpha_{m}-\alpha_{m}^{\star}).$

Note that $a_{m}$, $b_{m}$ ($m=0,1,2,...$) and $\bar{\eta}(0)$ are real-
valued. The associated first-order interfacial energy density per spanwise
wavelength is proportional to

$e_{0}=\frac{k}{2\pi}\int_{0}^{2\pi}\int_{0}^{2\pi/k}\left[\bar{\eta}(0)e^{\mathrm{i}kz}+c.c\right]^{2}\mathrm{d}z\mathrm{d}\theta=2\pi{\bm{a}}^{T}\bm{a},$
(20)

where $\mathcal{R}_{\text{int}}^{0}(0)=1$ was implied, and where we defined

$\bm{a}=[a_{N},a_{N-1},...,a_{1},\sqrt{2}a_{0},b_{N},b_{N-1},...,b_{1}].$ (21)

Let us define $\bar{\eta}_{m}(t)$ as the evolution at time $t$ of each term
$\bar{\eta}_{m}(0)=e^{\mathrm{i}m\theta}$ of the inital condition (19). Thanks
to the linearity of the evolution equation (18), the interface shape at $t=T$
reads simply

$\bar{\eta}(T)=\sum_{m=-N}^{N}\alpha_{m}\bar{\eta}_{m}(T).$ (22)

The corresponding interfacial energy density per wavelength is proportional to

$e(T)=\frac{k}{2\pi}\int_{0}^{2\pi}\int_{0}^{2\pi/k}|\bar{\eta}(T)e^{\mathrm{i}kz}+c.c|^{2}\mathrm{d}z\mathcal{R}_{\text{int}}^{0}(T)\mathrm{d}\theta.$
(23)

We further show in appendix D that $e(T)=2\bm{a}^{T}\mathsfbi{E}(T)\bm{a}$
where $\mathsfbi{E}(T)$ is a real-valued, symmetric, strictly positive
definite matrix of size $(2N+1)\times(2N+1)$. In this manner, the (optimal)
transient gain $G(T)$ defined as

$G(T)=\max_{\bm{a}}\frac{e(T)}{e_{0}}=\pi^{-1}\frac{\bm{a}^{T}\mathsfbi{E}(T)\bm{a}}{\bm{a}^{T}\bm{a}}$
(24)

is simply the largest eigenvalue of $\mathsfbi{E}(T)$ divided by $\pi$, and
the associated eigenvector provides directly the Fourier mode coefficients of
the optimal initial condition.

###  II.5 Numerical method

Figure 3:  The numerical domain used for computing the base flow and the
linear stability analysis. $\Omega_{xy}$ denotes the liquid bulk. The
boundaries of the numerical domain are denoted by
$\partial\Omega_{xy}=\partial\Sigma_{\mathrm{w}}\cup\partial\Sigma_{\mathrm{int}}$:
$\partial\Sigma_{\mathrm{w}}$ represents the interior wall of the tube with
the radius of $\beta$, and $\partial\Sigma_{\mathrm{int}}$ represents the gas-
liquid interface. The cross-section of the interface is initially a circle of
unit radius, concentric with the tube (sketched with the dashed line).

The base flow calculation, linear stability, and transient growth analyses are
performed numerically employing the finite element package COMSOL
Multiphysics${}^{\text{TM}}$. A triangular mesh is generated on the two-
dimensional domain shown in figure 3. In the following sections, the area
increment in the bulk cross-section is denoted by $\mathrm{d}A_{\Omega_{xy}}$.
On the boundary $j$, the increment of the surface area is denoted by
$\mathrm{d}A_{\Sigma_{j}}$, and the increment of arc length is denoted by
$\mathrm{d}s$. The grid size is controlled by the vertex densities on the
boundaries $\partial\Sigma_{\mathrm{w}}$ and $\partial\Sigma_{\mathrm{int}}$.
The variational formulations of the base flow equations (2)-(6), linear
stability equations (12), and linearised Navier-Stokes equations (18) are
discretised spatially using quadratic (P2) Lagrange elements for
${\mathbf{u}}^{0}$, $\tilde{{\mathbf{u}}}$, ${\mathbf{\bar{}}{u}}$,
$\tilde{\eta}$ and $\bar{\eta}$ as well as the base flow interface geometry,
and linear (P1) Lagrange elements for $p^{0}$, $\tilde{p}$, and $\bar{p}$.
While the interface conditions are applied in the built-in modules of COMSOL
Multiphysics${}^{\text{TM}}$, the linearised condition (16) is applied by
means of the Lagrange multipliers of quadratic (P2) order. The employed
discretisation yields approximately 500’000 degrees of freedom for the base
flow, linear stability, and transient growth analyses.

The time-dependent base flow is computed using the laminar two-phase flow
module coupled with the moving mesh module. The numerical time step is set by
the Backward differentiation formula with maximum differentiation order of 2.
The solver is initialised by the Backward Euler consistent initialization with
an initial step fraction of $10^{-9}$. At each time step, Newton’s method is
used to solve the non-linear equations, where the relative tolerance for
convergence is set to $10^{-6}$. Following the computed base flow, and after
extracting the geometric characteristics of the base interface at each time
step, the generalised eigenvalue problem (12) is solved using the shift-invert
Arnoldi method (for more details on the development of the variational
formulation, the implementation of the linear stability eigenvalue problem and
their corresponding boundary conditions see appendix B).

The computation time associated with obtaining the base flow for a given set
of parameters, followed by the stability analysis for $\sim$20 values of $k$,
is of the order of a few hours on a single Intel core at 3.6 GHz. For the
transient growth analysis, the computation of the propagator matrices with
$-5\leq m\leq 5$ over 10 values of $k$ takes a few days; the calculation of
the optimal transient gain for the same $k$ values and $\sim$12 values of the
temporal horizon $T$ takes approximately one hour. Both the base flow and
stability analysis model are validated with the analytical solutions available
in the literature for the flow of liquid coating the inside of a vertical tube
(for more details about the series of validation tests see appendix E).

##  III Results

###  III.1 Linear stability analysis

In this section, we present the linear stability characteristics associated
with the drainage of a liquid film coating the inside of a horizontal circular
tube. The results of the frozen base flow linear stability analysis in terms
of different parameters are presented hereafter. First, an overview of the
stability characteristics of the flow is presented in §III.1.1. Then, the
influence of $Bo$ and $Oh$ is demonstrated in §III.1.2-III.1.3, and the
stability diagram is presented in §III.1.4. Additional comments on the effect
of $Bo$ and on the most unstable wavenumber are given in §III.1.5-III.1.6. The
validity of the frozen base flow assumption is discussed in §III.1.7.

####  III.1.1 Stability of draining film at different instants of its
evolution

Figure 4: Linear stability analysis of the base flow presented in figure 2;
(a) Three least stable eigenmodes for $t=0^{+},k=0.55$; colour map presents
the axial eigenvelocity, black arrows show the in-plane eigenvelocity field,
and black dashed line shows the superposition of the eigeninterface (with
arbitrary amplitude) onto the base interface; (b) Three-dimensional render of
the same perturbed interfaces; (c) The dispersion curve of the same three
modes at $t=0^{+}$; (d) Temporal variation of the dispersion curve of the
unstable mode (mode 1); The green line marks the large time limit;
$Oh\rightarrow\infty,\ \delta/\beta=0.3,\ Bo=0.05$.

Figure 4(a-b) shows the three least stable modes of the base flow presented in
figure 2 at the initial time, $t=0^{+}$, and figure 4(c) presents the
corresponding dispersion curves (growth rate $\sigma_{r}$ as a function of the
axial wavenumber $k$). Only one unstable mode (mode 1) is detected, whose
dispersion curves at later time instants are shown in figure 4(d). The
characteristics of this mode are symmetry with respect to the vertical mid-
plane, strong interface modulation at the bottom of the bubble, and weak
interface modulation at the top of the bubble. The fact that this mode is
unstable at $k=0$ is reminiscent of the instability properties of a purely
viscous liquid thread in absence of inertia (see Eggers & Villermaux [17]).
This mode resembles the structure of the unstable interface observed in the
experiments of Duclaux et al. [16]. The maximal growth rate of the unstable
mode initially decreases before increasing and eventually saturating to a
value inferior to the initial one (see the green line in figure 4(d)), as the
flow evolves further. The corresponding wavenumber decreases monotonously.

####  III.1.2 Effect of the Bond ($Bo$) number

Figure 5: Influence of $Bo$ on the base flow and its long-time linear stability for $Oh\rightarrow\infty,\ \delta/\beta=0.3,\ Bo=\\{0.05,0.25,0.35,0.65\\}$; (a) Evolution of the relative thickness at the north pole; dashed line present the exponential fit $H_{N}^{0}/\delta=a_{1}\exp[-t/T_{1}]+a_{2}\exp[-t/T_{2}]$ and the values of $\\{a_{1},a_{2},T_{1},T_{2}\\}$ are given in table 1; inset shows the variation of $t_{f}$ as a function of $Bo$; (b) Bubble interface at $t=t_{f}$ when the relative thickness at the north pole reaches $10\%$; dashed line shows the tube wall; (c) Dispersion curve of the unstable eigenmode, mode 1, for $Bo=\\{0.05,0.25,0.65\\}$. $Bo$ | $a_{1}$ | $T_{1}$ | $a_{2}$ | $T_{2}$
---|---|---|---|---
0.05 | 0.986 | 5.07 | 0.014 | 51.5
0.25 | 0.831 | 4.75 | 0.178 | 52.9
0.35 | 0.788 | 4.74 | 0.217 | 58.0
0.65 | 0.731 | 5.07 | 0.261 | 103
Table 1: Fitting parameters associated with the relative pole thickness
$H_{N}^{0}/\delta=a_{1}\exp[-t/T_{1}]+a_{2}\exp[-t/T_{2}]$ for different
values of $Bo$ presented in figure05(a).

In this section, the influence of $Bo$ on the base flow and its stability is
illustrated. Figure 5(a-b) show for several $Bo$ values the characteristics of
the north pole evolution until $t=t_{f}$ when the relative thickness at the
north pole, $H_{N}^{0}/\delta$, diminishes below $10\%$ for an initial thick
film $\delta/\beta=0.3$. The relative pole thickness decays exponentially with
two distinct time scales, presented in table 1. Increasing $Bo$ results in an
exponential increase of $t_{f}$ (figure 5(a)). Figure 5(b) presents for the
same flows the bubble interface at $t=t_{f}$. This figure evidences the
horizontal bubble widening as it rises towards the upper wall. Low surface
tension or large gravitational effects, i.e. high $Bo$, results in more
deviation of the interface cross-section from its initial circular shape. Note
that as $Bo$ increases, the location where the gap between the interface and
the tube is thinnest moves from the north pole to the sides.

The influence of $Bo$ on the dispersion curve of mode 1 is presented in figure
5(d). The evolution of the dispersion curves for small $Bo$ values (here
$Bo=\\{0.05,0.25\\}$) is similar to that demonstrated in §III.1.1. Increasing
$Bo$ reduces the maximal growth rate of the unstable mode and the range of
unstable wavenumbers, $k$. Further increasing $Bo$ results in the suppression
of the instability after a finite time (see $Bo=0.65$). In this investigation,
no other unstable or marginally stable mode was detected while varying $Bo$.

Figure 6: Influence of $Oh$; (a) Evolution of the relative thickness at the
north pole; (b) The relative film thickness variation as a function of the
velocity at the north pole; (c) Interface at $t=t_{f}$ when the relative
thickness at the north pole reaches $10\%$; dashed line ($--$) shows the tube
wall; (d) Dispersion curve at $t=0^{+}$; (e) Dispersion curve at $t=t_{f}$;
$Bo=0.05,\ \delta/\beta=0.3,\ Oh=\\{0.01,0.05,0.25,\infty\\}$ (corresponding
to $0\leq({Bo}/{Oh})^{2}\delta^{4}\leq\mathcal{O}(1)$).

####  III.1.3 Effect of the Ohnesorge ($Oh$) number

In the present study, all of the aforementioned results were limited to the
quasi-inertialess regime, where
$\left({Bo}/{Oh}\right)^{2}\delta^{4}\rightarrow 0$. One can interpret
$\left({Bo}/{Oh}\right)^{2}\delta^{4}=Re\,\delta^{2}$ where the Reynolds
number, $Re$, is defined upon the same length scale and drainage velocity
scale given in (1). In this section, we study the base flow and its stability
in presence of inertia by decreasing $Oh$ while keeping the rest of the flow
parameters identical to those presented in §II.2 and §III.1.1. Figure 6(a-b)
depict the influence of inertia on the base flow evolution from the initial
time $t=0^{+}$ to a terminal time $t=t_{f}$ identical to that presented in
§III.1.2. Compared to the quasi-inertialess drainage, the inertial bubble
accelerates and decelerates with some delay. Eventually, the relative film
thickness at the north pole and its velocity converge to that of the quasi-
inertialess flow. The delays in the bubble acceleration and deceleration
increase by decreasing $Oh$. However, the bubble interfaces in all cases are
identical at the terminal time (figure 6(c)).

Figure 6(d-e) present the effect of $Oh$ on the dispersion curve of the
unstable mode at $t=0^{+}$ and $t=t_{f}$, respectively. Including inertia
(finite $Oh$) affects the initial dispersion relation of the unstable mode in
two ways. The first effect is associated with $k=0$, where the immediate
bubble drift in the inertialess flow was seen to induce a non-zero growth
rate. The addition of inertia, no matter small, yields a zero growth rate, a
regularisation by inertia also observed for the Rayleigh-Plateau instability
of a liquid thread (see Eggers & Villermaux [17]). The second effect is linked
to the shape of the dispersion curve, which is similar to that of the
classical thin-film Rayleigh-Plateau instability in the vicinity of a wall and
in the absence of gravity, as obtained by a lubrication analysis
$\sigma_{r}\propto k^{2}-k^{4}$ [10, 20], for $Oh>0.25$. In the limit of small
$Oh$, inertia diminishes the growth rate of the perturbations without altering
the instability cut-off compared to the inertialess limit.

Figure 7: Linear stability stability diagram; (a) in $\delta$ vs $Bo$ plane;
the blue and red bullets are obtained from the present numerical stability
analysis; three base interfaces at $t=t_{f}$ are shown for $\delta/\beta=0.3$
(corresponding to $\delta=0.43$) and $Bo=\\{0.05,0.25,0.65\\};$ (b) the same
diagram shown in $\delta$ vs $Bo/(1-Bo)^{2}$ plane; the green dashed line
shows the analytical stability margin predicted by Duclaux et al. [16], the
black squares and the passing dashed line show their transition experimental
data and their best linear fit, $\delta=0.26Bo/(1-Bo)^{2}$, respectively; the
crosses show their experimental data; the continuous black line shows the best
linear fit to the present linear analysis: $\delta=0.14Bo/(1-Bo)^{2}$;

####  III.1.4 Asymptotic linear stability diagram

The $\\{\delta,Bo\\}$ space is investigated to follow the transition of the
interface from unstable to stable. Figure 7(a) presents the stability diagram
based on the dispersion curve of the most unstable mode obtained from the
linear stability analysis at large time $t=t_{f}$ (it is verified that
evaluation of the stability at $t=\\{0.75t_{f},7t_{f}\\}$ does not alter the
obtained marginal curve). For each $\delta$, increasing $Bo$ results in
interface stabilisation. Figure 7(b) presents the same linear stability phase
diagram in $\delta$ versus $Bo/(1-Bo)^{2}$ plane. Such an abscissa is
suggested by Duclaux et al. [16] who proposed an analytical dispersion
relation using the lubrication theory and assuming a circular base interface.
Note that their definition of $Bo$ is different from ours and the abscissa is
adapted to the present study. Even though their analysis suggests the
stability threshold as $\delta=12Bo/(1-Bo)^{2}$, their experiments suggest
$\delta=0.26Bo/(1-Bo)^{2}$ for the stability transition. The present numerical
stability analysis suggests the best fit to the separatrix as
$\delta=0.14Bo/(1-Bo)^{2}$. Therefore, our stability analysis yields a
significant improvement, although a discrepancy remains. Note that the
definition by Duclaux et al. [16] of the transition regime is not the
appearance of growing perturbations but rather that of a visible top/down
asymmetry. It is therefore reassuring to observe their transition points
(black squares) fall in the unstable region predicted by our stability
analysis (see figure 7(b)). However, for a few points, the linear analysis
predicts instability while the experiments report a stable interface,
resulting in a small relative error in $\delta$.

####  III.1.5 Why does increasing Bond number stabilise the flow?

Figure 8: Decoupled effects of $Bo$ on the linear stability of the flow at
different scales for $\delta/\beta=0.3$: (a) the maximal growth rate obtained
by isolating $Bo$ for the base flow ($Bo_{b}=\\{0.01,0.65\\}$) and varying it
for perturbations ($Bo_{p}$); (b) the maximal growth rate obtained by varying
$Bo$ for the base flow ($Bo_{b}$) and isolating it for the unstable
perturbations ($Bo_{p}=\\{0.01,0.65\\})$; the varying base interfaces are
shown; In both panels, thin dashed lines show the tube wall and thin solid
lines show the base interfaces, the thick solid lines represent a stable
regime, and pentagrams mark the stabilisation thresholds, and the green colour
presents the analytical prediction by Duclaux et al. [16].

To shed light on the reasons for the present improvement in the prediction of
the stability threshold, we recall from §II.2 that the present study exhibits
for a high $Bo$ a remarkable deformation of the rising bubble (also shown in
figure 7(a)). The interface deformation is disregarded in the analysis of
Duclaux et al. [16]. Yet, their analysis predicts the suppression of the
instability. To elucidate the stabilising role of the Bond number, we take
benefit of the linear instability formalism to decouple artificially the
effect of $Bo$ into (i) its capacity to deform the base state as a result of
drainage and capillary forces, and (ii) its explicit role in the perturbation
equations (12) and (II.3). As detailed below, we will therefore distinguish
the ”perturbation Bond number ” $Bo_{p}$ and the ”base Bond number ” $Bo_{b}$.

  1. 1.

First, we isolate the effect of the perturbation Bond $Bo_{p}$ on the
stability properties by fixing the base state according to two extreme values
of the base Bond number $Bo_{b}$. In other words, we artificially prevent the
Bond number from deforming the base state but let it act on the perturbations.
We consider the two limits of low ($Bo_{b}=0.01$) and high ($Bo_{b}=0.65$)
base Bond numbers for $\delta/\beta=0.3$. The former is unstable and the
latter stable according to our linear analysis (figure 7), and the
stabilisation threshold is $Bo\approx 0.575$. Fixing the bubble interface to
that obtained at $t_{f}$, the effect of $Bo_{p}$ is investigated by varying it
over a wide range. Remember that the base flow at $t_{f}$ is quasi-stagnant.
The variation of the maximal growth rate, $\sigma_{r}^{max}$, is presented in
figure 8(a) along with its corresponding prediction from the dispersion
relation of Duclaux et al. [16] (green) and the ”matching” value obtained when
$Bo_{b}=Bo_{p}$ (black squares). The results show that the growth rate for the
bubble with a circular cross-section ($Bo_{b}=0.01$, blue circles) diminishes
with $Bo_{p}$. However, the interface does not stabilise (the growth rate does
not cross zero) even for a rather large perturbation Bond number,
$Bo_{p}=10^{2}$. In contrast, the deformed interface of $Bo_{b}=0.65$ (red
crosses) stabilises as $Bo_{p}$ exceeds a threshold (red pentagram) inferior
to that of the ”matching” case $Bo_{b}=Bo_{p}$ (black pentagram), and far
larger than that of Duclaux et al. [16] (green pentagram). Moreover, for small
$Bo_{b}$ where the bubble interface is almost circular, the isolated growth
rate (blue circle) is close to the ”matching” value (black square). We observe
that the present linear analysis does not agree with the growth rate predicted
by Duclaux et al. [16] even in the low $Bo$ limit, where the interface is
circular. It is seen that the base bubble deformation state plays a
significant role: a more deformed bubble (larger $Bo_{b}$) is less unstable so
that a lower $Bo_{p}$ number is needed to eventually stabilise the flow than
when $Bo_{p}=Bo_{b}$. Remember that, in contrast, a quasi-cylindrical bubble
(smaller $Bo_{b}$) is found more unstable than the ”matching” case
$Bo_{p}=Bo_{b}$, and is actually never stabilised. This demonstrates the
essential role of the perturbation Bond number $Bo_{p}$.

  2. 2.

Next, we artificially fix the perturbation Bond number $Bo_{p}$ to the same
extreme values ($Bo_{p}=0.01$ and $Bo_{p}=0.65$) and let the base Bond number
$Bo_{b}$ influence the base state. Not surprisingly, a larger perturbation
Bond number ($Bo_{p}=0.65$, red crosses) is more stabilising than the
”matching” case (black squares), while a smaller perturbation Bond number
($Bo_{p}=0.01$, blue circles) is more destabilising. The associated critical
$Bo_{b}$ numbers (blue and red pentagrams) bracket the ”matching” critical
Bond number (black pentagram). This demonstrates the role of the Bond number
on the base state and, in turn, on the stability properties.

To conclude, decoupling the effects of $Bo$ on the base flow and perturbations
suggests that its influence on the base state (as revealed by the bubble
deformation at $t_{f}$) is essential to stabilise the flow. Without base state
deformation, the instability would persist at large Bond numbers.

This remark contrasts with what Duclaux et al. [16] concluded, namely that the
whole flow would become stable as soon as any region would become stable (in
this case the north pole). The mismatch between their analytical threshold and
experiments (as well as our analysis) could therefore be possibly explained by
the circular bubble cross-section assumption in their derivation of the linear
stability characteristics. Another source of discrepancy may be related to
their dispersion relation, where the growth rate is non-uniform in $\theta$,
whereas in principle the normal mode assumption forbids any dependency of the
growth rate on the spatial coordinates.

####  III.1.6 Maximal unstable wavenumber

The maximal wavenumber, associated with the maximal growth rate of the
unstable mode 1 is shown in figure 9 for a wide range of parameters.
Increasing $Bo$ decreases the optimal wavelength and $k_{max}\rightarrow 0$ on
the stability boundary. While Duclaux et al. [16] propose the fit
$k_{max}^{2}=0.5-1.25Bo$, the present linear analysis evidences a similar
trend but with different fitting parameters that depend on the initial film
thickness. These fitting parameters are presented in figure 9(b). Note that
considering the universal relation of Duclaux et al. [16], one expects a
unique stabilisation threshold of $Bo=0.4$ for all film thickness values in
contrast to their phase diagram.

Figure 9: Maximal wavenumber ($k_{max}$) obtained from the linear stability
analysis: (a) $k_{max}^{2}$ vs $Bo$ for different initial film thickness
values; the black dashed lines show the linear fit to the present numerical
data, $k_{max}^{2}=c_{0}(\delta)-c_{1}(\delta)Bo$, and the green dashed line
shows the prediction by Duclaux et al. [16]; (b) the prefactors obtained from
the fit shown in panel (a).

####  III.1.7 Validity of the frozen frame assumption

Figure 10: Frozen frame analysis validity check for the unstable modes
presented in figure 5(c): (a) Temporal evolution of $1/\max(\sigma_{r})$; (b)
Comparison of the evolution rate between the base and perturbed flows.

We should recall from §II.3 that the key assumption behind the frozen frame
approach for the linear stability analysis is that the perturbations evolve
much faster than the base flow [39]. To assess the validity of this
assumption, one can consider the characteristic time scales of the base flow
evolution presented in table 1. After the initial bubble drift, the base flow
decelerates mainly with two distinct time scales; a fast time scale $T_{1}\sim
5$ and a slow time scale $T_{2}>50$. On the other hand, the relevant time
scale for the unstable perturbation’s evolution can be computed at different
time steps as $1/\max(\sigma_{r})$, which is presented in figure 10(a). To
indicate how fast the base flow interface evolves compared to the unstable
perturbation,
$|\partial_{t}H_{N}^{0}/\partial_{t}{\eta^{1}}|\sim\frac{\delta\left(a_{1}/T_{1}\exp[-t/T_{1}]+a_{2}/T_{2}\exp[-t/T_{2}]\right)}{\max(\sigma_{r})\exp\left[\max(\sigma_{r})t\right]}$
is presented in figure 10(b). The numerator is calculated from the double-
exponent fit of the base flow, and the denominator is calculated from the
temporal growth of the most unstable mode. A look at this indicator reveals
that the base flow and perturbations may evolve at similar rates at small
times, and even in some cases that the base flow may evolve faster than the
perturbation. Nevertheless, at large times, when the drainage has decayed, the
unstable perturbation grows faster than the base flow. Thus, the frozen frame
assumption is only valid for large times, and loses its rigor in the early
stages. We recall that the criterion to distinguish between stable and
unstable modes in figure 7 is the growth rate of the modes at large times,
sufficient for the frozen frame assumption to hold. However, it is crucial to
properly take into account the base flow evolution at early times and to
account for possible nonmodal mechanisms. For this reason, we carry out a
transient growth analysis whose results are presented and discussed in the
next section.

###  III.2 Transient growth analysis

Figure 11: The parameters $\delta/\beta=0.3$ and $Bo=0.65$ have been selected.
(a) Isocontours of the transient gain $G$ defined in (24) in the $k-T$ plane,
the temporal horizon $T\in[0,t_{f}]$ (only the region where $G\geq 1$ is
shown). (b) Along the thin dashed lines in (a), comparison of the transient
gain (full line) with its approximation using the dispersion curve only
(dotted line), for $k=0.2$ and $k=0.5$ (maximum gain highlighted by a circle
and a diamond, respectively).

Figure 12: (a) Transient gains for $\delta=0.43$ (corresponding to
$\delta/\beta=0.3$), and $Bo=0.05$ (significantly unstable), $Bo=0.55$
(slightly unstable), $Bo=0.575$ (stable) and $Bo=0.65$ (stable). Bullet (resp.
diamond) markers highlight the largest gain for linearly stable (resp.
unstable) parameters. Inset: zoom highlighting the largest gains for
$Bo=0.575$ and $Bo=0.65$. The wavenumber $k$ leading to the largest growth has
been selected for all curves. (b) Phase diagram obtained from the transient
growth analysis. Colours stand for $G_{\text{max}}=\max_{k,T}G$ with
$k\in[0.1,0.9]$ and $T\in[0,t_{f}]$. The thin dotted line corresponds to
$\delta/\beta=0.3$. The blue (resp. red) dashed-dotted line links the last
unstable (resp. first stable) points for a given $\delta$ and by increasing
$Bo$.

We present in what follows the results of the transient growth analysis,
conducted mainly in the stable part of the parameter space in figure 7.
Indeed, a flow can be linearly (asymptotically) stable but transiently amplify
an initial perturbation so significantly that it is destabilised through a
subcritical transition.

The transient gain for $\delta/\beta=0.3$ and $Bo=0.65$ is shown in figure
11(a) as a function of the wavenumber $k$ and the temporal horizon $T$. Longer
waves are subject to larger maximum transient gain but require more time
before reaching it. The maximum attainable transient gain over all considered
$k$ and $T$ is $G_{\text{max}}=2.64$, reached for the lowest considered $k$.
This maximum gain value is of order unity, such that, for the considered
parameters at least, a small-amplitude initial perturbation is not expected to
be sufficiently amplified by the flow to trigger nonlinearities.

The extent to which the transient growth is driven by nonmodal mechanisms, and
may or may not be estimated directly from dispersion curves is investigated in
figure 11(b), where the transient gains are compared with the integrated
growth rate $\exp\left[\int_{0}^{T}2\sigma_{r}(k,t)dt\right]$ for $k=0.2$ and
$k=0.5$ (thin dashed lines). We recall that the growth rate $\sigma_{r}(k,t)$
for $\delta/\beta=0.3$ and $Bo=0.65$ is shown in figure 5(c). For $k=0.5$ in
figure 11(b), the curves are qualitatively similar but present some small
quantitative discrepancies: the maximum transient gain reaches a slightly
larger value at a larger temporal horizon, $T=5$ against $T=2.5$. Thus,
between these two times, the energy of the perturbation can grow although the
frozen base flow is linearly stable. For $k=0.2$ the transient gain is also
systematically larger, however the two curves present little differences,
implying this transient gain to be driven mostly by modal mechanisms.

In figure 12(a), the gain for $\delta/\beta=0.3$ and $Bo=0.65$ is compared to
those of three additional Bond numbers: $Bo=0.575$, $Bo=0.55$ and $Bo=0.05$,
which correspond to a stable, a slightly and a strongly unstable flow,
respectively. The wavenumber associated with the largest gain was selected.
For $Bo=0.05$, the maximum growth rate remains large (see figure 5(c)) and the
gain quickly aligns on the fast exponential growth; specifically, the slope of
the line marked with the black diamond is two times the maximum growth rate.
For $Bo=0.55$, the maximum growth rate rapidly converges towards very small
values but does not become negative, thus the gain also grows exponentially
for large times. Overall, results shown in figure 12(a) suggest that large
gains can only be attained through instability.

To assess the generality of this conclusion, we report
$G_{\text{max}}\doteq\max_{k,T}G$ for $k\in[0.1,0.9]$ and $T\in[0,t_{f}]$ in
figure 12(b), in the same $\\{\delta,Bo\\}$ parameter space as in figure 7.
The considered points were voluntarily chosen in the stable regime and the
closest to the separatrix (i.e., the first red dot for a given $\delta$ in
figure 7(b)). For all the considered parameters, $G_{\text{max}}$ does not
exceed the small value of $G_{\text{max}}=3.8$. In this sense, the linear
stability analysis seems to provide a sufficiently complete description of
this flow, since linear stability seems to imply small transient gains;
conversely, linear instability seems to be a necessary and sufficient
condition for large gains. Moreover, as we saw for $\delta/\beta=0.3$ and
$Bo=0.65$, and we have checked these conclusions to be true for all the other
points in figure 12(b), $G_{\text{max}}$ is reasonably well characterised
using solely the information about the growth rate. In conclusion, nonmodal
mechanisms are believed to have little influence on the flow, at least for the
parameters considered in this paper.

##  IV Summary and conclusion

In this work, we have studied the draining flow of a viscous liquid film
coating the inner wall of a horizontal tube. First, the temporal evolution of
the axially invariant base flow with an initially stagnant film of uniform
thickness was computed numerically. In the absence of inertia, the base flow
exhibits an instantaneous upward drift of the core bubble, followed by an
exponential decay of the flow as the interface approaches the tube upper wall.

Next, the stability of the evolving base flow was investigated by means of a
linear stability analysis under the frozen frame assumption. It was verified
that this assumption holds at large times. One unstable mode was observed
which features the characteristics of the Rayleigh-Plateau mode: axial
interface undulations, horizontal symmetry, and vertical asymmetry evidencing
stronger modulation at the bottom. These features match the interface shape
observed in the experiments of Duclaux et al. [16]. The maximal growth rate of
this mode decreases as the bubble approaches the tube upper wall.

A parametric study was then conducted in the space of dimensionless parameters
$\\{Bo,Oh,\beta\\}$. This study suggests that increasing $Bo$, equivalent to
weakening surface tension in comparison with gravity, results in a
horizontally wider bubble interface that rises more slowly than quasi-
cylindrical bubbles observed at small Bond numbers. The bubble slowdown
coincides with the displacement of the location of minimum film thickness from
the north pole to the sides. With increasing $Bo$, the deformed interface
becomes less unstable, and the Rayleigh-Plateau instability is suppressed
above a critical $Bo$ value. We also demonstrated that inertial forces,
achieved at finite $Oh$, affect the flow mainly at the onset of the drainage
such that the bubble drifts upward more smoothly, and the drainage is delayed.
However, inertia does not alter the shape of the interface, nor its linear
stability regime, and has a minor influence on the maximal growth rate at
large times.

Finally, a stability diagram was sketched by investigating the linear
stability of the deformed interface in the limit of large times for various
$\\{Bo,\beta\\}$, confirming the stabilising effect of larger Bond numbers and
thinner films. By relaxing assumptions of past studies (e.g. circular base
interface), the present linear study showed an interesting improvement on
existing theoretical results for the transition between stable and unstable
interfaces. A slight discrepancy with the experimental data of Duclaux et al.
[16] remains, resulting in a few stable experimental conditions predicted to
be unstable by our linear stability analysis. Comparing the experimental
protocol of Duclaux et al. [16] with the assumptions of the present study
suggests some possibilities for the observed slight mismatch. Firstly, the
most unstable wavenumber diminishes with $Bo$, which results in the promotion
of very long wavelengths. Capturing wavelengths longer than the finite length
of the experimental apparatus is not possible, which may result in the
stabilization of perturbations that would be unstable in an infinitely long
tube. Additionally, in the experiments, the gaseous core is pushed into a wet
tube by a syringe pump; any residual axial velocity, left-right asymmetry, or
axial variations in the deposited film thickness [6] may all have some
influence on the stability. Finally, we note that linear analysis cannot
capture the influence of non-linearities on the final pattern. Non-linearities
may quickly saturate the growth of the linear modes to a very small amplitude
[23] which may make it difficult to observe the instability experimentally.

Lastly, a transient growth analysis was conducted so as to relax the frozen
base flow assumption. It conclusively demonstrated the small importance of
nonmodal mechanisms, both because the transient gains were systematically of
order unity in the linearly stable region of the parameter space, and because
they were rather well predicted from the growth rate of the leading
eigenvalue. This justifies the relevance of the linear stability analysis for
this flow when compared with experimental data.

The present study, based on the full Navier-Stokes equations and capable of
handling complex interface geometries, paves the way for a wealth of future
investigations. For instance, it would be of interest to investigate the
effect of inclining the tube: as gravity becomes non-orthogonal to the tube
axis, one can expect a competition between transverse drainage and
longitudinal advection, possibly resulting in a transition from absolute to
convective instability. Another exciting perspective is that of a film coating
the outside of an inclined tube: in addition to the aforementioned
competition, this configuration offers the possibility of rich nonlinear
dynamics such as pinch-off. Finally we recall that the present study has
focused on the regime of small and intermediate Bond numbers, where the
Rayleigh-Taylor instability is suppressed [42, 3]. A natural extension of our
work should explore larger Bond numbers. We conjecture the existence of an
interval of stable Bond numbers, with smaller $Bo$ unstable to the Rayleigh-
Plateau instability (like in the present study) and larger $Bo$ unstable to
the Rayleigh-Taylor instability [42, 3, 5].

##  Appendix A Derivation of the interface boundary conditions

In this section, the derivation of the interface boundary conditions for the
perturbed flow is elaborated. These conditions should be imposed on the
perturbed interface, i.e. on
$r=\mathcal{R}_{\text{int}}^{0}+\epsilon\eta^{1}$, while $\eta^{1}$ is already
a part of the problem unknowns. By using the Taylor expansion, that is,
projecting radially on the base interface, i.e. on
$r=\mathcal{R}_{\text{int}}^{0}(\theta,t)$, any flow quantity at the perturbed
interface can be readily approximated. This projection is referred to as
flattening and for an arbitrary function $f(r,\theta,z,t)$ can be expressed as

$f|_{(r=\mathcal{R}_{\text{int}}^{0}+\epsilon\eta^{1},\theta,z,t)}=f|_{(r=\mathcal{R}_{\text{int}}^{0},\theta,z,t)}+\epsilon\eta^{1}\
\partial_{r}f|_{(r=\mathcal{R}_{\text{int}}^{0},\theta,z,t)}+\mathcal{O}(\epsilon^{2}).$
(25)

By substituting the decomposed state vector of (7), into the interface
conditions (5)-(6), then using the ansatz of (8), and applying the
aforementioned flattening, we can formulate these conditions as a set of
equivalent constraints on the boundary of the base interface. The linearised
form of the kinematic condition (5) writes

$\partial_{t}\left(\mathcal{R}_{\text{int}}^{0}+\epsilon\eta^{1}\right)+\left({\mathbf{u}}^{0}+\epsilon{\mathbf{u}}^{1}\right)\cdot\nabla\left(\mathcal{R}_{\text{int}}^{0}+\epsilon\eta^{1}\right)=\left({\mathbf{u}}^{0}+\epsilon{\mathbf{u}}^{1}\right){\mathbf{\cdot}e}_{r}\quad\text{
at }\ r=\mathcal{R}_{\text{int}}^{0}+\epsilon\eta^{1},$ (26)

where the gradient vector in the Cylindrical coordinate can be expressed as
$\nabla=\left(\partial_{r},1/r\partial_{\theta},\partial_{z}\right)^{T}$.
Applying (25) to (26) and using the ansatz of (8) readily results in (11).

The linearised dynamic condition (6) writes

$\left(\underline{\underline{\tau}}^{0}+\epsilon\underline{\underline{\tau}}^{1}\right)\
\left({\mathbf{n}}^{0}+\epsilon{\mathbf{n}}^{1}\right)=\left(\kappa^{0}+\epsilon\kappa^{1}\right)\left({\mathbf{n}}^{0}+\epsilon{\mathbf{n}}^{1}\right)\quad\text{
at }\ r=\mathcal{R}_{\text{int}}^{0}+\epsilon\eta^{1},$ (27)

Applying (25) to (27) and using the ansatz of (8) readily results in (16). In
order to express interface conditions in the Cartesian coordinates, the terms
which are expressed in the Cylindrical coordinates should be transformed by
employing the Jacobian transformations as

$\displaystyle{\mathbf{e}}_{r}$ $\displaystyle=\cos{\theta}\
{\mathbf{e}}_{x}+\sin{\theta}\ {\mathbf{e}}_{y},$
$\displaystyle{\mathbf{e}}_{\theta}$ $\displaystyle=-\sin{\theta}\
{\mathbf{e}}_{x}+\cos{\theta}\ {\mathbf{e}}_{y},$ $\displaystyle\partial_{r}$
$\displaystyle=\cos{\theta}\ {\partial}_{x}+\sin{\theta}\ {\partial}_{y},$
$\displaystyle\partial_{\theta}$
$\displaystyle=\frac{{\mathbf{t}}^{0}\cdot\nabla_{s}}{{\mathbf{t}}^{0}\cdot\nabla_{s}\theta},$
(28)

where ${\mathbf{t}}^{0}$ denotes the unit tangent vector, and
$\nabla_{s}=\nabla-{\mathbf{n}}^{0}\left({\mathbf{n}}^{0}\cdot\nabla\right)$
is the tangential derivative on the base interface. Both conditions (11) and
(16) include the normal vector and the curvature of the perturbed interface
whose formulation is given in appendix C. For further details concerning the
numerical implementation of the boundary conditions, see appendix B.

##  Appendix B Variational formulation of the linear stability analysis and
implementation of boundary conditions

Implementation of the numerical scheme and development of the variational
formulation associated with the governing equations presented in section II
are elaborated in this appendix, recalling that the numerical domain is shown
in figure 3. To develop the variational form of (12), firstly the normal mode
of (8) is applied to the system of equations (9)-(11). Then it is internally
multiplied by the vector of the test functions
$\psi=[\psi_{p},\psi_{{\mathbf{u}}},\psi_{\eta}]$, where
$\psi_{{\mathbf{u}}}=[\psi_{u_{x}},\psi_{u_{y}},\psi_{u_{z}}]$. The resulting
scalar product is integrated on $\Omega_{xy}$, which in the linear order gives

$\displaystyle\biggl{\\{}\iint_{\Omega_{xy}}\psi_{p}^{\star}\left(\tilde{\nabla}\cdot\tilde{{\mathbf{u}}}\right)\
\mathrm{d}A_{\Omega_{xy}}$ $\displaystyle+$
$\displaystyle\iint_{\Omega_{xy}}\psi_{{\mathbf{u}}}^{\star}\cdot\left(\left(\frac{Bo}{Oh}\right)^{2}\delta^{4}\left(\sigma\tilde{\mathbf{u}}+{\mathbf{u}}^{0}\cdot\tilde{\nabla}\tilde{\mathbf{u}}+\tilde{\mathbf{u}}\cdot\nabla{\mathbf{u}}^{0}\right)-\tilde{\nabla}\cdot\underline{\underline{\tilde{\tau}}}\right)\
\mathrm{d}A_{\Omega_{xy}}$ $\displaystyle+$
$\displaystyle\int_{\partial\Sigma_{\mathrm{int}}}\psi_{\eta}^{\star}\left(\sigma\tilde{\eta}+\left(-\partial_{r}u^{0}_{r}+\frac{\partial_{r}u^{0}_{\theta}\
\partial_{\theta}\mathcal{R}_{\text{int}}^{0}}{\mathcal{R}_{\text{int}}^{0}}-\frac{u^{0}_{\theta}\
\partial_{\theta}\mathcal{R}_{\text{int}}^{0}}{\left(\mathcal{R}_{\text{int}}^{0}\right)^{2}}\right)\tilde{\eta}+\frac{u^{0}_{\theta}}{\mathcal{R}_{\text{int}}^{0}}\partial_{\theta}\tilde{\eta}+\
\frac{\partial_{\theta}\mathcal{R}_{\text{int}}^{0}}{\mathcal{R}_{\text{int}}^{0}}\tilde{u}_{\theta}-{\tilde{u}_{r}}\right)\
\mathrm{d}s\biggr{\\}}$ $\displaystyle+$ $\displaystyle\text{c.c.}=0.$ (29)

It should be noted that in a complex system, the applied scalar product is
Hermitian, defined as
$\left\langle{\mathbf{a},b}\right\rangle={\mathbf{a}^{\star}\cdot b}$ where
the superscript $\star$ denotes the complex conjugate. In the last line of
this system of equations, kinematic condition (11) is used to define
$\tilde{\eta}$ only on $\partial\Sigma_{\mathrm{int}}$. After integrating by
part,
$\psi_{{\mathbf{u}}}^{\star}\cdot\left(\tilde{\nabla}\cdot\underline{\underline{\tilde{\tau}}}\right)=\tilde{\nabla}\cdot\left(\underline{\underline{\tilde{\tau}}}\psi_{{\mathbf{u}}}^{\star}\right)-tr\left(\underline{\underline{\tilde{\tau}}}^{T}\left(\tilde{\nabla}\psi_{{\mathbf{u}}}\right)^{\star}\right)$
, and then applying the Gauss’s theorem,
$\iint_{\Omega_{xy}}\tilde{\nabla}\cdot\left(\underline{\underline{\tilde{\tau}}}\psi_{{\mathbf{u}}}^{\star}\right)\
\mathrm{d}A_{\Omega_{xy}}=\int_{\partial\Omega_{xy}}\left(\underline{\underline{\tilde{\tau}}}\psi_{{\mathbf{u}}}^{\star}\right)\cdot{\mathbf{n}}^{0}\
\mathrm{d}s$, (B) implies

$\displaystyle\biggl{\\{}\iint_{\Omega_{xy}}\psi_{p}^{\star}\left(\tilde{\nabla}\cdot\tilde{{\mathbf{u}}}\right)\
\mathrm{d}A_{\Omega_{xy}}$ $\displaystyle+$
$\displaystyle\iint_{\Omega_{xy}}\psi_{{\mathbf{u}}}^{\star}\cdot\left(\left(\frac{Bo}{Oh}\right)^{2}\delta^{4}\left(\sigma\tilde{\mathbf{u}}+{\mathbf{u}}^{0}\cdot\tilde{\nabla}\tilde{\mathbf{u}}+\tilde{\mathbf{u}}\cdot\nabla{\mathbf{u}}^{0}\right)\right)\
\mathrm{d}A_{\Omega_{xy}}$ $\displaystyle+$
$\displaystyle\iint_{\Omega_{xy}}tr\left(\underline{\underline{\tilde{\tau}}}^{T}\left(\tilde{\nabla}\psi_{{\mathbf{u}}}\right)^{\star}\right)\
\mathrm{d}A_{\Omega_{xy}}$ $\displaystyle+$
$\displaystyle\int_{\partial\Omega_{xy}}-\left(\underline{\underline{\tilde{\tau}}}\psi_{{\mathbf{u}}}^{\star}\right)\cdot{\mathbf{n}}^{0}\
\mathrm{d}s$ $\displaystyle+$
$\displaystyle\int_{\partial\Sigma_{\mathrm{int}}}\psi_{\eta}^{\star}\left(\sigma\tilde{\eta}+\left(-\partial_{r}u^{0}_{r}+\frac{\partial_{r}u^{0}_{\theta}\
\partial_{\theta}\mathcal{R}_{\text{int}}^{0}}{\mathcal{R}_{\text{int}}^{0}}-\frac{u^{0}_{\theta}\
\partial_{\theta}\mathcal{R}_{\text{int}}^{0}}{\left(\mathcal{R}_{\text{int}}^{0}\right)^{2}}\right)\tilde{\eta}+\frac{u^{0}_{\theta}}{\mathcal{R}_{\text{int}}^{0}}\partial_{\theta}\tilde{\eta}+\
\frac{\partial_{\theta}\mathcal{R}_{\text{int}}^{0}}{\mathcal{R}_{\text{int}}^{0}}\tilde{u}_{\theta}-{\tilde{u}_{r}}\right)\
\mathrm{d}s\biggr{\\}}$ $\displaystyle+$ $\displaystyle\text{c.c.}=0.$ (30)

$\underline{\underline{\tilde{\tau}}}$ is symmetric, thus
$\left(\underline{\underline{\tilde{\tau}}}\psi_{{\mathbf{u}}}^{\star}\right)\cdot{\mathbf{n}}^{0}=\left(\underline{\underline{\tilde{\tau}}}{\mathbf{n}}^{0}\right)\cdot\psi_{{\mathbf{u}}}^{\star}$.
Using the dynamic condition (16) and the fact that
$\psi_{{\mathbf{u}}}|_{\partial\Sigma_{\mathrm{f}}}=0$ (because of the no-slip
condition on the solid wall), the variational form of (12) implies

$\displaystyle\biggl{\\{}\iint_{\Omega_{xy}}\psi_{p}^{\star}\left(\tilde{\nabla}\cdot\tilde{{\mathbf{u}}}\right)\
\mathrm{d}A_{\Omega_{xy}}$ (31) $\displaystyle+$
$\displaystyle\iint_{\Omega_{xy}}\psi_{{\mathbf{u}}}^{\star}\cdot\left(\left(\frac{Bo}{Oh}\right)^{2}\delta^{4}\sigma\tilde{\mathbf{u}}\right)\
\mathrm{d}A_{\Omega_{xy}}$ (32) $\displaystyle+$
$\displaystyle\iint_{\Omega_{xy}}\psi_{{\mathbf{u}}}^{\star}\cdot\left(\left(\frac{Bo}{Oh}\right)^{2}\delta^{4}\left({\mathbf{u}}^{0}\cdot\tilde{\nabla}\tilde{\mathbf{u}}+\tilde{\mathbf{u}}\cdot\nabla{\mathbf{u}}^{0}\right)\right)\
\mathrm{d}A_{\Omega_{xy}}$ (33) $\displaystyle+$
$\displaystyle\iint_{\Omega_{xy}}tr\left(\underline{\underline{\tilde{\tau}}}^{T}\left(\tilde{\nabla}\psi_{{\mathbf{u}}}\right)^{\star}\right)\
\mathrm{d}A_{\Omega_{xy}}$ (34) $\displaystyle+$
$\displaystyle\int_{\partial\Sigma_{\mathrm{int}}}\left(\underline{\underline{\tau}}^{0}\
\tilde{\mathbf{n}}+\tilde{\eta}\,\partial_{r}\underline{\underline{\tau}}^{0}\
{\mathbf{n}}^{0}-\left(\kappa_{0}\tilde{\mathbf{n}}+\tilde{\kappa}{\mathbf{n}}^{0}\right)\right)\cdot\psi_{{\mathbf{u}}}^{\star}\
\mathrm{d}s$ (35) $\displaystyle+$
$\displaystyle\int_{\partial\Sigma_{\mathrm{int}}}\psi_{\eta}^{\star}\left(\sigma\tilde{\eta}\right)\
\mathrm{d}s$ (36) $\displaystyle+$
$\displaystyle\int_{\partial\Sigma_{\mathrm{int}}}\psi_{\eta}^{\star}\left(\left(-\partial_{r}u^{0}_{r}+\frac{\partial_{r}u^{0}_{\theta}\
\partial_{\theta}\mathcal{R}_{\text{int}}^{0}}{\mathcal{R}_{\text{int}}^{0}}-\frac{u^{0}_{\theta}\
\partial_{\theta}\mathcal{R}_{\text{int}}^{0}}{\left(\mathcal{R}_{\text{int}}^{0}\right)^{2}}\right)\tilde{\eta}+\frac{u^{0}_{\theta}}{\mathcal{R}_{\text{int}}^{0}}\partial_{\theta}\tilde{\eta}+\
\frac{\partial_{\theta}\mathcal{R}_{\text{int}}^{0}}{\mathcal{R}_{\text{int}}^{0}}\tilde{u}_{\theta}-{\tilde{u}_{r}}\right)\
\mathrm{d}s\biggr{\\}}$ (37) $\displaystyle+$ $\displaystyle\text{c.c.}=0.$
(38)

This variational equation can be readily implemented and solved in COMSOL
Multiphysics${}^{\text{TM}}$. It is sufficient to solve the first part (in
$\\{\\}$) and the c.c. is known consequently. The matrix representation of
(31)-(38) is shown in figure 13.

Figure 13: Matrix representation of the variational system (31)-(38), solved
in COMSOL Multiphysics${}^{\text{TM}}$; blue represents the implementation of
(9)-(10); white represents the implementation of the no-slip boundary
condition on the solid wall; green represents the implementation of the
dynamic boundary condition (16); a represents the implementation of the
kinematic condition (11).

##  Appendix C Characterization of an arbitrary interface

In this section, we present the geometrical characterization of an arbitrary
interface parameterised in Cylindrical coordinates. The aim is to develop the
characteristics of a three-dimensional interface, as well as the linear
perturbations applied to this interface. The key properties of interest are
the normal vector and the local curvature.

###  C.1 Normal vector

We recall that the unit normal vector of a linearly perturbed interface can be
decomposed as ${\mathbf{n}}={\mathbf{n}}^{0}+\epsilon{\mathbf{n}}^{1}$ which
reads

${\textbf{n}}^{0}=\begin{pmatrix}n^{0}_{r}\\\ n^{0}_{\theta}\\\
n^{0}_{z}\end{pmatrix}=A^{-1/2}\begin{pmatrix}1\\\
-\frac{1}{\mathcal{R}_{\text{int}}^{0}}\partial_{\theta}{\mathcal{R}_{\text{int}}^{0}}\\\
\partial_{z}{\mathcal{R}_{\text{int}}^{0}}\end{pmatrix},$ (39)

${\mathbf{n}}^{1}=\begin{pmatrix}n^{1}_{r}\\\ n^{1}_{\theta}\\\
n^{1}_{z}\end{pmatrix}=\begin{pmatrix}B_{r}\eta^{1}+C_{r}\partial_{\theta}\eta^{1}+D_{r}\partial_{z}\eta^{1}\\\
B_{\theta}\eta^{1}+C_{\theta}\partial_{\theta}\eta^{1}+D_{\theta}\partial_{z}\eta^{1}\\\
B_{z}\eta^{1}+C_{z}\partial_{\theta}\eta^{1}+D_{z}\partial_{z}\eta^{1}\end{pmatrix},$
(40)

where

$\displaystyle A$
$\displaystyle=1+\left(\frac{1}{\mathcal{R}_{\text{int}}^{0}}\partial_{\theta}{\mathcal{R}_{\text{int}}^{0}}\right)^{2}+\left(\partial_{z}{\mathcal{R}_{\text{int}}^{0}}\right)^{2},$
$\displaystyle B_{r}$
$\displaystyle=\frac{A^{-3/2}}{\left({\mathcal{R}_{\text{int}}^{0}}\right)^{3}}\left(\partial_{\theta}{\mathcal{R}_{\text{int}}^{0}}\right)^{2},$
$\displaystyle B_{\theta}$
$\displaystyle=A^{-3/2}\left(-\frac{1}{\left({\mathcal{R}_{\text{int}}^{0}}\right)^{4}}\left(\partial_{\theta}{\mathcal{R}_{\text{int}}^{0}}\right)^{3}+\frac{A}{\left({\mathcal{R}_{\text{int}}^{0}}\right)^{2}}\partial_{\theta}{\mathcal{R}_{\text{int}}^{0}}\right),$
$\displaystyle B_{z}$
$\displaystyle=-\frac{A^{-3/2}}{\left({\mathcal{R}_{\text{int}}^{0}}\right)^{3}}\left(\partial_{\theta}{\mathcal{R}_{\text{int}}^{0}}\right)^{2}\
\partial_{z}{\mathcal{R}_{\text{int}}^{0}},$ $\displaystyle C_{r}$
$\displaystyle=-\frac{A^{-3/2}}{\left({\mathcal{R}_{\text{int}}^{0}}\right)^{2}}\partial_{\theta}{\mathcal{R}_{\text{int}}^{0}},$
$\displaystyle C_{\theta}$
$\displaystyle=-\frac{A^{-3/2}}{{\mathcal{R}_{\text{int}}^{0}}}\left(1+\left(\partial_{z}{\mathcal{R}_{\text{int}}^{0}}\right)^{2}\right),$
$\displaystyle C_{z}$
$\displaystyle=\frac{A^{-3/2}}{\left({\mathcal{R}_{\text{int}}^{0}}\right)^{2}}\partial_{\theta}{\mathcal{R}_{\text{int}}^{0}}\
\partial_{z}{\mathcal{R}_{\text{int}}^{0}},$ $\displaystyle D_{r}$
$\displaystyle=-A^{-3/2}\partial_{z}{\mathcal{R}_{\text{int}}^{0}},$
$\displaystyle D_{\theta}$
$\displaystyle=\frac{A^{-3/2}}{{\mathcal{R}_{\text{int}}^{0}}}\partial_{z}{\mathcal{R}_{\text{int}}^{0}}\
\partial_{\theta}{\mathcal{R}_{\text{int}}^{0}},$ $\displaystyle D_{z}$
$\displaystyle=A^{-3/2}\left(\left(\partial_{z}{\mathcal{R}_{\text{int}}^{0}}\right)^{2}-A\right).$
(41)

For a base interface of the form $r=\mathcal{R}_{\text{int}}^{0}(t,\theta)$
and a perturbation ansatz (8), the normal vector of the base interface can be
further simplified as

${\textbf{n}}^{0}=\begin{pmatrix}n^{0}_{r}\\\ n^{0}_{\theta}\\\
n^{0}_{z}\end{pmatrix}=A^{-1/2}\begin{pmatrix}1\\\
-\frac{1}{\mathcal{R}_{\text{int}}^{0}}\partial_{\theta}{\mathcal{R}_{\text{int}}^{0}}\\\
0\end{pmatrix},$ (42)

and simplification of the linearised perturbation of the normal vector
(40)-(C.1) implies

${\mathbf{\tilde{}}{n}}=\begin{pmatrix}\tilde{n}_{r}\\\ \tilde{n}_{\theta}\\\
\tilde{n}_{z}\end{pmatrix}=\begin{pmatrix}B_{r}\tilde{\eta}+C_{r}\partial_{\theta}\tilde{\eta}\\\
B_{\theta}\tilde{\eta}+C_{\theta}\partial_{\theta}\tilde{\eta}\\\
\mathrm{i}kD_{z}\tilde{\eta}\end{pmatrix},$ (43)

where

$\displaystyle A$
$\displaystyle=1+\left(\frac{1}{\mathcal{R}_{\text{int}}^{0}}\partial_{\theta}{\mathcal{R}_{\text{int}}^{0}}\right)^{2},$
$\displaystyle B_{r}$
$\displaystyle=\frac{A^{-3/2}}{\left({\mathcal{R}_{\text{int}}^{0}}\right)^{3}}\left(\partial_{\theta}{\mathcal{R}_{\text{int}}^{0}}\right)^{2},$
$\displaystyle B_{\theta}$
$\displaystyle=A^{-3/2}\left(-\frac{1}{\left({\mathcal{R}_{\text{int}}^{0}}\right)^{4}}\left(\partial_{\theta}{\mathcal{R}_{\text{int}}^{0}}\right)^{3}+\frac{A}{\left({\mathcal{R}_{\text{int}}^{0}}\right)^{2}}\partial_{\theta}{\mathcal{R}_{\text{int}}^{0}}\right),$
$\displaystyle C_{r}$
$\displaystyle=-\frac{A^{-3/2}}{\left({\mathcal{R}_{\text{int}}^{0}}\right)^{2}}\partial_{\theta}{\mathcal{R}_{\text{int}}^{0}},$
$\displaystyle C_{\theta}$
$\displaystyle=-\frac{A^{-3/2}}{{\mathcal{R}_{\text{int}}^{0}}},$
$\displaystyle D_{z}$ $\displaystyle=-A^{-1/2}.$ (44)

###  C.2 Curvature

We recall that the local curvature of a linearly perturbed interface can be
decomposed as
${\mathbf{\kappa}}={\mathbf{\kappa}}^{0}+\epsilon{\mathbf{\kappa}}^{1}$ which
reads

$\kappa^{0}=\frac{1}{\mathcal{R}_{\text{int}}^{0}}\left(n^{0}_{r}+\partial_{\theta}n^{0}_{\theta}\right)+\partial_{z}n^{0}_{z},$
(45) $\displaystyle\kappa^{1}$
$\displaystyle=-\frac{\eta^{1}}{\left({\mathcal{R}_{\text{int}}^{0}}\right)^{2}}\left(n^{0}_{r}+\partial_{\theta}n^{0}_{\theta}\right)+\frac{1}{{\mathcal{R}_{\text{int}}^{0}}}\left(n^{1}_{r}+\partial_{\theta}n^{1}_{\theta}\right)+\partial_{z}n^{1}_{z}$
$\displaystyle=E\eta^{1}+F\partial_{\theta}\eta^{1}+G\partial_{z}\eta^{1}+D_{z}\partial^{2}_{zz}\eta^{1}+\frac{C_{\theta}}{{\mathcal{R}_{\text{int}}^{0}}}\partial^{2}_{\theta\theta}\eta^{1},$
(46)

where

$\displaystyle E$
$\displaystyle=\frac{1}{\mathcal{R}_{\text{int}}^{0}}\left(B_{r}+\partial_{\theta}B_{\theta}\right)-\frac{1}{\left({\mathcal{R}_{\text{int}}^{0}}\right)^{2}}\left(n^{0}_{r}+\partial_{\theta}n^{0}_{\theta}\right)+\partial_{z}B_{z},$
$\displaystyle F$
$\displaystyle=\frac{1}{\mathcal{R}_{\text{int}}^{0}}\left(C_{r}+\partial_{\theta}C_{\theta}+B_{\theta}+D_{\theta}\partial_{z}\right)+\partial_{z}C_{z}+C_{z}\partial_{z},$
$\displaystyle G$ $\displaystyle=B_{z}+\partial_{z}D_{z}.$ (47)

Note that the subscript $r$ in (39)-(C.2) does not imply the real part. For a
base interface of the form $r=\mathcal{R}_{\text{int}}^{0}(t,\theta)$ and a
perturbation ansatz (8), the curvature of the base interface can be further
simplified as

$\kappa^{0}=\frac{1}{\mathcal{R}_{\text{int}}^{0}}\left(n^{0}_{r}+\partial_{\theta}n^{0}_{\theta}\right),$
(48)

and simplification of the linearised perturbation of the curvature (C.2)-(C.2)
implies

$\displaystyle\tilde{\kappa}$
$\displaystyle=-\frac{\tilde{\eta}}{\left({\mathcal{R}_{\text{int}}^{0}}\right)^{2}}\left(n^{0}_{r}+\partial_{\theta}n^{0}_{\theta}\right)+\frac{1}{{\mathcal{R}_{\text{int}}^{0}}}\left(\tilde{n}_{r}+\partial_{\theta}\tilde{n}_{\theta}\right)+\mathrm{i}k\tilde{n}_{z}$
$\displaystyle=E\tilde{\eta}+F\partial_{\theta}\tilde{\eta}-k^{2}D_{z}\tilde{\eta}+\frac{C_{\theta}}{{\mathcal{R}_{\text{int}}^{0}}}\partial^{2}_{\theta\theta}\tilde{\eta},$
(49)

where

$\displaystyle E$
$\displaystyle=\frac{1}{\mathcal{R}_{\text{int}}^{0}}\left(B_{r}+\partial_{\theta}B_{\theta}\right)-\frac{1}{\left({\mathcal{R}_{\text{int}}^{0}}\right)^{2}}\left(n^{0}_{r}+\partial_{\theta}n^{0}_{\theta}\right),$
$\displaystyle F$
$\displaystyle=\frac{1}{\mathcal{R}_{\text{int}}^{0}}\left(C_{r}+\partial_{\theta}C_{\theta}+B_{\theta}\right).$

##  Appendix D Derivation of a simplified expression for the energy density of
the transient response.

The interfacial energy density per wavelength is proportional to

$\begin{split}e(T)&=\frac{k}{2\pi}\int_{0}^{2\pi}\int_{0}^{2\pi/k}|\bar{\eta}(T)e^{\mathrm{i}kz}+c.c|^{2}\mathrm{d}z\mathcal{R}_{\text{int}}^{0}(T)\mathrm{d}\theta\\\
&=2\int_{0}^{2\pi}|\bar{\eta}(T)|^{2}\mathcal{R}_{\text{int}}^{0}(T)\mathrm{d}\theta\\\
&=2\sum_{m=-N}^{N}\sum_{n=-N}^{N}\alpha_{m}^{\star}\alpha_{n}\int_{0}^{2\pi}\bar{\eta}_{m}(T)^{\star}\bar{\eta}_{n}(T)\mathcal{R}_{\text{int}}^{0}(T)\mathrm{d}\theta\\\
&=2\bm{\alpha}^{H}\mathsfbi{A}(T)\bm{\alpha}\end{split}$ (50)

where the strictly positive definite and Hermitian matrix $\mathsfbi{A}(T)$ is
such that

$[\mathsfbi{A}(T)]_{mn}=\int_{0}^{2\pi}\bar{\eta}_{m}(T)^{\star}\bar{\eta}_{n}(T)\mathcal{R}_{\text{int}}^{0}(T)\mathrm{d}\theta,$
(51)

for $-N\leq m\leq N$ and $-N\leq n\leq N$, and where we defined

$\begin{split}\bm{\alpha}=[\alpha_{N}^{\star},\alpha_{N-1}^{\star},...,\alpha_{1}^{\star},\alpha_{0},\alpha_{1},...,\alpha_{N-1},\alpha_{N}]^{T}.\end{split}$
(52)

The leading eigenvector of $\mathsfbi{A}(T)$ has no particular reason to have
its $N$ first elements equal to the complex conjugate of its last $N$. Thus,
it cannot correspond directly to the optimal set of $\bm{\alpha}$, that should
satisfy this last constraint. For this reason we introduce the matrix

$\begin{split}\mathsfbi{M}^{-1}=\begin{bmatrix}\mathsfbi{I}&\mathsfbi{O}&\mathsfbi{P}\\\
\mathsfbi{O}^{T}&\sqrt{2}&\mathsfbi{O}^{T}\\\
-\mathrm{i}\mathsfbi{I}&\mathsfbi{O}&\mathrm{i}\mathsfbi{P}\end{bmatrix},\quad\text{with}\quad\mathsfbi{I}_{ij}=\delta_{ij},\quad\mathsfbi{P}_{ij}=\delta_{N-i+1,j},\quad\text{and}\quad\mathsfbi{O}_{i}=0.\end{split}$
(53)

for $1\leq i,j\leq N$. Namely, $\mathsfbi{I}$ is the identity matrix of size
$N\times N$, $\mathsfbi{P}$ is the identity matrix mirrored around its
vertical axis (it contains ones on the diagonal from the bottom-left to the
top-right and zeros everywhere else) and $\mathsfbi{O}$ is a vectors of zeros
of size $N\times 1$. In this manner, we have directly
$\bm{a}=\mathsfbi{M}^{-1}\bm{\alpha}$, such that :

$\begin{split}e(T)&=2\bm{\alpha}^{H}\mathsfbi{A}(T)\bm{\alpha}=2\bm{a}^{T}\mathsfbi{M}^{H}\mathsfbi{A}(T)\mathsfbi{M}\bm{a}=2\bm{a}^{T}\Re[\mathsfbi{M}^{H}\mathsfbi{A}(T)\mathsfbi{M}]\bm{a}\end{split}$
(54)

where we used that $\bm{a}$ and $e(T)$ are a real-valued. Therefore, defining
$\mathsfbi{E}(T)=\Re[\mathsfbi{M}^{H}\mathsfbi{A}(T)\mathsfbi{M}]$ leads to
the desired result.

##  Appendix E Validation of the numerical model

Figure 14: Numerical model validation; (a) Base flow: Two snapshots of the
liquid film thickness, $H^{0}=\beta-\mathcal{R}_{\text{int}}^{0}$, in the
right half-plane, at time $t=\\{2\beta,4\beta\\}$. Solid lines present the
solution from the present study, and the dashed lines present the solution
obtained by Balestra et al. [3]; $Oh=10.05,Bo=1960.2,\delta/\beta=0.01$; (b)
Stability analysis: dispersion curves of the two least stable modes associated
with the viscous film coating inside a vertical tube, namely $|m|=\\{0,1\\}$.
The continuous lines present the analytical solution obtained from the Stokes
equations, and the circle represents the results from the present model;
$Oh\rightarrow\infty,Bo=1,\delta=0.1$.

The developed numerical scheme is validated hereafter. Several measures are
taken to ensure the correspondence of the model, based on the asymptotic
limits where an analytical solution may exist.

###  E.1 Base flow model

The present base flow model is validated with Balestra et al. [3] who studied
a similar base flow in the limit of thin film, $\delta\ll 1$, and small
surface tension, $Bo\gg 1$, by employing lubrication equations [33]. Figure
14(a) shows that the present model gives a solution of the base flow in full
agreement with the solution of Balestra et al. [3].

###  E.2 Linear stability analysis model

The present linear stability model is validated with the analytical solution
that Camassa et al. [10] presented for the gravity-driven flow of a viscous
film that coats the interior side of a vertical tube (where gravity points in
$z$ direction in figure 1). The corresponding base flow is parallel and can be
expressed in cylindrical coordinates as

$u^{0}_{z}=\frac{\delta^{-2}}{2}\left(\frac{\beta^{2}-r^{2}}{2}+\ln{\frac{r}{\beta}}\right),\quad
p^{0}=-1,\quad\mathcal{R}_{\text{int}}^{0}=1,$ (55)

For the linear stability analysis, Camassa et al. [10] employed an
approximation of the Stokes equations for long jets, referred to as the long-
wavelength approximation [36] and compared the results with the analytical
solution, in terms of Bessel functions, obtained by solving the full Stokes
equations [22]. Thanks to the axisymmetry of the base flow, Camassa et al.
[10] considered a perturbation as in (8) with the Fourier ansatz exponent of
$\exp[\sigma t+\mathrm{i}kz+\mathrm{i}m\theta]$ where $m$ denotes the
azimuthal wavenumber. Both of the aforementioned equations match in the limit
of thin liquid film, $\beta\rightarrow 1$, where the long-wavelength
approximation gives a dispersion relation of
$\sigma=S\left(k^{2}-k^{4}\right)-\mathrm{i}k$ with $S=S(Bo,\delta)$ being
constant for the axisymmetric perturbation $m=0$. All of the helical
perturbations, $|m|>0$, are known to be linearly stable for such a flow [37].
Figure 14(b) presents the full agreement between the present linear stability
model and the analytical solution for a relatively thin film thickness
$\delta=0.1$. It should be noted that despite the axisymmetric nature of the
validated case, this presented validation holds also for an arbitrary
interface. For this aim, the geometrical symmetry in the numerical reference
frame is broken by setting the origin of the coordinates system at an
arbitrary location inside the bubble, $(x,y)=(0.2,0.7)$.

###  E.3 Grid independency

Figure 15: Mesh convergence proof for
$Oh\rightarrow\infty,Bo=0.05,\delta/\beta=0.3$: (a) base flow;
$H_{N}^{0}/\delta$ vs $u_{y}^{0}$; (b) linear stability analysis;
$\sigma_{r}^{max}$ vs $N_{dof}$. All of the results presented in this
manuscript are obtained for M3.

A convergence study for the base flow evolution and the linear stability of
the unstable eigenvalue is presented in figure 15, for
$\\{Oh\rightarrow\infty,Bo=0.05,\delta/\beta=0.3\\}$. Mesh resolution is
controlled by setting the number of divisions on the solid wall and interface
boundaries. Mesh convergence is already attained for the presented grids. All
of the presented results in the manuscript are obtained employing M3.

## References

  * [1]
  * Augello et al. [2018] Augello, L., Fani, A. & Gallaire, F. 2018 The influence of the entry region on the instability of a coflowing injector device. J. Phys.: Condens. Matter 30, 284003, https://doi.org/10.1088/1361-648X/aac790.
  * Balestra et al. [2016] Balestra, G., Brun, P.-T. & Gallaire, F. 2016 Rayleigh-Taylor instability under curved substrates: An optimal transient growth analysis. Phys. Rev. Fluids 1(8), 083902, https://doi.org/10.1103/PhysRevFluids.1.083902.
  * Balestra et al. [2018] Balestra, G., Kofman, N., Brun, P.-T., Scheid, B. & Gallaire, F. 2018 Three-dimensional Rayleigh–Taylor instability under a unidirectional curved substrate. J. Fluid Mech. 837, 19–47, https://doi.org/10.1017/jfm.2017.817.
  * Balestra et al. [2018] Balestra, G., Nguyen, D. M.-P. & Gallaire, F. 2018 Rayleigh-Taylor instability under a spherical substrate. Phys. Rev. Fluids 3 (8), 084005, https://link.aps.org/doi/10.1103/PhysRevFluids.3.084005.
  * Balestra et al. [2018] Balestra, G., Zhu, L. & Gallaire, F. 2018 Viscous Taylor droplets in axisymmetric and planar tubes: From Bretherton’s theory to empirical models. Microfluid. Nanofluid. 22, 67, https://doi.org/10.1007/s10404-018-2084-y.
  * Benilov [2006] Benilov, E. S. 2006 Does surface tension stabilize a liquid film inside a rotating horizontal cylinder? Part 2: Multi-dimensional disturbances. Stud. Appl. Maths 116, 1–20, https://doi.org/10.1111/j.1467-9590.2005.00331.x.
  * Benilov et al. [2005] Benilov, E. S., Kopteva, N. & O’Brien, S. B. G. 2005 Does surface tension stabilize a liquid film inside a rotating horizontal cylinder?. Q. J. Mech. Appl. Maths 58, 158–-200, https://doi.org/10.1093/qjmamj/hbi004.
  * Bian et al. [2010] Bian, S., Tai, C-F., Halpern, D., Zheng, Y. & Grotberg, J.B. 2016 Experimental study of flow fields in an airway closure model. J. Fluid Mech. 647, 391–402, https://doi.org/10.1017/S0022112010000091.
  * Camassa et al. [2014] Camassa, R., Ogrosky, H. R. & Olander J. 2014 Viscous film flow coating the interior of a vertical tube. Part 1. Gravity-driven flow. J. Fluid Mech. 745, 682–715, https://doi.org/10.1017/jfm.2014.90.
  * Camassa et al. [2016] Camassa, R., Marzuola, J., Ogrosky, H. R. & Vaughn, N. 2016 Traveling waves for a model of gravity-driven film flows in cylindrical domains. Phys. D 333, 254–265, https://doi.org/10.1016/j.physd.2015.12.003.
  * Camassa et al. [2017] Camassa, R., Ogrosky, H. R. & Olander J. 2017 Viscous film-flow coating the interior of a vertical tube. Part 2. Air-driven flow. J. Fluid Mech. 825, 1056–1090, https://doi.org/10.1017/jfm.2017.409.
  * Del Guercio et al. [2014] Del Guercio, G. and Cossu, C. and Pujals, G. 2014 Optimal streaks in the circular cylinder wake and suppression of the global instability. Journal of Fluid Mechanics, 752, 572-588.
  * Ding et al. [2018] Ding, Z., Liu, R., Wong, T. N. & Yang, C. 2018 Absolute instability induced by Marangoni effect in thin liquid film flows on vertical cylindrical surfaces. Chem. Eng. Sc., 177, 261–269, https://doi.org/10.1016/j.ces.2017.11.039.
  * Dobson & Chato [1998] Dobson, M. K. & Chato, J. C. 1998 Condensation in Smooth Horizontal Tubes. ASME. J. Heat Transfer, 120(1), 193-–213, https://doi.org/10.1115/1.2830043.
  * Duclaux et al. [2006] Duclaux, V., Clanet, C. & Quéré, D. 2006 The effects of gravity on the capillary instability in tubes. J. Fluid Mech., 556, 217–226, https://doi.org/10.1017/S0022112006009505.
  * Eggers & Villermaux [2008] Eggers, J. & Villermaux, E. 2008 Physics of liquid jets. Rep. Prog. Phys., 71, 036601, http://stacks.iop.org/RoPP/71/036601.
  * Fermigier et al. [1987] Fermigier, M., Limat, L., Westfreid, J. E., Boudinet, P. & Quilliet, C 1992 Two-dimensional patterns in Rayleigh–Taylor instability of a thin layer. J. Fluid Mech., 236, 349–-383, https://doi.org/10.1017/S0022112092001447.
  * Frenkel et al. [1987] Frenkel, A. L., Babchin, A. J., Levich, B. G., Shlang, T. & Sivashinsky, G. I. 1987 Annular flows can keep unstable films from breakup: Nonlinear saturation of capillary instability. J. Colloid Interface Sci., 115 1, 225–233, https://doi.org/10.1016/0021-9797(87)90027-0.
  * Gallaire & Brun [2017] Gallaire, F. & Brun, P.-T.  2017 Fluid dynamic instabilities: theory and application to pattern forming in complex media. Phil. Trans. R. Soc. A, 375 (2093), 20160155, https://doi.org/10.1098/rsta.2016.0155.
  * Goldsmith & Mason [1963] Goldsmith, H. L. & Mason, S. G. 1963 The flow of suspensions through tubes. II. Single large bubbles. J. Colloid Interface Sci., 18(3), 237–261, https://doi.org/10.1016/0095-8522(63)90015-1.
  * Goren [1962] Goren, S. 1962 The instability of an annular thread of fluid. J. Fluid Mech., 12(2), 309–319, https://doi.org/10.1017/S002211206200021X.
  * Halpern & Grotberg [2003] Halpern, D. & Grotberg, J. 2003 Nonlinear saturation of the Rayleigh instability due to oscillatory flow in a liquid-lined tube. J. Fluid Mech., 492, 251–270, https://doi.org/10.1017/S0022112003005573.
  * Heil et al. [2008] Heil, M., Hazel, A. L. & Smith, J. A. 2008 The mechanics of airway closure. Respir. Physiol. Neurobiol., 163(1–3), 214–221, https://doi.org/10.1016/j.resp.2008.05.013.
  * Hu & Cubaud [2020] Hu, X. & Cubaud, T. 2020 From droplets to waves: Periodic instability patterns in highly viscous microfluidic flows. J. Fluid Mech. 887, A27, https://doi.org/10.1017/jfm.2019.1009.
  * Joseph et al. [1997] Joseph, D. D., Bai, R., Chen, K. & Renardy, Y. Y. 1997 Core–annular flows. Ann. Rev. Fluid Mech. 29, 65–90, https://doi.org/10.1017/jfm.2019.1009.
  * Joseph & Renardy [1993] Joseph, D. & Renardy, Y. 1993 Fundamentals of Two-Fluid Dynamics, Part II. Springer.
  * Levy et al. [2014] Levy, R., Hill, D. B, Forest, M. G. & Grotberg, J. B. 2014 Pulmonary Fluid Flow Challenges for Experimental and Mathematical Modeling. Integr. Comp. Biol., 54(6), 985-–1000, https://doi.org/10.1093/icb/icu107.
  * Liu & Ding [2017] Liu, R. & Ding, Z. 2017 Stability of viscous film flow coating the interior of a vertical tube with a porous wall. Phys. Rev. E 95(5), 053101, https://link.aps.org/doi/10.1103/PhysRevE.95.053101.
  * Ogrosky [2021] Ogrosky, H. R. 2021 Impact of viscosity ratio on falling two-layer viscous film flow inside a tube. Phys. Rev. Fluids, 6, 104005, https://doi.org/10.1103/PhysRevFluids.6.104005.
  * Ogrosky [2021] Ogrosky, H. R. 2021 Linear stability and nonlinear dynamics in a long-wave model of film flows inside a tube in the presence of surfactant. J. Fluid Mech., 908, A23, https://doi.org/10.1017/jfm.2020.878.
  * O’Neill & Mudawar [2020] O’Neill, L. & Mudawar, I. 2020 Review of two-phase flow instabilities in macro- and micro-channel systems. Intl J. Heat Mass transfer, 157, 119738, https://doi.org/10.1016/j.ijheatmasstransfer.2020.119738.
  * Oron et al. [1997] Oron, A., Davis, S. H. & Bankoff, S. G. 1997 Long-scale evolution of thin liquid films. Rev. Mod. Phys. 69(931), 931–980, https://doi.org/10.1103/RevModPhys.69.931.
  * Plateau [1873] Plateau, J. A. F. 1873 Statique expérimentale et théorique des liquides soumis aux seules forces moléculaires. Gauthier-Villars.
  * Revankar & Pollock [2005] Revankar, S. T. & Pollock, D.  2005 Laminar film condensation in a vertical tube in the presence of noncondensable gas. Appl. Math. Model., 29(4), 341–359, https://doi.org/10.1016/j.apm.2004.09.010.
  * Reynolds [1886] Reynolds, O. 1886 IV. On the theory of lubrication and its application to Mr. Beauchamp tower’s experiments, including an experimental determination of the viscosity of olive oil. Phil Trans. R. Soc. Lond., 177, 157-234.
  * Rayleigh [1878] Rayleigh, L. 1878 On the instability of jets. Proc. Lond. Math. Soc., s1-10, 4-13, https://doi.org/10.1112/plms/s1-10.1.4.
  * Rayleigh [1882] Rayleigh, L. 1878 Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density. Proc. Lond. Math. Soc., s1-14, 170–177, https://doi.org/10.1112/plms/s1-14.1.170.
  * Tan & Homsy [1986] Tan, C. T. & Homsy, G. 1986 Stability of miscible displacements in porous media: Rectilinear flow. Phys. Fluids, 29(11), 3549–3556.
  * Taylor [1950] Taylor, G. I. 1950 The instability of liquid surfaces when accelerated in a direction perpendicular to their planes. I. Proc. R. Soc. Lond., A201, 192-–196, https://doi.org/10.1098/rspa.1950.0052.
  * Teng et al. [1999] Teng, H., Cheng, P. & Zhao, T. S. 1999 Instability of condensate film and capillary blocking in small-diameter-thermosyphon condensers. Heat Mass Transfer, 42, 3071–3083, https://doi.org/10.1016/S0017-9310(98)90375-1.
  * Trinh et al. [2014] Trinh, P. H., Kim, H., Hammoud, N., Howell, P. D., Chapman, S. J. & Stone, H. A. 2014 Curvature suppresses the Rayleigh-Taylor instability. Phys. Fluids 26(5), 051704, https://doi.org/10.1063/1.4876476.