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

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.

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.

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.