# Rayleigh-Taylor instability in variable density swirling flows

Bastien Di Pierro, Malek Abid

► To cite this version:

Bastien Di Pierro, Malek Abid. Rayleigh-Taylor instability in variable density swirling flows. The European Physical Journal B: Condensed Matter and Complex Systems, 2012, 85 (69), pp.1-8. 10.1140/epjb/e2012-20540-6 . hal-00701221

HAL Id: hal-00701221

https://hal.science/hal-00701221

Submitted on 29 Mar 2021

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Rayleigh-Taylor instability in variable density swirling flows

B. Dipierro and M. Abidª

IRPHE – UMR 6594, Technopole de Château-Gombert, 49 rue Joliot Curie, B.P. 146, 13384 Marseille Cedex 13, France

**Abstract.** Using linear instability theory and nonlinear dynamics, the Rayleigh-Taylor instability of variable density swirling flows is studied. It is found that the flow topology could be predicted, when the instability sets in, using a function $\chi$ dependent on density and axial and azimuthal velocities. It is shown that even when the inner axial-flow is heavier than the outer one (a favorable case for the development of the Rayleigh-Taylor instability thanks to the centrifugal force) the instability is not necessarily Rayleigh-Taylor-dominated. It is also shown that when the Rayleigh-Taylor instability develops, it is helical.

## 1 Introduction

Swirling flows are important for many application devices and as a fundamental problem considering their relevance to aircraft trailing vortices, vortical transport of momentum and energy in meteorology, and vortex bursting (or breakdown) [1–4]. Hence, they are widely studied in the literature ([1–7], among many others) but the physical mechanisms of their instabilities, when density variations are present, are not generally discussed.

In a variable density swirling flow (for example a columnar vortex with an axial flow, a configuration that will receive a particular attention in this study) there are ingredients for the development of two fundamental instabilities: Rayleigh-Taylor [8–10] (due to density variations) and Taylor-Couette [11–13] (promoted by differential rotation). Intuitively, the instability of the flow is a mixing of Rayleigh-Taylor and Taylor-Couette instabilities. Therefore, a natural question (from a physical point of view) is how to formulate the instability problem to account for the contribution of each of these latter instabilities to the general, three-dimensional, instability of a variable density swirling flow? Note that the key instability dynamics of Rayleigh-Taylor and Taylor-Couette can both be represented with a minimum of two dimensions. The former is characterized by “bubbles” for positively buoyant fluid and “spikes” for negatively buoyant fluid, the latter being axisymmetric with the organization of the flow as toroidal patterns.

The objective of the present study is to answer the above question based on the recent asymptotic analysis by [14]. To this end, since Taylor experiments [11] of a fluid between two rotating cylinders have indicated that the disturbance could be axisymmetric when the instability

occurs, let us define two particular instabilities: an instability with an axisymmetric centrifugal-force that could be identified with an axisymmetric Taylor-Couette instability, and a two-dimensional “orthogonal state” to the previous one, with no axisymmetry of the centrifugal force, that could be identified with a two-dimensional Rayleigh-Taylor instability (see below for more precision about the orthogonality). We show that the nonlinear evolution of the three-dimensional instability is a mixing of these two (pure) states with a coupling term of axial and azimuthal velocities. This coupling term is responsible for the three-dimensionality of the flow, giving it a helical form.

## 2 Linear asymptotics

Consider an inviscid and incompressible swirling flow with a basic velocity field (in a cylindrical coordinate system $(r, \theta, x)$ with the basis $(\mathbf{e}_r, \mathbf{e}_\theta, \mathbf{e}_x$); $(x, y, z)$ being the usual cartesian coordinate system),

$$ \mathbf{U}_0(\mathbf{x}) = (0, V(r), W(r))_{(\mathbf{e}_r, \mathbf{e}_\theta, \mathbf{e}_x)}, \quad r = \sqrt{y^2 + z^2} \quad (1) $$

and a density field, $\rho_b(r)$. Let $f(r)$ be the amplitude, $k$ the axial wavenumber, $m$ the azimuthal wavenumber and $\omega$ the complex angular frequency of a perturbation to the flow. For a helical perturbation of the form $f(r) \exp(i(m\theta + kx - \omega t))$, in the limit of large wavenumbers, it is shown in [14] that the leading order (in $\epsilon = O(1/k) = O(1/m)$) of the perturbation growth rate $\sigma \equiv \Im(\omega)$ is equal to $\chi(m, k, r_0)$ where the function $\chi$ is defined as:

$$ \chi^2(m, k, r) = \frac{m^2 G^2 + k^2 r^2 H^2 + 2mkr\Omega W'}{m^2 + k^2 r^2}, \quad (2) $$

<footnote>ª e-mail: abid@irphe.univ-mrs.fr</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** (a) A schematic of a cylindrical Rayleigh-Taylor instability. The flow acceleration, $\Gamma$, is oriented outward from the vortex core, and if the vortex is heavier than the ambient fluid (i.e. $\nabla \rho$ is oriented from the outside to the vortex core), this leads to a Rayleigh-Taylor unstable configuration. (b) Axisymmetric Taylor-Couette instability.</figcaption></figure>

and $r_0$ is the location where $\chi$ is extremal. In the $\chi$ expression, $\Omega = V/r$ is the local angular velocity, $H$ is defined as: $H^2 = G^2 - \Phi$, $G$ (the analogue of the Brunt-Väisäla frequency) is defined as: $G^2 = -(V^2/r)\rho_b'/\rho_b$ and $\Phi = r^{-3}d(r^2V^2)/dr = 2\Omega\Theta$ (the Rayleigh discriminant). The prime, $'$, is the derivative with respect to the $r$ variable.

It is also shown in [14] that for the particular case ($m \gg 1, k = 0$), the leading order (in $\epsilon = O(1/m)$) of the growth rate is equal to $G(r_0)$ ($r_0$ is the location where $G$ is extremal), and for the particular case ($m = 0, k \gg 1$) the leading order (in $\epsilon = O(1/k)$) of the growth rate is equal to $H(r_0)$ ($r_0$ is the location where $H$ is extremal). In the present paper we note that these two particular cases could be related to two classically well known instabilities. In the case with ($m \neq 0, k = 0$), the perturbation grows only in the azimuthal direction leading eventually to the formation of finger-like structures of heavy fluid penetrating the light one, as explained and sketched in Figure 1a. This is the Rayleigh-Taylor instability where the centrifugal acceleration replaces the gravity. In the case with ($m = 0, k \neq 0$), the perturbation grows axisymmetrically, and is wavy in the axial direction, leading to the organization of the flow in toroidal structures characteristic of the Taylor-Couette instability (Fig. 1b). Therefore, we will call the two orthogonal eigenstates ($m \neq 0, k = 0$) and ($m = 0, k \neq 0$) Rayleigh-Taylor and Taylor-Couette, respectively. By inspection of the $\chi$ expression (Eq. (2)), clearly the three-dimensional instability of a variable density swirling flow is a mixing of both Rayleigh-Taylor instability (with a weight equal to $m^2/(m^2 + k^2r_0^2)$), Taylor-Couette instability (with a weight equal to $k^2r_0^2/(m^2 + k^2r_0^2)$) and an instability due to the coupling of rotation and the axial shear (with a weight equal to $2mkr_0/(m^2 + k^2r_0^2)$). Note that the coupling term is responsible for the three-dimensionalization of the flow instability. Note also that for a fixed $k$ and

$m \to \infty$ the function $\chi \to G$ and the instability should be Rayleigh-Taylor-dominated. For, $m$ fixed and $k \to \infty$ the function $\chi \to H$ and the instability should be Taylor-Couette-dominated.

The linear asymptotical-results of Di Pierro and Abid [14] are valid for large azimuthal and axial wavenumbers. It will be shown here, numerically, that they work well even for low values of $m$ and $k$ and that they are also valid in the nonlinear regime.

## 3 Nonlinear study

We now study the instability behavior beyond the point that linearizing assumptions are valid, and use numerical simulation to study this non-linear regime. The governing equations are:

$$ \frac{\partial \rho}{\partial t} + \mathbf{u} \cdot \nabla \rho = 0, \quad (3) $$

$$ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} + \frac{\nabla P}{\rho} - \frac{1}{Re\rho}\Delta\mathbf{u} = 0, \quad (4) $$

$$ \nabla \cdot \mathbf{u} = 0, \quad (5) $$

where $\rho$ is the variable (in space and time) density field, $Re = \bar{\rho}Ua/\mu$ is the Reynolds number with $\mu$ the dynamic viscosity, $\bar{\rho}$ a density scale (chosen here as the density at infinity), $U$ a characteristic velocity (chosen as the velocity on the axis) and $a$ the characteristic length scale (chosen as the vortex-core length). These equations are the full variable density Navier-Stokes equations and there is no Boussinesq approximation [15,16]. Therefore large density variations are allowed. To appreciate the analogy with a Rayleigh-Taylor system note that for $Re = \infty$ any basic velocity (Eq. (1)), $U_0(x)$, and density, $\rho_b(r)$, are equilibrium solutions of the previous equations (see the appendix) provided that $\nabla P = \rho_b(V^2/r)e_r$ (to be compared with an equilibrium of a planar density-stratification in the gravity field, $g = -ge_x$, for which $\nabla P = -\rho_b ge_x$).

To study this possible competition, we solve the variable density, constant dynamical viscosity, Navier-Stokes equations using numerical simulations. The time is made nondimensional using the time scale $a/U$. We use periodic boundary conditions and the initial condition is a superposition of the basic flow (Eq. (1)) and a perturbation flow obtained as an eigenmode of the linearized Navier-Stokes equations (about the basic flow) [17]. We seek the nonlinear behavior of the solution (after the time of the saturation of the linear instability, noted $t_s$ and defined in Fig. 2). The solution is obtained numerically using a pseudo-spectral (Fourier) method for the spatial discretization and a second order Runge-Kutta scheme for the time discretization. Aliased modes are removed using the “2/3” rule [18,19]. The computation is solved in a cartesian domain with a volume equal to $(4\pi a)^3$. The use of a pseudo-spectral algorithm in a cubic (or parallelepipedic) box, for a flow that is initially axisymmetric, is possible for flows with an exponential decay of their vorticity
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2. Plot against time of the energy, $E_{mk}$, of the perturbation modes a) ($m=1, k=2$) for the experimental profile (Exp) and b) ($m=3, k=1$) for the Batchelor-q-Vortex (BqV). Solid-line: the linear instability-theory prediction. Dot-dashed line: the energy growth obtained using the numerical solution of the variable density Navier-Stokes equations. The saturation time is defined as the time when the exponential growth, predicted by linear instability theory, is no longer valid due to nonlinearities of Navier-Stokes equations: for $t < t_s^i$ (i denoting one of the two profiles) the dynamics is linear, and for $t > t_s^i$ it is nonlinear. Note that for $t < t_s^i$ the growth rate of the perturbation predicted using linear instability theory (solid line) is well recovered in the numerical simulation (dot-dashed line): this is a validation of the numerical procedure.</figcaption></figure>

with the radial distance. The advantages of this method (compared to a method using a cylindrical coordinate system) are: (a) this pseudo-spectral method is easy to implement and parallelize; (b) there is no singularity at the origin of the coordinate system, hence there is no need to implement the pole conditions (or equations) [20]. In the literature, many studies have been done with a pseudo-

spectral algorithm in a parallelipipedic box, for a flow that is initially axisymmetric, principally for jets, swirling jets and wakes [21–25]. Note that these numerical studies are done in the constant density case. In the present work, we extended the numerical method to the variable density case, as described below.

Using standard notations, at time $t_n = n\delta t$ ($n \in \mathbb{N}$), where $\delta t$ is the time-step, equations (3)–(5) are solved as follows. Let $\mathbb{P}_\rho = (\mathbb{I}_d - (1/\rho)\nabla\mathbb{L}_\rho^{-1}\nabla\cdot)$ where $\mathbb{I}_d$ is the identity operator, $\mathbb{L}_\rho^{-1}$ is the inverse (in the space of zero mean functions) of $\mathbb{L}_\rho \equiv \nabla\cdot(1/\rho\nabla)$ and $\nabla\cdot$ is the divergence operator. Note that $\mathbb{P}_\rho$ is an exact projector, i.e., $\mathbb{P}_\rho^2 = \mathbb{P}_\rho$, and that $\nabla\cdot\mathbb{P}_\rho\mathbf{v} = 0$, for any differentiable vector field $\mathbf{v}$ and density profile $\rho$ (i.e., $\mathbb{P}_\rho$ is a projector onto the space of divergence-free vector fields). The equation

$$ \left(\frac{\partial \mathbf{u}}{\partial t}\right)^n = -(\mathbf{u}^n \cdot \nabla)\mathbf{u}^n + \frac{1}{Re\rho^n} \Delta \mathbf{u}^n - \frac{\nabla P^n}{\rho^n} \\ \equiv \mathbf{F}(\mathbf{u}^n, \rho^n) - \frac{\nabla P^n}{\rho^n} \qquad (6) $$

is advanced in time using a modified, fully explicit, second order Runge-Kutta scheme with two intermediate projection steps:

$$ \mathbf{K}_1^* = \mathbf{F}(\mathbf{u}^n, \rho^n), \qquad (7) $$

$$ \mathbf{K}_1^n = \mathbb{P}_\rho \mathbf{K}_1^*, \qquad (8) $$

$$ \mathbf{K}_2^* = \mathbf{F}(\mathbf{u}^n + \delta t \mathbf{K}_1^n, \rho^n), \qquad (9) $$

$$ \mathbf{K}_2^n = \mathbb{P}_\rho \mathbf{K}_2^* \qquad (10) $$

$$ \mathbf{u}^{n+1} = \mathbf{u}^n + \frac{1}{2} \delta t (\mathbf{K}_1^n + \mathbf{K}_2^n). \qquad (11) $$

Note that this time integration is chosen to be fully explicit because all the terms in $\mathbf{F}(\cdot)$ are nonlinear (especially the diffusive term: $(\Delta \mathbf{u}^n/(\rho^n Re))$). Note also that the projection step needs the inversion of $\mathbb{L}_\rho$ that will be detailed when the pressure is calculated.

The density is calculated with a classical second-order Runge-Kutta scheme,

$$ \left(\frac{\partial \rho}{\partial t}\right)^n = -\mathbf{u}^n \cdot \nabla \rho^n \qquad (12) \\ \equiv G(\mathbf{u}^n, \rho^n), $$

$$ \rho^{n+1} = \rho^n + \frac{1}{2}\delta t(G(\mathbf{u}^n, \rho^n) + G(\mathbf{u}^n, \rho^n + \delta t G(\mathbf{u}^n, \rho^n))). \qquad (14) $$

Consistency of the numerical scheme with variable density Navier-Stokes equations leads to the following equation for the pressure:

$$ \mathbb{L}_\rho P^n = \nabla \cdot \left( \frac{\nabla P^n}{\rho^n} \right) = \frac{1}{2} \nabla \cdot (\mathbf{K}_1^* + \mathbf{K}_2^*), \qquad (15) $$
with $\mathbf{K}_1^*$ and $\mathbf{K}_2^*$ expressed above. Noting that the modified Runge-Kutta scheme, for the velocity, needs the computation of $\mathbb{P}_{\rho}\mathbf{K}_1^*$ and $\mathbb{P}_{\rho}\mathbf{K}_2^*$ separately, equation (15) is written as Poisson equations,

$$ \Delta P_i^n = \frac{\nabla \rho^n \cdot \nabla P_i^n}{\rho^n} + \rho^n \nabla \cdot \mathbf{K}_i^* \equiv h_i^n, \quad i = 1, 2, \quad (16) $$

and, if needed, the pressure is computed as $P^n = (P_1^n + P_2^n)/2$. The solution of equation (16), for the “partial” pressure $P_i$, gives $\mathbb{L}_{\rho}^{-1}\nabla \cdot \mathbf{K}_i^*$ that is necessary in the computation of $\mathbb{P}_{\rho}\mathbf{K}_i^*$. To this end, equation (16) is solved iteratively. Let $k \in \mathbb{N}$ be the current iteration index. Hence, $P_{i,k}^n$ is computed as:

$$ \begin{aligned} P_{i,k}^n &= \Delta^{-1} \left( \frac{\nabla \rho^n \cdot \nabla P_{i,k-1}^n}{\rho^n} + \rho^n \nabla \cdot \mathbf{K}_i^* \right) \\ &\equiv \Delta^{-1} h_{i,k-1}^n, \end{aligned} \quad (17) $$

and the iterative process is stopped when the residual, $||\Delta P_{i,k}^n - h_{i,k-1}^n||$, is smaller than a desired accuracy (for a given norm). Note that $\Delta^{-1}$ is the laplacian inverse for zero mean functions. The number of iterations needed to converge is very dependent on the initial value $P_{i,0}^n$. Therefore, to accelerate slightly this process, the previous converged value of $P_i$ is used to initialize the next one: $P_{i,0}^{n+1} = P_{i,c}^n$, where $P_{i,c}^n$ is the converged “partial” pressure at the previous time. In practice, we find that six iterations are sufficient to obtain a velocity divergence $O(10^{-12})$.

The overall numerical procedure is validated as follows: (a) we check the spectral accuracy by verifying the exponential decay of the power spectrum density (at high wavenumbers) of each computed field (velocity, pressure and density); (b) we check that the numerical time-error, $\epsilon_t$, is of order two when decreasing the time-step $\delta t$, i.e., $\epsilon_t = O((\delta t)^2)$ in the limit $\delta t \to 0$; (c) we compute the growth rates of flow perturbations at early times (when nonlinearities are small: $t \ll t_s^i$, $i$ denoting one of the two profiles used, see below and Fig. 2) and compare them to growth rates obtained using diagonalization of linearized Navier-Stokes equations [17]: we check that differences between growth rates are less than 5%. For example, for the radial velocity (which is initially zero), $u_r$, computed using the Navier-Stokes equations, let us define the energy in a helical structure as:

$$ E_{mk}(t) = \left| \int_{\mathcal{D}} d\mathbf{x} u_r^2(\mathbf{x}, t) e^{-i(kx+m\theta)} \right|, \quad (18) $$

where $\mathcal{D}$ is the volume of the computational box. In Figure 2, the energy $E_{mk}$ of the mode ($m = 1, k = 2$) is plotted as a function of time (for $s = 2$ and using the experimental profile as a basic flow, see below for its definition). The energy growth-rate predicted using the linear instability theory (the continuous line) is 0.67. The energy growth-rate obtained using the numerical solution of Navier-Stokes equations (the dot-dashed line) is: 0.66 for $t < t_s^{exp}$.

We study the competition between Rayleigh-Taylor and Taylor-Couette instabilities for two particular flows. The first one is the well known Batchelor $q$ vortex (Bqv) [26], that we extend to the variable density case:

$$ \begin{aligned} V(r)/W(0) &= q \frac{1 - \exp(-(r/a)^2)}{(r/a)}, \\ W(r)/W(0) &= \exp(-(r/a)^2), \\ \rho_b(r)/\rho_\infty &= 1 + (s-1) \exp(-(r/a)^2), \end{aligned} \quad (19) $$

where $q$ is the swirl number and $s$ is the density ratio between the vortex center and the ambient fluid. The second one is a fit of experimentally obtained vortices with axial flows [27], that we extend also to the variable density case:

$$ \begin{aligned} V(r)/W(0) &= q(r/a) \exp(-(r/a)^4), \\ W(r)/W(0) &= \frac{1}{1 + (\exp((r/a)^2 \log(2)) - 1)^3}, \\ \rho_b(r)/\rho_\infty &= 1 + \frac{s-1}{1 + (\exp((r/a)^2 \log(2)) - 1)^3}. \end{aligned} \quad (20) $$

In what follows the length scale is chosen as the vortex-core length $a$, the velocity scale as the axial velocity at $r = 0$ ($U = W(0)$) and the density scale as the density at infinity $\bar{\rho} = \rho_\infty \equiv \rho_b(r \to \infty)$. Two, new, nondimensional numbers could be introduced: the Atwood number $At = (\rho_b(0) - \rho_\infty)/(\rho_b(0)+\rho_\infty) = (s-1)/(s+1)$, and the Taylor number $Ta = \Omega^2 a^4/(\mu/\rho_\infty)^2 = q^2 Re^2$ (with $\Omega = qU/a$ the angular-velocity scale).

The azimuthal velocity of the Bqv has an algebraic decay. To accommodate this velocity profile in the periodic domain of computation we remark that it is divergence free and its vorticity ($\omega_0 = \nabla \times U_0$) has an exponential decay. Therefore, we compute the velocity as: $U_0 = -\Delta^{-1}\nabla \times \omega_0$, which is automatically periodic if the operator $\nabla$ is the gradient on the space of periodic functions.

In practice, if $L_y$ and $L_z$ are the lateral dimensions of the computational box, we choose $L_y \gg a$ and $L_z \gg a$ (see Fig. 3) to ensure that the vortex neighbors, introduced by periodicity, do not affect the growth rates of the linear instability. We used $L_x = L_y = L_z \simeq 14a$. Two resolutions (192³ and 256³) are used to ensure convergence. For azimuthal modes, the characteristic scale is $2\pi a$ (since they grow near the vortex core, $r = a$, where the gradients are high). Therefore, with lateral dimensions $L_y = L_z \simeq 14a$, we have $(2\pi a/14a)256 \simeq 114$ points. If we consider, for example, that 5 points per wavelength are sufficient to accurately represent a wave [28], we are able to accurately represent azimuthal modes up to $m \simeq 22$. This is enough for the physics we are considering (instability, transition to turbulence and not turbulence).
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Geometry of the computational box in the lateral plane $(y, z)$. The lateral lengths are chosen such that $L_y \gg a$ and $L_z \gg a$ to ensure that the vortex neighbors, introduced by the periodization process, do not affect the growth rates of the linear instability.</figcaption></figure>

The reynolds number is chosen as $Re = 1000$, that is (initially) the self-advection is 1000 times more important than the viscous diffusion. This value of $Re$ is selected as the highest possible value that retains numerical stability, at the used resolutions. It is interesting to mention that this value of Reynolds number is close to ones used in a variable-density jet experiment [29], and a constant-density swirling-jet experiment [27].

Note that these profiles are exact solutions of Euler equations. When these profiles are used as initial conditions for the Navier-Stokes equations, the corresponding velocity profiles diffuse (see the appendix). For example, when $s = 1$, the axial velocity of the Bqv diffuses such that $a = a_0\sqrt{1 + 4t/Re}$, where $a_0$ is the initial vortex core length (indeed the axisymmetric diffusion equation has a solution of the form $\exp(-(r/a)^2)/a$ provided that $a(da/dt) = 2/Re$). Note also that the azimuthal velocity of the Bqv is stable to axisymmetric perturbations, however that of experimental profiles is unstable to this kind of perturbations [30,31].

For the Bqv, in the linear framework (using the function $\chi$), we find that $\chi$ is dominated by $G$ if $s > 1$ for all $m, k$ and $q$. Furthermore, the coupling term is not identically zero. Therefore, we expect the development of a three-dimensional Rayleigh-Taylor instability and the formation of helical spikes in the nonlinear regime. This is confirmed by the results presented, for example, in Figure 4. The evolution of the Bqv (for $s = 2$ and $q = -0.5$), initially disturbed with the mode $k = 1$ and $m = 3$, is

shown. At the beginning, three helical fingers grew radially. The heavy fluid, initially localized in the vortex core, is gradually transported far from the axis, owing to these helical structures. When nonlinearities become important, fingers roll-up into spikes and at the end of the simulation, the heavy fluid is localized around $r \sim 2$. The simulation is repeated with other unstable $(m, k)$ modes, and the same flow topology is observed. Hence, the helical Rayleigh-Taylor instability is dominant for the Bqv with $s > 1$. It is interesting to mention that for higher values of Re, vorticity would be greater and one would expect tighter roll-up of spikes. Therefore, for the late evolution of this roll-up a time-varying Reynolds number, based on the evolving length-scale of the spikes, would be more appropriate.

For $s < 1$, for the Bqv, we find that $\chi$ is dominated by the coupling term $\Omega W'$. The topology of the flow in the nonlinear regime is different from the previously described Rayleigh-Taylor instability. It is shown in Figure 5 for $s = 0.5$ and $q = -0.5$ and for the same mode ($k = 1$ and $m = 3$), for example. Helical structures appear and there is no spikes formation.

For the experimental profile we find that the function $\chi$ is principally dominated by $H$ for all $q, s, m$ and $k$ values. We therefore expect that in the nonlinear regime the flow develops a topology similar to that resulting from the Taylor-Couette instability. This is confirmed in Figure 6. We see that disturbances, even in the nonlinear stage, are concentrated near $r = 1$, as if the flow was bounded (note that a confined dynamics may be also observed when important density stratifications are present in Rayleigh-Taylor systems [32,33]). Furthermore, although the fluid topology is helical it presents parts of toroidal structures similar to Taylor vortices. Note that even if $s > 1$ (a favorable case for the development of the Rayleigh-Taylor instability), the instability is not Rayleigh-Taylor dominated.

## 4 Conclusion

We now summarize our main findings. The three-dimensional instability of swirling flows, submitted to density variations, is formulated as a competition among three distinct instabilities: Rayleigh-Taylor, Taylor-Couette and an instability due to the coupling between the rotation and the axial shear. We have presented the results of a linear asymptotic analysis, valid for incompressible inviscid fluids in the limit of large wavenumbers. We have demonstrated here that the analysis is also valid for nonlinear development at low wavenumbers, as confirmed by numerical simulations. It is found that the flow topology could be predicted, when the instability sets in, using the function $\chi$ dependent on density and axial and azimuthal velocities. As an application it is shown that even when the inner axial-flow is heavier than the outer one (a favorable case for the development of the Rayleigh-Taylor instability thanks to the centrifugal force) the instability is not necessarily Rayleigh-Taylor-dominated.
<figure><img src="image_6.png" /><figcaption>Fig. 4. Nonlinear evolution of a variable density Bqv, with $q = -0.5$ ($Ta = 2.5 \times 10^5$) and $s = 2$ ($At = 0.33$). The initial “interface” disturbance has $k = 1$, $m = 3$ and an initial amplitude $O(10^{-4})$ times that of the initial velocity. (a)–(d) Contour plot of the density at $x = 0$ for $t/t_s^{Bqv} = 0, 0.6, 1.1$ and $1.8$ respectively ($t_s^{Bqv}$ is the linear-instability time-saturation). (d) Isosurface of the density at a medium level for $t/t_s^{Bqv} = 1.5$. Three helical spikes are ejected radially. Only the central part (near the vortex) of the computational box is shown. Note that spatial scales are in unit of $a$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Same conditions as in Figure 4 but with $s = 0.5$ ($At = -0.33$). Helical unstable structures due to the coupling of rotation and axial shear appear. Only the central part (near the vortex) of the computational box is shown.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Nonlinear evolution of a variable density experimental profile with $q = -0.5$ ($Ta = 2.5 \times 10^5$) and $s = 2$ ($At = 0.33$). The initial “interface” disturbance has $k = 2$, $m = 1$ and an initial amplitude $O(10^{-4})$ times that of the initial velocity. (a)–(d) Contour plot of the density at $y = 0$, for $t/t_s^{Exp} = 0, 0.5, 1$ and $2$ respectively ($t_s^{Exp}$ is the linear-instability time-saturation). (e) Isosurface of the density at a medium level for $t/t_s^{Exp} = 1.5$. Although the fluid topology is helical, it presents parts of toroidal structures. Note, here, that even if $s > 1$ (a favorable case for the development of the Rayleigh-Taylor instability), the instability is not Rayleigh-Taylor-dominated. Only the central part (near the vortex) of the computational box is shown. Recall that spatial scales are in unit of a.</figcaption></figure>

## Appendix A

In a cylindrical coordinate system $(r, \theta, x)$ the variable density Navier-Stokes equations for velocity, $\mathbf{u} = (u_r, u_\theta, u_x)$, and density, $\rho$, are:

$$ \frac{\partial \rho}{\partial t} + (\mathbf{u} \cdot \nabla)\rho = 0, \qquad (A.1) $$

$$ \frac{\partial u_r}{\partial t} + (\mathbf{u} \cdot \nabla)u_r - \frac{u_\theta^2}{r} = -\frac{1}{\rho} \frac{\partial P}{\partial r} + \frac{1}{Re\rho} \left( \Delta u_r - \frac{2}{r^2} \frac{\partial u_\theta}{\partial \theta} - \frac{u_r}{r^2} \right), \qquad (A.2) $$

$$ \frac{\partial u_\theta}{\partial t} + (\mathbf{u} \cdot \nabla)u_\theta + \frac{u_r u_\theta}{r} = -\frac{1}{\rho r} \frac{\partial P}{\partial \theta} + \frac{1}{Re\rho} \left( \Delta u_\theta + \frac{2}{r^2} \frac{\partial u_r}{\partial \theta} - \frac{u_\theta}{r^2} \right), \qquad (A.3) $$

$$ \frac{\partial u_x}{\partial t} + (\mathbf{u} \cdot \nabla)u_x = -\frac{1}{\rho} \frac{\partial P}{\partial x} + \frac{1}{Re\rho} \Delta u_x, \qquad (A.4) $$

$$ \frac{1}{r} \frac{\partial (ru_r)}{\partial r} + \frac{1}{r} \frac{\partial u_\theta}{\partial \theta} + \frac{\partial u_x}{\partial x} = 0, \qquad (A.5) $$

where

$$ (\mathbf{u} \cdot \nabla) f = u_r \frac{\partial f}{\partial r} + \frac{u_\theta}{r} \frac{\partial f}{\partial \theta} + u_x \frac{\partial f}{\partial x}, $$

and

$$ \Delta f = \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial f}{\partial r} \right) + \frac{1}{r^2} \frac{\partial^2 f}{\partial \theta^2} + \frac{\partial^2 f}{\partial x^2}. $$

For any stationary (axisymmetric) velocity, $\mathbf{U} = (u_r = 0, u_\theta = V(r), u_x = W(r))$, and density, $\rho = \rho_b(r)$, profiles the following terms vanish:

$$ \frac{\partial u_r}{\partial t}, \frac{\partial u_\theta}{\partial t}, \frac{\partial u_x}{\partial t}, (\mathbf{u} \cdot \nabla)u_r, (\mathbf{u} \cdot \nabla)u_\theta, (\mathbf{u} \cdot \nabla)u_x, \frac{u_r u_\theta}{r}, $$

and the continuity equation (A.5) is automatically verified.

Therefore $(\mathbf{U}, \rho_b)$ are exact solutions of Euler equations ($Re = \infty$) provided that:

$$ \frac{V^2}{r} = \frac{1}{\rho_b} \frac{\partial P}{\partial r}. $$

In the constant density case (choose for example $\rho = 1$), if the swrling-flow velocity $\mathbf{U}$ is used as an initial condition for the Navier-Stokes equations, the velocity components $u_\theta$ and $u_x$ will diffuse according to:

$$ \frac{\partial u_\theta}{\partial t} = \frac{1}{Re} \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u_\theta}{\partial r} \right), \quad u_\theta(r, t=0) = V(r), $$

and

$$ \frac{\partial u_x}{\partial t} = \frac{1}{Re} \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u_x}{\partial r} \right), \quad u_x(r, t=0) = W(r). $$

The pressure is given by:

$$ \frac{u_\theta^2}{r} = \frac{\partial P}{\partial r}, $$

and there is no generation of radial velocity.
In the numerical simulations, we checked that when a swirling jet is used as an initial condition there is no generation of azimuthal modes *m* = 4 and *m* = 8 by periodicity. That is to say the swirling jet only diffuses and there is no generation of radial velocity by periodicity. This is also a validation that the values of lateral lengths of the computational box are sufficient and the ghost swirling-flows have no effects on the flow of the computational box.

We would like to acknowledge valuable comments of anonymous referees that ameliorated the presentation of the present work. The "Centre Informatique National de l'enseignement supérieur" (CINES) is also acknowledged for computational facility.

## References

1. L.N. Howard, A.S. Gupta, J. Fluid Mech. **14**, 463 (1962)

2. M. Lessen, P.W. Singh, F. Paillet, J. Fluid Mech. **63**, 753 (1974)

3. M. Lessen, F. Paillet, J. Fluid Mech. **65**, 769 (1974)

4. K. Stewartson, Phys. Fluids **25**, 1953 (1982)

5. K.S. Eckhoff, L. Storesletten, J. Fluid Mech. **89**, 401 (1978)

6. S. Leibovich, K. Stewartson, J. Fluid Mech. **126**, 335 (1983)

7. S. Leblanc, A. LeDuc, J. Fluid Mech. **537**, 433 (2005)

8. J.W.S. Rayleigh, Proc. Lond. Math. Soc. **14**, 170 (1883)

9. G. Taylor, Proc. R. Soc. Lond. A **201**, 192 (1950)

10. D.H. Sharp, Physica D **12**, 3 (1984)

11. G.I. Taylor, Philos. Trans. R. Soc. Lond. **223**, 289 (1923)

12. C.D. Andereck, S.S. Liu, H.L. Swinney, J. Fluid Mech. **164**, 155 (1986)

13. J. Abshagen, A. Schultz, G. Pfister, Lecture Notes in Physics **476**, 63 (1996)

14. B.D. Pierro, M. Abid, Phys. Rev. E **82**, 046312 (2010)

15. J.P. Bell, D.L. Marcus, J. Comput. Phys. **101**, 334 (1992)

16. A.S. Almgren, J.B. Bell, P. Colella, L.H. Howell, M.L. Welcome, J. Comput. Phys. **142**, 1 (1998)

17. M. Abid, J. Fluid Mech. **605**, 19 (2008)

18. G.S. Patterson Jr., S.A. Orszag, Phys. Fluids **14**, 2538 (1971)

19. S.A. Orszag, J. Atmos. Sci. **28**, 1074 (1971)

20. K. Mohseni, T. Colonius, J. Comput. Phys. **157**, 787 (2000)

21. M. Abid, M.E. Brachet, Phys. Fluids A **5**, 2582 (1993)

22. P. Brancher, J.M. Chomaz, P. Huerre, Phys. Fluids **6**, 1768 (1994)

23. M. Abid, M.E. Brachet, Phys. Fluids **10**, 469 (1998)

24. I. Delbende, J.M. Chomaz, P. Huerre, J. Fluid Mech. **335**, 229 (1998)

25. I. Delbende, M. Rossi, Phys. Fluids **17**, 044103 (2005)

26. G.K. Batchelor, J. Fluid Mech. **20**, 645 (1964)

27. F. Gallaire, J.M. Chomaz, J. Fluid Mech. **494**, 223 (2003)

28. L.N. Threfethen, *Spectral Methods in Matlab* (Siam, 2000), pp. 66–67

29. S. Boujemaa, M. Amielh, M. Chauve, C.R. Acad. Sci. **332**, 933 (2004)

30. L.D. Landau, E.M. Lifshitz, *Fluid Mechanics* (Pergamon Press, 1987), p. 99

31. B.J. Bayly, Phys. Fluids **31**, 56 (1988)

32. J.W. Jacobs, S.B. Dalziel, J. Fluid Mech. **542**, 251 (2005)

33. A. Scagliarini, L. Biferale, M. Sbragaglia, K. Sugiyama, F. Toschi, Phys. Fluids **22**, 055101 (2010)