Multiplicity of states in Taylor-Couette flow of a dense granular gas

Molecular dynamics simulations with a purely repulsive Lennard-Jones potential and a normal damping force is used to simulate the granular flow in the annular region between two differentially-rotating cylinders, called the Taylor-Couette flow. The flow transition from the azimuthally-invariant Circular Couette flow (CCF) to the Taylor-vortex flow (TVF) is studied by increasing the rotation rate (ωi) of the inner cylinder, with the outer cylinder being kept stationary. Multiplicity of states, highlighting the hysteretic nature of the “CCF ↔ TVF” transition, is observed over a wide range of rotation rates. The onset of Taylor vortices is quantified in terms of the maximum radial velocity and the net circulation per vortex.


Introduction
The flow between two concentric, independently rotating cylinders, called the Taylor-Couette flow (TCF), has attracted considerable attention since the foundational work of Couette [1] and Mallock [2] which was followed up by Taylor's seminal linear stability analysis [3] on incompressible TCF, confirming the genesis of Taylor vortices beyond a minimum rotation speed of the inner cylinder. Experiments of Coles [4] showed the multiplicity of patterns at specified rotation-rates of the inner cylinder, and ascribed this phenomena to initial conditions and the acceleration protocol employed to reach the final state. A zoo of patterns, along with the existence of multiple states, was uncovered by Andereck et al. [5] in an incompressible TCF with independently rotating inner and outer cylinders. In the case of small aspect-ratio (Γ = h/δ ∼ O(1), where h and δ = R o − R i are the height of two cylinders and the annular gap between them, respectively) Taylor-Couette (TC) systems, the Taylor vortices that break the mid-plane axial symmetry are admissible as verified in experiments [6]. Anomalous Taylor vortices, with outflow near the top and bottom end-walls, have also been observed [7] experimentally.
The bulk of the scientific literature on TCF deals with incompressible fluids, but there have been limited works studying the effect of compressibility [8,9] and rarefaction [10] on this flow. Molecular Dynamics (MD) simulations have been successfully employed to study the TCF [11][12][13]. While the former work [11] demonstrated that the critical Reynolds/Taylor number for the onset of Taylor vortices and the nature of bifurcation can be quantitatively predicted by small-scale MD simulations of a Lennard-Jones fluid, the compressibility effects on the flow dynamics were not analysed. The compressibility ef- * e-mail: g_nandu@cb.amrita.edu * * e-mail: meheboob@jncasr.ac.in fects and the microscopic dissipation can significantly affect the flow patterns; for example, pronounced variations of density along both radial and axial directions were observed [13] away from the critical point. After a brief overview of our recent work [13] that pertains to a TCcell with small aspect ratio (Γ ≤ 4), the MD simulation results of a TC-cell with an aspect ratio Γ = 3 is discussed. This geometry should, intuitively, admit Taylor vortices that break the mid-plane axial symmetry close to the critical point. This is followed by discussions regarding a relatively tall TC-cell with Γ = 10 where simulations are carried out to further explore the bifurcation scenario and the multiplicity of states in a dense granular gas undergoing TCF.

Simulation details
The open source LAMMPS MD code [14], compiled on a mixed CPU-GPU architecture, with modifications to model the Taylor-Couette flow with the 'reflecting' endwall conditions [13] is employed for the simulations. The particles interact with one another with a Lennard-Jones(LJ) potential resulting in a force of the form F e (r) = 24 r 2 σ , where r is the distance between two particle centres, σ is the effective diameter of the particles and is the depth of the potential well -the cutoff distance, r c = 2 1/6 σ, results in a purely 'repulsive' potential known as the Weeks-Chandler-Andersen (WCA) potential. A damping force of the form, F d = −mγ n v n implements a velocity-dependent normal restitution force as a function of mass of the particle (m), the normal component of the relative velocity (v n ) and the damping coefficient (γ n ). The parameters chosen for the potential are σ = = m = 1 and the damping coefficient is set to γ n = 5 (thereby qualifying the system to be a dissipative/granular gas). critical-point isoline Figure 1: Phase diagram of various patterns in the (ω i , Γ)-plane; the velocity-vectors plots in the (r, z)-plane at selected values of (ω i , Γ) are shown on side panels. Adapted from Ref. [13].
For the results discussed herein, the gap width is δ = R o − R i = 75σ − 50σ = 25σ and the radius ratio is η = R i /R o = 2/3; the temperatures of the inner and outer walls (T w ) are set to 1. Reflecting boundary conditions are enforced at the top and bottom boundaries. The height h of the cylinder is varied such that Γ = h/δ ∈ (0.5, 10); for example, results will be presented for a geometry corresponding to Γ = h/δ = 3 (i.e. h = 75σ) and a relatively tall TC-cell with Γ = 10 (i.e. h = 250σ). The latter system corresponds to the number of particles being N = 1.25 × 10 6 for an averaged reduced density of ρ av = 0.5, placing the simulations in the domain of a 'dense' granular gas. In all simulations, the outer cylinder is stationary (i.e. ω o = 0) and only the inner-cylinder is rotated with angular velocity ω i . The particle positions are time integrated via a velocity-Verlet integrator with a time-step of ∆t = 0.005 [13].
Since the Taylor vortices represent azimuthally invariant states, the simulation data are averaged over azimuthal (θ) directions at any radial and axial positions. Spatial binning is performed with d r = d z = 1.25σ, along the radial (r) and axial (z) directions, respectively. The instantaneous reduced density at any time t is related to the number of particles in the bin N bin (t) and the volume of the bin, V bin , as ρ = σ 3 N bin (t)/V bin . The characteristic length and time scales are δ and ω −1 i resulting in dimensionless length z * = z/δ, time t * = tω i and velocity v * = v/ω i δ.
The simulation protocol consists of two stages: (i) an equilibration phase of 10 5 time steps during which the inner cylinder is uniformly accelerated from rest to the prescribed angular velocity, ω i and (ii) a production phase of 2 × 10 6 time steps during which the hydrodynamic fields are spatially averaged over the radial and axial bins and temporally averaged over every 2000 time-steps. For example, the primary bifurcation leads to a single vortex at Γ = 1, two, three and four vortices, stacked along the axial direction, at Γ = 2, 3 and 4, respectively. It is seen from Fig. 2 that the critical rotation rate (ω cr i ) for "CCF → T V F" transition is nearly constant ω cr ≈ 0.054 at Γ ≥ 1; however, ω cr i decreases sharply at Γ < 1, implying that the flow is more stable in a shorter TC cell having its height lesser than its annular gap (with h < δ).

Taylor vortices and the role of microscopic damping
Based on the speed of the inner-cylinder v i = ω i R i and the gap-width δ = R o − R i , the flow Reynolds number can be estimated from Re = v i δ/µ, where = ρ p φ = (πρ/6)ρ p is the mass-density of the fluid, ρ p is the material density of particles. The shear viscosity, µ, of a dense hard-sphere fluid is given by [15] (1) Note that g 0 (φ) = (1−φ/2)/(1−φ) 3 is the radial distribution function and f 2 (φ) is a dimensionless function of particle volume fraction (φ). Using (1), the Reynolds number can be expressed as with T * = T/(ω i δ) being the dimensionless temperature, and we have set k B = 1 = m as in MD simulations. Referring to Fig. 1, the mean/global temperature at the onset of Taylor vortices is T * = T * av ≈ 0.28 (i.e. at ω i = ω cr ≈ 0.054) at Γ = 4 [13]. Substituting the above value of T * , along with f 2 (φ av ) ≈ 1.1029 (evaluated at a mean volume fraction of φ = φ av = πρ av /6 ≈ 0.2618) and R i /σ = 50 in Eq. (2), we find a critical Reynolds number of Re cr ≈ 85.6 that holds for a granular gas (γ n = 5). For the same system without microscopic damping (γ n = 0), we found ω cr ≈ 0.07, and T * = T * av ≈ 0.403 [13], resulting in a critical Reynolds number of Re cr (γ n = 0) ≈ 71.4 < Re cr (γ n 0).
The above estimates do not depend on the aspect-ratio (Γ = h/δ) of the TC-cell as long as Γ > 1. Therefore, we conclude that the microscopic dissipation (γ n > 0) stabilizes this flow, thereby delaying the onset of Taylor vortices. It should be noted that the estimate of Re cr (γ n = 0) ≈ 71.4 from MD simulations is very close to the linearstability prediction of Re cr = 76.4 [16] for an incompressible TCF with a radius ratio of η = 2/3. Returning to the phase-diagram in Fig. 1, we note that the primary transition at Γ = 3 yields a 3-vortex system at ω i > ω cr ≈ 0.054. The evolution of vortices/rolls with time can be understood from Fig. 2; while the main panel shows the temporal evolution of the system temperature (defined as T * s = ( r z T (r, z)/n r n z )/(δω i ) 2 ) for ω i = 0.06 (green line) and ω i = 0.070 (red line), the insets and the side panels display the corresponding meridional velocity fields. For both cases, it is seen that the initial pattern is a 4-roll configuration which transitions, over time, with the coarsening of the top-most roll which is followed by the enlargement of the remaining 3 rolls to span the axial height (h). Note that, in case of ω i = 0.070, the 4-roll to 3-roll transition is accompanied by a distinct reduction in energy, a trait that is not discernible in the case of ω i = 0.060. For Γ = 3 where one would intuitively assume an asymmetric pattern (3 rolls) as the system transitions from CCF to TVF, it is seen that the transition first occurs to a symmetric pattern (4 rolls) which then transitions to an asymmetric (3 roll) state at large enough time when the system equilibrates to a statistical steady statethe 3-roll state seems to be a lower energy state at ω i ∼ ω cr and hence stable.

Secondary bifurcation and results for Γ = 10
Extensive results presented in Ref. [13] confirmed that the secondary transition from the TVF state, resulting in different number of rolls (at ω i > ω cr ; see Fig. 2), is in fact hysteretic. For example, at Γ = 4, the primary bifurcation results in a 4-roll state at ω i = 0.04 which persists over ω i ∈ (0.054, 0.3) beyond which the 6-roll state takes over; moreover, the 5-roll state was found to be stable over 0.18 ≤ ω i ≤ 0.37. Therefore, the 4-, 5-and 6-roll states can coexist with each other (viz. Fig. 11 in Ref. [13]) over a finite range of ω i at Γ = 4 in TCF of a dissipative granular gas; the underlying transition between 4 ↔ 5-rolls and 5 ↔ 6-rolls was found to occur via a saddle-node bifurcation. Interestingly, in the absence of microscopic dissipation (γ n = 0), the onset of 4-and 5-roll states was found to occur simultaneously at ω i ≈ 0.054; therefore, the 4-and 5-roll states coexist with each other over ω cr ≤ ω i ≤ 0.4 for a dense "molecular" gas undergoing TCF.
With parameter values as in Fig. 2, new simulations were conducted in a taller TC-cell with Γ = 10. The visualizations of the rolls in terms of the specific angular momentum, L = ρv θ r , are displayed in Fig. 3 at different values of ω i . With increasing ω i , we find "11 → 10 → 11 → 12 → 12 * → 12 → 13 → 10 → 12 * → 14 * " rolls; with some states consisting of 'outward' jets near both end-walls; note that the even-roll states with inward and outward end-wall jets are equivalent solutions belonging to the upper and lower branches of perfect pitchfork bifurcation for reflective end-wall conditions. The outward jets, transporting the angular momentum from the rotating inner cylinder to the stationary outer cylinder are clearly visible.
The above bifurcation sequence among differentnumbered roll-states appears to be smooth when the maximum radial velocity,  is used as an order parameter. The solid line in figure 4(a) is a fit through the data points upto ω i ≤ 0.25: with ω cr ≈ 0.054 and β = 1/2. The critical rotation rate for the onset of Taylor vortices is approximately the same as that for Γ = 1, 2, 3 and 4 as marked by the neutral curve in figure 1. However, the power-law exponent in Eq. (6) is β ≈ 1/2 which can be compared with its value of β(Γ = 2) ≈ 0.585 and β(Γ = 4) ≈ 0.675 for smaller aspect-ratio systems [13]. This indicates that the wellknown square-root scaling, ∆v r = α(ω i /ω cr − 1) 1/2 , for pitchfork bifurcations would hold in MD simulations of TCF if the system-size is sufficiently large. The nature of secondary bifurcations from even to oddnumbered vortices is difficult to ascertain from figure 4(a) since the solution branches for different roll-numbers are closely located in the (∆v r , ω i )-plane. However, it is clear from Fig. 4(b) that the odd-and even-numbered Taylor vortices can co-exist with each other; note that the 11-(squares) and 13-roll (triangle) states appear to lie on different branches in the (ζ, ω i )-plane. Note that these results were also obtained from individual runs by slowly increasing the rotation rate of the inner-cylinder to the desired ω i and allowing the system to attain a statistical steady state over a long time. As discussed in Ref. [13], a numerical continuation scheme should be used to track a particular roll-state and the associated limit/turning points (if any) with increasing/decreasing ω i for Γ = 10 too. This requires time-consuming simulations, which will be addressed in the future. The present miscroscopic model is a toy model of a granular gas. A more realistic microscopic model [12,17] must incorporate normal and tangential restitutions as well as Coulomb friction; this along with athermal walls should be implemented to analyse "granular" TCF (without gravity) so as to make concrete conclusions on the role of microscopic dissipations on Taylor vortices and subsequent secondary transitions. However what is interesting is that this simplistic toy model of a granular gas correctly predicts the transition from CCF and is sufficient to capture hysteretic phenomena such as the multiplicity of states.