Interfacial behaviour in two-fluid Taylor–Couette flow

Interfacial behaviour in two-fluid Taylor–Couette flow Summary The flow of a system of two viscous fluids between two concentric counter-rotating cylinders is discussed. A simple theory is presented that describes the evolution of shape of the interface between the fluids when they have near equal densities and identical viscosities. This suggests that the interface is neutrally stable, but that after sufficient time there are nevertheless points on the profile at which the curvature becomes very large. As a consequence, the interface develops cusp-like portions in its profile. A novel spectral method is developed for this problem in which the interface is represented as a region of finite width and over which the density changes rapidly but smoothly. The results confirm the general predictions of the asymptotic theory for rotation in a horizontal plane but when the rotation occurs vertically additional features develop in the flow. 1. Introduction The Kelvin–Helmholtz instability is a classic example of an interfacial motion that arises when there is a discontinunity in velocity across the interface of two inviscid fluids. Many systems exhibit the phenomenon which often initially appears as waves that gradually wind up to form large spirals. Kelvin–Helmholtz mechanisms can be observed in many practical circumstances including the flow of air over water and are thought to be operative within Jupiter’s red spot. Many textbooks describe the physics behind the instability, including the account given in (1). In a seminal paper Moore (2) proved that when an infinitesimally thin interface between two inviscid fluids begins to deform under the instability, then at certain locations the curvature becomes unbounded within a finite time $$t_c$$. This realisation has serious implications for the accurate tracking of the evolution of an interface as some sort of approximate numerical device is required to follow the development of the flow beyond time $$t_c$$. Most techniques effectively replace the mathematical vortex sheet produced at a thin interface with some form of interfacial zone of finite width. Perhaps the most popular approach is based on the idea of a ‘vortex-blob’, as pioneered by Krasny (3), although a spectral method combined with smoothing has recently been developed by Chen and Forbes (4). Baker et al. (5) have demonstrated that interfacial curvature singularities may arise in the Rayleigh–Taylor instability, in which a heavier fluid overlies a lighter one. In that problem, Forbes (6) has shown numerically that the curvature singularities in inviscid Rayleigh–Taylor flow are replaced by small regions of intense vorticity in the corresponding viscous problem. Curvature singularities along interfaces in viscous fluids are known to occur in a variety of situations and their role in the break-up of jets has been discussed by Eggers and Villermaux (7). Somewhat earlier, Buckmaster (8) had shown that an inviscid bubble within a viscous fluid undergoing straining flow can develop pointed ends which denote locations of infinite curvature. Moreover rising bubbles in non-Newtonian fluids may generate cusp singularities on their trailing side, and this intriguing phenomenon has been studied both theoretically and experimentally by Joseph et al. (9). Recently, Forbes et al. (10) studied the interfacial flow between two viscous fluids that obey the Stokes equations. Numerical schemes based on a new spectral series method detected what seemed to be an effective curvature singularity after a finite time. It is not difficult to appreciate the origin of this phenomenon if we consider the simple model of two viscous fluids moving in the $$x$$–$$y$$ plane within a horizontal channel with walls $$y=\pm 1$$. If the walls move in opposite directions but with the same speed $$c$$, the interface $$y=0$$ is subjected to a viscous shear. For two fluids with near equal densities the situation is that of a simple Couette shear flow for which the horizontal and vertical velocity components are just $$u(x,y)=-cy$$ and $$v(x,y)=0$$. If the interface is perturbed according to a periodic disturbance of amplitude $$A$$ and wavenumber $$k$$, then at time $$t = 0$$ it consists simply of the marker particles $$(x , y) = \left( \xi , A \cos k\xi \right)$$ where $$\xi$$ plays the role of a parameter. After time $$t$$, the initial interface has undergone shear given by the Couette flow so that its shape becomes $$(x , y) = \left( \xi - Act \cos k\xi , A \cos k\xi \right)$$ with an associated curvature   \begin{equation} \kappa = -\frac { A k^2 \cos k\xi} {\bigl[ \bigl(1 + Act k \sin k\xi \bigr)^2 + A^2 k^2\sin^2k\xi \bigr]^{3/2}}. \end{equation} (1.1) Although this curvature does not become precisely singular within finite time, inspection of the denominator shows that it can become very large. Indeed, near crests and troughs of the waveform, $$\sin k\xi\approx 0$$ and if it is also slightly negative then the whole denominator is small at times $$t$$ for which $$1 + Act k \sin k\xi\approx 0$$. This is only possible once $$t$$ is sufficiently large, but assuming this is the case, it is a routine exercise to deduce that when the denominator of (1.1) is least, the corresponding curvature $$\kappa=O(t^3)$$ for $$t\gg 1$$. In practice, accurate numerical methods for the solution of this problem also encounter large spikes in curvature along the interface, and these appear essentially indistinguishable from true singularities, for the reasons just described. Numerical schemes then fail at about the critical time unless particular precautions are taken. This behaviour motivated a spectral solution technique developed for this problem by Forbes et al. (10) who noted that the results bore a superficial similarity to purely inviscid interfaces subject to the Moore curvature singularity (2). The numerical scheme in (10) could only be continued past a critical time by invoking special smoothing techniques at the interface which then overturned and became multi-valued. Previous numerical schemes, such as those described by Pozrikidis (11) and Li and Renardy (12) do not mention this explicitly, although behaviour consistent with this effect is clearly evident in the numerical results of (12) in their Figure 3.3. The objective of this article then is to explore the evolution of interfacial instability in Taylor–Couette flow (that is, flow between two rotating concentric cylinders). First we sketch a simple analysis to demonstrate how, if the fluids are of similar densities, then after some time curvature spikes of a quickly growing size can be anticipated. We then develop a spectral technique which confirms that even though there might be an initially circular interface between the two fluids the large curvature spikes eventually appear. It is then of interest to explore the extent to which the numerical results support the prediction of the elementary asymptotic theory. There is a vast literature discussing Taylor–Couette flow from theoretical, numerical and experimental perspectives. Early experimental work was conducted by Mallock (13), (14) and Couette (15) while pioneering theoretical studies by Andereck et al. (16) uncovered the rich phase-space possessed by the flow. Careful numerical studies were conducted by Marcus (17), (18) and the effect on the stability of two-fluid flow by the inclusion of surfactant has been examined by Peng and Zhu (19). Numerous other papers have explored many stability aspects of Taylor–Couette flow and the reader desiring a more extensive and leisurely discussion are referred to the papers mentioned here and references therein. We structure the remainder of the article by first considering asymptotically the evolution that describes interfacial instability that exists on the boundary between two fluids of near equal densities. This analysis is discussed presently in section 2 and shows that, just as for the plane Couette problem, there are isolated points along the interface at which the curvature can become arbitrarily large. We emphasise that although these values are large, they are never infinite, so are not curvature singularities. Our predictions are confirmed by numerical solutions that are described in section 3. The article closes with a few remarks in section 4. 2. Taylor–Couette flow between rotating cylinders Consider the motion induced in two viscous fluids confined to the region between two concentric cylinders with inner radius $$a$$ and outer radius $$b$$. The smaller cylinder rotates with an angular speed $$\omega$$ while the outer one moves in the opposite direction with steady angular speed $$-\omega$$. Our fluids are supposed to possess identical viscosities and are distinguished by a small difference in density. Furthermore, they are separated by an interface positioned at the neutral radius $$r_n = \sqrt{2 a^2 b^2 / (a^2 + b^2)}$$ at which, although the angular speed is zero, viscous shear is nevertheless present. A periodic perturbation in the azimuthal angle $$\theta$$ is superimposed on the interface, and a sketch of the situation we envisage is shown in Fig. 1. Fig. 1. View largeDownload slide Definition diagram for Taylor–Couette flow between rotating concentric (2.2) for a mode with $$q = 4$$. In this example $$b/a = 2$$, the amplitude is $$\varepsilon = 0.15$$ and the interface is shown at dimensionless time $$\omega t = 1$$ Fig. 1. View largeDownload slide Definition diagram for Taylor–Couette flow between rotating concentric (2.2) for a mode with $$q = 4$$. In this example $$b/a = 2$$, the amplitude is $$\varepsilon = 0.15$$ and the interface is shown at dimensionless time $$\omega t = 1$$ The Navier–Stokes equations expressed in cylindrical polar coordinates admit a simple steady-state solution for the azimuthal component $$v$$ of velocity given by   \begin{equation} v(r , \theta ) = - \frac {\omega(b^2+a^2)}{(b^2 - a^2)} \biggl[ r - \frac{r_n^2}{r} \biggr] , \end{equation} (2.1) where $$r_n$$ is the neutral radius (1). In the special case when both fluids have equal density, so that there is a single effective fluid, the interface is simply a line of marker points in the rotational shear flow (2.1). If a perturbation is made to the shape of the interface, so that initially   \begin{equation} r(\xi ) = r_n + \varepsilon a \cos ( q\xi ),\qquad\quad \theta (\xi ) = \xi , \nonumber \end{equation} where $$q\in\mathbb{N}$$ and the radial amplitude of the disturbance is $$\varepsilon a$$, then after a time $$t$$ the Taylor–Couette flow (2.1) deforms the initial interface to the new shape $$\left( r(\xi ) , \theta (\xi , t) \right)$$ given by   \begin{eqnarray} r ( \xi ) & = & r_n + \varepsilon a \cos ( q\xi ) \nonumber \\ \theta ( \xi ,t ) & = & \xi - \frac {(b^2+a^2)\omega t}{(b^2 - a^2)} \biggl[ 1 - \frac {r_n^2}{r^2 (\xi )} \biggr] . \end{eqnarray} (2.2) The curvature $$\kappa$$ of an interface defined in polar form in terms of the parameter $$\xi$$, becomes   \begin{equation} \kappa = \frac {r^2 (\theta')^3 + 2 (r')^2 \theta' - r \bigl[ r'' \theta' - \theta'' r' \bigr]} {\bigl[ r^2 (\theta')^2 + (r')^2 \bigr]^{3/2}} \,, \end{equation} (2.3) where a dash denotes differentation with respect to $$\xi$$. It is of course technically possible to substitute the expressions (2.2) for the two functions $$r(\xi)$$ and $$\theta (\xi ,t)$$ although this is not necessary for our purposes. Rather, we merely need to note that   \begin{equation} r^2(\theta')^2+(r')^2=\left[ r_n^2+O(\varepsilon)\right]\left[ 1+\frac{2bq\omega(b^2+a^2)}{(b^2-a^2)r_n}\varepsilon t \sin q\xi+O(\varepsilon^2)\right]^2+a^2q^2\varepsilon^2\sin^2(q\xi)\,; \end{equation} (2.4) it is evident that this expression can never exactly vanish but as time evolves so this function can be made very small at locations where   \begin{equation} \varepsilon\sin q\xi\approx -\frac{(b^2-a^2)r_n}{2bq\omega (a^2+b^2)}t^{-1} \nonumber \end{equation} so that the denominator in (2.3) becomes $$O(t^{-3})$$. At such points, it is easily verified that the numerator of (2.3) remains $$O(1)$$ so that overall $$\kappa=O(t^3)$$. Thus, just as for the plane Couette situation, large spikes in the interfacial curvature develop after a certain time elapses, close to crests and troughs of the wave profile. A typical evolution of the curvature distribution is depicted in Fig. 2, which illustrates the curvature of a mode with azimuthal wavenumber $$q = 3$$. The outer cylinder has been assumed to have twice the radius of the inner one, in which case the location of the neutral radius is $$r_n=\sqrt{1.6}a\approx 1.265a$$. The interface is plotted in Fig. 2a and its structure is reminiscent of cylindrical Kelvin–Helmholtz instabilities computed by Forbes and Cosgrove (20). Solutions are shown at the two times $$\omega t = 2$$ and $$\omega t = 4$$ and it is observed that as time progresses so the interface becomes increasingly distorted. Indeed, by the later time $$\omega t=4 $$ the troughs and crests of the wave profiles have developed sharp cusps suggesting the formation of regions of high curvature. This conjecture is confirmed by Fig. 2b, in which the wave profile is shown again, but now with a third axis added to indicate the size of the curvature $$\kappa$$ according to the formula (2.3). On each of the three lobes evident in (a), there are two locations at which the magnitude of the curvature becomes extremely large. Even at a time as modest as $$\omega t=4$$ the curvatures calculated from the formula (2.3) have spikes of significant amplitudes and so the curvature values in Fig. 2b are only shown where $$|\kappa| a < 100$$ to aid visualisation. Fig. 2. View largeDownload slide Results for the wave profiles for Taylor–Couette flow between concentric cylinders with $$b=2a$$. Results are shown for the third mode wave $$q=3$$, amplitude parameter $$\varepsilon = 0.15$$ and the two times $$\omega t = 2$$ (dashed line) and $$\omega t = 4$$ (solid line). The lower panel depicts the interface curvature at the later time Fig. 2. View largeDownload slide Results for the wave profiles for Taylor–Couette flow between concentric cylinders with $$b=2a$$. Results are shown for the third mode wave $$q=3$$, amplitude parameter $$\varepsilon = 0.15$$ and the two times $$\omega t = 2$$ (dashed line) and $$\omega t = 4$$ (solid line). The lower panel depicts the interface curvature at the later time 3. Numerical solution for Taylor–Couette flow To assess the applicability of the asymptotic predictions developed above some independent numerical tests were conducted. Rather than deal with a sharp interface of infinitesimal thickness, a Boussinesq approximation was adopted in which the two-fluid flow is modelled by a single fluid with a density that varies continuously. A feature of this profile is that it incorporates a narrow zone over which the density changes rapidly, but nevertheless smoothly. The use of an abrupt jump in density across the interface in the assumed initial condition has been tried but found to generate small oscillations in the numerical results at later times. Consequently, our infinitesimal interface is replaced with a narrow interfacial zone of finite width, which may be thought of as an effective mixing region for the two fluids. The density $$\rho$$ of the fluid between the rotating cylinders $$a < r < b$$ was taken to be $$\rho = \rho_1 + \bar{\rho}$$ , where $$\rho_1$$ is the reference density in the fluid between the inner cylinder at $$r = a$$ and the interface, and the function $$\bar{\rho} (r, \theta , t)$$ gives the perturbation to that reference density. The continuity equation for the velocity vector q of the fluid is split into an incompressible component   \begin{equation} {\rm div }\,\, {\bf q} = 0 , \end{equation} (3.1) associated with the constant density $$\rho_1$$, and a transport component   \begin{equation} \frac {\partial \bar{\rho}}{\partial t} + {\bf q}\cdot \nabla \bar{\rho} = \sigma \nabla^2 \bar{\rho} \end{equation} (3.2) for the perturbation. In this latter equation, the constant $$\sigma$$ denotes a diffusion coefficient related to a Prandtl number, as discussed by Farrow and Hocking (21). In terms of cylindrical polar coordinates, the velocity vector is $${\bf q} = u\, {\bf e_r} + v\, {\bf e_{\theta}}$$ in which the dimensionless vectors $${\bf e_r}$$ and $${\bf e_{\theta}}$$ are unit vectors in the radial and azimuthal directions respectively. It is possible to satisfy the incompressibility condition (3.1) exactly, by using a streamfunction $$\Psi ( r, \theta ,t)$$ and requiring that the velocity components $$u$$ and $$v$$ obey the conditions   \begin{equation} u = \frac {1}{r} \frac {\partial \Psi}{\partial \theta} \quad {\rm and} \quad v = - \frac {\partial \Psi}{\partial r} . \end{equation} (3.3) It is convenient now to consider the vorticity vector $${\rm curl }\,{\bf q}$$, since for planar flow it has only the single component $$\zeta$$ normal to the plane of the motion, given by   \begin{equation} \zeta = \frac {1}{r} \frac {\partial \bigl( r v \bigr)}{\partial r} - \frac {1}{r} \frac {\partial u}{\partial \theta} \nonumber \end{equation} when written in cylindrical polar coordinates. It is straightforward to show that the vorticity component $$\zeta$$ is related to the streamfunction $$\Psi$$ through a Poisson equation   \begin{equation} \zeta = - \biggl[ \frac {1}{r} \frac {\partial}{\partial r} \biggl( r \frac {\partial \Psi}{\partial r} \biggr) + \frac {1}{r^2} \frac {\partial^2 \Psi}{\partial \theta^2} \biggr] = - \nabla^2 \Psi . \end{equation} (3.4) The conservation of linear momentum for this viscous fluid is expressed by the Navier–Stokes equations. The pressure can be eliminated by taking the vector curl and this leaves the scalar vorticity equation   \begin{equation} \frac {\partial \zeta}{\partial t} + u \frac {\partial \zeta}{\partial r} + \frac {v}{r} \frac {\partial \zeta}{\partial \theta} = \frac {\mu}{\rho_1} \biggl[ \frac {\partial^2 \zeta}{\partial r^2} + \frac {1}{r} \frac {\partial\zeta}{\partial r} + \frac {1}{r^2} \frac {\partial^2 \zeta}{\partial \theta^2} \biggr] . \end{equation} (3.5) This development is classical, and the general form of the vorticity equation is given in the text by Batchelor (22, page 267). It is assumed in (3.5) that any body forces on the fluid, such as gravity, are orthogonal to the plane of the fluid motion and so do not appear; the motion therefore occurs in horizontal planes. These equations are to be solved subject to the boundary conditions   \begin{eqnarray} u = 0 , \quad v = a \omega \quad &{\rm on }\,\quad r = a \nonumber \\ u = 0 , \quad v = - b \omega \quad &{\rm on }\,\quad r = b . \end{eqnarray} (3.6) The numerical solution of these Boussinesq equations is accomplished here using a novel spectral-type method. Initially, a representation is needed for the streamfunction $$\Psi$$ introduced in (3.3) and which satisfies boundary conditions (3.6). A simple function of this type is the asymptotic solution in section 2, for which the streamfunction is simply   \begin{equation} \Psi (r, \theta ,t) = \frac {\omega}{b^2 - a^2} \biggl[ \frac {1}{2} \bigl( a^2 + b^2 \bigr) \bigl( r^2 - a^2 \bigr) - 2a^2 b^2 \ln \biggl( \frac{r}{a} \biggr) \biggr] . \nonumber \end{equation} The vorticity $$\zeta$$ calculated from (3.4) is just a constant for this function, and so represents a trivial solution to (3.5). Moreover, it satisfies the boundary conditions (3.6) identically. A general expression for the streamfunction between the two rotating cylinders is therefore chosen to take the form   \begin{eqnarray} \Psi (r, \theta ,t) & = & \frac {\omega}{b^2 - a^2} \biggl[ \frac {1}{2} \bigl( a^2 + b^2 \bigr) \bigl( r^2 - a^2 \bigr) - 2a^2 b^2 \ln \biggl( \frac{r}{a} \biggr) \biggr] \nonumber \\ & + & \sum_{n=0}^N \sum_{m=1}^{M+2} W_{m,n}(r) \bigl[ P_{mn}(t) \cos (n\theta ) + Q_{mn}(t) \sin (n\theta ) \bigr] , \end{eqnarray} (3.7) in which the functions $$W_{m,n}(r)$$ must satisfy the homogeneous boundary conditions   \begin{equation} W_{m,n}(a) = W_{m,n}(b) = 0 , \end{equation} (3.8) since $$u = 0$$ at the cylinders $$r = a$$, $$b$$ as required by (3.6). It simplifies calculations considerably if the Laplacian $$\nabla^2 \Psi$$ has the same mathematical form as $$\Psi$$ itself and, for this to be the case, it is necessary to choose the functions $$W_{m,n} (r)$$ to be   \begin{equation} W_{mn}(r) = J_n \left( \alpha_{nm} \frac {r}{a} \right) Y_n \left( \alpha_{nm} \right) - J_n \left(\alpha_{nm} \right) Y_n \left( \alpha_{nm} \frac {r}{a} \right), \end{equation} (3.9) in which $$J_n$$ and $$Y_n$$ are respectively the Bessel functions of the first and second kind. Notice that the first condition $$W_{m,n}(a) = 0$$ in (3.8) is satisfied identically by these radial basis functions (3.9) and the second requirement $$W_{m,n}(b) = 0$$ is met if   \begin{equation} J_n \bigl( \alpha_{nm} A \bigr) Y_n \bigl( \alpha_{nm} \bigr) - J_n \bigl( \alpha_{nm} \bigr) Y_n \bigl( \alpha_{nm} A \bigr) = 0 \quad \text{with } \quad A = \frac {b}{a}; \end{equation} (3.10) this serves to define the constants $$\alpha_{nm}$$. The expression on the left of equation (3.10) is often referred to as a Bessel function cross product (see Abramowitz and Stegun (23, page 374)), and it is known that for real $$n$$ there are infinitely many real simple solutions to (3.10). Therefore, the notation $$\alpha_{nm}$$ used in the definition (3.9) denotes the $$m$$-th root of (3.10) with integer order $$n$$ fixed. Abramowitz and Stegun (23, formula 9.5.28) give good estimates for the $$m$$-th root, and these have been used here as a basis from which to calculate the zeros $$\alpha_{nm}$$. It has been found that these roots may be calculated in the correct location and order using the simple bisection algorithm, and accordingly this was the method implemented here. Recently, however, Horsley (24) has developed an efficient algorithm for computing zeros of Bessel function cross products and has proved identities related to their analytical structure. It still remains to satisfy the two conditions $$v (a , \theta , t) = a \omega$$ and $$v (b , \theta , t) = - b \omega$$ in (3.6), and this is achieved using the additional terms at orders $$m = M+1$$ and $$m = M+2$$ in the representation (3.7). The velocity component $$v$$ is calculated from (3.7) using (3.3) and collecting the coefficients of the $$\cos (n\theta )$$ terms gives the two equations   \begin{eqnarray} W_{M+1,n}^{\prime} (a) P_{M+1,n} (t) + W_{M+2,n}^{\prime} (a) P_{M+2,n} (t) & = & - \sum_{m=1}^M W_{m,n}^{\prime} (a) P_{m,n} (t)\,, \nonumber \\ W_{M+1,n}^{\prime} (b) P_{M+1,n} (t) + W_{M+2,n}^{\prime} (b) P_{M+2,n} (t) & = & - \sum_{m=1}^M W_{m,n}^{\prime} (b) P_{m,n} (t)\,, \end{eqnarray} (3.11) for the coefficients $$P_{M+1,n}$$ and $$P_{M+2,n}$$, at each value of $$n = 1 , \dots , N$$. Similar equations for $$Q_{M+1,n}$$ and $$Q_{M+2,n}$$ are obtained from the $$\sin (n\theta )$$ terms. This gives the final representation of the streamfunction to be   \begin{eqnarray} \Psi (r, \theta ,t) & = & \frac {\omega}{b^2 - a^2} \biggl[ \frac {1}{2} \bigl( a^2 + b^2 \bigr) \bigl( r^2 - a^2 \bigr) - 2a^2 b^2 \ln \biggl( \frac{r}{a} \biggr) \biggr] \nonumber \\ & + & \sum_{m=1}^M P_{m0} (t) \bigl[ W_{m,0}(r) - S_{m0} W_{M+1,0}(r) + T_{m0} W_{M+2,0}(r) \bigr] \nonumber \\ & + & \sum_{n=1}^N \sum_{m=1}^M \bigl[ P_{mn}(t) \cos (n\theta ) + Q_{mn}(t) \sin (n\theta ) \bigr] \nonumber \\ & & \quad \quad \times \bigl[ W_{m,n}(r) - S_{mn} W_{M+1,n}(r) + T_{mn} W_{M+2,n}(r) \bigr] . \end{eqnarray} (3.12) The remaining constants in this expression are defined to be   \begin{eqnarray} S_{mn} & = & \frac {\alpha_{n,M+2} \Delta_{M+2,n} - \alpha_{n,m} \Delta_{m,n}} {\alpha_{n,M+2} \Delta_{M+2,n} - \alpha_{n,M+1} \Delta_{M+1,n}}\,, \nonumber \\ T_{mn} & = & \frac {\alpha_{n,M+1} \Delta_{M+1,n} - \alpha_{n,m} \Delta_{m,n}} {\alpha_{n,M+2} \Delta_{M+2,n} - \alpha_{n,M+1} \Delta_{M+1,n}} . \end{eqnarray} (3.13) These results can be deduced using the identity   \begin{equation} W_{j,n}^{\prime} (a) W_{k,n}^{\prime} (b) = - \frac {2 \alpha_{n,k}}{\pi a^2} \Delta_{k,n} \nonumber \end{equation} which may in turn be obtained using standard recurrence relations for the derivatives of Bessel functions (23, page 361) and the definition (3.10) of the zeros $$\alpha_{n,m}$$. The constants $$\Delta_{m,n}$$ are defined to be   \begin{equation} \Delta_{m,n} = J_n \bigl( \alpha_{nm} \bigr) Y_{n+1} \bigl( \alpha_{nm} A \bigr) - J_{n+1} \bigl( \alpha_{nm} A \bigr) Y_n \bigl( \alpha_{nm} \bigr), \end{equation} (3.14) where, again, $$A = b/a$$. With this form (3.12) for the streamfunction and the sets of constants in (3.13) and (3.14), it is now possible to verify directly that the boundary conditions (3.6) for the azimuthal velocity component $$v$$ are satisfied at the two rotating cylinders, making use of the Wronskian for Bessel functions. Next we turn to the density perturbation function $$\bar{\rho}$$ which must vanish on the inner cylinder $$r=a$$ and take the value $$\rho_2-\rho_1$$ on the outer cylinder $$r=b$$. These requirements suggest that an appropriate form for the density perturbation is   \begin{eqnarray} \bar{\rho} (r,\theta ,t) & = & \bigl( \rho_2 - \rho_1 \bigr) \frac {\ln (r/a)}{\ln (b/a)} + \sum_{m=1}^M A_{m0}(t) W_{m0}(r) \nonumber \\ & + & \sum_{n=1}^N \sum_{m=1}^M \bigl[ A_{mn}(t) \cos (n\theta ) + B_{mn}(t) \sin (n\theta ) \bigr] W_{mn}(r) , \end{eqnarray} (3.15) where the functions $$W_{mn}(r)$$ are as defined in (3.9) and vanish on $$r=a$$ and $$r=b$$ by virtue of the conditions (3.10). Since the functions $$W_{mn}(r)$$ are linear combinations of Bessel functions of $$r$$, they must satisfy the Bessel differential equation, and obey the orthogonality condition   \begin{equation} \int_a^b r W_{mn}(r) W_{kn}(r) \, dr = 0 \quad \text{if } m \neq k \end{equation} (3.16) for fixed order $$n$$. Unfortunately, it does not appear possible to evaluate this integral in closed form in the case when $$m = k$$, and so the constants   \begin{equation} {\cal C}_{mn} = \int_a^b r \bigl[ W_{mn}(r) \bigr]^2 \, dr \end{equation} (3.17) were calculated numerically to very high accuracy. We chose to initiate our calculation with the smooth density perturbation profile   \begin{equation} \bar{\rho} (r,\theta ,0 ) = \frac {1}{2} \bigl( \rho_2 - \rho_1 \bigr) \biggl[ 1 + \tanh \bigl( \beta [r - r_n-\varepsilon a\cos(q\theta ) ] \bigr) \biggr] , \nonumber \end{equation} in which the constant $$\beta$$ controls the interface thickness and $$\bar\rho\to 0$$ as $$r\to a$$ and $$\bar\rho\to\rho_2-\rho_1$$ as $$r\to b$$ In practice the parameter $$\beta$$ was chosen to be large so that the main variation in $$\bar\rho$$ occurs over a narrow interval surrounding $$r=r_n$$. It was found that good results were found with $$\beta \approx 20 / a$$. The integer $$q$$ is the mode number of the assumed initial perturbation. A system of differential equations for the unknown coefficient functions $$P_{mn}(t)$$, $$Q_{mn}(t)$$, $$A_{mn}(t)$$ and $$B_{mn}(t)$$ can be derived using Fourier analysis of the vorticity equation (3.5) and the equation (3.2) for the perturbation density. From the vorticity equation it follows that   \begin{eqnarray} \frac {d P_{k \ell}(t)}{d t} & = & - \frac {\mu}{\rho_1 a^2} P_{k \ell}(t) \alpha_{\ell ,k}^2 \nonumber \\ & - & \frac {(2-\delta_{\ell})a^2}{2\pi \alpha_{\ell ,k}^2 {\cal C}_{k \ell}} \int_a^b \int_{-\pi}^{\pi} \biggl( ru \frac {\partial\zeta}{\partial r} + v \frac {\partial\zeta}{\partial\theta} \biggr) \cos (\ell\theta ) W_{k \ell}(r) \, d\theta \, dr \nonumber \\ & & \quad k = 1 , 2 , \dots , M , \quad \ell =0, 1 , 2 , \dots , N , \end{eqnarray} (3.18) where $$\delta_\ell=1$$ if $$\ell=0$$ and zero otherwise. A similar set of differential equations is obtained for the other coefficients $$Q_{k \ell}(t)$$ in the representation (3.12), except that the functions $$\cos (\ell\theta )$$ in the integrand are replaced by $$\sin (\ell\theta )$$. The transport equation (3.2) for the density leads to the system   \begin{eqnarray} \frac {d A_{k \ell}(t)}{d t} & = & - \frac {\sigma}{a^2} A_{k \ell}(t) \alpha_{\ell ,k}^2 \nonumber \\ & - & \frac {2-\delta_\ell}{2\pi {\cal C}_{k \ell}} \int_a^b \int_{-\pi}^{\pi} \biggl( ru \frac {\partial\bar{\rho}}{\partial r} + v \frac {\partial\bar{\rho}}{\partial\theta} \biggr) \cos (\ell\theta ) W_{k \ell}(r) \, d\theta \, dr \nonumber \\ & & \quad k = 1 , 2 , \dots , M , \quad \ell = 0, 1 , 2 , \dots , N . \end{eqnarray} (3.19) Once more, the corresponding differential equation for the coefficients $$B_{k \ell}(t)$$ is obtained by replacing $$\cos (\ell\theta )$$ with $$\sin (\ell\theta )$$ in the integrand in (3.19). The system (3.18)–(3.19) was integrated forward in time with the adaptive Runge–Kutta–Fehlberg package ode45 supplied as part of the MATLAB suite of routines proving well suited for this purpose. The fluid was taken to be initially at rest which is ensured by taking all the $$P_{mn}(0)=Q_{mn}(0)=0$$. Following (2.2), an interface perturbed at the $$q$$-th wavenumber was defined to be $${\cal R} ( \theta ) = r_n + \varepsilon a \cos ( q\theta ).$$ Fourier analysis of the initial profile for $$\bar{\rho}$$ based on the representation (3.15) at time $$t = 0$$ then shows that   \begin{equation} A_{m0}(0) = \frac {1}{2 \pi {\cal C}_{m0}} \int_a^b \int_{-\pi}^{\pi} \biggl[ \bar{\rho} (r, \theta ,0) - \bigl( \rho_2 - \rho_1 \bigr) \frac {\ln (r/a)}{\ln (b/a)} \biggr] r W_{m 0}(r) \, d\theta \, dr \nonumber \end{equation} and   \begin{equation} A_{mn}(0) = \frac {1}{\pi {\cal C}_{mn}} \int_a^b \int_{-\pi}^{\pi} \bar{\rho} (r, \theta ,0) \cos (n\theta ) r W_{m n}(r) \, d\theta \, dr . \nonumber \end{equation} A similar formula for $$B_{mn}(0)$$ can be deduced and has the same form, except that the term $$\cos (n\theta )$$ is replaced with $$\sin (n\theta )$$ in the integrand. Some sample results expressed in dimensionless units are illustrated in Fig. 3. Here all lengths are written in terms of the inner cylinder radius $$a$$ while the reference for time is the inverse rotational frequency $$1 / \omega$$ and the density scale is $$\rho_1$$, the density of the inner fluid. This re-scaling throws up three further dimensionless parameters; a density ratio $$D = {\rho_2}/{\rho_1}$$, a parameter $$S=a^2\omega/\sigma$$ related to a Prandtl number and a Reynolds number $$R_e = {\rho_1 a^2 \omega} /\mu$$ based on the viscosity $$\mu$$ of the fluids. For definiteness we set the values $$S=R_e=10^4$$ although their precise values are not crucial. For completeness, a Froude number $$F$$ was defined based on the component of the gravitational acceleration vector $${\bf g}$$ acting across the plane of rotation; if the unit vector $${\bf n}$$ denotes the normal to the rotation plane, then the Froude number $$F$$ is obtained from the relation Fig. 3. View largeDownload slide Interfacial profiles at two different dimensionless times (a) $$\omega t = 2$$ and (b) $$\omega t = 4$$, for rotational Taylor–Couette flow between concentric cylinders. Results are shown in dimensionless form, for a third-mode wave with $$q = 3$$, an amplitude parameter $$\varepsilon = 0.15$$ and relative cylinder radii $$b/a = 2$$. The dashed line in each diagram is the interface shape predicted by the asymptotic solution (2.2) Fig. 3. View largeDownload slide Interfacial profiles at two different dimensionless times (a) $$\omega t = 2$$ and (b) $$\omega t = 4$$, for rotational Taylor–Couette flow between concentric cylinders. Results are shown in dimensionless form, for a third-mode wave with $$q = 3$$, an amplitude parameter $$\varepsilon = 0.15$$ and relative cylinder radii $$b/a = 2$$. The dashed line in each diagram is the interface shape predicted by the asymptotic solution (2.2)   \begin{equation} \frac {1}{F^2} = \frac {\| {\bf g} \times {\bf n} \|}{a \omega^2} . \end{equation} (3.20) Clearly this parameter plays no role in the results of Fig. 3, since rotation has been assumed to occur in a horizontal plane, but this condition will be relaxed later. Contours of the dimensionless density perturbation $$\bar{\rho} / \rho_1$$ are depicted in Fig. 3, and the two regions of different density are evident in both diagrams. The interfacial zone is quite narrow at the earlier time $$\omega t = 2$$ but has widened slightly by the time $$\omega t = 4$$ shown in Fig. 3b. As time progresses, the sinusoidal fingers around the neutral radius become increasingly distorted by the counter-rotation of the inner and outer cylinders to the extent that eventually overhanging fingers of each fluid appear in the other. These fingers continue to elongate and thin, as is evident in Fig. 3b. The numerical algorithm does not continue easily for times much larger than the value $$\omega t = 4$$ shown in this diagram, since slender fingers of fluid apparently break up rapidly into a complicated mixed region that is increasingly difficult to resolve numerically. The results shown in Fig. 3 have been obtained using $$M = 31$$ and $$N = 41$$ Fourier coefficients, with $$121$$ numerical points in the radial coordinate and $$161$$ points azimuthally. The results are well converged, having been compared carefully against predictions generated with fewer Fourier coefficients. The dashed lines superposed on Fig. 3 denote the interfacial shapes predicted by the asymptotic solution (2.2) at these two different times. The agreement with the numerical solution is clearly excellent. The asymptotic interface lies precisely along the boundary of the two different fluids in Fig. 3a when $$\omega t = 2$$, and this property is maintained to the later time $$\omega t = 4$$ in Fig. 3b, despite the interfacial zone having widened a little in the interim. This gives strong confidence in the usefulness of the asymptotic development in section 2. To conclude this presentation of the numerical results, a situation is now considered in which the rotation takes place in a vertical plane, so that the acceleration of gravity is now directed across the rotating fluids. In this case, the Froude number $$F$$ defined in (3.20) now has a finite value. While this situation might perhaps be slightly contrived from the point of view of the physics, it is nevertheless interesting because, in addition to the Kelvin–Helmholtz instability caused by shear at the interface, there will be the possibility that Rayleigh–Taylor type instabilities come into play. This is likely since heavier fluid will overlie lighter fluid, and the two will attempt to exchange positions. Thus, as time increases, it is perhaps to be expected that the less dense fluid will accumulate at the top and the heavier fluid will gravitate towards the bottom, destroying the three-fold rotational symmetry in Fig. 3. The inclusion of gravitational effects requires that the additional term   \begin{equation} - \frac {g}{\rho_1} \biggl[ \cos\theta \frac {\partial \bar{\rho}}{\partial r} - \frac {\sin\theta}{r} \frac {\partial \bar{\rho}}{\partial \theta} \biggr] \nonumber \end{equation} be added to the right-hand side of the vorticity equation (3.5) that governs the behaviour of the vorticity component $$\zeta$$. This is a term that accounts for buoyancy between the different density fluid layers. It must also be included in (3.18) for the coefficients $$P_{mn}(t)$$ to give relations of the form   \begin{eqnarray} \frac {d P_{k \ell}(t)}{d t} & = & - \frac {\mu}{\rho_1 a^2} P_{k \ell}(t) \alpha_{\ell ,k}^2 \nonumber \\ & - & \frac {a^2}{\pi \alpha_{\ell ,k}^2 {\cal C}_{k \ell}} \int_a^b \int_{-\pi}^{\pi} r {\cal H} (r, \theta , t) \cos (\ell\theta ) W_{k \ell}(r) \, d\theta \, dr \nonumber \\ & & \quad k = 1 , 2 , \dots , M , \quad \ell = 1 , 2 , \dots , N , \nonumber \end{eqnarray} in which   \begin{equation} {\cal H} (r, \theta , t) = \biggl( u \frac {\partial\zeta}{\partial r} + \frac {v}{r} \frac {\partial\zeta}{\partial\theta} \biggr) + \frac {g}{\rho_1} \biggl( \cos\theta \frac {\partial \bar{\rho}}{\partial r} - \frac {\sin\theta}{r} \frac {\partial \bar{\rho}}{\partial \theta} \biggr) . \nonumber \end{equation} A similar set of modified differential equations applies for the coefficients $$Q_{mn}(t)$$. Figure 4 shows some sample results for $$F = 1/2$$, in which $$F = \sqrt{a \omega^2 /g}$$ is the Froude number defined in (3.20). Once again, the density ratio is $$D = 1.2$$ and the ratio of the cylinder radii is $$b/a = 2$$. The first diagram in part (a) presents density contours at the dimensionless time $$\omega t = 0.75$$ for a solution with $$q = 3$$ and amplitude $$\varepsilon = 0.15$$. At this relatively early time, the three waves in the solution are clearly visible and have only just started to overturn. Some smaller wavelength ripples along the interface may be visible in some regions, and for comparison, the asymptotic interface (2.2), valid for rotation in a horizontal plane, has been overlaid on the diagram. The agreement with the numerical solution at time $$\omega t = 0.75$$ is still reasonable, although the three-fold rotational symmetry of the initial disturbance has begun to disappear. Figure 4b shows the density profiles at the later time $$\omega t = 1.25$$. Again the asymptotic result (2.2) is shown for comparison, and demonstrates the extent to which the rotational symmetry has been destroyed by the additional gravitational effect. In particular, jets of the lighter inner fluid are now beginning to penetrate further into the region at the top of the diagram, which is as expected. There are small-amplitude shorter waves riding on the back of the jet on the right-hand side of the figure, and these are a combination of Kelvin–Helmholtz and Rayleigh–Taylor instabilities that ultimately contribute to mixing of the two fluids on short length scales, making the numerical solution beyond this time increasingly difficult. Fig. 4. View largeDownload slide Interfacial profiles at two different dimensionless times (a) $$\omega t = 0.75$$ and (b) $$\omega t = 1.25$$, for rotational Taylor–Couette flow between concentric cylinders. Results are shown in dimensionless form, for a third-mode wave with $$q = 3$$ and amplitude parameter $$\varepsilon = 0.15$$. The density ratio is $$D = 1.2$$, the cylinder radii are $$b/a = 2$$, and the Froude number is $$F = 1/2$$. The dashed line in each diagram is the interface shape predicted by the asymptotic solution (2.2) for flow in a horizontal plane Fig. 4. View largeDownload slide Interfacial profiles at two different dimensionless times (a) $$\omega t = 0.75$$ and (b) $$\omega t = 1.25$$, for rotational Taylor–Couette flow between concentric cylinders. Results are shown in dimensionless form, for a third-mode wave with $$q = 3$$ and amplitude parameter $$\varepsilon = 0.15$$. The density ratio is $$D = 1.2$$, the cylinder radii are $$b/a = 2$$, and the Froude number is $$F = 1/2$$. The dashed line in each diagram is the interface shape predicted by the asymptotic solution (2.2) for flow in a horizontal plane Equation (3.12) was evaluated to give the streamfunction $$\Psi$$ at the same values of the parameters as used for Fig. 4, and when $$\omega t = 1.25$$. Some representative contours are plotted in Fig. 5. For steady flow, the velocity vector $${\bf q}$$ is parallel to these level curves of $$\Psi$$, and so the contours in Fig. 5 give at least an approximate description of the locations of the streamlines. As expected, many of these are circles, reflecting the obvious fact that the dominant flow pattern is simply a result of the rotation of the two cylinders. Nevertheless, there is also a small recirculating zone on the right side of the picture, showing that the jet at about that location in Fig. 4b is associated with the development of significant vorticity there. This is confirmed by an examination of the vorticity function $$\zeta$$, which is computed from the spectral representation (3.12) combined with the knowledge that $$\zeta=-\nabla^2\Psi$$. Fig. 5. View largeDownload slide Streamlines for the case shown in Fig. 4, at the later time $$\omega t = 1.25$$ Fig. 5. View largeDownload slide Streamlines for the case shown in Fig. 4, at the later time $$\omega t = 1.25$$ Some comments on the accuracy of the numerical solutions to this Boussinesq model are appropriate here. A demonstration of the accuracy of these results has already been given in Fig. 3, where the close agreement between the asymptotic theory and the Boussinesq viscous results with $$M = 31$$, $$N = 41$$ Fourier coefficients has been noted. As a further check on the convergence of the numerical results with increasing numbers of Fourier coefficients, Fig. 6 shows a comparison between results obtained with $$M = 21$$, $$N = 31$$ coefficients and the solution for the same case generated with $$M = 31$$, $$N = 41$$ coefficients. These solutions are for the same parameter set as used in Figs. 4 and 5, and at the dimensionless time $$\omega t = 1$$. There is very close agreement between these two results, as is evident from the figures, although the more accurate results in Fig. 6b are capable of resolving some features at a finer length scale, as is to be expected. Thus, although the broad shapes of the three interfacial fingers are the same, it is possible to see small-scale Kelvin–Helmholtz waves on the two fingers at the bottom and the right of the diagram in Fig. 6b, that are not quite so precisely resolved in Fig. 6a. Fig. 6. View largeDownload slide A comparison of results for the density using the same parameters as in Fig. 4, at dimensionless time $$\omega t = 1$$, with (a) $$M = 21$$, $$N = 31$$ Fourier coefficients and (b) $$M = 31$$, $$N = 41$$ coefficients Fig. 6. View largeDownload slide A comparison of results for the density using the same parameters as in Fig. 4, at dimensionless time $$\omega t = 1$$, with (a) $$M = 21$$, $$N = 31$$ Fourier coefficients and (b) $$M = 31$$, $$N = 41$$ coefficients A further example of the effects of rotation under gravity is presented in Fig. 7. The parameters are again the same as in Fig. 4, except that the Froude number has been increased to the value $$F = 1$$; this may be thought of either as increasing the rotation rate of the cylinders or else mitigating the effect of gravity. A solution is illustrated at the dimensionless time $$\omega t = 3$$. This represents about the latest time, for this Froude number, at which reliable numerical results could be achieved. The reason for this is apparent from Fig. 7 and is associated with the very fine length scales that develop as the two fluids mix across their interface. Again, the solution was started with an interfacial disturbance with $$q = 3$$, but by time $$\omega t = 3$$ in Fig. 7, the bottom-most finger has detached from the main body of the inner fluid, and has almost completely disappeared through mixing. The remaining two fingers are towards the top of the diagram, as expected, and they, too, have developed into complicated structures, apparently with detached portions, and noticeable small-amplitude instability waves ride along the overall elongated finger structures of the jets. Fig. 7. View largeDownload slide Interfacial profile at dimensionless time $$\omega t = 3$$, for rotational Taylor–Couette flow between concentric cylinders. Results are shown in dimensionless form, for a third-mode wave with $$q = 3$$ and amplitude parameter $$\varepsilon = 0.15$$ The density ratio is $$D = 1.2$$, the radii are $$b/a = 2$$ and the Froude number is $$F = 1$$ Fig. 7. View largeDownload slide Interfacial profile at dimensionless time $$\omega t = 3$$, for rotational Taylor–Couette flow between concentric cylinders. Results are shown in dimensionless form, for a third-mode wave with $$q = 3$$ and amplitude parameter $$\varepsilon = 0.15$$ The density ratio is $$D = 1.2$$, the radii are $$b/a = 2$$ and the Froude number is $$F = 1$$ 4. Conclusions When two flowing inviscid fluids are separated by an interface of infinitesimal thickness, the interface is unstable. According to linearised theory, any small-amplitude disturbance to the interface will grow exponentially with time if the two fluids move with different speeds. Forbes and Cosgrove (20) recently showed that a similar situation exists in cylindrical geometry, where two adjacent fluids move under the effects of a line vortex. An interface disturbance that grows exponentially, however, eventually violates the linearisation assumption that disturbances are small, so that non-linear effects ultimately become important. Moore (2) has shown how non-linearity results in the formation of a curvature singularity at the interface, for inviscid fluids, and within finite time. The Moore curvature singularity is a subtle and complicated effect, and Moore’s asymptotic analysis (2) is a difficult calculation. Intuitively, it might be anticipated that the addition of fluid viscosity would so damp the interface that Moore’s curvature singularity now could no longer form. This is indeed the case; however, pathological behaviour at a narrow interface between viscous fluids can nevertheless still occur, in spite of the fact that there can be no velocity jump at the interface, and even when the interfacial disturbance starts from a smooth differentiable shape at $$t = 0$$. This is arguably a somewhat surprising fact and appears not to have been commented upon explicitly in the literature. Furthermore, unlike the complicated asymptotic analysis leading to Moore’s (2) result for the curvature singularity at an interface between inviscid fluids, in the viscous case a very simple analysis reveals this behaviour. For viscous flow of two fluids in a parallel-walled channel, Forbes et al. (10) demonstrated that a careful numerical solution revealed somewhat strange behaviour in the curvature of a narrow interface, and suggested that curvature singularity might be the cause. A significant aim of the present article has therefore been to explore whether similar behaviour at the interface might again be possible, and we have presented a simple analytical demonstration that, only after a certain time has elapsed, curvature spikes of rapidly-growing magnitude are indeed possible. In that sense, they appear numerically to be very similar to the genuine singularity of Moore (2) for inviscid fluids, even although the curvature magnitude in this viscous case does not become infinite within finite time. Perhaps the main aim of this present article has been to develop the somewhat novel spectral technique for solving the non-linear problem described in Section 3, and to assess the extent to which the simple asymptotic analysis of section 2 can represent the solution behaviour. There is an initially circular interface between the two fluids, and the asymptotic analysis shows that large curvatures can develop, even when the initial disturbance to the interface is quite smooth. The straightforward analysis of section 2 is sufficient to reveal this effect, and the accuracy of this asymptotic formula (2.2) has indeed been confirmed by the numerical solution of section 3. Nevertheless, there is an important difference between the spectral numerical solutions of section 3 and the analysis of section 2. The former is based on Boussinesq theory, in which the exact mathematical interface is replaced by a narrow zone of finite width, over which density changes rapidly, but smoothly. Thus, although there is close agreement with the asymptotic solutions in Fig. 3, the numerical results can nevertheless not produce the same near-singular behaviour of the asymptotic solutions. The important conclusion, then, is that the Moore singularity for inviscid fluids and the curvature spikes that develop with viscous fluids are both consequences of the assumption of an infinitessimally thin interface; the success of ‘vortex blob’ methods (3), (5) for numerical solution of the inviscid problem is not that they mimic viscosity, but rather that they create a diffuse interface of finite width, which prevents the Moore singularity from developing. Furthermore, as the curvature spikes discussed in this article are able to be predicted using quite simple arguments, it seems likely that they could occur quite generally at narrow interfaces between flowing viscous fluids. The phenomenon therefore appears worthy of further study. In addition, it may be possible to use simple results of the type presented in Section 2 as the basis of more complete asymptotic analyses. Acknowledgements We are grateful to the referees for numerous helpful suggestions that helped us clarify many aspects of the article. Discussion with Rhys A. Paul and Stephen J. Walters is gratefully acknowledged. This work is associated with Australian Research Council Discovery Grant DP140100094. References 1. Drazin P. G. and Reid W. H. Hydrodynamic Stability , 2nd edn. ( Cambridge University Press, Cambridge 2004). Google Scholar CrossRef Search ADS   2. Moore D. W. The spontaneous appearance of a singularity in the shape of an evolving vortex sheet, Proc. Roy. Soc. London A  365 ( 1979) 105– 119. Google Scholar CrossRef Search ADS   3. Krasny R. Desingularization of periodic vortex sheet roll-up, J. Comput. Phys.  65 ( 1986) 4. Chen M. J. and Forbes L. K. Accurate methods for computing inviscid and viscous Kelvin–Helmholtz instability, J. Comput. Phys.  230 ( 2011) 1499– 1515. Google Scholar CrossRef Search ADS   5. Baker G. Caflisch R. E. and Siegel M. Singularity formation during Rayleigh–Taylor instability, J. Fluid Mech.  252 ( 1993) 51– 78. Google Scholar CrossRef Search ADS   6. Forbes L. K., The Rayleigh-Taylor instability for inviscid and viscous fluids, J. Engin. Math.  65 ( 2009) 273– 290. Google Scholar CrossRef Search ADS   7. Eggers J. and Villermaux E. Physics of liquid jets, Rep. Prog. Phys.  71 ( 2008) 036601. Google Scholar CrossRef Search ADS   8. Buckmaster. J. D. Pointed bubbles in slow viscous flow, J. Fluid Mech.  55 ( 1972) 385– 400. Google Scholar CrossRef Search ADS   9. Joseph D. D. Nelson J. Renardy M. and Renardy Y. Two-dimensional cusped interfaces, J. Fluid Mech.  223 ( 1991) 383– 409. Google Scholar CrossRef Search ADS   10. Forbes L. K. Paul R. A. Chen M. J. and Horsley D. E. Kelvin–Helmholtz creeping flow at the interface between two viscous fluids, ANZIAM J.  56 ( 2015) 317– 358. Google Scholar CrossRef Search ADS   11. Pozrikidis C. Instability of two-layer creeping flow in a channel with parallel-sided walls, J. Fluid Mech.  351 ( 1997) 139– 165. Google Scholar CrossRef Search ADS   12. Li J. and Renardy Y. Numerical study of flows of two immiscible liquids at low Reynolds number, SIAM Rev.  42 ( 2000) 417– 439. Google Scholar CrossRef Search ADS   13. Mallock A. Determination of the viscosity of water, Proc. R. Soc. Lond. A  45 ( 1888) 126– 132. Google Scholar CrossRef Search ADS   14. Mallock A. Experiments on fluid viscosity, Phil. Trans. R. Soc. Lond. A  187 ( 1896) 41– 56. Google Scholar CrossRef Search ADS   15. Couette M. Études sur le frottement des liquides, Ann. Chim. Phys.  6 ( 1890) 433– 510. 16. Andereck C. D. Liu S. S. and Swinney H. L. Flow regimes in a circular Couette system with independently rotating cylinders, J. Fluid Mech.  164 ( 1986) 155– 183. Google Scholar CrossRef Search ADS   17. Marcus P. S. Simulation of Taylor–Couette flow. Part 1. Numerical methods and comparison with experiment, J. Fluid Mech.  146 ( 1984) 45– 64. Google Scholar CrossRef Search ADS   18. Marcus P. S. Simulation of Taylor–Couette flow. Part 2. Numerical results for wavy-vortex flow with one travelling wave, J. Fluid Mech.  146 ( 1984) 65– 113. Google Scholar CrossRef Search ADS   19. Peng J. and Zhu K.-Q. Linear instability of two-fluid Taylor–Couette flow in the presence of surfactant, J. Fluid Mech.  651 ( 2010) 357– 385. Google Scholar CrossRef Search ADS   20. Forbes L. K. and Cosgrove J. M., A line vortex in a two-fluid system, J. Engin. Math.  84 ( 2014), 181– 199. Google Scholar CrossRef Search ADS   21. Farrow D. E. and Hocking G. C. A numerical model for withdrawal from a two-layer fluid, J. Fluid Mech.  549 ( 2006) 141– 157. Google Scholar CrossRef Search ADS   22. Batchelor G. K. An Introduction to Fluid Mechanics. ( Cambridge University Press 1967). Google Scholar CrossRef Search ADS   23. Abramowitz M. and Stegun I. A. (eds), Handbook of Mathematical Functions . ( Dover, New York 1972). 24. Horsley D. E. Bessel phase functions: calculation and application, Numerische Math.  136 ( 2017) 679– 702. Google Scholar CrossRef Search ADS   © The Author, 2017. Published by Oxford University Press; all rights reserved. For Permissions, please email: journals.permissions@oup.com http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png The Quarterly Journal of Mechanics and Applied Mathematics Oxford University Press

Interfacial behaviour in two-fluid Taylor–Couette flow

Loading next page...
 
/lp/ou_press/interfacial-behaviour-in-two-fluid-taylor-couette-flow-i2EQ5Kmox8
Publisher
Oxford University Press
Copyright
© The Author, 2017. Published by Oxford University Press; all rights reserved. For Permissions, please email: journals.permissions@oup.com
ISSN
0033-5614
eISSN
1464-3855
D.O.I.
10.1093/qjmam/hbx025
Publisher site
See Article on Publisher Site

Abstract

Summary The flow of a system of two viscous fluids between two concentric counter-rotating cylinders is discussed. A simple theory is presented that describes the evolution of shape of the interface between the fluids when they have near equal densities and identical viscosities. This suggests that the interface is neutrally stable, but that after sufficient time there are nevertheless points on the profile at which the curvature becomes very large. As a consequence, the interface develops cusp-like portions in its profile. A novel spectral method is developed for this problem in which the interface is represented as a region of finite width and over which the density changes rapidly but smoothly. The results confirm the general predictions of the asymptotic theory for rotation in a horizontal plane but when the rotation occurs vertically additional features develop in the flow. 1. Introduction The Kelvin–Helmholtz instability is a classic example of an interfacial motion that arises when there is a discontinunity in velocity across the interface of two inviscid fluids. Many systems exhibit the phenomenon which often initially appears as waves that gradually wind up to form large spirals. Kelvin–Helmholtz mechanisms can be observed in many practical circumstances including the flow of air over water and are thought to be operative within Jupiter’s red spot. Many textbooks describe the physics behind the instability, including the account given in (1). In a seminal paper Moore (2) proved that when an infinitesimally thin interface between two inviscid fluids begins to deform under the instability, then at certain locations the curvature becomes unbounded within a finite time $$t_c$$. This realisation has serious implications for the accurate tracking of the evolution of an interface as some sort of approximate numerical device is required to follow the development of the flow beyond time $$t_c$$. Most techniques effectively replace the mathematical vortex sheet produced at a thin interface with some form of interfacial zone of finite width. Perhaps the most popular approach is based on the idea of a ‘vortex-blob’, as pioneered by Krasny (3), although a spectral method combined with smoothing has recently been developed by Chen and Forbes (4). Baker et al. (5) have demonstrated that interfacial curvature singularities may arise in the Rayleigh–Taylor instability, in which a heavier fluid overlies a lighter one. In that problem, Forbes (6) has shown numerically that the curvature singularities in inviscid Rayleigh–Taylor flow are replaced by small regions of intense vorticity in the corresponding viscous problem. Curvature singularities along interfaces in viscous fluids are known to occur in a variety of situations and their role in the break-up of jets has been discussed by Eggers and Villermaux (7). Somewhat earlier, Buckmaster (8) had shown that an inviscid bubble within a viscous fluid undergoing straining flow can develop pointed ends which denote locations of infinite curvature. Moreover rising bubbles in non-Newtonian fluids may generate cusp singularities on their trailing side, and this intriguing phenomenon has been studied both theoretically and experimentally by Joseph et al. (9). Recently, Forbes et al. (10) studied the interfacial flow between two viscous fluids that obey the Stokes equations. Numerical schemes based on a new spectral series method detected what seemed to be an effective curvature singularity after a finite time. It is not difficult to appreciate the origin of this phenomenon if we consider the simple model of two viscous fluids moving in the $$x$$–$$y$$ plane within a horizontal channel with walls $$y=\pm 1$$. If the walls move in opposite directions but with the same speed $$c$$, the interface $$y=0$$ is subjected to a viscous shear. For two fluids with near equal densities the situation is that of a simple Couette shear flow for which the horizontal and vertical velocity components are just $$u(x,y)=-cy$$ and $$v(x,y)=0$$. If the interface is perturbed according to a periodic disturbance of amplitude $$A$$ and wavenumber $$k$$, then at time $$t = 0$$ it consists simply of the marker particles $$(x , y) = \left( \xi , A \cos k\xi \right)$$ where $$\xi$$ plays the role of a parameter. After time $$t$$, the initial interface has undergone shear given by the Couette flow so that its shape becomes $$(x , y) = \left( \xi - Act \cos k\xi , A \cos k\xi \right)$$ with an associated curvature   \begin{equation} \kappa = -\frac { A k^2 \cos k\xi} {\bigl[ \bigl(1 + Act k \sin k\xi \bigr)^2 + A^2 k^2\sin^2k\xi \bigr]^{3/2}}. \end{equation} (1.1) Although this curvature does not become precisely singular within finite time, inspection of the denominator shows that it can become very large. Indeed, near crests and troughs of the waveform, $$\sin k\xi\approx 0$$ and if it is also slightly negative then the whole denominator is small at times $$t$$ for which $$1 + Act k \sin k\xi\approx 0$$. This is only possible once $$t$$ is sufficiently large, but assuming this is the case, it is a routine exercise to deduce that when the denominator of (1.1) is least, the corresponding curvature $$\kappa=O(t^3)$$ for $$t\gg 1$$. In practice, accurate numerical methods for the solution of this problem also encounter large spikes in curvature along the interface, and these appear essentially indistinguishable from true singularities, for the reasons just described. Numerical schemes then fail at about the critical time unless particular precautions are taken. This behaviour motivated a spectral solution technique developed for this problem by Forbes et al. (10) who noted that the results bore a superficial similarity to purely inviscid interfaces subject to the Moore curvature singularity (2). The numerical scheme in (10) could only be continued past a critical time by invoking special smoothing techniques at the interface which then overturned and became multi-valued. Previous numerical schemes, such as those described by Pozrikidis (11) and Li and Renardy (12) do not mention this explicitly, although behaviour consistent with this effect is clearly evident in the numerical results of (12) in their Figure 3.3. The objective of this article then is to explore the evolution of interfacial instability in Taylor–Couette flow (that is, flow between two rotating concentric cylinders). First we sketch a simple analysis to demonstrate how, if the fluids are of similar densities, then after some time curvature spikes of a quickly growing size can be anticipated. We then develop a spectral technique which confirms that even though there might be an initially circular interface between the two fluids the large curvature spikes eventually appear. It is then of interest to explore the extent to which the numerical results support the prediction of the elementary asymptotic theory. There is a vast literature discussing Taylor–Couette flow from theoretical, numerical and experimental perspectives. Early experimental work was conducted by Mallock (13), (14) and Couette (15) while pioneering theoretical studies by Andereck et al. (16) uncovered the rich phase-space possessed by the flow. Careful numerical studies were conducted by Marcus (17), (18) and the effect on the stability of two-fluid flow by the inclusion of surfactant has been examined by Peng and Zhu (19). Numerous other papers have explored many stability aspects of Taylor–Couette flow and the reader desiring a more extensive and leisurely discussion are referred to the papers mentioned here and references therein. We structure the remainder of the article by first considering asymptotically the evolution that describes interfacial instability that exists on the boundary between two fluids of near equal densities. This analysis is discussed presently in section 2 and shows that, just as for the plane Couette problem, there are isolated points along the interface at which the curvature can become arbitrarily large. We emphasise that although these values are large, they are never infinite, so are not curvature singularities. Our predictions are confirmed by numerical solutions that are described in section 3. The article closes with a few remarks in section 4. 2. Taylor–Couette flow between rotating cylinders Consider the motion induced in two viscous fluids confined to the region between two concentric cylinders with inner radius $$a$$ and outer radius $$b$$. The smaller cylinder rotates with an angular speed $$\omega$$ while the outer one moves in the opposite direction with steady angular speed $$-\omega$$. Our fluids are supposed to possess identical viscosities and are distinguished by a small difference in density. Furthermore, they are separated by an interface positioned at the neutral radius $$r_n = \sqrt{2 a^2 b^2 / (a^2 + b^2)}$$ at which, although the angular speed is zero, viscous shear is nevertheless present. A periodic perturbation in the azimuthal angle $$\theta$$ is superimposed on the interface, and a sketch of the situation we envisage is shown in Fig. 1. Fig. 1. View largeDownload slide Definition diagram for Taylor–Couette flow between rotating concentric (2.2) for a mode with $$q = 4$$. In this example $$b/a = 2$$, the amplitude is $$\varepsilon = 0.15$$ and the interface is shown at dimensionless time $$\omega t = 1$$ Fig. 1. View largeDownload slide Definition diagram for Taylor–Couette flow between rotating concentric (2.2) for a mode with $$q = 4$$. In this example $$b/a = 2$$, the amplitude is $$\varepsilon = 0.15$$ and the interface is shown at dimensionless time $$\omega t = 1$$ The Navier–Stokes equations expressed in cylindrical polar coordinates admit a simple steady-state solution for the azimuthal component $$v$$ of velocity given by   \begin{equation} v(r , \theta ) = - \frac {\omega(b^2+a^2)}{(b^2 - a^2)} \biggl[ r - \frac{r_n^2}{r} \biggr] , \end{equation} (2.1) where $$r_n$$ is the neutral radius (1). In the special case when both fluids have equal density, so that there is a single effective fluid, the interface is simply a line of marker points in the rotational shear flow (2.1). If a perturbation is made to the shape of the interface, so that initially   \begin{equation} r(\xi ) = r_n + \varepsilon a \cos ( q\xi ),\qquad\quad \theta (\xi ) = \xi , \nonumber \end{equation} where $$q\in\mathbb{N}$$ and the radial amplitude of the disturbance is $$\varepsilon a$$, then after a time $$t$$ the Taylor–Couette flow (2.1) deforms the initial interface to the new shape $$\left( r(\xi ) , \theta (\xi , t) \right)$$ given by   \begin{eqnarray} r ( \xi ) & = & r_n + \varepsilon a \cos ( q\xi ) \nonumber \\ \theta ( \xi ,t ) & = & \xi - \frac {(b^2+a^2)\omega t}{(b^2 - a^2)} \biggl[ 1 - \frac {r_n^2}{r^2 (\xi )} \biggr] . \end{eqnarray} (2.2) The curvature $$\kappa$$ of an interface defined in polar form in terms of the parameter $$\xi$$, becomes   \begin{equation} \kappa = \frac {r^2 (\theta')^3 + 2 (r')^2 \theta' - r \bigl[ r'' \theta' - \theta'' r' \bigr]} {\bigl[ r^2 (\theta')^2 + (r')^2 \bigr]^{3/2}} \,, \end{equation} (2.3) where a dash denotes differentation with respect to $$\xi$$. It is of course technically possible to substitute the expressions (2.2) for the two functions $$r(\xi)$$ and $$\theta (\xi ,t)$$ although this is not necessary for our purposes. Rather, we merely need to note that   \begin{equation} r^2(\theta')^2+(r')^2=\left[ r_n^2+O(\varepsilon)\right]\left[ 1+\frac{2bq\omega(b^2+a^2)}{(b^2-a^2)r_n}\varepsilon t \sin q\xi+O(\varepsilon^2)\right]^2+a^2q^2\varepsilon^2\sin^2(q\xi)\,; \end{equation} (2.4) it is evident that this expression can never exactly vanish but as time evolves so this function can be made very small at locations where   \begin{equation} \varepsilon\sin q\xi\approx -\frac{(b^2-a^2)r_n}{2bq\omega (a^2+b^2)}t^{-1} \nonumber \end{equation} so that the denominator in (2.3) becomes $$O(t^{-3})$$. At such points, it is easily verified that the numerator of (2.3) remains $$O(1)$$ so that overall $$\kappa=O(t^3)$$. Thus, just as for the plane Couette situation, large spikes in the interfacial curvature develop after a certain time elapses, close to crests and troughs of the wave profile. A typical evolution of the curvature distribution is depicted in Fig. 2, which illustrates the curvature of a mode with azimuthal wavenumber $$q = 3$$. The outer cylinder has been assumed to have twice the radius of the inner one, in which case the location of the neutral radius is $$r_n=\sqrt{1.6}a\approx 1.265a$$. The interface is plotted in Fig. 2a and its structure is reminiscent of cylindrical Kelvin–Helmholtz instabilities computed by Forbes and Cosgrove (20). Solutions are shown at the two times $$\omega t = 2$$ and $$\omega t = 4$$ and it is observed that as time progresses so the interface becomes increasingly distorted. Indeed, by the later time $$\omega t=4 $$ the troughs and crests of the wave profiles have developed sharp cusps suggesting the formation of regions of high curvature. This conjecture is confirmed by Fig. 2b, in which the wave profile is shown again, but now with a third axis added to indicate the size of the curvature $$\kappa$$ according to the formula (2.3). On each of the three lobes evident in (a), there are two locations at which the magnitude of the curvature becomes extremely large. Even at a time as modest as $$\omega t=4$$ the curvatures calculated from the formula (2.3) have spikes of significant amplitudes and so the curvature values in Fig. 2b are only shown where $$|\kappa| a < 100$$ to aid visualisation. Fig. 2. View largeDownload slide Results for the wave profiles for Taylor–Couette flow between concentric cylinders with $$b=2a$$. Results are shown for the third mode wave $$q=3$$, amplitude parameter $$\varepsilon = 0.15$$ and the two times $$\omega t = 2$$ (dashed line) and $$\omega t = 4$$ (solid line). The lower panel depicts the interface curvature at the later time Fig. 2. View largeDownload slide Results for the wave profiles for Taylor–Couette flow between concentric cylinders with $$b=2a$$. Results are shown for the third mode wave $$q=3$$, amplitude parameter $$\varepsilon = 0.15$$ and the two times $$\omega t = 2$$ (dashed line) and $$\omega t = 4$$ (solid line). The lower panel depicts the interface curvature at the later time 3. Numerical solution for Taylor–Couette flow To assess the applicability of the asymptotic predictions developed above some independent numerical tests were conducted. Rather than deal with a sharp interface of infinitesimal thickness, a Boussinesq approximation was adopted in which the two-fluid flow is modelled by a single fluid with a density that varies continuously. A feature of this profile is that it incorporates a narrow zone over which the density changes rapidly, but nevertheless smoothly. The use of an abrupt jump in density across the interface in the assumed initial condition has been tried but found to generate small oscillations in the numerical results at later times. Consequently, our infinitesimal interface is replaced with a narrow interfacial zone of finite width, which may be thought of as an effective mixing region for the two fluids. The density $$\rho$$ of the fluid between the rotating cylinders $$a < r < b$$ was taken to be $$\rho = \rho_1 + \bar{\rho}$$ , where $$\rho_1$$ is the reference density in the fluid between the inner cylinder at $$r = a$$ and the interface, and the function $$\bar{\rho} (r, \theta , t)$$ gives the perturbation to that reference density. The continuity equation for the velocity vector q of the fluid is split into an incompressible component   \begin{equation} {\rm div }\,\, {\bf q} = 0 , \end{equation} (3.1) associated with the constant density $$\rho_1$$, and a transport component   \begin{equation} \frac {\partial \bar{\rho}}{\partial t} + {\bf q}\cdot \nabla \bar{\rho} = \sigma \nabla^2 \bar{\rho} \end{equation} (3.2) for the perturbation. In this latter equation, the constant $$\sigma$$ denotes a diffusion coefficient related to a Prandtl number, as discussed by Farrow and Hocking (21). In terms of cylindrical polar coordinates, the velocity vector is $${\bf q} = u\, {\bf e_r} + v\, {\bf e_{\theta}}$$ in which the dimensionless vectors $${\bf e_r}$$ and $${\bf e_{\theta}}$$ are unit vectors in the radial and azimuthal directions respectively. It is possible to satisfy the incompressibility condition (3.1) exactly, by using a streamfunction $$\Psi ( r, \theta ,t)$$ and requiring that the velocity components $$u$$ and $$v$$ obey the conditions   \begin{equation} u = \frac {1}{r} \frac {\partial \Psi}{\partial \theta} \quad {\rm and} \quad v = - \frac {\partial \Psi}{\partial r} . \end{equation} (3.3) It is convenient now to consider the vorticity vector $${\rm curl }\,{\bf q}$$, since for planar flow it has only the single component $$\zeta$$ normal to the plane of the motion, given by   \begin{equation} \zeta = \frac {1}{r} \frac {\partial \bigl( r v \bigr)}{\partial r} - \frac {1}{r} \frac {\partial u}{\partial \theta} \nonumber \end{equation} when written in cylindrical polar coordinates. It is straightforward to show that the vorticity component $$\zeta$$ is related to the streamfunction $$\Psi$$ through a Poisson equation   \begin{equation} \zeta = - \biggl[ \frac {1}{r} \frac {\partial}{\partial r} \biggl( r \frac {\partial \Psi}{\partial r} \biggr) + \frac {1}{r^2} \frac {\partial^2 \Psi}{\partial \theta^2} \biggr] = - \nabla^2 \Psi . \end{equation} (3.4) The conservation of linear momentum for this viscous fluid is expressed by the Navier–Stokes equations. The pressure can be eliminated by taking the vector curl and this leaves the scalar vorticity equation   \begin{equation} \frac {\partial \zeta}{\partial t} + u \frac {\partial \zeta}{\partial r} + \frac {v}{r} \frac {\partial \zeta}{\partial \theta} = \frac {\mu}{\rho_1} \biggl[ \frac {\partial^2 \zeta}{\partial r^2} + \frac {1}{r} \frac {\partial\zeta}{\partial r} + \frac {1}{r^2} \frac {\partial^2 \zeta}{\partial \theta^2} \biggr] . \end{equation} (3.5) This development is classical, and the general form of the vorticity equation is given in the text by Batchelor (22, page 267). It is assumed in (3.5) that any body forces on the fluid, such as gravity, are orthogonal to the plane of the fluid motion and so do not appear; the motion therefore occurs in horizontal planes. These equations are to be solved subject to the boundary conditions   \begin{eqnarray} u = 0 , \quad v = a \omega \quad &{\rm on }\,\quad r = a \nonumber \\ u = 0 , \quad v = - b \omega \quad &{\rm on }\,\quad r = b . \end{eqnarray} (3.6) The numerical solution of these Boussinesq equations is accomplished here using a novel spectral-type method. Initially, a representation is needed for the streamfunction $$\Psi$$ introduced in (3.3) and which satisfies boundary conditions (3.6). A simple function of this type is the asymptotic solution in section 2, for which the streamfunction is simply   \begin{equation} \Psi (r, \theta ,t) = \frac {\omega}{b^2 - a^2} \biggl[ \frac {1}{2} \bigl( a^2 + b^2 \bigr) \bigl( r^2 - a^2 \bigr) - 2a^2 b^2 \ln \biggl( \frac{r}{a} \biggr) \biggr] . \nonumber \end{equation} The vorticity $$\zeta$$ calculated from (3.4) is just a constant for this function, and so represents a trivial solution to (3.5). Moreover, it satisfies the boundary conditions (3.6) identically. A general expression for the streamfunction between the two rotating cylinders is therefore chosen to take the form   \begin{eqnarray} \Psi (r, \theta ,t) & = & \frac {\omega}{b^2 - a^2} \biggl[ \frac {1}{2} \bigl( a^2 + b^2 \bigr) \bigl( r^2 - a^2 \bigr) - 2a^2 b^2 \ln \biggl( \frac{r}{a} \biggr) \biggr] \nonumber \\ & + & \sum_{n=0}^N \sum_{m=1}^{M+2} W_{m,n}(r) \bigl[ P_{mn}(t) \cos (n\theta ) + Q_{mn}(t) \sin (n\theta ) \bigr] , \end{eqnarray} (3.7) in which the functions $$W_{m,n}(r)$$ must satisfy the homogeneous boundary conditions   \begin{equation} W_{m,n}(a) = W_{m,n}(b) = 0 , \end{equation} (3.8) since $$u = 0$$ at the cylinders $$r = a$$, $$b$$ as required by (3.6). It simplifies calculations considerably if the Laplacian $$\nabla^2 \Psi$$ has the same mathematical form as $$\Psi$$ itself and, for this to be the case, it is necessary to choose the functions $$W_{m,n} (r)$$ to be   \begin{equation} W_{mn}(r) = J_n \left( \alpha_{nm} \frac {r}{a} \right) Y_n \left( \alpha_{nm} \right) - J_n \left(\alpha_{nm} \right) Y_n \left( \alpha_{nm} \frac {r}{a} \right), \end{equation} (3.9) in which $$J_n$$ and $$Y_n$$ are respectively the Bessel functions of the first and second kind. Notice that the first condition $$W_{m,n}(a) = 0$$ in (3.8) is satisfied identically by these radial basis functions (3.9) and the second requirement $$W_{m,n}(b) = 0$$ is met if   \begin{equation} J_n \bigl( \alpha_{nm} A \bigr) Y_n \bigl( \alpha_{nm} \bigr) - J_n \bigl( \alpha_{nm} \bigr) Y_n \bigl( \alpha_{nm} A \bigr) = 0 \quad \text{with } \quad A = \frac {b}{a}; \end{equation} (3.10) this serves to define the constants $$\alpha_{nm}$$. The expression on the left of equation (3.10) is often referred to as a Bessel function cross product (see Abramowitz and Stegun (23, page 374)), and it is known that for real $$n$$ there are infinitely many real simple solutions to (3.10). Therefore, the notation $$\alpha_{nm}$$ used in the definition (3.9) denotes the $$m$$-th root of (3.10) with integer order $$n$$ fixed. Abramowitz and Stegun (23, formula 9.5.28) give good estimates for the $$m$$-th root, and these have been used here as a basis from which to calculate the zeros $$\alpha_{nm}$$. It has been found that these roots may be calculated in the correct location and order using the simple bisection algorithm, and accordingly this was the method implemented here. Recently, however, Horsley (24) has developed an efficient algorithm for computing zeros of Bessel function cross products and has proved identities related to their analytical structure. It still remains to satisfy the two conditions $$v (a , \theta , t) = a \omega$$ and $$v (b , \theta , t) = - b \omega$$ in (3.6), and this is achieved using the additional terms at orders $$m = M+1$$ and $$m = M+2$$ in the representation (3.7). The velocity component $$v$$ is calculated from (3.7) using (3.3) and collecting the coefficients of the $$\cos (n\theta )$$ terms gives the two equations   \begin{eqnarray} W_{M+1,n}^{\prime} (a) P_{M+1,n} (t) + W_{M+2,n}^{\prime} (a) P_{M+2,n} (t) & = & - \sum_{m=1}^M W_{m,n}^{\prime} (a) P_{m,n} (t)\,, \nonumber \\ W_{M+1,n}^{\prime} (b) P_{M+1,n} (t) + W_{M+2,n}^{\prime} (b) P_{M+2,n} (t) & = & - \sum_{m=1}^M W_{m,n}^{\prime} (b) P_{m,n} (t)\,, \end{eqnarray} (3.11) for the coefficients $$P_{M+1,n}$$ and $$P_{M+2,n}$$, at each value of $$n = 1 , \dots , N$$. Similar equations for $$Q_{M+1,n}$$ and $$Q_{M+2,n}$$ are obtained from the $$\sin (n\theta )$$ terms. This gives the final representation of the streamfunction to be   \begin{eqnarray} \Psi (r, \theta ,t) & = & \frac {\omega}{b^2 - a^2} \biggl[ \frac {1}{2} \bigl( a^2 + b^2 \bigr) \bigl( r^2 - a^2 \bigr) - 2a^2 b^2 \ln \biggl( \frac{r}{a} \biggr) \biggr] \nonumber \\ & + & \sum_{m=1}^M P_{m0} (t) \bigl[ W_{m,0}(r) - S_{m0} W_{M+1,0}(r) + T_{m0} W_{M+2,0}(r) \bigr] \nonumber \\ & + & \sum_{n=1}^N \sum_{m=1}^M \bigl[ P_{mn}(t) \cos (n\theta ) + Q_{mn}(t) \sin (n\theta ) \bigr] \nonumber \\ & & \quad \quad \times \bigl[ W_{m,n}(r) - S_{mn} W_{M+1,n}(r) + T_{mn} W_{M+2,n}(r) \bigr] . \end{eqnarray} (3.12) The remaining constants in this expression are defined to be   \begin{eqnarray} S_{mn} & = & \frac {\alpha_{n,M+2} \Delta_{M+2,n} - \alpha_{n,m} \Delta_{m,n}} {\alpha_{n,M+2} \Delta_{M+2,n} - \alpha_{n,M+1} \Delta_{M+1,n}}\,, \nonumber \\ T_{mn} & = & \frac {\alpha_{n,M+1} \Delta_{M+1,n} - \alpha_{n,m} \Delta_{m,n}} {\alpha_{n,M+2} \Delta_{M+2,n} - \alpha_{n,M+1} \Delta_{M+1,n}} . \end{eqnarray} (3.13) These results can be deduced using the identity   \begin{equation} W_{j,n}^{\prime} (a) W_{k,n}^{\prime} (b) = - \frac {2 \alpha_{n,k}}{\pi a^2} \Delta_{k,n} \nonumber \end{equation} which may in turn be obtained using standard recurrence relations for the derivatives of Bessel functions (23, page 361) and the definition (3.10) of the zeros $$\alpha_{n,m}$$. The constants $$\Delta_{m,n}$$ are defined to be   \begin{equation} \Delta_{m,n} = J_n \bigl( \alpha_{nm} \bigr) Y_{n+1} \bigl( \alpha_{nm} A \bigr) - J_{n+1} \bigl( \alpha_{nm} A \bigr) Y_n \bigl( \alpha_{nm} \bigr), \end{equation} (3.14) where, again, $$A = b/a$$. With this form (3.12) for the streamfunction and the sets of constants in (3.13) and (3.14), it is now possible to verify directly that the boundary conditions (3.6) for the azimuthal velocity component $$v$$ are satisfied at the two rotating cylinders, making use of the Wronskian for Bessel functions. Next we turn to the density perturbation function $$\bar{\rho}$$ which must vanish on the inner cylinder $$r=a$$ and take the value $$\rho_2-\rho_1$$ on the outer cylinder $$r=b$$. These requirements suggest that an appropriate form for the density perturbation is   \begin{eqnarray} \bar{\rho} (r,\theta ,t) & = & \bigl( \rho_2 - \rho_1 \bigr) \frac {\ln (r/a)}{\ln (b/a)} + \sum_{m=1}^M A_{m0}(t) W_{m0}(r) \nonumber \\ & + & \sum_{n=1}^N \sum_{m=1}^M \bigl[ A_{mn}(t) \cos (n\theta ) + B_{mn}(t) \sin (n\theta ) \bigr] W_{mn}(r) , \end{eqnarray} (3.15) where the functions $$W_{mn}(r)$$ are as defined in (3.9) and vanish on $$r=a$$ and $$r=b$$ by virtue of the conditions (3.10). Since the functions $$W_{mn}(r)$$ are linear combinations of Bessel functions of $$r$$, they must satisfy the Bessel differential equation, and obey the orthogonality condition   \begin{equation} \int_a^b r W_{mn}(r) W_{kn}(r) \, dr = 0 \quad \text{if } m \neq k \end{equation} (3.16) for fixed order $$n$$. Unfortunately, it does not appear possible to evaluate this integral in closed form in the case when $$m = k$$, and so the constants   \begin{equation} {\cal C}_{mn} = \int_a^b r \bigl[ W_{mn}(r) \bigr]^2 \, dr \end{equation} (3.17) were calculated numerically to very high accuracy. We chose to initiate our calculation with the smooth density perturbation profile   \begin{equation} \bar{\rho} (r,\theta ,0 ) = \frac {1}{2} \bigl( \rho_2 - \rho_1 \bigr) \biggl[ 1 + \tanh \bigl( \beta [r - r_n-\varepsilon a\cos(q\theta ) ] \bigr) \biggr] , \nonumber \end{equation} in which the constant $$\beta$$ controls the interface thickness and $$\bar\rho\to 0$$ as $$r\to a$$ and $$\bar\rho\to\rho_2-\rho_1$$ as $$r\to b$$ In practice the parameter $$\beta$$ was chosen to be large so that the main variation in $$\bar\rho$$ occurs over a narrow interval surrounding $$r=r_n$$. It was found that good results were found with $$\beta \approx 20 / a$$. The integer $$q$$ is the mode number of the assumed initial perturbation. A system of differential equations for the unknown coefficient functions $$P_{mn}(t)$$, $$Q_{mn}(t)$$, $$A_{mn}(t)$$ and $$B_{mn}(t)$$ can be derived using Fourier analysis of the vorticity equation (3.5) and the equation (3.2) for the perturbation density. From the vorticity equation it follows that   \begin{eqnarray} \frac {d P_{k \ell}(t)}{d t} & = & - \frac {\mu}{\rho_1 a^2} P_{k \ell}(t) \alpha_{\ell ,k}^2 \nonumber \\ & - & \frac {(2-\delta_{\ell})a^2}{2\pi \alpha_{\ell ,k}^2 {\cal C}_{k \ell}} \int_a^b \int_{-\pi}^{\pi} \biggl( ru \frac {\partial\zeta}{\partial r} + v \frac {\partial\zeta}{\partial\theta} \biggr) \cos (\ell\theta ) W_{k \ell}(r) \, d\theta \, dr \nonumber \\ & & \quad k = 1 , 2 , \dots , M , \quad \ell =0, 1 , 2 , \dots , N , \end{eqnarray} (3.18) where $$\delta_\ell=1$$ if $$\ell=0$$ and zero otherwise. A similar set of differential equations is obtained for the other coefficients $$Q_{k \ell}(t)$$ in the representation (3.12), except that the functions $$\cos (\ell\theta )$$ in the integrand are replaced by $$\sin (\ell\theta )$$. The transport equation (3.2) for the density leads to the system   \begin{eqnarray} \frac {d A_{k \ell}(t)}{d t} & = & - \frac {\sigma}{a^2} A_{k \ell}(t) \alpha_{\ell ,k}^2 \nonumber \\ & - & \frac {2-\delta_\ell}{2\pi {\cal C}_{k \ell}} \int_a^b \int_{-\pi}^{\pi} \biggl( ru \frac {\partial\bar{\rho}}{\partial r} + v \frac {\partial\bar{\rho}}{\partial\theta} \biggr) \cos (\ell\theta ) W_{k \ell}(r) \, d\theta \, dr \nonumber \\ & & \quad k = 1 , 2 , \dots , M , \quad \ell = 0, 1 , 2 , \dots , N . \end{eqnarray} (3.19) Once more, the corresponding differential equation for the coefficients $$B_{k \ell}(t)$$ is obtained by replacing $$\cos (\ell\theta )$$ with $$\sin (\ell\theta )$$ in the integrand in (3.19). The system (3.18)–(3.19) was integrated forward in time with the adaptive Runge–Kutta–Fehlberg package ode45 supplied as part of the MATLAB suite of routines proving well suited for this purpose. The fluid was taken to be initially at rest which is ensured by taking all the $$P_{mn}(0)=Q_{mn}(0)=0$$. Following (2.2), an interface perturbed at the $$q$$-th wavenumber was defined to be $${\cal R} ( \theta ) = r_n + \varepsilon a \cos ( q\theta ).$$ Fourier analysis of the initial profile for $$\bar{\rho}$$ based on the representation (3.15) at time $$t = 0$$ then shows that   \begin{equation} A_{m0}(0) = \frac {1}{2 \pi {\cal C}_{m0}} \int_a^b \int_{-\pi}^{\pi} \biggl[ \bar{\rho} (r, \theta ,0) - \bigl( \rho_2 - \rho_1 \bigr) \frac {\ln (r/a)}{\ln (b/a)} \biggr] r W_{m 0}(r) \, d\theta \, dr \nonumber \end{equation} and   \begin{equation} A_{mn}(0) = \frac {1}{\pi {\cal C}_{mn}} \int_a^b \int_{-\pi}^{\pi} \bar{\rho} (r, \theta ,0) \cos (n\theta ) r W_{m n}(r) \, d\theta \, dr . \nonumber \end{equation} A similar formula for $$B_{mn}(0)$$ can be deduced and has the same form, except that the term $$\cos (n\theta )$$ is replaced with $$\sin (n\theta )$$ in the integrand. Some sample results expressed in dimensionless units are illustrated in Fig. 3. Here all lengths are written in terms of the inner cylinder radius $$a$$ while the reference for time is the inverse rotational frequency $$1 / \omega$$ and the density scale is $$\rho_1$$, the density of the inner fluid. This re-scaling throws up three further dimensionless parameters; a density ratio $$D = {\rho_2}/{\rho_1}$$, a parameter $$S=a^2\omega/\sigma$$ related to a Prandtl number and a Reynolds number $$R_e = {\rho_1 a^2 \omega} /\mu$$ based on the viscosity $$\mu$$ of the fluids. For definiteness we set the values $$S=R_e=10^4$$ although their precise values are not crucial. For completeness, a Froude number $$F$$ was defined based on the component of the gravitational acceleration vector $${\bf g}$$ acting across the plane of rotation; if the unit vector $${\bf n}$$ denotes the normal to the rotation plane, then the Froude number $$F$$ is obtained from the relation Fig. 3. View largeDownload slide Interfacial profiles at two different dimensionless times (a) $$\omega t = 2$$ and (b) $$\omega t = 4$$, for rotational Taylor–Couette flow between concentric cylinders. Results are shown in dimensionless form, for a third-mode wave with $$q = 3$$, an amplitude parameter $$\varepsilon = 0.15$$ and relative cylinder radii $$b/a = 2$$. The dashed line in each diagram is the interface shape predicted by the asymptotic solution (2.2) Fig. 3. View largeDownload slide Interfacial profiles at two different dimensionless times (a) $$\omega t = 2$$ and (b) $$\omega t = 4$$, for rotational Taylor–Couette flow between concentric cylinders. Results are shown in dimensionless form, for a third-mode wave with $$q = 3$$, an amplitude parameter $$\varepsilon = 0.15$$ and relative cylinder radii $$b/a = 2$$. The dashed line in each diagram is the interface shape predicted by the asymptotic solution (2.2)   \begin{equation} \frac {1}{F^2} = \frac {\| {\bf g} \times {\bf n} \|}{a \omega^2} . \end{equation} (3.20) Clearly this parameter plays no role in the results of Fig. 3, since rotation has been assumed to occur in a horizontal plane, but this condition will be relaxed later. Contours of the dimensionless density perturbation $$\bar{\rho} / \rho_1$$ are depicted in Fig. 3, and the two regions of different density are evident in both diagrams. The interfacial zone is quite narrow at the earlier time $$\omega t = 2$$ but has widened slightly by the time $$\omega t = 4$$ shown in Fig. 3b. As time progresses, the sinusoidal fingers around the neutral radius become increasingly distorted by the counter-rotation of the inner and outer cylinders to the extent that eventually overhanging fingers of each fluid appear in the other. These fingers continue to elongate and thin, as is evident in Fig. 3b. The numerical algorithm does not continue easily for times much larger than the value $$\omega t = 4$$ shown in this diagram, since slender fingers of fluid apparently break up rapidly into a complicated mixed region that is increasingly difficult to resolve numerically. The results shown in Fig. 3 have been obtained using $$M = 31$$ and $$N = 41$$ Fourier coefficients, with $$121$$ numerical points in the radial coordinate and $$161$$ points azimuthally. The results are well converged, having been compared carefully against predictions generated with fewer Fourier coefficients. The dashed lines superposed on Fig. 3 denote the interfacial shapes predicted by the asymptotic solution (2.2) at these two different times. The agreement with the numerical solution is clearly excellent. The asymptotic interface lies precisely along the boundary of the two different fluids in Fig. 3a when $$\omega t = 2$$, and this property is maintained to the later time $$\omega t = 4$$ in Fig. 3b, despite the interfacial zone having widened a little in the interim. This gives strong confidence in the usefulness of the asymptotic development in section 2. To conclude this presentation of the numerical results, a situation is now considered in which the rotation takes place in a vertical plane, so that the acceleration of gravity is now directed across the rotating fluids. In this case, the Froude number $$F$$ defined in (3.20) now has a finite value. While this situation might perhaps be slightly contrived from the point of view of the physics, it is nevertheless interesting because, in addition to the Kelvin–Helmholtz instability caused by shear at the interface, there will be the possibility that Rayleigh–Taylor type instabilities come into play. This is likely since heavier fluid will overlie lighter fluid, and the two will attempt to exchange positions. Thus, as time increases, it is perhaps to be expected that the less dense fluid will accumulate at the top and the heavier fluid will gravitate towards the bottom, destroying the three-fold rotational symmetry in Fig. 3. The inclusion of gravitational effects requires that the additional term   \begin{equation} - \frac {g}{\rho_1} \biggl[ \cos\theta \frac {\partial \bar{\rho}}{\partial r} - \frac {\sin\theta}{r} \frac {\partial \bar{\rho}}{\partial \theta} \biggr] \nonumber \end{equation} be added to the right-hand side of the vorticity equation (3.5) that governs the behaviour of the vorticity component $$\zeta$$. This is a term that accounts for buoyancy between the different density fluid layers. It must also be included in (3.18) for the coefficients $$P_{mn}(t)$$ to give relations of the form   \begin{eqnarray} \frac {d P_{k \ell}(t)}{d t} & = & - \frac {\mu}{\rho_1 a^2} P_{k \ell}(t) \alpha_{\ell ,k}^2 \nonumber \\ & - & \frac {a^2}{\pi \alpha_{\ell ,k}^2 {\cal C}_{k \ell}} \int_a^b \int_{-\pi}^{\pi} r {\cal H} (r, \theta , t) \cos (\ell\theta ) W_{k \ell}(r) \, d\theta \, dr \nonumber \\ & & \quad k = 1 , 2 , \dots , M , \quad \ell = 1 , 2 , \dots , N , \nonumber \end{eqnarray} in which   \begin{equation} {\cal H} (r, \theta , t) = \biggl( u \frac {\partial\zeta}{\partial r} + \frac {v}{r} \frac {\partial\zeta}{\partial\theta} \biggr) + \frac {g}{\rho_1} \biggl( \cos\theta \frac {\partial \bar{\rho}}{\partial r} - \frac {\sin\theta}{r} \frac {\partial \bar{\rho}}{\partial \theta} \biggr) . \nonumber \end{equation} A similar set of modified differential equations applies for the coefficients $$Q_{mn}(t)$$. Figure 4 shows some sample results for $$F = 1/2$$, in which $$F = \sqrt{a \omega^2 /g}$$ is the Froude number defined in (3.20). Once again, the density ratio is $$D = 1.2$$ and the ratio of the cylinder radii is $$b/a = 2$$. The first diagram in part (a) presents density contours at the dimensionless time $$\omega t = 0.75$$ for a solution with $$q = 3$$ and amplitude $$\varepsilon = 0.15$$. At this relatively early time, the three waves in the solution are clearly visible and have only just started to overturn. Some smaller wavelength ripples along the interface may be visible in some regions, and for comparison, the asymptotic interface (2.2), valid for rotation in a horizontal plane, has been overlaid on the diagram. The agreement with the numerical solution at time $$\omega t = 0.75$$ is still reasonable, although the three-fold rotational symmetry of the initial disturbance has begun to disappear. Figure 4b shows the density profiles at the later time $$\omega t = 1.25$$. Again the asymptotic result (2.2) is shown for comparison, and demonstrates the extent to which the rotational symmetry has been destroyed by the additional gravitational effect. In particular, jets of the lighter inner fluid are now beginning to penetrate further into the region at the top of the diagram, which is as expected. There are small-amplitude shorter waves riding on the back of the jet on the right-hand side of the figure, and these are a combination of Kelvin–Helmholtz and Rayleigh–Taylor instabilities that ultimately contribute to mixing of the two fluids on short length scales, making the numerical solution beyond this time increasingly difficult. Fig. 4. View largeDownload slide Interfacial profiles at two different dimensionless times (a) $$\omega t = 0.75$$ and (b) $$\omega t = 1.25$$, for rotational Taylor–Couette flow between concentric cylinders. Results are shown in dimensionless form, for a third-mode wave with $$q = 3$$ and amplitude parameter $$\varepsilon = 0.15$$. The density ratio is $$D = 1.2$$, the cylinder radii are $$b/a = 2$$, and the Froude number is $$F = 1/2$$. The dashed line in each diagram is the interface shape predicted by the asymptotic solution (2.2) for flow in a horizontal plane Fig. 4. View largeDownload slide Interfacial profiles at two different dimensionless times (a) $$\omega t = 0.75$$ and (b) $$\omega t = 1.25$$, for rotational Taylor–Couette flow between concentric cylinders. Results are shown in dimensionless form, for a third-mode wave with $$q = 3$$ and amplitude parameter $$\varepsilon = 0.15$$. The density ratio is $$D = 1.2$$, the cylinder radii are $$b/a = 2$$, and the Froude number is $$F = 1/2$$. The dashed line in each diagram is the interface shape predicted by the asymptotic solution (2.2) for flow in a horizontal plane Equation (3.12) was evaluated to give the streamfunction $$\Psi$$ at the same values of the parameters as used for Fig. 4, and when $$\omega t = 1.25$$. Some representative contours are plotted in Fig. 5. For steady flow, the velocity vector $${\bf q}$$ is parallel to these level curves of $$\Psi$$, and so the contours in Fig. 5 give at least an approximate description of the locations of the streamlines. As expected, many of these are circles, reflecting the obvious fact that the dominant flow pattern is simply a result of the rotation of the two cylinders. Nevertheless, there is also a small recirculating zone on the right side of the picture, showing that the jet at about that location in Fig. 4b is associated with the development of significant vorticity there. This is confirmed by an examination of the vorticity function $$\zeta$$, which is computed from the spectral representation (3.12) combined with the knowledge that $$\zeta=-\nabla^2\Psi$$. Fig. 5. View largeDownload slide Streamlines for the case shown in Fig. 4, at the later time $$\omega t = 1.25$$ Fig. 5. View largeDownload slide Streamlines for the case shown in Fig. 4, at the later time $$\omega t = 1.25$$ Some comments on the accuracy of the numerical solutions to this Boussinesq model are appropriate here. A demonstration of the accuracy of these results has already been given in Fig. 3, where the close agreement between the asymptotic theory and the Boussinesq viscous results with $$M = 31$$, $$N = 41$$ Fourier coefficients has been noted. As a further check on the convergence of the numerical results with increasing numbers of Fourier coefficients, Fig. 6 shows a comparison between results obtained with $$M = 21$$, $$N = 31$$ coefficients and the solution for the same case generated with $$M = 31$$, $$N = 41$$ coefficients. These solutions are for the same parameter set as used in Figs. 4 and 5, and at the dimensionless time $$\omega t = 1$$. There is very close agreement between these two results, as is evident from the figures, although the more accurate results in Fig. 6b are capable of resolving some features at a finer length scale, as is to be expected. Thus, although the broad shapes of the three interfacial fingers are the same, it is possible to see small-scale Kelvin–Helmholtz waves on the two fingers at the bottom and the right of the diagram in Fig. 6b, that are not quite so precisely resolved in Fig. 6a. Fig. 6. View largeDownload slide A comparison of results for the density using the same parameters as in Fig. 4, at dimensionless time $$\omega t = 1$$, with (a) $$M = 21$$, $$N = 31$$ Fourier coefficients and (b) $$M = 31$$, $$N = 41$$ coefficients Fig. 6. View largeDownload slide A comparison of results for the density using the same parameters as in Fig. 4, at dimensionless time $$\omega t = 1$$, with (a) $$M = 21$$, $$N = 31$$ Fourier coefficients and (b) $$M = 31$$, $$N = 41$$ coefficients A further example of the effects of rotation under gravity is presented in Fig. 7. The parameters are again the same as in Fig. 4, except that the Froude number has been increased to the value $$F = 1$$; this may be thought of either as increasing the rotation rate of the cylinders or else mitigating the effect of gravity. A solution is illustrated at the dimensionless time $$\omega t = 3$$. This represents about the latest time, for this Froude number, at which reliable numerical results could be achieved. The reason for this is apparent from Fig. 7 and is associated with the very fine length scales that develop as the two fluids mix across their interface. Again, the solution was started with an interfacial disturbance with $$q = 3$$, but by time $$\omega t = 3$$ in Fig. 7, the bottom-most finger has detached from the main body of the inner fluid, and has almost completely disappeared through mixing. The remaining two fingers are towards the top of the diagram, as expected, and they, too, have developed into complicated structures, apparently with detached portions, and noticeable small-amplitude instability waves ride along the overall elongated finger structures of the jets. Fig. 7. View largeDownload slide Interfacial profile at dimensionless time $$\omega t = 3$$, for rotational Taylor–Couette flow between concentric cylinders. Results are shown in dimensionless form, for a third-mode wave with $$q = 3$$ and amplitude parameter $$\varepsilon = 0.15$$ The density ratio is $$D = 1.2$$, the radii are $$b/a = 2$$ and the Froude number is $$F = 1$$ Fig. 7. View largeDownload slide Interfacial profile at dimensionless time $$\omega t = 3$$, for rotational Taylor–Couette flow between concentric cylinders. Results are shown in dimensionless form, for a third-mode wave with $$q = 3$$ and amplitude parameter $$\varepsilon = 0.15$$ The density ratio is $$D = 1.2$$, the radii are $$b/a = 2$$ and the Froude number is $$F = 1$$ 4. Conclusions When two flowing inviscid fluids are separated by an interface of infinitesimal thickness, the interface is unstable. According to linearised theory, any small-amplitude disturbance to the interface will grow exponentially with time if the two fluids move with different speeds. Forbes and Cosgrove (20) recently showed that a similar situation exists in cylindrical geometry, where two adjacent fluids move under the effects of a line vortex. An interface disturbance that grows exponentially, however, eventually violates the linearisation assumption that disturbances are small, so that non-linear effects ultimately become important. Moore (2) has shown how non-linearity results in the formation of a curvature singularity at the interface, for inviscid fluids, and within finite time. The Moore curvature singularity is a subtle and complicated effect, and Moore’s asymptotic analysis (2) is a difficult calculation. Intuitively, it might be anticipated that the addition of fluid viscosity would so damp the interface that Moore’s curvature singularity now could no longer form. This is indeed the case; however, pathological behaviour at a narrow interface between viscous fluids can nevertheless still occur, in spite of the fact that there can be no velocity jump at the interface, and even when the interfacial disturbance starts from a smooth differentiable shape at $$t = 0$$. This is arguably a somewhat surprising fact and appears not to have been commented upon explicitly in the literature. Furthermore, unlike the complicated asymptotic analysis leading to Moore’s (2) result for the curvature singularity at an interface between inviscid fluids, in the viscous case a very simple analysis reveals this behaviour. For viscous flow of two fluids in a parallel-walled channel, Forbes et al. (10) demonstrated that a careful numerical solution revealed somewhat strange behaviour in the curvature of a narrow interface, and suggested that curvature singularity might be the cause. A significant aim of the present article has therefore been to explore whether similar behaviour at the interface might again be possible, and we have presented a simple analytical demonstration that, only after a certain time has elapsed, curvature spikes of rapidly-growing magnitude are indeed possible. In that sense, they appear numerically to be very similar to the genuine singularity of Moore (2) for inviscid fluids, even although the curvature magnitude in this viscous case does not become infinite within finite time. Perhaps the main aim of this present article has been to develop the somewhat novel spectral technique for solving the non-linear problem described in Section 3, and to assess the extent to which the simple asymptotic analysis of section 2 can represent the solution behaviour. There is an initially circular interface between the two fluids, and the asymptotic analysis shows that large curvatures can develop, even when the initial disturbance to the interface is quite smooth. The straightforward analysis of section 2 is sufficient to reveal this effect, and the accuracy of this asymptotic formula (2.2) has indeed been confirmed by the numerical solution of section 3. Nevertheless, there is an important difference between the spectral numerical solutions of section 3 and the analysis of section 2. The former is based on Boussinesq theory, in which the exact mathematical interface is replaced by a narrow zone of finite width, over which density changes rapidly, but smoothly. Thus, although there is close agreement with the asymptotic solutions in Fig. 3, the numerical results can nevertheless not produce the same near-singular behaviour of the asymptotic solutions. The important conclusion, then, is that the Moore singularity for inviscid fluids and the curvature spikes that develop with viscous fluids are both consequences of the assumption of an infinitessimally thin interface; the success of ‘vortex blob’ methods (3), (5) for numerical solution of the inviscid problem is not that they mimic viscosity, but rather that they create a diffuse interface of finite width, which prevents the Moore singularity from developing. Furthermore, as the curvature spikes discussed in this article are able to be predicted using quite simple arguments, it seems likely that they could occur quite generally at narrow interfaces between flowing viscous fluids. The phenomenon therefore appears worthy of further study. In addition, it may be possible to use simple results of the type presented in Section 2 as the basis of more complete asymptotic analyses. Acknowledgements We are grateful to the referees for numerous helpful suggestions that helped us clarify many aspects of the article. Discussion with Rhys A. Paul and Stephen J. Walters is gratefully acknowledged. This work is associated with Australian Research Council Discovery Grant DP140100094. References 1. Drazin P. G. and Reid W. H. Hydrodynamic Stability , 2nd edn. ( Cambridge University Press, Cambridge 2004). Google Scholar CrossRef Search ADS   2. Moore D. W. The spontaneous appearance of a singularity in the shape of an evolving vortex sheet, Proc. Roy. Soc. London A  365 ( 1979) 105– 119. Google Scholar CrossRef Search ADS   3. Krasny R. Desingularization of periodic vortex sheet roll-up, J. Comput. Phys.  65 ( 1986) 4. Chen M. J. and Forbes L. K. Accurate methods for computing inviscid and viscous Kelvin–Helmholtz instability, J. Comput. Phys.  230 ( 2011) 1499– 1515. Google Scholar CrossRef Search ADS   5. Baker G. Caflisch R. E. and Siegel M. Singularity formation during Rayleigh–Taylor instability, J. Fluid Mech.  252 ( 1993) 51– 78. Google Scholar CrossRef Search ADS   6. Forbes L. K., The Rayleigh-Taylor instability for inviscid and viscous fluids, J. Engin. Math.  65 ( 2009) 273– 290. Google Scholar CrossRef Search ADS   7. Eggers J. and Villermaux E. Physics of liquid jets, Rep. Prog. Phys.  71 ( 2008) 036601. Google Scholar CrossRef Search ADS   8. Buckmaster. J. D. Pointed bubbles in slow viscous flow, J. Fluid Mech.  55 ( 1972) 385– 400. Google Scholar CrossRef Search ADS   9. Joseph D. D. Nelson J. Renardy M. and Renardy Y. Two-dimensional cusped interfaces, J. Fluid Mech.  223 ( 1991) 383– 409. Google Scholar CrossRef Search ADS   10. Forbes L. K. Paul R. A. Chen M. J. and Horsley D. E. Kelvin–Helmholtz creeping flow at the interface between two viscous fluids, ANZIAM J.  56 ( 2015) 317– 358. Google Scholar CrossRef Search ADS   11. Pozrikidis C. Instability of two-layer creeping flow in a channel with parallel-sided walls, J. Fluid Mech.  351 ( 1997) 139– 165. Google Scholar CrossRef Search ADS   12. Li J. and Renardy Y. Numerical study of flows of two immiscible liquids at low Reynolds number, SIAM Rev.  42 ( 2000) 417– 439. Google Scholar CrossRef Search ADS   13. Mallock A. Determination of the viscosity of water, Proc. R. Soc. Lond. A  45 ( 1888) 126– 132. Google Scholar CrossRef Search ADS   14. Mallock A. Experiments on fluid viscosity, Phil. Trans. R. Soc. Lond. A  187 ( 1896) 41– 56. Google Scholar CrossRef Search ADS   15. Couette M. Études sur le frottement des liquides, Ann. Chim. Phys.  6 ( 1890) 433– 510. 16. Andereck C. D. Liu S. S. and Swinney H. L. Flow regimes in a circular Couette system with independently rotating cylinders, J. Fluid Mech.  164 ( 1986) 155– 183. Google Scholar CrossRef Search ADS   17. Marcus P. S. Simulation of Taylor–Couette flow. Part 1. Numerical methods and comparison with experiment, J. Fluid Mech.  146 ( 1984) 45– 64. Google Scholar CrossRef Search ADS   18. Marcus P. S. Simulation of Taylor–Couette flow. Part 2. Numerical results for wavy-vortex flow with one travelling wave, J. Fluid Mech.  146 ( 1984) 65– 113. Google Scholar CrossRef Search ADS   19. Peng J. and Zhu K.-Q. Linear instability of two-fluid Taylor–Couette flow in the presence of surfactant, J. Fluid Mech.  651 ( 2010) 357– 385. Google Scholar CrossRef Search ADS   20. Forbes L. K. and Cosgrove J. M., A line vortex in a two-fluid system, J. Engin. Math.  84 ( 2014), 181– 199. Google Scholar CrossRef Search ADS   21. Farrow D. E. and Hocking G. C. A numerical model for withdrawal from a two-layer fluid, J. Fluid Mech.  549 ( 2006) 141– 157. Google Scholar CrossRef Search ADS   22. Batchelor G. K. An Introduction to Fluid Mechanics. ( Cambridge University Press 1967). Google Scholar CrossRef Search ADS   23. Abramowitz M. and Stegun I. A. (eds), Handbook of Mathematical Functions . ( Dover, New York 1972). 24. Horsley D. E. Bessel phase functions: calculation and application, Numerische Math.  136 ( 2017) 679– 702. Google Scholar CrossRef Search ADS   © The Author, 2017. Published by Oxford University Press; all rights reserved. For Permissions, please email: journals.permissions@oup.com

Journal

The Quarterly Journal of Mechanics and Applied MathematicsOxford University Press

Published: Feb 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Unlimited reading

Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere.

Stay up to date

Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates.

Organize your research

It’s easy to organize your research with our built-in tools.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

Monthly Plan

  • Read unlimited articles
  • Personalized recommendations
  • No expiration
  • Print 20 pages per month
  • 20% off on PDF purchases
  • Organize your research
  • Get updates on your journals and topic searches

$49/month

Start Free Trial

14-day Free Trial

Best Deal — 39% off

Annual Plan

  • All the features of the Professional Plan, but for 39% off!
  • Billed annually
  • No expiration
  • For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588

$360/year

billed annually
Start Free Trial

14-day Free Trial