# Taylor bubble motion in stagnant and flowing liquids in vertical pipes. Part
II: Linear stability analysis

H. A. Abubakar1,2  O. K. Matar1<EMAIL_ADDRESS>1Department of
Chemical Engineering, Imperial College London, London SW7 2AZ, UK 2Department
of Chemical Engineering, Ahmadu Bello University, Zaria 810107, Nigeria

###### Abstract

In this study, we examine the linear stability of an axisymmetric Taylor
bubble moving steadily in a flowing liquid enclosed in a circular tube.
Linearisation is performed about axisymmetric base states obtained in Part I
of this study by Abubakar & Matar (2021). The stability is characterised by
the dimensionless inverse viscosity $\left(Nf\right)$, Eötvös
$\left(Eo\right)$, and Froude numbers $\left(U_{m}\right)$, the latter being
based on the centreline liquid velocity. The analysis shows that there exist
regions of $(Nf,Eo,U_{m})$ space within which the bubble is unstable and
assumes an asymmetric shape. To elucidate the mechanisms underlying the
instability, an energy budget analysis is carried out which reveals that
perturbation growth is driven by the bubble pressure for $Eo\geq 100$, and by
the tangential interfacial stress for $Eo<100$. Examples of the asymmetric
bubble shapes and their associated flow fields are also provided near the
onset of instability for a wide range of $Nf$, $Eo$, and $U_{m}$.

##  1 Introduction

Slug formation in gas-liquid flows is characterised by intermittent flow of
large gas bubbles, separated by liquid masses. In vertical circular tubes,
these large bubbles are bullet-shaped with diameter approximately equal to
that of the tube, and are known as Taylor bubbles. Because of the pseudo-
periodic nature of the motion of these bubbles, the study of the behaviour of
a single Taylor bubble is considered as a paradigm for understanding the slug
flow regime in vertical tubes. For this reason, extensive experimental
(Griffith & Wallis, 1961; Moissis & Griffith, 1962; White & Beardmore, 1962;
Nicklin et al., 1962; Campos & Guedes de Carvalho, 1988; Bugg & Saad, 2002;
Nogueira et al., 2006; Llewellin et al., 2012; Rana et al., 2015; Pringle et
al., 2015; Fershtman et al., 2017), and a number of theoretical studies
(Dumitrescu, 1943; Brown, 1965; Collins et al., 1978; Funada et al., 2005;
Fabre, 2016) have been carried out to examine the steady bubble shape, rise
speed, and velocity field surrounding the bubble. In addition, numerical
simulations have also been performed a large proportion of which are for
axisymmetric bubbles (Mao & Dukler, 1990, 1991; Bugg & Saad, 2002; Taha & Cui,
2006; Lu & Prosperetti, 2009; Kang et al., 2010), with comparatively fewer
studies focusing on the fully three-dimensional case (Lizarraga-Garcia et al.,
2017; Anjos et al., 2014; Taha & Cui, 2006).

The shape of the Taylor bubble ‘nose’ plays an important role in determining
the rise speed, an important feature for developing predictive models for the
slug flow regime (Mao & Dukler, 1990). A Taylor bubble rising in a stagnant,
or upward-flowing liquid, is generally axisymmetric and moves at constant
speed in an inertia-dominated regime. This is not the case for Taylor bubbles
moving in downward-flowing liquids, however, as experiments have confirmed the
existence of a critical liquid velocity beyond which the bubble shape loses
axisymmetry (Fershtman et al., 2017; Fabre & Figueroa-Espinoza, 2014; Polonsky
et al., 1999; Martin, 1976; Nicklin et al., 1962; Griffith & Wallis, 1961). An
example of asymmetric bubble shapes in downward-flowing liquids is shown in
Figure 1 in which it is seen that the bubble nose becomes distorted, and in an
attempt to avoid the fast-moving fluid at the centre of the tube, the bubble
moves closer to the tube wall (Nicklin et al., 1962), rising faster than it
would have done had it remained axisymmetric (Polonsky et al., 1999; Martin,
1976). It was also noted by Martin (1976) that for a downward-flowing liquid,
a stable axisymmetric Taylor bubble can only be observed in tubes where
surface tension effects are dominant, which is typical of small diameter tubes
characterised by low Eötvös numbers; furthermore, the absolute value of the
downward liquid velocity at which a bubble loses its axisymmetry decreases
with increasing tube diameter.

Motivated by the aforementioned observations, Lu & Prosperetti (2006) carried
out a linear stability analysis of a Taylor bubble moving in a downward-
flowing liquid using potential flow theory and with negligible surface
tension. They demonstrated that an axisymmetric Taylor bubble rising in a
liquid with a laminar velocity profile subjected to irrotational perturbations
becomes unstable beyond a critical negative liquid velocity. Following on from
the theoretical investigation of Lu & Prosperetti (2006) and the numerical
simulation of Figueroa-Espinoza & Fabre (2011), Fabre & Figueroa-Espinoza
(2014) performed experimental investigations in tubes of diameters 20, 40, and
80 mm, complemented by numerical simulations using the boundary element method
of Ha Ngoc & Fabre (2006). They showed that the radius of curvature of the
bubble nose plays a key role in the stability of the Taylor bubble and that
the onset of instability is dependent on the Eötvös number.

(a)

(b)

Figure 1: Taylor bubble shapes in a stagnant, (a), and downward-flowing
liquid, (b), reproduced from Fabre & Figueroa-Espinoza (2014)

In this paper, we examine the stability of an axisymmetric steadily moving
Taylor bubble in stagnant and flowing liquids, with particular attention given
to the transition of the bubble shape from symmetric to asymmetric. To the
best of our knowledge, the study of Lu & Prosperetti (2006) represents the
first attempt in the literature to understand the mechanism governing this
transition using linear theory. However, experimental studies (Fabre &
Figueroa-Espinoza, 2014) have shown that the onset of the instability is
dependent on surface tension and by neglecting it, Lu & Prosperetti (2006)
overestimate this onset. Moreover, existing studies of this transition have
focused on the inertia-dominated regime, necessitating the need for
investigating the parameter spaces where both viscous and surface tension
effects are important. Thus, we carry out a linear stability analysis to
understand how the forces acting on the bubble, characterised by the
dimensionless inverse viscosity, Eötvös, and Froude numbers (the latter being
based on the centreline liquid velocity) affect the loss of bubble
axisymmetry. In addition, an energy budget analysis is carried out to
determine the dominant, perturbation energy-producing terms that drive the
instability. In a companion paper to the present work by Abubakar & Matar
(2021), we computed the steady state solutions for an axisymmetric Taylor
bubble rising steadily in stagnant and flowing liquids for different values of
the aforementioned dimensionless parameters. These solutions serve as the base
states for the linear stability analysis carried out herein. We relate the
conclusions drawn from the steady state analysis by Abubakar & Matar (2021) to
the mechanisms governing the instability.

The rest of this paper is organised as follows. In Section 2, we provide
details of the governing equations, weak formulation, normal modes analysis,
and validation of the numerical technique used to carry out the computations.
A discussion of the results of the linear stability, and of the energy budget
analyses is provided in Sections 3 and 4, respectively. Finally, in Section 5,
concluding remarks and perspectives for future research on remaining open
questions related to Taylor bubble motion are provided.

##  2 Problem formulation

###  2.1 Governing equations and boundary conditions

Figure 2:  Three-dimensional axisymmetric Taylor bubble moving with a steady
speed $U_{b}$ through a liquid which is either stagnant or flowing (upwards or
downwards) in a vertically-aligned tube of diameter $D$.

We consider a situation in which the two-dimensional axisymmetric steady state
solution is revolved in the azimuthal direction to form a three-dimensional
Taylor bubble, as shown in Figure 2. Here, we have adopted a cylindrical polar
coordinate system $(r,\theta,z)$ in which $r$, $\theta$, and $z$ denote the
radial, azimuthal, and axial coordinates, respectively. In the steady state
analysis of Abubakar & Matar (2021), the density, $\rho_{g}$, and viscosity,
$\mu_{g}$, of the gas phase are assumed to be very small compared to their
liquid counterparts, $\rho$ and $\mu$, respectively. Hence, the dynamics in
the gas phase is approximated by a constant pressure ${P_{b}}$, its influence
being restricted to the interface separating the phases designated by
$\Gamma_{b}$, and, consequently, only the liquid flow field and the $P_{b}$
need to be determined.

The flow is governed by the Navier-Stokes and continuity equations which are
non-dimensionalised by the characteristic length, velocity, and pressure
scales, $D,\sqrt{gD},\;\mbox{and}\;\rho gD,$ respectively, where $D$ is the
tube diameter, and $g$ is the gravitational accleration. These equations are
cast in a frame-of-reference translating with the velocity of the steadily-
rising bubble nose $\mathbf{u}_{b}=-U_{b}\mathbf{i}_{z}$, in which $U_{b}$ is
the constant rise speed, and given by

$\frac{\partial\mathbf{u}}{\partial
t}+\left(\mathbf{u}\cdot\nabla\right)\mathbf{u}-\nabla\cdot\mathbf{T}=\mathbf{0},\qquad\text{in}\qquad\Omega\left(t\right)$
(1) $\nabla\cdot\mathbf{u}=0,\qquad\text{in}\qquad\Omega\left(t\right)$ (2)

where $\Omega$ denotes the domain of interest,
$\mathbf{u}=(u_{r},u_{\theta},{u_{z}})$ where $u_{r}$, $u_{\theta}$, and
$u_{z}$ are the components of the liquid velocity vector in the moving frame-
of-reference $\mathbf{u}$ in the $r$, $\theta$, and $z$ directions,
respectively, and $t$ denotes time; $\mathbf{T}$ is the total stress tensor
given by

$\mathbf{T}=-{p}\mathbf{I}+2{Nf}^{-1}\mathbf{E}(\mathbf{u}),$ (3)

in which $p$ represents the dynamic pressure,
$\mathbf{E}=(\nabla\mathbf{u}+\nabla\mathbf{u}^{T})/2$ is the rate of
deformation tensor, and $\nabla=\mathbf{i}_{r}\frac{\partial}{\partial
r}+\mathbf{i}_{\theta}\frac{1}{r}\frac{\partial}{\partial\theta}+\mathbf{i}_{z}\frac{\partial}{\partial
z}$ is the gradient operator in cylindrical polar coordinates where
$\mathbf{i}_{r}$, $\mathbf{i}_{\theta}$, and ${\mathbf{i}_{z}}$ are the unit
vectors in the radial, azimuthal, and axial directions, respectively, and
$\mathbf{I}$ is the unit tensor. In equation $\mathbf{T}$, the inverse
viscosity number $Nf$ is defined as follows

$Nf=\frac{\rho\left(gD^{3}\right)^{\frac{1}{2}}}{\mu}.$ (4)

The boundary of the domain, $\Gamma$, is divided into
$\Gamma_{in},\;\Gamma_{out},\;\Gamma_{wall}$, and $\Gamma_{b}$ as shown in
Figure 2, with the subscripts ‘in’, ‘out’, ‘wall’, and ‘b’ representing the
inlet, outlet, wall, and bubble boundaries, respectively. At the wall, no-slip
and no-penetration boundary conditions are imposed:

$\mathbf{u}=-\mathbf{u}_{b}\quad\mathrm{on\quad\Gamma_{wall}}.$ (5)

At the inlet, prescribed values, $\mathbf{u}_{in}$, are specified for the
velocity components:

$\mathbf{u}=\mathbf{u}_{in}-\mathbf{u}_{b}\quad\mathrm{on\quad\Gamma_{in}}.$
(6)

At the domain outlet, we impose the following conditions:

$\displaystyle\mathrm{\mathbf{n}\cdot\mathbf{T}\cdot\mathbf{n}}=0,$ (7a)
$\displaystyle\mathrm{\left(\mathbf{I}-\mathbf{n}\otimes\mathbf{n}\right)}\cdot\mathbf{u}=\mathbf{0},$
(7b)

where $\mathbf{n}$ is the unit normal vector to the boundary. Finally, at the
interface $\Gamma_{b}$, we set

$\displaystyle\mathbf{n}\cdot\mathbf{T}\cdot{\mathbf{n}}+P_{b}-z-{Eo}^{-1}\kappa=0,$
(8) $\displaystyle\mathbf{n}\cdot\mathbf{T}\times\mathbf{n}=\mathbf{0},$ (9)
$\displaystyle\frac{d\mathbf{r}_{b}}{dt}\cdot\mathbf{n}-\mathbf{u}\cdot\mathbf{n}=0,$
(10)

where $\kappa$ is the curvature of the interface, $\mathbf{r}_{b}(t)$
represents the position vector of all the points on the portion of the
boundary that corresponds to the interface $\Gamma_{b}$, and ${Eo}$ is the
dimensionless Eötvös number expressed by

$Eo=\frac{\rho gD^{2}}{\gamma},$ (11)

in which $\gamma$ denotes the (constant) surface tension. Equations
$\left(\ref{eq:normal_stress_bc}\right)$-$\left(\ref{eq:kinematic_bc}\right)$
correspond to the normal stress, tangential stress, and kinematic boundary
conditions, respectively. Note that gravity appears in
$\left(\ref{eq:normal_stress_bc}\right)$ as $z$ because the hydrostatic
component of the pressure has been subtracted from the total pressure, leaving
only the hydrodynamic part.

###  2.2 Weak forms and perturbation equations

We begin the model development from the weak forms of the momentum and
continuity equations (see appendix A):

$\displaystyle\int_{V}\left\\{\frac{\partial{\mathbf{u}}}{\partial
t}\cdot\mathbf{\Phi}+\left[\left({\mathbf{u}}\cdot\nabla\right)\mathbf{u}\right]\cdot\mathbf{\Phi}+2{Nf}^{-1}\mathbf{E}\left({\mathbf{u}}\right):\mathbf{E}\left(\mathbf{\Phi}\right)-{p}\left(\nabla\ldotp\mathbf{\Phi}\right)\right\\}dV,$
$\displaystyle\quad{}-\int_{A_{b}}\left\\{\left[{Eo}^{-1}\kappa+\mathrm{z}-\mathrm{P}_{\mathrm{b}}\right]\mathbf{n}\cdot\mathbf{\Phi}\right\\}dA_{b}=0,$
(12)
$\int_{V}\left\\{\left(\nabla\ldotp\mathbf{u}\right)\varphi\right\\}dV=0.$
(13)

Also, we have the weak form of the kinematic boundary condition:

$\int_{A_{b}}\left\\{\left[{\frac{d\mathbf{r}_{b}}{dt}\cdot\mathbf{n}-\mathbf{u}\cdot\mathbf{n}}\right]\xi\right\\}dA_{b}=0.$
(14)

In equations (12)-(14), $\mathbf{\Phi}$, $\varphi$, and $\xi$ denote the test
functions for the velocity vector, pressure, and interface deformation
magnitude, respectively.

Similar to domain perturbation, we assume that the three-dimensional base flow
domain is perturbed by the addition of infinitesimal deformation field
$\tilde{\mathbf{x}}$ to its position vector and that equations (12)-(14) are
valid on the three-dimensional perturbed domain. A similar approach has been
adopted in the three-dimensional linear stability analysis of coating flow
problems by Carvalho & Scriven (1999); Christodoulou & Scriven (1988). The
perturbed three-dimensional domain can then be linearised around the base
state three-dimensional axisymmetric domain, the deformation field restricted
to the interface, just as it is expected in the classical linear stability
approach; finally, the linearisation of the perturbed flow field variables
about the base state solution is carried out. This approach to the derivation
of the linear stability model can be seen as an extension of the total
linearisation method used in solving two-dimensional viscous free boundary
problems (Kruyt et al., 1988; Cuvelier & Schulkes, 1990) to three dimensions.

Let the position vector of the perturbed domain be written as

$\mathbf{r}=\mathbf{r}^{0}+\epsilon\tilde{\mathbf{x}},$ (15)

where $\mathbf{r}^{0}=(r^{0},\theta^{0},z^{0})$ represents the position vector
of the unperturbed three-dimensional base flow domain,
$\tilde{\mathbf{x}}=(\tilde{x}_{r},\tilde{x}_{\theta},\tilde{x}_{z})$ is a
deformation field defined over the entire base flow domain, and $\epsilon\ll
1$ to signify the infinitesimally small nature of the applied perturbations.
The linearised elemental volume of the perturbed three-dimensional domain is
given as (Cairncross et al., 2000; Carvalho & Scriven, 1999)

$\displaystyle dV=rdrd\theta dz$ $\displaystyle=$
$\displaystyle(1+\nabla\cdot\tilde{\mathbf{x}})r^{0}dr^{0}d\theta^{0}dz^{0}$
(16) $\displaystyle=$
$\displaystyle(1+\nabla\cdot\tilde{\mathbf{x}})d\Omega^{0}d\theta^{0};$

thus, we can relate the elemental volume in the perturbed three-dimensional
domain to the base flow two-dimensional axisymmetric domain,
$d\Omega^{0}=r^{0}dr^{0}dz^{0}$. Similarly, an elemental area on the perturbed
interface in the three-dimensional domain, $dA_{b}$, can be related to base
flow length of arc, $\Gamma^{0}_{b}$ in the two-dimensional axisymmetric
domain:

$dA_{b}=\left(1+\nabla_{s}\cdot\tilde{\mathbf{x}}\right)d\Gamma^{0}_{b}d\theta^{0},$
(17)

where $\nabla_{s}$ is the surface gradient operator; the interface terms can
be linearised as follows

$\displaystyle\mathbf{n}$
$\displaystyle=\mathbf{n}^{0}+\epsilon\tilde{\mathbf{n}},$ (18a)
$\displaystyle\kappa$ $\displaystyle=\kappa^{0}+\epsilon\tilde{\kappa},$ (18b)
$\displaystyle\mathbf{\Phi}$
$\displaystyle=\mathbf{\Phi}+\epsilon\left(\tilde{\mathbf{x}}\cdot\nabla\right){\mathbf{\Phi}},$
(18c) $\displaystyle\mathbf{u}$
$\displaystyle=\mathbf{u}+\epsilon\left(\tilde{\mathbf{x}}\cdot\nabla\right){\mathbf{u}};$
(18d)

here, $\mathbf{n}^{0}$ and $\kappa^{0}$ are the base state interface normal
vector and curvature, and $\tilde{\mathbf{n}}$ and $\tilde{\kappa}$ represent
the normal vector and curvature perturbations, respectively (see Appendix B):

$\tilde{\mathbf{n}}=-\mathbf{t}^{0}\left(\mathbf{n}^{0}\cdot\frac{d\tilde{\mathbf{x}}}{ds^{0}}\right)-\frac{\mathbf{n}^{0}}{r^{0}}\cdot\frac{\partial\tilde{\mathbf{x}}}{\partial\theta^{0}},\mathbf{i}_{\theta}$
(19)
$\tilde{\kappa}=\frac{1}{r^{0}}\frac{d}{ds^{0}}\left[r^{0}\left(\mathbf{n}^{0}\cdot\frac{d\tilde{\mathbf{x}}}{ds^{0}}\right)\right]+2\left(\mathbf{t}^{0}\cdot\frac{d\tilde{\mathbf{x}}}{ds^{0}}\right)\left(\mathbf{t}^{0}\cdot\frac{d\mathbf{n}^{0}}{ds^{0}}\right)\\\
+\frac{\mathbf{n}^{0}}{{r^{0}}^{2}}\cdot\frac{\partial^{2}\tilde{\mathbf{x}}}{\partial{\theta^{0}}^{2}}+\frac{\tilde{x}_{r}n_{r}^{0}}{{r^{0}}^{2}}-\frac{d\mathbf{n}^{0}}{ds^{0}}\cdot\frac{d\tilde{\mathbf{x}}}{ds^{0}},$
(20)

where $\frac{d}{ds^{0}}=\mathbf{t}\cdot\nabla$ is the derivative along the arc
length $s$ on the base state interface. Substitution into equations (12)-(14)
of equations (15)-(18), together with the flow field perturbations

$\displaystyle\mathbf{u}=\mathbf{u}^{0}+\epsilon\tilde{\mathbf{u}},$ (21a)
$\displaystyle{p}={p}^{0}+\epsilon\tilde{{p}},$ (21b)

followed by neglecting all terms of order $\epsilon^{2}$ respectively yields
the following leading order momentum, continuity, and kinematic condition
equations

$\displaystyle\int_{0}^{2\pi}\left\\{\int_{\Omega^{0}}\left\\{\left[\left({\mathbf{u}^{0}}\cdot\nabla\right)\mathbf{u}^{0}\right]\cdot\mathbf{\Phi}+2{Nf}^{-1}\mathbf{E}\left({\mathbf{u}^{0}}\right):\mathbf{E}\left(\mathbf{\Phi}\right)-{p}^{0}\left(\nabla\ldotp\mathbf{\Phi}\right)\right\\}d\Omega^{0}\right.$
$\displaystyle\left.\quad{}-\int_{\Gamma_{b}^{0}}\left\\{\left[{Eo}^{-1}\kappa^{0}+z^{0}-P_{b}^{0}\right]\mathbf{n}^{0}\cdot\mathbf{\Phi}\right\\}d\Gamma_{b}^{0}\right\\}d\theta^{0},$
(22)
$\displaystyle\int_{0}^{2\pi}\int_{\Omega^{0}}\left\\{\left\\{\left(\nabla\ldotp\mathbf{u}^{0}\right)\varphi\right\\}\left[1+\nabla\cdot\tilde{\mathbf{x}}\right]\right\\}d\Omega^{0}d\theta^{0},$
(23)
$\displaystyle\int_{0}^{2\pi}\int_{\Gamma_{b}^{0}}\left\\{\left\\{\left[{\frac{d\mathbf{r}_{b}^{0}}{dt}\cdot\mathbf{n}^{0}-\mathbf{u}^{0}\cdot\mathbf{n}^{0}}\right]\xi\right\\}\left[1+\nabla_{s}\cdot\tilde{\mathbf{x}}\right]\right\\}d\Gamma_{b}^{0}d\theta^{0}.$
(24)

It is also possible to write the following equations at $O(\epsilon)$ to yield
equations that feature the perturbation variables: the momentum conservation
equation,

$\displaystyle\int_{0}^{2\pi}\int_{\Omega^{0}}\left\\{\frac{\partial\tilde{\mathbf{u}}}{\partial
t}\cdot\mathbf{\Phi}+{\left[\left({\mathbf{u}^{0}}\cdot\nabla\right)\tilde{\mathbf{u}}+\left(\tilde{\mathbf{u}}\cdot\nabla\right)\mathbf{u}^{0}\right]\cdot\mathbf{\Phi}}+2{Nf}^{-1}\mathbf{E}\left(\tilde{\mathbf{u}}\right):\mathbf{E}\left(\mathbf{\Phi}\right)\right\\}d\Omega^{0}d\theta^{0}$
$\displaystyle\quad{}-\int_{0}^{2\pi}\int_{\Omega^{0}}\left\\{\tilde{p}\left(\nabla\ldotp\mathbf{\Phi}\right)\right\\}d\Omega^{0}d\theta^{0}$
$\displaystyle\quad{}-\int_{0}^{2\pi}\int_{\Gamma_{b}^{0}}\left\\{\left[{Eo}^{-1}\tilde{\kappa}+\tilde{z}\right]\mathbf{n}^{0}\cdot\mathbf{\Phi}\right\\}d\Gamma_{b}^{0}d\theta^{0}$
$\displaystyle\quad{}+{\int_{0}^{2\pi}\int_{\Gamma^{0}_{b}}\tilde{\mathbf{x}}\cdot\mathbf{n}^{0}\left\\{{\left[\left({\mathbf{u}^{0}}\cdot\nabla\right)\mathbf{u}^{0}\right]\cdot\mathbf{\Phi}}+2{Nf}^{-1}\mathbf{E}\left(\mathbf{u}^{0}\right):\mathbf{E}\left(\mathbf{\Phi}\right)-p^{0}\left(\nabla\ldotp\mathbf{\Phi}\right)\right\\}d\Gamma^{0}_{b}d\theta^{0}}$
$\displaystyle\quad{}-\int_{0}^{2\pi}\int_{\Gamma_{b}^{0}}\left\\{\left[{Eo}^{-1}\kappa^{0}+z^{0}-P_{b}^{0}\right]\left[\tilde{\mathbf{n}}\cdot\mathbf{\Phi}+\left[\left(\tilde{\mathbf{x}}\cdot\nabla\right){\mathbf{\Phi}}\right]\cdot\mathbf{n}^{0}+\left(\nabla_{s}\cdot\tilde{\mathbf{x}}\right)\mathbf{n}^{0}\cdot\mathbf{\Phi}\right]\right\\}d\Gamma_{b}^{0}d\theta^{0}$
$\displaystyle\quad{}=0,$ (25)

where $\mathbf{u}^{0}$, $p^{0}$, and $P_{b}^{0}$ represent the base flow
solutions for the variables; and $\tilde{\mathbf{u}}$ and $\tilde{p}$ denote
the perturbations to the flow field variables, and the last two lines of (25)
are due to the linearisation of the domain and boundary terms of equation (12)
where we have used Gauss’s divergence theorem to restrict the deformation to
the interface as it is expected in classical linear stability formulation; the
continuity equation,

$\int_{0}^{2\pi}\int_{\Omega^{0}}\left\\{\left(\nabla\ldotp\tilde{\mathbf{u}}\right)\varphi\right\\}d\Omega^{0}d\theta^{0}=0,$
(26)

and the kinematic condition:

$\int_{0}^{2\pi}\int_{\Gamma_{b}^{0}}\left\\{\left[{\frac{d\tilde{\mathbf{x}}}{dt}\cdot\mathbf{n}^{0}-\tilde{\mathbf{u}}\cdot\mathbf{n}^{0}-\mathbf{u}^{0}\cdot\tilde{\mathbf{n}}-\left[\left(\tilde{\mathbf{x}}\cdot\nabla\right)\mathbf{u}^{0}\right]\cdot\mathbf{n}^{0}}\right]\xi\right\\}d\Gamma_{b}^{0}d\theta^{0}=0.$
(27)

Simplifying equations (25)-(27) further by substituting for
$\tilde{\mathbf{n}}$ and $\tilde{\kappa}$ using equations (19) and (20),
respectively, and taking the deformation field to be of the form
$\tilde{\mathbf{x}}=\tilde{h}\mathbf{n}^{0}$, since the deformation has been
restricted to the interface and making use of the relations

$\displaystyle\mathbf{u}^{0}=\left(\mathbf{u}^{0}\cdot\mathbf{n}^{0}\right)\mathbf{n}^{0}$
$\displaystyle+\left(\mathbf{u}^{0}\cdot\mathbf{t}^{0}\right)\mathbf{t}^{0},$
(28a)
$\displaystyle{\nabla}=\left(\mathbf{I}-\mathbf{n}^{0}\otimes\mathbf{n}^{0}\right)\cdot\mathbf{\nabla}$
$\displaystyle+\left(\mathbf{n}^{0}\otimes\mathbf{n}^{0}\right)\cdot{\nabla},$
(28b)

equations (25)-(27), after some algebra, can be expressed as follows

$\displaystyle\int_{0}^{2\pi}\int_{\Omega^{0}}\left\\{\frac{\partial\tilde{\mathbf{u}}}{\partial
t}\cdot\mathbf{\Phi}+{\left[\left({\mathbf{u}^{0}}\cdot\nabla\right)\tilde{\mathbf{u}}+\left(\tilde{\mathbf{u}}\cdot\nabla\right)\mathbf{u}^{0}\right]\cdot\mathbf{\Phi}}+2{Nf}^{-1}\mathbf{E}\left(\tilde{\mathbf{u}}\right):\mathbf{E}\left(\mathbf{\Phi}\right)\right\\}d\Omega^{0}d\theta$
$\displaystyle\quad{}-\int_{0}^{2\pi}\int_{\Omega^{0}}\left\\{\tilde{p}\left(\nabla\ldotp\mathbf{\Phi}\right)\right\\}d\Omega^{0}d\theta$
$\displaystyle\quad{}-\int_{0}^{2\pi}\int_{\Gamma_{b}^{0}}{Eo}^{-1}\left\\{-\frac{d\tilde{h}}{ds}\left[\mathbf{n}\cdot\frac{d\mathbf{\Phi}}{ds}-\kappa_{a}\left(\mathbf{t}\cdot\mathbf{\Phi}\right)\right]+\left[\tilde{h}\left(\kappa_{a}^{2}+\kappa_{b}^{2}\right)+\frac{1}{r^{2}}\frac{\partial^{2}\tilde{h}}{\partial\theta^{2}}\right]\mathbf{n}\cdot\mathbf{\Phi}\right\\}d\Gamma_{b}^{0}d\theta$
$\displaystyle\quad{}-\int_{0}^{2\pi}\int_{\Gamma_{b}^{0}}\left\\{\tilde{h}n_{z}\left(\mathbf{n}\cdot\mathbf{\Phi}\right)\right\\}d\Gamma_{b}^{0}d\theta$
$\displaystyle\quad{}+\int_{0}^{2\pi}\int_{\Gamma^{0}_{b}}\tilde{h}\left\\{{\left(\mathbf{u}^{0}\cdot\mathbf{t}\right)\left[\frac{d\mathbf{u}^{0}}{ds}\cdot\mathbf{\Phi}\right]}+\left[-p^{0}+2Nf^{-1}\left(\mathbf{t}\cdot\frac{d\mathbf{u}^{0}}{ds}\right)\right]\left(\mathbf{t}\cdot\frac{d\mathbf{\Phi}}{ds}\right)\right.$
$\displaystyle\quad{}\left.+\left[-p^{0}+2Nf^{-1}\frac{u_{r}^{0}}{r}\right]\left(\frac{\Phi_{r}}{r}+\frac{1}{r}\frac{\partial\Phi_{\theta}}{\partial\theta}\right)\right\\}d\Gamma^{0}_{b}d\theta$
$\displaystyle\quad{}+\int_{0}^{2\pi}\int_{\Gamma_{b}^{0}}\left\\{\left[{Eo}^{-1}\kappa+z-P_{b}^{0}\right]\left[\left(\mathbf{t}\cdot\mathbf{\Phi}\right)\frac{d\tilde{h}}{ds}+\frac{\Phi_{\theta}}{r}\frac{\partial\tilde{h}}{\partial\theta}+\tilde{h}\kappa\left(\mathbf{n}\cdot\mathbf{\Phi}\right)\right]\right\\}d\Gamma_{b}^{0}d\theta=0,$
(29)
$\int_{0}^{2\pi}\int_{\Omega^{0}}\left\\{\left(\nabla\ldotp\tilde{\mathbf{u}}\right)\varphi\right\\}d\Omega
d\theta=0,$ (30)
$\int_{0}^{2\pi}\int_{\Gamma_{b}^{0}}\left\\{\left[{\frac{d\tilde{h}}{dt}-\tilde{\mathbf{u}}\cdot\mathbf{n}+\left(\mathbf{u}^{0}\cdot\mathbf{t}\right)\frac{d\tilde{h}}{ds}-\tilde{h}\left(\frac{d\mathbf{u}^{0}}{dn}\cdot\mathbf{n}\right)}\right]\xi\right\\}d\Gamma_{b}^{0}d\theta=0,$
(31)

where $\frac{d}{dn}=\left(\mathbf{n}\cdot\nabla\right)$ is the derivative in
the normal direction. In equations (29)-(31), we have suppressed the use of
the superscript $`0^{\prime}$ to designate base state quantities for the unit
tangent and normal vectors for the sake of brevity.

###  2.3 Normal modes

Let us take the following normal mode forms for the perturbation variables:

$\displaystyle\tilde{\mathbf{u}}\left(r,\theta,z,t\right)$
$\displaystyle=\hat{\mathbf{u}}\left(r,z\right)e^{\left(\mathfrak{i}m\theta+\beta
t\right)},$ (32a) $\displaystyle\tilde{p}\left(r,\theta,z,t\right)$
$\displaystyle=\hat{p}\left(r,z\right)e^{\left(\mathfrak{i}m\theta+\beta
t\right)},$ (32b) $\displaystyle\tilde{h}\left(s,\theta,t\right)$
$\displaystyle=\hat{h}\left(s\right)e^{\left(\mathfrak{i}m\theta+\beta
t\right)},$ (32c)

and their corresponding test functions as

$\displaystyle\mathbf{\Phi}\left(r,\theta,z\right)$
$\displaystyle=\bar{\mathbf{\Phi}}\left(r,z\right)e^{\left(-\mathfrak{i}m\theta\right)},$
(33a) $\displaystyle{\varphi}\left(r,\theta,z\right)$
$\displaystyle=\bar{\varphi}\left(r,z\right)e^{\left(-\mathfrak{i}m\theta\right)},$
(33b) $\displaystyle\xi\left(s,\theta\right)$
$\displaystyle=\bar{\xi}\left(s\right)e^{\left(-\mathfrak{i}m\theta\right)},$
(33c)

where $\hat{\mathbf{u}}$, $\hat{p}$, and $\hat{h}$ are complex functions of
space representing the amplitude of the velocity, pressure, and interface
deformation perturbations, respectively; $m$ is a dimensionless (integer) wave
number in the azimuthal direction $\theta$; $\beta=\beta_{R}+i\beta_{I}$ is
the complex growth rate which can be decomposed into its real $\beta_{R}$ and
imaginary $\beta_{I}$ parts denoting the temporal growth rate and frequency,
respectively: if $\beta_{R}$ is positive (negative), the disturbance grows
(decays) exponentially in time and the base flow is linearly unstable
(stable); if $\beta_{R}$ is zero, the disturbance is neutrally stable.
Substituting equations (32) and (33) into (29)-(31), separating the momentum
equation into its components, yields the following equations governing the
normal mode evolution of the perturbations as a function of $Nf$, $Eo$,
$U_{m}$, and $m$:

$\displaystyle\int_{\Omega^{0}}\left\\{\beta\hat{u}_{r}{\bar{\Phi}_{r}}+\left[\hat{u}_{r}\frac{\partial{u}_{r}^{0}}{\partial
r}+\hat{u}_{z}\frac{\partial{u}_{r}^{0}}{\partial
z}+{u}_{r}^{0}\frac{\partial\hat{u}_{r}}{\partial
r}+{u}_{z}^{0}\frac{\partial\hat{u}_{r}}{\partial
z}\right]{\bar{\Phi}_{r}}+Nf^{-1}\left[2\frac{\partial\hat{u}_{r}}{\partial
r}\frac{\partial\bar{\Phi}_{r}}{\partial r}\right.\right.$
$\displaystyle\left.\left.\quad{}+\left(2+m^{2}\right)\frac{\hat{u}_{r}\bar{\Phi}_{r}}{r^{2}}+3\mathfrak{i}m\frac{\hat{u}_{\theta}\bar{\Phi}_{r}}{r^{2}}-\mathfrak{i}m\frac{\bar{\Phi}_{r}}{r}\frac{\partial\hat{u}_{\theta}}{\partial
r}+\frac{\partial\bar{\Phi}_{r}}{\partial
z}\left(\frac{\partial\hat{u}_{r}}{\partial
z}+\frac{\partial\hat{u}_{z}}{\partial r}\right)\right]\right.$
$\displaystyle\left.\quad{}-\hat{p}\left(\frac{\partial\bar{\Phi}_{r}}{\partial
r}+\frac{{\bar{\Phi}_{r}}}{r}\right)\right\\}d\Omega^{0}$
$\displaystyle\quad{}-\int_{\Gamma_{b}^{0}}{Eo}^{-1}\left\\{-\frac{d\hat{h}}{ds}\left[n_{r}\frac{d\bar{\Phi}_{r}}{ds}-\kappa_{a}\left(t_{r}\bar{\Phi}_{r}\right)\right]+\hat{h}\left[\kappa_{a}^{2}+\kappa_{b}^{2}-\frac{m^{2}}{r^{2}}\right]n_{r}\bar{\Phi}_{r}\right\\}d\Gamma_{b}^{0}$
$\displaystyle\quad{}-\int_{\Gamma_{b}^{0}}\left\\{\hat{h}n_{z}\left(n_{r}\bar{\Phi}_{r}\right)\right\\}d\Gamma_{b}^{0}$
$\displaystyle\quad{}+\int_{\Gamma^{0}_{b}}\hat{h}\left\\{{\left({u}^{0}_{r}{t}_{r}\right)\left[\frac{d{u}^{0}_{r}}{ds}\bar{\Phi}_{r}\right]}+\left[-p^{0}+2Nf^{-1}\left({t}_{r}\frac{d{u}^{0}_{r}}{ds}\right)\right]\left({t}_{r}\frac{d\bar{\Phi}_{r}}{ds}\right)\right.$
$\displaystyle\left.\quad{}+\left[-p^{0}+2Nf^{-1}\frac{u_{r}^{0}}{r}\right]\left(\frac{\bar{\Phi}_{r}}{r}\right)\right\\}d\Gamma^{0}_{b}$
$\displaystyle\quad{}+\int_{\Gamma_{b}^{0}}\left\\{\left[{Eo}^{-1}\kappa+z-P_{b}^{0}\right]\left[\left({t}_{r}\bar{\Phi}_{r}\right)\frac{d\hat{h}}{ds}+\hat{h}\kappa\left({n}_{r}\bar{\Phi}_{r}\right)\right]\right\\}d\Gamma_{b}^{0}=0,$
(34)
$\displaystyle\int_{\Omega^{0}}\left\\{\beta\hat{u}_{\theta}{\bar{\Phi}_{\theta}}+\left[{u}_{r}^{0}\frac{\partial\hat{u}_{\theta}}{\partial
r}+{u}_{z}^{0}\frac{\partial\hat{u}_{\theta}}{\partial
z}+\frac{{u}_{r}^{0}\hat{u}_{\theta}}{r}\right]{\bar{\Phi}_{\theta}}+Nf^{-1}\left[\left(1+2m^{2}\right)\frac{\hat{u}_{\theta}\bar{\Phi}_{\theta}}{r^{2}}+\frac{\partial\hat{u}_{\theta}}{\partial
z}\frac{\partial\bar{\Phi}_{\theta}}{\partial z}\right.\right.$
$\displaystyle\quad{}\left.\left.+\frac{\partial\hat{u}_{\theta}}{\partial
r}\frac{\partial\bar{\Phi}_{\theta}}{\partial
r}-\left(\frac{\hat{u}_{\theta}}{r}\frac{\partial\bar{\Phi}_{\theta}}{\partial
r}+\frac{\bar{\Phi}_{\theta}}{r}\frac{\partial\hat{u}_{\theta}}{\partial
r}\right)+\mathfrak{i}m\left(\frac{\hat{u}_{r}}{r}\frac{\partial\bar{\Phi}_{\theta}}{\partial
r}+\frac{\hat{u}_{z}}{r}\frac{\partial\bar{\Phi}_{\theta}}{\partial
z}\right)-3\mathfrak{i}m\frac{\hat{u}_{r}\bar{\Phi}_{\theta}}{r^{2}}\right]\right.$
$\displaystyle\quad{}\left.-p\left(-\mathfrak{i}m\frac{\bar{\Phi}_{\theta}}{r}\right)\right\\}d\Omega^{0}+\int_{\Gamma^{0}_{b}}\hat{h}\left\\{\left[-p^{0}+2Nf^{-1}\frac{u_{r}^{0}}{r}\right]\left(-\mathfrak{i}m\frac{\bar{\Phi}_{\theta}}{r}\right)\right\\}d\Gamma^{0}_{b}$
$\displaystyle\quad{}+\int_{\Gamma_{b}^{0}}\left\\{\hat{h}\left[{Eo}^{-1}\kappa+z-P_{b}^{0}\right]\left[\mathfrak{i}m\frac{\bar{\Phi}_{\theta}}{r}\right]\right\\}d\Gamma_{b}^{0}=0,$
(35)
$\displaystyle\int_{\Omega^{0}}\left\\{\beta\hat{u}_{z}{\bar{\Phi}_{z}}+\left[\hat{u}_{r}\frac{\partial{u}_{z}^{0}}{\partial
r}+\hat{u}_{z}\frac{\partial{u}_{z}^{0}}{\partial
z}+{u}_{r}^{0}\frac{\partial\hat{u}_{z}}{\partial
r}+{u}_{z}^{0}\frac{\partial\hat{u}_{z}}{\partial
z}\right]{\bar{\Phi}_{z}}+Nf^{-1}\left[2\frac{\partial\hat{u}_{z}}{\partial
z}\frac{\partial\bar{\Phi}_{z}}{\partial
z}+m^{2}\frac{\hat{u}_{z}\bar{\Phi}_{z}}{r^{2}}\right.\right.$
$\displaystyle\quad{}\left.\left.-\mathfrak{i}m\frac{\bar{\Phi}_{z}}{r}\frac{\partial\hat{u}_{\theta}}{\partial
z}+\frac{\partial\bar{\Phi}_{z}}{\partial
r}\left(\frac{\partial\hat{u}_{r}}{\partial
z}+\frac{\partial\hat{u}_{z}}{\partial
r}\right)\right]-\hat{p}\left(\frac{\partial\bar{\Phi}_{z}}{\partial
z}\right)\right\\}d\Omega^{0}$
$\displaystyle\quad{}-\int_{\Gamma_{b}^{0}}{Eo}^{-1}\left\\{-\frac{d\hat{h}}{ds}\left[n_{z}\frac{d\bar{\Phi}_{z}}{ds}-\kappa_{a}\left(t_{z}\bar{\Phi}_{z}\right)\right]+\hat{h}\left[\kappa_{a}^{2}+\kappa_{b}^{2}-\frac{m^{2}}{r^{2}}\right]n_{z}\bar{\Phi}_{z}\right\\}d\Gamma_{b}^{0}$
$\displaystyle\quad{}-\int_{\Gamma_{b}^{0}}\left\\{\hat{h}n_{z}\left(n_{z}\bar{\Phi}_{z}\right)\right\\}d\Gamma_{b}^{0}$
$\displaystyle\quad{}+\int_{\Gamma^{0}_{b}}\hat{h}\left\\{{\left({u}^{0}_{z}{t}_{z}\right)\left[\frac{d{u}^{0}_{z}}{ds}\bar{\Phi}_{z}\right]}+\left[-p^{0}+2Nf^{-1}\left({t}_{z}\frac{d{u}^{0}_{z}}{ds}\right)\right]\left({t}_{z}\frac{d\bar{\Phi}_{z}}{ds}\right)\right\\}d\Gamma^{0}_{b}$
$\displaystyle\quad{}+\int_{\Gamma_{b}^{0}}\left\\{\left[{Eo}^{-1}\kappa+z-P_{b}^{0}\right]\left[\left({t}_{z}\bar{\Phi}_{z}\right)\frac{d\hat{h}}{ds}+\hat{h}\kappa\left({n}_{z}\bar{\Phi}_{z}\right)\right]\right\\}d\Gamma_{b}^{0}=0,$
(36) $\int_{\Omega^{0}}\left\\{\left[\frac{\partial\hat{u}_{r}}{\partial
r}+\frac{\hat{u}_{r}}{r}-\mathfrak{i}m\frac{\hat{u}_{\theta}}{r}+\frac{\partial\hat{u}_{z}}{\partial
z}\right]\bar{\varphi}\right\\}d\Omega^{0}=0,$ (37)
$\int_{\Gamma_{b}^{0}}\left\\{\left[{\beta\hat{h}-\hat{\mathbf{u}}\cdot\mathbf{n}+\left(\mathbf{u}^{0}\cdot\mathbf{t}\right)\frac{d\hat{h}}{ds}-\hat{h}\left(\frac{d\mathbf{u}^{0}}{dn}\cdot\mathbf{n}\right)}\right]\bar{\xi}\right\\}d\Gamma_{b}^{0}=0.$
(38)

The combined finite element forms for the perturbations equations (34)-(38)
can be recast as a generalised eigenvalue problem

$\beta\mathbf{B}\mathbf{y}=\mathbf{J}\mathbf{y},$ (39)

with $\beta$ being the eigenvalue, $\mathbf{B}$ the mass matrix, $\mathbf{y}$
the eigenfunctions, and $\mathbf{J}$ the Jacobian matrix respectively given by

$\displaystyle{\mathbf{B}=\begin{bmatrix}{\mathbf{B}}_{r,r}&0&0&0&0\\\
0&{\mathbf{B}}_{\theta,\theta}&0&0&0\\\ 0&0&{\mathbf{B}}_{z,z}&0&0\\\
0&0&0&0&0\\\
0&0&0&0&{\mathbf{B}}_{h}\end{bmatrix},}\quad{\mathbf{y}=\begin{bmatrix}\hat{\mathbf{v}}_{r}\\\
\hat{\mathbf{v}}_{\theta}\\\ \hat{\mathbf{v}}_{z}\\\ \hat{\mathbf{p}}\\\
\hat{\mathbf{h}}\end{bmatrix},}\quad\quad$
$\displaystyle{\mathbf{J}=-\begin{bmatrix}\left({\mathbf{C}}_{r,r}+{\mathbf{K}}_{r,r}\right)&{\mathbf{K}}_{r,\theta}&\left({\mathbf{C}}_{r,z}+{\mathbf{K}}_{r,z}\right)&-{\mathbf{Q}}_{r}^{T}&{\mathbf{M}}_{r,h}\\\
{\mathbf{K}}_{r,\theta}^{T}&\left({\mathbf{C}}_{\theta,\theta}+{\mathbf{K}}_{\theta,\theta}\right)&{\mathbf{K}}_{z,\theta}^{T}&-{\mathbf{Q}}_{\theta}^{T}&{\mathbf{M}}_{\theta,h}\\\
\left({\mathbf{C}}_{z,r}+{\mathbf{K}}_{r,z}^{T}\right)&{\mathbf{K}}_{z,\theta}&\left({\mathbf{C}}_{z,z}+{\mathbf{K}}_{z,z}\right)&-{\mathbf{Q}}_{z}^{T}&{\mathbf{M}}_{z,h}\\\
-{\mathbf{Q}}_{r}&-{\mathbf{Q}}_{\theta}&-{\mathbf{Q}}_{z}&0&0\\\
-{\mathbf{X}}_{r}&0&-{\mathbf{X}}_{z}&0&{\mathbf{X}}_{h}\end{bmatrix}}$

where $\mathbf{B}_{i,i}\quad\forall\;i=r,\theta,z$ are coefficient matrices
for the $i$ component of velocity in the $i$ component of momentum equation.
$\mathbf{C}_{i,j}\quad\text{and}\quad\mathbf{K}_{i,j}\quad\forall\;i=r,\theta,z$
are the convective and viscous coefficient matrices for the $j$ component of
velocity in the $i$ component of momentum equation, respectively.
$\mathbf{M}_{i,h}\quad\forall\;i=r,\theta,z$ are coefficient matrices for
interface deformation magnitude in the $i$ component of momentum equation.
$\mathbf{Q}_{j}\quad\forall\;j=r,\theta,z$ are coefficient matrices for the
$j$ component of velocity in the continuity equation, whose Hermitian
transpose correspond to the coefficient matrices of pressure in the $j$
component of momentum equation. $\mathbf{X}_{j}\quad\forall\;j=r,z$ are
coefficient matrices for the $j$ component of velocity in the kinematic
boundary condition. $\mathbf{B}_{h}\quad\text{and}\quad\mathbf{X}_{h}$ are the
coefficient matrices of interface deformation magnitude in the kinematic
boundary condition for the mass and Jacobian matrices, respectively. The
operator $\left(\cdot\right)^{T}$ denotes the Hermitian transpose. The full
expressions for the coefficient matrices in the mass and Jacobian matrices are
given in Appendix C.

The boundary conditions at the inlet, wall, and outlet reduce to the following
conditions on the perturbations:

$\hat{\mathbf{u}}=\mathbf{0}\quad\text{on}\quad\Gamma_{in}^{0}\;\text{and}\;\Gamma_{wall}^{0},$
(40)
$\mathbf{n}\cdot\hat{\mathbf{T}}\cdot\mathbf{n}=0\quad\text{and}\quad\left(\mathbf{I}-\mathbf{n}\otimes\mathbf{n}\right)\cdot\hat{\mathbf{u}}=0\quad\text{on}\quad\Gamma_{out}^{0},$
(41)

where the tensor $\mathbf{\hat{T}}$ is expressed by

$\displaystyle\hat{\mathbf{T}}=-\hat{p}+2Nf^{-1}\hat{\mathbf{E}}\left(\hat{\mathbf{u}}\right);\quad\quad\hat{\mathbf{E}}\left(\hat{\mathbf{u}}\right)=\frac{1}{2}\left[\nabla\hat{\mathbf{u}}+{\nabla\hat{\mathbf{u}}}^{T}\right].$

We stress that while it is customary to impose additional conditions along the
axis of symmetry $\Gamma_{sym}^{0}$, we did not apply any such conditions in
this case because the model equations were written around the perturbed three-
dimensional domain and then linearised before integrating out the $\theta$
dependence.

###  2.4 Numerical method and validation

####  2.4.1 Linear algebra

The linear stability of the steady state solutions as a parametric function of
the system dimensionless groups is determined by solving a generalized,
asymmetric matrix eigenvalue problem given by equation (39). The asymmetric
nature of the problem, as can be seen by inspection of the Jacobian matrix, is
due to the convective and interface deformation contributions to the matrix.
In addition, it can also be seen that the mass matrix is singular, with
several rows having identically zero entries, which can be attributed to the
absence of time derivative terms in the continuity equation and though not
obvious, due to the requirement that the perturbation vectors satisfy
homogeneous essential boundary conditions (Carvalho & Scriven, 1999;
Natarajan, 1992). The implication of having singularity in the mass matrix is
that the number of true eigenvalues is smaller than the dimension of the
problem, with the number of missing eigenvalues equal to the number of
algebraic equations having identically zero row in the matrix (Carvalho &
Scriven, 1999). We carry out a shift-and-invert transformation (Christodoulou
& Scriven, 1988) to map these missing, or so-called ‘infinite’, eigenvalues to
zero wherein the original problem is transformed to

$\displaystyle\left(\mathbf{J}-\upsilon\mathbf{B}\right)^{-1}\mathbf{B}\mathbf{y}=\tau\mathbf{y},$
(42)

where the shift $\upsilon$ is a complex constant, $\tau$ is the eigenvalue of
the new problem and is related to the eigenvalue of the original problem
$\beta$ by

$\tau=\frac{1}{\left(\beta-\upsilon\right)}.$ (43)

We solve this eigenvalue problem by using an iterative Arnoldi method
available in ARPACK (Lehoucq et al., 1997), which can be called within
FreeFem++, using the standard Taylor-Hood element for the flow field variables
and piecewise quadratic element for interface deformation magnitude as in the
steady state simulations. The accuracy of the converged leading eigenvalue is
confirmed by ensuring that the residual
$\left|\mathbf{J}\mathbf{y}-\beta\mathbf{B}\mathbf{y}\right|$ is always less
than $1\times 10^{-10}$.

####  2.4.2 Validation

We test the validity of our theoretical and numerical procedure by examining
the stability of a spherical bubble of fixed volume in a stagnant liquid with
negligible gravitational and boundary effects. The bubble is stable under
these conditions and its motion is governed by an analytical solution
(Prosperetti, 1980; Miller & Scriven, 1968). We compare our numerical results
for the eigenvalues with this solution given in Prosperetti (1980) for small
amplitude normal mode perturbations. The characteristic scales used for the
non-dimensionalisation of space, velocity, and pressure in the governing
equations are $R$, $\sqrt{\gamma/(\rho R)}$, and $\gamma/R$, respectively,
where $R$ is the bubble radius, so that the validation problem is
parameterised by the Ohnesorge number, $Oh\mu/\sqrt{\rho R\gamma}$.

Based on the scaling above, the dimensionless form of the characteristic
equation for the bubble oscillations reads (Prosperetti, 1980)

$\displaystyle\left[\mathcal{H}^{\left(1\right)}_{m-\frac{1}{2}}\left(X^{*}\right)\right]\beta^{2}+Oh\left[4m(m+2)^{2}-2(m+2)(2m+1)(\mathcal{H}^{(1)}_{m-\frac{1}{2}}\left(X^{*}\right)+2)\right]\beta$
$\displaystyle+(m+1)(m-1)(m+2)(\mathcal{H}^{(1)}_{m-\frac{1}{2}}\left(X^{*}\right)+2)=0,$
(44)

where $X^{*}$ is a rescaled growth rate, and
$\mathcal{H}^{\left(1\right)}_{j}(X^{*})$ is a Hankel function of the first
kind:

$X^{*}=\left[\frac{\beta}{Oh}\right]^{\frac{1}{2}},~{}~{}{\rm
and}~{}~{}~{}\mathcal{H}^{(1)}_{j}(X^{*})=\frac{X^{*}H^{(1)}_{j+1}(X^{*})}{H^{(1)}_{j}(X^{*})}.$
(45)

For a fixed value of $Oh$, we solve iteratively for $\beta$. The initial guess
used is the solution to the following equation (Prosperetti, 1980)

$\beta^{2}-2Oh\left[(m+2)(2m+1)\right]\beta+(m+1)(m-1)(m+2)=0.$ (46)

Once the solution for the first eigenvalue is obtained, we use the associated
$X^{*}$ for the previous $Oh$ as the initial guess for the next value of $Oh$.
We implemented the solution steps in MATLAB and generated the analytical
solution for $0\leq Oh\leq 1$.

At steady state, in the absence of gravity and since the liquid surrounding
the bubble is stagnant $(\mathbf{u}=0)$, the governing equations reduce to

$\displaystyle\nabla p$ $\displaystyle=0,$ (47) and the normal stress boundary
condition to $\displaystyle-p+P_{b}$
$\displaystyle=\kappa\;\text{on}\;{\Gamma}_{b}.$ (48)

Equation (47) implies that pressure field in the liquid phase surrounding the
bubble is a constant, $P_{a}$, so that the bubble pressure becomes

$P_{b}=\kappa+P_{a}\;~{}~{}~{}\text{on}\;{\Gamma}_{b}.$ (49)

For the linear stability analysis, the value of $P_{a}$ was set to zero
without loss of generality.

We solve the modified forms of the perturbation equations (34)-(38) using the
base state solutions computed as set out above. Figures 3(a) and 3(b)
respectively show excellent agreement between the real and imaginary parts of
the eigenvalues computed and the analytical solution of (44) as a function of
the Ohnesorge number for four different azimuthal wavenumbers.

(a)

(b)

Figure 3: Validation of the theoretical and numerical procedure for an
oscillating bubble in the absence of gravitational and boundary effects.
Comparison between the amplification rates, (a), and oscillation frequencies,
(b), from the analytical solution given by equation (44) (coloured continuous
solid line) and our numerically-generated growth rates (coloured markers), for
modes m = 2, 3, 4 and 5.

##  3 Linear stability results

In this section, we provide a discussion of the linear stability results
starting with the dependence of the growth rate $\beta_{R}$ obtained from the
leading eigenvalues associated primarily with the first two modes, $m=1$ and
$m=2$, as a parametric function of $Nf$, $Eo$, and $U_{m}$. The asymmetric
bubble shape near instability onset, associated with the most dangerous linear
mode, is also discussed, and a stability map is plotted which clearly
demarcates the stability boundary as a function of the system parameters.

###  3.1 Dominant modes of instability

The linear stability of axisymmetric steady states associated with $40\leq
Nf\leq 100$ and $20\leq Eo\leq 300$ was investigated for downward liquid flow
characterised by $U_{m}<0$; these base states had been computed and reported
in the companion paper to the present work by Abubakar & Matar (2021). In
Figure 4, we show the dependence of the growth rate $\beta_{R}$ on $U_{m}$ for
modes $m=1$, $m=2$, and $m=3$ for $Nf=40,60,80,100$, and $Eo=20,180,300$. For
all the cases shown in this figure, it is evident that $m=1$ is the most
unstable mode, which corresponds to a deflection of the bubble away from the
axis of symmetry and occurs over a well-defined range of negative $U_{m}$
values, $\Delta U_{m}$, for which $\beta_{R}>0$ indicating the presence of a
linear instability. For sufficiently large $Eo$ and $Nf$, the $m=2$ is also
unstable, though it remains sub-dominant to the $m=1$ mode, as illustrated in
Figure 4(l), for instance, for the $Nf=100$ and $Eo=300$ case. Modes
associated with $m=0$ and $m\geq 3$ have $\beta_{R}\leq 0$ for all values of
$Nf$, $Eo$, and $U_{m}$ studied, and play no role in the transition to linear
instability. Furthermore, for all the cases examined, the eigenvalues
associated with the $m=1$ and $m=2$ modes are real.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 4: Growth rate, $\beta_{R}$, as a function of $U_{m}$ for $Nf=40$, 60,
80, and 100, shown in (a)-(c), (d)-(f), (g)-(i), and, (j)-(l), respectively,
with $Eo=20$ in (a), (d), (g), and (j), $Eo=180$ in (b), (e), (h), and (k),
and $Eo=300$ in (c), (f), (i), and (l). The results are shown for the modes
$m=1,2$, and $3$.

From Figure 4, it is seen that for $Nf=40,60,80,100$, increasing $Eo$ is
accompanied by a decrease in the magnitude of $U_{m}$ required for instability
and a widening of $\Delta U_{m}$ though this trend appears to saturate at
large $Eo$. Moreover, in Figure 5, which depicts the variation of $\beta_{R}$
with $Eo$ and with $Nf=80$ held constant, the critical $Eo$ for which the
$m=1$ mode is destabilised is reduced threefold as $U_{m}$ is varied from -0.2
to -0.55. Thus, the results presented in Figures 4 and 5 demonstrate that
increasing the velocity of the downward-flowing liquid and/or decreasing the
relative significance of surface tension forces is destabilising.

Inspection of Figure 4 also reveals that decreasing $Nf$ for $Eo=180$ and
$Eo=300$, that is, for weak surface tension, appears to have little effect on
the critical $U_{m}$, $\Delta U_{m}$, and the magnitude of $\beta_{R}$ for the
most dangerous mode, $m=1$. In contrast, for $Eo=20$, decreasing $Nf$ leads to
a substantial decrease in the critical $U_{m}$ value and is therefore strongly
destabilising indicating that viscous effects gain in significance as the
relative importance of surface tension increases with decreasing $Eo$ for
sufficiently low $Eo$. Furthermore, from Figure 4 (d,a), (e,b) and (f,c), it
is also seen that decreasing $Nf$ from $Nf=60$ to $Nf=40$ also leads to a
large reduction in the critical $U_{m}$ even at high $Eo$ values for
sufficiently small $Nf$.

(a)

(b)

(c)

Figure 5: Growth rate, $\beta_{R}$, as a function of $Eo$, with $Nf=80$ and
$U_{m}=-0.20$, (a), $U_{m}=-0.35$, (b) and $U_{m}=-0.55$, (c). The results are
shown for $m=1,2$, and $3$.

We offer an explanation of the trends highlighted above by focusing on the
bubble nose and appealing to the steady state results of Abubakar & Matar
(2021) who had noted that the frontal radius of curvature of the nose,
$R_{F}$, in stagnant and downward-flowing liquids increases with $Eo$ then
saturates for $Eo\gtrsim 100$ and is weakly-dependent on $Nf$ for $Nf\gtrsim
60$. For $Nf<60$, $R_{F}$ is reduced with decreasing $Nf$. For $Eo<100$,
$R_{F}$ exhibits a turning point in $Eo$ for all $Nf$ values studied with a
well-defined cross-over $Eo$ value below which the magnitude of the $R_{F}$
minima increase with decreasing $Nf$. These results demonstrate that bubble
noses become flatter with decreasing and increasing $Nf$ for sufficiently
small and large $Eo$, respectively. Abubakar & Matar (2021) have also shown
that increasingly negative $U_{m}$ has a similar effect leading to flatter
bubble noses regardless of the value of $Nf$ and $Eo$; this is attributed to
the increase in the normal stress exerted on the bubble nose relative to that
in a stagnant liquid due to the commensurate increase in the opposing inertial
force in the downward liquid flow. The results of Abubakar & Matar (2021) for
the nose curvature and its dependence on $Nf$, $Eo$, and $Um$ thus appear to
mirror the linear stability trends presented in Figures 4-5. In particular,
the significant destabilisation of the bubble with $Eo$ increasing from 20 to
180 and its subsequent saturation, the destabilising effect of $Nf$ with
decreasing $Nf$ at $Eo=20$ (see figure 4(a,j)), the weak dependence on the
stability characteristics for $Nf=40,60,80,100$ for $Eo=180,300$ (see Figure
4(b,k) and (c,l)) can be correlated to the parametric dependence of the nose
curvature $R_{F}$ on $Nf$, $Eo$, and $U_{m}$.

###  3.2 Asymmetric bubble shapes

We now study the influence of the parameters $Nf$, $Eo$, and $U_{m}$ on the
shapes of the eigenfunctions focusing on those associated with the interfacial
deformation in order to highlight the base state bubble regions targeted by
the instability. For every point on the three-dimensional axisymmetric base
state interface with position vector $\left(r^{0},\theta^{0},z^{0}\right)$ in
cylindrical coordinates, from (15) and (32c), and recalling that
$\tilde{\mathbf{x}}=\tilde{h}\mathbf{n}^{0}$, the corresponding deformed
interface points in Cartesian coordinates can be constructed:

$\displaystyle x$
$\displaystyle=r^{0}\cos\left(\theta^{0}\right)+\epsilon\left[h_{R}n_{r}\cos\left(m\theta^{0}\right)-h_{I}n_{r}\sin\left(m\theta^{0}\right)\right]\cos\left(\theta^{0}\right)\quad\forall\;\theta^{0}\;\in\;\left[0,2\pi\right]$
(50a) $\displaystyle y$
$\displaystyle=r^{0}\sin\left(\theta^{0}\right)+\epsilon\left[h_{R}n_{r}\cos\left(m\theta^{0}\right)-h_{I}n_{r}\sin\left(m\theta^{0}\right)\right]\sin\left(\theta^{0}\right)\quad\forall\;\theta^{0}\;\in\;\left[0,2\pi\right]$
(50b) $\displaystyle z$
$\displaystyle=z^{0}+\epsilon\left[h_{R}n_{z}\cos\left(m\theta^{0}\right)-h_{I}n_{z}\sin\left(m\theta^{0}\right)\right]\quad\forall\;\theta^{0}\;\in\;\left[0,2\pi\right]$
(50c)

where $h_{R}$ and $h_{I}$ denote the real and imaginary parts of the interface
deformation in the normal direction; $n_{r}$ and $n_{z}$ remain the radial and
axial components of the unit normal to the base state interface, respectively.
In our discussion below of the three-dimensional bubble shape immediately
following the transition to instability, we assign a value to the parameter
$\epsilon$, which signifies the formally infinitesimal size of the
perturbation, to enhance the visualisation of the results.

Figures 6 and 7 show the influence of the parameters $Nf$, $Eo$ and $U_{m}$ on
the eigenfunctions for the interface deformation $\hat{h}$. It is seen that
the the nose, film, and bottom regions of the bubble are targeted to varying
degrees; the precise definitions of these regions are in Abubakar & Matar
(2021). For the dominant eigenmode $m=1$, the peaks in $\hat{h}$ coincide with
the nose and bottom regions for high and low $Eo$, respectively. Around the
bottom region, the observed peaks in $\hat{h}$ are either due to the tail
structure at high $Eo$ (see the middle and right column in Figure 6), or the
undulation in the film region close to the bubble bottom at low $Eo$ (see the
left column in Figure 6) though in the latter case we note that the $m=1$ mode
is linearly stable for the parameters used to generate these results
($U_{m}=-0.2$, $Eo=20$, and $Nf=40,60,80,100$). For eigenmode $m=2$, the peak
in $\hat{h}$ coincides with the bottom region except at higher magnitude of
downward liquid flow velocity, $U_{m}$, where similar peaks are seen in the
nose region, as shown in Figure 7(c). In Figures 6 and 7, we also show
enlarged views of the variation of $\hat{h}$ with the arc length $s$ for the
bottom region to demonstrate that these boundary layer-like regions in
$\hat{h}$ have been resolved adequately.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 6: Interface deformation eigenfunctions $\hat{h}$ for eigenmodes $m=1$
and $m=2$ as a function of base state axial position of the interface,
$z^{0}$, with $U_{m}=-0.20$ and $Nf=40$, 60, 80, and 100, shown in (a)-(c),
(d)-(f), (g)-(i), and (j)-(l), respectively, with $Eo=20$ in (a), (d), (g),
and (j), $Eo=180$ in (b), (e), (h), and (k), and $Eo=300$ in (c), (f), (i),
and (l). The axisymmetric base state bubble shape is also shown (coloured
purple) as a reference in order to highlight the regions targeted by the
instability. The insets depict enlarged views of $\hat{h}$ varying with the
arc length $s$ for $m=1$ and $m=2$ in the bubble bottom region.

In Figures 8 and 9, we show the effects of $U_{m}$, $Eo$, and $Nf$ on the
three-dimensional bubble shapes obtained by adding the interface deformation
associated with mode $m=1$ to the base state taking the value of
$\epsilon=0.05$ for clarity of presentation; these correspond to the shapes
one might expect to observe experimentally at the onset of instability. Also
shown in Figures 8 and 9 are two-dimensional projections of the shapes in the
$(r,z)$ plane which highlight the deviations from the base state within the
framework of linear theory. The results shown in Figure 8 are for the same
parameter values used to generate Figure 7, which correspond to the large
$Eo$, weak surface tension limit; panel (d) of this figure also depicts the
analogous results for the deformed bubble shape when $m=2$. These results
illustrate the asymmetry of the nose for Taylor bubble motion in downward
flowing liquids for negligible surface tension and are reminiscent of the
asymmetric shape shown in Figure 1(b). In Figure 9, on the other hand, the
results are associated with $Eo=20$ at which surface tension effects are
significant. Here, it is clearly seen that in addition to asymmetries in the
nose region, the instability also targets the undulation in the bottom region.
Figure 9(d) provides a clear demonstration that the asymmetry is most
pronounced in this region for the fastest downward flowing liquid case; in
contrast, the bubble nose remains essentially axisymmetric in this case.

(a)

(b)

(c)

Figure 7: Interface deformation eigenfunctions $\hat{h}$ for eigenmodes $m=1$
and $m=2$ as a function of base state axial position of the interface,
$z^{0}$, for $Nf=100$, $Eo=220$, and with $U_{m}=-0.25$, (a), $U_{m}=-0.40$,
(b), and $U_{m}=-0.55$, (c). The axisymmetric base state bubble shape is also
shown (coloured purple) in order to highlight the regions targeted by the
instability. The insets depict enlarged views of $\hat{h}$ varying with the
arc length $s$ for $m=1$ and $m=2$ in the bubble bottom region.

(a)

(b)

(c)

(d)

Figure 8: Three-dimensional bubble shapes and their two-dimensional
projections (solid lines) in the $(r,z)$ plane obtained by adding the
interface deformation for $m=1$ to the base state (dashed lines) with
$\epsilon=0.05$ for $U_{m}=-0.25,-0.40,-0.55$, shown in (a)-(c), respectively,
with $Nf=100$ and $Eo=220$. In (d), we show the analogous shape for $m=2$ with
$U_{m}=-0.55$, $Nf=100$, and $Eo=220$.

###  3.3 Stability maps

We show in Figure 10 stability maps in $(U_{m},Eo)$ with $Nf$ varying
parametrically that depict the boundaries demarcating regions of linear
instability for Taylor bubbles moving in downward flowing liquids
characterised by $U_{m}<0$. In each case, examples of the three-dimensional,
asymmetric bubble shape at the onset of instability is also shown. The general
trend observed is that for $Eo>60$, the magnitude of the critical $U_{m}$
decreases with increasing $Eo$ for a fixed $Nf$, saturating for large $Eo$ and
$Nf$, beyond $Eo\gtrsim 100$ and $Nf\gtrsim 60$. For $Eo<60$, there is a
turning point in the stability map that highlights the fact that an increase
in the magnitude of the downward flow velocity is needed to overcome the
strong surface tension forces at sufficiently low $Eo$; furthermore, the
critical $U_{m}$ magnitude increases with $Nf$ reflecting the destabilising
effect of viscous stresses over this range of $Eo$. The results presented in
Figure 10 are consistent with the trends discussed in the previous sections,
and, in particular, those associated with the normalised frontal radius of
curvature of the nose, $R_{F}$, presented in Abubakar & Matar (2021); this
provides a further indication that the dependence of $R_{F}$ on $U_{m}$, $Eo$,
and $Nf$ controls the linear stability characteristics of the bubble motion.

We also show in Figure 10 the curve for the critical $U_{m}$ for which the
axisymmetric base state has $U_{b}=0$. It is noticeable that for $Eo\gtrsim
100$ for all $Nf$ studied, this curve is in the linearly unstable region
indicating that bubbles whose motion has been arrested due to a downward
flowing liquid over this range of $Eo$ and $Nf$ cannot have an axisymmetric
shape. In contrast, in the complementary range of $Eo$, the critical $U_{m}$
curve for such bubbles is outside the unstable region implying that they can
sustain an axisymmetric shape despite the downward liquid flow. Examples of
these cases have been observed experimentally; see, for instance, the
axisymmetric Taylor bubble shown to have been held stationary in downward
liquid flow by Nigmatulin (2001), for $Nf=6087.38$ and $Eo=33.06$. We note
that the $Eo$ value is within the range of the linearly stable region shown in
Figure 10(e). Though the experimental $Nf$ value is outside of the range we
studied, the fact that the linear stability boundaries appear to saturate at
large $Nf$ suggests that our results can still provide a reasonable indication
of the behaviour observed experimentally.

##  4 Energy budget analysis

In this section, we analyse how energy is transferred from the base flow to
the perturbations by studying the growth of the perturbation kinetic energy
(Hu & Joseph, 1989; Hooper & Boyd, 1983). By investigating the contribution of
the mechanisms of different physical origin that account for energy
production, one can identify the dominant ones that drive instability
(Boomkamp & Miesen, 1996). This analysis has been used to study destabilising
mechanisms in parallel two-phase flows (Sahu et al., 2009; Selvam et al.,
2007; Sahu et al., 2007; Ó Náraigh et al., 2011) and their classification
(Boomkamp & Miesen, 1996).

###  4.1 Energy balance formulation

To derive the equation for the growth of the disturbance kinetic energy, one
multiplies the continuous forms of the momentum perturbation equations for the
velocity components with their corresponding complex conjugates, integrates
over the domain, adds the resulting equations, and simplifies as appropriate.
Our perturbation equations for the velocity components, however, were derived
from the weak form of the continuous momentum equations written around the
perturbed domain; the development of the weak form of the energy equation must
therefore follow the same strategy. Thus, we obtain the energy equation from
the derived perturbation equations for the velocity components by setting the
test functions for the latter to equal the complex conjugate of the velocity
perturbations, $\bar{{\Phi}}_{r}=\hat{{u}}_{r}^{*}$,
$\bar{{\Phi}}_{\theta}=\hat{{u}}_{\theta}^{*}$,
$\bar{{\Phi}}_{z}=\hat{{u}}_{z}^{*}$, followed by necessary simplifications:

$\displaystyle\beta_{R}$
$\displaystyle\int_{\Omega^{0}}\left\\{\left|\hat{u}_{r}\right|^{2}+\left|\hat{u}_{\theta}\right|^{2}+\left|\hat{u}_{z}\right|^{2}\right\\}d\Omega^{0}+\int_{\Omega^{0}}\left\\{\left|\hat{u}_{r}\right|^{2}\frac{\partial
u_{r}^{0}}{\partial r^{0}}+\left|\hat{u}_{z}\right|^{2}\frac{\partial
u_{z}^{0}}{\partial
z}+\frac{\left|\hat{u}_{\theta}\right|^{2}u_{r}^{0}}{r}\right.$
$\displaystyle\left.+\left(\frac{\partial u_{r}^{0}}{\partial
z}+\frac{\partial u_{z}^{0}}{\partial
r^{0}}\right)\mathbb{R}\left\\{\hat{u}_{r}\hat{u}_{z}^{*}\right\\}\right\\}d{\Omega^{0}}+\int_{\Omega}Nf^{-1}\left\\{2\left|\frac{\partial\hat{u}_{r}}{\partial
r}\right|^{2}+2\left|\frac{\partial\hat{u}_{z}}{\partial z}\right|^{2}\right.$
$\displaystyle\left.+\frac{1}{r^{2}}\left[\left(2+m^{2}\right)\left|\hat{u}_{r}\right|^{2}+\left(1+2m^{2}\right)\left|\hat{u}_{\theta}\right|^{2}+m^{2}\left|\hat{u}_{z}\right|^{2}\right]+\left|\frac{\partial\hat{u}_{\theta}}{\partial
r}\right|^{2}+\left|\frac{\partial\hat{u}_{\theta}}{\partial
z}\right|^{2}\right.$
$\displaystyle\left.+\left|\frac{\partial\hat{u}_{r}}{\partial
z}+\frac{\partial\hat{u}_{z}}{\partial
r}\right|^{2}-\frac{1}{r}\frac{\partial\left|\hat{u}_{\theta}\right|^{2}}{\partial
r}+\frac{6m}{{r}^{2}}\mathbb{I}m\left\\{\hat{u}_{r}\hat{u}_{\theta}^{*}\right\\}+\frac{2m}{r}\mathbb{I}m\left(\hat{u}_{z}^{*}\frac{\partial\hat{u}_{\theta}}{\partial
z}+\hat{u}_{r}^{*}\frac{\partial\hat{u}_{\theta}}{\partial
r}\right)\right\\}d{\Omega^{0}}$
$\displaystyle-\int_{\Gamma_{b}^{0}}Eo^{-1}\left\\{-\frac{d\hat{h}}{ds}\left[\mathbf{n}\cdot\frac{d\hat{\mathbf{u}}^{*}}{ds}-\kappa_{a}\left(\mathbf{t}\cdot\hat{\mathbf{u}}^{*}\right)\right]+\hat{h}\left[\kappa_{a}^{2}+\kappa_{b}^{2}-\frac{m^{2}}{r^{2}}\right]\mathbf{n}\cdot\hat{\mathbf{u}}^{*}\right\\}d\Gamma_{b}^{0}$
$\displaystyle-\int_{\Gamma_{b}^{0}}\left\\{\hat{h}n_{z}\left(\mathbf{n}\cdot\hat{\mathbf{u}}^{*}\right)\right\\}d\Gamma_{b}^{0}$
$\displaystyle+\int_{\Gamma^{0}_{b}}\hat{h}\left\\{{\left(\mathbf{u}^{0}\cdot\mathbf{t}\right)\left[\frac{d\mathbf{u}^{0}}{ds}\cdot\hat{\mathbf{u}}^{*}\right]}+\left[-p^{0}+2Nf^{-1}\left(\mathbf{t}\cdot\frac{d\mathbf{u}^{0}}{ds}\right)\right]\left(\mathbf{t}\cdot\frac{d\hat{\mathbf{u}}^{*}}{ds}\right)\right.$
$\displaystyle\left.+\left[-p^{0}+2Nf^{-1}\frac{u_{r}^{0}}{r}\right]\left(\frac{\hat{u}_{r}^{*}}{r}-\mathfrak{i}m\frac{\hat{u}_{\theta}^{*}}{r}\right)\right\\}d\Gamma^{0}_{b}$
$\displaystyle+\int_{\Gamma_{b}^{0}}\left\\{\left[{Eo}^{-1}\kappa+z-P_{b}^{0}\right]\left[\left(\mathbf{t}\cdot\hat{\mathbf{u}}^{*}\right)\frac{d\hat{h}}{ds}+\hat{h}\left(\mathfrak{i}m\frac{\hat{u}_{\theta}^{*}}{r}\right)+\hat{h}\kappa\left(\mathbf{n}\cdot\hat{\mathbf{u}}^{*}\right)\right]\right\\}d\Gamma_{b}^{0}=0$
(51)

where the symbol $|.|$ represents the magnitude of a complex function;
$\mathbb{R}$ and $\mathbb{I}$ denote the real and imaginary part of a complex
function, respectively. Equation (51) is the energy budget formulation that
governs the evolution of the disturbance kinetic equation.

(a)

(b)

(c)

(d)

Figure 9: Three-dimensional bubble shapes and their two-dimensional
projections (solid lines) in the $(r,z)$ plane obtained by adding the
interface deformation for $m=1$ to the base state (dashed lines) with
$\epsilon=0.05$ for $U_{m}=-0.35,-0.55,-0.80,-1.10$, and $Nf=40,60,80,100$,
shown in (a)-(d), respectively, with $Eo=20$.

(a)

(b)

(c)

(d)

(e)

Figure 10: Stability maps depicting the boundaries demarcating (shown by the
full circles) the regions of linear instability in $(U_{m},Eo)$ space
characterised by a transition from axisymmetric to asymmetric bubble shapes in
downward liquid flow with $U_{m}<0$ for (a) $Nf=40$, (b) $Nf=60$, (c) $Nf=80$,
(d) $Nf=100$ (e) $Nf=120$. The dashed lines represent the curves of critical
$U_{m}$ for which the axisymmetric base state corresponds to an arrested
bubble in downward liquid flow with $U_{b}=0$.

Following Boomkamp & Miesen (1996), we express the energy balance equation as

$\dot{E}=REY+DIS+INT,$ (52)

where $\dot{E}$ corresponds to the time range of change of the perturbation
kinetic energy given by the following relation

$\dot{E}=\beta_{R}\int_{\Omega^{0}}\left\\{\left|\hat{u}_{r}\right|^{2}+\left|\hat{u}_{\theta}\right|^{2}+\left|\hat{u}_{z}\right|^{2}\right\\}d\Omega^{0}\equiv\beta_{R}~{}KIN,$
(53)

wherein $KIN$ represents the total kinetic energy associated with the
perturbation velocity field, which equals $\dot{E}$ when multiplied by the
growth rate $\beta_{R}$ (which is positive for an unstable flow). We also
introduce the following definitions for the terms $REY$ and $DIS$ that appear
on the right-hand-side of equation (52) (Boomkamp & Miesen, 1996):

$REY=-\int_{\Omega^{0}}\left\\{\left|\hat{u}_{r}\right|^{2}\frac{\partial
u_{r}^{0}}{\partial r^{0}}+\left|\hat{u}_{z}\right|^{2}\frac{\partial
u_{z}^{0}}{\partial
z}+\frac{\left|\hat{u}_{\theta}\right|^{2}u_{r}^{0}}{r}+\left(\frac{\partial
u_{r}^{0}}{\partial z}+\frac{\partial u_{z}^{0}}{\partial
r^{0}}\right)\mathbb{R}\left\\{\hat{u}_{r}\hat{u}_{z}^{*}\right\\}\right\\}d{\Omega^{0}},\\\
$ (54) $\displaystyle DIS=$
$\displaystyle-\int_{\Omega}Nf^{-1}\left\\{2\left|\frac{\partial\hat{u}_{r}}{\partial
r}\right|^{2}+2\left|\frac{\partial\hat{u}_{z}}{\partial
z}\right|^{2}+\frac{1}{r^{2}}\left[\left(2+m^{2}\right)\left|\hat{u}_{r}\right|^{2}+\left(1+2m^{2}\right)\left|\hat{u}_{\theta}\right|^{2}\right.\right.$
$\displaystyle\left.\left.+m^{2}\left|\hat{u}_{z}\right|^{2}\right]+\left|\frac{\partial\hat{u}_{\theta}}{\partial
r}\right|^{2}+\left|\frac{\partial\hat{u}_{\theta}}{\partial
z}\right|^{2}+\left|\frac{\partial\hat{u}_{r}}{\partial
z}+\frac{\partial\hat{u}_{z}}{\partial
r}\right|^{2}-\frac{1}{r}\frac{\partial\left|\hat{u}_{\theta}\right|^{2}}{\partial
r}+\frac{6m}{{r}^{2}}\mathbb{I}m\left\\{\hat{u}_{r}\hat{u}_{\theta}^{*}\right\\}\right.$
$\displaystyle\left.+\frac{2m}{r}\mathbb{I}m\left(\hat{u}_{z}^{*}\frac{\partial\hat{u}_{\theta}}{\partial
z}+\hat{u}_{r}^{*}\frac{\partial\hat{u}_{\theta}}{\partial
r}\right)\right\\}d{\Omega}.$ (55)

Here, $REY$ denotes the rate of energy transfer by the Reynolds stress from
the base flow to the disturbed flow, and $DIS$ represents the rate of viscous
dissipation of energy of the disturbed flow. We also provide a breakdown for
$INT$, the rate of work done by the velocity and stress disturbances in
deforming the interface (Boomkamp & Miesen, 1996)

$\displaystyle INT=NOR+TAN,$ (56a) $\displaystyle NOR=TEN+HYD+BUB,$ (56b)

which we have decomposed into its normal, $NOR$, and tangential, $TAN$,
components with $NOR$ further subdivided into $TEN$, $HYD$, and $BUB$,
representing work done at the interface against surface tension, gravity, and
bubble pressure, respectively. Expressions for $TAN$, $NOR$, and $TEN$, $HYD$
and $BUB$, are respectively provided by

$\displaystyle TAN=$
$\displaystyle\mathbb{R}\left\\{-\int_{\Gamma^{0}_{b}}\hat{h}\left\\{{\left(\mathbf{u}^{0}\cdot\mathbf{t}\right)\left[\frac{d\mathbf{u}^{0}}{ds}\cdot\hat{\mathbf{u}}^{*}\right]}+\left[-p^{0}+2Nf^{-1}\left(\mathbf{t}\cdot\frac{d\mathbf{u}^{0}}{ds}\right)\right]\left(\mathbf{t}\cdot\frac{d\hat{\mathbf{u}}^{*}}{ds}\right)\right.\right.$
$\displaystyle\left.\left.+\left[-p^{0}+2Nf^{-1}\frac{u_{r}^{0}}{r}\right]\left(\frac{\hat{u}_{r}^{*}}{r}-\mathfrak{i}m\frac{\hat{u}_{\theta}^{*}}{r}\right)\right\\}d\Gamma^{0}_{b}+\int_{\Gamma_{b}^{0}}Eo^{-1}\left\\{\frac{d\hat{h}}{ds}\left[\kappa_{a}\left(\mathbf{t}\cdot\hat{\mathbf{u}}^{*}\right)\right]\right\\}d\Gamma_{b}^{0}\right.$
$\displaystyle\left.-\int_{\Gamma_{b}^{0}}\left\\{\left[{Eo}^{-1}\kappa+z-P_{b}^{0}\right]\left[\left(\mathbf{t}\cdot\hat{\mathbf{u}}^{*}\right)\frac{d\hat{h}}{ds}+\hat{h}\left(\mathfrak{i}m\frac{\hat{u}_{\theta}^{*}}{r}\right)\right]\right\\}d\Gamma_{b}^{0}\right\\},$
(57) $\displaystyle NOR=$
$\displaystyle\mathbb{R}\left\\{\int_{\Gamma_{b}^{0}}Eo^{-1}\left\\{-\frac{d\hat{h}}{ds}\left[\mathbf{n}\cdot\frac{d\hat{\mathbf{u}}^{*}}{ds}\right]+\hat{h}\left[\kappa_{a}^{2}+\kappa_{b}^{2}-\frac{m^{2}}{r^{2}}\right]\mathbf{n}\cdot\hat{\mathbf{u}}^{*}\right\\}d\Gamma_{b}^{0}\right.$
$\displaystyle\left.+\int_{\Gamma_{b}^{0}}\left\\{\hat{h}n_{z}\left(\mathbf{n}\cdot\hat{\mathbf{u}}^{*}\right)\right\\}d\Gamma_{b}^{0}-\int_{\Gamma_{b}^{0}}\left\\{\left[{Eo}^{-1}\kappa+z-P_{b}^{0}\right]\left[\hat{h}\kappa\left(\mathbf{n}\cdot\hat{\mathbf{u}}^{*}\right)\right]\right\\}d\Gamma_{b}^{0}\right\\},$
(58)
$TEN=\mathbb{R}\left\\{\int_{\Gamma_{b}^{0}}Eo^{-1}\left\\{-\frac{d\hat{h}}{ds}\left[\mathbf{n}\cdot\frac{d\hat{\mathbf{u}}^{*}}{ds}\right]+\hat{h}\left[\kappa_{a}^{2}+\kappa_{b}^{2}-\kappa^{2}-\frac{m^{2}}{r^{2}}\right]\mathbf{n}\cdot\hat{\mathbf{u}}^{*}\right\\}d\Gamma_{b}^{0}\right\\},$
(59)
$HYD=\mathbb{R}\left\\{\int_{\Gamma_{b}^{0}}\left\\{\hat{h}\left(n_{z}-z\kappa\right)\left(\mathbf{n}\cdot\hat{\mathbf{u}}^{*}\right)\right\\}d\Gamma_{b}^{0}\right\\},$
(60)
$BUB=\mathbb{R}\left\\{\int_{\Gamma_{b}^{0}}\left\\{P_{b}^{0}\left[\hat{h}\kappa\left(\mathbf{n}\cdot\hat{\mathbf{u}}^{*}\right)\right]\right\\}d\Gamma_{b}^{0}\right\\}.$
(61)

We normalised the energy terms using $KIN$, so the energy balance equation
becomes

$\beta_{R}=REY^{*}+DIS^{*}+TEN^{*}+HYD^{*}+BUB^{*}+TAN^{*},$ (62)

where the asterisk designates the normalisation by $KIN$. In Table 1, we
demonstrate for eigenmode $m=1$, $Nf=80$, and $Eo=140$, and various $U_{m}$
values, that the difference between the growth rate $\beta_{R}$ on the left-
hand-side of (62) computed from the linear stability analysis and the sum of
the energy terms on the right-hand-side of (62), $SUM$, denoted as $DIF$ in
Table 1 is negligibly small. These results inspire confidence in our procedure
for computing the terms in (62).

Table 1: Balance of energy distribution among the various energy terms in equation (62) for $Nf=80$, $Eo=140$, and mode $m=1$. $U_{m}$ | $\beta_{R}$ | $REY^{*}$ | $DIS^{*}$ | $TEN^{*}$ | $HYD^{*}$ | $BUB^{*}$ | $TAN^{*}$ | $SUM$ | $DIF$
---|---|---|---|---|---|---|---|---|---
$-0.40$ | $0.2824$ | $-0.0059$ | $-2.6795$ | $-0.1070$ | $0.3066$ | $5.9845$ | $-3.21631$ | $0.2824$ | $2.78\times 10^{-8}$
$-0.30$ | $0.1128$ | $0.1224$ | $-2.6494$ | $0.0528$ | $0.5397$ | $3.1162$ | $-1.0690$ | $0.1128$ | $1.85\times 10^{-7}$
$-0.20$ | $-0.6354$ | $0.2656$ | $-2.6537$ | $0.2919$ | $0.8280$ | $0.9208$ | $0.2840$ | $-0.6354$ | $4.2.0\times 10^{-8}$
$-0.10$ | $-0.2442$ | $0.3260$ | $-2.6008$ | $0.5245$ | $1.2070$ | $-0.5126$ | $0.8117$ | $-0.2442$ | $2.23\times 10^{-6}$
$0.00$ | $-0.4474$ | $0.2269$ | $-2.1745$ | $0.5412$ | $1.1340$ | $-0.9108$ | $0.7358$ | $-0.4473$ | $1.40\times 10^{-5}$
$0.10$ | $-0.6404$ | $0.0413$ | $-1.6299$ | $0.4210$ | $0.6429$ | $-0.5639$ | $0.4481$ | $-0.6404$ | $-7.67\times 10^{-6}$

(a)

(b)

(c)

Figure 11: Breakdown of the energy budget for the eigenmode $m=1$ as a
function of $U_{m}$ for $Eo=20$, 180, and 300, shown in (a)-(c), respectively,
with $Nf=80$. In each panel, the vertical dashed line marks the $U_{m}$ value
for which $\beta_{R}=0$.

###  4.2 Energy analysis results

From the linear stability analysis results for downward liquid flow presented
in section 3, the $m=1$ mode was identified as being the most unstable one.
Here, we examine the contribution of each term in equation (62) in order to
elucidate their roles in driving instability and to identify the most dominant
destabilising mechanism. In Figure 11, the energy analysis results for $Nf=80$
with $Eo=20,180,$ and $300$ are shown. These parameter values are chosen to
correspond to those used to generate a representative subset of the results
presented in section 3. In the discussion below, the asterisk decoration which
appears in equation (62) is suppressed for the sake of brevity.

It is seen clearly in Figure 11(a) that for $Eo=20$, $TAN$ is overwhelmingly
the dominant mechanism, followed by $BUB$, over the range of $U_{m}$ for which
$\beta_{R}>0$. Over the remainder of the $U_{m}$ range studied in which
$\beta_{R}<0$, $TAN$ decreases monotonically but remains destabilising while
$BUB$ becomes stabilising. $HYD$, $REY$, and $TEN$ are destabilising over this
range while, as expected, $DIS$ is stabilising. Upon increasing $Eo$ to $180$
and $300$, as shown in Figure 11(b) and 11(c), respectively, the dominant
mechanism over the unstable range of $U_{m}$ switches to $BUB$ followed in
relative dominance by $HYD$, which is marginally destabilising, while $TAN$ is
strongly stabilising. These results suggest that the dominant mechanism
depends on the relative significance of surface tension forces. In the case of
negligible surface tension, characterised by high $Eo$ values, the
destabilising mechanism is related to the bubble pressure which indicates that
the origin of the instability is in the gas phase. At low $Eo$, the
instability originates in the liquid phase with the energy provided by the
tangential stress component.

It is instructive to split $TAN$, given by equation (57), into its constituent
components based on the base state groups that supply energy to the
perturbations:

$\displaystyle
TAN_{ut}=\mathbb{R}\left\\{-\int_{\Gamma^{0}_{b}}\hat{h}{\left(\mathbf{u}^{0}\cdot\mathbf{t}\right)\left[\frac{d\mathbf{u}^{0}}{ds}\cdot\hat{\mathbf{u}}^{*}\right]}\right\\},$
(63a) $\displaystyle
TAN_{strs}=\mathbb{R}\left\\{-\int_{\Gamma^{0}_{b}}\hat{h}\left\\{\left[-p^{0}+2Nf^{-1}\left(\mathbf{t}\cdot\frac{d\mathbf{u}^{0}}{ds}\right)\right]\left(\mathbf{t}\cdot\frac{d\hat{\mathbf{u}}^{*}}{ds}\right)\right.\right.$
$\displaystyle\left.\left.\qquad{}\qquad{}\qquad{}+\left[-p^{0}+2Nf^{-1}\frac{u_{r}^{0}}{r}\right]\left(\frac{\hat{u}_{r}^{*}}{r}-\mathfrak{i}m\frac{\hat{u}_{\theta}^{*}}{r}\right)\right\\}\right\\},$
(63b) $\displaystyle
TAN_{ts}=\mathbb{R}\left\\{\int_{\Gamma_{b}^{0}}Eo^{-1}\left\\{\frac{d\hat{h}}{ds}\left[\kappa_{a}\left(\mathbf{t}\cdot\hat{\mathbf{u}}^{*}\right)\right]-\kappa\left[\left(\mathbf{t}\cdot\hat{\mathbf{u}}^{*}\right)\frac{d\hat{h}}{ds}+\hat{h}\left(\mathfrak{i}m\frac{\hat{u}_{\theta}^{*}}{r}\right)\right]\right\\}d\Gamma_{b}^{0}\right\\},$
(63c) $\displaystyle
TAN_{g}=\mathbb{R}\left\\{\int_{\Gamma_{b}^{0}}-z\left[\left(\mathbf{t}\cdot\hat{\mathbf{u}}^{*}\right)\frac{d\hat{h}}{ds}+\hat{h}\left(\mathfrak{i}m\frac{\hat{u}_{\theta}^{*}}{r}\right)\right]d\Gamma_{b}^{0}\right\\},$
(63d) $\displaystyle
TAN_{pb}=\mathbb{R}\left\\{\int_{\Gamma_{b}^{0}}P_{b}\left[\left(\mathbf{t}\cdot\hat{\mathbf{u}}^{*}\right)\frac{d\hat{h}}{ds}+\hat{h}\left(\mathfrak{i}m\frac{\hat{u}_{\theta}^{*}}{r}\right)\right]d\Gamma_{b}^{0}\right\\}.$
(63e)

The terms $TAN_{ut}$, $TAN_{strs}$, $TAN_{ts}$, $TAN_{g}$ and $TAN_{pb}$
denote the contributions to $TAN$ due to streaming tangential velocity,
tangential stress, surface tension, gravity and bubble pressure on the
interface as captured by the base state terms $\mathbf{u}\cdot\mathbf{t}$,
$\left[-p^{0}+2Nf^{-1}\left(\mathbf{t}\cdot\frac{d\mathbf{u}^{0}}{ds}\right)\right]+\left[-p^{0}+2Nf^{-1}\frac{u_{r}^{0}}{r}\right]$,
$Eo$, $z$ and $P_{b}$ in the expressions, respectively. The base state
contribution to $TAN_{strs}$ is related to the tangential stress because it
can be obtained by taking the double dot product of the stress tensor and the
tangential projection operator,
$\left(\mathbf{I}-\mathbf{n}\otimes\mathbf{n}\right)$.

In Figure 12, we plot the dependence of the constituents of $TAN$ on $U_{m}$
for $m=1$, $Eo=20$, 180, and 300, and $Nf=80$. Inspection of Figure 12(a)
reveals that in the $Eo=20$ case, for which surface tension effects are
important, the major contributor to $TAN$ corresponds to $TAN_{strs}$ and
exerts a destabilising influence on the bubble motion.

(a)

(b)

(c)

Figure 12: Breakdown of $TAN$ for the eigenmode $m=1$ into its constituent
components given by equations (63) as a function of $U_{m}$ for $Eo=20$, 180,
and 300, shown in (a)-(c), respectively, with $Nf=80$. In each panel, the
vertical dashed line marks the $U_{m}$ value for which $\beta_{R}=0$.

As can also be seen in Figure 12a, although $TAN_{strs}$ remains
destabilising, it gives way to $TAN_{g}$ as the dominant contributor to $TAN$
with decreasing magnitude of $U_{m}$, while $TAN_{pb}$ is sufficiently
stabilising so as to render $\beta_{R}<0$; the contributions of $TAN_{ut}$ and
$TAN_{ts}$ are relatively negligible and they play an insignificant role in
the bubble stability.

In Figure 12(b,c) generated for $Eo=180$ and 300 for which surface tension
effects are weak, the dominant destabilising contribution to $TAN$ is due to
$TAN_{g}$ with the sub-dominant $TAN_{strs}$ and $TAN_{pb}$ exerting a
stabilising influence over the majority of the $U_{m}$ range investigated. The
reversal in the role of $TAN_{strs}$ as we cross over from relatively low to
high $Eo$ values shown in Figure 12(b,c) is consistent with the results
discussed in the previous sections which indicated that viscous effects are
destabilising (stabilising) for low (high) $Eo$. This is further illustrated
in Figure 13 in which we plot the breakdown of the energy budget (see Figure
13a) and the constituents of $TAN$ (see Figure 13b) as a function of $Eo$ for
$m=1$ with $U_{m}=-0.2$ and $Nf=80$. It is clear that $TAN_{strs}$ switches
roles in the $Eo$ interval $\left(60,100\right]$ and $TAN$ exhibits a similar
behaviour over a somewhat larger $Eo$ range.

(a)

(b)

Figure 13: Breakdown of the energy budget (see equation (62)), (a), and the
$TAN$ constituents (see equation (63)), (b), with $Eo$, for $m=1$ with $Nf=80$
and $U_{m}=-0.20$. The vertical dashed line marks the $Eo$ value for which
$\beta_{R}=0$.

Lastly, we show in Figure 14a breakdown of the energy budget and of the $TAN$
constituents as a function of $Nf$ for $m=1$ with $Eo=300$ and $U_{m}=-0.25$.
It is seen clearly in Figure 14a that for this large $Eo$ case, $BUB$ provides
the dominant destabilising contribution with $TAN$ and $DIS$ inducing
stability. Inspection of Figure 14b reveals that although $TAN_{g}$ is
destabilising over the range of $Nf$ studied, $TAN_{pb}$ is also destabilising
for $Nf<60$; this acts to reduce the stabilising effect associated with the
increase in viscous effects and reduction in $Nf$.

(a)

(b)

Figure 14: Breakdown of the energy budget (see equation (62)), (a), and the
$TAN$ constituents (see equation (63)), (b), with $Nf$, for $m=1$ with
$Eo=300$ and $U_{m}=-0.25$.

We now establish a connection with the work of Lu & Prosperetti (2006) who
concluded that it is the normal component of gravity on the interface that
drives the transition to asymmetric Taylor bubble shape. It is worth
mentioning, however, that the analysis of Lu & Prosperetti (2006) was carried
out locally around the nose region under the assumptions that the effects of
viscosity and surface tension are negligible. In contrast, our analysis shows
that the dominant destabilising mechanisms depend on the relative significance
of surface tension characterised by $Eo$: for low and high $Eo$, the
tangential stress, $TAN$ (with $TAN_{strs}$ being the main contributor) and
the work done at the interface against the bubble pressure, $BUB$, are chiefly
responsible for instability, respectively. A look at Figures 13(a) and 14(a)
shows that the energy term due to normal component of gravity on the
interface, $HYD$, is an increasing function of $Nf$ and $Eo$. It is plausible
that at very high $Nf$ and $Eo$ $HYD$ may overtake $BUB$ as the most dominant
destabilising energy term. It is therefore likely that the mechanisms
governing the instability in regimes of extremely negligible surface tension
and viscosity is different from that identified in this investigation.

##  5 Summary and conclusions

We have examined the linear stability of Taylor bubbles in stagnant and
flowing liquids in vertical pipes focusing on the case of downward liquid
flow. The base state, characterised by constant bubble and axisymmetric
shapes, was computed by Abubakar & Matar (2021) as a function of the Eötvös
and inverse viscosity numbers, $Eo$ and $Nf$, and the (centreline) speed of
the downward flowing liquid, $U_{m}$. A finite element linear stability model
was derived using the concepts of domain perturbation and total linearisation
method presented in Carvalho & Scriven (1999) and Kruyt et al. (1988),
respectively. The model was validated by comparing its predictions with
analytical results for the growth rate and frequency for the normal mode
small-amplitude oscillation of a spherical bubble in an unbounded stagnant
liquid under the assumption of negligible gravity.

Our linear stability framework was then used to examine the stability of the
base states obtained by Abubakar & Matar (2021). Our results demonstrated that
the leading unstable mode corresponds to $m=1$, where $m$ is the azimuthal
wavenumber of the applied perturbation. We constructed stability maps showing
the dependence of the critical magnitude of $U_{m}$ on $Eo$, with $Nf$ varying
parametrically, which demarcate the regions in $(U_{m},Eo)$ space wherein the
flow is linearly unstable. At low $Eo$, for which surface tension effects are
significant, the instability targets an undulation in the bottom region of the
bubble with the three-dimensional bubble shape exhibiting an asymmetric bulge
in that region. For weak surface tension effects, characterised by high $Eo$,
the most unstable mode corresponds to a deflection of the bubble nose away
from the axis of symmetry. The stability maps also show the locus of points
for which the bubble are stationary in a downward flowing liquid and
highlights the regions in parameter space in which they are linearly stable or
unstable resulting in axisymmetric and asymmetric shapes, respectively.

To elucidate the origins of the transition to linear instability and
asymmetric bubble shapes, an energy budget analysis was performed to analyse
the contribution of various physical mechanisms to the production of
perturbation energy. This analysis showed that the major contribution to
energy production that drives the instability comes from the bubble pressure
and the tangential stress for high and low $Eo$ values, respectively. The
insights gained from the energy analysis, and the trends observed in the
linear stability characteristics, were used to establish clear connections to
the influence of $Eo$, $Nf$, and $U_{m}$ on the curvature of the bubble nose,
which plays a crucial role in the stability of the flow.

Declaration of Interests. The authors report no conflict of interest.

### Acknowledgements

This work was supported by the Engineering $\&$ Physical Sciences Research
Council UK through the EPSRC MEMPHIS (grant no. EP/K003976/1) and PREMIERE
(grant no. EP/T000414/1) Programme Grants, and a Nigerian government PTDF
scholarship for HA. OKM also acknowledges funding from PETRONAS and the Royal
Academy of Engineering for a Research Chair in Multiphase Fluid Dynamics, and
the PETRONAS Centre for Engineering of Multiphase Systems. We also acknowledge
HPC facilities provided by the Research Computing Service (RCS) of Imperial
College London for the computing time.

##  Appendix A Weak formulations

The transformation of the governing equations into their weak forms involves
three steps: multiplying the governing equations for each variables with their
corresponding test functions and integrating over the domain, integrating by
part to reduce the order of integration, and, finally, incorporating the
boundary conditions into the resulting relations. Before proceeding to derive
the weak forms of the equations, it is important to define the necessary
functional spaces to which the solution and test functions must belong
(Heinrich & Pepper, 1999):

  1. 1.

The $L^{2}\left(\Omega\right)$ space: This is a space of functions
$f\left(\mathbf{r}\right)$ defined in $\Omega$ that are square integrable over
$\Omega$:

$L^{2}\left(\Omega\right)=\left\\{f\left(\mathbf{r}\right)|\int_{\Omega}\left(f\left(\mathbf{r}\right)\right)^{2}d\Omega<\infty\right\\}.$
(64)

  2. 2.

The $L^{2}_{0}\left(\Omega\right)$ subspace: This is a subspace of
$L^{2}\left(\Omega\right)$ defined in $\Omega$ such that for functions defined
in $L^{2}\left(\Omega\right)$, the following equation is satisfied:

$L^{2}_{0}\left(\Omega\right)=\left\\{f|f\in L^{2}\left(\Omega\right)\mbox{
and }\int_{\Omega}f\left(\mathbf{r}\right)=\mathrm{0}\right\\}.$ (65)

  3. 3.

The Sobolev space $H^{1}\left(\Omega\right)$: this is a space of functions
$f\left(\mathbf{r}\right)$ defined in $\Omega$ such that both the function and
all its first partial derivatives are in $L^{2}\left(\Omega\right)$

$H^{1}\left(\Omega\right)=\left\\{f\left(\mathbf{r}\right)|\int_{\Omega}\left[|f|^{2}+|\nabla
f|^{2}\right]d\Omega<\infty\right\\}.$ (66)

  4. 4.

The Sobolev subspace $H^{1}_{0}\left(\Omega\right)$: this is a subspace of the
Sobolev space $H^{1}\left(\Omega\right)$ for which the functions defined in
space $H^{1}$ vanish on the portions of the boundary of $\Omega$ where
Dirichilet boundary conditions are imposed ( i.e
$\Gamma_{D}=\Gamma_{in}+\Gamma_{wall}$ ):

$H^{1}_{0}\left(\Omega\right)=\left\\{f|f\in H^{1}\left(\Omega\right)\mbox{
and }f\left(\mathbf{r}\right)=\mathbf{0}\mbox{ if
}\mathbf{r}\in\Gamma_{D}\right\\}.$ (67)

Let $\mathbf{\Phi}\in H^{1}_{0}$ and $\varphi\in L^{2}_{0}$ be the test
functions corresponding to $\mathbf{u}\in H^{1}$ and $\mathrm{p}\in L^{2}$,
respectively. Next we take the inner product of equation (1) with
$\mathbf{\Phi}$ , multiply equation (2) with $\varphi$ and integrate the
equations over the domain:

$\int_{\Omega}\left\\{\frac{\partial{\mathbf{u}}}{\partial
t}\cdot\mathbf{\Phi}+\left[\left({\mathbf{u}}\cdot\nabla\right)\mathbf{u}\right]\cdot\mathbf{\Phi}-\left[\nabla\cdot\mathbf{T}\right]\cdot\mathbf{\Phi}\right\\}d\Omega=0,$
(68)
$\int_{\Omega}\left\\{\left(\nabla\ldotp\mathbf{u}\right)\varphi\right\\}d\Omega=0.$
(69)

Equations (68) and (69) are the weighted residual forms of the momentum and
continuity equations, respectively. Integrating the last term on the left-
hand-side of (68) by parts,

$\int_{\Omega}\left\\{\left[\nabla\cdot\mathbf{T}\right]\cdot\mathbf{\Phi}\right\\}d\Omega=\int_{\Gamma}\left\\{\mathbf{n}\cdot\mathbf{T}\cdot\mathbf{\Phi}\right\\}d\Gamma_{b}-\int_{\Omega}\left\\{\mathbf{T}:\nabla\left(\mathbf{u}\right)\right\\}d\Omega,$
(70)

where

$\Gamma=\Gamma_{b}+\Gamma_{in}+\Gamma_{wall}+\Gamma_{out}.$

Enforcing the outlet boundary condition and taking into consideration that
$\mathbf{\Phi}\in H^{1}_{0}$, hence $\mathbf{\Phi}$ are zero where essential
boundary conditions are imposed, we are left with

$\Gamma=\Gamma_{b}.$

Making use of (3), equation (68) becomes

$\int_{\Omega}\left\\{\frac{\partial{\mathbf{u}}}{\partial
t}\cdot\mathbf{\Phi}+\left[\left({\mathbf{u}}\cdot\nabla\right)\mathbf{u}\right]\cdot\mathbf{\Phi}+2{Nf}^{-1}\mathbf{E}\left({\mathbf{u}}\right):\mathbf{E}\left(\mathbf{\Phi}\right)-{p}\left(\nabla\ldotp\mathbf{\Phi}\right)\right\\}d\Omega\\\
-\int_{\Gamma_{b}}\left\\{\mathbf{n}\cdot\mathbf{T}\cdot\mathbf{\Phi}\right\\}d\Gamma_{b}=0.$
(71)

The traction term in the last term on the left-hand-side of (71) can be
decomposed into its normal and tangential components (Pozrikidis, 2011):

$\mathbf{n}\cdot{\mathbf{T}}=\left[\mathbf{n}\cdot{\mathbf{T}}\cdot\mathbf{n}\right]\mathbf{n}+\mathbf{n}\times\left[\mathbf{n}\cdot{\mathbf{T}}\times\mathbf{n}\right],$
(72)

thereby allowing the incorporation of normal stress condition (8) and
tangential stress condition (9) into (71) to give

$\displaystyle\int_{\Omega}\left\\{\frac{\partial{\mathbf{u}}}{\partial
t}\cdot\mathbf{\Phi}+\left[\left({\mathbf{u}}\cdot\nabla\right)\mathbf{u}\right]\cdot\mathbf{\Phi}+2{Nf}^{-1}\mathbf{E}\left({\mathbf{u}}\right):\mathbf{E}\left(\mathbf{\Phi}\right)-{p}\left(\nabla\ldotp\mathbf{\Phi}\right)\right\\}d\Omega$
$\displaystyle\quad{}-\int_{\Gamma_{b}}\left\\{\left[{Eo}^{-1}\kappa+\mathrm{z}-\mathrm{P}_{\mathrm{b}}\right]\mathbf{n}\cdot\mathbf{\Phi}\right\\}d\Gamma_{b}=0.$
(73)

Equations (69) and (73) are the weak forms of the governing equations.

##  Appendix B Curvature linearisation

For a three-dimensional axisymmetric surface, the interface location in any
$(r,z)$ plane is sufficient for computing the curvature of the surface.
Consider that the interface in any such plane is spanned by a curve with the
coordinates of any point on the curve being $(r,z)$. In addition, let the
interface be parametrised by length of arc $s$ so that the position vector of
any point on the interface is given as

$\mathbf{r}=r(s)\mathbf{i}_{r}+z(s)\mathbf{i}_{z}.$ (74)

The total curvature at any given point on the interface is defined as

$\kappa=-\nabla_{s}\cdot\mathbf{n},$ (75)

where $r\left(s\right)$ and $z\left(s\right)$ are the radial and axial
coordinates of the points on the interface, respectively; $\mathbf{n}$ remains
the unit normal to the interface and $\nabla_{s}$ is the surface tangential
gradient operator which is given as

$\nabla_{s}=\left(\mathbf{I}-\mathbf{n}\otimes\mathbf{n}\right)\nabla.$ (76)

For a three-dimensional axisymmetric surface, (76) simplifies to

$\nabla_{s}=t_{r}\frac{d}{ds}\mathbf{i}_{r}+\frac{1}{r}\frac{\partial}{\partial\theta}\mathbf{i}_{\theta}+t_{z}\frac{d}{ds}\mathbf{i}_{z},$
(77)

and equation (75) becomes

$\kappa=-\left(\mathbf{t}\cdot\frac{d\mathbf{n}}{ds}+\frac{n_{r}}{r}\right),$
(78)

where $\mathbf{t}$ is the unit tangent vector to the interface with components
$t_{r}$, $t_{\theta}$ and $t_{z}$ in the radial, azimuthal and axial
directions, respectively; $\frac{d}{ds}=\left(\mathbf{t}\cdot\nabla\right)$ is
an operator that denotes the derivative in the tangential direction; $n_{r}$
is the radial component of the unit normal vector, $\mathbf{n}$.

Let us imagine that the deformed interface can be expressed as a summation of
the undeformed interface and a very small deformation. Thus, the deformed
interface can be written as

$\mathbf{r}=\mathbf{r}^{0}+{\mathbf{x}},$ (79)

where
$\mathbf{r}^{0}=r^{0}\mathbf{i}_{r}+\theta^{0}\mathbf{i}_{\theta}+z^{0}\mathbf{i}_{z}$
and $\mathbf{r}=r\mathbf{i}_{r}+\theta\mathbf{i}_{\theta}+z\mathbf{i}_{z}$ are
the undeformed (i.e base) and deformed (i.e perturbed) interface position
vectors, respectively;
$\mathbf{x}=x_{r}\mathbf{i}_{r}+x_{\theta}\mathbf{i}_{\theta}+x_{z}\mathbf{i}_{z}$,
as mentioned in the previous section, is the interface deformation vector and
is taken to be very small in magnitude. Linearisation of the unit normal to,
and elemental arc length on, the deformed interface about the undeformed
interface give (Ramanan & Engelman, 1996; Kruyt et al., 1988; Weatherburn,
1927)

$\mathbf{n}=\mathbf{n}^{0}-\mathbf{n}^{0}\times\nabla_{s}\times\mathbf{x},$
(80) $ds=\left(1+\mathbf{t}^{0}\cdot\frac{d\mathbf{x}}{ds^{0}}\right)ds^{0},$
(81)
$\frac{d}{ds}=\left(1-\mathbf{t}^{0}\cdot\frac{d\mathbf{x}}{ds^{0}}\right)\frac{d}{ds^{0}}.$
(82)

On further simplification of (80)

$\mathbf{n}=\mathbf{n}^{0}-\mathbf{t}^{0}\left(\mathbf{n}^{0}\cdot\frac{d\mathbf{x}}{ds^{0}}\right)-\frac{\mathbf{n}^{0}}{r^{0}}\cdot\frac{\partial\mathbf{x}}{\partial\theta^{0}}\mathbf{i}_{\theta}.$
(83)

Substituting (19) together with (79) and (82) into (78), the linearised
curvature neglecting all terms of nonlinear in $x$ gives

$\kappa=\kappa^{0}+\kappa^{1},$ (84)

with

$\kappa^{0}=-\left(\mathbf{t}^{0}\cdot\frac{d\mathbf{n}^{0}}{ds^{0}}+\frac{n_{r}^{0}}{r^{0}}\right),$
(85)
$\kappa^{1}=\frac{1}{r^{0}}\frac{d}{ds^{0}}\left[r^{0}\left(\mathbf{n}^{0}\cdot\frac{d\mathbf{x}}{ds^{0}}\right)\right]+2\left(\mathbf{t}^{0}\cdot\frac{d\mathbf{x}}{ds^{0}}\right)\left(\mathbf{t}^{0}\cdot\frac{d\mathbf{n}^{0}}{ds^{0}}\right)+\frac{\mathbf{n}^{0}}{{r^{0}}^{2}}\cdot\frac{\partial^{2}\mathbf{x}}{\partial{\theta^{0}}^{2}}+\frac{x_{r}n_{r}^{0}}{{r^{0}}^{2}}\\\
-\frac{d\mathbf{n}^{0}}{ds^{0}}\cdot\frac{d\mathbf{x}}{ds^{0}}.$ (86)

When (20) is further simplified by allowing the deformation vector to be of
the form

$\mathbf{x}=h\mathbf{n}^{0},$ (87)

it results in

$\kappa^{1}=\frac{1}{r^{0}}\frac{d}{ds^{0}}\left(r^{0}\frac{dh}{ds^{0}}\right)+h\left[{\kappa^{0}}^{2}_{a}+{\kappa^{0}}^{2}_{b}+\frac{1}{{r^{0}}^{2}}\frac{\partial^{2}h}{\partial{\theta^{0}}^{2}}\right],$
(88)

which is the same as the expression derived for curvature perturbation in
Chireux et al. (2015), albeit through a different and longer route. In
arriving at (88), we have used the following Frenet-Serret relations

$\displaystyle\frac{d\mathbf{t}}{ds}$ $\displaystyle=\kappa\mathbf{n},$ (89a)
$\displaystyle\frac{d\mathbf{n}}{ds}$ $\displaystyle=-\kappa\mathbf{t}.$ (89b)

Equations (86) and (88) are the expressions for curvature deformation and can
be used for an axisymmetric deformation by setting the term containing
derivative with respect to the azimuthal coordinate to zero to obtain

$\kappa^{1}=\frac{1}{r^{0}}\frac{d}{ds^{0}}\left(r^{0}\frac{dh}{ds^{0}}\right)+h\left[{\kappa^{0}}^{2}_{a}+{\kappa^{0}}^{2}_{b}\right].$
(90)

In equations (80)-(88), $\mathbf{n}^{0}$ is the unit normal vector to the
undeformed interface and $n_{r}^{0}$, $n_{\theta}^{0}$ and $n_{z}^{0}$ are its
component in the radial, azimuthal and axial directions, respectively;
$\mathbf{t}^{0}$ is the unit tangent vector to the undeformed interface with
components $t_{r}^{0}$, $t_{\theta}^{0}$ and $t_{z}^{0}$ in the radial,
azimuthal and axial directions, respectively; $ds$ and $ds^{0}$ are the
elemental arc length for the deformed and undeformed interfaces, respectively;
$\kappa^{0}$ is the curvature of the undeformed surface and $\kappa^{1}$ is
the addition to the undeformed interface curvature (also referred to as
curvature perturbation in the context of linear stability analysis) due to
linearisation of the deformed interface about the undeformed interface; $h$ is
the magnitude of the interface deformation in the direction normal to the
undeformed interface. ${\kappa^{0}}_{a}$ and ${\kappa^{0}}_{b}$ are the two
principal curvatures of the undeformed interface (85) corresponding to the
curvature in the $r-z$ and $r-\theta$ planes, respectively, defined as

$\displaystyle{\kappa^{0}}_{a}$
$\displaystyle=\mathbf{t}^{0}\cdot\frac{d\mathbf{n}^{0}}{ds^{0}},$ (91a)
$\displaystyle{\kappa^{0}}_{b}$ $\displaystyle=\frac{n_{r}^{0}}{r^{0}}.$ (91b)

##  Appendix C Galerkin finite element approximation

In appendix A, we derived the weak forms of the governing equations in the
continuous domain. Let us divide the domain, $\Omega$ into smaller subsets,
finite element, so that each element occupies a sub-domain $\Omega^{e}$ with
boundary $\Gamma^{e}$. The discretised domain is therefore an assemblage of
$n_{e}$ finite elements that make up the domain:

$\int_{\Omega^{h}}\left\\{\cdot\right\\}d\Omega^{h}=\sum_{k=1}^{ne}\int_{\Omega^{e}_{k}}\left\\{\cdot\right\\}d\Omega^{e}_{k},$
(92)
$\int_{\Gamma^{h}}\left\\{\cdot\right\\}d\Gamma^{h}=\sum_{k=1}^{ne}\int_{\Gamma^{e}_{k}}\left\\{\cdot\right\\}d\Gamma^{e}_{k},$
(93)

where $\Omega^{h}$ represents the discretised domain with boundary
$\Gamma^{h}$, and $ne$ is the number of elements in the discretised domain.

Figure 15: Two-dimensional triangulated domain with linear approximations on
the subdomains thereby leading to three nodes.

The geometry of the finite element could be a triangle or quadrilateral in a
two-dimensional domain and in three-dimensional domain, it could be a
tetrahedral or hexahedral. On each element, the unknown variables of the
problem, known as trial functions, are approximated as a linear combination of
unknown parameters and piecewise polynomial functions, basis (or shape )
functions. In addition, the test functions corresponding to each variables are
set to equal to the basis functions used in approximating the variables.
Depending on the geometry of an element and the order of polynomial used in
approximating the unknown variable on the element, $n_{v}$ number of
nodes/unknown parameters can be associated with the element as shown in Figure
15. Thus, on each element, the trial and test functions are approximated as

$\displaystyle\mathrm{u_{r}\left(\mathbf{x},t\right)}$
$\displaystyle=\sum_{k=1}^{n_{v}}\psi^{k}\left(\mathbf{x}\right)v_{r}^{k}\left(t\right)=\mathbf{\Psi}^{T}\mathbf{v}_{r},$
(94a) $\displaystyle\mathrm{u}_{\theta}\left(\mathbf{x},t\right)$
$\displaystyle=\sum_{k=1}^{n_{v}}\psi^{k}\left(\mathbf{x}\right)v_{\theta}^{k}\left(t\right)=\mathbf{\Psi}^{T}\mathbf{v}_{\theta},$
(94b) $\displaystyle\mathrm{u}_{z}\left(\mathbf{x},t\right)$
$\displaystyle=\sum_{k=1}^{n_{v}}\psi^{k}\left(\mathbf{x}\right)v_{z}^{k}\left(t\right)=\mathbf{\Psi}^{T}\mathbf{v}_{z},$
(94c) $\displaystyle\mathrm{p}\left(\mathbf{x},t\right)$
$\displaystyle=\sum_{l=1}^{n_{p}}\lambda^{l}\left(\mathbf{x}\right)\mathfrak{p}^{l}\left(t\right)=\mathbf{\Lambda}^{T}\mathbf{p},$
(94d) $\displaystyle h$
$\displaystyle=\sum_{k=1}^{n_{h}}\eta_{k}\left(\mathbf{r}_{b}\right)\hslash_{k}\left(t\right)=\text{\boldmath${\eta}$}^{T}\mathbf{h},$
(94e)
$\displaystyle\Phi_{r}=\psi^{k},\quad\Phi_{\theta}=\psi^{k},\quad\Phi_{z}=\psi^{k}\quad\forall\quad
k=1,\cdots,n_{v};$ (95a) $\displaystyle\varphi=\lambda^{l}\quad\forall\quad
l=1,\cdots,n_{p};$ (95b) $\displaystyle\xi=\eta_{m}\quad\forall\quad
m=1,\cdots,n_{h}.$ (95c)

In equations (94) and (95a), $\psi$, $\lambda$ and $\eta$ are the shape
functions for the velocity components, pressure, and interface deformation
magnitude, respectively; $n_{v}$, $n_{p}$ and $n_{h}$ are the number of nodes
for the velocity components, pressure, and interface deformation magnitude,
respectively, and $v_{j}\;\forall\;j=r,\;\theta,\;z$ , $\mathfrak{p}$ and
$\hslash$ are the corresponding unknown nodal parameters of the velocity
components, pressure, and interface deformation magnitude.
$\mathbf{v}_{i}\;\forall\;i=r,\;\theta\;z$, $\mathbf{p}$ and $\mathbf{h}$ are
column vectors of the unknown nodal parameters for the velocity components,
pressure, and interface deformation magnitude, respectively; and
$\mathbf{\Psi}$, $\mathbf{\Lambda}$ and ${\eta}$ are the column vectors of the
shape functions for the velocity components, pressure, and interface
deformation magnitude, respectively.

It is seen that the shape functions for the velocity components are the same
but differ from that of the pressure. This is to avoid having an over-
constrained system of discrete equations. In fact, the shape function used for
pressure must be at least one order lower than that of the velocity field
(Reddy & Gartling, 2010).

Using equations (94) and (95a), we can derive, for example, the finite element
model for the continuity equation by writing the weak form of it (37) over an
element and substituting (94) and (95a) to get

$\left[\int_{\Omega^{e}}\left\\{\mathbf{\Lambda}\left(\frac{\partial\mathbf{\Psi}^{T}}{\partial
r}+\frac{\mathbf{\Psi}^{T}}{r}\right)\right\\}d\Omega^{e}\right]\mathbf{v}_{r}+\left[\int_{\Omega^{e}}\left\\{\mathfrak{i}m{\mathbf{\Lambda}}\frac{{\mathbf{\Psi}}^{T}}{r}\right\\}d\Omega^{e}\right]\mathbf{v}_{\theta}\\\
+\left[\int_{\Omega^{e}}\left\\{\mathbf{\Lambda}\frac{\partial\mathbf{\Psi}^{T}}{\partial
z}\right\\}d\Omega^{e}\right]\mathbf{v}_{z}=\mathbf{0},$ (96)

which can be compactly written as

$\mathbf{Q}^{T}\mathbf{v}=\mathbf{0},$ (97)

where $\mathbf{v}=\begin{bmatrix}\mathbf{v}_{r}\\\ \mathbf{v}_{\theta}\\\
\mathbf{v}_{z}\end{bmatrix}$ and $\mathbf{Q}=\begin{bmatrix}\mathbf{Q}_{r}\\\
\mathbf{Q}_{\theta}\\\ \mathbf{Q}_{z}\end{bmatrix}$ .
$\mathbf{Q}_{i}\quad\forall\;i=r,\;\theta\;z$ are the coefficient matrices for
the velocity components. Similar compact relations to (97) can be derived for
the momentum equations and the kinematic boundary condition. The combined
formed of which is given in equation (39) and full expressions for the terms
are:
continuity equation

$\displaystyle{\mathbf{Q}}_{r}=\int_{\Omega^{0,e}}\left\\{\mathbf{\Lambda}\left(\frac{\partial{\mathbf{\Psi}}^{T}}{\partial
r}+\frac{{\mathbf{\Psi}}^{T}}{r}\right)\right\\}d\Omega^{0,e},$ (98)
$\displaystyle{\mathbf{Q}}_{\theta}=\int_{\Omega^{0,e}}\left\\{\mathfrak{i}m{\mathbf{\Lambda}}\frac{{\mathbf{\Psi}}^{T}}{r}\right\\}d\Omega^{0,e},$
(99)
$\displaystyle{\mathbf{Q}}_{z}=\int_{\Omega^{0,e}}\left\\{\mathbf{\Lambda}\frac{\partial{\mathbf{\Psi}}^{T}}{\partial
z}\right\\}d\Omega^{0,e},$ (100)
$\displaystyle\mathbf{Q}_{p}=\int_{\Omega^{0,e}}\left\\{\varepsilon\mathbf{\Lambda}\mathbf{\Lambda}^{T}\right\\}d\Omega^{0,e};$
(101)

momentum equation: growth rate part

$\displaystyle{\mathbf{B}}_{i,i}$
$\displaystyle=\int_{\Omega^{0,e}}\left\\{{\mathbf{\Psi}}{\mathbf{\Psi}}^{T}\right\\}d\Omega^{0,e}\quad\forall\;i=r,\theta,z,$
(102)

momentum equation: convective part

$\displaystyle\mathbf{C}_{i,i}$
$\displaystyle=\int_{\Omega^{0,e}}\left\\{\mathbf{\Psi}\left[\left({\mathbf{u}^{0}}\cdot\nabla\right)\mathbf{\Psi}^{T}+\delta_{r,i}\mathbf{\Psi}^{T}\frac{\partial
u_{r}^{0}}{\partial r}+\delta_{z,i}\mathbf{\Psi}^{T}\frac{\partial
u_{z}^{0}}{\partial z}\right]\right\\}d\Omega^{0,e}\;\forall\;i=r,z,$ (103)
$\displaystyle{\mathbf{C}}_{\theta,\theta}$
$\displaystyle=\int_{\Omega^{0,e}}\left\\{{\mathbf{\Psi}}\left[\left({\mathbf{u}^{0}}\cdot\nabla\right){\mathbf{\Psi}}^{T}+\frac{u_{r}^{0}{\mathbf{\Psi}}^{T}}{r}\right]\right\\}d\Omega^{0,e},$
(104) $\displaystyle\mathbf{C}_{r,z}$
$\displaystyle=\int_{\Omega^{0,e}}\left\\{\mathbf{\Psi}\mathbf{\Psi}^{T}\frac{\partial
u_{r}^{0}}{\partial z}\right\\}d\Omega^{0,e},$ (105)
$\displaystyle\mathbf{C}_{z,r}$
$\displaystyle=\int_{\Omega^{0,e}}\left\\{\mathbf{\Psi}\mathbf{\Psi}^{T}\frac{\partial
u_{z}^{0}}{\partial r}\right\\}d\Omega^{0,e};$ (106)

momentum equation: viscous part

$\displaystyle{\mathbf{K}}_{i,i}$
$\displaystyle=\int_{\Omega^{0,e}}Nf^{-1}\left\\{2\frac{\partial{\mathbf{\Psi}}}{\partial
i}\frac{\partial{\mathbf{\Psi}}^{T}}{\partial
i}+\left(2\delta_{r,i}+m^{2}\right)\frac{{\mathbf{\Psi}}{\mathbf{\Psi}}^{T}}{r^{2}}+\delta_{r,i}\frac{\partial{\mathbf{\Psi}}}{\partial
z}\frac{\partial{\mathbf{\Psi}}^{T}}{\partial
z}+\delta_{z,i}\frac{\partial{\mathbf{\Psi}}}{\partial
r}\frac{\partial{\mathbf{\Psi}}^{T}}{\partial r}\right.$
$\displaystyle\qquad{}+\left.\delta_{z,i}\frac{\partial{\mathbf{\Psi}}}{\partial
r}\frac{\partial{\mathbf{\Psi}}^{T}}{\partial
r}\right\\}d\Omega^{0}\;\forall\;i=r,z,$ (107a)
$\displaystyle{\mathbf{K}}_{\theta,\theta}$
$\displaystyle=\int_{\Omega^{0,e}}Nf^{-1}\left\\{\left(1+2m^{2}\right)\frac{{\mathbf{\Psi}}{\mathbf{\Psi}}^{T}}{r^{2}}+\frac{\partial{\mathbf{\Psi}}}{\partial
z}\frac{\partial{\mathbf{\Psi}}^{T}}{\partial
z}+\frac{\partial{\mathbf{\Psi}}}{\partial
r}\frac{\partial{\mathbf{\Psi}}^{T}}{\partial r}\right.$
$\displaystyle\qquad{}-\left.\left(\frac{{\mathbf{\Psi}}^{T}}{r}\frac{\partial{\mathbf{\Psi}}}{\partial
r}+\frac{{\mathbf{\Psi}}}{r}\frac{\partial{\mathbf{\Psi}}^{T}}{\partial
r}\right)\right\\}d\Omega^{0,e},$ (107b) $\displaystyle\mathbf{K}_{r,z}$
$\displaystyle=\int_{\Omega^{0,e}}Nf^{-1}\left\\{\frac{\partial\mathbf{\Psi}}{\partial
z}\frac{\partial\mathbf{\Psi}^{T}}{\partial r}\right\\}d\Omega^{0,e},$ (107c)
$\displaystyle{\mathbf{K}}_{i,\theta}$
$\displaystyle=\int_{\Omega^{0,e}}Nf^{-1}\left\\{\delta_{r,i}3\mathfrak{i}m\frac{{\mathbf{\Psi}}{\mathbf{\Psi}}^{T}}{r^{2}}-\mathfrak{i}m\frac{{\mathbf{\Psi}}}{r}\frac{\partial{\mathbf{\Psi}}^{T}}{\partial
i}\right\\}d\Omega^{0,e}\;\forall\;i=r,z;$ (107d)

momentum equation: interface deformation part

$\displaystyle{\mathbf{M}}_{i,h}=$
$\displaystyle-\int_{\Gamma_{b}^{0,e}}{Eo}^{-1}\left\\{-\left[n_{i}\frac{d{\mathbf{\Psi}}}{ds}-\kappa_{a}\left(t_{i}{\mathbf{\Psi}}\right)\right]\frac{d\text{\boldmath${{\eta}}$}^{T}}{ds}+n_{i}\left[\kappa_{a}^{2}+\kappa_{b}^{2}-\frac{m^{2}}{r^{2}}\right]{\mathbf{\Psi}}\text{\boldmath${{\eta}}$}^{T}\right\\}d\Gamma_{b}^{0,e}$
$\displaystyle-\int_{\Gamma_{b}^{0}}\left\\{n_{z}n_{i}{\mathbf{\Psi}}\text{\boldmath${\eta}$}^{T}\right\\}d\Gamma_{b}^{0}$
$\displaystyle+\int_{\Gamma^{0}_{b}}\left\\{\left\\{{\left({u}^{0}_{i}{t}_{i}\right)\left[\frac{d{u}^{0}_{i}}{ds}{\mathbf{\Psi}}\right]}+\left[-p^{0}+2Nf^{-1}\left({t}_{i}\frac{d{u}^{0}_{i}}{ds}\right)\right]\left({t}_{i}\frac{d{\mathbf{\Psi}}}{ds}\right)\right\\}\text{\boldmath${\eta}$}^{T}\right\\}d\Gamma^{0}_{b}$
$\displaystyle+\int_{\Gamma_{b}^{0}}\left\\{\left[{Eo}^{-1}\kappa+z-P_{b}^{0}\right]\left[\left({t}_{i}{\mathbf{\Psi}}\right)\frac{d\text{\boldmath${\eta}$}^{T}}{ds}+\kappa\left({n}_{i}{\mathbf{\Psi}}\text{\boldmath${\eta}$}^{T}\right)\right]\right\\}d\Gamma_{b}^{0}\quad\forall\;i=r,z,$
(108a) $\displaystyle{\mathbf{M}}_{\theta,h}=$
$\displaystyle\int_{\Gamma^{0}_{b}}\left\\{\left[\left(-p^{0}+2Nf^{-1}\frac{u_{r}^{0}}{r}\right)-\left({Eo}^{-1}\kappa+z-P_{b}^{0}\right)\right]\left(-\mathfrak{i}m\frac{{\mathbf{\Psi}}\text{\boldmath${\eta}$}^{T}}{r}\right)\right\\}d\Gamma^{0}_{b};$
(108b)

kinematic boundary condition

$\displaystyle{\mathbf{X}}_{i}$
$\displaystyle=\int_{\Gamma_{b}^{0,e}}n_{i}\left\\{\text{\boldmath${\eta}$}\mathbf{\Psi}^{T}\right\\}d\Gamma_{b}^{0,e}\quad\;\forall\;i=r,z,$
(109a) $\displaystyle{\mathbf{X}}_{h}$
$\displaystyle=\int_{\Gamma_{b}^{0,e}}\left\\{\left(\mathbf{u}^{0}\cdot\mathbf{t}\right)\text{\boldmath${\eta}$}\frac{d\text{\boldmath${\eta}$}^{T}}{ds}-\left(\frac{d\mathbf{u}^{0}}{dn}\cdot\mathbf{n}\right)\text{\boldmath${\eta}$}\text{\boldmath${\eta}$}^{T}\right\\}d\Gamma_{b}^{0,e},$
(109b) $\displaystyle{\mathbf{B}}_{h}$
$\displaystyle=\int_{\Gamma_{b}^{0,e}}\left\\{\text{\boldmath${\eta}$}\text{\boldmath${\eta}$}^{T}\right\\}d\Gamma_{b}^{0,e}.$
(109c)

## References

  * Abubakar & Matar (2021) Abubakar, H. A. & Matar, O. K. 2021 Taylor bubble motion in stagnant and flowing liquids in vertical pipes. part i: Steady-states. Submitted to J. Fluid Mech. .
  * Anjos et al. (2014) Anjos, G., Mangiavacchi, N., Borhani, N. & Thome, J. R. 2014 3D ALE finite-element method for two-phase flows with phase change. Heat Transfer Engineering 35 (5), 537–547.
  * Boomkamp & Miesen (1996) Boomkamp, P.A.M. & Miesen, R.H.M. 1996 Classification of instabilities in parallel two-phase flow. Int. J. Multiphase Flow 22, 67–88.
  * Brown (1965) Brown, R.A.S. 1965 The mechanics of large gas bubbles in tubes I. Bubble velocities in stagnant liquids. Can. J. Chem. Eng 43, 217–223.
  * Bugg & Saad (2002) Bugg, J. D. & Saad, G. A. 2002 The velocity field around a Taylor bubble rising in a stagnant viscous fluid: Numerical and experimental results. Int. J. Multiphase Flow 28, 791–803.
  * Cairncross et al. (2000) Cairncross, R.A., Schunk, P.R., Baer, T.A., Rao, R.R. & Sackinger, P.A. 2000 A finite element method for free surface flows of incompressible fluids in three dimensions. Part I. Boundary fitted mesh motion. Int. J. Numer. Meth. Fluids 33, 375–403.
  * Campos & Guedes de Carvalho (1988) Campos, J.B.L.M. & Guedes de Carvalho, J.R.F. 1988 An experimental study of the wake of gas slugs rising in liquids. J. Fluid Mech. 196, 27–37.
  * Carvalho & Scriven (1999) Carvalho, M.S. & Scriven, L.E. 1999 Three-dimensional stability analysis of free surface flows: application to forward deformable roll coating. J. Comput. Phys. 151, 534–562.
  * Chireux et al. (2015) Chireux, V., Fabre, D., Risso, F. & Tordjeman, P. 2015 Oscillations of a liquid bridge resulting from the coalescence of two droplets. Physics of Fluids 27, 062103.
  * Christodoulou & Scriven (1988) Christodoulou, C.N. & Scriven, L.E. 1988 Finding leading modes of a viscous free surface flow: An asymmetric generalized eigenproblem. J. Sci. Comput. 3 (4), 355–406.
  * Collins et al. (1978) Collins, R., De Moraes, F., Davidson, J. & Harrison, D. 1978 The motion of a large gas bubble rising through liquid flowing in a tube. J. Fluid Mech 89, 497–514.
  * Cuvelier & Schulkes (1990) Cuvelier, C. & Schulkes, R.M.S.M. 1990 Some numerical methods for the computation of capillary free boundaries governed by the Navier-Stokes equations. SIAM Rev. 32 (3), 355–423.
  * Dumitrescu (1943) Dumitrescu, D.T. 1943 Strömung an einer Luftblase im senkrechten Rohr. Z. Angew. Math. Mech 23 (3), 139–149.
  * Fabre (2016) Fabre, J. 2016 A long bubble rising in still liquid in a vertical channel: a plane inviscid solution. J. Fluid Mech. 794, R4.
  * Fabre & Figueroa-Espinoza (2014) Fabre, J. & Figueroa-Espinoza, B. 2014 Taylor bubble rising in a vertical pipe against laminar or turbulent downward flow: symmetric to asymmetric shape transition. J. Fluid Mech. 755, 485–502.
  * Fershtman et al. (2017) Fershtman, A., Babin, V., Barnea, D. & Shemer, L. 2017 On shapes and motion of an elongated bubble in downward liquid pipe flow. Physics of Fluids 29, 112103.
  * Figueroa-Espinoza & Fabre (2011) Figueroa-Espinoza, B. & Fabre, J. 2011 Taylor bubble moving in a flowing liquid in vertical channel: transition from symmetric to asymmetric shape. J. Fluid Mech. 679, 432–454.
  * Funada et al. (2005) Funada, T., Joseph, D., Maehara, T. & Yamashita, S. 2005 Ellipsoidal model of the rise of a Taylor bubble in a round tube. Int. J. Multiph. Flow 31, 473–491.
  * Griffith & Wallis (1961) Griffith, P. & Wallis, G. B. 1961 Two phase slug flow. ASME: J. Heat Transfer 83, 07–320.
  * Ha Ngoc & Fabre (2006) Ha Ngoc, H. & Fabre, J. 2006 A boundary element method for calculating the shape and velocity of two-dimensional long bubble in stagnant and flowing liquid. Engng Anal. Bound. Elem. 30, 539–552.
  * Heinrich & Pepper (1999) Heinrich, J.C. & Pepper, D.W. 1999 Intermediate finite element method : fluid flow and heat transfer applications. Taylor & Francis Group.
  * Hooper & Boyd (1983) Hooper, A.P. & Boyd, W.G.C. 1983 Shear-flow instability at the interface between two viscous fluids. J. Fluid Mech. 128, 507–528.
  * Hu & Joseph (1989) Hu, H.H & Joseph, D.D. 1989 Lubricated pipelining: stability of core-annular flow. Part 2. J. Fluid Mech. 205, 359–396.
  * Kang et al. (2010) Kang, C.W., Quan, S.P. & Lou, J. 2010 Numerical study of a Taylor bubble rising in stagnant liquids. Phys. Rev. E. 81, 1539–3755.
  * Kruyt et al. (1988) Kruyt, N.P., Cuvelier, C., Segal, A. & Van Der Zanden, J. 1988 A total linearization method for solving viscous free boundary flow problems by the finite-element method. Int. J. Numer. Methods Fluids 8, 351–363.
  * Lehoucq et al. (1997) Lehoucq, R.B., Sorensen, D.C. & Yang, C. 1997 ARPACK user’s guide: solution of large scale eigenvalue problems with implicitly restarted Arnoldi methods. SIAM .
  * Lizarraga-Garcia et al. (2017) Lizarraga-Garcia, E., Buongiorno, J., Al-Safran, E. & Lakehal, D. 2017 A broadly applicable unified closure relation for Taylor bubble rise velocity in pipes with stagnant liquid. Int. J. Multiph. Flow 89, 345–358.
  * Llewellin et al. (2012) Llewellin, E.W., Del Bello, E., Taddeucci, J., Scarlato, P. & Lane, S.J. 2012 The thickness of the falling film of liquid around a Taylor bubble. Proc. R. Soc. A 468, 1041–1064.
  * Lu & Prosperetti (2006) Lu, X & Prosperetti, A. 2006 Axial stability of Taylor bubbles. J. Fluid Mech. 568, 173–192.
  * Lu & Prosperetti (2009) Lu, X & Prosperetti, A. 2009 A numerical study of Taylor bubbles. Ind. Eng. Chem. Res. 48, 242–252.
  * Mao & Dukler (1990) Mao, Z.S. & Dukler, A.E. 1990 The motion of Taylor bubbles in vertical tubes I. A numerical simulation for the shape and rise velocity of Taylor bubbles in stagnant and flowing liquids. J. Comput. Phys. 91, 2055–2064.
  * Mao & Dukler (1991) Mao, Z.S. & Dukler, A.E. 1991 The motion of Taylor bubbles in vertical tubes II. Experimental data and simulations for laminar and turbulent flow. Chem. Eng. Sci. 46, 132–160.
  * Martin (1976) Martin, C.S. 1976 Vertically downward two-phase slug flow. Trans. ASME J. Fluids Engng 98, 715–722.
  * Miller & Scriven (1968) Miller, C.A. & Scriven, L.E. 1968 The oscillations of a fluid droplet immersed in another fluid. J . Fluid Mech. 32, 417–435.
  * Moissis & Griffith (1962) Moissis, R. & Griffith, P. 1962 Entrance effect in a two-phase slug flow. J. Heat Transf. 84, 29–38.
  * Natarajan (1992) Natarajan, R. 1992 An Arnoldi-based iterative scheme for nonsymmetric matrix pencils arising in finite element stability problems. J. Comput. Phys. 100, 128–142.
  * Nicklin et al. (1962) Nicklin, D., Wilkes, J. & Davidson, J. 1962 Two-phase flow in vertical tubes. Trans. Inst. Chem. Engrs 40, 61–68.
  * Nigmatulin (2001) Nigmatulin, T.R. 2001 Surface of a Taylor bubble in vertical cylindrical flows. Dokl. Phys. 46, 803–805.
  * Nogueira et al. (2006) Nogueira, S., Riethmuller, M.L., Campos, J.B.L.M. & Pinto, A.M.F.R. 2006 Flow patterns in the wake of a Taylor bubble rising through vertical columns of stagnant and flowing Newtonian liquids: an experimental study. Chem. Eng. Sci. 61, 7199–7212.
  * Ó Náraigh et al. (2011) Ó Náraigh, L., Spelt, P.D.M., Matar, O. K. & Zaki, T.A. 2011 Interfacial instability in turbulent flow over a liquid film in a channel. Int. J. Multiphase Flow 37, 812–830.
  * Polonsky et al. (1999) Polonsky, S., Shemer, L. & Barnea, D. 1999 The relation between the Taylor bubble motion and the velocity field ahead of it. Int. J. Multiphase Flow 25, 957–975.
  * Pozrikidis (2011) Pozrikidis, C. 2011 Introduction to theoretical and computational fluid dynamics. Taylor & Francis Group.
  * Pringle et al. (2015) Pringle, C.C.T., Ambrose, S., Azzopardi, B.J. & Rust, A.C. 2015 The existence and behaviour of large diameter Taylor bubbles. Int. J. Multiphase Flow 72, 318–323.
  * Prosperetti (1980) Prosperetti, A. 1980 Normal-mode analysis for the oscillations of a viscous liquid drop in an immiscible liquid. Journal de Mécanique 19 (1), 149–181.
  * Ramanan & Engelman (1996) Ramanan, N. & Engelman, M.S. 1996 An algorithm for simulation of steady free surface flows. Int. J. Numer. Methods Fluids 22, 103–120.
  * Rana et al. (2015) Rana, B.K., Das, A.K. & Das, P.K. 2015 Mechanism of bursting Taylor bubbles at free surfaces. Langmuir 31, 9870–9881.
  * Reddy & Gartling (2010) Reddy, J.N. & Gartling, D.K. 2010 The finite element method in heat transfer and fluid dynamics.. Taylor & Francis Group.
  * Sahu et al. (2009) Sahu, K.C., Ding, H., Valluri, P. & Matar, O.K. 2009 Linear stability analysis and numerical simulation of miscible two-layer channel flow. Physics of Fluids 21, 042104.
  * Sahu et al. (2007) Sahu, K.C., Valluri, P., Spelt, P.D.M. & Matar, O.K. 2007 Linear instability of pressure-driven channel flow of a newtonian and Herschel-Bulkley fluid. Physics of Fluids 19, 122101\.
  * Selvam et al. (2007) Selvam, B., Merk, S., Govindarajan, R. & Meiburg, E. 2007 Stability of miscible core-annular flows with viscosity stratification. J. Fluid Mech. 592, 23–49.
  * Taha & Cui (2006) Taha, T. & Cui, Z.F. 2006 CFD modeling of slug flow in vertical tubes. Chem. Eng. Sci. 61, 676–687.
  * Weatherburn (1927) Weatherburn, C.E. 1927 Differential geometry of three dimensions. UK: Cambridge University press, part 127, pp 251.
  * White & Beardmore (1962) White, E.T. & Beardmore, R.H. 1962 The velocity of rise of single cylindrical air bubbles through liquids contained in vertical tubes. Chem. Eng. Sci. 17, 351–361.