Transient Dynamics in Counter-Rotating Stratified Taylor–Couette Flow

Abstract

This study focuses on the investigation of stratified Taylor–Couette flow (STCF) using non-modal analysis, which has received relatively limited attention compared to other shear flows. The dynamics of perturbations under different temperature conditions are explored, and their patterns of amplification are analyzed. The study highlights the correlation between flow configurations, emphasizing the similarity in transient dynamics despite different speed ratios. The subcritical effects of thermal stratification on disturbance dynamics are examined, considering the interplay between viscous and buoyancy effects counteracted by strong centrifugal forces. It is found that increasing the wall temperature beyond a critical value leads to buoyancy forces dominating, resulting in a linear increase in the amplification factor. The research reveals significant deviations from previous results, indicating the significant role of temperature stratification.

1. Introduction

The transition from laminar to turbulent fluid flow has been of interest for over a century. Beyond the laminar state, fluid flows can exhibit several states before they reach a turbulent state. These states have played a role in making significant contributions in industries such as automobile, civil, aerospace, and others. To address the instabilities arising from flow disturbances, a multitude of mathematical techniques have been developed within hydrodynamic stability theory [1,2].

Numerical and analytical methods, including bifurcation theory and other suitable approaches, have been employed to identify these flow patterns [3,4]. Among the various techniques used, the transient growth method has gained extensive utilization in the exploration of subcritical transitions in fluid flow, proving to be highly successful over several decades. This method has been invaluable in capturing and elucidating phenomena associated with stability problems in fluid flows that cannot be effectively addressed by linear stability theory alone [5,6,7,8,9].

The plane Couette and pipe Poiseuille flows are known to be stable for all Reynolds numbers within the theoretical mathematical framework of linear stability analysis, but in practice, these flows become turbulent for sufficiently large Reynolds numbers, as accounted for by Romanov [10], Davey [11], Drazin and Reid [12] and recent studies by others [13,14]. In practice, the plane Poiseuille flow does exhibit an instability at a critical Reynolds number Re${}_c\approx$ 2300, but theoretically this phenomenon is observed at a much larger value, Re${}_c\approx$ 5772, at which the so-called Tollmien–Schlichting wave depicting an unstable mode appears [15,16]. The inconsistency in results has been attributed to the non-normality of the linearized Navier–Stokes operator. It is known that the non-normality is responsible for the spike (i.e., transient growth) in the growth rate within the linear regime that eventually triggers the transition from laminar to turbulent flow. This behavior is widely known as a subcritical transition or bypass transition [5,17,18]. The non-normal route to turbulence has raised many concerns in the literature from different schools of thought [7,19,20], conversely, Reddy and Henningson [21] objectively addressed the issues with insightful results for two- and three-dimensional plane Couette and Poiseuille flows. The utilization of the transient growth method for investigating subcritical bifurcations has witnessed widespread adoption and significant advancements over the last two decades, with both theoretical and experimental justifications [8,22,23,24].

The circular Couette flow, often known as Taylor–Couette flow (TCF), has been widely studied experimentally and theoretically [25,26,27,28,29]. The TCF problem has been known to be one of the pioneering problems of hydrodynamic stability for several decades, since the early experimental and theoretical works of Taylor [30]. In 1965, Coles [31] experimentally investigated the transition of viscous flow in concentric cylinders and reported that at certain Reynolds numbers, where the linear theory predicted laminar flow, there exists spiral turbulence in the counter-rotating flow [31]. In the following year, Van Atta [32] carried out similar experiments on the spiral transition of the flow and obtain related results. These results were confirmed experimentally later, in 1989, by Hegseth et al. [33] and later, in 2002, by Prigent et al. [34].

The non-normality problem of the TCF linear operator was initially investigated by Grebhardt and Grossmann in 1992 [26]. However, Hristova et al. [35] conducted a comprehensive investigation of this phenomenon using pseudospectral analysis. In 2002, Meseguer [36] examined the transient energy growth associated with the parametric regime reported by Coles [31] and observed a significant amplification of the transient energy in the counter-rotating regime, which aligned well with Coles’s [31] experimental findings. It was discovered that the non-axisymmetric modes experience greater amplification and they could potentially trigger the subcritical transition to turbulence [19,36]. Maretzke et al. [37] further expanded on this perspective through numerical and analytical methods. The authors explored the transient energy growth rates of the modes across all TCF regimes and established a universal Re${}^{2/3}$ scaling. Additionally, it was established that the linear stability and transient growth rate remained unaffected by the ratio of the cylinders’ rotation rates in the limit where there was no axial perturbation dependency.

However, there has been a significant lack of attention given to the application of non-modal analysis in the context of stratified Taylor–Couette flow (STCF). Thus, the primary objective of this research is to extend the existing experimental work conducted by Cole [31] in the domain of thermal convection by employing non-modal analysis techniques. Building upon the configuration reported by Meseguer [36], we introduce radial heating and proceed with a comprehensive investigation of the subcritical phenomena associated with counter-rotating flows. This article specifically focuses on the optimal transient growth rate observed in counter-rotating stratified Taylor–Couette flow, with a particular emphasis on non-axisymmetric perturbations. To provide a comprehensive understanding of the topic, Section 2 presents a concise overview of the governing equations underlying the Taylor–Couette problem. The derived model is then linearized in Section 3. Section 4 offers an in-depth discussion of the results obtained from our investigation. Finally, the overall conclusions drawn from this study are summarized in Section 5.

2. Mathematical Formulation

2.1. Model

We examine the behavior of an incompressible fluid with kinematic viscosity $\mu$, and a density $\rho ={\rho }_o+{\rho }^{\prime }$, which is contained between two concentric cylinders of infinite height, where ${\rho }_o$ is the average density of the fluid and ${\rho }^{\prime }$ is the variation in the density. The inner cylinder has a radius $r_i$, and rotates at an angular velocity ${\Omega }_i>0$, while the outer cylinder has a radius $r_o$, and rotates at an angular velocity ${\Omega }_o<0$. Additionally, we assume a negative temperature gradient with an inner temperature $T_i=\overline{T}+\Delta T/2$, and an outer temperature $T_o=\overline{T}-\Delta T/2$, as depicted in Figure 1, where $\overline{T}$ is the ambient temperature and $\Delta T$ is the temperature difference between the two cylinders. Thus, the model is governed by the incompressible Navier–Stokes equations for rapidly rotating flows, as described by Lopez et al. [28,38,39]:

$${\rho }_o({\partial }_t\mathbf{v}+(\mathbf{v}·\nabla )\mathbf{v})=-\nabla \Pi +\mu {\nabla }^2\mathbf{v}+{\rho }^{\prime }\nabla \Pi +{\rho }^{\prime }(\mathbf{v}·\nabla )\mathbf{v},$$

(1)

$$\nabla ·\mathbf{v}=0,\mathrm{and}$$

(2)

$${\partial }_{t}T+(\mathbf{v}·\nabla )T=\kappa {\nabla }^{2}T,$$

(3)

where $\kappa$ is the thermal diffusivity of the fluid; $\Pi =p+{\rho }_o\Phi$ is composed of the gradient terms pressure p, and constant gravitational potential ${\rho }_o\Phi$.

The boundary conditions are given as:

$$\mathbf{v}·{\mathbf{e}}_{\theta }={\Omega }_i\mathrm{on}r=r_i\mathrm{and}\mathbf{v}·{\mathbf{e}}_{\theta }={\Omega }_o\mathrm{on}r=r_o.$$

(4)

Equation (1) is defined in the stationary reference frame and the last term, (${\rho }^{\prime }(\mathbf{v}·\nabla )\mathbf{v}$), describes the centrifugal buoyancy contribution of the flow. This term is necessary in order to accurately account for strong vortex dynamics or flows with rapidly rotating walls.

2.2. Boussinesq Approximation

The system is defined by a cylindrical coordinate system $(r,\theta ,z)$. We proceed with the Boussinesq assumption by ignoring density variations in the flow, and we only keep density contributions from the gravitational and centrifugal terms. Furthermore, we assume a constant gravitational acceleration, g, in the vertical direction, z. By assuming that the temperature has a linearly decreasing effect on the density, we have:

$$\rho ={\rho }_o(1-\beta T)$$

(5)

where $\beta =-(1/{\rho }_o){\partial \rho /\partial T|}_{\rho ={\rho }_o}$ is the coefficient of thermal expansion. Since $\rho ={\rho }_o+{\rho }^{\prime }$, we have:

$${\rho }^{\prime }=-{\rho }_o\beta T.$$

(6)

Now, applying (6) to the buoyancy term in (1) and using $\Phi =gz$, we have:

$${\partial }_t\mathbf{v}+(\mathbf{v}·\nabla )\mathbf{v}=-\nabla \Pi +\nu {\nabla }^2\mathbf{v}-\beta gT{\mathbf{e}}_z-\beta T(\mathbf{v}·\nabla )\mathbf{v}.$$

(7)

2.3. Nondimensionalization

In this section, we proceed by converting Equations (2), (3), and (7) to dimensionless forms by using $\Delta T$, the gap-width $d=r_o-r_i$, $d^2/\nu$, and ${(\nu /d)}^2$ as the temperature, length, viscous time, and reduced pressure $p/{\rho }_o$ scaling factors, respectively. The inner $r_i$ and outer $r_o$ radii are expressed in dimensionless form as $\eta /(1-\eta )$ and $1/(1-\eta )$, respectively. With these scaling factors we obtained the new dimensionless parameters listed in

Table 1. Definitions of the dimensionless parameters.

Parameter	Expression
Grashof number	$\mathrm{Gr}=\beta g\Delta T{d}^3/{\nu }^2$
Inner Reynolds number	${\mathrm{Re}}_i={\Omega }_i{r}_{i}d/\nu$
Outer Reynolds number	${\mathrm{Re}}_o={\Omega }_o{r}_{o}d/\nu$
Prandtl number	$\mathrm{Pr}=\nu /\kappa$
Radius ratio	$\eta =r_i/r_o$
Relative density	$ϵ=\beta \Delta T$

. In addition, we introduce a dummy variable ${\rm Y}\in \{0,1\}$ into the equation:

$${\partial }_t\mathbf{v}+(\mathbf{v}·\nabla )\mathbf{v}=-\nabla p+{\nabla }^2\mathbf{v}+{\rm Y}\mathrm{Gr}T{\mathbf{e}}_z+{\rm Y}ϵT(\mathbf{v}·\nabla )\mathbf{v},$$

(8)

$$\nabla ·\mathbf{v}=0,$$

(9)

$${\partial }_{t}T+(\mathbf{v}·\nabla )T=\frac{1}{\mathrm{Pr}}{\nabla }^{2}T,$$

(10)

where ${\mathbf{e}}_z$ is the unit vector in the axial direction z.

If ${\rm Y}$ is 1, then (8) describes the STCF system and if it is 0, then it defines the TCF system without stratification. Thus, we use these parameters and the dimensionless Equations (8)–(10) in subsequent sections for our analysis and discussion.

2.4. Basic Equations

Given the nondimensional equations, let the velocity be expressed as $\mathbf{v}=u{\mathbf{e}}_r+v{\mathbf{e}}_{\theta }+w{\mathbf{e}}_z$ in cylindrical coordinates, where u, v, and w are the radial, azimuthal, and axial components of the velocity, respectively, and ${\mathbf{e}}_r$ and ${\mathbf{e}}_{\theta }$ are the unit vectors in the r, $\theta$, and z directions. The basic velocity $\mathbf{U}$, becomes:

$$\mathbf{U}=U_b{\mathbf{e}}_r+V_b{\mathbf{e}}_{\theta }+W_b{\mathbf{e}}_z$$

(11)

and the basic components of the velocity $U_b$, $V_b$, and $W_b$ correspond to the respective radial, azimuthal, and axial basic profiles. Thus, we proceed by deriving expressions for the basic profiles by assuming that $\mathbf{U}$ is independent of time t, and the azimuthal $\theta$ and axial z directions. In addition, in order to have a constant pressure gradient in the axial direction, we impose the zero mass flux constraint,

$${\mathit{\int }}_{r_i}^{r_o}r{W}_b\left(r\right)dr=0$$

(12)

and thus, with these assumptions we obtain the following analytical solutions:

$$U_b\left(r\right) = 0$$

(13)

$$V_b\left(r\right) = Ar+\frac{B}{r}$$

(14)

$$\begin{gathered} W_b\left(r\right) = \frac{\mathrm{Gr}{\rm Y}}{{(1-\eta )}^2}[C({\left(r(1-\eta )\right)}^2-{\eta }^2) \\ \hspace{1em} \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +\left(C(1-{\eta }^2)+\frac{1}{8}(1-{\left(r(1-\eta )\right)}^2)\right)\left(2T_b\left(r\right)-1\right)] \end{gathered}$$

(15)

$$T_b\left(r\right) = \frac{ln\left(\right(r/\eta \left)\right(1-\eta \left)\right)}{ln\left(\eta \right)}+\frac{1}{2}$$

(16)

where $V_b$, $W_b$, and $T_b$ remain the same solutions introduced in the related literature. The constants A, B and C are given by the expressions:

$$A=\frac{{\mathrm{Re}}_o-\eta {\mathrm{Re}}_i}{1+\eta },B=\frac{\eta ({\mathrm{Re}}_i-\eta {\mathrm{Re}}_o)}{(1-\eta )(1-{\eta }^2)},$$

(17)

and

$$C=\frac{({\eta }^2-3)({\eta }^2-1)+4ln\left(\eta \right)}{16({\eta }^2-1)((1+{\eta }^2)ln\left(\eta \right)+1-{\eta }^2)}.$$

(18)

3. Linearization

Furthermore, we suppose that the basic state is subjected to an infinitesimal perturbation:

$$\mathbf{v}=\mathbf{U}+\delta \tilde{\mathbf{v}},$$

(19)

$$p=P_b+\delta \tilde{p},$$

(20)

$$T=T_b+\delta \tilde{T},$$

(21)

where $\delta$ is assumed to be a very small constant and $\tilde{\mathbf{v}}$, $\tilde{p}$ and $\tilde{T}$ are the residues (fluctuations) of the velocity, pressure and temperature, respectively. We substitute Equations (19)–(21), into the dimensionless governing Equations (8)–(10). We then proceed by collecting the O($\delta$) terms and ignore higher-order terms. Thus, the linearized equations become:

$${\partial }_t\tilde{\mathbf{v}}=-\nabla \tilde{p}+{\nabla }^2\tilde{\mathbf{v}}+{\rm Y}G\tilde{T}{\mathbf{e}}_z+({\rm Y}ϵ\tilde{T}-1)((\mathbf{U}·\nabla )\tilde{\mathbf{v}}+(\tilde{\mathbf{v}}·\nabla )\mathbf{U}),$$

(22)

$$\nabla ·\tilde{\mathbf{v}}=0,$$

(23)

$${\partial }_t\tilde{T}=\frac{1}{\mathrm{Pr}}{\nabla }^2\tilde{T}-(\mathbf{U}·\nabla )\tilde{T}-(\tilde{\mathbf{v}}·\nabla )T_b.$$

(24)

Our objective is to describe the transition of Taylor–Couette flow from laminar to turbulent flow. In particular, we seek to identify any differences between isothermal and stratified Taylor–Couette flow. To facilitate our analysis, and expecting harmonic motion, we utilize the Fourier transform in the azimuthal $\theta$ and axial z directions, while we employ the Gauss–Lobatto collocation point Chebyschev polynomial expansion in the radial direction r; see Figure 1. Based on these considerations, we propose the following functional forms:

$$\begin{array}{ccc}\hfill \tilde{\mathbf{v}}& \approx & \widehat{\mathbf{v}}(r,t)e^{i(n\theta +kz)},\tilde{p}\approx \widehat{p}(r,t)e^{i(n\theta +kz)},\tilde{T}\approx \widehat{T}(r,t)e^{i(n\theta +kz)}\hfill \end{array}$$

where the pressure $\tilde{p}$, temperature $\tilde{T}$, and velocities $\tilde{u}$, $\tilde{v}$, and $\tilde{w}$ are functions of the radial, azimuthal, and axial coordinate directions. Furthermore, $n\in \mathbb{N}$ and $k\in \mathbb{R}$ are the azimuthal and axial wavenumbers, respectively. The frequencies of the disturbance are characterized by the wavenumbers. The wavelengths in the homogeneous directions $\theta$ and z are $L_{\theta }=2\pi /n$ and $L_z=2\pi /k$, respectively. Thus, we proceed by substituting the above functional forms into the linearized Equations (22)–(24), which after further simplification become:

$$\begin{array}{cc}\hfill \frac{\partial \widehat{u}}{\partial t}=& -\frac{\partial \widehat{p}}{\partial r}+\mathcal{F}(n,k,{\rm Y})\widehat{u}-\left[\frac{2in}{r^2}-\frac{2V_b}{r}(1-{\rm Y}ϵT_b)\right]\widehat{v}-\frac{{\rm Y}ϵV_b^2}{r}\widehat{T},\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill \frac{\partial \widehat{v}}{\partial t}=& -\frac{in\widehat{p}}{r}+\mathcal{F}(n,k,{\rm Y})\widehat{v}+\left[\frac{2in}{r^2}-(1-{\rm Y}ϵT_b)\left(\frac{V_b}{r}+\frac{\partial V_b}{\partial r}\right)\right]\widehat{u},\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill \frac{\partial \widehat{w}}{\partial t}=& -ik\widehat{p}+\left[\mathcal{F}(n,k,{\rm Y})+\frac{1}{r^2}\right]\widehat{w}-\frac{\partial W_b}{\partial r}(1-{\rm Y}ϵT_b)\widehat{u}+{\rm Y}\mathrm{Gr}\widehat{T},\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill D_{+}\widehat{u}+& \frac{in}{r}\widehat{v}+ik\widehat{w}=0,\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill \frac{\partial \widehat{T}}{\partial t}=& \frac{1}{\mathrm{Pr}}\left[D_{+}D-\frac{n^2}{r^2}-k^2-i\mathrm{Pr}\left(\frac{n{V}_b}{r}+k{W}_b\right)\right]\widehat{T}-\frac{\partial T_b}{\partial r}\widehat{u},\hfill \end{array}$$

(29)

where

$$\mathcal{F}(n,k,{\rm Y})=D_{+}D-\frac{n^2+1}{r^2}-k^2-i\left(\frac{n{V}_b}{r}+k{W}_b\right)(1-ϵ{\rm Y}{T}_b),$$

$D=d/dr$, and $D_{+}=D+1/r$. Manipulating Equations (25)–(29), we can represent the solenoidal system of equations in a compact form as an initial value problem:

$$\frac{\partial }{\partial t}\mathbf{B}\widehat{\mathbf{q}}=\mathbf{A}\widehat{\mathbf{q}},$$

(30)

where $\mathbf{A}$ and $\mathbf{B}$ are defined in Appendix B. By assuming a solution of the form

$$\widehat{\mathbf{v}}=\breve{\mathbf{v}}\left(r\right)e^{\lambda t},\widehat{p}=\breve{p}\left(r\right)e^{\lambda t},\mathrm{and}\widehat{T}=\breve{T}\left(r\right)e^{\lambda t}.$$

(31)

We transform the initial value problem to a generalized eigenvalue problem:

$$\lambda \widehat{\mathbf{q}}=\mathbf{L}\widehat{\mathbf{q}}\mathrm{and}\mathbf{L}={\mathbf{B}}^{-1}\mathbf{A}.$$

(32)

The eigenvalues are complex and define the temporal linear stability of the flow. That is, if the real parts of all eigenvalues are negative, the flow is linearly stable and the amplitude of the perturbations will decay in time. Furthermore, if there exists an eigenvalue $\lambda$ whose real part is positive, then the flow is unstable and the amplitude of the perturbation will grow in time. We assume that the perturbation of the velocity and temperature of the fluid motion must vanish at the walls:

$$\breve{\mathbf{v}}\left(r_i\right)=\breve{\mathbf{v}}\left(r_o\right)=\breve{T}\left(r_i\right)=\breve{T}\left(r_o\right)=0.$$

(33)

4. Numerical Results

The amplification factor, often denoted as $G_0$, is a fundamental concept in the study of transient growth in fluid dynamics. It quantifies the maximum amplification of perturbations or disturbances in a given flow system over a certain time interval. Transient growth refers to the temporary amplification of perturbations in a flow before they eventually decay or become insignificant. This phenomenon is particularly relevant in the subcritical transition to turbulence, where small disturbances can trigger the onset of turbulent behavior. The amplification factor, $G_0$, represents the maximum energy amplification that can be achieved by the optimal initial perturbations within a specified time frame. It characterizes the efficiency of energy transfer from the mean flow to the perturbations. A larger value of $G_0$ indicates a higher potential for transient growth and the presence of stronger amplification mechanisms in the flow. The computation of $G_0$, as defined in Appendix A, involves solving an eigenvalue problem or performing a non-modal analysis, depending on the nature of the flow and the mathematical framework employed. By determining the optimal initial conditions and analyzing the linearized equations governing the flow, $G_0$ can be evaluated. In addition, the structures and modes that contribute most significantly to the amplification can be identified. In our study, we employed well-established approaches from the literature to compute the amplification factor $G_0$. These approaches have been widely used and validated in previous research, providing reliable and effective methods for evaluating the maximum amplification of disturbances in fluid flows. By utilizing these established techniques, we ensure the accuracy and robustness of our calculations, allowing us to gain valuable insights into the transient growth and stability characteristics of the flow system under investigation [5,40,41].

The configuration is chosen so that our results can compare well with earlier studies of Taylor–Couette flow (TCF) without stratification [35,36,37]. Hereafter, we apply stratification to the problem to investigate the behavior of the transient energy growth of the flow. Thus, we continue with an equivalent configuration used by Meseguer [36] for the study of energy transient growth in the TCF problem. We choose fixed values of $\eta =0.881$, $ϵ=0.067$, and $\mathrm{Pr}=68$, and various values of Gr along with various wavenumber pairs ($n\in [0,15]$, $k\in [0,15]$). In addition, we consider specific counter-rotating pairs of values for Re${}_i$ and Re${}_o$ from the set of the most significant pairs reported by Meseguer [36] for the TCF case. Furthermore, for proper justification of this investigation, we consider the symmetry introduced due to the periodic assumption made for the azimuthal (i.e., the SO(2)-symmetry because of the invariance of the system with respect to the azimuthal rotation about the z-axis) and axial (i.e., the O(2)-symmetry because of the invariance of the system with respect to the reflections and translations of the axial axis) coordinate directions.

Due to the fact that no experimental work exists for the choice of configuration to which we can compare our results, we have introduced a parameter ${\rm Y}$ in the governing dimensionless equation with values chosen to be either 0 or 1. The value of ${\rm Y}$ is set to 0 in order to compare our findings with already existing results obtained for TCF [35,37]. The results, as indicated in Figure 2, are in perfect agreement with the results reported by Hristova et al. [35], Meseguer [36], and recently by works of Maretzke et al. [37]. In addition, we verify that the eigenvalues of the STCF case are in perfect agreement with Lopez [38], and McFadden et al. [38].

The nondimensional scheme we used is similar to Meseguer [36] and Maretzke et al. [37], but is slightly different from the scales used by Hristova et al. [35]. Thus, we made the conversions ${\mathrm{Re}}_i^{\mathrm{alt}}=2{\mathrm{Re}}_i$ and ${\mathrm{Re}}_o^{\mathrm{alt}}=2{\mathrm{Re}}_o/\eta$ in computing the results with those reported by Hristova et al. [35] as performed by others [36,37]. This is because Hristova et al. [35] used a scale of $d/2$ to nondimensionalize the length and fix the angular rotation speed in defining the problem.

Meseguer [36] and Maretzke et al. [37] reported that close to the region where ${\mathrm{Re}}_i=\eta {\mathrm{Re}}_o$ (i.e., the cylinders rotate at the same speed) no transient growth is observed. However, it is observed that the situation is different when the temperature increases beyond a certain threshold. There exists an amplification of the perturbations at specific times for different relative counter-rotating speeds in the presence of temperature variations. We continue by applying radial heating to the TCF problem and investigate the effect of the thermal stratification for various Grashof, Gr, numbers with specific counter-rotating pairs Re${}_i$, Re${}_o$ and wavenumbers. Figure 3 shows the amplification factor G against t for the four different configurations C${}_1$ to C${}_4$ for various values of Gr.

We find that G has a minimum around Gr = 3250 for configuration C${}_1$, around Gr = 4800 for configuration C${}_2$, around Gr = 5750 for configuration C${}_3$, and around Gr = 7150 for configuration C${}_4$. These results are verified explicitly in Figure 4b. Figure 4 shows that as the temperature increases, the amplification of the perturbations monotonically decay, before reaching a turning point where it enters a new phase and begins increasing monotonically with temperature.

The growth rate shows a linear relationship with Gr in all configurations, as evident from Figure 4a. The growth rate remains consistently negative across all values of Gr, indicating that the fluid flow is linearly stable within the considered range of flows and implying that buoyancy does not induce a linear instability. Figure 4a clearly demonstrates an increase in the growth rate as Gr increases, which aligns with the proportional relationship between $G_0$ and Gr, as illustrated. We observe in Figure 4b that initially there is an observable downward trend in the linear decay of $G_0$ as Gr increases until it reaches a critical value Gr. Beyond this critical value, a linear increase in $G_0$ is observed as Gr continues to rise. By examining the decay phase boundaries, the order of $G_0$ for each configuration can be represented by the expression $G_0\left(C_4\right)>G_0\left(C_3\right)>G_0\left(C_2\right)>G_0\left(C_1\right)$. Similarly, for the increasing phase, the order of $G_0$ can be characterized as $G_0\left(C_4\right)<G_0\left(C_3\right)<G_0\left(C_2\right)<G_0\left(C_1\right)$.

An intriguing observation depicted in Figure 4b is the relationship between $C_1$ and the other configurations, namely, $C_2$, $C_3$, and $C_4$, with respect to their ratio and the magnitude of their growth rates. Notably, there is a bigger difference in the order of growth rates between $C_1$ and $C_2$ when compared to the difference between $C_2$ and $C_3$ across all Gr values. The consistency in the proportionality of their speed ratio suggests that their growth rates should exhibit similar differences, in the same order of magnitude. This observation could be due to the influence of buoyancy-induced contributions captured in the amplification reflected by $G_0$.

The underlying behavior of the fluid motion can be attributed to an internal subcritical phenomenon because of the induced thermal stratification on the dynamics of the disturbances. During the decay phases of the amplification factor, it is evident that the combination of viscous and buoyancy influences the amplitude of the perturbations.

In Figure 5, the amplification factor G is depicted to illustrate its growth as Gr increases. Each row in the sequence of plots represents the evolving behavior of G as Gr increases. The first column has the values Gr $=$ 2250, 3800, 3750, and 5150 and the last column has the values Gr $=$ 4250, 5800, 7750, and 9150. They have Gr values slightly distant from the turning points, Gr $=$ 3250, 4800, 5750, 7150 showcasing a single peak. Additionally, the plots in the second and fourth columns demonstrate two distinct behaviors in their evolution. Notably, the second column of plots display the emergence of a new crest on the wave front. This new crest continues to grow until it surpasses the initial crest in amplitude, eventually becoming the sole peak after the turning point for sufficiently high Gr values. The third column has used a value of Gr where both crescents exhibit similar amplitudes. The initial crest gradually merges into the latter one, and with further increases in Gr, the former completely disappears due to induced buoyancy disturbances.

In Figure 6, we illustrate contour plots of the amplified perturbations in the radial velocity $\breve{u}$ at the value of Gr when $G=G_0$.

The sequence of contour plots follows that of Figure 5 for C${}_1$ to C${}_4$. Notice $t_{-2}$, $t_{-1}$, $t_0$, $t_1$, and $t_2$ follow an arithmetic progression, where $t_0$ is defined as the time when $G=G_0$. We see in Figure 6, that at $G=G_0$ we have a reversal in the rotational direction of the perturbed flow.

5. Conclusions

Over the past twenty years, researchers have widely utilized non-normal analysis to investigate the linear transient characteristics of wall-bounded flows, particularly shear flows. However, there has been a relatively limited focus on applying this approach to examine the behavior of stratified Taylor–Couette flow (STCF).

This study aims to explore the dynamics of perturbations under different temperature conditions and analyze their patterns of amplification. Here we emphasize the correlation between various flow configurations, highlighting their similarity in transient dynamics despite differing speed ratios. Additionally, we investigate the subcritical effects of thermal stratification on the dynamics of the disturbances. The interplay between viscous and buoyancy effects, counteracted by strong centrifugal forces, is found to influence the growth and decay of the amplification factor. As the temperature increases beyond a turning point, buoyancy forces become dominant, leading to a linear increase in the amplification factor. The findings of this study reveal that incorporating temperature into the analysis significantly alters the expected outcomes, deviating considerably from the results obtained by Meseguer [36] and Coles [31]. This supports the notion that buoyancy, induced by the stratification of temperature variation, plays an important role. The transient analysis suggests the presence of complex dynamical phenomena specific to transient dynamics, which cannot be captured by solely focusing on the growth rates of the eigenvalues. An important result of our work is that we have identified the possibility of a two-dimensional subcritical bifurcation in STCF.

Transient growth of linearly stable disturbances is believed to play an important role in the subcritical transition of laminar boundary layers and the self-sustained nature of boundary layer fluctuations in a fully turbulent flow that the modal analysis might not capture within its operational framework. Indeed, our work shows that this is possible for the case of Taylor–Couette flow that has the additional influence of stratification.

In fact, we have performed direct numerical simulations (DNS) and have established that two-dimensional subcritical states are possible for the ubiquitous flow considered in this manuscript. Differences between the linearized theory of the present work and DNS are attributed to nonlinear advection mechanisms. A further investigation of this phenomenon is ongoing.

The transient growth phenomenon refers to an algebraic amplification of small-amplitude disturbances prior to an exponential decay farther downstream, that the full linearized modal mechanism does not capture. Therefore, transient growth has been proposed as a likely mechanism behind laminar–turbulent transition scenarios that cannot be elucidated by the classical paradigm of hydrodynamic instabilities. In fact, our work tackles the case where competition between hydrodynamic and convective instabilities is also present.

We show that the competition between convective and purely hydrodynamic instabilities through transient growth might play an important role in the self-generation of turbulence in fully turbulent wall shear flows through the subcritical generation of nonlinear states as predicted by the fully nonlinear Navier–Stokes equations [39]. This does not mean that transient growth is an alternative to the sequential approach to the turbulent regime, but rather an aid to generating additional branches to such a deterministic approach to turbulence that the sequence of bifurcation approach (SBA) offers. There is of course the alternative approaches by [42] that attempt to explain fully developed turbulence.