Rayleigh-Taylor instability in variable density swirling flows

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 χ 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.


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][2][3][4]. Hence, they are widely studied in the literature ( [1][2][3][4][5][6][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][9][10] (due to density variations) and Taylor-Couette [11][12][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 instabila e-mail: abid@irphe.univ-mrs.fr ity 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 threedimensionality of the flow, giving it a helical form.

Linear asymptotics
Consider an inviscid and incompressible swirling flow with a basic velocity field (in a cylindrical coordinate system (r, θ, x) with the basis (e r , e θ , e x ); (x, y, z) being the usual cartesian coordinate system), and a density field, ρ b (r). Let f (r) be the amplitude, k the axial wavenumber, m the azimuthal wavenumber and ω the complex angular frequency of a perturbation to the flow. For a helical perturbation of the form f (r) exp(i(mθ + kx − ωt)), in the limit of large wavenumbers, it is shown in [14] that the leading order (in = O(1/k) = O(1/m)) of the perturbation growth rate σ ≡ (ω) is equal to χ(m, k, r 0 ) where the function χ is defined as: and r 0 is the location where χ is extremal. In the χ expression, Ω = V/r is the local angular velocity, H is defined as: The prime, , is the derivative with respect to the r variable.
It is also shown in [14] that for the particular case (m 1, k = 0), the leading order (in = 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 1) the leading order (in = 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 = 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 = 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 = 0, k = 0) and (m = 0, k = 0) Rayleigh-Taylor and Taylor-Couette, respectively. By inspection of the χ 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 2 r 2 0 )), Taylor-Couette instability (with a weight equal to k 2 r 2 0 /(m 2 + k 2 r 2 0 )) and an instability due to the coupling of rotation and the axial shear (with a weight equal to 2mkr 0 /(m 2 + k 2 r 2 0 )). Note that the coupling term is responsible for the three-dimensionalization of the flow instability. Note also that for a fixed k and m → ∞ the function χ → G and the instability should be Rayleigh-Taylor-dominated. For, m fixed and k → ∞ the function χ → 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.

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: where ρ is the variable (in space and time) density field, Re =ρU a/μ is the Reynolds number with μ the dynamic viscosity,ρ 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 = ∞ any basic velocity (Eq. (1)), U 0 (x), and density, ρ b (r), are equilibrium solutions of the previous equations (see the appendix) provided that ∇P = ρ 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 ∇P = −ρ 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πa) 3 . The use of a pseudo-spectral algorithm in a cubic (or parallelipipedic) box, for a flow that is initially axisymmetric, is possible for flows with an exponential decay of their vorticity 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 i s (i denoting one of the two profiles) the dynamics is linear, and for t > t i s it is nonlinear. Note that for t < t i s 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.
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][22][23][24][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δt (n ∈ N), where δt is the time-step, equations (3)-(5) are solved as follows. Let P ρ = ( where I d is the identity operator, L −1 ρ is the inverse (in the space of zero mean functions) of L ρ ≡ ∇·(1/ρ∇) and ∇· is the divergence operator. Note that P ρ is an exact projector, i.e., P 2 ρ = P ρ , and that ∇ · P ρ v = 0, for any differentiable vector field v and density profile ρ (i.e., P ρ is a projector onto the space of divergence-free vector fields). The equation is advanced in time using a modified, fully explicit, second order Runge-Kutta scheme with two intermediate projection steps: Note that this time integration is chosen to be fully explicit because all the terms in F(.) are nonlinear (especially the diffusive term: (Δu n /(ρ n Re))). Note also that the projection step needs the inversion of L ρ that will be detailed when the pressure is calculated. The density is calculated with a classical second-order Runge-Kutta scheme, ≡ G(u n , ρ n ), + G(u n , ρ n + δtG(u n , ρ n )).
Consistency of the numerical scheme with variable density Navier-Stokes equations leads to the following equation for the pressure: with K * 1 and K * 2 expressed above. Noting that the modified Runge-Kutta scheme, for the velocity, needs the computation of P ρ K * 1 and P ρ K * 2 separately, equation (15) is written as Poisson equations, and, if needed, the pressure is computed as P n = (P n 1 + P n 2 )/2. The solution of equation (16), for the "partial " pressure P i , gives L −1 ρ ∇ · K * i that is necessary in the computation of P ρ K * i . To this end, equation (16) is solved iteratively. Let k (∈ N) be the current iteration index. Hence, P n i,k is computed as: and the iterative process is stopped when the residual, ||ΔP n i,k − h n i,k−1 ||, is smaller than a desired accuracy (for a given norm). Note that Δ −1 is the laplacian inverse for zero mean functions. The number of iterations needed to converge is very dependent on the initial value P n i,0 . Therefore, to accelerate slightly this process, the previous converged value of P i is used to initialize the next one: P n+1 i,0 = P n i,c , where P n i,c 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 follow: (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, t , is of order two when decreasing the time-step δt, i.e., t = O((δt) 2 ) in the limit δt → 0; (c) we compute the growth rates of flow perturbations at early times (when nonlinearities are small: t t i s , 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: where 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 en-ergy growth-rate obtained using the numerical solution of Navier-Stokes equations (the dot-dashed line) is: 0.66 for t < t Exp s . 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: 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: In what follows the length scale is chosen as the vortexcore length a, the velocity scale as the axial velocity at r = 0 (U = W (0)) and the density scale as the density at infinityρ = ρ ∞ ≡ ρ b (r → ∞). Two, new, nondimensional numbers could be introduced: the Atwood number At = (ρ b (0)−ρ ∞ )/(ρ b (0)+ρ ∞ ) = (s−1)/(s+1), and the Taylor number T a = Ω 2 a 4 /(μ/ρ ∞ ) 2 = q 2 Re 2 (with Ω = 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 (ω 0 = ∇ × U 0 ) has an exponential decay. Therefore, we compute the velocity as: which is automatically periodic if the operator ∇ 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 a and L z 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 14a. Two resolutions (192 3 and 256 3 ) are used to ensure convergence. For azimuthal modes, the characteristic scale is 2πa (since they grow near the vortex core, r = a, where the gradients are high). Therefore, with lateral dimensions L y = L z 14a, we have (2πa/14a)256 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 22. This is enough for the physics we are considering (instability, transition to turbulence and not turbulence). 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 constantdensity swirling-jet experiement [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 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 χ), we find that χ 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 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 ∼ 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 χ is dominated by the coupling term Ω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 χ 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.

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 χ 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.   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.