2.1. Boundary conditions

As in the work of Hirshfeld & Rapaport (Reference Hirshfeld and Rapaport1998), the inner and outer cylinders are treated as thermal walls with a wall temperature $T_w$. Any particle that collides with the inner or outer cylindrical wall has all memory of its velocity erased and is assigned a random velocity in addition to a mean velocity corresponding to the tangential velocity of the cylinder ($v_\theta =\omega _iR_i$). Referring to figure 1$(b)$ and using the transformation rule between Cartesian ($x,y,z$) and polar ($r, \theta , z$) coordinates, a reflected particle from the cylindrical walls is assigned the following velocity:

(2.2)\begin{equation} (c'_x, c'_y, c'_z )= \sqrt{\frac{k_B T_w}{m}} \,(v_{Gx}, v_{Gy}, v_{Gz} ) + \omega_iR_i(-\sin\theta, \cos\theta, 0), \end{equation}

where $v_{Gi}$ is sampled from the Gaussian distribution

(2.3)\begin{equation} p(v_{Gi}) = \left(\frac{m}{2{\rm \pi} k_BT_w}\right)^{1/2}\exp\left(- \frac{mv_{Gi}^2}{2k_B T_w}\right), \end{equation}

with $k_B$ being the Boltzmann constant. If the assigned velocity is such that the particle leaves the simulation domain, then the directional sense of $v_{Gx}$ and $v_{Gy}$ is inverted ($v_{Gx}=-v_{Gx}$ and $v_{Gy}=-v_{Gy}$).

The axial boundary condition is implemented using the following rule: any particle that moves outside the top or bottom endwalls, by a distance $s_z$ during time integration, is put back into the box by the same distance and the $z$-component of the velocity of that particle is inverted such that

(2.4)\begin{equation} \boldsymbol{c}' \boldsymbol{\cdot} \hat{\boldsymbol{k}} = - \boldsymbol{c} \boldsymbol{\cdot} \hat{\boldsymbol{k}}, \end{equation}

where $\hat {\boldsymbol {k}}$ is the unit vector along the axial direction; note that the reflecting particle's $x$- and $y$-component velocities remain unchanged.

The above thermal-wall condition (2.2) and reflective endwall condition (2.4) are similar to those employed in the previous MD simulations (Hirshfeld & Rapaport Reference Hirshfeld and Rapaport1998; Trevelyan & Zaki Reference Trevelyan and Zaki2016) of the TCF; similar boundary conditions were also used by Stefanov & Cercignani (Reference Stefanov and Cercignani1993) to study the onset of vortices in Monte Carlo simulations of the TCF of a rarefied gas. The role of stationary endwalls, resulting in no-slip velocity, will be briefly analysed for few specific cases in § 5.

2.2. Simulation methodology with CPU and GPU

The open-source LAMMPS MD code (Plimpton Reference Plimpton1995) is employed to integrate the equations of motion in the Cartesian coordinate system using a velocity–Verlet integrator. To speed up computations, LAMMPS has been built from source with the graphics processing unit (GPU) package enabled and optimized for performance on a GPU with 3072 CUDA (Compute Unified Device Architecture) cores and 12 GB GPU random-access memory (RAM). The central processing unit (CPU) computational component consists of a system total of 64 GB CPU RAM equally distributed amongst two sockets, each housing a 2.4 GHz Haswell architecture processor having eight cores, of which 14 ($2\times 7$) cores are employed for computations. The average scale-up obtained by this ‘CPU + GPU’ combination for the family of runs corresponding to $\varGamma =2$ and $\delta =25 \sigma$ is found to be approximately $2.3$ times that of the CPU-only runs employing 14 cores. The scale-up performance is seen to decay with both increasing $\omega _i$ and increasing number of particles. The latter can be attributed to the more frequent neighbour list rebuilds, which can occur in the CPU only, owing to the nature of the potential used.

All simulations have been carried out for a mean reduced density of $\rho _{av}=0.5$ as in previous studies on TCF (Hirshfeld & Rapaport Reference Hirshfeld and Rapaport1998; Trevelyan & Zaki Reference Trevelyan and Zaki2016); the aspect ratio ($\varGamma = h/\delta = 0.5, 0.75, 1.0, 1.33, 2, 4, 10$) was varied by changing the height ($h$) of the cylinders. In most simulations, the outer and inner cylinder radii are taken as $R_o=75\sigma$ and $R_i=50\sigma$, such that the gap width is $\delta =(R_o-R_i)=25\sigma$, the radius ratio is $\eta =R_i/R_o=2/3$ and the curvature is $\kappa =\delta /R_i=1/2$. The latter two quantities ($\eta , \kappa$) suggest that our simulations correspond to ‘wide-gap’ TC geometry. The total number of particles within the annular gap is calculated from

(2.5)\begin{equation} N = \rho_{av}{\rm \pi} \left(\frac{R_o^2}{\sigma^2}\right) (1-\eta^2 ) \varGamma\left(\frac{\delta}{\sigma}\right), \end{equation}

which varied from $60\,000$ to $1.25\times 10^6$, depending on the value of $\varGamma \in (0.5,10)$. At smaller values of $\varGamma < 1$, a four-fold increase of $N$ (by keeping $\eta$, $\kappa$ and $\varGamma$ fixed) did not affect the observed patterns.

The simulation procedure consists of two stages: (i) an equilibration run that consist of $10^{5}$ time steps during which the inner cylinder is uniformly accelerated from rest to the required angular velocity and (ii) a production run consisting of another $2\times 10^6$ time steps during which the (instantaneous) hydrodynamic fields are calculated in each bin at every step, and then time-averaging over these bin-averaged properties is carried out over every $2000$ time steps as detailed in § 2.3. Most simulations have been carried out with a time step of ${\rm \Delta} t=0.005$; a few simulations were checked with smaller time steps of ${\rm \Delta} t=0.0025$ (at $\varGamma =2$ and $4$) and $0.001$ (at $\varGamma =4/3$), which did not alter the results.

2.3. Hydrodynamic fields and the statistical steady state of TVF

To calculate hydrodynamic fields from the instantaneous particle positions and velocities, the simulation domain is binned along the radial and axial directions with bin sizes of $\textrm {d}z= \textrm {d}r =1.25\sigma$. Since we will be presenting results on CCF and TVF, both of which being azimuthally invariant states, the hydrodynamic fields will be calculated by averaging particle data over all $\theta$. Note that the mapping from the Cartesian $(x,y,z)$ to the polar $(r, \theta , z)$ coordinate system is done via the planar rotation matrix to obtain the velocity components of a particle ($c_r, c_\theta , c_z$) in polar coordinates.

The instantaneous ‘reduced’ density is defined as $\rho ={\sigma ^3 N(t)}/{V_{bin}}$, where $V_{bin}$ is the volume of the bin and $N(t)$ is the instantaneous number of particles in a bin; the reduced density is related to the particle volume fraction ($\phi$) via $\phi =({\rm \pi} /6)\rho$. For a bin in the ($r, z$) plane, the local reduced density is calculated from

(2.6)\begin{equation} \rho(r,z)=\frac{{\sigma^3}}{{\rm \pi}(R_{bo}^2 - R_{bi}^2)\,\textrm{d}z} \left\langle\sum_{i=1}^{N}{\delta(r-r_i)\delta(z-z_i)}\right\rangle_\theta , \end{equation}

where $R_{bo}$ and $R_{bi}$ denote the outer and inner radii of the bin whose axial height is $\textrm {d}z$, and ($r_i, z_i$) denote the coordinates of all particles residing in this particular bin. The subscript $\theta$ on the angular brackets in (2.6) implies a spatial averaging over all particles in the $\theta$-direction, in addition to temporal averaging over a few thousand snapshots of the system in the steady state.

The instantaneous velocity of a group of particles within a bin, at any instant in time, is defined as ${\sum _{i=1}^{N}\boldsymbol {c}}/{N(t)}$, where $\boldsymbol {c}$ denotes the instantaneous velocity vector of the particles in that bin. Thus, the local hydrodynamic velocity is calculated from

(2.7)\begin{equation} \boldsymbol{v}(r,z)=\frac{\left\langle {\sum_{i=1}^{N} \boldsymbol{c_{i}}(t) \delta(r-r_i)\delta(z-z_i)} \right\rangle_\theta}{\left\langle\sum_{i=1}^{N}{\delta(r-r_i) \delta(z-z_i)}\right\rangle_\theta} , \end{equation}

where $\boldsymbol {c}_i$ is the instantaneous velocity of the $i$th particle in the bin, and the angular brackets denote an averaging over time as well as $\theta$.

The local kinetic temperature in a bin, which is a measure of the fluctuation velocities of the particles within that bin, is defined as

(2.8)\begin{equation} T(r,z)=\frac{\left\langle m \sum_{i=1}^{N} \left({\boldsymbol{c}_i(t)-\boldsymbol{v}(r,z)}\right) \boldsymbol{\cdot} ({\boldsymbol{c}_i(t)-\boldsymbol{v}(r,z)})\delta(r-r_i)\delta(z-z_i) \right\rangle_\theta}{3k_B \left\langle\sum_{i=1}^{N}{\delta(r-r_i)\delta(z-z_i)}\right\rangle_\theta} . \end{equation}

For the simulations considered here, $k_B=m=1$ and ${\boldsymbol C}=(\boldsymbol {c}-\boldsymbol {v})$ denotes the peculiar velocity of the particles within the bin during the time frame of interest. Using the gap width ($\delta =R_o-R_i$) and the inverse of the inner-cylinder rotation rate $(\omega _i^{-1})$ as the length and time scales, respectively, the following dimensionless fields are defined:

(2.9a–d)\begin{equation} z^* = \frac{z}{\delta}, \quad t^* = {\omega_i}t, \quad {\boldsymbol v}^* = \frac{\boldsymbol v}{ \omega_i\delta }, \quad T^* = \frac{T}{(\omega_i\delta)^2}. \end{equation}

To ascertain whether the system reached a statistical steady state, we monitor the global temperature defined as

(2.10)\begin{equation} T_s=\frac{\sum^{r}\sum^{z} T(r,z)}{n_r n_z }, \end{equation}

where $n_r$ and $n_z$ denote the number of bins along the radial and axial directions, respectively. The constancy of the global temperature $T_s$ with time is taken as the statistical steady state of the flow – see the main panels in figure 2. The inset of each panel in figure 2 displays a pair of counter-rotating rolls at the steady state for parameter values of $(a)$$\gamma _n=5$ and $(b)$$\gamma _n=0$ (i.e. with zero microscopic damping force), with $\omega _i=0.1$ and $\varGamma =h/\delta =2$. For both cases, there is an ‘outward’ jet at the mid-height of the cylinder and two ‘inward’ jets near the two endwalls – the resulting pattern is called TVF or Taylor rolls (Taylor Reference Taylor1923). All results presented below correspond to a damping coefficient $\gamma _n=5$ (Rapaport Reference Rapaport1998, Reference Rapaport2004) but the robustness of our reported results is verified in appendix A by performing simulations with $\gamma _n=0$.

Figure 2.

Main panel: temporal evolution of the global kinetic temperature (dimensionless), equation (2.10), of the fluid. Inset: snapshot of steady Taylor vortices, displaying radial–axial velocity vectors. Parameter values are $(a)$$\gamma _n=5$ (microscopic damping coefficient) and $(b)$$\gamma _n=0$, with $\varGamma =2$ and $\omega _i=0.1$.

3.1. Determination of the critical rotation rate and Reynolds number

To identify the location of the primary bifurcation point $\omega _i=\omega ^{cr}$ for the onset of TVF, the simulations have been performed for a range of $\omega _i$ for specified values of the aspect ratio $\varGamma$. Based on the radial velocity variations such as in figure 3$(a)$, the maximum radial velocity,

(3.1)\begin{equation} {\rm \Delta} v_r =\max _{0\leq z^* \leq \Gamma} [v_r (z)], \end{equation}

is chosen as an order parameter to identify the transition from CCF to TVF. The variation of ${\rm \Delta} v_r$ with $\omega _i$ is displayed in figure 3$(b)$. The critical rotation rate $\omega _i^{cr}$ for transition from CCF to TVF is obtained by fitting all data for ${\rm \Delta} v_r$ using the following function:

(3.2)\begin{equation} {\rm \Delta} v_r = \alpha(\omega_i - \omega^{cr})^{\beta} , \end{equation}

with $\alpha \approx 2.05$ and $\beta \approx 0.565$. Estimating the bifurcation point (figure 3$b$) at which ${\rm \Delta} v_r\approx 0$ yields a value of $\omega ^{cr}\approx 0.054$ for the present system with $\varGamma =2$ and $\omega _o=0$; this critical rotation rate holds for all $\varGamma >1$ (see § 3.2).

Based on the inner-cylinder speed $v_{i}=\omega _{i} R_i$ as the velocity scale and the gap width $\delta =R_o-R_i$ as the length scale, the Reynolds number can be estimated from

(3.3)\begin{equation} Re = \frac{\varrho v_{i} \delta}{\mu} , \end{equation}

where $\varrho =\rho _p\phi =({\rm \pi} \rho /6)\rho _p$ is the mass density of the fluid, $\rho _p$ is the material density of particles and $\mu$ is the shear viscosity. The expression of the latter for a ‘dense’ hard-sphere fluid is given by Chapman & Cowling (Reference Chapman and Cowling1970) and Saha & Alam (Reference Saha and Alam2016) as

(3.4)\begin{equation} \frac{\mu}{\varrho \sigma\sqrt{(k_B/m)T}} = \frac{5\sqrt{\rm \pi}}{96\phi g_0}\left(1 + \frac{8}{5} \phi g_0\right)^2 + \frac{8}{5\sqrt{\rm \pi}}\phi g_0 =f_2(\phi), \end{equation}

where $g_0(\phi ) =(1-\phi /2)/(1-\phi )^3$ is the radial distribution function and $f_2(\phi )$ is a dimensionless function of particle volume fraction ($\phi$). Combining (3.3) and (3.4) we arrive at the following expression for the Reynolds number:

(3.5)\begin{equation} Re = \frac{\varrho \omega_{i}R_i \delta}{\varrho \sigma \sqrt{(k_B/m)T} \, f_2(\phi) } \equiv \frac{(R_i/\sigma)}{f_2(\phi)\sqrt{T^*}}, \end{equation}

with $T^*=T/(\omega _i\delta )$ being the dimensionless temperature, and we have set $k_B=1=m$ as in the present simulations. The mean/global temperature at the onset of Taylor vortices is $T^*=T_{av}^*\approx 0.28$ (i.e. at $\omega _i=\omega ^{cr}\approx 0.054$). Substituting the value of $T^*$, along with $f_2(\phi _{av})\approx 1.1029$ (evaluated at a mean volume fraction of $\phi =\phi _{av}={\rm \pi} \rho _{av}/6\approx 0.2618$) and $R_i/\sigma =50$ in (3.5), we find a critical Reynolds number of $Re^{cr}\approx 85.6$.

On the other hand, for the same system without microscopic damping ($\gamma _n=0$, see appendix A), we have $\omega ^{cr}=0.07$ and $T^*=T_{av}^*\approx 0.4028$, resulting in a critical Reynolds number of $Re^{cr}(\gamma _n=0)\approx 71.4 < Re^{cr}(\gamma _n\neq 0)$. The above estimates do not depend on the aspect ratio ($\varGamma =h/\delta$) of the TC cell as long as $\varGamma > 1$. Overall, the microscopic dissipation ($\gamma _n>0$) increases the value of the critical Reynolds number, thereby delaying the onset of Taylor vortices, which is expected. Note that our estimate of $Re^{cr}(\gamma _n=0)\approx 71.4$ is close to the predicted value of $Re^{cr}=76.4$ (Sparrow, Munro & Jonsson Reference Sparrow, Munro and Jonsson1964) for an incompressible TCF with a radius ratio of $\eta =2/3$.

3.2. Phase diagram of patterns: symmetric and asymmetric vortices

Following the procedure described in § 3.1, the critical rotation rate $\omega ^{cr}$ for the bifurcation from CCF to TVF has been determined for different values of the aspect ratio $\varGamma$. The neutral/critical locus of $\omega _i=\omega ^{cr}$ is marked in figure 4 by the dashed and dot-dashed lines. It is seen that the critical rotation rate is nearly constant ($\omega ^{cr}\approx 0.054$) for $\varGamma \geq 1$, and the data points for $\varGamma \in (0.2, 1)$ can be fitted via

(3.6)\begin{equation} \omega^{cr} = \frac{0.051}{\varGamma}, \quad \forall \ \varGamma <1, \end{equation}

marked by the dashed line in figure 4. Equation (3.6) implies that the critical rotation rate increases for short cylinders or, in other words, the transition from CCF to TVF is delayed in short cylinders with $\varGamma < 1$.

Figure 4.

Phase diagram of patterns in the ($\varGamma , \omega _i$) plane. Different symbols denote different numbers of Taylor vortices as shown in the legend (bottom right corner); the empty symbols denote a sense of rotation opposite to the one observed in the case of the corresponding filled symbols. Each symbol at ($\varGamma , \omega _i$) corresponds to a steady-state pattern that has been obtained from an individual run (such as figure 2) by uniformly accelerating the inner cylinder to a specified $\omega _i$ in $10^5$ time steps (${\rm \Delta} t$) from the rest state, and subsequently allowing the simulations to run for $2\times 10^6{\rm \Delta} t$ to calculate hydrodynamic fields; see § 2.2 for details on simulation protocols. The vector visualizations corresponding to points A, B, D and F have been scaled by factors of $15$, $2$, $5$ and $2$, respectively, to improve clarity.

The main panel of figure 4 summarizes all patterns for selected aspect ratios $\varGamma \in (0.2,4)$ that have been obtained from long-time simulations (such as in figure 2) by varying the inner-cylinder rotation rate $\omega _i \in (0.05,0.4)$ that approximately spans a range of inner-cylinder Reynolds number of $Re=\omega _iR_i\delta /(\mu /\varrho )\in (50,\ 500)$. The Taylor vortices having different numbers of rolls are marked in figure 4 by different symbols at each $\varGamma$. The steady-state snapshots with increasing $\omega _i>\omega ^{cr}$, depicting transitions from one pattern to another, are shown on the left and right panels of figure 4 in terms of the radial–axial velocity-vector plots. For example, the snapshots I and J represent a transition from a two-roll to a three-roll TVF (marked by diamonds and triangles, respectively, in figure 4) that occurs at $\omega _i\approx 0.27$ at $\varGamma =2$. A transition from a four-roll to a five-roll state can be ascertained from the snapshots K and L displayed on the lower right panels that occurs at $\varGamma =4$.

The snapshots A $\to$ B and C $\to$ D represent transition from (i) an ‘unsteady’ state to a ‘steady’ single-roll state and (ii) a single-roll state to an ‘asymmetric’ two-roll state at $\varGamma =0.5$ and $0.75$, respectively. When the aspect ratio is increased to $\varGamma =1$, a transition from a single-roll state to a ‘symmetric’ two-roll state is found to occur; see the snapshots E $\to$ F. On the other hand, the snapshots G $\to$ H represent a transition between two possible two-roll configurations, with one having ‘outward’ flows (G) near the endwalls and the other having ‘inward’ flows (H) near the endwalls. The latter transition will be analysed in detail in § 5.1.

Among the patterns identified in figure 4, the odd-numbered or asymmetric rolls break the $Z_2$-symmetry (reflection about the mid-height $z=h/2$). Such asymmetric states, along with the midplane symmetric Taylor vortices with outward flow at the endwalls, are collectively known as ‘anomalous modes’ according to the classification of Benjamin (Reference Benjamin1978a,Reference Benjaminb). The latter-type vortices have been extensively studied in incompressible TCF with no-slip, free-slip and asymmetric boundary conditions at the endwalls (Benjamin & Mullin Reference Benjamin and Mullin1981; Cliffe Reference Cliffe1983; Lücke et al. Reference Lücke, Mihelcic, Wingerath and Pfister1984; Nakamura et al. Reference Nakamura, Toya, Yamashita and Ueki1989; Mullin & Blohm Reference Mullin and Blohm2001). In the context of the present findings, pertinent questions are these: Are the present asymmetric roll states similar to those known in incompressible fluids? How does the compressibility of the fluid affect the primary and secondary bifurcation scenario? Does the bifurcation between different types of roll states occur via supercritical or subcritical transitions? These questions are addressed below.

4.2. Single-vortex and asymmetric two-vortex states: $\varGamma \leq 1$

Here we discuss results pertaining to short cylinders with aspect ratios of $\varGamma \leq O(1)$; some of the observed patterns are marked in the phase diagram in figure 4. The bifurcation diagram for $\varGamma =1$ is shown in the middle panel of figure 12 in the (${\rm \Delta} v_r, \omega _i$) plane. The primary bifurcation from the CCF state yields a ‘single-vortex’ state for which three snapshots are displayed on the left ($\omega _i=0.065$) and top ($\omega _i=0.15$ and $0.3$) panels of figure 12. We use the terminology ‘single vortex’ to refer to the existence of one vortex spanning the axial height of the TC cell. Each snapshot represents streamline patterns (at the steady state) in the ($r, z$) plane. This has been drawn by employing the ‘surface-LIC’ (surface line integral convolution) module of the ParaView software (Ahrens, Geveci & Law Reference Ahrens, Geveci and Law2005). Note that the single-vortex state is slightly off-centred, indicating its geometrical asymmetry, which persists over the range of $\omega _i$ studied. The bottom row and the right column of figure 12 display streamline patterns that belong to the lower branch of the bifurcation diagram. An asymmetric pair of vortices is found at $\omega _i=0.085$ (the bottom left panel) that persists up to a rotation rate of $\omega _i=0.1$. Both these asymmetric two-roll states consist of a stronger vortex lying below a weaker one. With increasing rotation rate, the asymmetric two-vortex state degenerates into a symmetric two-roll state; see the snapshots for $\omega _i=0.25$ (the bottom right panel), $\omega _i=0.3$ (middle right panel) and $\omega _i=0.4$ (top right panel) in figure 12.

Figure 12.

Evolution of Taylor rolls with increasing $\omega _i$ along the upper and lower branches of the bifurcation diagram (middle panel) for $\varGamma =1$. The middle panel represents the bifurcation diagram in the (${\rm \Delta} v_r, \omega _i$) plane, with the circles and diamonds representing single-roll and two-roll states, respectively, obtained from simulations in which the inner cylinder is started from the rest state to a specified $\omega _i$ as described in § 2.2. The triangles on the upper and lower branches were obtained from two separate ‘continuation’ runs in which the inner-cylinder rotation $\omega _i$ is (i) increased from $\omega _i=0.25$ (upper branch) to $0.35$ and (ii) decreased from $\omega _i=0.17$ (lower branch) to $0.06$ in steps of ${\rm \Delta} \omega _i=0.01$. The ramping rate is $|\textrm {d}\omega _i/\textrm {d}t|= 2\times 10^{-5}$ for both cases. Refer to figure 21 for the protocol on ‘$\omega _i$ continuation’ runs.

The above sequence of pattern evolution can be understood from the bifurcation diagram (viz. the middle panel in figure 12) on which the upper and lower branch solutions are connected via two limit points through an unstable branch (marked by the dashed line). It is clear that the underlying bifurcation is mediated via a hysteretic transition, resulting in (i) ‘1-roll $\rightarrow$ 2-rolls’ transition with increasing $\omega _i$ from the CCF, with the transition to the two-roll state occurring at the right limit point $\omega _i\approx 0.3$, and (ii) ‘symmetric 2-rolls $\to$ asymmetric 2-rolls $\to$ 1-roll’ transition with decreasing $\omega _i$ from the symmetric two-roll state, with the transition to the single-roll state occurring at the left limit point $\omega _i\approx 0.08$. We have verified that the above bifurcation sequence also holds at smaller values of $\varGamma =0.75$ for which the bifurcation diagram looks similar to that in figure 12. These findings on the overall bifurcation scenario for the single-vortex mode and the impact of stationary endwalls on them will be compared and contrasted with related works on incompressible TCF (Benjamin & Mullin Reference Benjamin and Mullin1981; Furukawa et al. Reference Furukawa, Watanabe, Toya and Nakamura2002) in § 5.2.

The hydrodynamic fields for $\varGamma =1$ are displayed in figure 13. Note that the radial profiles of the azimuthal velocity $v_\theta (r)$ look structurally similar (with the formation of a plateau around the mid-gap at $\omega _i\geq 0.15$) to those for $\varGamma =2$ (figure 6$c$) and $\varGamma =4$ (figure 9$b$). The asymmetric $v_\theta (z)$ profiles for $\omega _i=0.065$ and $0.15$ in figure 13(d) represent the single-vortex mode (viz. figure 12). The remaining profiles for $\omega _i=0.1$, $0.3$ and $0.4$ represent two-roll modes and confirm the existence of a mid-height ‘inward’ jet and two ‘outward’ jets near two walls in each case. The key signature of compressible Taylor vortices is strong density stratification along both radial and axial directions that persists also in short cylinders with $\varGamma \leq O(1)$ as is evident from figure 13$(b{,}e)$. The specific angular momentum field ${\mathcal {L}}(r)$ in figure 13$(c)$ indicates the emergence of its non-monotonic radial variation with increasing $\omega _i\geq 0.25$, similar to those uncovered for $\varGamma =2$ (figure 7$b$) and $\varGamma =4$ (figure 9$c$).

Figure 13.

Radial (top row) and axial (bottom row) variations of $(a{,}d)$ azimuthal velocity, $(b{,}e)$ reduced density and $(c{,}f)$ specific angular momentum for $\varGamma =1$. The corresponding snapshots can be found in figure 12, which refer to single-roll ($\omega _i=0.065, 0.15$), asymmetric two-roll ($\omega _i=0.1$) and symmetric two-roll ($\omega _i=0.25, 0.4$) states.

About the last issue, the origin of the non-monotonic angular momentum profiles (viz. figures 7$b$, 9$c$ and 13$c$) can be tied to the underlying density stratification, as we explain below. Recall that we have found strong radial variation of density (see figures 6$a$, 9$a$ and 13$b$), with the density maximum shifting towards the stationary outer cylinder with increasing $\omega _i$. As explained in § 4.1, the constancy of the shear stress across the Couette gap immediately leads to the formation of a plateau in the azimuthal velocity (figures 6$c$, 9$c$ and 13$a$) around the mid-gap with increasing $\omega _i$. Therefore, the coupling of the density and azimuthal velocity via the shear stress is likely to be responsible for the non-monotonic radial variation of the angular momentum (${\mathcal {L}}(r)=\langle \rho v_\theta r \rangle$) in the compressible TCF of a dense gas. Collectively, the signature of the fluid compressibility is implicated as stratifications in all hydrodynamic fields, and the structural characteristics of the hydrodynamic profiles with increasing $\omega _i$ look similar at all $\varGamma$ (irrespective of the number of vortices along the axial direction).