Energy Budget Characterisation of the Optimal Disturbance in Stratified Shear Flow

: Stratified Taylor–Couette flow (STCF) undergoes transient growth. Recent studies have shown that there exists transient amplification in the linear regime of counter-rotating STCF. The kinetic budget of the optimal transient perturbation is analysed numerically to simulate the interaction of the shear production (SP), buoyancy flux (BP)


Introduction
The Taylor-Couette system consists of a thin layer of a viscous fluid confined between two vertical rotating coaxial cylinders.The stability conditions for co-rotating cylinders was first examined by [1].Counter-rotating Taylor-Couette flow refers to the flow of a fluid confined between two coaxial and independently rotating cylinders in opposite directions.Stratified Taylor-Couette flow (STCF) introduces an additional layer of complexity by incorporating density variations in the fluid.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 [2][3][4][5][6][7][8][9][10].Understanding the dynamics of stratified Taylor-Couette flow holds relevance in many practical applications, including geophysics, oceanography, and industrial processes in the oil and gas sector.The combination of the horizontal angular shear induced by the cylinders and buoyancy effects due to density stratification gives rise to intriguing patterns, transitions, and stability regimes.In recent years, numerous studies have been conducted to better understand the underlying mechanisms and characteristics of this flow regime [11,12].This flow configuration has been extensively studied for decades and in recent times, primarily in the absence of stratification [13][14][15][16][17]. Conversely, in the presence of stratification, the flow becomes more complex with more degrees of freedom [18,19].
Although linear stability analysis predicts plane Couette and pipe Poiseuille flows to be stable, they exhibit turbulence in practical scenarios for sufficiently large Reynolds numbers, as indicated by [20][21][22], and recent studies [23,24].For plane Poiseuille flow, an instability is predicted to occur beyond a Reynolds number around 5772; in practice, the critical Reynolds number is found to be around 2300.
The inconsistency in these results is attributed to the non-normality of the linearised Navier-Stokes operator, leading to transient growth in the linear regime and triggering the transition from laminar to turbulent flow [25][26][27].
A recent work by Godwin et al. [28] focused 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 analysed the transient stability of STCF by considering small perturbations to 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 it does not capture the physical mechanism that triggers the transient bifurcation phenomena at Grashof number Gr ≈ Gr c .As in [28], the critical value Gr c is the value of Gr when the amplification factor reaches a local minimum.
There is a preference to considering these ubiquitous flows within the modal analysis framework, but this analysis predicts linear stability for all values of the Reynolds number R, which is contrary to experimental observations.Transient growth exists when the eigenvectors of the system are non-normal.To investigate transient growth, non-modal analysis is necessary, and this is the reason that a study based on non-modal analysis is presented here.Additionally, Couette flow is asymptotically stable according to linear theory, and since non-modal analysis determines the largest possible growth of a perturbation in a finite time interval, also called optimal growth, it is hoped that non-modal analysis will bridge the gap between theory and experiment.The initial disturbance yielding optimal growth is called an optimal initial condition, and in our study, we have also allowed for thermal disturbances.The latter ones can be removed in a limiting process, allowing therefore the possibility of connection between linear and nonlinear regimes.
The primary goal of this study is to explore the physical processes responsible for triggering a 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 [28].
Since non-modal analysis is a relatively new method of analysis, it is used here as an important tool to predict instabilities.Nowadays, it is understood that a perturbation in a shear flow can experience significant transient energy growth.This growth is responsible for the initial linear amplification of disturbances and is mainly responsible for the appearance of a subcritical transition to turbulence.Non-modal effects can therefore explain the discrepancy between the observed experimental results in incompressible shear flows and the theoretically predicted critical Reynolds number for a linear instability in wall-bounded shear flows.We therefore need to use non-modal analysis to gain further insight into the linear stability of Couette flow.
It is reported that as Gr increases, the maximum amplification factor initially decays before it eventually starts growing after reaching a critical value Gr c .Assuming other parameters of the model are kept constant, the value Gr c changes as a function of different speed ratios of the cylinders.The observed decay/growth pattern of G 0 is attributed to the interplay between the induced shear and buoyancy.In this study, the mechanism that triggers this observed pattern is investigated.
In Section 1, a brief overview of the investigated problem is presented.Section 2 focuses on deriving, linearising, and numerically discretising 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 summarises the conclusions drawn from this study.

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 Ω i > 0 and is held at the temperature T i = T + ∆T/2, where T is the ambient temperature and ∆T is the temperature difference between the two cylinders.The outer cylinder has radius r o , rotates with an angular velocity of Ω o < 0 and is held at the temperature T o = T − ∆T/2.We assume that the density ρ varies linearly with the temperature such that where β = −(1/ρ o )∂ρ/∂T| ρ=ρ o is the coefficient of thermal expansion.We assume that the dynamic viscosity of the fluid µ and thermal conductivity κ are constant, whilst we define ν = µ/ρ o as the kinematic viscosity.We notice we can write 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 (bold here denotes a vector) where v, Ť, x, ť and p are the dimensionless velocity, temperature, spatial coordinate, time and pressure, with d = r o − r i .By defining η = r i /r o we can express the dimensionless inner and outer radii as η/(1 − η) and 1/(1 − η), respectively.By dropping the accents for convenience, the resulting dimensionless equations become where e z is the unit vector in the vertical axial direction, and the Grashof number, Prandtl number and relative density are given by , and ϵ = β∆T.
The dimensionless boundary conditions are given as where e θ is the unit vector in the azimuthal direction, and the inner and outer Reynolds numbers are given by Due to the cylindrical boundaries, we perform our calculations in cylindrical polar coordinates.The steady base state velocity takes the form In order to have a constant pressure gradient in the axial direction, we impose the zero mass flux constraint, namely The well-known analytical base state solution to this problem is given by where the constants A, B and C are given by and Now, we consider small perturbations to the base state in the form v = U + δ ṽ, p = P b + δ p, and where δ is assumed to be a small constant, P b is the base state pressure, not presented here for convenience, and ṽ is the perturbation to the velocity.Substituting Equation ( 12) into the dimensionless Equations ( 2)-( 4) and linearising in δ, we obtain ∇ • ṽ = 0, ( To perform a stability analysis on this system, we seek perturbations of the form ṽ = v(r, t)e i(nθ+kz) + c.c., p = p(r, t)e i(nθ+kz) + c.c., T = T(r, t)e i(nθ+kz) + c.c.
where n ∈ N and k ∈ R are the azimuthal and axial wavenumbers, respectively.Substituting these expansions into the linearized Equations ( 13)-( 15) and simplifying yields where By letting q = [ û, v, ŵ, p, T] T , we can write Equations ( 16)- (20) as where 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 and By assuming a solution of the form v = v(r)e λt , p = p(r)e λt , and T = T(r)e λt , ( we transform the initial value problem to a generalized eigenvalue problem: The scalar eigenvalue λ is complex and defines the temporal stability of the flow.That is, if the real part of λ is negative, the flow is stable and the amplitude of the perturbations will decay in time.Conversely, if the real part of λ 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: In other words, we impose a homogeneous boundary condition for the velocity and temperature perturbations at the respective walls of the cylinders.

Perturbation Energy Budget
The asymptotic stability approach employed in perturbation analysis can effectively capture the length and time scales of unstable modes [29,30].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 [25,[31][32][33].Exploring the processes influencing the kinetic energy distribution of perturbations can provide valuable insights [34][35][36].
The presence of kinetic energy in thermal stratified shear flow significantly impacts the perturbation's decay or growth process [35,37].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.
To calculate the kinetic energy of the optimal perturbation, we multiply the linearized Equations (A1)-(A4) by the complex conjugate of the optimal perturbation and evaluate the average over the entire domain in azimuthal and axial coordinating directions.Using the singular value decomposition method we computed the largest (negative) eigenvalues of equation (25).The associated eigenvectors to these eigenvalues are the optimal perturbations.One notes that eigenvalues remain negative as this problem is predicted to be linearly stable even though full nonlinear simulations and experiments show this problem to be unstable.This process yields the following energy budget equation: where KE is the kinetic energy of the perturbation, SP is the shear production, BP is the buoyancy production, EF is the energy flux, CKEV is the convergence of kinetic energy of the optimal perturbation due to viscosity and VD is the viscosity dissipation.All of these terms are defined in Appendix A.

Results and Discussion
Before discussing our results, it is 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 a 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 the energy budget quantities within the framework of optimal transient perturbations particularly in relation to the amplification factor.
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 four different configurations are given.Furthermore, the critical value Gr c for each configuration is also included.The analysis is carried out over a range of Grashof number values along with the following values: Pr = 68, ϵ = 0.067 and η = 0.881, with values of Re i and Re o given in Table 1.These values are chosen to allow a comparison of the current results with those in the literature and correspond to the experimental values in [38].The results obtained for the above parameter values are illustrative of the results we obtained from a larger number of numerical simulations.To obtain sufficiently accurate results, 128 Chebyshev nodes were used to compute the energy quantities.
The least stable modes for each configuration are first identified.The energy quantities are computed using the optimal disturbances as inputs for each of the selected Gr values.These values are chosen for illustrative purposes to show how the energy quantities oscillate between the cylinders depending on the value of Gr.For a more detailed discussion on how the amplification factor depends on Gr, please see Figure 4 of [28].
Figure 1 illustrates the energy quantities for the C 1 configuration for two values of Gr, namely Gr = 2500 (blue) and Gr = 3500 (red).The curves are presented at the time when the amplification factor is optimal.
In Figure 1, the blue curve corresponds to Gr < Gr c whilst the red curve corresponds to Gr > Gr c .All of the blue curves in Figure 1 show peaks concentrated around the outer cylinder, whilst the red curves show peaks concentrated around the inner cylinder.Notice that the peaks of BVP are at least 10 3 times smaller than the peaks of the other six quantities.When BP is positive, buoyancy causes the warmer, less dense fluid to rise, while the cold dense fluid sinks.We found that configurations C 2 , C 3 and C 4 produced very similar results to those of C 1 .
Figure 1.The optimal perturbation energy (temperature) budget quantities for the configuration C 1 ; see Table 1.The blue curve corresponds to Gr = 2500, whilst the red curve corresponds to Gr = 3500.The curves are presented at the time when the amplification factor is optimal.
In Figure 2, we illustrate contour plots in the r-z plane of the perturbations to the base state temperature at the optimal wavenumber for θ = 0.68.The values of Gr above each contour plot in the left column correspond to values of Gr less than Gr c , the middle column corresponds to values of Gr within 5% of Gr c , whilst the right column corresponds to values of Gr greater than Gr c .The rows from top to bottom correspond to the configurations The left column in Figure 2 shows that the extrema in the perturbed temperature occurs towards the outer cylinder for values of Gr < Gr c .Conversely, the right column in Figure 2 shows that the extrema in the perturbed temperature occurs towards the inner cylinder for values of Gr > Gr c .The perturbations to the angular and vertical velocity were very similar to that of the temperature comtour plots and thus are not presented here.These results are consistent with the energy budget quantities shown in Figure 1.
In Figure 3, we illustrate contour plots in the r − z plane of the perturbations to the base state radial velocity at the optimal wavenumber for θ = 0.68.The values of Gr above each contour plot in the left column correspond to values of Gr less than Gr c , the middle column corresponds to values of Gr within 5% of Gr c , whilst the right column corresponds to values of Gr greater than Gr c .The rows from top to bottom correspond to the configurations C 1 -C 4 .The left column in Figure 3 shows that the extrema in the perturbed radial velocity occurs towards the outer cylinder for values of Gr < Gr c .Conversely, the right column in Figure 3 shows that the extrema in the perturbed radial velocity occurs towards the inner cylinder for values of Gr > Gr c .Again, these results are consistent with the energy budget quantities shown in Figure 1.We notice that there are alternating convective rolls arranged vertically above each other, which is consistent with the vortices which appear as secondary states in the Taylor-Couette flow.

Conclusions
In this study, we have revisited the stratified Taylor-Couette flow via the non-modal stability analysis the singular value decomposition method and via computing the largest negative eigenvalues of equation (25).Sequentially we used the energy budget equation to perform our calculations.Our results are presented for values of the Grashof number less than, around, and greater than the critical Grashof number, where the critical Grashof number has been defined as the value which minimises the optimal amplification factor.We found that when Gr < Gr c , then the perturbations to the base state flow are focused around the outer wall; conversely, when Gr > Gr c , then the perturbations to the base state flow are focused around the inner wall.
The complex dynamics of mixed convection flow with a particular focus on shear instability and its manifestation in the presence of stratified conditions were investigated with 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.
Our 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 SP and BF near the outer cylinder wall induce 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 towards 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.
Our analysis also highlighted the role of buoyancy and shear drag in triggering shorttime and long-time instability.We found that 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 an 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.
The perturbation functions change shape even for small Gr increments as shown in Figures 2 and 3 due to the proximity of the optimal amplification factor.Small changes indicate the change of stability and the possibility of existence of nonlinear solutions.Sometimes shifted to the inner wall and sometimes to the outer wall the flow meanders as it rearranges itself for effective mass transport.This work is slightly distinctive to modal analysis where the eigenvalues indicate the critical mode(s) and indicates an alternative transfer to nonlinearity.
Furthermore, our findings have shown that the viscous processes CKEV and VD 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.

Figure 2 .
Figure 2. Contours of the perturbations to the base state temperature of the optimal perturbation in the r-z plane for θ = 0.68.Configurations C 1 , C 2 , C 3 , and C 4 are organized in rows from top to bottom, for selected Gr values, as indicated above each contour plot.

Figure 3 .
Figure 3. Contours of the perturbations to the base state radial velocity of the optimal perturbation in the r-z plane for θ = 0.68.Configurations C 1 , C 2 , C 3 , and C 4 are organized in rows from top to bottom, for selected Gr values, as indicated above each contour plot.

Author Contributions: 2 ∂
Conceptualization, L.E.G.; methodology, L.E.G. and S.C.G.; software, L.E.G.; validation, L.E.G. and S.C.G.; formal analysis, L.E.G. and S.C.G.; investigation, L.E.G. and S.C.G.; resources, L.E.G. and S.C.G.; writing-original draft preparation, L.E.G.; writing-review and editing, L.E.G., P.M.J.T., T.A. and S.C.G.; supervision, P.M.J.T. and S.C.G.; project administration, P.M.J.T. and S.C.G.; funding acquisition, L.E.G. and S.C.G.All authors have read and agreed to the published version of the manuscript.Funding: This research was funded by RISE Horizon 2020 ATM2BT, Atomistic to Molecular Turbulence, Grant No. 824022, PTDF ED/PHD/GLE/826/16 scholarship and DTI EPSRC grant, Aston University sponsorship.each of the first three equations by the associated velocity component and then take the sum of these equations we can integrate over the axial and azimuthal directions and assume periodic boundary conditions in both those directions to give If we multiply Equation (A4) by T and average as before, we obtain 1 T2 ∂t = − < ȗ T > ∂T b ∂r + 1 Pr < T∇ 2 T > .