Energy Budget Characterisation of the Optimal Disturbance in Stratified Shear Flow

1. Introduction

Counter rotating Stratified Taylor Couette Flow (STCF) is a complex fluid dynamics phenomenon that has attracted significant attention from researchers due to its relevance in various industrial and environmental applications. In recent years, numerous studies have been conducted to better understand the underlying mechanisms and characteristics of this flow regime [1,2]. Counter rotating Taylor Couette Flow refers to the flow of a fluid confined between two coaxial and independently rotating cylinders in opposite directions. This flow configuration has been extensively studied for decades and in recent times, primarily in the absence of stratification [3,4,5,6,7]. Conversely, in the presence of stratification the flow becomes more complex with more degrees of freedom in the parametric study of the dynamic behaviour [8,9]. Stratification refers to the presence of density gradients in the fluid, which can occur due to variations in temperature, salinity, or other properties.The buoyancy effect, resulting from the density differences caused by stratification, introduces a new dimension to the dynamics of STCF. It interacts with the rotational forces, and its impact on the stability of the flow has been the subject of considerable research interest [10,11,12,13,14,15,16,17,18].

A recent work by Godwin et al. [19] focussed specifically on the short-time instability of STCF. To explore the short-time instability of STCF, the authors conducted a detailed numerical investigation employing transient growth analysis techniques. The authors analyzed the transient stability of STCF by considering small perturbations superimposed on the base flow. By systematically varying the Grashof number for various model configurations, they were able to quantify the influence of these parameters on the transient instability characteristics. The study sheds light on the transient instability induced by buoyancy in counter rotating STCF within the bounds of the linear regime considered, but does not capture the physical mechanism that triggers the transient bifurcation phenomena at $Gr\approx G{r}_c$. The value $G{r}_c$ is the critical value for which the transient bifurcation occurs with respect to the model configuration.

The primary goal of this study is to explore the physical processes responsible for triggering transient bifurcation, with specific objectives outlined as follows: developing a linear model for STCF and conducting a numerical evaluation using the Chebyshev collocation method, determining the least stable eigenvalue for each configuration, identifying the maximum amplification factor and optimal perturbation for the least stable mode, and ultimately deriving and evaluating the energy budget quantities associated with the optimal perturbation. This paper primarily focuses on the last objective, as the other aspects have already been satisfactorily addressed in previous work [19].

In Section 1, a brief overview of the investigated problem is presented to facilitate a thorough understanding of the topic. Section 2 focuses on deriving, linearizing, and numerically discretizing the governing equations of the STCF model. The discussion of optimal transient budget quantities is covered in Section 3, with detailed derivations provided in Appendix A. Section 4 delves into a comprehensive analysis of the investigation’s results. Finally, Section 5 summarizes the conclusions drawn from this study.

2. Mathematical Model

In this study, we consider an incompressible fluid between two concentric rotating heated vertical cylinders of infinite height. The inner cylinder has radius $r_i$, rotates with an angular velocity of ${\Omega }_i>0$ and is held at the temperature $T_i=\overline{T}+\Delta T/2$. The outer cylinder has radius $r_o$, rotates with an angular velocity of ${\Omega }_o<0$ and is held at the temperature $T_o=\overline{T}-\Delta T/2$. Here, $\overline{T}$ is the ambient temperature, and $\Delta T$ is the temperature difference between the two cylinders. We assume that the density $\rho$ varies linearly with the temperature such that

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

(1)

where $\beta =-(1/{\rho }_o){\partial \rho /\partial T|}_{\rho ={\rho }_o}$ is the coefficient of thermal expansion. We assume that the dynamic viscosity of the fluid $\mu$ and thermal conductivity $\kappa$ are constant, whilst we define $\nu =\mu /{\rho }_o$ as the kinematic viscosity. We notice we can write $\rho ={\rho }_o+{\rho }^{\prime }$ where ${\rho }^{\prime }=-{\rho }_o\beta (T-\overline{T})$.

The presence of temperature gradients with the rotating concentric cylinders induces complex flow patterns and thermal effects. The variation in the density, due to the temperature gradients, induces various buoyancy-driven flow regimes. The equations are nondimensionalised using

$$\mathbf{v}=\frac{\nu }{d}\check{\mathbf{v}},T=\overline{T}+(\Delta T)\check{T},\mathbf{x}=d\check{\mathbf{x}},t=\frac{d^2}{\nu }\check{t},p=p_a+\frac{{\rho }_o{\nu }^2}{d^2}\check{p},$$

where $\check{\mathbf{v}}$, $\check{T}$, $\check{\mathbf{x}}$, $\check{t}$ and $\check{p}$ are the dimensionless velocity, temperature, spatial coordinate, time and pressure, with $d=r_o-r_i$. By defining $\eta =r_i/r_o$ we can express the dimensionless inner and outer radii as $\eta /(1-\eta )$ and $1/(1-\eta )$, respectively. By dropping the accents for convenience, the resulting dimensionless equations become

$${\mathbf{v}}_t+(\mathbf{v}·\nabla )\mathbf{v}=-\nabla p+{\nabla }^2\mathbf{v}+GrT{\mathbf{e}}_z+ϵT(\mathbf{v}·\nabla )\mathbf{v},$$

(2)

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

(3)

$$T_t+(\mathbf{v}·\nabla )T=\frac{1}{Pr}{\nabla }^{2}T,$$

(4)

where ${\mathbf{e}}_z$ is the unit vector in the vertical axial direction and the Grashof number, Prandtl number and relative density are given by

$$Gr=\frac{\beta g\Delta T{d}^3}{{\nu }^2},Pr=\frac{\nu }{\kappa },andϵ=\beta \Delta T.$$

The dimensionless boundary conditions are given as:

$$\mathbf{v}={Re}_i{\mathbf{e}}_{\theta }onr=\frac{\eta }{1-\eta }and\mathbf{v}={Re}_o{\mathbf{e}}_{\theta }onr=\frac{1}{1-\eta },$$

(5)

where ${\mathbf{e}}_{\theta }$ is the unit vector in the azimuthal direction and the inner and outer Reynolds numbers are given by

$${Re}_i=\frac{{\Omega }_i{r}_{i}d}{\nu }and{Re}_o=\frac{{\Omega }_o{r}_{o}d}{\nu }.$$

Due to the cylindrical boundaries, we shall change to cylindrical polar coordinates. The steady base state velocity takes the form

$$\mathbf{U}\left(r\right)=V_b{\mathbf{e}}_{\theta }+W_b{\mathbf{e}}_z.$$

(6)

In order to have a constant pressure gradient in the axial direction, we impose the zero mass flux constraint, namely

$${\int }_{\frac{\eta }{1-\eta }}^{\frac{1}{1-\eta }}r{W}_{b}dr=0.$$

(7)

The well known analytical base state solution to this problem is given by

$$\begin{array}{ccc}\hfill V_b& =& Ar+\frac{B}{r},T_b=\frac{ln\left(\right(r/\eta \left)\right(1-\eta \left)\right)}{ln\left(\eta \right)}+\frac{1}{2},\hfill \\ \hfill W_b& =& \frac{Gr}{{(1-\eta )}^2}[C({\left(r(1-\eta )\right)}^2-{\eta }^2)\hfill \end{array}$$

(8)

$$\begin{array}{cc}+& \left(C(1-{\eta }^2)+\frac{1}{8}(1-{\left(r(1-\eta )\right)}^2)\right)\left(2T_b\left(r\right)-1\right)],\hfill \end{array}$$

(9)

where the constants A, B and C are given by

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

(10)

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)}.$$

(11)

Now, we consider small perturbations to the base state in the form

$$\mathbf{v}=\mathbf{U}+\delta \tilde{\mathbf{v}},p=P_b+\delta \tilde{p},andT=T_b+\delta \tilde{T}$$

(12)

where $\delta$ is assumed to be a small constant, $P_b$ is the base state pressure, not presented here for convenience, and $\tilde{\mathbf{v}}$ is the perturbation to the velocity. Substituting equation (12) into the dimensionless equations (2) to (4) and linearising in $\delta$ we obtain

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

(13)

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

(14)

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

(15)

To perform a stability analysis on this system we seek perturbations of the form

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

where $n\in \mathbb{N}$ and $k\in \mathbb{R}$ are the azimuthal and axial wavenumbers, respectively. Substituting these expansions into the linearized equations (13) - (15) and simplifying yields

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

(16)

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

(17)

$$\begin{array}{ccc}\hfill \frac{\partial \widehat{w}}{\partial t}& =& \left(\mathcal{F}-\mathcal{H}\right)\widehat{w}-ik\widehat{p}-\frac{\partial W_b}{\partial r}(1-ϵT_b)\widehat{u}+Gr\widehat{T},\hfill \end{array}$$

(18)

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

(19)

$$\begin{array}{ccc}\hfill \frac{\partial \widehat{T}}{\partial t}& =& \frac{1}{Pr}\left[\mathcal{F}-iPr\left(\frac{n{V}_b}{r}+k{W}_b\right)\right]\widehat{T}-\frac{\partial T_b}{\partial r}\widehat{u},\hfill \end{array}$$

(20)

where ${\displaystyle D=\frac{\partial }{\partial r}}$, ${\displaystyle D_{+}=D+\frac{1}{r}}$ and

$$\mathcal{F}=D_{+}D-\frac{n^2}{r^2}-k^{2}and\mathcal{H}=i\left(\frac{n{V}_b}{r}+k{W}_b\right)(1-ϵT_b).$$

By letting $\widehat{\mathbf{q}}={[\widehat{u},\widehat{v},\widehat{w},\widehat{p},\widehat{T}]}^T$ we can write equations (16)-() as

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

(21)

where

$$A=\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1\end{array}\right]$$

(22)

and

$$B=\left[\begin{array}{ccccc}\mathcal{F}-\mathcal{H}-\frac{1}{r^2}& \frac{2V_b}{r}(1-ϵT_b)-\frac{2in}{r^2}& 0& -D& -\frac{ϵV_b^2}{r}\\ \frac{2in}{r^2}-(1-ϵT_b)\left(\frac{V_b}{r}+\frac{\partial V_b}{\partial r}\right)& \mathcal{F}-\mathcal{H}-\frac{1}{r^2}& 0& -\frac{in}{r}& 0\\ -\frac{\partial W_b}{\partial r}(1-ϵT_b)& 0& \mathcal{F}-\mathcal{H}& -ik& Gr\\ D_{+}& \frac{in}{r}& ik& 0& 0\\ -\frac{\partial T_b}{\partial r}& 0& 0& 0& \frac{{\displaystyle \mathcal{F}}}{Pr}-i\left(\frac{n{V}_b}{r}+k{W}_b\right)\end{array}\right].$$

(23)

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},and\widehat{T}=\breve{T}\left(r\right)e^{\lambda t},$$

(24)

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

$$\lambda \mathbf{A}\widehat{\mathbf{q}}=\mathbf{B}\widehat{\mathbf{q}},\Rightarrow |\lambda \mathbf{A}-\mathbf{B}|=0.$$

(25)

The scalar eigenvalue $\lambda$ is complex and defines the temporal stability of the flow. That is if the real part of $\lambda$ is negative the flow is stable and the amplitude of the perturbations will decay in time. Conversely, if the real part of $\lambda$ is positive the flow is unstable and the amplitude of the perturbation will grow asymptotically in time. Further more, in order to define the boundary conditions, 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)=\mathbf{0}and\breve{T}\left(r_i\right)=\breve{T}\left(r_o\right)=0.$$

(26)

In other words, we impose homogeneous boundary condition for the velocity and temperature perturbations at the respective walls of the cylinders.

3. Perturbation Energy Budget

The asymptotic stability approach employed in perturbation analysis can effectively capture the length and time scales of unstable modes [20,21]. However, this method lacks an intuitive physical interpretation of the internal processes driving the instability. To gain a deeper understanding of the physical underpinnings of the instability, it is advantageous to adopt a more intuitive approach [22,23,24,25]. Exploring the processes influencing the kinetic energy distribution of perturbations can provide valuable insights [26,27,28].

The presence of kinetic energy in thermal stratified shear flow significantly impacts the perturbation’s decay or growth process [27,29]. The perturbation process is characterized by the conservation or transformation of energy across different states. For perturbations to grow, an energy exchange must occur. In the case of shear flow, perturbations grow by accessing kinetic energy. Investigating the quantities of energy sources interacting to form perturbation kinetic energy allows for a further understanding of the nature of induced instability—whether it is induced by shear or buoyancy—due to transient growing modes. When induced buoyancy is greater than zero, buoyant fluid rises near the wall, and dense fluid sinks. This process is also influenced by shear quantity. If the shear contribution of the energy source is less than zero, it implies that shear opposes induced convection due to buoyancy.

To calculate the kinetic energy of the optimal perturbation, we multiply the linearized equations (A1)-() by the complex conjugate of the optimal perturbation and evaluate the average over the entire domain in azimuthal and axial coordinating directions. This process yields the energy budget equation:

$$\frac{\mathsf{d}KE}{\mathsf{d}t}=\mathrm{SP}+\mathrm{BP}+\mathrm{EF}+\mathrm{CKEV}+\mathrm{VD}.$$

(27)

Please see appendix A for the derivation.

4. Result

Before delving into the discussion of our results, it’s important to acknowledge that in mixed convection flow, shear instability can manifest when there are relatively moving shear layers of fluid. The perturbations that arise from this instability may experience damping due to stable stratification because the perturbation’s energy must be expended to counteract the force of gravity induced by the stratification. Consequently, the growth rate of these perturbations will decrease, and the instability will be suppressed if the stratification is sufficiently strong. These intricate processes are elucidated through the examination of energy budget quantities within the framework of optimal transient perturbations, particularly in relation to the amplification factor, as elaborated upon in the subsequent sections.

The kinetic budget quantities of the optimal growth perturbation are numerical investigated for different stratified STCF configurations. In Table 1 the parameter values for the 4 different configurations are given.

The analysis is carried out with fixed values for Pr=68, $ϵ$=0.08, $\eta$ = 0.881, Re, and selected Gr. To obtain sufficiently accurate results, 128 Chebyshev nodes were used to compute the energy quantities.

The results are computed in a systematic manner. The least stable modes for each configuration are first identified followed by the optimal perturbation. The energy quantities are computed using the optimal disturbances as inputs for the selected Gr values. These values are chosen for illustrative purposes to show that the energy quantities oscillate between the cylinders before and after the turning points. For a more detailed discussion of these turning points please see Figure 4 in Godwin et al. [19].

Figure 1 illustrates energy quantities prior to the turning point for configurations $C_1$ and $C_2$.

The terms KE, SP, BP, EF, CKEV, VD and BVP are defined in appendix A. We found that configurations $C_3$ and $C_4$ produced very similar results. The values of Gr are deliberately selected to be significantly smaller than the turning point, that is Gr $\ll {Gr}_c$. In the proximity of the outer cylinder wall, induced SP and BF draw KE from the mean flow into the perturbation, causing energy to radially flux inwards towards the center of the gap between the cylinders through EF. This energy, when spread by EF, converges/diverges at specific radial locations close to the outer cylinder wall. Moreover, Figure 1 reveals a sharp peak in KE near the outer cylinder wall due to the high values of SP and BF. The positive implication of BF indicates that buoyant fluid rises while dense fluid sinks near the outer cylinder wall. Conversely, the buoyancy flux resulting from the buoyancy variance in the heat equation is negatively amplified, albeit with a very small magnitude of about ${10}^{-5}$.

Figure 2 depicts energy quantities around the turning point for configurations $C_1$ and $C_2$.

We found that configurations $C_3$ and $C_4$ produced very similar results. Depending on the Gr values around this point, the kinetic energy (KE) exhibits a concentrated distribution of significant peaks either towards the outer or inner cylinder, with a minor bump on the tail side of the amplitude. No specific value was identified while varying the Gr values across the transition from the decaying to the growth phase of the amplification factor, $G_0$, due to transient discontinuity around the turning point. Additionally, the energy flux (EF) displays converging/diverging peaks, with one peak being more pronounced than the other in relation to the Gr value near that point. This suggests that EF evenly spreads the energy from the center to the cylinder walls.

Figure 3 illustrates energy quantities beyond the turning point for configurations $C_1$ and $C_2$.

We found that configurations $C_3$ and $C_4$ produced very similar results. The values of Gr have been selected to be significantly larger than the turning point, denoted as Gr $\gg {Gr}_c$. Near the wall of the inner cylinder, the induced SP and BF extract KE from the mean flow into the perturbation, causing the energy to flux radially inwards towards the center of the inner cylinders through EF. The spreading of energy by EF leads to convergence/divergence at specific radial locations near the wall, $r_i$. Additionally, Figure 3 also shows a sharp peak in KE near $r_i$ due to the high concentration of SP and BF. The negative implication of BF implies that shear-induced instability must exert more work against gravity to grow by suppressing buoyant fluid and lifting dense fluid near $r_i$. Conversely, the buoyancy flux resulting from the buoyancy variance in the heat equation is positively amplified, albeit with a very small magnitude of about ${10}^{-5}$.

The opposite signs of the amplitude of BP and SP may trigger showtime or long-time instability depending on the proportionality of magnitude of each amplification. We see that has because the disparity between the SP and BP values differed significantly it tends to cause a short time amplification. However, for a sufficiently high value of BP, the short time amplification dies as shown in Figure 4.

This long-time behaviour is what is capture by the asymptotic analysis of the eigenvalue. Although, it can be tempting to conclude that buoyancy is the trigger for the short time instability, but we argue otherwise, that it this the shear drag as it opposes the buoyancy to triggers the short time instability as evident in the long time run captured in Figure 4.

The viscous processes (CKEV and VD) is unable to generate the optimal perturbation kinetic energy (KE) akin to the energy flux (EF). Instead, these viscous terms disperse energy without adding to the overall KE, and they vanish at the boundaries. The dissipative terms (VD) exhibit a negative definiteness, leading to the decay of KE. It’s worth noting that viscosity functions as an energy sink for the perturbation, compelling the velocity at the boundary to gradually diminish to zero. This process results in a distortion of the perturbation, causing the wall-normal velocity ($u_r$) of the perturbation to incline against the transient background shear shown in Figure 4, ultimately yielding a positive shear production (SP) amplification. While viscosity doesn’t directly contribute to KE, it indirectly influences it by configuring the perturbation in a manner that enables it to extract KE from the background mean flow.

In Figure 5 we illustrate the vertical vorticity for the four configurations.

In Figure 5, the left column (for values of Gr before the turning point), shows that the largest variations in the vertical vorticity are in the region near the outer cylinder, however, the right column (for values of Gr after the turning point), shows that the largest variations in the vertical vorticity are in the region near the inner cylinder.

In Figure 6 we illustrate the horizontal velocity in the z plane.

Figure 6, the left column (for values of Gr before the turning point), shows that the largest variations in the horizontal velocity are in the region near the outer cylinder, however, the right column (for values of Gr after the turning point), shows that the largest variations in the horizontal velocity are in the region near the inner cylinder.

In Figure 7 we illustrate the vertical velocity in the azimuthal plane.

Figure 7, the left column (for values of Gr before the turning point), shows that the largest variations in the vertical velocity are in the region near the outer cylinder, however, the right column (for values of Gr after the turning point), shows that the largest variations in the vertical velocity are in the region near the inner cylinder.

5. Conclusion

The complex dynamics of mixed convection flow with a particular focus on shear instability and its manifestation in the presence of stratified conditions is investigated further in this paper. The analysis involves a systematic numerical investigation of kinetic budget quantities for optimal growth perturbations in various stratified configurations of Taylor-Couette flow. The key findings shed light on the interplay of different energy components and their effects on the transient growth.

The results reveal distinct behaviours of energy quantities before, around, and beyond the turning point in the amplification factor. Prior to the turning point, the values of shear production (SP) and buoyancy flux (BF) near the outer cylinder wall induces significant kinetic energy (KE) peaks. This positive amplification of BF indicates a buoyant fluid rise, while the negative amplification suggests shear-induced instability opposing gravity. Around the turning point, the kinetic energy distribution exhibits peaks either toward the outer or inner cylinder, emphasizing the transient nature of the amplification process. Beyond the turning point, the energy quantities show similar amplification peaks, with SP and BF near the inner cylinder wall causing significant KE peaks.

The analysis also highlights the role of buoyancy and shear drag in triggering short-time and long-time instability. The disparity between BF and SP values significantly influences the nature of amplification, with a sufficiently high shear production leading to the decay of short-time amplification. This long-time behavior is captured by asymptotic analysis of the eigenvalue. Contrary to initial assumptions, the study argues that shear drag, rather than buoyancy, is the primary trigger for short-time instability.

Furthermore, the viscous processes (CKEV and VD) are found to disperse energy without contributing to the overall KE, while dissipative terms (VD) lead to the decay of KE. Viscosity acts as an energy sink, gradually diminishing the velocity at the boundary and influencing the perturbation’s configuration to extract KE from the background mean flow.

In summary, our analysis provides insight into the intricate processes governing mixed convection flow with shear instability, emphasizing the role of buoyancy, shear drag, and viscosity in shaping the flow dynamics. These findings contribute to a deeper understanding of the underlying mechanisms in complex fluid systems.