2.1 Equations of motion

The governing equations for the fluid velocity, $\boldsymbol{u}_{j}$ , pressure, $p_{j}$ , density, $\unicode[STIX]{x1D70C}_{j}$ , and viscosity, $\unicode[STIX]{x1D707}_{j}$ , in each fluid layer are the conservation of mass equation, an equation of motion and an incompressibility condition given respectively by

(2.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{j}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{j}\boldsymbol{u}_{j})=0, & \displaystyle\end{eqnarray}$$

(2.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{j}\frac{\text{D}\boldsymbol{u}_{j}}{\text{D}t}=-\unicode[STIX]{x1D735}p_{j}-\unicode[STIX]{x1D70C}_{j}\unicode[STIX]{x1D734}\times (\unicode[STIX]{x1D734}\times \boldsymbol{x})-2\unicode[STIX]{x1D70C}_{j}\unicode[STIX]{x1D734}\times \boldsymbol{u}_{j}+\unicode[STIX]{x1D707}_{j}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{j}+2\unicode[STIX]{x1D65A}_{j}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D707}_{j}+\unicode[STIX]{x1D707}_{j}\unicode[STIX]{x1D735}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{j}), & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(2.1c ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}\unicode[STIX]{x1D70C}_{j}}{\text{D}t}=0, & \displaystyle\end{eqnarray}$$

$j=1,2$

$\unicode[STIX]{x1D65A}_{j}=(\unicode[STIX]{x1D735}\boldsymbol{u}_{j}+\unicode[STIX]{x1D735}\boldsymbol{u}_{j}^{\mathsf{T}})/2$

$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{j}=0$

$\unicode[STIX]{x1D6FA}$

$a$

$\unicode[STIX]{x1D70C}_{0}=(\unicode[STIX]{x1D70C}_{1}+\unicode[STIX]{x1D70C}_{2})/2$

$\unicode[STIX]{x1D707}_{0}=(\unicode[STIX]{x1D707}_{1}+\unicode[STIX]{x1D707}_{2})/2$

(2.2a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{j}^{\prime }}{\unicode[STIX]{x2202}t^{\prime }}+\unicode[STIX]{x1D735}^{\prime }\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{j}^{\prime }\boldsymbol{u}_{j}^{\prime })=0, & \displaystyle\end{eqnarray}$$

(2.2b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}\boldsymbol{u}_{j}^{\prime }}{\text{D}t^{\prime }}=-\frac{1}{\unicode[STIX]{x1D70C}_{j}^{\prime }}\unicode[STIX]{x1D735}^{\prime }p_{j}^{\prime }+r^{\prime }\hat{\boldsymbol{r}}-2\hat{\boldsymbol{z}}\times \boldsymbol{u}_{j}^{\prime }+\frac{1}{Re}\frac{\unicode[STIX]{x1D707}_{j}^{\prime }}{\unicode[STIX]{x1D70C}_{j}^{\prime }}\unicode[STIX]{x1D6FB}^{\prime 2}\boldsymbol{u}_{j}^{\prime }+\frac{1}{Re}\frac{2}{\unicode[STIX]{x1D70C}_{j}^{\prime }}\unicode[STIX]{x1D65A}_{j}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735}^{\prime }\unicode[STIX]{x1D707}_{j}^{\prime }, & \displaystyle\end{eqnarray}$$

(2.2c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}^{\prime }\boldsymbol{\cdot }\boldsymbol{u}_{j}^{\prime }=0, & \displaystyle\end{eqnarray}$$

$Re=\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D6FA}a^{2}/\unicode[STIX]{x1D707}_{0}$

$\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D6FA}^{2}a^{2}$

$\hat{\boldsymbol{r}}$

$\hat{\boldsymbol{z}}$

$Ek=\unicode[STIX]{x1D707}_{0}/(\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D6FA}a^{2})$

Hydrostatic solutions to the governing system (2.2) are denoted by a superscript ‘ $\ast$ ’ and are given by

(2.3a-c ) $$\begin{eqnarray}\displaystyle \boldsymbol{u}_{j}^{\ast }=\mathbf{0},\quad \unicode[STIX]{x1D70C}_{j}^{\ast }=\unicode[STIX]{x1D70C}_{j}^{\ast }(r),\quad p_{j}^{\ast }=p_{0}+\int _{0}^{r}\unicode[STIX]{x1D70C}_{j}^{\ast }(\unicode[STIX]{x1D709})\unicode[STIX]{x1D709}\,\text{d}\unicode[STIX]{x1D709}, & & \displaystyle\end{eqnarray}$$

for a reference pressure $p_{0}$ at the origin and a radially varying initial density distribution $\unicode[STIX]{x1D70C}_{j}^{\ast }=\unicode[STIX]{x1D70C}_{j}^{\ast }(r)$ . In the hydrostatic configuration the form of $\unicode[STIX]{x1D707}_{j}$ is arbitrary.

We consider linear perturbations to the hydrostatic initial condition of the form

(2.4a,b ) $$\begin{eqnarray}\displaystyle \boldsymbol{u}_{j}=\boldsymbol{u}_{j}^{\ast }+\unicode[STIX]{x1D716}\boldsymbol{U}_{j}(r,\unicode[STIX]{x1D703},t),\quad \unicode[STIX]{x1D70C}_{j}=\unicode[STIX]{x1D70C}_{j}^{\ast }(r)+\unicode[STIX]{x1D716}\unicode[STIX]{x1D70E}_{j}(r)\text{e}^{\text{i}(m\unicode[STIX]{x1D703}+\unicode[STIX]{x1D714}t)}, & & \displaystyle\end{eqnarray}$$

(2.4c,d ) $$\begin{eqnarray}\displaystyle p_{j}=p_{j}^{\ast }+\unicode[STIX]{x1D716}P_{j}(r)\text{e}^{\text{i}(m\unicode[STIX]{x1D703}+\unicode[STIX]{x1D714}t)},\quad \unicode[STIX]{x1D707}_{j}=\unicode[STIX]{x1D707}_{j}^{\ast }(r)+\unicode[STIX]{x1D716}\unicode[STIX]{x1D702}_{j}(r)\text{e}^{\text{i}(m\unicode[STIX]{x1D703}+\unicode[STIX]{x1D714}t)}, & & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D716}\ll 1$

$m\in \mathbb{N}$

$\unicode[STIX]{x1D713}$

(2.5) $$\begin{eqnarray}\displaystyle \boldsymbol{U}_{j}=\left(\frac{1}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{j}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}},-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{j}}{\unicode[STIX]{x2202}r},0\right), & & \displaystyle\end{eqnarray}$$

and (2.2c ) is automatically satisfied in each layer. Writing the streamfunction as $\unicode[STIX]{x1D713}_{j}=\unicode[STIX]{x1D719}_{j}(r)\text{e}^{\text{i}(m\unicode[STIX]{x1D703}+\unicode[STIX]{x1D714}t)}$ we have

(2.6) $$\begin{eqnarray}\displaystyle \boldsymbol{u}_{j}=\unicode[STIX]{x1D716}\left(\frac{\text{i}m\unicode[STIX]{x1D719}_{j}}{r},-\frac{\text{d}\unicode[STIX]{x1D719}_{j}}{\text{d}r},0\right)\text{e}^{\text{i}(m\unicode[STIX]{x1D703}+\unicode[STIX]{x1D714}t)}. & & \displaystyle\end{eqnarray}$$

The solution growth rate is controlled by $\text{Im}(\unicode[STIX]{x1D714})$ , its precession is controlled by $\text{Re}(\unicode[STIX]{x1D714})$ and the azimuthal wavenumber of the perturbation to the interface is $m$ . The linearised mass conservation (2.2a ) equation gives

(2.7a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}_{j}=-\frac{m\unicode[STIX]{x1D719}_{j}}{\unicode[STIX]{x1D714}r}\frac{\text{d}\unicode[STIX]{x1D70C}_{j}^{\ast }}{\text{d}r}. & & \displaystyle\end{eqnarray}$$

The linearised radial equation of motion gives

(2.7b ) $$\begin{eqnarray}\displaystyle & & \displaystyle -\frac{m\unicode[STIX]{x1D714}\unicode[STIX]{x1D719}_{j}}{r}=-\frac{1}{\unicode[STIX]{x1D70C}_{j}^{\ast }}\frac{\text{d}P_{j}}{\text{d}r}+\frac{\unicode[STIX]{x1D70E}_{j}r}{\unicode[STIX]{x1D70C}_{j}^{\ast }}-2\frac{\text{d}\unicode[STIX]{x1D719}_{j}}{\text{d}r}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\frac{\text{i}m}{Re}\frac{\unicode[STIX]{x1D707}_{j}^{\ast }}{\unicode[STIX]{x1D70C}_{j}^{\ast }}\left\{\frac{1}{r}\frac{\text{d}}{\text{d}r}\left[r\frac{\text{d}}{\text{d}r}\left(\frac{\unicode[STIX]{x1D719}_{j}}{r}\right)\right]-\frac{(1+m^{2})\unicode[STIX]{x1D719}_{j}}{r^{3}}+\frac{2}{r^{2}}\frac{\text{d}\unicode[STIX]{x1D719}_{j}}{\text{d}r}+\frac{2}{\unicode[STIX]{x1D707}_{j}^{\ast }}\frac{\text{d}\unicode[STIX]{x1D707}_{j}^{\ast }}{\text{d}r}\left[\frac{1}{r}\frac{\text{d}\unicode[STIX]{x1D719}_{j}}{\text{d}r}-\frac{\unicode[STIX]{x1D719}_{j}}{r^{2}}\right]\right\},\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

and the linearised azimuthal equation of motion is

(2.7c ) $$\begin{eqnarray}\displaystyle -\text{i}\unicode[STIX]{x1D714}\frac{\text{d}\unicode[STIX]{x1D719}_{j}}{\text{d}r} & = & \displaystyle -\frac{\text{i}mP_{j}}{\unicode[STIX]{x1D70C}_{j}^{\ast }r}-2\frac{\text{i}m\unicode[STIX]{x1D719}_{j}}{r}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{Re}\frac{\unicode[STIX]{x1D707}_{j}^{\ast }}{\unicode[STIX]{x1D70C}_{j}^{\ast }}\left\{-\frac{1}{r}\frac{\text{d}}{\text{d}r}\left(r\frac{\text{d}^{2}\unicode[STIX]{x1D719}_{j}}{\text{d}r^{2}}\right)+\frac{(1+m^{2})}{r^{2}}\frac{\text{d}\unicode[STIX]{x1D719}_{j}}{\text{d}r}\right.\nonumber\\ \displaystyle & & \displaystyle -\left.\frac{2m^{2}\unicode[STIX]{x1D719}_{j}}{r^{3}}-\frac{1}{\unicode[STIX]{x1D707}_{j}^{\ast }}\frac{\text{d}\unicode[STIX]{x1D707}_{j}^{\ast }}{\text{d}r}\left[\frac{\text{d}^{2}\unicode[STIX]{x1D719}_{j}}{\text{d}r^{2}}-\frac{1}{r}\frac{\text{d}\unicode[STIX]{x1D719}_{j}}{\text{d}r}+\frac{m^{2}\unicode[STIX]{x1D719}_{j}}{r^{2}}\right]\right\}.\end{eqnarray}$$

Eliminating $P_{j}$ and $\unicode[STIX]{x1D70E}_{j}$ from (2.7b ) using (2.7a ) and (2.7c ) yields a fourth order, one-dimensional linear ordinary differential equation, of an Orr–Sommerfeld type, for perturbations to the system given by

(2.8) $$\begin{eqnarray}\displaystyle & & \displaystyle \text{i}\unicode[STIX]{x1D714}\left\{\left(\unicode[STIX]{x1D719}_{j}^{\prime \prime }+\frac{\unicode[STIX]{x1D719}_{j}^{\prime }}{r}-\frac{m^{2}\unicode[STIX]{x1D719}_{j}}{r^{2}}\right)+\frac{\unicode[STIX]{x1D70C}_{j}^{\ast \prime }}{\unicode[STIX]{x1D70C}_{j}^{\ast }}\left(\unicode[STIX]{x1D719}_{j}^{\prime }+\frac{m(m-2\unicode[STIX]{x1D714})\unicode[STIX]{x1D719}_{j}}{\unicode[STIX]{x1D714}^{2}r}\right)\right\}\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{1}{Re}\frac{\unicode[STIX]{x1D707}_{j}^{\ast }}{\unicode[STIX]{x1D70C}_{j}^{\ast }}\left\{\unicode[STIX]{x1D719}_{j}^{\prime \prime \prime \prime }+\frac{2\unicode[STIX]{x1D719}_{j}^{\prime \prime \prime }}{r}-(1+2m^{2})\left[\frac{\unicode[STIX]{x1D719}_{j}^{\prime \prime }}{r^{2}}-\frac{\unicode[STIX]{x1D719}_{j}^{\prime }}{r^{3}}\right]+\frac{m^{2}(m^{2}-4)\unicode[STIX]{x1D719}_{j}}{r^{4}}\right.\nonumber\\ \displaystyle & & \displaystyle \qquad \left.+\,\frac{\unicode[STIX]{x1D707}_{j}^{\ast \prime }}{\unicode[STIX]{x1D707}_{j}^{\ast }}\left(2\unicode[STIX]{x1D719}_{j}^{\prime \prime \prime }+\frac{\unicode[STIX]{x1D719}_{j}^{\prime \prime }}{r}-\frac{(1+2m^{2})\unicode[STIX]{x1D719}_{j}^{\prime }}{r^{2}}+\frac{3m^{2}\unicode[STIX]{x1D719}_{j}}{r^{3}}\right)+\frac{\unicode[STIX]{x1D707}_{j}^{\ast \prime \prime }}{\unicode[STIX]{x1D707}_{j}^{\ast }}\left(\unicode[STIX]{x1D719}_{j}^{\prime \prime }-\frac{\unicode[STIX]{x1D719}_{j}^{\prime }}{r}+\frac{m^{2}\unicode[STIX]{x1D719}_{j}}{r^{2}}\right)\right\},\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where a prime symbol now indicates differentiation with respect to $r$ .

2.2 Boundary conditions

In both the inviscid and viscous cases we require finite velocities on the axis $r=0$ and a no-penetration condition (no-normal-velocity condition) on $r=1$ such that

(2.9a,b ) $$\begin{eqnarray}\displaystyle |\boldsymbol{u}(r=0)|<\infty ,\quad \boldsymbol{u}(r=1)\boldsymbol{\cdot }\hat{\boldsymbol{r}}=0. & & \displaystyle\end{eqnarray}$$

In the case of a viscous outer fluid we will also require a no-slip condition on $r=1$ , viz.

(2.10) $$\begin{eqnarray}\displaystyle \boldsymbol{u}(r=1)\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D73D}}=0. & & \displaystyle\end{eqnarray}$$

These three conditions become, in terms of (2.6)

(2.11a-d ) $$\begin{eqnarray}\displaystyle \lim _{r\rightarrow 0}\left|\frac{\unicode[STIX]{x1D719}_{1}}{r}\right|<\infty ,\quad |\unicode[STIX]{x1D719}_{1}^{\prime }(0)|<\infty ,\quad \unicode[STIX]{x1D719}_{2}(1)=0,\quad \unicode[STIX]{x1D719}_{2}^{\prime }(1)=0, & & \displaystyle\end{eqnarray}$$

where the first two conditions are the velocity regularity condition, the third is the no-penetration condition and the final condition is the no-slip condition that applies when the outer fluid is viscous.

The condition of stress continuity at the interface of two fluids with surface tension is given dimensionally by

(2.12) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E5}\{\unicode[STIX]{x1D748}\boldsymbol{\cdot }\hat{\boldsymbol{n}}\}=\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{\boldsymbol{n}})\hat{\boldsymbol{n}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6E5}\{\cdot \}$ indicates the jump in a quantity from the outer fluid 2 to the inner fluid 1 across the interface $\mathscr{S}$ , $\unicode[STIX]{x1D748}$ is the stress tensor, $\hat{\boldsymbol{n}}$ is a unit normal vector at the interface directed from fluid 1 into fluid 2, and $\unicode[STIX]{x1D6FE}$ is the coefficient of surface tension. Non-dimensionalising (and dropping the prime notation immediately) we have

(2.13) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E5}\{\unicode[STIX]{x1D748}\boldsymbol{\cdot }\hat{\boldsymbol{n}}\}=\frac{1}{We}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{\boldsymbol{n}})\hat{\boldsymbol{n}}, & & \displaystyle\end{eqnarray}$$

where $We=\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D6FA}^{2}a^{3}/\unicode[STIX]{x1D6FE}$ is a non-dimensional Weber number representing the ratio of inertial to curvature effects.

Taking the interface, $\mathscr{S}$ , to be defined by

(2.14) $$\begin{eqnarray}\displaystyle \mathscr{S}:=r-(r_{0}+\unicode[STIX]{x1D716}\text{e}^{\text{i}(m\unicode[STIX]{x1D703}+\unicode[STIX]{x1D714}t)})=0, & & \displaystyle\end{eqnarray}$$

then the unit vector pointing from fluid 1 into fluid 2 is given by

(2.15) $$\begin{eqnarray}\displaystyle \hat{\boldsymbol{n}}=\frac{\unicode[STIX]{x1D735}\mathscr{S}}{|\unicode[STIX]{x1D735}\mathscr{S}|}=(1+O(\unicode[STIX]{x1D716}^{2}))\hat{\boldsymbol{r}}+\left(-\frac{\text{i}m}{r}\unicode[STIX]{x1D716}\text{e}^{\text{i}(m\unicode[STIX]{x1D703}+\unicode[STIX]{x1D714}t)}+O(\unicode[STIX]{x1D716}^{2})\right)\hat{\unicode[STIX]{x1D73D}}. & & \displaystyle\end{eqnarray}$$

Hence, on $\mathscr{S}$

(2.16) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{\boldsymbol{n}}=(1+O(\unicode[STIX]{x1D716}^{2}))\hat{\boldsymbol{r}}+\left(-\frac{\text{i}m}{r_{0}}\unicode[STIX]{x1D716}\text{e}^{\text{i}(m\unicode[STIX]{x1D703}+\unicode[STIX]{x1D714}t)}+O(\unicode[STIX]{x1D716}^{2})\right)\hat{\unicode[STIX]{x1D73D}}, & \displaystyle\end{eqnarray}$$

(2.17) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{\boldsymbol{n}}=\frac{1}{r_{0}}+\unicode[STIX]{x1D716}\frac{m^{2}-1}{r_{0}^{2}}\text{e}^{\text{i}(m\unicode[STIX]{x1D703}+\unicode[STIX]{x1D714}t)}+O(\unicode[STIX]{x1D716}^{2}). & \displaystyle\end{eqnarray}$$

The stress tensor in each fluid layer is given in non-dimensional terms by

(2.18) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D748}_{j}=-p_{j}\unicode[STIX]{x1D644}+\frac{2\unicode[STIX]{x1D707}_{j}}{Re}\unicode[STIX]{x1D65A}_{j}\quad \Rightarrow \quad \unicode[STIX]{x1D6E5}\{\unicode[STIX]{x1D748}\boldsymbol{\cdot }\hat{\boldsymbol{n}}\}=\unicode[STIX]{x1D6E5}\{-p\hat{\boldsymbol{n}}\}+\frac{2}{Re}\unicode[STIX]{x1D6E5}\{\unicode[STIX]{x1D707}\unicode[STIX]{x1D65A}\boldsymbol{\cdot }\hat{\boldsymbol{n}}\}. & & \displaystyle\end{eqnarray}$$

We consider the two terms on the right-hand side separately. Taylor expanding about $r=r_{0}$ we have that

(2.19) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E5}\{-p\hat{\boldsymbol{n}}\}=-\left[p^{\ast }+\unicode[STIX]{x1D716}\left(\left.\frac{\text{d}p^{\ast }}{\text{d}r}\right|_{r=r_{0}}+P\right)\text{e}^{\text{i}(m\unicode[STIX]{x1D703}+\unicode[STIX]{x1D714}t)}+O(\unicode[STIX]{x1D716}^{2})\right]_{-}^{+}\hat{\boldsymbol{n}}, & & \displaystyle\end{eqnarray}$$

where the jump on the right-hand side is across $r=r_{0}$ (as distinct from the jump across the interface). Substituting in both the hydrostatic condition $\text{d}p_{j}^{\ast }/\text{d}r=\unicode[STIX]{x1D70C}_{j}^{\ast }r$ and the expression for $\hat{\boldsymbol{n}}$ on $\mathscr{S}$ , equation (2.16), we may rewrite this last expression as

(2.20) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E5}\{-p\hat{\boldsymbol{n}}\} & = & \displaystyle -\{[p^{\ast }+\unicode[STIX]{x1D716}(\unicode[STIX]{x1D70C}^{\ast }r_{0}+P)\text{e}^{\text{i}(m\unicode[STIX]{x1D703}+\unicode[STIX]{x1D714}t)}]_{-}^{+}+O(\unicode[STIX]{x1D716}^{2})\}\hat{\boldsymbol{r}}\nonumber\\ \displaystyle & & \displaystyle +\,\left\{\frac{\unicode[STIX]{x1D716}\text{i}m}{r_{0}}[p^{\ast }]_{-}^{+}\,\text{e}^{\text{i}(m\unicode[STIX]{x1D703}+\unicode[STIX]{x1D714}t)}+O(\unicode[STIX]{x1D716}^{2})\right\}\hat{\unicode[STIX]{x1D73D}}.\end{eqnarray}$$

For the flow under consideration, described in terms of (2.6), the rate of strain tensor $\unicode[STIX]{x1D65A}$ is given by

(2.21) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D65A}=\frac{\unicode[STIX]{x1D716}}{2r^{2}}\left(\begin{array}{@{}cc@{}}2\,\text{i}\,m(r\unicode[STIX]{x1D719}^{\prime }-\unicode[STIX]{x1D719}) & -r^{2}\unicode[STIX]{x1D719}^{\prime \prime }+r\unicode[STIX]{x1D719}^{\prime }-m^{2}\unicode[STIX]{x1D719}\\ -r^{2}\unicode[STIX]{x1D719}^{\prime \prime }+r\unicode[STIX]{x1D719}^{\prime }-m^{2}\unicode[STIX]{x1D719} & -2\,\text{i}\,m(r\unicode[STIX]{x1D719}^{\prime }-\unicode[STIX]{x1D719})\end{array}\right)\text{e}^{\text{i}(m\unicode[STIX]{x1D703}+\unicode[STIX]{x1D714}t)}. & & \displaystyle\end{eqnarray}$$

Hence,

(2.22) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E5}\{\unicode[STIX]{x1D707}\unicode[STIX]{x1D65A}\boldsymbol{\cdot }\hat{\boldsymbol{n}}\} & = & \displaystyle \left\{\frac{\unicode[STIX]{x1D716}\text{i}m}{r_{0}^{2}}[\unicode[STIX]{x1D707}^{\ast }(r_{0}\unicode[STIX]{x1D719}^{\prime }-\unicode[STIX]{x1D719})]_{-}^{+}\,\text{e}^{\text{i}(m\unicode[STIX]{x1D703}+\unicode[STIX]{x1D714}t)}+O(\unicode[STIX]{x1D716}^{2})\right\}\hat{\boldsymbol{r}}\nonumber\\ \displaystyle & & \displaystyle +\left\{\frac{\unicode[STIX]{x1D716}}{2r_{0}^{2}}[\unicode[STIX]{x1D707}^{\ast }(-r_{0}^{2}\unicode[STIX]{x1D719}^{\prime \prime }+r_{0}\unicode[STIX]{x1D719}^{\prime }-m^{2}\unicode[STIX]{x1D719})]_{-}^{+}\,\text{e}^{\text{i}(m\unicode[STIX]{x1D703}+\unicode[STIX]{x1D714}t)}+O(\unicode[STIX]{x1D716}^{2})\right\}\hat{\unicode[STIX]{x1D73D}}.\end{eqnarray}$$

Combining (2.16), (2.17), (2.18), (2.20) and (2.22) with the stress continuity condition (2.13) we have at $O(1)$ and $O(\unicode[STIX]{x1D716})$ in the normal direction and $O(\unicode[STIX]{x1D716})$ in the tangential direction, respectively, the following conditions for stress continuity across $\mathscr{S}$

(2.23a ) $$\begin{eqnarray}\displaystyle & \displaystyle [p^{\ast }]_{-}^{+}=-\frac{1}{We}\frac{1}{r_{0}}, & \displaystyle\end{eqnarray}$$

(2.23b ) $$\begin{eqnarray}\displaystyle & \displaystyle \left[\unicode[STIX]{x1D70C}^{\ast }r_{0}+P-\frac{2\unicode[STIX]{x1D707}^{\ast }}{Re}\frac{\text{i}m}{r_{0}^{2}}(r_{0}\unicode[STIX]{x1D719}^{\prime }-\unicode[STIX]{x1D719})\right]_{-}^{+}=-\frac{1}{We}\frac{m^{2}-1}{r_{0}^{2}}, & \displaystyle\end{eqnarray}$$

(2.23c ) $$\begin{eqnarray}\displaystyle & \displaystyle [\unicode[STIX]{x1D707}^{\ast }(r_{0}^{2}\unicode[STIX]{x1D719}^{\prime \prime }-r_{0}\unicode[STIX]{x1D719}^{\prime }+m^{2}\unicode[STIX]{x1D719})]_{-}^{+}=0, & \displaystyle\end{eqnarray}$$

The kinematic condition at the interface, $\mathscr{S}$ , is that the fluid at the interface should move with the velocity of the interface, hence

(2.24) $$\begin{eqnarray}\displaystyle \frac{\text{D}}{\text{D}t}(r-[r_{0}+\unicode[STIX]{x1D716}\text{e}^{\text{i}(m\unicode[STIX]{x1D703}+\unicode[STIX]{x1D714}t)}])|_{\mathscr{S}}=0\quad \Rightarrow \quad u|_{\mathscr{S}}\sim \unicode[STIX]{x1D716}\text{i}\unicode[STIX]{x1D714}\text{e}^{\text{i}(m\unicode[STIX]{x1D703}+\unicode[STIX]{x1D714}t)}. & & \displaystyle\end{eqnarray}$$

The linearised kinematic condition in terms of (2.6) gives

(2.25) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{j}(r_{0})=\frac{\unicode[STIX]{x1D714}r_{0}}{m},\quad j=1,2. & & \displaystyle\end{eqnarray}$$

This condition may be also seen to match the normal fluid velocities at the interface. In the case of two viscous fluids, the tangential fluid velocities at the interface are also forced to match, and this condition is satisfied at $O(\unicode[STIX]{x1D716})$ when

(2.26) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{1}^{\prime }(r_{0})=\unicode[STIX]{x1D719}_{2}^{\prime }(r_{0}). & & \displaystyle\end{eqnarray}$$

Equations (2.11), (2.23), (2.25) and (2.26) are the complete set of boundary conditions for the problem.

2.3 Numerical simulation

A number of numerical simulations were performed using a volume-of-fluid method. The method was implemented using modifications of the ‘interFoam’ and ‘twoLiquidMixingLayer’ solvers, available as part of the OpenFOAM distribution (Weller et al. Reference Weller, Tabor, Jasak and Fureby1998). The standard implementation of the velocity equation was modified to account for the rotation by including both the centrifugal term and the Coriolis term, otherwise the solvers were unchanged. An extensive review and discussion of the implementation of the interFoam solver is provided in Deshpande, Anumolu & Trujillo (Reference Deshpande, Anumolu and Trujillo2012). The solver uses a volume-of-fluid approach to solve the continuity equation

(2.27) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}\boldsymbol{u})=0, & & \displaystyle\end{eqnarray}$$

and momentum equation

(2.28) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D70C}\boldsymbol{u})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}\boldsymbol{u}\boldsymbol{u}) & = & \displaystyle -\unicode[STIX]{x1D735}p+[\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D707}\unicode[STIX]{x1D735}\boldsymbol{u})+\unicode[STIX]{x1D735}\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D707}]\nonumber\\ \displaystyle & & \displaystyle -\,\int _{\unicode[STIX]{x1D6E4}}[\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{\boldsymbol{n}}]\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{x}_{s})\hat{\boldsymbol{n}}\,\text{d}\unicode[STIX]{x1D6E4}(\boldsymbol{x}_{s}),\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}$ is the surface tension coefficient, $\unicode[STIX]{x1D6FF}$ is the Dirac delta function in this instance and $\unicode[STIX]{x1D6E4}$ denotes the interface between the two phases. The momentum equation was modified by including the rotational terms, $-\unicode[STIX]{x1D70C}\,\unicode[STIX]{x1D734}\times (\unicode[STIX]{x1D734}\times \boldsymbol{x})-2\unicode[STIX]{x1D70C}\,\unicode[STIX]{x1D734}\times \boldsymbol{u}$ on the right-hand side of (2.28). A pressure-implicit with splitting of operators (PISO) (Issa Reference Issa1986) based predictor–corrector method is then used to first predict an updated velocity field and then correct and update the pressure and velocity fields, enforcing incompressibility. The circular domain was meshed with a 5-block mesh and typically contained $4.8\times 10^{5}$ cells. Time was stepped forward using the second-order Crank–Nicolson scheme, limited by the Courant–Friedrichs–Lewy (CFL) number. The schemes for calculating spatial gradients and the Laplacian were second order; the divergence was calculated using a limited second-order upwind scheme (see van Leer Reference van Leer1974).

3.1 Two-layer stable and unstable inviscid solutions: Rayleigh approximation

The simplest configuration we consider consists of two layers of inviscid fluid that have a density difference. In this case $\unicode[STIX]{x1D707}_{j}^{\ast }(r)=0$ and $\unicode[STIX]{x1D70C}_{j}^{\ast }(r)=\unicode[STIX]{x1D70C}_{j}$ is constant in each layer such that $\text{d}\unicode[STIX]{x1D70C}_{j}^{\ast }/\text{d}r=0$ . The initial hydrostatic conditions for the velocity and density are

(3.1a-c ) $$\begin{eqnarray}\displaystyle \boldsymbol{u}_{j}^{\ast }=\mathbf{0},\quad \unicode[STIX]{x1D70C}_{j}^{\ast }=\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D70C}_{1}, & r<r_{0},\\ \unicode[STIX]{x1D70C}_{2}, & r>r_{0},\end{array}\right.\quad p_{j}^{\ast }=p_{0}+\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D70C}_{1}r^{2}/2, & r<r_{0},\\ (\unicode[STIX]{x1D70C}_{1}-\unicode[STIX]{x1D70C}_{2})r_{0}^{2}/2+\unicode[STIX]{x1D70C}_{2}r^{2}/2, & r>r_{0},\end{array}\right. & & \displaystyle\end{eqnarray}$$

where $r=r_{0}\in (0,1)$ is the location of the initial unperturbed interface and $p_{0}$ is a reference pressure at the origin. The $O(1)$ pressure continuity condition (2.23a ) has been applied. In the special case of equal volumes of fluid in each layer $r_{0}=2^{-1/2}$ .

The form of $\unicode[STIX]{x1D70C}_{j}^{\ast }$ and $\unicode[STIX]{x1D707}_{j}^{\ast }$ leads to the simplification of (2.8), removing all the viscous effects leaving a one-dimensional Laplace equation that may be considered a Rayleigh equation. Specifically, we define the linear differential operator ${\mathcal{L}}$ and have

(3.2) $$\begin{eqnarray}\displaystyle {\mathcal{L}}[\unicode[STIX]{x1D719}]:=\unicode[STIX]{x1D719}_{j}^{\prime \prime }+\frac{\unicode[STIX]{x1D719}_{j}^{\prime }}{r}-\frac{m^{2}\unicode[STIX]{x1D719}_{j}}{r^{2}}=0, & & \displaystyle\end{eqnarray}$$

for $\unicode[STIX]{x1D714}\neq 0$ , $j=1,2$ , together with the pressure perturbation

(3.3) $$\begin{eqnarray}\displaystyle P_{j}=-\unicode[STIX]{x1D70C}_{j}\left(2\unicode[STIX]{x1D719}_{j}-\frac{\unicode[STIX]{x1D714}r}{m}\unicode[STIX]{x1D719}_{j}^{\prime }\right). & & \displaystyle\end{eqnarray}$$

Equation (3.2) yields power-law solutions for $m\geqslant 1$ given by

(3.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{j}=c_{j1}r^{-m}+c_{j2}r^{m},\quad j=1,2, & & \displaystyle\end{eqnarray}$$

with four unknown constants. The dispersion relation may now be found by fitting appropriate matching and boundary conditions.

We enforce the velocity regularity condition on $r=0$ and the no-penetration condition on $r=1$ . We then apply the kinematic condition and the stress continuity condition on the perturbed interface $r=r_{0}+\unicode[STIX]{x1D716}\exp \{\text{i}(m\unicode[STIX]{x1D703}+\unicode[STIX]{x1D714}t)\}$ as described in § 2.2. These five conditions with only four free constants then lead to the dispersion relation. The velocity regularity condition on $r=0$ implies that $|\unicode[STIX]{x1D719}_{1}/r|<\infty$ as $r\rightarrow 0$ and so we take $c_{11}=0$ . The no-penetration condition on $r=1$ implies that $\unicode[STIX]{x1D719}_{2}(1)=0$ and so we take $c_{22}=-c_{21}$ .

Applying the kinematic condition (2.25) then gives

(3.5a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{1}=\frac{\unicode[STIX]{x1D714}r_{0}}{m}\left(\frac{r}{r_{0}}\right)^{m},\quad \unicode[STIX]{x1D719}_{2}=\frac{\unicode[STIX]{x1D714}r_{0}}{m}\left(\frac{r}{r_{0}}\right)^{m}\left(\frac{1-r^{-2m}}{1-r_{0}^{-2m}}\right). & & \displaystyle\end{eqnarray}$$

The continuity of stress at the interface now gives the dispersion relation. In the present inviscid case, in the absence of surface tension, we have continuity of pressure across $\mathscr{S}$ . The tangential stress condition (2.23c ) is automatically satisfied as $\unicode[STIX]{x1D707}_{j}^{\ast }=0$ , and the leading-order pressure continuity condition (2.23a ) is satisfied by the hydrostatic solution (3.1). The remaining $O(\unicode[STIX]{x1D716})$ pressure continuity condition, (2.23b ), simplifies to

(3.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{1}\left[r_{0}-2\unicode[STIX]{x1D719}_{1}(r_{0})+\frac{\unicode[STIX]{x1D714}r_{0}}{m}\unicode[STIX]{x1D719}_{1}^{\prime }(r_{0})\right]=\unicode[STIX]{x1D70C}_{2}\left[r_{0}-2\unicode[STIX]{x1D719}_{2}(r_{0})+\frac{\unicode[STIX]{x1D714}r_{0}}{m}\unicode[STIX]{x1D719}_{2}^{\prime }(r_{0})\right]. & & \displaystyle\end{eqnarray}$$

Substituting in (3.5) leads to the dispersion relation

(3.7) $$\begin{eqnarray}\displaystyle \mathscr{A}(\unicode[STIX]{x1D714}^{2}-2\unicode[STIX]{x1D714}+m)=\frac{(1+\mathscr{A})\unicode[STIX]{x1D714}^{2}}{1-r_{0}^{2m}}, & & \displaystyle\end{eqnarray}$$

where we have used an Atwood number, $\mathscr{A}$ , as in Scase et al. (Reference Scase, Baldwin and Hill2017b ), defined by

(3.8) $$\begin{eqnarray}\displaystyle \mathscr{A}=\frac{\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1}}{\unicode[STIX]{x1D70C}_{2}+\unicode[STIX]{x1D70C}_{1}}\in [-1,1]. & & \displaystyle\end{eqnarray}$$

We can further simplify (3.7) by introducing a modified Atwood number, $\mathscr{A}^{\ast }$ , given by

(3.9) $$\begin{eqnarray}\displaystyle \mathscr{A}^{\ast }=\mathscr{A}\left(\frac{1-r_{0}^{2m}}{1+\mathscr{A}r_{0}^{2m}}\right)\in [-1,1], & & \displaystyle\end{eqnarray}$$

that has the same sign as $\mathscr{A}$ but has a shape factor based on the initial position of the interface, the azimuthal wavenumber and the Atwood number. In terms of the modified Atwood number we may rewrite (3.7) as

(3.10) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}^{2}+\mathscr{A}^{\ast }(2\unicode[STIX]{x1D714}-m)=0. & & \displaystyle\end{eqnarray}$$

In the special case $\mathscr{A}=\mathscr{A}^{\ast }=0$ , that occurs when there is no density contrast, we can immediately show that $\unicode[STIX]{x1D714}=0$ and there is no growth or precession of the perturbation.

Since in (2.14) we choose the form of the interface to be $r=r_{0}+\unicode[STIX]{x1D716}\exp \{\text{i}(m\unicode[STIX]{x1D703}+\unicode[STIX]{x1D714}t)\}$ , a negative imaginary part in $\unicode[STIX]{x1D714}$ corresponds to growth of the interface. The dispersion relations (3.7) and (3.10) have solutions with negative imaginary parts when their discriminants are negative. This condition is met (except for the special case $m=1$ , $\mathscr{A}=-1$ ) when $\mathscr{A}^{\ast }<0$ , i.e. $\mathscr{A}\in [-1,0)$ . We have unstable growth when the Atwood number is negative, i.e. when the density of the fluid in the inner layer is greater than that in the outer layer, as may be anticipated.

For the special case $\mathscr{A}=\mathscr{A}^{\ast }=-1$ , that occurs when there is no fluid in the outer layer, we find $\unicode[STIX]{x1D714}=1\pm \text{i}\sqrt{m-1}$ , i.e. except for $m=1$ all modes grow in amplitude. When $m=1$ the inner fluid layer, the only fluid layer present, remains circular, but is displaced, precessing about the origin. Note that $m=1$ is the only perturbation mode that gives rise to a non-zero fluid velocity at the origin of the system and may be considered a centre of mass oscillation.

The solution to (3.10) is

(3.11) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}=-\mathscr{A}^{\ast }\pm \{\mathscr{A}^{\ast }(\mathscr{A}^{\ast }+m)\}^{1/2}. & & \displaystyle\end{eqnarray}$$

When the Atwood number is negative then perturbations initially grow exponentially in time with precessional rate $-\mathscr{A}^{\ast }$ and growth rate $\{-\mathscr{A}^{\ast }(\mathscr{A}^{\ast }+m)\}^{1/2}$ . The stable dispersion relation, (3.7) with $\mathscr{A}>0$ , necessarily has one positive root and one negative root, denoted $\unicode[STIX]{x1D714}_{\infty }^{\pm }$ , as a result of $m(1+\mathscr{A}r_{0}^{2m})>0$ . One wave precesses clockwise, while the other precesses in an anticlockwise direction.

In the limit of large azimuthal wavenumber, $m\rightarrow \infty$ , then $r_{0}^{2m}\rightarrow 0$ and

(3.12) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}\sim -\mathscr{A}\pm \{\mathscr{A}(\mathscr{A}+m)\}^{1/2}, & & \displaystyle\end{eqnarray}$$

independently of $r_{0}$ , as the curvature of the interface is not ‘felt’ by the system for the high azimuthal wavenumbers since the azimuthal wavelength is too short compared to $r_{0}$ . The classical two-dimensional non-rotating dispersion relation due to Rayleigh (Reference Rayleigh1883) and Taylor (Reference Taylor1950) in a domain of width $2\unicode[STIX]{x03C0}r_{0}$ for perturbations with wavenumber $k=m/r_{0}$ is $\unicode[STIX]{x1D714}^{2}=g\mathscr{A}m/r_{0}$ (using dimensional quantities) and so we observe that the rotating solution at large $m$ behaves as the classical solution, with no curvature of the interface and an effective gravitational acceleration given by $g\sim r_{0}\unicode[STIX]{x1D6FA}^{2}$ .

We can gain further insight into the physical mechanisms at work by proceeding dimensionally. Distinguishing between the centrifugal term in the governing equation of motion (2.1b ), $\unicode[STIX]{x1D734}\times (\unicode[STIX]{x1D734}\times \boldsymbol{x})$ and the Coriolis term, $2\hat{\unicode[STIX]{x1D734}}\times \boldsymbol{u}$ , by introducing a ‘hat’ on the rotation vector in the Coriolis term, the dimensional dispersion relation (3.10) is given by

(3.13) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}^{2}+\mathscr{A}^{\ast }(2\hat{\unicode[STIX]{x1D6FA}}\unicode[STIX]{x1D714}-\unicode[STIX]{x1D6FA}^{2}m)=0, & & \displaystyle\end{eqnarray}$$

where $\mathscr{A}^{\ast }=\mathscr{A}(a^{2m}-r_{0}^{2m})/(a^{2m}+\mathscr{A}r_{0}^{2m})$ and $\hat{\unicode[STIX]{x1D6FA}}\equiv \unicode[STIX]{x1D6FA}$ : the hat notation is used only to highlight terms originating from the Coriolis force. The constant term in (3.13) is now identifiable with the centrifugal term from the equation of motion (2.1b ), while the linear term is associated with the Coriolis term. The discriminant of (3.13) is proportional to $\mathscr{A}^{\ast 2}\hat{\unicode[STIX]{x1D6FA}}^{2}+\mathscr{A}^{\ast }\unicode[STIX]{x1D6FA}^{2}m$ where the first term can now be associated with Coriolis effects and the second term with centrifugal effects. In order for a perturbation to grow in time the discriminant must be negative. However, the first term of the discriminant is always positive and so Coriolis effects can never lead to unstable growth. Any growth that occurs is due to the second term of the discriminant which is associated with the centrifugal term in (2.1b ) and hence any growth is therefore ‘centrifugally driven’. For any configuration where we have centrifugally driven growth, the Coriolis term must slow the growth rate but can never prevent growth as a result of $|\mathscr{A}^{\ast 2}\hat{\unicode[STIX]{x1D6FA}}^{2}|<|\mathscr{A}^{\ast }\unicode[STIX]{x1D6FA}^{2}m|$ .

It may also be observed from the dimensional form of the modified Atwood number that the limit $\mathscr{A}^{\ast }\rightarrow \mathscr{A}$ as $m\rightarrow \infty$ may also be realised by considering the outer fluid to be unbounded, i.e. taking $a\rightarrow \infty$ keeping $m$ and $r_{0}$ fixed. The conclusion of the equivalence of the present system and the classical two-dimensional system in the limit of high azimuthal wavenumber, $m$ , therefore also follows for the limit of an unbounded outer fluid layer, $a\rightarrow \infty$ .

The dispersion relation of Tao et al. (Reference Tao, He, Ye and Busse2013) is recovered in the limit of: a hydrostatic initial interface, $g=R\unicode[STIX]{x1D6FA}^{2}$ (in their notation); and $-kR=m$ where $m$ is our azimuthal wavenumber and $k$ is their Cartesian wavenumber. The minus sign accounts for the difference in the sense of our $\hat{\unicode[STIX]{x1D73D}}$ and their $\hat{\boldsymbol{x}}$ . Finally $\mathscr{A}^{\ast }$ is replaced by $-\mathscr{A}$ , due to neglecting the curvature of the interface and taking an unbounded outer domain ( $r_{0}/a\rightarrow 0$ ), the sign change is due to the difference in definition of the Atwood number between the two treatments.

3.2 Surface tension between immiscible fluids

Returning to non-dimensional quantities, we now consider the effect of surface tension acting between two inviscid immiscible fluid layers. The pressure field in the hydrostatic solution is modified from the previous case considered in § 3.1, due to the finite Weber number in (2.23a ), and is now given by

(3.14) $$\begin{eqnarray}\displaystyle p_{j}^{\ast }=p_{0}+\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D70C}_{1}r^{2}/2+1/(We\,r_{0}), & r<r_{0},\\ (\unicode[STIX]{x1D70C}_{1}-\unicode[STIX]{x1D70C}_{2})r_{0}^{2}/2+\unicode[STIX]{x1D70C}_{2}r^{2}/2, & r>r_{0}\end{array}\right. & & \displaystyle\end{eqnarray}$$

while the hydrostatic, velocity and density fields remain as before. The kinematic conditions are as in (2.25), hence (3.5) remains unchanged. The $O(\unicode[STIX]{x1D716})$ normal stress continuity condition at the interface (2.23b ) is

(3.15) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{1}\left(r_{0}-2\unicode[STIX]{x1D719}_{1}(r_{0})+\frac{\unicode[STIX]{x1D714}r_{0}}{m}\unicode[STIX]{x1D719}_{1}^{\prime }(r_{0})\right)=\unicode[STIX]{x1D70C}_{2}\left(r_{0}-2\unicode[STIX]{x1D719}_{2}(r_{0})+\frac{\unicode[STIX]{x1D714}r_{0}}{m}\unicode[STIX]{x1D719}_{2}^{\prime }(r_{0})\right)+\frac{1}{We}\frac{m^{2}-1}{r_{0}^{2}}. & & \displaystyle\end{eqnarray}$$

This modifies (3.10) to give

(3.16) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}^{2}+\mathscr{A}^{\ast }(2\unicode[STIX]{x1D714}-m)=\frac{\mathscr{A}^{\ast }}{We}\frac{m(m^{2}-1)}{2\mathscr{A}r_{0}^{3}}=:\mathscr{A}^{\ast }S(We;\mathscr{A},m,r_{0}), & & \displaystyle\end{eqnarray}$$

defining the surface tension adjustment factor $S$ . It follows immediately from the form of (3.16) that surface tension has the greatest effect when the wavenumber, $m$ , is high and when the mean radius of the interface, $r_{0}$ , is small, as expected on physical grounds. The result (3.16) agrees with the results for a rotating column of liquid (Hocking & Michael Reference Hocking and Michael1959) and a rotating inviscid drop (Patzek et al. Reference Patzek, Basaran, Benner and Scriven1995), up to a change of frame of reference and non-dimensionalisation, in the special case $\mathscr{A}=-1$ .

The solution to (3.16), modified compared to (3.11), is

(3.17) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}=-\mathscr{A}^{\ast }\pm \{\mathscr{A}^{\ast }(\mathscr{A}^{\ast }+m+S)\}^{1/2}. & & \displaystyle\end{eqnarray}$$

Inspection of (3.17) and comparison with (3.11) shows that the effect of the surface tension can be interpreted as modifying the azimuthal wavenumber. In a stable configuration the effect is to increase the apparent azimuthal wavenumber, as $S>0$ , enhancing the frequency of oscillation, as might be anticipated since the surface tension applies an additional restorative force on the interface. In an unstable configuration the effect is to decrease the apparent azimuthal wavenumber, as $S<0$ , inhibiting the growth of the instability. The surface tension affects the constant term in the quadratic dispersion relation (3.17) and so influences the contribution of the centrifugal forcing to the system. It is therefore possible for the surface tension not only to influence the growth rate of unstable modes but indeed completely stabilise an otherwise unstable mode if the surface tension is large enough. It follows from (3.17) that for a naturally unstable mode, $m$ , the surface tension is able to completely stabilise the mode when $S\leqslant -(\mathscr{A}^{\ast }+m)$ . Alternatively, we observe that the surface tension is able to completely suppress the growth of high frequency modes above a critical value $m^{\ast }$ given by the implicit relationship $-m^{\ast }=\mathscr{A}^{\ast }+S(We;\mathscr{A},m^{\ast },r_{0})$ or

(3.18) $$\begin{eqnarray}\displaystyle \mathscr{A}(1-r_{0}^{2m^{\ast }})+(1+\mathscr{A}r_{0}^{2m^{\ast }})m^{\ast }\left[1+\frac{1}{We}\frac{(m^{\ast 2}-1)}{2\mathscr{A}r_{0}^{3}}\right]=0. & & \displaystyle\end{eqnarray}$$

When $r_{0}$ is close to the outer boundary, say $r_{0}=1-\unicode[STIX]{x1D700}$ , $\unicode[STIX]{x1D700}\ll 1$ , (n.b.  $\unicode[STIX]{x1D700}$ is different to the small parameter $\unicode[STIX]{x1D716}$ in (2.4), it is assumed that $0\ll \unicode[STIX]{x1D716}\ll \unicode[STIX]{x1D700}\ll 1$ ) then we find

(3.19) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}\sim m\left\{\pm \left(1+\frac{1}{We}\frac{m^{2}-1}{2\mathscr{A}}\right)^{1/2}\left(\frac{2\mathscr{A}}{1+\mathscr{A}}\unicode[STIX]{x1D700}\right)^{1/2}-\frac{2\mathscr{A}}{1+\mathscr{A}}\unicode[STIX]{x1D700}+O(\unicode[STIX]{x1D700}^{3/2})\right\}, & & \displaystyle\end{eqnarray}$$

and we see that a given mode with wavenumber $m$ is stable for $We>-(m^{2}-1)/2\mathscr{A}$ (cf. (3.18) as $r_{0}\rightarrow 1$ ).

3.3 Diffusion of the interface between miscible fluids

It is well known in the study of Kelvin–Helmholtz instability that diffusion of the interface between the two fluid layers can lead to some modes of instability being suppressed. Here we investigate whether a diffuse fluid layer between the inner and outer fluids causes some modes of centrifugally driven Rayleigh–Taylor instability to be suppressed similarly. We begin by considering the form of a diffuse interface between the two layers. We consider a two-layer stratification with a sharp interface (as in (3.1), with $\unicode[STIX]{x1D70C}^{\ast }=\unicode[STIX]{x1D70C}_{1}$ in the inner layer and $\unicode[STIX]{x1D70C}^{\ast }=\unicode[STIX]{x1D70C}_{2}$ in the outer layer) that is subject to diffusion of density governed by the dimensional diffusion equation

(3.20) $$\begin{eqnarray}\displaystyle \frac{\text{D}\unicode[STIX]{x1D70C}^{\ast }}{\text{D}t}=\unicode[STIX]{x1D705}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D70C}^{\ast }, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D705}$ is a constant diffusivity. A natural scaling for this problem in isolation is to take $t=a^{2}\unicode[STIX]{x1D705}^{-1}t^{\prime }$ , however to enable comparison with the temporal development of the instability we choose $t=\unicode[STIX]{x1D6FA}^{-1}t^{\prime }$ . Therefore, non-dimensionalising as in § 2 and dropping the prime notation we obtain

(3.21) $$\begin{eqnarray}\displaystyle \frac{\text{D}\unicode[STIX]{x1D70C}^{\ast }}{\text{D}t}=\frac{1}{Pe}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D70C}^{\ast }, & & \displaystyle\end{eqnarray}$$

where $Pe=\unicode[STIX]{x1D6FA}a^{2}/\unicode[STIX]{x1D705}$ is a non-dimensional Péclet number that relates the advection time scale to the diffusion time scale. We consider the diffusion of density in the hydrostatic fluid such that $\unicode[STIX]{x1D70C}^{\ast }=\unicode[STIX]{x1D70C}^{\ast }(r,t)$ , $\boldsymbol{u}^{\ast }=\mathbf{0}$ and hence

(3.22) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\ast }}{\unicode[STIX]{x2202}t}=\frac{1}{Pe}\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\ast }}{\unicode[STIX]{x2202}r}\right),\quad r\in (0,1),\,t\in [0,\infty ). & & \displaystyle\end{eqnarray}$$

The initial condition is as in (3.1) and we enforce Neumann boundary conditions $\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\ast }/\unicode[STIX]{x2202}r=0$ on $r=0,1$ . By separation of variables a series solution to the above problem can be found and is given by

(3.23) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}^{\ast }=1+\mathscr{A}-2\mathscr{A}r_{0}\left\{r_{0}+2\mathop{\sum }_{n=1}^{\infty }\frac{\text{J}_{1}(\unicode[STIX]{x1D706}_{n}r_{0})\text{J}_{0}(\unicode[STIX]{x1D706}_{n}r)}{\unicode[STIX]{x1D706}_{n}\text{J}_{0}(\unicode[STIX]{x1D706}_{n})^{2}}\text{e}^{-\unicode[STIX]{x1D706}_{n}^{2}t/Pe}\right\}, & & \displaystyle\end{eqnarray}$$

where $\text{J}$ is a Bessel function of the first kind and $\unicode[STIX]{x1D706}_{n}$ for $n=1,2,\ldots$ are the zeros of $\text{J}_{1}$ , i.e. $\unicode[STIX]{x1D706}_{1}\approx 3.83$ , $\unicode[STIX]{x1D706}_{2}\approx 7.02$ etc. It can be seen that in the long-time limit the density is constant everywhere equal to the mean density. The ‘thickness’ of the diffuse layer, to be denoted by $\unicode[STIX]{x1D6FF}$ , is subjective, but one method is to define the thickness by

(3.24) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}=(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})\left(\left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\ast }}{\unicode[STIX]{x2202}r}\right|_{r=r_{0}}\right)^{-1}, & & \displaystyle\end{eqnarray}$$

that follows from fitting a piecewise linear density profile that has the same gradient as the full solution at $r=r_{0}$ . An approximate scaling for $\unicode[STIX]{x1D6FF}$ can be found as follows. If we consider the diffuse layer to be small, such that the rapid changes in $\unicode[STIX]{x1D70C}^{\ast }$ take place in a narrow region $r=r_{0}+\unicode[STIX]{x1D700}x$ , where $\unicode[STIX]{x1D700}\ll 1$ (and again, $\unicode[STIX]{x1D700}$ is an arbitrary small parameter, not related to $\unicode[STIX]{x1D716}$ in § 2) and the time scale for diffusion is taken such that $t=\unicode[STIX]{x1D700}^{2}Pe\,\unicode[STIX]{x1D70F}$ then (3.22) becomes

(3.25) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\ast }}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}=\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D70C}^{\ast }}{\unicode[STIX]{x2202}x^{2}}+O(\unicode[STIX]{x1D700}), & & \displaystyle\end{eqnarray}$$

which may therefore be considered Cartesian at leading order as the curvature effects are negligible. This form of the diffusion equation accepts a similarity solution

(3.26) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}^{\ast }=\frac{\unicode[STIX]{x1D70C}_{1}+\unicode[STIX]{x1D70C}_{2}}{2}+\frac{(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})}{2}\text{erf}\left(\frac{x}{2\sqrt{\unicode[STIX]{x1D70F}}}\right), & & \displaystyle\end{eqnarray}$$

where $\text{erf}(x)$ is the error function $2\unicode[STIX]{x03C0}^{-1/2}\int _{0}^{x}\exp \{-\unicode[STIX]{x1D709}^{2}\}\,\text{d}\unicode[STIX]{x1D709}$ and $\unicode[STIX]{x1D70C}^{\ast }$ satisfies the boundary conditions $\unicode[STIX]{x1D70C}^{\ast }\rightarrow \unicode[STIX]{x1D70C}_{1}$ as $x\rightarrow -\infty$ , $\unicode[STIX]{x1D70C}^{\ast }\rightarrow \unicode[STIX]{x1D70C}_{2}$ as $x\rightarrow \infty$ , and the initial conditions (3.1). Therefore, at early times we may approximate the diffuse layer given by (3.23) as

(3.27) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}^{\ast }=1+\mathscr{A}\text{erf}\left(\frac{r-r_{0}}{2\sqrt{t/Pe}}\right). & & \displaystyle\end{eqnarray}$$

Figure 2 shows the solution of (3.20) as given in (3.23) (thick solid line) compared to the early time approximation (3.27) (thin solid line, indistinguishable from the full solution at the scale shown). The form of solution in (3.27) leads to a simple approximation for $\unicode[STIX]{x1D6FF}$ , as defined in (3.24), giving

(3.28) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}\sim 2\sqrt{\unicode[STIX]{x03C0}t/Pe}+O(\unicode[STIX]{x1D700}), & & \displaystyle\end{eqnarray}$$

which is also shown in figure 2.

Figure 2.

A comparison of the approximations to the diffuse layer for $\mathscr{A}=1/2$ , $t/Pe=10^{-3}$ , $r_{0}=2^{-1/2}$ . The exact solution (3.23) is shown as solid bold. The similarity solution (3.27) is plotted but cannot be distinguished by eye at this scale from the exact solution. The piecewise power-law approximation (3.29) is shown as a solid line. The piecewise linear approximation is shown dashed for comparison. Dotted lines indicate the width of the diffuse layer $\unicode[STIX]{x1D6FF}$ , given in (3.28).

We consider the effect of diffusion on a system with large Péclet number, where diffusion is slow compared to the growth of an instability, or the oscillation of an interfacial wave. This would correspond to an experimental situation where the interface had been allowed to diffuse before the experiment was conducted, but during the experiment the effects of further diffusion could be neglected. As such we seek to approximate the effect of diffusion by constructing a piecewise density profile. Examining (2.8) we see that if we can approximate the diffuse interface by a piecewise density $\unicode[STIX]{x1D70C}^{\ast }$ such that $\unicode[STIX]{x1D70C}^{\ast }$ is of the form $\unicode[STIX]{x1D70C}^{\ast }=\unicode[STIX]{x1D6FD}r^{\unicode[STIX]{x1D6FC}}$ for some constants $\unicode[STIX]{x1D6FD}$ and $\unicode[STIX]{x1D6FC}$ , then the form of (2.8) is unchanged, and power-law solutions are still admitted. We define $r_{0-}=r_{0}-\unicode[STIX]{x1D6FF}/2$ , the inner edge of the diffuse layer, and $r_{0+}=r_{0}+\unicode[STIX]{x1D6FF}/2$ , the outer edge of the diffuse layer and identify a three-layer system $(j\in \{1,2,3\})$ where

(3.29) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}^{\ast }=\left\{\begin{array}{@{}lcl@{}}\unicode[STIX]{x1D70C}_{1}, & 0\leqslant r<r_{0-}; & j=1\\ \unicode[STIX]{x1D6FD}r^{\unicode[STIX]{x1D6FC}}, & r_{0-}\leqslant r<r_{0+}; & j=2\\ \unicode[STIX]{x1D70C}_{2}, & r_{0+}\leqslant r\leqslant 1; & j=3\end{array}\right. & & \displaystyle\end{eqnarray}$$

and where $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ are chosen to make the density continuous and are given by

(3.30a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}=\log \left(\frac{1-\mathscr{A}}{1+\mathscr{A}}\right)\left[\log \left(\frac{r_{0-}}{r_{0+}}\right)\right]^{-1},\quad \unicode[STIX]{x1D6FD}=\frac{1-\mathscr{A}}{r_{0-}^{\unicode[STIX]{x1D6FC}}}. & & \displaystyle\end{eqnarray}$$

It follows that

(3.31) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D70C}^{\ast \prime }}{\unicode[STIX]{x1D70C}^{\ast }}=\left\{\begin{array}{@{}ll@{}}0, & 0\leqslant r<r_{0-}\\ \unicode[STIX]{x1D6FC}/r, & r_{0-}\leqslant r<r_{0+}\\ 0, & r_{0+}\leqslant r\leqslant 1.\end{array}\right. & & \displaystyle\end{eqnarray}$$

The corresponding initial hydrostatic pressure field is

(3.32) $$\begin{eqnarray}\displaystyle p^{\ast }=p_{0}+\left\{\begin{array}{@{}ll@{}}\displaystyle \frac{\unicode[STIX]{x1D70C}_{1}}{2}r^{2}, & 0\leqslant r<r_{0-}\\ \displaystyle \frac{\unicode[STIX]{x1D70C}_{1}}{2}r_{0-}^{2}+\frac{\unicode[STIX]{x1D6FD}}{2+\unicode[STIX]{x1D6FC}}[r^{2+\unicode[STIX]{x1D6FC}}-r_{0-}^{2+\unicode[STIX]{x1D6FC}}], & r_{0-}\leqslant r<r_{0+}\\ \displaystyle \frac{\unicode[STIX]{x1D70C}_{1}}{2}r_{0-}^{2}+\frac{\unicode[STIX]{x1D6FD}}{2+\unicode[STIX]{x1D6FC}}[r_{0+}^{2+\unicode[STIX]{x1D6FC}}-r_{0-}^{2+\unicode[STIX]{x1D6FC}}]-\frac{\unicode[STIX]{x1D70C}_{2}}{2}r_{0+}^{2}+\frac{\unicode[STIX]{x1D70C}_{2}}{2}r^{2}, & r_{0+}\leqslant r\leqslant 1.\end{array}\right. & & \displaystyle\end{eqnarray}$$

The governing equation (2.8) is as in the two-layer no-diffusion case in the inner ( $j=1$ ) and outer ( $j=3$ ) layers and so we have

(3.33a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{1}(r)=c_{11}r^{-m}+c_{12}r^{m},\quad \unicode[STIX]{x1D719}_{3}(r)=c_{31}r^{-m}+c_{32}r^{m}. & & \displaystyle\end{eqnarray}$$

The governing equation (2.8) in the transitional ( $j=2$ ) layer is

(3.34) $$\begin{eqnarray}\displaystyle \left(\unicode[STIX]{x1D719}_{2}^{\prime \prime }+\frac{\unicode[STIX]{x1D719}_{2}^{\prime }}{r}-\frac{m^{2}\unicode[STIX]{x1D719}_{2}}{r^{2}}\right)+\frac{\unicode[STIX]{x1D6FC}}{r}\left(\unicode[STIX]{x1D719}_{2}^{\prime }+\frac{m(m-2\unicode[STIX]{x1D714})\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x1D714}^{2}r}\right)=0. & & \displaystyle\end{eqnarray}$$

This yields

(3.35) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{2}^{\prime \prime }+(1+\unicode[STIX]{x1D6FC})\frac{\unicode[STIX]{x1D719}_{2}^{\prime }}{r}-\left(m^{2}-\frac{m(m-2\unicode[STIX]{x1D714})\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D714}^{2}}\right)\frac{\unicode[STIX]{x1D719}_{2}}{r^{2}}=0, & & \displaystyle\end{eqnarray}$$

with the power-law solution

(3.36) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{2}(r)=c_{21}r^{-\unicode[STIX]{x1D6FC}/2-\unicode[STIX]{x1D712}}+c_{22}r^{-\unicode[STIX]{x1D6FC}/2+\unicode[STIX]{x1D712}},\quad \text{where}~\unicode[STIX]{x1D712}=\frac{1}{2}\left\{\unicode[STIX]{x1D6FC}^{2}+4\left(m^{2}-\frac{m(m-2\unicode[STIX]{x1D714})\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D714}^{2}}\right)\right\}^{1/2}. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

We apply the velocity regularity condition at $r=0$ and no-penetration condition at $r=1$ to find $c_{11}=0$ and $c_{32}=-c_{31}$ . The boundary conditions are as derived in § 2.2, but are applied at $r_{0-}$ and $r_{0+}$ . These may be shown to be that $\unicode[STIX]{x1D719}_{1}=\unicode[STIX]{x1D719}_{2}$ and $\unicode[STIX]{x1D719}_{1}^{\prime }=\unicode[STIX]{x1D719}_{2}^{\prime }$ at $r=r_{0-}$ , and $\unicode[STIX]{x1D719}_{2}=\unicode[STIX]{x1D719}_{3}$ and $\unicode[STIX]{x1D719}_{2}^{\prime }=\unicode[STIX]{x1D719}_{3}^{\prime }$ at $r=r_{0+}$ . Enforcing these matching conditions, to write $c_{12}$ , $c_{22}$ at $r_{0-}$ and $c_{22}$ and $c_{31}$ at $r_{0+}$ in terms of $c_{12}$ leaves two expressions for the ratio $c_{21}/c_{22}$ which, combined, yield the dispersion relation

(3.37) $$\begin{eqnarray}\displaystyle \frac{2(\unicode[STIX]{x1D712}+m)+\unicode[STIX]{x1D6FC}}{2(\unicode[STIX]{x1D712}-m)-\unicode[STIX]{x1D6FC}}\left[\frac{(2(\unicode[STIX]{x1D712}-m)-\unicode[STIX]{x1D6FC})r_{0+}^{2m}-2(\unicode[STIX]{x1D712}+m)+\unicode[STIX]{x1D6FC}}{(2(\unicode[STIX]{x1D712}+m)+\unicode[STIX]{x1D6FC})r_{0+}^{2m}-2(\unicode[STIX]{x1D712}-m)-\unicode[STIX]{x1D6FC}}\right]=\left(\frac{r_{0-}}{r_{0+}}\right)^{2\unicode[STIX]{x1D712}}. & & \displaystyle\end{eqnarray}$$

For small diffuse layers where $\unicode[STIX]{x1D6FF}\ll 1$ , and $\unicode[STIX]{x1D6FF}$ may be approximated using (3.28) we have from (3.37), that

(3.38) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D712}^{2} & {\sim} & \displaystyle \left(\frac{\unicode[STIX]{x1D6FC}_{0}r_{0}}{2\unicode[STIX]{x1D6FF}}\right)^{2}\left\{1+\frac{4m}{\mathscr{A}^{\ast }\unicode[STIX]{x1D6FC}_{0}r_{0}}\unicode[STIX]{x1D6FF}\right.\nonumber\\ \displaystyle & & \displaystyle \left.-\left[\frac{1}{6r_{0}^{2}}-\left(\frac{2m}{\unicode[STIX]{x1D6FC}_{0}r_{0}}\right)^{2}\left(1+\frac{1+\mathscr{A}}{\mathscr{A}^{3}}\frac{(2\mathscr{A}+\unicode[STIX]{x1D6FC}_{0})[1-\mathscr{A}(1-r_{0}^{2m})]}{(1-r_{0}^{2m})^{2}}\right)\right]\unicode[STIX]{x1D6FF}^{2}+O(\unicode[STIX]{x1D6FF}^{3})\right\},\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}_{0}=\log [(1-\mathscr{A})/(1+\mathscr{A})]$ . It follows that

(3.39) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}\sim \unicode[STIX]{x1D714}_{\infty }\left\{1-\left(\frac{m}{\unicode[STIX]{x1D6FC}_{0}r_{0}}\right)\left(\frac{\mathscr{A}^{\ast }\unicode[STIX]{x1D714}_{\infty }}{\mathscr{A}^{\ast }+\unicode[STIX]{x1D714}_{\infty }}\right)\frac{1+\mathscr{A}}{2\mathscr{A}^{3}}\frac{(2\mathscr{A}+\unicode[STIX]{x1D6FC}_{0})[1-\mathscr{A}(1-r_{0}^{2m})]}{(1-r_{0}^{2m})^{2}}\unicode[STIX]{x1D6FF}+O(\unicode[STIX]{x1D6FF}^{2})\right\}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D714}_{\infty }$ is the zero-diffusion solution from (3.11). It follows immediately that as the thickness of the diffuse layer, $\unicode[STIX]{x1D6FF}$ , tends to zero, the zero-diffusion solution $\unicode[STIX]{x1D714}_{\infty }$ is recovered.

Figure 3.

Comparison of the effect of the various models of the diffuse interface on $\unicode[STIX]{x1D719}$ . The no-diffusion case is the solid piecewise line comprising two solid thin line segments that are joined at the black circle. The piecewise diffuse case is the thick black line, the solution is in three segments that coincide at the white circles. $\unicode[STIX]{x1D719}_{1}$ continues dot-dashed after the first white circle, $\unicode[STIX]{x1D719}_{2}$ outside the diffuse region is shown dashed and $\unicode[STIX]{x1D719}_{3}$ is shown dotted to the left of the second white circle. The numerical solutions corresponding to the exact solution of the diffusion equation and the similarity solution of the diffusion equation are shown as a thin solid line and a dashed line, but the difference cannot be seen by eye at this scale. The parameters chosen were: $\mathscr{A}=-1/2$ , $t=10^{-3}Pe$ , $m=5$ , $r_{0}=2^{-1/2}$ . The corresponding values of $\unicode[STIX]{x1D714}$ were: for the non-diffuse case $\unicode[STIX]{x1D714}=0.49-1.49\text{i}$ ; for the piecewise diffusion solution $\unicode[STIX]{x1D714}=0.40-1.35\text{i}$ ; for the exact and similarity solution diffusion solutions $\unicode[STIX]{x1D714}=0.37-1.31\text{i}$ .

Figure 3 compares the solutions $\unicode[STIX]{x1D719}$ for the various approximations to the diffuse layer with the sharp zero-diffusion solutions. The zero-diffusion solution (3.5) is shown as a thin solid line where the transition between $\unicode[STIX]{x1D719}_{1}$ and $\unicode[STIX]{x1D719}_{2}$ takes place at the solid black circle at a cusp. For the chosen parameters, $\mathscr{A}=-1/2$ , $m=5$ , $r_{0}=2^{-1/2}$ , the dispersion relation (3.7) gives $\unicode[STIX]{x1D714}=0.49-1.49\text{i}$ for the most unstable mode. We then considered the change in $\unicode[STIX]{x1D714}$ if the interface was allowed to diffuse over a time $t=10^{-3}Pe$ prior to any perturbation to the system. We can use the exact solution for $\unicode[STIX]{x1D70C}^{\ast }$  (3.23) and then numerically solve the appropriate form of (2.8) given by

(3.40) $$\begin{eqnarray}\displaystyle \left(\unicode[STIX]{x1D719}^{\prime \prime }+\frac{\unicode[STIX]{x1D719}^{\prime }}{r}-\frac{m^{2}\unicode[STIX]{x1D719}}{r^{2}}\right)+\frac{\unicode[STIX]{x1D70C}^{\ast \prime }}{\unicode[STIX]{x1D70C}^{\ast }}\left(\unicode[STIX]{x1D719}^{\prime }+\frac{m(m-2\unicode[STIX]{x1D714})\unicode[STIX]{x1D719}}{\unicode[STIX]{x1D714}^{2}r}\right)=0, & & \displaystyle\end{eqnarray}$$

across the whole domain to find $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D714}$ . This solution is shown as a thin solid line, where the kinematic condition has been enforced at $r=r_{0}$ indicated by the lower black circle. The solution to (3.40) subject to $\unicode[STIX]{x1D719}(0)=0$ , $\unicode[STIX]{x1D719}(1)=0$ was calculated separately in $[0,r_{0}]$ and $[r_{0},1]$ with the kinematic condition enforced for each solution at $r=r_{0}$ . The eigenvalue $\unicode[STIX]{x1D714}$ is chosen to ensure that $\unicode[STIX]{x1D719}$ is at least class $C^{2}[0,1]$ . The method of solution was to formulate the problem in terms of Chebyshev polynomials (see Driscoll, Hale & Trefethen (Reference Driscoll, Hale and Trefethen2014) for details) and the method generalises well to the fully viscous case, as is discussed in § 3.5. The same method can be implemented using the similarity solution (3.27), and this is shown as a dashed line, though it is hard to distinguish from the exact solution at this scale. These numerical solutions yield an eigenvalue $\unicode[STIX]{x1D714}=0.37-1.31\text{i}$ for the most unstable mode, showing that the diffusion inhibits the most unstable mode. Finally, the piecewise solution is shown as the thick bold line and comprises three segments that transition at the white circles located on $r=r_{0\pm }$ . The $\unicode[STIX]{x1D719}_{1}$ solution is plotted, dot-dashed, to the right of the left-hand white circle though it is no longer part of the solution. Similarly, $\unicode[STIX]{x1D719}_{2}$ is shown to the left of $r_{0-}$ and the right of $r_{0+}$ dashed. The final part of the solution, $\unicode[STIX]{x1D719}_{3}$ , is shown to the left of the right-hand white circle, dotted, before it becomes part of the solution. The predicted most unstable mode in this case has eigenvalue $\unicode[STIX]{x1D714}=0.40-1.35\text{i}$ .

3.4 Viscous fluid layers

We now consider the case of two fluid layers separated by a sharp interface with differing, but constant, densities and viscosities. We begin by considering the fluids to be miscible with no surface tension acting at the interface. In each fluid layer $j=1,2$ , the governing equation (2.8) simplifies as

(3.41) $$\begin{eqnarray}\displaystyle & & \displaystyle \text{i}\unicode[STIX]{x1D714}\left\{\left(\unicode[STIX]{x1D719}_{j}^{\prime \prime }+\frac{\unicode[STIX]{x1D719}_{j}^{\prime }}{r}-\frac{m^{2}\unicode[STIX]{x1D719}_{j}}{r^{2}}\right)\right\}\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{1}{Re}\frac{\unicode[STIX]{x1D707}_{j}^{\ast }}{\unicode[STIX]{x1D70C}_{j}^{\ast }}\left\{\unicode[STIX]{x1D719}_{j}^{\prime \prime \prime \prime }+\frac{2\unicode[STIX]{x1D719}_{j}^{\prime \prime \prime }}{r}-(1+2m^{2})\left[\frac{\unicode[STIX]{x1D719}_{j}^{\prime \prime }}{r^{2}}-\frac{\unicode[STIX]{x1D719}_{j}^{\prime }}{r^{3}}\right]+\frac{m^{2}(m^{2}-4)\unicode[STIX]{x1D719}_{j}}{r^{4}}\right\}.\end{eqnarray}$$

This can be rewritten in terms of the differential operator, ${\mathcal{L}}$ , defined in (3.2) as

(3.42) $$\begin{eqnarray}\displaystyle \text{i}\unicode[STIX]{x1D714}{\mathcal{L}}[\unicode[STIX]{x1D719}_{j}]=\frac{1}{Re}\frac{\unicode[STIX]{x1D707}_{j}^{\ast }}{\unicode[STIX]{x1D70C}_{j}^{\ast }}{\mathcal{L}}^{2}[\unicode[STIX]{x1D719}_{j}]. & & \displaystyle\end{eqnarray}$$

This fourth-order linear ordinary differential equation can be seen to factorise and so we find solutions are given by solutions of the two second-order ordinary differential equations

(3.43a,b ) $$\begin{eqnarray}\displaystyle {\mathcal{L}}[\unicode[STIX]{x1D719}_{j}]=0,\quad \text{and}\quad \frac{1}{Re}\frac{\unicode[STIX]{x1D707}_{j}^{\ast }}{\unicode[STIX]{x1D70C}_{j}^{\ast }}{\mathcal{L}}[\unicode[STIX]{x1D719}_{j}]=\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D719}_{j}. & & \displaystyle\end{eqnarray}$$

Hence, the general solution to (3.41) is given by

(3.44) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{j}=c_{j1}r^{-m}+c_{j2}r^{m}+c_{j3}\text{J}_{m}(q_{j}\unicode[STIX]{x1D714}^{1/2}r)+c_{j4}\text{Y}_{m}(q_{j}\unicode[STIX]{x1D714}^{1/2}r), & & \displaystyle\end{eqnarray}$$

where $j=1,2$ corresponds to the inner and outer fluid layers respectively, $\text{Y}$ is a Bessel function of the second kind and $q_{j}=(1-\text{i})\sqrt{Re\,\unicode[STIX]{x1D70C}_{j}/2\unicode[STIX]{x1D707}_{j}}$ . The pressure perturbation is given by

(3.45) $$\begin{eqnarray}\displaystyle P_{j}(r)=-\unicode[STIX]{x1D70C}_{j}\left(2\unicode[STIX]{x1D719}_{j}-\frac{\unicode[STIX]{x1D714}r}{m}\unicode[STIX]{x1D719}_{j}^{\prime }\right)+\frac{\unicode[STIX]{x1D707}_{j}}{Re}\frac{\text{i}r}{m}\left\{\unicode[STIX]{x1D719}_{j}^{\prime \prime \prime }+\frac{\unicode[STIX]{x1D719}_{j}^{\prime \prime }}{r}-\frac{(1+m^{2})\unicode[STIX]{x1D719}_{j}^{\prime }}{r^{2}}+\frac{2m^{2}\unicode[STIX]{x1D719}_{j}}{r^{3}}\right\}. & & \displaystyle\end{eqnarray}$$

Velocity regularity at $r=0$ forces $c_{11}=0$ and $c_{14}=0$ . The no-slip and no-penetration conditions at $r=1$ can be enforced by taking

(3.46a ) $$\begin{eqnarray}\displaystyle & \displaystyle c_{21}=-\frac{q_{2}\sqrt{\unicode[STIX]{x1D714}}}{2m}(\text{J}_{m+1}(q_{2}\sqrt{\unicode[STIX]{x1D714}})c_{23}+\text{Y}_{m+1}(q_{2}\sqrt{\unicode[STIX]{x1D714}})c_{24}), & \displaystyle\end{eqnarray}$$

(3.46b ) $$\begin{eqnarray}\displaystyle & \displaystyle c_{22}=-(\text{J}_{m}(q_{2}\sqrt{\unicode[STIX]{x1D714}})c_{23}+\text{Y}_{m}(q_{2}\sqrt{\unicode[STIX]{x1D714}})c_{24})+\frac{q_{2}\sqrt{\unicode[STIX]{x1D714}}}{2m}(\text{J}_{m+1}(q_{2}\sqrt{\unicode[STIX]{x1D714}})c_{23}+\text{Y}_{m+1}(q_{2}\sqrt{\unicode[STIX]{x1D714}})c_{24}). & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Substituting into the stress continuity conditions (2.23) using (3.45), the interfacial jump conditions are, in the absence of surface tension, in the tangential and normal directions respectively

(3.47a ) $$\begin{eqnarray}\displaystyle & \displaystyle [\unicode[STIX]{x1D707}^{\ast }(r_{0}^{2}\unicode[STIX]{x1D719}^{\prime \prime }-r_{0}\unicode[STIX]{x1D719}^{\prime }+m^{2}\unicode[STIX]{x1D719})]_{-}^{+}=0, & \displaystyle\end{eqnarray}$$

(3.47b ) $$\begin{eqnarray}\displaystyle & \displaystyle \left[\unicode[STIX]{x1D70C}^{\ast }\left(r_{0}-2\unicode[STIX]{x1D719}+\frac{\unicode[STIX]{x1D714}r_{0}}{m}\unicode[STIX]{x1D719}^{\prime }\right)+\frac{\unicode[STIX]{x1D707}^{\ast }}{Re}\frac{\text{i}r_{0}}{m}\left\{\unicode[STIX]{x1D719}^{\prime \prime \prime }-\frac{3m^{2}}{r_{0}^{3}}(r_{0}\unicode[STIX]{x1D719}^{\prime }-\unicode[STIX]{x1D719})\right\}\right]_{-}^{+}=0, & \displaystyle\end{eqnarray}$$

In the special case that the dynamic viscosity of each layer is equal, then $\unicode[STIX]{x1D707}_{j}=1$ for $j=1,2$ and we can make use of the continuity of fluid velocity condition, that $\unicode[STIX]{x1D719}_{j}$ and $\unicode[STIX]{x1D719}_{j}^{\prime }$ are continuous across the interface, to simplify (3.47) as

(3.48a ) $$\begin{eqnarray}\displaystyle & \displaystyle [\unicode[STIX]{x1D719}^{\prime \prime }]_{-}^{+}=0, & \displaystyle\end{eqnarray}$$

(3.48b ) $$\begin{eqnarray}\displaystyle & \displaystyle \left[\unicode[STIX]{x1D70C}^{\ast }\{m+\unicode[STIX]{x1D714}(\unicode[STIX]{x1D719}^{\prime }-2)\}+\frac{\text{i}}{Re}\unicode[STIX]{x1D719}^{\prime \prime \prime }\right]_{-}^{+}=0. & \displaystyle\end{eqnarray}$$

(3.49a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{1}(r_{0})=\frac{\unicode[STIX]{x1D714}r_{0}}{m}, & \displaystyle\end{eqnarray}$$

(3.49b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{2}(r_{0})=\frac{\unicode[STIX]{x1D714}r_{0}}{m}, & \displaystyle\end{eqnarray}$$

(3.49c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{1}^{\prime }(r_{0})=\unicode[STIX]{x1D719}_{2}^{\prime }(r_{0}), & \displaystyle\end{eqnarray}$$

$c_{12}$

$c_{13}$

$c_{21}$

$c_{22}$

$\unicode[STIX]{x1D714}$

$\unicode[STIX]{x1D719}_{1}(r_{0})$

$\unicode[STIX]{x1D719}_{1}^{\prime }(r_{0})$

$\unicode[STIX]{x1D719}_{2}(r_{0})$

$\unicode[STIX]{x1D719}_{2}^{\prime }(r_{0})$

Figure 4.

The parameters for the simulations shown are $\mathscr{A}=5/6$ , $m=20$ , $r_{0}=2^{-1/2}$ and where viscous $Re=2\unicode[STIX]{x03C0}\times 10^{6}$ . (a) The predicted normalised amplitude of a perturbation to the interface against time shown for the viscous analytical prediction (bold solid), the inviscid theoretical prediction (thin solid) and the numerical simulation (white data points). The amplitude comprises the superposition of two waves, one travelling in the positive $\unicode[STIX]{x1D703}$ direction and the other in the negative $\unicode[STIX]{x1D703}$ direction. (b) The plot shows how the eigenvalues of the two modes of solution are modified as the viscosity of the two fluid layers is increased. The inviscid solution ( $\unicode[STIX]{x1D714}_{\infty }^{\pm }$ , white data points) is given by (3.7). As the Reynolds number decreases the positive solution approaches zero, whereas the negative solution asymptotes to a finite real value (dashed vertical line) and infinite imaginary value. The real value of the asymptote is given by (A 4). (c) The behaviour of the negative solution as the Reynolds number is varied. In the small Reynolds number limit the solution behaves as $\unicode[STIX]{x1D714}\sim \unicode[STIX]{x1D706}+\text{i}\unicode[STIX]{x1D705}Re^{-1}$ . The real part of the solution is shown in bold and the imaginary part is shown as a thin line. The power-law behaviour of the solutions is indicated by the dashed lines. (d) The behaviour of the positive solution as the Reynolds number is varied. In the small Reynolds number limit the solution behaves as $\unicode[STIX]{x1D714}\sim \text{i}\unicode[STIX]{x1D705}Re$ (see (3.50)).

Figure 4 shows the behaviour of stable configurations of viscous fluid layers of uniform but differing density and a sharp interface. The viscosity of the fluid layers is uniform and equal. For the plots shown the Atwood number is $\mathscr{A}=5/6$ , the azimuthal wavenumber is $m=20$ , the initial interface position is $r_{0}=2^{-1/2}$ , giving equal fluid volume in each layer. Figure 4(a) is a plot of the amplitude of the interfacial disturbance that propagates around the interface with time. The white circles are taken from a low viscosity numerical simulation with $Re=2\unicode[STIX]{x03C0}\times 10^{6}$ . The data points are well modelled by the inviscid solution (3.11) which comprises a superposition of two counter-propagating interfacial waves (thin solid line). Closer agreement is found by using the full viscous solution that better captures the slow decay in amplitude of the perturbations with time (thick solid line). Figure 4(b) shows the behaviour of the eigenvalue $\unicode[STIX]{x1D714}$ as the viscosity of the fluid layers varies. The inviscid solutions, $\unicode[STIX]{x1D714}_{\infty }^{\pm }$ , given by (3.11), are shown as white circles. The positive solution $\unicode[STIX]{x1D714}_{\infty }^{+}\approx 3.33$ corresponds to a wave that precesses in a clockwise direction about the interface whereas the negative solution $\unicode[STIX]{x1D714}_{\infty }^{-}\approx -5.00$ corresponds to a wave that precesses in an anticlockwise direction about the interface. As the viscosity of the fluid layers is increased the precession rates of the solutions decreases, the waves are less able to propagate around the interface. The positive solution tends to zero as $Re\rightarrow 0$ representing a stationary perturbation. The negative solution degenerates as its imaginary part tends to positive infinity, meaning the solution decays instantaneously.

We can examine the behaviour of the viscous solutions in the limit of small Reynolds number. The root associated with $\unicode[STIX]{x1D714}_{\infty }^{+}$ behaves as

(3.50) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714} & {\sim} & \displaystyle \frac{\text{i}\mathscr{A}m}{2}\left\{1-\frac{r_{0}^{2(m-1)}}{2}[(m^{2}(1-r_{0}^{2})^{2}+2r_{0}^{2})(1-\unicode[STIX]{x1D702})+m(1-r_{0}^{4})(1+\unicode[STIX]{x1D702})]-\unicode[STIX]{x1D702}r_{0}^{4m}\right\}\nonumber\\ \displaystyle & & \displaystyle \times \,\frac{r_{0}^{2}}{m^{2}-1}\{1+\unicode[STIX]{x1D702}r_{0}^{2(m-1)}(m^{2}(1-r_{0}^{2})^{2}+2r_{0}^{2})+\unicode[STIX]{x1D702}^{2}r_{0}^{4m}\}^{-1}Re+O(Re^{2}),\end{eqnarray}$$

where $\unicode[STIX]{x1D702}=(\unicode[STIX]{x1D707}_{2}-\unicode[STIX]{x1D707}_{1})/(\unicode[STIX]{x1D707}_{2}+\unicode[STIX]{x1D707}_{1})$ is the viscosity contrast (cf. (18) Alvarez-Lacalle et al. Reference Alvarez-Lacalle, Ortín and Casademunt2004, noting that they have labelled their fluid layers the other way round to the present authors). We can see that this mode is stable for positive Atwood number, $\mathscr{A}>0$ , and unstable for negative Atwood number, $\mathscr{A}<0$ . When the Atwood number is negative, the magnitude of the growth rate is determined largely by the system Reynolds number, $Re$ , the Atwood number, $\mathscr{A}$ , and the azimuthal wavenumber, $m$ . The dependence on the viscosity contrast $\unicode[STIX]{x1D702}$ is multiplied by terms involving $r_{0}$ to powers of $2(m-1)$ or higher. As $r_{0}\in (0,1)$ this indicates therefore, for moderate values of $m$ , the dependence of the growth rate on $\unicode[STIX]{x1D702}$ may be weak and the growth rate for most modes of instability may depend largely only the viscosity of the most viscous layer, but it does not matter whether this layer is the inner layer or the outer layer.

For the negative root associated with $\unicode[STIX]{x1D714}_{\infty }^{-}$ , we can show $\unicode[STIX]{x1D714}\sim \unicode[STIX]{x1D706}+\text{i}\unicode[STIX]{x1D705}Re^{-1}$ where the expressions for the coefficients $\unicode[STIX]{x1D705}$ and $\unicode[STIX]{x1D706}$ are unwieldy, but are given in appendix A. Figures 4(c) and 4(d) show the asymptotic behaviours of the real and imaginary parts of the eigenvalue $\unicode[STIX]{x1D714}$ as the Reynolds number is varied. Both solutions approach the inviscid solution like $Re^{-1/2}$ as $Re\rightarrow \infty$ . In the high viscosity limit as $Re\rightarrow 0$ the positive solution tends to zero as shown in (3.50), the negative solution tends to a constant real part and a singular imaginary part as described in appendix A.

In order for the system at low Reynolds number to support Saffman–Taylor instability (Saffman & Taylor Reference Saffman and Taylor1958) we would require there to exist an $\unicode[STIX]{x1D702}=\unicode[STIX]{x1D702}^{\ast }\in (-1,1)$ such that $\unicode[STIX]{x1D714}=0$ in (3.50), the neutral stability case, i.e. there would exist for $m\geqslant 2$ an $\unicode[STIX]{x1D702}^{\ast }$ that satisfies

(3.51) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702}^{\ast }=\frac{2-r_{0}^{2(m-1)}[m^{2}(1-r_{0}^{2})^{2}+m(1-r_{0}^{4})+2r_{0}^{2}]}{2r_{0}^{4m}-r_{0}^{2(m-1)}[m^{2}(1-r_{0}^{2})^{2}-m(1-r_{0}^{4})+2r_{0}^{2}]}\in [-1,1], & & \displaystyle\end{eqnarray}$$

however no such $\unicode[STIX]{x1D702}^{\ast }$ exists, and so the system does not support pure Saffman–Taylor instability at low Reynolds number.

Figure 5.

The plots show the behaviour of the eigenvalue $\unicode[STIX]{x1D714}$ with changes in the system Reynolds number, $Re$ , for $\mathscr{A}=-5/6$ , $m=20$ , $r_{0}=2^{-1/2}$ . The variation with the viscosity contrast is also shown: the bold solid lines are $\unicode[STIX]{x1D702}=0$ , no contrast, matched viscosities; the thin solid lines are $\unicode[STIX]{x1D702}=5/6$ , the outer layer is 11 times more viscous than the inner layer; the dashed lines are $\unicode[STIX]{x1D702}=-5/6$ , the inner layer is 11 times more viscous than the outer layer. The stable root, left, has the same asymptotic behaviour as the negative root in figure 4(b), diverging to a positive infinite imaginary part as $Re^{-1}$ in the limit $Re\rightarrow 0$ . In the limit of large Reynolds number the solution tends toward the inviscid solution shown by the white circle. The unstable solution tends to zero as $Re\rightarrow 0$ and tends to the inviscid solution as $Re\rightarrow \infty$ , shown by the white circle. The asymptotic expression (3.50) agrees well with the numerical solution for near the $Re\rightarrow 0$ limit. While the stable branch is sensitive to the viscosity contrast, $\unicode[STIX]{x1D702}$ , the unstable branch is insensitive for all $Re$ ; the differences between the three lines cannot be observed easily at the scale shown.

We now turn our attention to viscous fluid layers separated by an initially sharp interface with an unstable density stratification. We consider the unstable complement of figure 4 where, as before, $m=20$ , $r_{0}=2^{-1/2}$ , but we reverse the sign of the Atwood number such that the system is unstable and $\mathscr{A}=-5/6$ . Figure 5 shows the behaviour of solutions to the system as the system Reynolds number is varied for three different values of the viscosity contrast, $\unicode[STIX]{x1D702}$ . The inviscid solutions (3.11) are shown as white circles and it can be seen that as $Re\rightarrow \infty$ all the solutions tend toward their inviscid limit. The numerical solutions tend toward the inviscid solution as $Re^{-1/2}$ as in figure 4. The stable, decaying, solutions are shown on the left and have the same behaviour asymptotically as $Re\rightarrow 0$ as the negative solution in figure 4, specifically that $\unicode[STIX]{x1D714}\sim \unicode[STIX]{x1D706}+\text{i}\unicode[STIX]{x1D705}Re^{-1}$ . The unstable solutions determine the growth of the perturbation at the interface and can be seen to tend to zero, the stationary solution, as $Re\rightarrow 0$ . The solution near $Re\rightarrow 0$ is well approximated by (3.50). For the solutions shown the viscosity contrasts were: $\unicode[STIX]{x1D702}=0$ , for the bold solid lines; $\unicode[STIX]{x1D702}=5/6$ for the thin solid lines, and $\unicode[STIX]{x1D702}=-5/6$ for the dashed lines. A contrast of $\unicode[STIX]{x1D702}=5/6$ corresponds to the outer layer being 11 times more viscous than the inner layer, whereas a contrast of $\unicode[STIX]{x1D702}=-5/6$ corresponds to the inner layer being 11 times more viscous than the outer layer. It can be seen that while the stable solution is sensitive to the viscosity contrast, the unstable solution, which controls the Rayleigh–Taylor growth of the perturbation, is insensitive to the value of $\unicode[STIX]{x1D702}$ across all values of $Re$ . This supports our earlier conjecture that the growth rate is determined primarily by the system Reynolds number and it does not matter whether the outer layer is more viscous than the inner layer or vice versa. In all cases we observe that as the system Reynolds number increases the growth rate of the perturbation is reduced.

Figure 6 shows snapshots of a numerical simulation of two fluids of equal volume and equal viscosity, where $\mathscr{A}=-1/2$ , $r_{0}=2^{-1/2}$ , $Re=2\unicode[STIX]{x03C0}\times 10^{2}$ , $\unicode[STIX]{x1D702}=0$ and the initial perturbation has azimuthal wavenumber $m=45$ . The initial amplitude of the perturbation was $\unicode[STIX]{x1D716}=4\times 10^{-3}$ . The four snapshots are from $t=0.57$ to $t=6.22$ , approximately one complete revolution of the system. The boundary of the domain is shown as a thick solid line and the initial position of the interface is the inner dashed circle in each image. The images are in the rotating frame of reference. The large dotted cross with a white circle at one end indicates the position of a fixed point in the ‘laboratory’, or fixed, inertial frame. This fixed point appears to move in a clockwise direction. As may be observed, the initial perturbation starts to grow. It can be seen in the first two images that the growth is dominated by motion in the radial direction, the effect of the Coriolis term is seen to be small at these times, again the growth of the instability is driven by centrifugal forces.

Figure 6.

Simulation of a centrifugally forced Rayleigh–Taylor instability. The parameters were $Re=2\unicode[STIX]{x03C0}\times 10^{2}$ , $\mathscr{A}=-1/2$ , $r_{0}=2^{-1/2}$ (shown dashed), $\unicode[STIX]{x1D702}=0$ and the initial perturbation was $m=45$ , $\unicode[STIX]{x1D716}=4\times 10^{-3}$ . The times shown are: (a) $t=0.57$ , (b)  $t=2.45$ , (c) $t=4.34$ and (d) $t=6.22$ .

Figure 7 is a comparison of various growth rates calculated from numerical simulations of the centrifugally driven Rayleigh–Taylor instability against the linear stability analysis predictions. For parameters $\mathscr{A}=-1/2$ , $r_{0}=2^{-1/2}$ , $m=45$ , $\unicode[STIX]{x1D702}=0$ and $Re=2\unicode[STIX]{x03C0}\times 10^{6}$ , the growth rate is indicated by the square data points. To a good approximation the data lie at early times on the straight line given by the inviscid approximation $\unicode[STIX]{x1D714}=0.50-4.72\text{i}$ . For a simulation at more moderate Reynolds numbers we expect the inviscid approximation to be a poor estimate of the growth rate. The white circle data points are from a simulation with $Re=2\unicode[STIX]{x03C0}\times 10^{2}$ . We may not neglect the effects of viscosity at this Reynolds number, as the growth rate is seen to be substantially lower than that of the high Reynolds number simulation. The theoretical prediction for this configuration is that $\unicode[STIX]{x1D714}=0.09-1.49\text{i}$ , and can be seen to match well the numerical data points (white circle). The remaining two sets of data are for configurations with a diffuse interface and with surface tension which are discussed in § 3.5.

Figure 7.

A comparison of the growth rates observed in numerical simulations compared to the theoretical predictions. The parameters for all simulations were $\mathscr{A}=-1/2$ , $m=45$ , $r_{0}=2^{-1/2}$ . The individual simulation runs had: (▫) low viscosity $Re=2\unicode[STIX]{x03C0}\times 10^{6}$ , no surface tension, no diffusion; (

) low viscosity $Re=2\unicode[STIX]{x03C0}\times 10^{6}$ , no surface tension, an initial diffuse interface with width $\unicode[STIX]{x1D6FF}=0.067$ ; (○) high viscosity $Re=2\unicode[STIX]{x03C0}\times 10^{2}$ , no surface tension, no diffusion; (▵) low viscosity $Re=2\unicode[STIX]{x03C0}\times 10^{6}$ , surface tension $We=5\times 10^{3}$ , no diffusion. The straight thin lines are the theoretically predicted growth rates.

3.5 Immiscible viscous layers with surface tension and miscible diffuse viscous fluid layers

Two classes of viscous flow are of practical interest. We may have a configuration where the two fluids are different, immiscible with different viscosities and surface tension between them. We may also consider the case of miscible fluids that have approximately equal viscosities, no surface tension between the layers, but a diffuse interface.

We first consider the case of two immiscible fluids with a constant density and constant viscosity contrast between them and surface tension acting at the interface. The governing version of the Orr–Sommerfeld equation (2.8) is as in § 3.4 and is given by (3.41), which accepts the same solution as before. However, the boundary conditions at the interface must be modified to account for the stress jump due to the surface tension, similarly to the inviscid case considered in § 3.2. The kinematic condition (3.49a ) and (3.49b ), the tangential velocity continuity condition (3.49c ) and the tangential stress condition (3.47a ) remain as before. However (3.47b ) is modified, following (2.23) to give

(3.52) $$\begin{eqnarray}\displaystyle \left[\unicode[STIX]{x1D70C}^{\ast }\left(r_{0}-2\unicode[STIX]{x1D719}+\frac{\unicode[STIX]{x1D714}r_{0}}{m}\unicode[STIX]{x1D719}\right)+\frac{\unicode[STIX]{x1D707}^{\ast }}{Re}\frac{\text{i}r}{m}\left\{\unicode[STIX]{x1D719}^{\prime \prime \prime }-\frac{3m^{2}}{r_{0}^{3}}(r_{0}\unicode[STIX]{x1D719}^{\prime }-\unicode[STIX]{x1D719})\right\}\right]_{-}^{+}=-\frac{1}{We}\frac{m^{2}-1}{r_{0}^{2}} & & \displaystyle\end{eqnarray}$$

(cf. the inviscid condition (3.15)) where (3.47a ) has been used to simplify (3.52). As in the miscible case, (3.49) and (3.47a ) may be used to find expressions for $c_{12}$ , $c_{13}$ , $c_{21}$ and $c_{22}$ in terms of $\unicode[STIX]{x1D714}$ and also then $\unicode[STIX]{x1D719}_{1}$ and $\unicode[STIX]{x1D719}_{2}$ and their gradients at $r_{0}$ . The final condition (3.52) yields the dispersion relation. The inclusion of the effect of surface tension in the right-hand side of (3.52) on the low Reynolds number asymptotic solution (3.50) is to introduce a correcting factor $1+(m^{2}-1)/(2We\,\mathscr{A}r_{0}^{3})=1+S/m$ into the expression for $\unicode[STIX]{x1D714}$ . This has the corollary that the surface tension may stabilise the mode with azimuthal wavenumber $m^{\ast }$ when $S=-m^{\ast }$ , or, equivalently, all modes $m>m^{\ast }$ are stabilised for $We=-(m^{\ast 2}-1)/(2\mathscr{A}r_{0}^{3})$ .

Figure 8.

The growth of modes of disturbance at an interface without surface tension (a,c) and with surface tension (b,d). Snap shots from the flows are shown in images (a) and (b) at $t=2.51$ . The surface tension in (b) and (d) corresponds to a cutoff azimuthal wavenumber $m^{\ast }=45$ for growth as given by (3.18), indicated by the bold vertical line. In (c) and (d) it can be seen that each individual mode evolves in time, indicated by circular data points, filled white at $t=0$ and getting darker with each time step $\unicode[STIX]{x0394}t=0.063$ . The predicted magnitude of each mode based on either (3.7) (c) or (3.17) (d), at the final time (black circle) is indicated by a white square. The final times are $t=0.44$ in (c) and $t=0.82$ in (d). It can be seen in (d) that modes that are above the cutoff wavenumber $m^{\ast }=45$ are suppressed at these early times, to a good approximation. The parameters for the flow were $\mathscr{A}=-0.5$ , and $We=5.8\times 10^{3}$ in (b) and (d).

Figure 8 shows a comparison between the behaviour of two viscous systems, one without surface tension at the interface (a), and one with surface tension at the interface (b). Both systems have equal fluid viscosities in each layer, i.e. $\unicode[STIX]{x1D702}=0$ , and an unstable density stratification, $\mathscr{A}=-1/2$ . The Reynolds number, as defined in § 2.1, in each simulation was $2\unicode[STIX]{x03C0}\times 10^{6}$ . The Reynolds number based on the speed of propagation of the perturbed front was less than $5\times 10^{3}$ , two orders of magnitude below that identified by Zhou (Reference Zhou2007) for the transition to the minimum state required for comparison to astrophysical-scale flows.

At $t=0$ a number of modes of perturbation between $m=5$ and $m=90$ were introduced at the interface with random amplitudes of order $4\times 10^{-3}$ . This initial perturbation is indicated by the white circles on figures 8(c) and 8(d). As time evolved and the instability developed, the amplitude of each mode was added to figures 8(c) and 8(d) as a progressively darker shaded data point for a given $m$ at time steps of $\unicode[STIX]{x0394}t=0.063$ . At the last time shown ( $t=0.44$ in (c) and $t=0.82$ in (d)) the last data point (black circle) is compared with its theoretical prediction, indicated by a white square. The simulation on the left (figure 8 a) had no surface tension acting at the interface, whereas the simulation on the right (figure 8 b) had surface tension acting at the interface and the effect can be immediately observed qualitatively. The small-scale instabilities apparent in the left hand simulation appear significantly suppressed in the right-hand simulation. The strength of the surface tension, $We=5800$ , was chosen such that the critical wavenumber, $m^{\ast }$ , given in (3.18) is approximately 45. Figure 8(d) shows that the behaviour of modes above $m=m^{\ast }$ is completely different to that below $m=m^{\ast }$ and their growth is suppressed at early times.

As the instability develops, some of the dense fluid moves towards the boundary of the domain, as is observed from the simulations. As the interface between the two fluids moves towards the boundary it may develop areas where its radius of curvature is greater than its initial value at $r=r_{0}$ . This therefore allows modes that were suppressed initially by the surface tension when the interface was nearer the centre of the system to develop as the interface approaches the boundary since the effects of surface tension are not felt as strongly in regions of lower interfacial curvature. We can return to figure 7 and consider the growth of a system where initially the perturbed mode is $m=45$ , but the surface tension $We=5\times 10^{3}$ should completely suppress the growth of the mode. The triangular data points show that initially the $m=45$ mode is unable to grow, but later as the interface has moved toward the boundary, due to instability at lower wavenumbers, the $m=45$ mode is able to grow.

We now turn to the case of two layers of miscible viscous fluid with matching viscosities whose interface has diffused over a period of time. The approach is identical to that followed in the inviscid case where we consider a diffuse layer thickness $\unicode[STIX]{x1D6FF}\sim 2\sqrt{\unicode[STIX]{x03C0}t/Pe}$ and a piecewise continuous density of the form (3.29). Assuming equal viscosities for the fluids we define the following quantity

(3.53) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D701}(r)=\frac{Re\,\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D714}\text{i}}{(2+\unicode[STIX]{x1D6FC})^{2}}r^{2+\unicode[STIX]{x1D6FC}}, & & \displaystyle\end{eqnarray}$$

and functions

(3.54) $$\begin{eqnarray}\displaystyle {\mathcal{F}}(c;\unicode[STIX]{x1D701}) & = & \displaystyle \text{}_{2}F_{3}\left(\left[\frac{\unicode[STIX]{x1D6FC}+2(c+\unicode[STIX]{x1D712})}{2(2+\unicode[STIX]{x1D6FC})},\frac{\unicode[STIX]{x1D6FC}+2(c-\unicode[STIX]{x1D712})}{2(2+\unicode[STIX]{x1D6FC})}\right]\!;\left[\frac{\unicode[STIX]{x1D6FC}}{2+\unicode[STIX]{x1D6FC}},\frac{\unicode[STIX]{x1D6FC}+2c}{2+\unicode[STIX]{x1D6FC}},\frac{\unicode[STIX]{x1D6FC}+2(c+1)}{2+\unicode[STIX]{x1D6FC}}\right]\!;\unicode[STIX]{x1D701}\right),\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(3.55) $$\begin{eqnarray}\displaystyle {\mathcal{G}}(c;\unicode[STIX]{x1D701}) & = & \displaystyle \text{}_{2}F_{3}\left(\left[\frac{\unicode[STIX]{x1D6FC}+2(c+\unicode[STIX]{x1D712})}{2(2+\unicode[STIX]{x1D6FC})},\frac{\unicode[STIX]{x1D6FC}+2(c-\unicode[STIX]{x1D712})}{2(2+\unicode[STIX]{x1D6FC})}\right]\!;\left[\frac{4+\unicode[STIX]{x1D6FC}}{2+\unicode[STIX]{x1D6FC}},\frac{\unicode[STIX]{x1D6FC}+2(c-1)}{2+\unicode[STIX]{x1D6FC}},\frac{\unicode[STIX]{x1D6FC}+2c}{2+\unicode[STIX]{x1D6FC}}\right]\!;\unicode[STIX]{x1D701}\right),\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $\text{}_{2}F_{3}$ is a hypergeometric function and $\unicode[STIX]{x1D712}$ is as defined in (3.36). We may write solutions for each region as

(3.56) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{1}=c_{11}r^{-m}+c_{12}r^{m}+c_{13}\text{J}_{m}(q_{1}\unicode[STIX]{x1D714}^{1/2}r)+c_{14}\text{Y}_{m}(q_{1}\unicode[STIX]{x1D714}^{1/2}r), & \displaystyle\end{eqnarray}$$

(3.57) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{2}=c_{21}r^{-m}{\mathcal{F}}(-m;\unicode[STIX]{x1D701})+c_{22}r^{m}{\mathcal{F}}(m;\unicode[STIX]{x1D701})+c_{23}r^{2-m}{\mathcal{G}}(2-m,\unicode[STIX]{x1D701})+c_{24}r^{2+m}{\mathcal{G}}(2+m,\unicode[STIX]{x1D701}), & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(3.58) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{3}=c_{31}r^{-m}+c_{32}r^{m}+c_{33}\text{J}_{m}(q_{2}\unicode[STIX]{x1D714}^{1/2}r)+c_{34}\text{Y}_{m}(q_{2}\unicode[STIX]{x1D714}^{1/2}r). & \displaystyle\end{eqnarray}$$

We then proceed as before, enforcing velocity regularity at $r=0$ which forces $c_{11}=0$ and $c_{14}=0$ , and no-slip, no-penetration conditions at $r=1$ which forces equivalent constraints to (3.46). This reduces the number of free constants to eight; two in the inner and outer layers, and four in the mid-layer. We now apply the kinematic condition and velocity continuity conditions at each interface (five conditions) and stress continuity at each interface (four conditions) to find the dispersion relation.

Figure 9 shows the effect of a diffuse interface on the real part of the solution $\unicode[STIX]{x1D719}$ . The piecewise solution comprises three segments that, as in the inviscid case, transition at the white circles. The unused parts of the solutions are shown, indicating their behaviour in the matching regions. The exact solution is shown as a thin solid line, the kinematic condition is enforced at $r=r_{0}$ , indicated by the black circle. The solution to (2.8) subject to $\unicode[STIX]{x1D719}(0)=\unicode[STIX]{x1D719}^{\prime }(0)=0$ , $\unicode[STIX]{x1D719}(1)=\unicode[STIX]{x1D719}^{\prime }(1)=0$ was calculated separately in $[0,r_{0}]$ and $[r_{0},1]$ with the kinematic condition enforced for each solution at $r=r_{0}$ . The free boundary condition $\unicode[STIX]{x1D719}^{\prime }(r_{0})$ in each domain and the eigenvalue $\unicode[STIX]{x1D714}$ are chosen to ensure that $\unicode[STIX]{x1D719}$ is at least class $C^{4}[0,1]$ . The method of solution, as in the inviscid case, was to formulate the problem in terms of Chebyshev polynomials. In figure 9(b) the comparison with the non-diffuse case is shown and it can be seen that the qualitative differences between the solutions are less marked than in the inviscid case since the viscosity has acted to smooth $\unicode[STIX]{x1D719}$ to some extent already.

Figure 9.

(a) The solution, $\unicode[STIX]{x1D719}$ , (bold solid line) to (2.8) for $m=12$ , $\mathscr{A}=-1/2$ , $r_{0}=2^{-1/2}$ , $Re=400$ , $\unicode[STIX]{x1D707}_{1}=1$ , $\unicode[STIX]{x1D707}_{2}=1$ and $\unicode[STIX]{x1D6FF}=2\sqrt{\unicode[STIX]{x03C0}t/Pe}\approx 0.11$ , where $t=10^{-3}Pe$ . The shown solution is the real part, and corresponds to the solution of the dispersion relation at $\unicode[STIX]{x1D714}=0.22-1.41\text{i}$ . The solution constitutes the matching of three analytic solutions; $\unicode[STIX]{x1D719}_{1}$ for $r\in [0,r_{0-}]$ which continues dot-dashed after the solution is matched to $\unicode[STIX]{x1D719}_{2}$ for $r\in [r_{0-},r_{0+}]$ . Outside this region $\unicode[STIX]{x1D719}_{2}$ is shown dashed. Finally in $r\in [r_{0+},1]$ , $\unicode[STIX]{x1D719}_{3}$ is used, but is shown dashed outside this range. The diffusion width $\unicode[STIX]{x1D6FF}$ is indicated, and the matching between solutions occurs at the white dots. The exact solution is shown as a thin solid line. The kinematic condition is enforced at the black circle. (b) Comparison with the zero-diffusion case (higher amplitude curve); again the transition between the two solutions is shown with a white dot, the solution corresponds to the solution of the dispersion relation at $\unicode[STIX]{x1D714}=0.29-1.62\text{i}$ .

3.6 Non-inertial flow and comparison with Hele-Shaw cell and porous media flow

A related problem to those considered above is of a two-dimensional rotating droplet, possibly lying within an unbounded fluid of different density and viscosity, rotating in a Hele-Shaw cell or a porous medium. This problem was considered by Schwartz (Reference Schwartz1989) and later by Alvarez-Lacalle et al. (Reference Alvarez-Lacalle, Ortín and Casademunt2004). If the outer fluid is considered unbounded and inertia is to be ignored then the non-dimensionalisation used above is not natural since $a\rightarrow \infty$ , but $\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D6FA}a^{2}/\unicode[STIX]{x1D707}_{0}\rightarrow 0$ . It is more appropriate therefore to take the radius of the ‘droplet’, $r_{0}$ , as the length scale in this case.

We make a Stokes flow approximation, taking $Re\rightarrow 0$ whereby inertial terms are ignored but the centrifugal term and the pressure term are balanced with the viscous term of the equation of motion (2.1b ). This requires a velocity scale $(\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D6FA}r_{0}^{2}/\unicode[STIX]{x1D707}_{0})\unicode[STIX]{x1D6FA}r_{0}$ and hence an implied time scale $\unicode[STIX]{x1D707}_{0}/(\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D6FA}^{2}r_{0}^{2})$ . We also require the pressure scale $\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D6FA}^{2}r_{0}^{2}$ , and define a Reynolds number $Re=(\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D6FA}r_{0}^{2}/\unicode[STIX]{x1D707}_{0})^{2}$ . The non-dimensional equations of motion are, assuming constant viscosity and density in each layer

(3.59a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{j}^{\prime }}{\unicode[STIX]{x2202}t^{\prime }}+\unicode[STIX]{x1D735}^{\prime }\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{j}^{\prime }\boldsymbol{u}_{j}^{\prime })=0, & \displaystyle\end{eqnarray}$$

(3.59b ) $$\begin{eqnarray}\displaystyle & \displaystyle Re\frac{\text{D}\boldsymbol{u}_{j}^{\prime }}{\text{D}t^{\prime }}=-\frac{1}{\unicode[STIX]{x1D70C}_{j}^{\prime }}\unicode[STIX]{x1D735}^{\prime }p_{j}^{\prime }+r^{\prime }\hat{\boldsymbol{r}}-2Re^{1/2}\hat{\boldsymbol{z}}\times \boldsymbol{u}_{j}^{\prime }+\frac{\unicode[STIX]{x1D707}_{j}^{\prime }}{\unicode[STIX]{x1D70C}_{j}^{\prime }}\unicode[STIX]{x1D6FB}^{\prime 2}\boldsymbol{u}_{j}^{\prime }, & \displaystyle\end{eqnarray}$$

(3.59c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}^{\prime }\boldsymbol{\cdot }\boldsymbol{u}_{j}^{\prime }=0. & \displaystyle\end{eqnarray}$$

$Re\rightarrow 0$

(3.60a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{j}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{j}\boldsymbol{u}_{j})=0, & \displaystyle\end{eqnarray}$$

(3.60b ) $$\begin{eqnarray}\displaystyle & \displaystyle \mathbf{0}=-\frac{1}{\unicode[STIX]{x1D70C}_{j}}\unicode[STIX]{x1D735}p_{j}+r\hat{\boldsymbol{r}}+\frac{\unicode[STIX]{x1D707}_{j}}{\unicode[STIX]{x1D70C}_{j}}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{j}, & \displaystyle\end{eqnarray}$$

(3.60c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{j}=0. & \displaystyle\end{eqnarray}$$

(3.61) $$\begin{eqnarray}\displaystyle 0=\unicode[STIX]{x1D719}_{j}^{\prime \prime \prime \prime }+\frac{2\unicode[STIX]{x1D719}_{j}^{\prime \prime \prime }}{r}-(1+2m^{2})\left[\frac{\unicode[STIX]{x1D719}_{j}^{\prime \prime }}{r^{2}}-\frac{\unicode[STIX]{x1D719}_{j}^{\prime }}{r^{3}}\right]+\frac{m^{2}(m^{2}-4)\unicode[STIX]{x1D719}_{j}}{r^{4}}={\mathcal{L}}^{2}[\unicode[STIX]{x1D719}_{j}]. & & \displaystyle\end{eqnarray}$$

This is a simplification of (2.8) in which the left-hand side, the inertial terms are zero and we have uniform viscosity in each layer. The corresponding pressure perturbation is given by

(3.62) $$\begin{eqnarray}\displaystyle P_{j}=\frac{\text{i}\unicode[STIX]{x1D707}_{j}}{m}\left((r\unicode[STIX]{x1D719}_{j}^{\prime \prime })^{\prime }-\frac{m^{2}+1}{r}\unicode[STIX]{x1D719}_{j}^{\prime }+\frac{2m^{2}\unicode[STIX]{x1D719}_{j}}{r^{2}}\right). & & \displaystyle\end{eqnarray}$$

This zero-inertia form of the Orr–Sommerfeld equation accepts power-law solutions of the form

(3.63) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{j}=c_{j1}r^{-m}+c_{j2}r^{-m+2}+c_{j3}r^{m}+c_{j4}r^{m+2},\quad j=1,2. & & \displaystyle\end{eqnarray}$$

The coefficients $c_{11}$ and $c_{12}$ are taken to be zero for velocity regularity at the origin. The no-slip and no-penetration conditions are now applied at $r=a$ such that $\unicode[STIX]{x1D719}_{2}(a)=0$ , $\unicode[STIX]{x1D719}_{2}^{\prime }(a)=0$ . The kinematic condition, which also enforces normal velocity continuity, is given now by $\unicode[STIX]{x1D719}_{1}(1)=\unicode[STIX]{x1D719}_{2}(1)=\unicode[STIX]{x1D714}/m$ . Tangential velocity continuity is ensured by setting $\unicode[STIX]{x1D719}_{1}^{\prime }(1)=\unicode[STIX]{x1D719}_{2}^{\prime }(1)$ . This leaves only one free constant, but the two stress continuity conditions to satisfy and so will yield a dispersion relation. The normal and tangential stress continuity conditions are, at order $\unicode[STIX]{x1D716}$ , respectively

(3.64) $$\begin{eqnarray}\displaystyle & \displaystyle [\unicode[STIX]{x1D707}^{\ast }\{\unicode[STIX]{x1D719}^{\prime \prime \prime }-3m^{2}(\unicode[STIX]{x1D719}^{\prime }-\unicode[STIX]{x1D719})\}-\text{i}m\unicode[STIX]{x1D70C}^{\ast }]_{-}^{+}=\frac{\text{i}m(m^{2}-1)}{We}, & \displaystyle\end{eqnarray}$$

(3.65) $$\begin{eqnarray}\displaystyle & \displaystyle [\unicode[STIX]{x1D707}^{\ast }\{\unicode[STIX]{x1D719}^{\prime \prime }-\unicode[STIX]{x1D719}^{\prime }+m^{2}\unicode[STIX]{x1D719}\}]_{-}^{+}=0, & \displaystyle\end{eqnarray}$$

where the jump occurs across $r=1$ , and the tangential condition has been used to simplify the normal condition. For $a\gg 1$ , $m>1$ , the dispersion relation is given by

(3.66) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714} & = & \displaystyle \frac{\text{i}\mathscr{A}m}{2}\left(1+\frac{1}{We}\frac{m^{2}-1}{2\mathscr{A}}\right)\nonumber\\ \displaystyle & & \displaystyle \times \left(1-\frac{a^{-2(m+1)}}{2}[(a^{2}-1)^{2}(1-\unicode[STIX]{x1D702})m^{2}+(a^{4}-1)(1+\unicode[STIX]{x1D702})m+2a^{2}(1-\unicode[STIX]{x1D702})]-a^{-4m}\unicode[STIX]{x1D702}\right)\nonumber\\ \displaystyle & & \displaystyle \times \,\{(m^{2}-1)[1+\unicode[STIX]{x1D702}(a^{-2(m+1)}[(a^{2}-1)^{2}m^{2}+2a^{2}]+a^{-4m}\unicode[STIX]{x1D702})]\}^{-1}\end{eqnarray}$$

(cf. (3.50) with the surface tension correction factor). The form of the dispersion relation shows immediately that solutions will either decay or grow, depending on whether the Atwood number is positive or negative respectively, but there are no precessional or travelling wave solutions possible as might be anticipated on physical grounds for a system with no inertia. There can be a balance between the stabilising effect of the surface tension and the Rayleigh–Taylor instability such that a perturbed interface neither grows nor decays and is stationary. This occurs when

(3.67) $$\begin{eqnarray}\displaystyle We=-\frac{m^{2}-1}{2\mathscr{A}}, & & \displaystyle\end{eqnarray}$$

which can only occur for $\mathscr{A}<0$ as would be expected again on physical grounds. If quantities are rewritten in terms of our original non-dimensionalisation (§ 2.1) and $\unicode[STIX]{x1D707}_{1}=\unicode[STIX]{x1D707}_{2}$ , i.e. $\unicode[STIX]{x1D702}=0$ , then the dispersion relation (3.66) is exactly the first term of (3.50) with the surface tension correction factor, as it must be.

We compare the results of Schwartz (Reference Schwartz1989) and Alvarez-Lacalle et al. (Reference Alvarez-Lacalle, Ortín and Casademunt2004) with the related non-inertial limit of the flows considered in §§ 3.1–3.5. In a Hele-Shaw cell the equations of motion are classically simplified under Stokes flow and lubrication approximations that respectively ignore any inertia in the flow and assume that gradients in the gap-width direction are much greater than gradients in the planar direction. For consideration of the flow in a Hele-Shaw cell or a porous media we return to (3.59) but instead of interpreting the equations in plane-polar coordinates we interpret them in cylindrical polar coordinates and apply a lubrication approximation whereby changes in the planar direction take place over much greater distances than changes in the gap-width $z$ -direction. That is, we apply a second scaling to (3.59) under the assumption that the dimensional gap width, $b$ , is small compared to the initial radius of the inner fluid layer. Hence, we take $\unicode[STIX]{x1D700}=b/r_{0}\ll 1$ . Under this rescaling we write the gradient operator and the velocity separately in their planar and vertical components: $\unicode[STIX]{x1D735}=\unicode[STIX]{x1D735}_{H}+\unicode[STIX]{x2202}\boldsymbol{\cdot }/\unicode[STIX]{x2202}z\hat{\boldsymbol{z}}$ , $\boldsymbol{u}=\boldsymbol{u}_{H}+w\hat{\boldsymbol{z}}$ . The pressure, $p$ , is scaled by $\unicode[STIX]{x1D700}^{-2}$ , and the system (3.59) is transformed to

(3.68a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{j}^{\prime }}{\unicode[STIX]{x2202}t^{\prime }}+\unicode[STIX]{x1D735}^{\prime }\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{j}^{\prime }\boldsymbol{u}_{j}^{\prime })=0, & \displaystyle\end{eqnarray}$$

(3.68b ) $$\begin{eqnarray}\displaystyle & \displaystyle Re\,\unicode[STIX]{x1D700}^{2}\frac{\text{D}\boldsymbol{u}_{Hj}^{\prime }}{\text{D}t^{\prime }}=-\frac{1}{\unicode[STIX]{x1D70C}_{j}^{\prime }}\unicode[STIX]{x1D735}_{H}^{\prime }p^{\prime }+r^{\prime }\hat{\boldsymbol{r}}-2Re^{1/2}\unicode[STIX]{x1D700}\,\hat{\boldsymbol{z}}\times \boldsymbol{u}_{Hj}^{\prime }+\frac{\unicode[STIX]{x1D707}_{j}^{\prime }}{\unicode[STIX]{x1D70C}_{j}^{\prime }}(\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D6FB}_{H}^{\prime 2}+\unicode[STIX]{x2202}_{z}^{2})\boldsymbol{u}_{Hj}^{\prime }, & \displaystyle\end{eqnarray}$$

(3.68c ) $$\begin{eqnarray}\displaystyle & \displaystyle Re\,\unicode[STIX]{x1D700}^{2}\frac{\text{D}w_{j}^{\prime }}{\text{D}t^{\prime }}=-\frac{1}{\unicode[STIX]{x1D700}^{2}}\frac{1}{\unicode[STIX]{x1D70C}_{j}^{\prime }}\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}z^{\prime }}+\frac{\unicode[STIX]{x1D707}_{j}^{\prime }}{\unicode[STIX]{x1D70C}_{j}^{\prime }}(\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D6FB}_{H}^{\prime 2}+\unicode[STIX]{x2202}_{z}^{2})w_{j}^{\prime }, & \displaystyle\end{eqnarray}$$

(3.68d ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}^{\prime }\boldsymbol{\cdot }\boldsymbol{u}_{j}^{\prime }=0. & \displaystyle\end{eqnarray}$$

(3.69a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{j}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{j}\boldsymbol{u}_{j})=0, & \displaystyle\end{eqnarray}$$

(3.69b ) $$\begin{eqnarray}\displaystyle & \displaystyle \mathbf{0}=-\frac{1}{\unicode[STIX]{x1D70C}_{j}}\unicode[STIX]{x1D735}_{H}p+r\hat{\boldsymbol{r}}+\frac{\unicode[STIX]{x1D707}_{j}}{\unicode[STIX]{x1D70C}_{j}}\frac{\unicode[STIX]{x2202}^{2}\boldsymbol{u}_{Hj}}{\unicode[STIX]{x2202}z^{2}}, & \displaystyle\end{eqnarray}$$

(3.69c ) $$\begin{eqnarray}\displaystyle & \displaystyle 0=\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}z}, & \displaystyle\end{eqnarray}$$

(3.69d ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{j}=0. & \displaystyle\end{eqnarray}$$

$p=p(r,\unicode[STIX]{x1D703},t)$

$z$

$z=0$

$z=\unicode[STIX]{x1D700}$

(3.70) $$\begin{eqnarray}\displaystyle \boldsymbol{u}_{j}=\frac{z(z-\unicode[STIX]{x1D700})}{2\unicode[STIX]{x1D707}_{j}}\unicode[STIX]{x1D735}_{H}\left(p_{j}-\frac{\unicode[STIX]{x1D70C}_{j}r^{2}}{2}\right), & & \displaystyle\end{eqnarray}$$

and we can define a vertically averaged velocity in each fluid layer

(3.71) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{j}=\frac{1}{\unicode[STIX]{x1D700}}\int _{0}^{\unicode[STIX]{x1D700}}\boldsymbol{u}_{j}\,\text{d}z=-\frac{\unicode[STIX]{x1D700}^{2}}{12\unicode[STIX]{x1D707}_{j}}\unicode[STIX]{x1D735}_{H}\left(p_{j}-\frac{\unicode[STIX]{x1D70C}_{j}r^{2}}{2}\right), & & \displaystyle\end{eqnarray}$$

a rotational version of Darcy’s law. As the pressure in each layer is independent of $z$ it follows from taking the curl of (3.71) that the vertically averaged velocity, $\boldsymbol{v}_{j}$ , is irrotational in each layer and so Alvarez-Lacalle et al. (Reference Alvarez-Lacalle, Ortín and Casademunt2004) proceed by taking a potential for the flow such that $\boldsymbol{v}_{j}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}_{j}$ . Hence, as a result of incompressibility, they seek to solve $\unicode[STIX]{x1D6FB}_{H}^{2}\unicode[STIX]{x1D711}_{j}=0$ in each layer. For consistency with notation in the previous sections we follow an equivalent method which is to express $\boldsymbol{v}_{j}$ in terms of a streamfunction, $\unicode[STIX]{x1D713}_{j}$ , since the velocity field is incompressible. The irrotationality of $\boldsymbol{v}_{j}$ then implies that $\unicode[STIX]{x1D6FB}_{H}^{2}\unicode[STIX]{x1D713}_{j}=0$ . We therefore seek normal mode solutions and solve Laplace’s equation for $\unicode[STIX]{x1D713}_{j}=\unicode[STIX]{x1D716}\unicode[STIX]{x1D719}_{j}(r)\exp \{\text{i}(m\unicode[STIX]{x1D703}+\unicode[STIX]{x1D714}t)\}$ (taking care to distinguish between $\unicode[STIX]{x1D700}$ the non-dimensional gap width of the Hele-Shaw cell, and $\unicode[STIX]{x1D716}$ the small perturbation to the background hydrostatic initial conditions). The pressure is then given by

(3.72) $$\begin{eqnarray}\displaystyle p_{j}=p_{j}^{\ast }-\unicode[STIX]{x1D716}\frac{12\unicode[STIX]{x1D707}_{j}}{\unicode[STIX]{x1D700}^{2}}\frac{\text{i}r}{m}\unicode[STIX]{x1D719}_{j}^{\prime }\,\text{e}^{\text{i}(m\unicode[STIX]{x1D703}+\unicode[STIX]{x1D714}t)},\quad \text{where}~p_{j}^{\ast }=p_{0j}+\frac{\unicode[STIX]{x1D70C}_{j}r^{2}}{2}, & & \displaystyle\end{eqnarray}$$

and $p_{0j}$ is a reference pressure. (Note that solving Laplace’s equation for $\unicode[STIX]{x1D713}$ leads to exactly the same equation for $\unicode[STIX]{x1D719}$ as the Rayleigh equation (3.2) in § 3.1 demonstrating the well-known phenomenon of three-dimensional viscous flow in a Hele-Shaw cell modelling two-dimensional inviscid flow.) The conditions of stress continuity are

(3.73a ) $$\begin{eqnarray}\displaystyle & \displaystyle [p^{\ast }]_{-}^{+}=-\frac{1}{We},\quad O(1)~\text{normal stress}, & \displaystyle\end{eqnarray}$$

(3.73b ) $$\begin{eqnarray}\displaystyle & \displaystyle \left[-\unicode[STIX]{x1D70C}^{\ast }+\frac{2\text{i}\unicode[STIX]{x1D707}^{\ast }}{m}\left\{\!\left(\!\frac{6}{\unicode[STIX]{x1D700}^{2}}+m^{2}\right)\unicode[STIX]{x1D719}^{\prime }-m^{2}\unicode[STIX]{x1D719}\!\right\}\!\right]_{-}^{+}=\frac{1}{We}(m^{2}-1),\quad O(\unicode[STIX]{x1D716})~\text{normal stress}, & \displaystyle\end{eqnarray}$$

(3.73c ) $$\begin{eqnarray}\displaystyle & \displaystyle [\unicode[STIX]{x1D707}^{\ast }(\unicode[STIX]{x1D719}^{\prime \prime }-\unicode[STIX]{x1D719}^{\prime }+m^{2}\unicode[STIX]{x1D719})]_{-}^{+}=0,\quad O(\unicode[STIX]{x1D716})~\text{tangential stress} & \displaystyle\end{eqnarray}$$

$r=1$

$\unicode[STIX]{x1D719}_{j}(1)=\unicode[STIX]{x1D714}/m$

$\unicode[STIX]{x1D719}_{1}^{\prime }(1)=\unicode[STIX]{x1D719}_{2}^{\prime }(1)$

The solutions to Laplace’s equation have only two free constants in each layer, and once the conditions of regularity at the origin and unboundedness in the outer layer are taken into account, only two free constants remain to satisfy the tangential velocity continuity condition, stress continuity conditions and kinematic condition. Alvarez-Lacalle et al. (Reference Alvarez-Lacalle, Ortín and Casademunt2004) enforce the kinematic condition, and thus normal velocity continuity, and normal stress continuity to find (in our notation)

(3.74) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{1}=\frac{\unicode[STIX]{x1D714}r^{m}}{m},\quad \unicode[STIX]{x1D719}_{2}=\frac{\unicode[STIX]{x1D714}r^{-m}}{m}, & & \displaystyle\end{eqnarray}$$

and hence find the dispersion relation

(3.75) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}=\frac{\text{i}\mathscr{A}m}{2}\left(1+\frac{1}{We}\frac{m^{2}-1}{2\mathscr{A}}\right)\left[\frac{6}{\unicode[STIX]{x1D700}^{2}}+m(m+\unicode[STIX]{x1D702})\right]^{-1}. & & \displaystyle\end{eqnarray}$$

We note that as in the plane-polar two-dimensional case, there can be a balance between the unstable density stratification and surface tension that leads to a stationary perturbation when $We=-(m^{2}-1)/2\mathscr{A}$ . The conditions of tangential stress and velocity continuity, are not enforced.

Figure 10.

Modified pressure contours in perturbed droplets rotating in porous media or Hele-Shaw cells. The initial droplet size is a perturbed circle of initially unit radius, and the cross is for scale (length 0.4 in each direction). The times shown are $t=15$ , 25, 30 and 35. There is no external fluid layer.

Schwartz (Reference Schwartz1989) considered the same problem in a porous media as well as in a Hele-Shaw cell. Figure 10 shows droplets with differing initial perturbations in porous media evolving in time, governed by Schwartz’s model; contours of $p-\unicode[STIX]{x1D70C}r^{2}/2$ are shown. Flows in Hele-Shaw cells of gap width $b$ are equivalent to flows in porous media, governed by Darcy’s law, with permeability $k=b^{2}/12$ . Therefore for a fluid droplet, with no outer fluid layer, in an unbounded porous media we have $\mathscr{A}=-1$ , $\unicode[STIX]{x1D702}=-1$ and it follows from (3.75) that

(3.76) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}=-\frac{\text{i}m}{2}\left(1-\frac{1}{We}\frac{m^{2}-1}{2}\right)\left[\frac{1}{2k^{\prime }}+m(m-1)\right]^{-1}, & & \displaystyle\end{eqnarray}$$

where $k^{\prime }=k/r_{0}^{2}$ is the non-dimensional permeability. The differences between (3.76) (that agrees with the result of Alvarez-Lacalle et al. (Reference Alvarez-Lacalle, Ortín and Casademunt2004)) and the result quoted in Schwartz (Reference Schwartz1989) are due to both applying normal stress continuity here, instead of pressure continuity, across the interface and not having a Coriolis term here for consistency with the narrow gap-width and Stokes flow approximations. The dispersion relation given in Schwartz (Reference Schwartz1989) is equal to that of Alvarez-Lacalle et al. (Reference Alvarez-Lacalle, Ortín and Casademunt2004) when $Re\rightarrow 0$ in the Schwartz solution and pressure continuity is enforced in the Alvarez-Lacalle et al. (Reference Alvarez-Lacalle, Ortín and Casademunt2004) solution (see their (13) as opposed to their (18)). The behaviour of the plane-polar two-dimensional system and the Hele-Shaw cell, porous media solutions can be seen to have some qualitatively similar behaviour. Comparing (3.75) with (3.66) in the limit $a\rightarrow \infty$ which gives

(3.77) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}=\frac{\text{i}\mathscr{A}m}{2}\left(1+\frac{1}{We}\frac{m^{2}-1}{2\mathscr{A}}\right)\frac{1}{m^{2}-1}, & & \displaystyle\end{eqnarray}$$

we see that the growth rates differ by a shape factor that depends on the azimuthal wavenumber and the viscosity contrast, but that the dependence on the Atwood number and the Weber number is identical. The key difference between the models is the interaction with the viscous term in the equation of motion and so it is consistent that the behaviour of the models is distinguished via a shape factor that depends on viscosity and on the size of a given perturbation.