2.1. Navier–Stokes equations under the Boussinesq approximation

In this study, we consider the Boussinesq approximation in which the reference density $\rho _{0}$ is assumed to be much larger than the density variation $\varrho$ (i.e. $\rho _{0}\gg \varrho$ ). The density variation $\varrho$ is assumed to satisfy a linear relation with the total temperature $\vartheta$ as $\varrho /\rho _{0}=-\alpha _{0} (\vartheta -\vartheta _{0})$ where $\alpha _{0}\gt 0$ is the thermal expansion coefficient and $\vartheta _{0}$ is the reference temperature. For velocity ${\boldsymbol{U}}=(U_{r},U_{\theta },U_{z})$ , temperature variation $\varTheta =\vartheta -\vartheta _{0}$ and associated pressure variation $P$ in cylindrical coordinates $(r,\theta ,z)$ , we consider the following continuity, momentum and energy equations:

(2.1) \begin{align} \nabla \cdot {\boldsymbol{U}}=0, \\[-24pt] \nonumber \end{align}

(2.2) \begin{align} \frac {\partial {\boldsymbol{U}}}{\partial t}+{\boldsymbol{U}}\cdot \nabla {\boldsymbol{U}}=-\frac {1}{\rho _{0}}\nabla P +\alpha _{0}g\varTheta \vec{\boldsymbol{e}}_{z}+\nu _{0}\nabla ^{2}{\boldsymbol{U}}, \\[-24pt] \nonumber \end{align}

(2.3) \begin{align} \frac {\partial \varTheta }{\partial t}+{\boldsymbol{U}}\cdot \nabla \varTheta =\kappa _{0}\nabla ^{2}\varTheta , \\[-2pt] \nonumber \end{align}

where $g$ is the gravitational acceleration in the vertical direction $z$ , $\nu _{0}$ is the kinematic viscosity, $\kappa _{0}$ is the thermal diffusivity and $\nabla ^{2}$ is the Laplacian operator. We consider cylindrical Couette flow in a stably stratified fluid as a base state with base velocity $\boldsymbol{U}_{B}= (0,V_{B}(r),0 )$ and base temperature $T_{B}(z)$ as

(2.4) \begin{align} V_{B}(r)=Ar+\frac {B}{r},\quad A=\varOmega _{i}\frac {\mu -\eta ^{2}}{1-\eta ^{2}},\quad B=\varOmega _{i}R_{i}^{2}\frac {1-\mu }{1-\eta ^{2}},\quad T_{B}(z)=\frac {\Delta T}{\Delta z}z, \end{align}

where $V_{B}(r)$ is the base azimuthal velocity, $A$ and $B$ are constants as a function of the angular velocities $\varOmega _{i}$ and $\varOmega _{o}$ and radii $R_{i}$ and $R_{o}$ (where the subscripts $i$ and $o$ denote the inner and outer cylinders, respectively), $\mu =\varOmega _{o}/\varOmega _{i}$ is the angular velocity ratio, $\eta =R_{i}/R_{o}$ is the radius ratio and $\Delta T$ is the temperature difference along the vertical distance $\Delta z$ . Throughout this paper, we consider only the case in which the outer cylinder is fixed and only the inner cylinder rotates (i.e. $\varOmega _{o}=0$ and $\mu =0$ ). The temperature gradient $\Delta T/\Delta z$ is assumed to be positive and constant so the fluid is stably stratified and the base temperature $T_{B}$ increases linearly with $z$ . The base pressure $P_{B}(r,z)$ balances the velocity $V_{B}$ and temperature $T_{B}$ by satisfying the following relations:

(2.5) \begin{align} -\frac {V_{B}^{2}}{r}=-\frac {1}{\rho _{0}}\frac {\partial P_{B}}{\partial r},\quad 0=-\frac {1}{\rho _{0}}\frac {\partial P_{B}}{\partial z}+\alpha _{0}g T_{B}. \end{align}

Equations (2.1)–(2.3) can be expressed in a non-dimensional form by considering $R_{i}$ as the reference length, $\varOmega _{i}^{-1}$ as the reference time, $R_{i}\varOmega _{i}$ as the reference velocity, $\rho _{0}R_{i}^{2}\varOmega _{i}^{2}$ as the reference pressure and $\Delta T$ as the reference temperature. The distance $\Delta z$ can be chosen arbitrarily without loss of generality; thus we choose $\Delta z=R_{i}$ for simplicity. The coordinates $z$ and $r$ are considered non-dimensional hereafter, and hence the base temperature can simply be expressed as $T_{B}(z)=z$ .

We consider perturbation with its velocity $\boldsymbol{u}={\boldsymbol{U}}-\boldsymbol{U}_{B}= (u_{r},u_{\theta },u_{z} )$ , pressure $p=P-P_{B}$ and temperature $T=\varTheta -T_{B}$ . By applying the base state and perturbation to (2.1)–(2.3) and subtracting the base-state equations (2.4)–(2.5), we obtain the Navier–Stokes equations for perturbation as follows:

(2.6) \begin{align} \frac {\partial u_{r}}{\partial r}+\frac {u_{r}}{r}+\frac {1}{r}\frac {\partial u_{\theta }}{\partial \theta }+\frac {\partial u_{z}}{\partial z}=0, \\[-26pt] \nonumber \end{align}

(2.7) \begin{align} \frac {\partial u_{r}}{\partial t}+\varOmega \frac {\partial u_{r}}{\partial \theta }-2\varOmega u_{\theta }+{N}_{r}=-\frac {\partial p}{\partial r}+\frac {1}{Re}\left (\nabla ^{2}u_{r}-\frac {u_{r}}{r^{2}}-\frac {2}{r^{2}}\frac {\partial u_{\theta }}{\partial \theta }\right ), \\[-26pt] \nonumber \end{align}

(2.8) \begin{align} \frac {\partial u_{\theta }}{\partial t}+\varOmega \frac {\partial u_{\theta }}{\partial \theta }+Zu_{r}+{N}_{\theta }=-\frac {1}{r}\frac {\partial p}{\partial \theta }+\frac {1}{Re}\left (\nabla ^{2}u_{\theta }-\frac {u_{\theta }}{r^{2}}+\frac {2}{r^{2}}\frac {\partial u_{r}}{\partial \theta }\right ), \\[-26pt] \nonumber \end{align}

(2.9) \begin{align} \frac {\partial u_{z}}{\partial t}+\varOmega \frac {\partial u_{z}}{\partial \theta }+{N}_{z}=-\frac {\partial p}{\partial z}+N^{2}T+\frac {1}{Re}\nabla ^{2}u_{z}, \\[-26pt] \nonumber \end{align}

(2.10) \begin{align} \frac {\partial T}{\partial t}+\varOmega \frac {\partial T}{\partial \theta }+u_{z}+{N}_{T}=\frac {1}{RePr}\nabla ^{2}T, \\[1pt] \nonumber \end{align}

where $\varOmega (r)=V_{B}/r$ is the base angular velocity, $Z(r)=(1/r)\mathrm {d}(r^{2}\varOmega )/\mathrm {d}r$ is the base axial vorticity and $N_{r}$ , $N_{\theta }$ , $N_{z}$ and $N_{T}$ are nonlinear terms:

(2.11) \begin{align} {N}_{r}=u_{r}\frac {\partial u_{r}}{\partial r}+\frac {u_{\theta }}{r}\frac {\partial u_{r}}{\partial \theta }+u_{z}\frac {\partial u_{r}}{\partial z}-\frac {u_{\theta }^{2}}{r}, \\[-26pt] \nonumber \end{align}

(2.12) \begin{align} {N}_{\theta }=u_{r}\frac {\partial u_{\theta }}{\partial r}+\frac {u_{\theta }}{r}\frac {\partial u_{\theta }}{\partial \theta }+u_{z}\frac {\partial u_{\theta }}{\partial z}+\frac {u_{r}u_{\theta }}{r}, \\[-26pt] \nonumber \end{align}

(2.13) \begin{align} {N}_{z}=u_{r}\frac {\partial u_{z}}{\partial r}+\frac {u_{\theta }}{r}\frac {\partial u_{z}}{\partial \theta }+u_{z}\frac {\partial u_{z}}{\partial z}, \\[-26pt] \nonumber \end{align}

(2.14) \begin{align} {N}_{T}=u_{r}\frac {\partial T}{\partial r}+\frac {u_{\theta }}{r}\frac {\partial T}{\partial \theta }+u_{z}\frac {\partial T}{\partial z}. \\[6pt] \nonumber \end{align}

The parameters $Re$ , $N$ and $Pr$ are the Reynolds number, non-dimensional Brunt–Väisälä frequency and Prandtl number, respectively, which are defined as

(2.15) \begin{align} Re=\frac {R_{i}^{2}\varOmega _{i}}{\nu _{0}},\quad N=\sqrt {\frac {\alpha _{0}g}{\varOmega _{i}^{2}}\frac {\Delta T}{\Delta z}},\quad Pr=\frac {\nu _{0}}{\kappa _{0}}. \end{align}

For convenience in comparison with other literature, we also use a conventional inner-cylinder-based Reynolds number $Re_{i}$ defined as

(2.16) \begin{align} Re_{i}=\frac {R_{i}\varOmega _{i}d}{\nu _{0}}=Re\left (\frac {1-\eta }{\eta }\right ), \end{align}

where $d=R_{o}-R_{i}$ is the gap size between the two cylinders.

2.2. Pseudo-spectral formulation for direct numerical simulation

In this study, spectral-based DNS are conducted to solve (2.6)–(2.10) numerically. To do so, we decompose the perturbation into the sum of modes using the following Fourier representation:

(2.17) \begin{align} \left ( \begin{array}{c} u_{r}\\ u_{\theta }\\ u_{z}\\ T\\ p \end{array} \right ) = \sum _{j=-M}^{M} \sum _{l=-K}^{K} \left ( \begin{array}{c} \tilde {u}_{\! jl}(r,t)\\ \tilde {v}_{\! jl}(r,t)\\ \tilde {w}_{\! jl}(r,t)\\ \tilde {T}_{\! jl}(r,t)\\ \tilde {p}_{\! jl}(r,t) \end{array} \right )\exp \left (\mathrm {i}m_{\! j}\theta +\mathrm {i}k_{l}z\right ), \end{align}

where $M$ and $K$ are the cut-off numbers of modes considered in the azimuthal direction $\theta$ and axial direction $z$ , respectively, $\tilde {u}_{\! jl}$ , $\tilde {v}_{\! jl}$ , $\tilde {w}_{\! jl}$ , $\tilde {p}_{\! jl}$ and $\tilde {T}_{\! jl}$ are the time-dependent mode shapes, $m_{j}=jm$ is the $j$ th multiple of the principal azimuthal wavenumber $m$ and $k_{l}=lk$ is the $l$ th multiple of the principal axial wavenumber $k$ . The above ansatz (2.17) is used for various problems such as the Taylor–Couette formulation in nsCouette code (López et al. Reference López, Feldmann, Rampp, Vela-Martín, Shi and Avila2020) or convection problems (Saltzman Reference Saltzman1962; Park et al. Reference Park, Moon, Seo and Baik2021a ). The formulation is also analogous to semi-linear models, which are applied to centrifugal instability of anti-cyclonic vortices (Yim et al. Reference Yim, Billant and Gallaire2020, Reference Yim, Billant and Gallaire2023). The semi-linear theory allows us to investigate directly nonlinear interaction between base flow and an instability mode, and the method is generalised by considering nonlinear interaction among multiple instability modes in low-order harmonics while neglecting the triad interaction leading to harmonics of orders higher than the cut-off numbers. For each mode with indices $j$ and $l$ , we apply the ansatz (2.17) to (2.6)–(2.10) and obtain

(2.18) \begin{align} \frac {\partial \tilde {u}_{\! jl}}{\partial r}+\frac {\tilde {u}_{\! jl}}{r}+\frac {\mathrm {i}m_{\! j}\tilde {v}_{\! jl}}{r}+\mathrm {i}k_{l}\tilde {w}_{\! jl}=0, \\[-26pt] \nonumber \end{align}

(2.19) \begin{align} \frac {\partial \tilde {u}_{\! jl}}{\partial t}+\mathrm {i}m_{\! j}\varOmega \tilde {u}_{\! jl}-2\varOmega \tilde {v}_{\! jl}+\tilde {N}_{r,jl}=-\frac {\partial \tilde {p}_{\! jl}}{\partial r}+\frac {1}{Re}\left (\tilde {\nabla }_{jl}^{2}\tilde {u}_{\! jl}-\frac {\tilde {u}_{\! jl}}{r^{2}}-\frac {2\mathrm {i}m_{\! j}\tilde {v}_{\! jl}}{r^{2}}\right ), \\[-26pt] \nonumber \end{align}

(2.20) \begin{align} \frac {\partial \tilde {v}_{\! jl}}{\partial t}+\mathrm {i}m_{\! j}\varOmega \tilde {v}_{\! jl}+Z\tilde {u}_{\! jl}+\tilde {N}_{\theta ,jl}=-\frac {\mathrm {i}m_{\! j}\tilde {p}_{\! jl}}{r}+\frac {1}{Re}\left (\tilde {\nabla }_{jl}^{2}\tilde {v}_{\! jl}-\frac {\tilde {v}_{\! jl}}{r^{2}}+\frac {2\mathrm {i}m_{\! j}\tilde {u}_{\! jl}}{r^{2}}\right ), \\[-26pt] \nonumber \end{align}

(2.21) \begin{align} \frac {\partial \tilde {w}_{\! jl}}{\partial t}+\mathrm {i}m_{\! j}\varOmega \tilde {w}_{\! jl}+\tilde {N}_{z,jl}=-\mathrm {i}k_{l}\tilde {p}_{\! jl}+N^{2}\tilde {T}_{\! jl}+\frac {1}{Re}\tilde {\nabla }_{jl}^{2}\tilde {w}_{\! jl}, \\[-26pt] \nonumber \end{align}

(2.22) \begin{align} \frac {\partial \tilde {T}_{\! jl}}{\partial t}+\mathrm {i}m_{\! j}\varOmega \tilde {T}_{\! jl}+\tilde {w}_{\! jl}+\tilde {N}_{T,jl}=\frac {1}{RePr}\tilde {\nabla }_{jl}^{2}\tilde {T}_{\! jl}, \end{align}

where $\tilde {N}_{r,jl}$ , $\tilde {N}_{\theta ,jl}$ , $\tilde {N}_{z,jl}$ and $\tilde {N}_{T,jl}$ are the terms convoluted from the nonlinear terms (2.11)–(2.14) using the Fourier transform and $\tilde {\nabla }_{jl}^{2}=\partial ^{2}/\partial r^{2}+(1/r)\partial /\partial r-k_{l}^{2}-m_{j}^{2}/r^{2}$ is the modal Laplacian operator. For the boundary conditions, we consider

(2.23) \begin{align} \tilde {u}_{\! jl}=\tilde {v}_{\! jl}=\tilde {w}_{\! jl}=\frac {\partial \tilde {T}_{\! jl}}{\partial r}=0, \end{align}

at both cylinders $r=1$ and $r=1/\eta$ .

An advantage of using the equations in a modal form is that (2.18)–(2.22) can further be simplified by eliminating the pressure $\tilde {p}_{\! jl}$ as follows:

(2.24) \begin{align} \mathcal {A}_{jl}\frac {\partial \tilde {\boldsymbol{q}}_{\! jl}}{\partial t}=\mathcal {B}_{jl}\tilde {\boldsymbol{q}}_{\! jl}+\tilde {\boldsymbol{N}}_{\! jl}, \end{align}

where $\tilde {\boldsymbol{q}}_{\! jl}$ is the vector of three variables, $\mathcal {A}_{jl}$ and $\mathcal {B}_{jl}$ are the differential operator matrices and $\tilde {\boldsymbol{N}}_{\! jl}$ is the vector comprised of the nonlinear terms, all of which depend on indices $j$ and $l$ such as whether $l=0$ or not and whether $j=0$ or not. We refer the reader to Appendix A for more details of these vectors and matrices. For spatial discretisation in the radial direction $r$ , a spectral method using the Chebyshev collocation points is considered (Antkowiak & Brancher Reference Antkowiak and Brancher2004; Park Reference Park2012). The vector $\tilde {\boldsymbol{N}}_{\! jl}$ is computed at each step using a pseudo-spectral method, similar to that in Deloncle et al. (Reference Deloncle, Billant and Chomaz2008). Equation (2.24) is written in a spectral form in the azimuthal and axial directions and we use the backward Fourier transform to obtain the vector $\boldsymbol{q}_{\! jl}$ in the physical space $(r,\theta ,z)$ . Then $\boldsymbol{q}_{\! jl}$ is used to compute the nonlinear term $\boldsymbol{N}_{\! jl}$ in the physical space and apply the forward Fourier transform to $\boldsymbol{N}_{\! jl}$ to compute $\tilde {\boldsymbol{N}}_{\! jl}$ in the spectral space $(r,m_{j},k_{l})$ . In this pseudo-spectral approach, we consider the numbers of collocation points $N_{\theta }=2M+1$ and $N_{z}=2K+1$ in the azimuthal and axial directions, respectively. The pseudo-spectral approach is similar to that in Guseva et al. (Reference Guseva, Willis, Hollerbach and Avila2015) who also consider the axial periodicity with the Fourier method but utilise a fractional time-stepping method by carefully computing an intermediate pressure using the influence-matrix method, which is especially effective in solving equations for the magnetic field. On the contrary, we avoid the use of the pressure $\tilde {p}_{\! jl}$ by considering different operator matrices $\mathcal {A}_{jl}$ and $\mathcal {B}_{jl}$ that depend on the indices $j$ and $l$ . Unlike Deloncle et al. (Reference Deloncle, Billant and Chomaz2008), the de-aliasing technique, which is required in turbulence simulations, is not implemented as the Reynolds number $Re_{i}$ considered in this work is not large enough to observe small-scale turbulence. For time advancement, we use a conventional semi-implicit scheme (see e.g. Kim et al. Reference Kim, Moin and Moser1987) with the Crank–Nicolson method for the linear term and the Adam–Bashforth method for the nonlinear term. For instance, at each step $n$ , we find the next step solution $\tilde {\boldsymbol{q}}_{\! jl}^{(n+1)}$ by solving numerically the discretised version of (2.24):

(2.25) \begin{align} \left (\mathcal {A}_{jl}-\frac {\Delta t}{2}\mathcal {B}_{jl}\right )\tilde {\boldsymbol{q}}_{\! jl}^{(n+1)}=\left (\mathcal {A}_{jl}+\frac {\Delta t}{2}\mathcal {B}_{jl}\right )\tilde {\boldsymbol{q}}_{\! jl}^{(n)}+\frac {\Delta t}{2}\left (3\tilde {\boldsymbol{N}}_{\! jl}^{(n)}-\tilde {\boldsymbol{N}}_{\! jl}^{(n-1)}\right ), \end{align}

where $\Delta t$ is the time step. In this study, the time step $\Delta t$ is fixed to $\Delta t=0.01$ , which is found to be sufficiently small for parameters considered in this study. It is verified that the Courant–Friedrichs–Lewy condition for perturbation velocity is satisfied and every DNS demonstrates no numerical divergence with this $\Delta t$ .

2.3. One-dimensional local linear stability analysis

By assuming that the perturbation is infinitesimally small, we can neglect nonlinear terms and perform 1-D local LSA using the normal mode:

(2.26) \begin{align} \left ( \begin{array}{c} u_{r}\\ u_{\theta }\\ u_{z}\\ T\\ p \end{array} \right ) = \left ( \begin{array}{c} \hat {u}(r)\\ \hat {v}(r)\\ \hat {w}(r)\\ \hat {T}(r)\\ \hat {p}(r) \end{array} \right )\exp \left (\mathrm {i}m\theta +\mathrm {i}kz-\mathrm {i}\omega t\right )+\text{c.c.}, \end{align}

where $\text{c.c.}$ denotes the complex conjugate, $\hat {u}$ , $\hat {v}$ , $\hat {w}$ , $\hat {p}$ and $\hat {T}$ are the mode shapes, $m$ is the azimuthal wavenumber, $k$ is the axial wavenumber and $\omega$ is the complex frequency $\omega =\omega _{r}+\mathrm {i}\omega _{i}$ , where $\omega _{r}={\textrm {Re}}(\omega )$ is the temporal frequency and $\omega _{i}={\textrm {Im}}(\omega )$ is the temporal growth rate. Throughout the paper, the non-dimensional wavenumber $k_{d}=kd/R_{i} = k(1-\eta )/\eta$ rescaled by the gap size $d$ is also used for convenience in comparison with other literature. The normal mode (2.26) is applied to (2.6)–(2.10) with the nonlinear terms neglected and, after manipulations, we can eliminate the pressure mode shape and obtain the following simplified eigenvalue problem:

(2.27) \begin{align} -\mathrm {i}\omega \mathcal {A}\hat {\boldsymbol{q}}=\mathcal {B}\hat {\boldsymbol{q}}, \end{align}

where $\hat {\boldsymbol{q}}= (\hat {u},\hat {v},\hat {T} )^{\mathrm {T}}$ and $\mathcal {A}$ and $\mathcal {B}$ are the operator matrices, which are essentially the same as $\mathcal {A}_{11}$ and $\mathcal {B}_{11}$ , respectively, the operator matrices in (2.24) with $j=l=1$ . To solve the eigenvalue problem (2.27), the MATLAB routine eig is used with the following boundary conditions imposed at both cylinders $r=1$ and $r=1/\eta$ :

(2.28) \begin{align} \hat {u}=\hat {v}=\hat {w}=\frac {\mathrm {d}\hat {T}}{\mathrm {d}r}=0. \end{align}

The Chebyshev spectral method is used for discretisation in the radial direction $r$ and a number of collocation points $N_{r}$ between 60 and 120 is considered in this study. The choice of $N_{r}$ depends on how large are the parameters such as the Reynolds number $Re$ or the Péclet number $Pe=RePr$ , how the eigenmodes are confined near the boundaries, etc. For more details of the LSA and numerical methods, we refer the reader to our previous work (e.g. Park Reference Park2012; Park & Billant Reference Park and Billant2013) that used the same code.

2.4. Two-dimensional bi-global linear stability analysis

For parameters considered in this study, 1-D LSA reveals that the axisymmetric mode with $m=0$ is the most unstable one for cylindrical Couette flow. As demonstrated in the next sections, centrifugal instability develops nonlinearly and saturates leading to axisymmetric Taylor vortices as a new base state. The new 2-D base state $\bar {\boldsymbol{Q}}(r,z)=(\bar {\boldsymbol{U}},\bar {T})$ is comprised of the new base velocity $\bar {\boldsymbol{U}}(r,z)= (\bar {U}(r,z),\bar {V}(r,z),\bar {W}(r,z) )$ and the new base temperature $\bar {T}(r,z)$ . For this Taylor vortex flow, we can analyse its secondary instability through 2-D bi-global LSA by considering a non-axisymmetric perturbation $\bar {\textbf {q}}$ with $m\neq 0$ expressed in the following ansatz:

(2.29) \begin{align} \bar {\textbf {q}}= \left ( \begin{array}{c} \bar {u}_{r}\\ \bar {u}_{\theta }\\ \bar {u}_{z}\\ \bar {T}\\ \bar {p} \end{array} \right ) = \left ( \begin{array}{c} \hat {u}_{m}(r,z)\\ \hat {v}_{m}(r,z)\\ \hat {w}_{m}(r,z)\\ \hat {T}_{m}(r,z)\\ \hat {p}_{m}(r,z) \end{array} \right )\exp \left (\mathrm {i}m\theta -\mathrm {i}\omega _{m} t\right )+\text{c.c.}, \end{align}

where $\hat {u}_{m}$ , $\hat {v}_{m}$ , $\hat {w}_{m}$ , $\hat {T}_{m}$ and $\hat {p}_{m}$ are the 2-D global mode shapes and $\omega _{m}$ is the complex frequency of the global mode. Similar to the 1-D local LSA, an eigenvalue problem can be formulated as

(2.30) \begin{align} -\mathrm {i}\omega _{m}\mathcal {A}_{m}\hat {\textbf {q}}_{m}=\mathcal {B}_{m}\hat {\textbf {q}}_{m}, \end{align}

where $\mathcal {A}_{m}$ and $\mathcal {B}_{m}$ are the operator matrices detailed in Appendix A.2 and $\hat {\textbf {q}}_{m}= (\hat {u}_{m},\hat {v}_{m},\hat {T}_{m} )^{\mathrm {T}}$ is the mode shape. The matrices $\mathcal {A}_{m}$ and $\mathcal {B}_{m}$ in the 2-D bi-global LSA have a much larger size than the size of the matrices $\mathcal {A}$ and $\mathcal {B}$ in (2.27). Thus we use the MATLAB routine eigs that computes only a few eigenvalues near a specified value using the Krylov–Schur algorithm (Stewart Reference Stewart2002). To solve the eigenvalue problem (2.30), we consider the following boundary conditions:

(2.31) \begin{align} \hat {u}_{m}=\hat {v}_{m}=\hat {w}_{m}=\frac {\partial \hat {T}_{m}}{\partial r}=0. \end{align}

The Chebyshev and Fourier spectral methods are used for discretisation in the radial and axial directions, respectively. The Fourier method imposes the periodic boundary condition in the axial direction $z$ given that the axisymmetric Taylor vortices as a new base state are also periodic in $z$ with wavelength $\lambda _{z}=2\pi /k$ .

3.1. Neutral stability curves

Figure 2(a) shows the growth rate $\omega _{i}$ of the axisymmetric mode ( $m=0$ ), which is found to be most unstable, versus the wavenumber $k_{d}$ for different parameter sets of $(N,Pr)$ at $\mu =0$ , $\eta =0.9$ and $Re_{i}=200$ . For every case presented in figure 2(a), the corresponding frequency $\omega _{r}$ is zero. It is known that centrifugal instability reaches its maximum $\omega _{i,max }$ as $k\rightarrow \infty$ in the inviscid limit $Re\rightarrow \infty$ as the growth rate scales as $\omega _{i,{inviscid}}=\omega _{i,max } -A_{0}/k^{2}$ (Billant & Gallaire Reference Billant and Gallaire2005; Park et al. Reference Park, Billant and Baik2017) while the viscous growth rate scales as $\omega _{i,{viscous}}=\omega _{i,{inviscid}}-A_{1}k^{2}/Re$ (Yim et al. Reference Yim, Billant and Ménesguen2016), where $A_{0}$ and $A_{1}$ are positive constants that depend on stratification and other parameters. Due to this characteristic varying with $k$ , the viscous growth rates are positive in a finite range of $k_{d}$ and each curve reaches its peak at a certain wavenumber $k_{d,max }$ . As the stratification increases from $N=0$ (grey) to $N=1$ (black) for $Pr=1$ , the growth-rate curves descend while the wavenumber $k_{d,max }$ at the peak of $\omega _{i}$ increases. For fixed $N=1$ , as the Prandtl number $Pr$ decreases from $Pr=1$ , the growth rate increases and, remarkably, the growth-rate curves for $Pr=10^{-2}$ and $10^{-4}$ overlap with the grey curve, the unstratified case with $N=0$ . This implies that the centrifugal instability in stratified and highly diffusive fluids with $Pr\ll 1$ behaves as the instability in unstratified fluids as the effect of stratification is suppressed by strong thermal diffusion. Another remarkable result in figure 2(a) is that the growth-rate curves overlap for $(N,Pr)=(10,10^{-2})$ and $(10^{2},10^{-4})$ , the cases with the same $P_{N}=N^{2} Pr=1$ . The growth rate for other values of $N$ and $Pr\ll 1$ is invariant if the rescaled parameter $P_{N}$ is the same. Such an invariant feature at the same $P_{N}$ but different $N$ and $Pr$ is similarly reported for both vertical and horizontal shear instabilities in vertically stratified fluids (Lignières et al. Reference Lignières1999; Park et al. Reference Park, Prat and Mathis2020, Reference Park, Prat and Mathis2021b ).

Figure 2.

(a) Growth-rate curves for various sets of $(N,Pr)$ at $\mu =0$ , $\eta =0.9$ , $Re_{i}=200$ and $m=0$ . (b) Neutral stability curves in the parameter space $(N,Re_{i})$ for different $Pr$ at $\mu =0$ , $\eta =0.9$ and $m=0$ . (c) The curves for different $Pr$ same as (b) but over a wider range of $N$ . (d) The curves same as (c) overlapped due to the rescaled parameter $P_{N}=N^{2}Pr$ on the abscissa. An additional thick grey line is a neutral stability curve obtained from the small- $Pr$ approximation.

Figure 2 $(b)$ display neutral stability curves, which denote the critical Reynolds number $Re_{i,c}$ at which the growth rate $\omega _{i}$ of the most unstable mode is zero, in the parameter space $(N,Re_{i})$ for different values of $Pr$ at $\mu =0$ , $\eta =0.9$ and $m=0$ . For $Pr\geqslant 1$ , the neutral stability curves increase rapidly with $N$ , a feature similarly found in Park et al. (Reference Park, Billant and Baik2017), and the $Pr=1$ case is more stable than other cases with $Pr\gt 1$ when $N\gt 3$ . This is expected as thermal dissipation is proportional to the diffusivity $\kappa _{0}$ ; thus lower diffusivity $\kappa _{0}$ (i.e. higher $Pr$ ) implies less thermal dissipation and more instability. In the range $0\lt N\lt 3$ , the situation is more complicated since the $Pr=10$ case is more stable than other $Pr$ cases for the axisymmetric perturbation. Although such high- $Pr$ dynamics at a moderate $N$ should be further investigated due to its great importance in geophysical and other contexts (e.g. oceanic flows where the analogous Schmidt number $Sc$ is around 700), the current study will only focus on the highly diffusive regime with low $Pr\leqslant 1$ . For $Pr\lt 1$ , the curves increase very slowly as $N$ increases and it is difficult to see in figure 2(b) the increase of the curves at $Pr=10^{-4}$ and $10^{-6}$ over the range $0\leqslant N\leqslant 20$ . Figure 2(c) displays the same neutral stability curves for $Pr\leqslant 1$ over a wider range of $N$ to see how slowly the neutral stability curves increase with $N$ as $Pr$ decreases. This increasing trend confirms that, while stratification enhances the stability of stratified Taylor–Couette flow, strong thermal diffusion destabilises by preventing the stabilising role of stratification and making the flow behave like an unstratified flow. In figure 2(d), we plot again the same neutral stability curves but over the rescaled parameter $P_{N}$ on the abscissa. All the curves now overlap each other, even for the case with $Pr=1$ , and can be described by a linear relation as $Re_{i,c}= 131.6+8.965P_{N}$ , where $Re_{i,c}=131.6$ is the critical Reynolds number for the unstratified case $N=0$ at $m=0$ and $\eta =0.9$ , the number agreeing with DiPrima et al. (Reference DiPrima, Eagles and Ng1984).

Figure 3.

(a,b) Real (solid) and imaginary (dashed) parts of the mode shape $\hat {v}(r)$ and rescaled mode $\hat {T}(r)/Pr$ for $Pr=1$ (red), $Pr=10^{-2}$ (blue) and $Pr=10^{-4}$ (grey) at $\mu =0$ , $\eta =0.9$ , $Re_{i}=200$ , $N=1$ , $m=0$ and $k_{d}=3.91$ . (c) Perturbation temperature $T(r,z)$ reconstructed from $\hat {T}(r)$ for $Pr=10^{-2}$ .

Figure 3 shows the real and imaginary parts of eigenmode shapes $\hat {v}(r)$ and $\hat {T}(r)$ for different $Pr$ for the wavenumber set $(m,k_{d})=(0,3.91)$ . The mode shapes are normalised by the maximum value of $\hat {v}$ and it is found that the mode shape $\hat {T}$ has a smaller amplitude than $\hat {v}$ for small $Pr$ . As the Prandtl number $Pr$ decreases, the amplitude of $\hat {T}$ decreases and scales as $O(Pr)$ . For a better comparison, figure 3(b) displays the mode shape $\hat {T}(r)$ divided by $Pr$ and we see that the rescaled $\hat {T}/Pr$ has the same mode shape for the cases with $Pr=10^{-2}$ and $Pr=10^{-4}$ . In figure 3(c), the perturbation temperature $T(r,z)=\hat {T}\exp (\mathrm {i}kz)+\hat {T}^{*}\exp (-\mathrm {i}kz)$ for $Pr=10^{-2}$ is plotted. Two in-phase waves, one near the inner cylinder and the other near the outer one, are clearly shown. The shape of this axisymmetric mode is different from that of a non-axisymmetric centrifugal instability mode, which is weakly sheared, or that of a SRI mode that is out of phase (Park et al. Reference Park, Billant and Baik2017, Reference Park, Billant, Baik and Seo2018). The perturbation temperature has one node (i.e. the zero crossing) around $r\simeq 1.053$ as it is the first mode which becomes the most unstable. Higher-order modes with a greater number of nodes are, on the other hand, found to be stable for the parameter set in figure 3 with $Pr=10^{-2}$ .

3.2. Small- $Pr$ approximation

The single dependence on $P_{N}$ for different $N$ and $Pr$ in the limit $Pr\ll 1$ can be understood by taking the small- $Pr$ approximation. Consider the Taylor expansions in terms of small $Pr$ as $\hat {\textbf {u}}= (\hat {u},\hat {v},\hat {w} )^{\mathrm {T}}=\hat {\textbf {u}}^{(0)}+Pr\hat {\textbf {u}}^{(1)}+O (Pr^{2} )$ and similarly for $\hat {T}$ and $\hat {p}$ as $\hat {T}=\hat {T}^{(0)}+Pr\hat {T}^{(1)}+\cdots$ and $\hat {p}=\hat {p}^{(0)}+Pr\hat {p}^{(1)}+\cdots$ . For the Péclet number $Pe$ being small, $Pe=RePr\ll 1$ (i.e. $Re\sim O(1)$ and $Pr\ll 1$ ), we obtain the leading-order equation for $\hat {T}^{(0)}$ as

(3.1) \begin{align} \frac {1}{RePr}\hat {\nabla }^{2}\hat {T}^{(0)}=0, \\[-6pt] \nonumber \end{align}

which leads to the analytic solution $\hat {T}^{(0)}(r)=a_{1}\mathrm {I}_{m}(kr)+a_{2}\mathrm {K}_{m}(kr)$ , where $a_{1}$ and $a_{2}$ are constants and $\mathrm {I}_{m}$ and $\mathrm {K}_{m}$ are modified Bessel functions of the first and second kinds, respectively (Abramowitz & Stegun Reference Abramowitz and Stegun1965). The constants $a_{1}$ and $a_{2}$ become zero if the no-flux conditions $\mathrm {d}\hat {T}^{(0)}/\mathrm {d}r=0$ are imposed at both cylinders; thus the leading-order solution simply becomes zero (i.e. $\hat {T}^{(0)}(r)=0$ ). At the next order, we obtain the following equations:

(3.2) \begin{align} \frac {\mathrm {d} \hat {u}^{(0)}}{\mathrm {d} r}+\frac {\hat {u}^{(0)}}{r}+\frac {\mathrm {i}m\hat {v}^{(0)}}{r}+\mathrm {i}k\hat {w}^{(0)}=0, \\[-24pt] \nonumber \end{align}

(3.3) \begin{align} -\mathrm {i}\omega \hat {u}^{(0)}+\mathrm {i}m\varOmega \hat {u}^{(0)}-2\varOmega \hat {v}^{(0)}=-\frac {\mathrm {d} \hat {p}^{(0)}}{\mathrm {d} r}+\frac {1}{Re}\left (\hat {\nabla }^{2}\hat {u}^{(0)}-\frac {\hat {u}^{(0)}}{r^{2}}-\frac {2\mathrm {i}m\hat {v}^{(0)}}{r^{2}}\right ), \\[-24pt] \nonumber \end{align}

(3.4) \begin{align} -\mathrm {i}\omega \hat {v}^{(0)}+\mathrm {i}m\varOmega \hat {v}^{(0)}+Z\hat {u}^{(0)}=-\frac {\mathrm {i}m\hat {p}^{(0)}}{r}+\frac {1}{Re}\left (\hat {\nabla }^{2}\hat {v}^{(0)}-\frac {\hat {v}^{(0)}}{r^{2}}+\frac {2\mathrm {i}m\hat {u}^{(0)}}{r^{2}}\right ), \\[-24pt] \nonumber \end{align}

(3.5) \begin{align} -\mathrm {i}\omega \hat {w}^{(0)}+\mathrm {i}m\varOmega \hat {w}^{(0)}=-\mathrm {i}k\hat {p}^{(0)}+P_{N}\hat {T}^{(1)}+\frac {1}{Re}\hat {\nabla }^{2}\hat {w}^{(0)}, \\[-24pt] \nonumber \end{align}

(3.6) \begin{align} \hat {w}^{(0)}=\frac {1}{Re}\hat {\nabla }^{2}\hat {T}^{(1)}, \\[-6pt] \nonumber \end{align}

where $P_{N}=N^{2}Pr$ . We see that the two parameters $N$ and $Pr$ representing thermal stratification and diffusion are simplified into a single parameter $P_{N}$ , as observed in other studies under the small-Péclet-number approximation (e.g. Lignières Reference Lignières1999; Park et al. Reference Park, Prat and Mathis2020). By considering the continuity equation and eliminating the pressure, (3.2)–(3.6) can further be simplified into the following eigenvalue problem:

(3.7) \begin{align} -\mathrm {i}\omega \mathcal {A}^{(0)}\hat {\boldsymbol{q}}^{(0)}=\mathcal {B}^{(0)}\hat {\boldsymbol{q}}^{(0)}, \end{align}

where $\hat {\boldsymbol{q}}^{(0)}$ is the mode shape vector and $\mathcal {A}^{(0)}$ and $\mathcal {B}^{(0)}$ are the operator matrices detailed in Appendix A. In figure 2(d), it is clearly shown that a neutral stability curve computed from (3.7) overlaps with other neutral stability curves when they are plotted over the rescaled Prandtl number $P_{N}=N^{2}Pr$ as the abscissa.

3.3. Perturbation energy analysis

The instability characteristics can also be understood by examining the evolution of perturbation energy. We first define the total perturbation energy $E(t)$ as

(3.8) \begin{align} E(t)=\frac {1}{2}\int _{0}^{L_{z}}\int _{0}^{2\pi }\int _{1}^{1/\eta }\left (u_{r}^{2}+u_{\theta }^{2}+u_{z}^{2}+N^{2}T^{2}\right )r\,\mathrm {d}r\,\mathrm {d}\theta \,\mathrm {d}z=\frac {1}{2}\left \lt \textbf {q}_{E}\,\textbf {:}\,\textbf {q}_{E}\right \gt , \end{align}

where $L_{z}$ is a vertical length of the domain assumed to be periodic (e.g. $L_{z}=2\pi /k$ if we consider one periodic length of a normal mode with the axial wavenumber $k$ ), the angle brackets denote the volume integral defined as $\left \lt \textbf {X}\right \gt =\int _{0}^{L_{z}}\int _{0}^{2\pi }\int _{1}^{1/\eta }\textbf {X} \,r\mathrm {d}r\mathrm {d}\theta \mathrm {d}z$ , $\textbf {q}_{E}= (u_{r},u_{\theta },u_{z},NT )^{\mathrm {T}}$ and $\textbf {:}$ denotes the Frobenius product (for more details, see also Park et al. Reference Park, Billant and Baik2017). For the momentum and energy equations (2.7)–(2.10), we multiply both sides by $ (u_{r}, u_{\theta }, u_{z}, N^{2}T )^{\mathrm {T}}$ and, after some manipulations, we obtain the following equation for perturbation energy:

(3.9) \begin{align} \frac {\partial E}{\partial t}=\left \lt -\frac {\mathrm {d}\varOmega }{\mathrm {d}r}\left (ru_{r}u_{\theta }\right )-\frac {1}{Re}\left [\nabla \textbf {u}\,\textbf {:}\,\nabla \textbf {u}+\frac {N^{2}}{Pr}\left (\nabla T\,\textbf {:}\,\nabla T\right )\right ]\right \gt , \end{align}

where the first term on the right-hand side within the angle brackets denotes the contribution from the mean angular shear $\mathrm {d}\varOmega /\mathrm {d}r$ and the second term corresponds to the kinetic and potential energy dissipation. The contribution from the mean angular shear is a main source of energy production if the momentum transfer term $u_{r}u_{\theta }$ is anti-correlated with the mean angular shear $\mathrm {d}\varOmega /\mathrm {d}r$ , a well-known mechanism for instability called the Orr mechanism (Orr Reference Orr1907). The momentum and thermal dissipation terms are always negative and thus stabilise the perturbation energy. We note that this mechanism is valid for both linear and nonlinear cases as the evolution equation (3.9) is derived from the nonlinear perturbation equations (2.7)–(2.10) and the nonlinear terms are cancelled out in the derivation process in which continuity is taken into account. If we apply the normal mode (2.26) to the evolution equation (3.9), we obtain the following expression for the growth rate:

(3.10) \begin{align} \omega _{i}=\frac {1}{\hat {E}}\left \lt -\frac {\mathrm {d}\varOmega }{\mathrm {d}r}\frac {r\left (\hat {u}^{*}\hat {v}+\hat {u}\hat {v}^{*}\right )}{2}-\frac {1}{Re}\left [\hat {\nabla }\hat {\textbf {u}}\,\textbf {:}\,\hat {\nabla }\hat {\textbf {u}}+\frac {N^{2}}{Pr}\left (\hat {\nabla }\hat {T}\,\textbf {:}\,\hat {\nabla }\hat {T}\right )\right ]\right \gt _{r}, \end{align}

where $*$ denotes the complex conjugate, $\lt \gt _{r}$ denotes the line integral over the radial coordinate $r$ as $\left \lt \hat {\textbf {X}}\right \gt _{r}=\int _{1}^{1/\eta }\hat {\textbf {X}} \,r\mathrm {d}r$ and $\hat {E}=\left \lt |\hat {u}|^{2}+|\hat {v}|^{2}+|\hat {w}|^{2}+|N\hat {T}|^{2}\right \gt _{r}$ .

Table 1 presents examples of the maximum growth rate $\omega _{i,max }$ and corresponding wavenumber $k_{d,max }$ for various parameter sets of $(N,Pr,m)$ at $\mu =0$ , $\eta =0.9$ and $Re_{i}=200$ . The table also details the production and dissipation terms $\mathcal {P}_{\Omega }$ , $\epsilon _{k}$ and $\epsilon _{p}$ , which are defined as

(3.11) \begin{align} \mathcal {P}_{\Omega }=\left \lt -\frac {\mathrm {d}\varOmega }{\mathrm {d}r}\frac {r\left (\hat {u}^{*}\hat {v}+\hat {u}\hat {v}^{*}\right )}{2}\right \gt _{r},\quad \epsilon _{k}=-\frac {1}{Re}\left \lt \hat {\nabla }\hat {\textbf {u}}:\hat {\nabla }\hat {\textbf {u}}\right \gt _{r},\quad \epsilon _{p}=-\frac {N^{2}}{RePr}\left \lt \hat {\nabla }\hat {T}\,\textbf {:}\,\hat {\nabla }\hat {T}\right \gt _{r}. \end{align}

The production and dissipation terms are calculated by implementing the eigenfunction into the expressions (3.11) and their sum shows a good agreement with the eigenvalue $\omega _{i,max }$ . For every case, the contribution from the thermal dissipation $\epsilon _{p}$ to the growth rate is small as $Pr\leqslant 1$ and the growth rate $\omega _{i}$ is mainly determined by a difference between the production $\mathcal {P}_{\Omega }$ and the viscous dissipation $\epsilon _{k}$ . For the cases with $Pr\ll 1$ , if we apply the small- $Pr$ approximation, we can re-express the thermal dissipation $\epsilon _{p}$ as

(3.12) \begin{align} \epsilon _{p}\simeq -\frac {P_{N}}{Re}\left \lt \hat {\nabla }\hat {T}_{1}\,\textbf {:}\,\hat {\nabla }\hat {T}_{1}\right \gt _{r}\sim O\left (\frac {P_{N}}{Re}\right ), \end{align}

where the dependence on the two parameters $N$ and $Pr$ is now expressed by the single parameter $P_{N}$ under the small- $Pr$ approximation. The scaling for thermal dissipation $\epsilon _{p}\sim O(P_{N}/Re)$ is based on the assumption $\hat {T}_{1}\sim O(1)$ under the small- $Pr$ approximation.

Table 1.

Values of the maximum growth rates $\omega _{i,max }$ with the corresponding wavenumber $k_{d,max }$ , production and dissipation terms for various $N$ , $Pr$ and $m$ at $\mu =0$ , $\eta =0.9$ and $Re_{i}=200$ .

3.4. Parametric investigations

Figure 4.

Neutral stability curves for different $m$ for (a) $Pr=1$ and (b) $Pr=10^{-4}$ at $\mu =0$ and $\eta =0.9$ . (c) Neutral stability curves for different $Pr$ over the rescaled parameter $P_{N}$ at $\mu =0$ , $\eta =0.9$ and $m=2$ . A thick grey line denotes the neutral stability curve from the small- $Pr$ approximation (SPA).

Figures 4(a) and 4(b) show neutral stability curves for axisymmetric and non-axisymmetric cases in the parameter space $(N,Re_{i})$ at $\mu =0$ and $\eta =0.9$ for $Pr=1$ and $10^{-4}$ , respectively. For both $Pr$ values, all the neutral stability curves increasing with $N$ ascend as $m$ increases; thus the lowest curves correspond to the axisymmetric case with $m=0$ . A difference between the two Pr values is that, for $Pr=1$ , neutral stability curves for non-axisymmetric cases ( $m\gt 0$ ) stay closer to the curve for the axisymmetric case ( $m=0$ ) as $N$ increases, while for $Pr=10^{-4}$ , the curves for $m\gt 0$ stay further from the curve for $m=0$ as $N$ increases. This implies that for strongly stratified fluids with $N\gg 0$ for $Pr=1$ , the axisymmetric mode will appear at the instability onset $Re_{i}=Re_{i,c}$ but then non-axisymmetric perturbations will also become unstable immediately after the onset; thus competition between the axisymmetric and non-axisymmetric modes will occur right above $Re_{i}\gt Re_{i,c}$ . On the other hand, for low $Pr$ and above the instability onset $Re_{i}\gt Re_{i,c}$ , the axisymmetric mode will become more dominant over non-axisymmetric modes as $N$ increases. Figure 4(c) displays neutral stability curves over the rescaled $P_{N}$ for different $Pr$ at $m=2$ . As similarly observed for the axisymmetric case $m=0$ in figure 2(c), the curves for $Pr\leqslant 10^{-2}$ overlap and agree with the prediction from the small- $Pr$ approximation. It is found that the line from the small- $Pr$ approximation for $m=2$ scales as $Re_{i,c}= 134.4+9.971P_{N}$ , where $Re_{c}=134.4$ as the $Re_{i}$ -axis intercept at $P_{N}=0$ corresponds to the critical Reynolds number for unstratified case $N=0$ . The slope 9.971 for $m=2$ is higher than the slope 8.965 of the $m=0$ case. It is verified that the critical Reynolds number and the slope increase with $m$ and thus they are at the lowest for $m=0$ . This implies that for low $Pr$ , the axisymmetric mode with $m=0$ is expected to be the most unstable one for $P_{N}=N^{2}Pr\geqslant 0$ .

Figure 5.

Neutral stability curves obtained from the small- $Pr$ approximation at $\mu =0$ .

Figure 5 shows neutral stability curves obtained from the small- $Pr$ approximation in the wide parameter space $(\eta ,Re_{i})$ for different values of $P_{N}$ at $\mu =0$ . The curve for $P_{N}=0$ agrees with the result in DiPrima et al. (Reference DiPrima, Eagles and Ng1984). The azimuthal wavenumber of the curves corresponds to $m=0$ as the axisymmetric perturbation is found to be more unstable than non-axisymmetric perturbations. An advantage of using the small- $Pr$ approximation is that, for $Pr\ll 1$ , the stability curves depend solely on a single parameter $P_{N}=N^{2}Pr$ instead of two parameters $N$ and $Pr$ ; thus parametric investigations become simplified. As $P_{N}$ increases, the curves ascend and, in particular, the wide-gap Taylor–Couette flow with small $\eta$ is strongly stabilised. It is also shown that the curves are less sensitive to the change in $P_{N}$ when $\eta$ is close to 1 (i.e. Taylor–Couette flow with a small gap).

This section has described how linear centrifugal instability of stratified Taylor–Couette flow is affected by strong thermal diffusion. It is shown that the stratification having a stabilising effect on the instability is suppressed by strong thermal diffusion. For the cases with $Pr\ll 1$ , the dependence on $N$ and $Pr$ is simplified further by a single rescaled parameter $P_{N}=N^{2}Pr$ as derived by the small- $Pr$ approximation. In § 4, we study how the instability modes develop nonlinearly in stratified and diffusive fluids. Nonlinear dynamical behaviours such as nonlinear saturation, secondary instability of the saturated state or transition to chaotic states are examined.

4.1. Nonlinear saturation to axisymmetric Taylor vortices

Figure 6.

(a) Time evolution of the total energy $E(t)$ (black) and representative modal energy components: $\tilde {E}_{01}$ (blue), $\tilde {E}_{11}$ (green) and $\tilde {E}_{00}$ (red) for case 1. The dashed lines denote 1-D LSA predictions on the growth of the modes with ( $m=0$ ) (blue) and $m=1$ (green). The green dotted line denotes a 2-D bi-global LSA prediction on the decay of the non-axisymmetric mode ( $m=1$ ). (b) Velocity field on the plane $(r,z)$ at $\theta =0$ and $t=500$ with contours denoting the total azimuthal velocity $U_{\theta }$ and vector plot denoting the transverse velocity field $(u_{r},u_{z})$ . (c) The total temperature profile $\varTheta (r,z)$ at $t=500$ .

Figure 6(a) shows an example of the time evolution of the total energy $E(t)$ for case 1 and a few representative examples of the modal energy $\hat {E}_{jl}$ for the axisymmetric mode with $(j,l)=(0,1)$ , the non-axisymmetric mode with $(j,l)=(1,1)$ and the mean-flow distortion with $(j,l)=(0,0)$ . It is recalled that the indices $j$ and $l$ are from the azimuthal wavenumber $m_{j}=jm$ and axial wavenumber $k_{l}=lk$ . At the initial stage, the perturbation energy is small enough and the total energy grows exponentially as $E(t)\sim \exp (2{\textrm {Im}}(\omega _{01})t)$ , where $\omega _{01}$ is the growth rate of the most unstable axisymmetric mode (i.e. $j=0$ and $l=1$ ). A good agreement between the LSA prediction and DNS is clearly shown in figure 6 $(a)$ for the growth of the $m=0$ and $m=1$ modes at early stage. Once the energy of the axisymmetric mode increases and saturates, the flow reaches an equilibrium state featuring axisymmetric Taylor vortices as shown in figure 6( $b$ ). In this saturation process, the mean-flow distortion also grows and its energy $\hat {E}_{00}$ saturates. While the energies of the axisymmetric mode and mean-flow distortion remain constant at large $t$ , the energy of the non-axisymmetric mode with $(j,l)=(1,1)$ decays as the Taylor vortices appear. It is not shown here but other modes with higher $m$ or $k$ also decay as the Taylor vortices appear. Using the Taylor vortices at saturation as a new base state, the bi-global LSA can predict the decay of non-axisymmetric modes with the index $l=1$ as we see a good agreement on the decay rate between the DNS and the bi-global LSA predictions. In figure 6( $c$ ), the total temperature $\varTheta (r,z)=T_{B}+T$ at the equilibrium state is displayed. For case 1, thermal diffusion is moderate with $Pr=1$ and thus the temperature perturbation $T$ can grow and have an amplitude comparable to $T_{B}$ at equilibrium. In this case, the total temperature $\varTheta$ no longer varies linearly but is affected and mixed by the Taylor vortices; i.e. the temperature $\varTheta$ increases (decreases) in the region where the direction of the velocity is upwards (downwards). The Taylor vortices in stratified fluids lead to baroclinicity with a radial temperature gradient; however, for case 1 with moderate $N=1$ and $Re_{i}=145$ , the overturning of the temperature does not occur and it does not become secondarily unstable. Supplementary movie 1 available at https://doi.org/10.1017/jfm.2025.261 demonstrates the nonlinear development of the Taylor vortices and perturbation through the saturation process. It is not shown here but with the same physical parameters of case 1, different sets of numerical parameters are tested (i.e. a higher-resolution case with $(N_{r},N_{\theta },N_{z})=(120,65,65)$ and another case with a longer domain $L_{z}=8\pi /k$ with resolution $(N_{r},N_{\theta },N_{z})=(60,33,129)$ ). The results are found to be the same qualitatively with the appearance of axisymmetric Taylor vortices without secondary instability and quantitatively with insignificant differences in the energy growth curves.

4.2. Bi-global linear stability analysis and secondary instability

Figure 7.

Growth-rate curves from 1-D LSA (dashed lines) and 2-D bi-global LSA (solid lines) for $(a)$ unstratified case $N=0$ with numbers denoting $m$ , $(b)$ $(N,Pr)=(1,1)$ , $(c)$ $(N,Pr)=(1,0.01)$ and $(d)$ $(N,Pr)=(1.5,0.01)$ . Black dashed lines denote the growth rate of the axisymmetric mode $m=0$ and other coloured curves in descending order (for $\mathcal {R}_{c}\lt 1$ ) denote the growth rate for $m=1$ and higher $m$ .

For the Taylor vortex flow, which is axisymmetric, we can explore its secondary instability by performing a 2-D bi-global LSA. To obtain the axisymmetric base state $\bar {\boldsymbol{Q}}(r,z)=(\bar {\boldsymbol{U}}(r,z),\bar {T}(r,z))$ for the bi-global LSA, 2-D DNS with $N_{\theta }=1$ (i.e. $M=0$ ) are conducted instead of 3-D DNS for computational efficiency. Another reason for using 2-D DNS is that the base state $\bar {\boldsymbol{Q}}$ does not become secondarily unstable but remains saturated, and thus we can use this steady and axially periodic base state in the bi-global LSA (see also Park et al. Reference Park, Hwang and Cossu2011). Figure 7 shows examples of the growth rates of the most unstable modes versus the Reynolds-number ratio $\mathcal {R}_{c}=Re_{i}/Re_{i,c}$ . Results are computed from the 1-D local LSA (dashed lines) and 2-D bi-global LSA (solid lines). Highly non-axisymmetric modes not appearing in figure 7 are more stable than the modes shown in figure 7. For every 1-D case in figure 7, the axisymmetric mode with $m=0$ becomes primarily unstable, as shown by black dashed lines. In both 1-D and 2-D LSA, the axial wavenumber $k=k_{c}$ of the axisymmetric mode at the onset of instability $\mathcal {R}_{c}=1$ . In the bi-global LSA, we use the Taylor vortices with axial wavenumber $k=k_{c}$ , the value at the onset of instability $Re_{i}=Re_{i,c}$ , and compute growth rates of non-axisymmetric modes with $m\gt 0$ by increasing $\mathcal {R}_{c}$ . For $\mathcal {R}_{c}\lt 1$ , growth rates of non-axisymmetric modes computed from the bi-global LSA are the same as those from the 1-D LSA since the axisymmetric Taylor vortices are not developed and the base state obtained from 2-D DNS is essentially cylindrical Couette flow (2.4). For the unstratified case with $N=0$ in figure 7(a), the axisymmetric mode becomes unstable at $Re_{i,c}\simeq 131.6$ . As the axisymmetric Taylor vortices appear for $Re_{i}\gt Re_{i,c}$ , the growth rates of non-axisymmetric modes for $m\geqslant 1$ are attenuated by this new base state as shown by coloured solid lines in figure 7(a). In the presence of the Taylor vortices, growth rates of weakly non-axisymmetric modes for $1\leqslant m\leqslant 4$ increase slowly with $\mathcal {R}_{c}$ and the non-axisymmetric mode with $m=1$ becomes unstable at $\mathcal {R}_{c}=1.075$ (i.e. $Re_{i}=Re_{i,2}\simeq 141.5$ , where $Re_{i,2}$ denotes the secondary critical Reynolds number at which a non-axisymmetric mode becomes secondarily unstable). When $\mathcal {R}_{c}=1.075$ , the growth rates of highly non-axisymmetric modes increase faster with $\mathcal {R}_{c}$ for the base flow case with Taylor vortices than the case with cylindrical Couette flow. This implies that highly non-axisymmetric modes are strongly promoted by the axisymmetric Taylor vortices for large $\mathcal {R}_{c}$ . However, the $m=1$ mode is the second unstable mode and it is difficult to predict in advance which non-axisymmetric mode becomes the next dominant mode for $Re_{i}\gt Re_{i,2}$ as nonlinear interactions involving the growth of the $m=1$ mode will lead to a new non-axisymmetric base state.

Characteristics of secondary instability are similar for a stratified case with $(N,Pr)=(1,1)$ in figure 7(b), the case where the axisymmetric mode becomes primarily unstable at a higher $Re_{i,c}\simeq 140.10$ and the non-axisymmetric mode with $m=1$ becomes secondarily unstable at a higher Reynolds-number ratio $\mathcal {R}_{c}=1.085$ (i.e. $Re_{i,2}\simeq 152.1$ ) than the unstratified case. For a highly diffusive case with $Pr=0.01$ in figure 7(c), characteristics of secondary instability change as weakly non-axisymmetric modes with $m=1$ and $2$ are stabilised by the Taylor vortices while a highly non-axisymmetric mode with $m=7$ becomes secondarily unstable at $\mathcal {R}_{c}=1.3$ (i.e. $Re_{i,2}=171.31$ ). This implies that the onset of secondary instability is delayed by strong thermal diffusion at $Pr=0.01$ and highly non-axisymmetric modes become dominant while weakly non-axisymmetric modes are suppressed. The onset of secondary instability is further delayed as stratification becomes stronger with $N=1.5$ as shown in figure 7(d), the case where the Taylor vortices become unstable by the $m=9$ mode at a higher ratio with $\mathcal {R}_{c}=1.466$ (i.e. $Re_{i,2}=193.25$ ).

The prominence of highly non-axisymmetric modes, which become dominant over weakly non-axisymmetric modes with low $m$ , can be understood by conducting an energetics analysis. Similar to (3.9), one can obtain an equation for the temporal evolution of perturbation energy from the linearised perturbation equations (A11)–(A15) as

(4.2) \begin{align} \frac {\partial \bar {E}}{\partial t}=\mathcal {P}_{\bar {U}}+\mathcal {P}_{\bar {V}}+\mathcal {P}_{\bar {W}}+\mathcal {P}_{\bar {\mathcal {T}}}-\mathcal {D}_{k}-\mathcal {D}_{p}, \end{align}

where $\mathcal {P}_{\bar {U},\bar {V},\bar {W},\bar {\mathcal {T}}}$ are the production terms induced by axisymmetric Taylor vortices with $(\bar {U},\bar {V},\bar {W},\bar {\mathcal {T}})$ and $\mathcal {D}_{K,P}$ are the viscous and thermal dissipation terms, respectively, all of which are defined as follows:

(4.3) \begin{eqnarray} &&\mathcal {P}_{\bar {U}}=-\left \lt \frac {\partial \bar {U}}{\partial r}\bar {u}^{2}+\frac {\bar {U}}{r}\bar {v}^{2}+\frac {\partial \bar {U}}{\partial z}\bar {u}\bar {w}\right \gt ,\quad \mathcal {P}_{\bar {V}}=-\left \lt \left (\frac {\partial \bar {V}}{\partial r}-\frac {\bar {V}}{r}\right )\bar {u}\bar {v}+\frac {\partial \bar {V}}{\partial z}\bar {v}\bar {w}\right \gt ,\nonumber \\ && \mathcal {P}_{\bar {W}}=-\left \lt \frac {\partial \bar {W}}{\partial r}\bar {u}\bar {w}+\frac {\partial \bar {W}}{\partial z}\bar {w}^{2}\right \gt ,\quad \mathcal {P}_{\bar {\mathcal {T}}}=-N^{2}\left \lt \frac {\partial \bar {\mathcal {T}}}{\partial r}\bar {u}\bar {T}+\frac {\partial \bar {\mathcal {T}}}{\partial z}\bar {w}\bar {T}\right \gt ,\nonumber \\ &&\mathcal {D}_{K}=\frac {1}{Re}\left \lt {\nabla }\bar {\textbf {u}}\,\textbf {:}\,{\nabla }\bar {\textbf {u}}\right \gt ,\quad \mathcal {D}_{P}=\frac {N^{2}}{RePr}\left \lt {\nabla }\bar {T}\,\textbf {:}\,{\nabla }\bar {T}\right \gt . \end{eqnarray}

We note that the dissipation terms $\mathcal {D}_{K,P}$ are always positive and play a stabilising role while the production terms $\mathcal {P}_{\bar {U},\bar {V},\bar {W},\bar {\mathcal {T}}}$ are not necessarily positive as they depend on the correlation between the base state $(\bar {U},\bar {V},\bar {W},\bar {\mathcal {T}})$ and perturbation variables $(\bar {u},\bar {v},\bar {w},\bar {T})$ . Applying the normal mode (2.29) to (4.2) and considering the normalisation $\hat {E}_{m}=1$ , where $\hat {E}_{m}$ is the modal energy defined as

(4.4) \begin{align} \hat {E}_{m}=\left \lt |\hat {u}_{m}|^{2}+|\hat {b}_{m}|^{2}+|\hat {w}_{m}|^{2}+N^{2}|\hat {T}_{m}|^{2}\right \gt _{rz},\,\left \lt \hat {\textbf {X}}_{m}\right \gt _{rz}=\int _{0}^{L_{z}}\int _{1}^{1/\eta }\hat {\textbf {X}}_{m} \,r\,\mathrm {d}r\,\mathrm {d}z, \end{align}

we obtain the following relation between the growth rate $\omega _{m,i}$ and the modal contribution terms:

(4.5) \begin{align} \omega _{m,i}=\hat {\mathcal {P}}_{m,\bar {U}}+\hat {\mathcal {P}}_{m,\bar {V}}+\hat {\mathcal {P}}_{m,\bar {W}}+\hat {\mathcal {P}}_{m,\bar {\mathcal {T}}}-\hat {\mathcal {D}}_{m,K}-\hat {\mathcal {D}}_{m,P}, \end{align}

where

(4.6) \begin{eqnarray} &&\hat {\mathcal {P}}_{m,\bar {U}}=-\left \lt \frac {\partial \bar {U}}{\partial r}|\hat {u}_{m}|^{2}+\frac {\bar {U}}{r}|\hat {v}_{m}|^{2}+\frac {\partial \bar {U}}{\partial z}\left (\frac {\hat {u}_{m}^{*}\hat {w}_{m}+\hat {u}_{m}\hat {w}^{*}_{m}}{2}\right )\right \gt _{rz},\nonumber \\ &&\hat {\mathcal {P}}_{m,\bar {V}}=-\left \lt \left (\frac {\partial \bar {V}}{\partial r}-\frac {\bar {V}}{r}\right )\left (\frac {\hat {u}_{m}^{*}\hat {v}_{m}+\hat {u}_{m}\hat {v}^{*}_{m}}{2}\right )+\frac {\partial \bar {V}}{\partial z}\left (\frac {\hat {v}_{m}^{*}\hat {w}_{m}+\hat {v}_{m}\hat {w}^{*}_{m}}{2}\right )\right \gt _{rz},\nonumber \\ && \hat {\mathcal {P}}_{m,\bar {W}}=-\left \lt \frac {\partial \bar {W}}{\partial r}\left (\frac {\hat {u}_{m}^{*}\hat {w}_{m}+\hat {u}_{m}\hat {w}^{*}_{m}}{2}\right )+\frac {\partial \bar {W}}{\partial z}|\hat {w}_{m}|^{2}\right \gt _{rz},\nonumber \\ && \hat {\mathcal {P}}_{m,\bar {\mathcal {T}}}=-N^{2}\left \lt \frac {\partial \bar {\mathcal {T}}}{\partial r}\left (\frac {\hat {u}_{m}^{*}\hat {T}_{m}+\hat {u}_{m}\hat {T}^{*}_{m}}{2}\right )+\frac {\partial \bar {\mathcal {T}}}{\partial z}\left (\frac {\hat {w}_{m}^{*}\hat {T}_{m}+\hat {w}_{m}\hat {T}^{*}_{m}}{2}\right )\right \gt _{rz},\nonumber \\ &&\hat {\mathcal {D}}_{m,K}=\frac {1}{Re}\left \lt \hat {\nabla }_{m}\hat {\textbf {u}}_{m}\,\textbf {:}\,\hat {\nabla }_{m}\hat {\textbf {u}}_{m}\right \gt _{rz},\quad \hat {\mathcal {D}}_{m,P}=\frac {N^{2}}{RePr}\left \lt \hat {\nabla }_{m}\hat {T}_{m}\,\textbf {:}\,\hat {\nabla }_{m}\hat {T}_{m}\right \gt _{rz}. \end{eqnarray}

Examples of the growth rate $\omega _{i}$ and instability contribution terms (4.6) computed from the eigenfunction $\hat {\textbf {q}}_{m}$ are presented in table 3. In the table, the growth rate $\omega _{m,i}$ computed from the eigenvalue problems (2.27) and (2.30) and the sum on the right-hand side of (4.5), which is obtained by integrating the eigenfunction computed from either (2.27) for 1-D LSA cases or (2.30) for 2-D LSA cases, are in good agreement with a very small difference of $O(10^{-4})$ . For 1-D LSA cases in which the base flow is cylindrical Couette flow, the axisymmetric mode with $m=0$ is most unstable due to the largest production $\hat {\mathcal {P}}_{m,\bar {V}}$ and the least total dissipation $\hat {\mathcal {D}}_{m}=\hat {\mathcal {D}}_{m,K}+\hat {\mathcal {D}}_{m,P}$ . This applies to both cases with $Pr=0.01$ and $Pr=1$ although the thermal dissipation $\hat {\mathcal {D}}_{m,P}$ is small when $Pr$ is small. As $m$ increases, the production $\hat {\mathcal {P}}_{m,\bar {V}}$ decreases while the dissipation $\hat {\mathcal {D}}_{m}$ increases, which leads to the decrease of the growth rate. For $m=1$ , we compare the 1-D LSA cases against the 2-D LSA cases where the base flow is 2-D axisymmetric Taylor vortices. It is clearly shown that the growth rate decreases as the flow becomes 2-D with decrease in the total production term $\hat {\mathcal {P}}_{m}=\hat {\mathcal {P}}_{m,\bar {U}}+\hat {\mathcal {P}}_{m,\bar {V}}+\hat {\mathcal {P}}_{m,\bar {W}}+\hat {\mathcal {P}}_{m,\bar {T}}$ and increase in the total dissipation $\hat {\mathcal {D}}_{m}$ . While other production terms such as $\hat {\mathcal {P}}_{m,\bar {U}}$ , $\hat {\mathcal {P}}_{m,\bar {W}}$ and $\hat {\mathcal {P}}_{m,\bar {T}}$ appear for 2-D cases due to the velocity and temperature gradients of the Taylor vortices in the radial and axial directions (see also figure 6 b,c), their contribution to the growth rate is smaller than the contribution from the azimuthal velocity $\bar {V}$ (i.e. $\hat {\mathcal {P}}_{m,\bar {V}}$ ). What is noteworthy in secondary instability of axisymmetric Taylor vortices is that the total dissipation $\hat {\mathcal {D}}_{m}$ increases monotonically with $m$ while the total production $\hat {\mathcal {P}}_{m}$ increases and then decreases as $m$ increases.

Table 3.

Growth rates $\omega _{m,i}$ and contribution terms computed from 1-D local LSA and 2-D bi-global LSA for $(Re_{i},N,Pr)=(200,1,0.01)$ and $(Re_{i},N,Pr)=(200,1,1)$ .

Figure 8.

(a,b) Growth rate $\omega _{m,i}$ (filled circles in a), production $\hat {\mathcal {P}}_{m}$ (open circles in b) and dissipation $\hat {\mathcal {D}}_{m}$ (crosses in b) versus azimuthal wavenumber $m$ for different $Pr$ : $Pr=10^{-6}$ (black), $10^{-4}$ (blue), $0.01$ (red) and $1$ (green) at $(Re_{i},N)=(200,1)$ . (c) Production and dissipation terms versus $m$ for $(Re_{i},N,Pr)=(200,1,0.01)$ .

Figure 8(a,b) demonstrates behaviours of the growth rate $\omega _{m,i}$ , production $\hat {\mathcal {P}}_{m}$ and dissipation $\hat {\mathcal {D}}_{m}$ with varying $m$ for different $Pr$ at $(Re_{i},N)=(200,1)$ . Except for $Pr=0.01$ , the growth rate $\omega _{m,i}$ increases first and then decreases as $m$ increases (figure 8 a). This is due to the fact that the production $\hat {\mathcal {P}}_{m}$ increases only for small $m$ and decreases as $m$ increases further while the dissipation $\hat {\mathcal {D}}_{m}$ increases exponentially with $m$ (figure 8 b). This results in the growth rate $\omega _{m,i}$ being positive only in a finite range of $m$ . The growth rates for $Pr=10^{-4}$ and $Pr=10^{-6}$ are very similar to each other and are higher than $\omega _{m,i}$ for $Pr=1$ . The growth rate behaviour for $Pr=0.01$ is more complicated as $\omega _{m,i}$ is negative for $m\leqslant 5$ before it shows a similar increasing/decreasing trend. Only highly non-axisymmetric modes in the range $6\leqslant m\leqslant 12$ are secondarily unstable for $Pr=0.01$ while non-axisymmetric modes for other $Pr$ cases are unstable in broader ranges of $m$ such as $1\leqslant m\leqslant 10$ for $Pr=1$ and $1\leqslant m\leqslant 13$ for $Pr=10^{-4}$ and $10^{-6}$ . This implies that a highly non-axisymmetric flow pattern, as further discussed in the next subsection, can appear for $Pr=0.01$ by secondary instability of Taylor vortices. Figure 8(c) describes variations of all the production and dissipation terms with $m$ at $(Re_{i},N,Pr)=(200,1,0.01)$ . The dominant production and dissipation terms are $\hat {\mathcal {P}}_{\bar {V}}$ and $\hat {\mathcal {D}}_{K}$ while other terms are small. The production term $\hat {\mathcal {P}}_{\bar {W}}$ becomes negative as $m\gt 1$ increases implying the stabilising role of the axial velocity $\bar {W}$ .

Figure 9 details how the growth rate, production and dissipation vary with $Pr$ for various non-axisymmetric modes at $Re_{i}=200$ and $N=1$ . In figure 9 $(a)$ , we see that, for low $m\leqslant 7$ , the growth rate $\omega _{m,i}$ decreases as $Pr$ decreases from $Pr=1$ and then it increases as $Pr$ further decreases from $Pr= 10^{-1}$ to $10^{-2}$ . For high $m\geqslant 8$ , the growth rate increases monotonically as $Pr$ decreases. Such behaviours are similar for the total production term $\hat {P}_{m}$ while the dissipation term $\hat {D}_{m}$ decreases overall monotonically as $Pr$ decreases, as shown in figure 9 $(b)$ . As revealed in figure 8 $(c)$ , the dominant contributions to the production and dissipation come from $\hat {P}_{m,\bar {V}}$ and $\hat {D}_{m,K}$ , respectively, which have the magnitude of $O(1)$ (figure 9 $c$ ). The $Pr$ behaviour of $\hat {P}_{m,\bar {V}}$ is similar to that of the total production $\hat {P}_{m}$ , which drives the growth rate. Although it is not straightforward to delineate the $Pr$ tendency on which non-axisymmetric mode becomes most unstable, one can infer from figure 9 $(c)$ that highly non-axisymmetric modes overall contribute to the secondary instability via its interaction with the azimuthal velocity $\bar {V}$ when $Pr$ is sufficiently low as $Pr\lt 10^{-2}$ , while the behaviour is more complicated in the range $10^{-2}\lt Pr\lt 1$ . In figure 9( $d$ ), we describe other contribution terms $\hat {P}_{m,\bar {U}}$ , $\hat {P}_{m,\bar {W}}$ and $\hat {D}_{m,P}$ , which are minor with a magnitude of $O(0.1)$ or less. The production $\hat {P}_{m,\bar {U}}$ is overall positive except for $m=1$ around $Pr=10^{-2}$ and behaves similar to $\hat {P}_{m,\bar {V}}$ while the contribution $\hat {P}_{m,\bar {W}}$ is overall negative except for the case with $m=1$ , which is positive, and overall decreases as $Pr$ decreases. As can be inferred from the expression in (4.6), thermal dissipation $\hat {D}_{m,P}$ overall decreases to zero as $Pr$ decreases.

Figure 9.

Variations with the Prandtl number $Pr$ for $(a)$ the growth rate $\omega _{m,i}$ , $(b)$ the total production $\hat {P}_{m}$ and dissipation $\hat {D}_{m}$ , $(c)$ the dominant production term $\hat {P}_{m,\bar {V}}$ and dissipation term $\hat {D}_{m,K}$ and $(d)$ other production terms $\hat {P}_{m,\bar {U}}$ , $\hat {P}_{m,\bar {W}}$ and dissipation term $\hat {D}_{m,P}$ for different $m=1,4,7,11$ and 14 at $Re_{i}=200$ and $N=1$ .

Figure 10.

Neutral stability curves from the 1-D LSA of the axisymmetric mode $m=0$ (solid lines denoting $Re_{i,c}$ ) and bi-global LSA of the non-axisymmetric modes (dashed lines denoting $Re_{i,2}$ ) for (a) $N=1$ and (b) $Pr=0.01$ . The numbers above the dashed lines indicate the azimuthal wavenumbers of the non-axisymmetric mode which becomes secondarily unstable.

Figure 10 displays different neutral stability curves obtained from the 1-D local LSA (solid lines) and 2-D bi-global LSA (dashed lines) to show how the critical Reynolds numbers $Re_{i,c}$ and $Re_{i,2}$ vary with the Prandtl number $Pr$ or the Brunt–Väisälä frequency $N$ . For $(N,Pr)=(1,1)$ in figure 10( $a$ ), primary centrifugal instability by an axisymmetric mode occurs at $Re_{i,c}=140.2$ while secondary instability occurs at $Re_{i,2}=152.0$ by a non-axisymmetric mode with $m=1$ . At $N=1$ , the critical Reynolds number $Re_{i,c}$ decreases monotonically as $Pr$ decreases while the secondary critical Reynolds number $Re_{i,2}$ increases as $Pr$ decreases from $Pr=1$ to $Pr=0.06$ and then $Re_{i,2}$ decreases monotonically from the peak at $Pr=0.06$ as $Pr$ further decreases. It is noteworthy that, in the range $10^{-3}\lt Pr\lt 1$ , the secondary instability is triggered by highly non-axisymmetric modes with $m\geqslant 4$ ; thus for $Re_{i}\gt Re_{i,2}$ , we expect to observe highly non-axisymmetric flow patterns in this range of $Pr$ . This highly non-axisymmetric pattern is further discussed in the following subsection. In the range $Pr\leqslant 10^{-3}$ , the critical Reynolds numbers $Re_{i,c}$ and $Re_{i,2}$ approach those of the unstratified case at $N=0$ (i.e. $Re_{i,c}=131.6$ and $Re_{i,2}=141.5$ ), the case with secondary instability triggered by a non-axisymmetric mode with $m=1$ . At $Pr=0.01$ , it is shown in figure 10 $(b)$ that the critical Reynolds number $Re_{i,c}$ does not change significantly in the range $0\leqslant N\leqslant 2$ (e.g. $Re_{i,c}=131.6$ at $N=0$ and $Re_{i,c}=132.0$ at $N=2$ ) while the secondary critical Reynolds number $Re_{i,2}$ increases monotonically with $N$ . The corresponding azimuthal wavenumber of the non-axisymmetric mode, which triggers secondary instability, also increases with $N$ .

Neutral stability curves in figure 10 delineate the regimes of axisymmetric and non-axisymmetric flow patterns in the parameter space $(Re_{i},Pr)$ or $(Re_{i},N)$ . The results are obtained by 1-D and 2-D LSA and we discuss in the next subsection about nonlinear simulation results how flow patterns change from axisymmetric Taylor vortices to non-axisymmetric wavy vortices via secondary instability or how the flow transitions to a chaotic and irregular state, a precursor stage prior to turbulence.

4.3. Transition of flow patterns

Figure 11.

(a) Temporal evolution of the total energy $E(t)$ (black) and modal energy components $\tilde {E}_{jl}$ for case 2: $\tilde {E}_{00}$ (red), $\tilde {E}_{01}$ (blue), $\tilde {E}_{11}$ (green), $\tilde {E}_{91}$ (yellow), and other energy components denoted by grey lines. (b) The corresponding instantaneous velocity profiles at different times on the plane $(r,z)$ at $\theta =0$ . The contours denote the azimuthal velocity $U_{\theta }(r,z)$ and the vector plot denotes the transverse velocity field $(U_{r},U_{z})$ .

Figure 12.

Time evolution of the vertical velocity $u_{z}$ in the $(x,y)$ plane at $z=0$ for case 2.

We now consider case 2 at $(Re_{i},N,Pr)=(200,1,0.01)$ , a case where the axisymmetric Taylor vortices become secondarily unstable by highly non-axisymmetric modes (see also figures 7 $c$ and 10 $a$ ). The perturbation energy plot in figure 11 $(a)$ shows that the axisymmetric mode becomes unstable and saturates quickly to the state of axisymmetric Taylor vortices (e.g. figure 11 $b$ at $t=100$ ). At this first saturation state, the energies of weakly non-axisymmetric modes like the $m=1$ mode (green line in figure 11 $a$ ) decay exponentially. In contrast, highly non-axisymmetric modes like the $m=9$ mode (yellow line in figure 11 $a$ ) become unstable as the axisymmetric Taylor vortices sustain and their amplitudes become comparable to that of the axisymmetric mode as $t\gt 200$ . The flow reaches a second saturation state in which the Taylor vortices start to oscillate vertically (figure 11 $b$ and supplementary movie 2) and become non-axisymmetric. The transition at which non-axisymmetric modes become dominant is clearly illustrated by the vertical velocity $u_{z}(x,y)$ plot for different time $t$ (figure 12 and supplementary movie 3). The most energetic non-axisymmetric mode with $m=1$ loses its energy at the first transition stage (e.g. $u_{z}$ at $t=100$ ) while other non-axisymmetric modes with higher $m$ gain their energies and become comparable in the transition stage $150\lt t\lt 200$ . Weakly and highly non-axisymmetric modes compete during the transition (e.g. $u_{z}$ at $t=150$ and $190$ ) and after $t\gt 200$ , the non-axisymmetric mode with $m=9$ becomes the second dominant mode after the axisymmetric mode $m=0$ forming the wavy Taylor vortices. This non-axisymmetric pattern in $u_{z}$ has a maximum amplitude around the centreline between the two cylinders and is different from the feature of SRI, which is triggered by two out-of-phase inertia–gravity waves confined near the two cylinders and has maxima near the cylinders (Park & Billant Reference Park and Billant2013; Park et al. Reference Park, Billant and Baik2017). Contours of the azimuthal vorticity $\omega _{\theta }=\partial u_{r}/\partial z-\partial u_{z}/\partial r$ in figure 13 $(a)$ clearly show the axisymmetric Taylor vortices at $t=100$ as the first saturation state and non-axisymmetric wavy Taylor vortices as the second saturation state. Various modal energies of the perturbation are expected to either saturate or decay after $t\gt 500$ and thus non-axisymmetric Taylor vortices sustain for large $t$ .

Figure 13.

$(a)$ Contours of the azimuthal vorticity $\omega _{\theta }$ at the boundary surfaces for case 2. $(b)$ Profiles of the averaged velocity $\bar {U}_{\theta }(r)$ and cylindrical Couette flow $V_{B}(r)$ for case 2. $(c)$ Nusselt number $Nu$ as a function of time $t$ for case 2 (black) and case 3 (blue).

Figure 13 $(b)$ shows profiles of the total averaged azimuthal velocity $\bar {U}_{\theta }(r,t)=V_{B}(r)+\tilde {v}_{00}(r,t)$ , the latter which is the azimuthal velocity of the mean-flow distortion denoting the azimuthal perturbation velocity averaged in the directions $\theta$ and $z$ . Compared with the profile of $\bar {U}_{\theta }$ at $t=100$ for the axisymmetric Taylor vortices at the first saturation, the distortion in the azimuthal velocity is reduced for non-axisymmetric wavy Taylor vortices at the second saturation. The reduced distortion in mean flow implies the reduction in the velocity gradient at the inner cylinder and the corresponding torque applied at the inner cylinder. Following Martínez Arias (Reference Martínez Arias2015), we define an axially averaged non-dimensional torque $G$ at the inner cylinder as

(4.7) \begin{align} G=\left (\frac {\eta Re_{i}}{1-\eta }\right )\frac {1}{2\pi L_{z}}\int _{0}^{2\pi }\int _{0}^{L_{z}}\left .\left (\frac {U_{\theta }}{r}-\frac {\partial U_{\theta }}{\partial r}\right )\right |_{r=1}\mathrm {d}\theta \,\mathrm {d}z,\, \end{align}

which gives the laminar torque $G_{\mathrm {lam}}$ for cylindrical Couette flow $U_{\theta }=V_{B}(r)$ as

(4.8) \begin{align} G_{\mathrm {lam}}=\frac {2\eta (1-\mu )Re_{i}}{(1+\eta )(1-\eta )^{2}}. \end{align}

We also define the Nusselt number $Nu$ as $Nu=G/G_{\mathrm {lam}}$ , which is the ratio between the transverse convective transport of angular velocity and the molecular transport of angular velocity as a measure of the non-dimensional angular momentum transfer (Dubrulle & Hersant Reference Dubrulle and Hersant2002; Eckhardt et al. Reference Eckhardt, Grossmann and Lohse2007; Martínez Arias Reference Martínez Arias2015). Figure 13 $(c)$ shows how the Nusselt number $Nu$ changes over time for case 2 (black line). Compared with the axisymmetric Taylor vortices in the range $50\leqslant t\leqslant 150$ having a constant $Nu=1.77$ during the first saturation process, the wavy Taylor vortices at the second saturation in the range $t\gt 200$ have a lower $Nu=1.68$ . This implies that the non-axisymmetric Taylor vortices after secondary instability lead to reduced angular momentum transfer compared with the axisymmetric Taylor vortices after primary centrifugal instability.

Figure 14.

(a) Temporal evolution of the total energy $E(t)$ (black) and modal energy components $\tilde {E}_{jl}$ for case 3: $\tilde {E}_{00}$ (red), $\tilde {E}_{01}$ (blue), $\tilde {E}_{11}$ (green), $\tilde {E}_{41}$ (brown), and other energy components denoted by grey lines. Corresponding instantaneous velocity field $u_{z}(x,y)$ at $z=0$ at (b) $t=200$ and (c) $t=500$ .

Figure 15.

(a) Evolution of the total temperature $\varTheta (r,z)$ at $\theta =0$ for case 3. (b) Corresponding spatio-temporal diagram of $\varTheta (t,z)$ measured at $r=1.05$ and $\theta =0$ .

Case 3 with $(N,Pr)=(1,1)$ considers the Reynolds number $Re_{i}=200$ at the ratio $\mathcal {R}_{2}=Re_{i}/Re_{i,2}=1.316$ , which is larger than the ratio $\mathcal {R}_{2}=1.168$ in case 2. In this case, it is observed that the secondary instability bypasses the saturation but leads to an irregular fluctuation as shown by the temporal variation of the Nusselt number $Nu$ in figure 13 $(c)$ and modal energies of perturbation in figure 14 $(a)$ . Primary centrifugal instability breaks down as $t\gt 100$ where other non-axisymmetric modes including the mode with $m=1$ become dominant with strong nonlinear modal interaction. For case 3, the most energetic mode varies over time and the flow exhibits a rather chaotic temporal fluctuation as shown by the energy evolution in figure 14 $(a)$ and by the velocity and temperature fields in the $(r,z)$ plane in supplementary movies 4 and 5. Although this irregular flow state does not show a small length-scale structure (see also figure 14 $b$ and supplementary movie 6), the levels of the modal energies of large-wavenumber modes are not negligible and turbulence is expected to occur for $(N,Pr)=(1,1)$ at high Reynolds number $Re_{i}\gt 200$ , the range to be further explored in a future study. The Prandtl number $Pr=1$ is not small and the temperature perturbation $T$ is comparable with the base temperature $T_{B}$ . Therefore, the fluctuation in the total temperature $\varTheta =T_{B}+T$ induced by the velocity fluctuation is visible as shown in figure 15 $(a)$ although the overturning of temperature, which is a potential source of baroclinic instability, does not occur. Chaotic fluctuation of the temperature $\varTheta$ after $t\gt 100$ is also clearly shown in the spatio-temporal diagram measured near the mid-point at $r=1.05$ (figure 15 $b$ ).

Figure 16.

$(a)$ Nusselt number $Nu$ versus Reynolds number $Re_{i}$ for $N=0$ (black), $(N,Pr)=(1,1)$ (blue), $(N,Pr)=(1,0.01)$ (red) and $(N,Pr)=(1,10^{-4})$ (green). Filled circles and squares indicate $Nu$ for the axisymmetric Taylor vortices and non-axisymmetric wavy vortices at saturation, respectively. Open squares with error bars indicate statistically averaged $Nu$ and the minimum and maximum of $Nu$ in the time interval considered (see also figure 13 $c$ ). $(b)$ Nusselt number $Nu$ same as $(a)$ but over the rescaled Reynolds-number ratio $\mathcal {R}_{c}$ . The grey line indicates the Nusselt number from (4.9). $(c)$ Temporal evolution of the total energy $E(t)$ for different $Pr$ at $Re_{i}=200$ and $N=1$ . Black line denotes the unstratified case with $N=0$ .

Figure 16 $(a)$ shows how the Nusselt number $Nu$ changes with the Reynolds number $Re_{i}$ for different sets of $(N,Pr)$ : $(1,1)$ (blue), $(1,0.01)$ (red), $(1,10^{-4})$ (green) and the unstratified case with $N=0$ (black). The Nusselt number $Nu$ is computed at the saturation state of the axisymmetric Taylor vortices (filled circles), which are computed through 3-D DNS with a numerical resolution the same as that in case 1, or for non-axisymmetric wavy vortices at the second saturation (filled squares) computed by 3-D DNS with a numerical resolution the same as that in case 2. For the flow in which its instability does not saturate and $Nu$ fluctuates like case 3, open squares with error bars are used to indicate the mean, minimum and maximum $Nu$ averaged in the time interval of fluctuation (see e.g. figure 13 $c$ ). In these cases with fluctuations, the same numerical resolution as in case 3 is used. It is clearly shown that $Nu$ for cases with the axisymmetric Taylor vortices increases fast with $Re_{i}$ and it increases slowly down as secondary instability occurs. The increasing trend of the Nusselt number $Nu$ is similar for unstratified cases with $N=0$ and stratified cases with $(N,Pr)=(1,10^{-4})$ . For both sets, secondary instability leading to a fluctuating flow state appears for $Re_{i}\gt 160$ and there is a jump in $Nu$ around $Re_{i}=180$ . The Nusselt number $Nu$ continues to increase similarly for both as $Re_{i}$ increases further for $Re_{i}\gt 180$ . For the cases with $(N,Pr)=(1,0.01)$ , the Nusselt number $Nu$ is higher than for other cases as secondary instability occurs at a later stage as $Re_{i,2}=171.31$ for $Pr=0.01$ . Up to $Re_{i}=210$ , the breakdown of secondary instability is not observed for $Pr=0.01$ . Unlike other cases with critical Reynolds numbers $Re_{i,c}=131.6{-}131.7$ , the centrifugal instability for $(N,Pr)=(1,1)$ occurs at a higher $Re_{i,c}=140.2$ and thus $Nu$ in this case is smaller than for other cases at the same $Re_{i}$ . Although the onset of primary centrifugal instability is delayed for $(N,Pr)=(1,1)$ , it is found that the breakdown of instability occurs early and the fluctuation in $Nu$ is larger than in other cases.

For a better comparison among the cases, we display in figure 16 $(b)$ the same $Nu$ against the Reynolds-number ratio $\mathcal {R}_{c}=Re_{i}/Re_{i,c}$ . For the axisymmetric Taylor vortices, DiPrima et al. (Reference DiPrima, Eagles and Ng1984) proposed the following scaling law for $Nu$ as a function of the ratio $\mathcal {R}_{c}$ :

(4.9) \begin{align} Nu=1+A\left (1-\frac {1}{\mathcal {R}_{c}^{2}}\right )+B\left (1-\frac {1}{\mathcal {R}_{c}^{2}}\right )^{2}, \end{align}

where $A=1.246$ and $B=-0.37$ are the constants for $\eta =0.9$ according to DiPrima et al. (Reference DiPrima, Eagles and Ng1984). It is clearly shown in figure 16 $(b)$ that $Nu$ for every case agrees well with the scaling law (4.9) in the range $\mathcal {R}_{c}\leqslant 1.1$ where the axisymmetric Taylor vortices are present. This implies that $Nu$ for the Taylor vortices, which are induced by primary centrifugal instability, is not strongly influenced by $Pr$ . When the flow is secondarily unstable or chaotic, the cases with $(N,Pr)=(1,0.01)$ for which secondary instability appears late as $\mathcal {R}_{c}\gt 1.3$ have higher $Nu$ than the scaling law (4.9) while other cases for which secondary instability appears early as $\mathcal {R}_{c}\gt 1.1$ have lower $Nu$ than the proposed scaling (4.9). This implies that the angular momentum transport represented by $Nu$ depends on $Pr$ once the flow becomes secondarily unstable or chaotic. However, as figure 16(a,b) only considers $N=1$ , more detailed investigations are required to confirm the relation between $Nu$ and $Pr$ in a wider parameter space of $(N,Re_{i})$ , an important topic to be further studied in the future to unravel the physics of turbulent angular momentum transport that depends on stratification and thermal diffusion.

Figure 17.

Instantaneous velocity contours for $U_{\theta }$ and vector plots for the transverse velocity field $(U_{r},U_{z})$ on the plane $(r,z)$ for different $Pr$ at $t=500$ , $\theta =0$ , $Re_{i}=200$ and $N=1$ . The rightmost panel corresponds to the unstratified case with $N=0$ .

As a summary of the nonlinear dynamics with varying $Pr$ , we display in figures 16 $(c)$ and 17 the temporal evolution of the total energy $E(t)$ and instantaneous velocity fields for different $Pr$ at $Re_{i}=200$ and $N=1$ . We see in figure 16 $(c)$ that, as $Pr$ is small, $Pr\lt 10^{-4}$ , the energy curves collapse and become identical to the curve of the unstratified case with $N=0$ . Instantaneous velocity fields obtained at $t=500$ from various DNS, which started from the same configuration of the initial condition described in the early part of § 4, are also identical between the stratified cases with $Pr=10^{-5}$ and $10^{-6}$ and unstratified case with $N=0$ . These results support the explanation in the introduction and figure 1 that the stratification effect is suppressed by strong thermal diffusion and the flow behaves as unstratified flow in both linear and nonlinear regimes. Nonetheless, there still remains a question as to whether our conclusion will remain valid for turbulent cases and if so, up to which Reynolds number $Re$ . Although DNS investigations like those by Prat & Lignières (Reference Prat and Lignières2013) verified this question for $Re$ of $O(100)$ , further investigations should be conducted at higher Reynolds numbers.