Stability Criteria to the Incompressible Inviscid Linear Fluid between Two Rotation Coaxial Cylinders

Hong Yang^*, Hao Liu^#

Department of Mathematical Sciences, Jinan University, Guangzhou, China.

DOI: 10.4236/jamp.2021.92020                

Abstract

The stability and instability phenomenon coupled with the rotation effect and the thermal convection between two concentric cylinders was studied. By means of the Normal-Modes method, the stability or instability criteria for the linearized system in terms of the oscillation frequency, the axial wavelength and the background thermal gradient are proved. Besides, some numerical simulation for the axisymmetric perturbations is presented.

Keywords

Stability or Instability Criteria, Linearized System, Rotation Flow

1. Introduction

The study of the hydrodynamic stability has a long history, one of the oldest problems considered is the stability and instability of shear flows, for example, Rayleigh [1], as well as by many other authors with new perspectives (see [2] and references therein). The stability or instability phenomenon for the fluid system has attracted renewed interest during the past decades, due to the mathematical challenges and many interesting physical phenomena they present [3]. Since the pioneering works of Couette and Taylor (see [4] [5] ), thousand of experimental, numerical and theoretical studies have considered different aspects of circular Couette flow [3] [6] [7]. The Taylor-Couette flow is a canonical and popular flow, it has also led to a very large number of studies and significant advances in the understanding of fluid stability [8] and transitions to turbulence [9] [10].

The Taylor-Couette flow is mainly to consider the sheared flow between two independently rotating concentric cylinders. However, Taylor-Couette flow with axial thermal stratification has received little attention [3] [11] [12], which combines horizontal shear and thermal convection, is of great interest in astrophysics, for example, the stratified Taylor-Couette flow serves as a model for instabilities in equatorial oceans [13] and [14]. The purpose of this paper is to study the Taylor-Couette flow with thermal convetion.

The study of the hydrodynamic stability for rotation flow (see Figure 1) begins with Rayleigh [1] (see also the overview in [15] ). In the work of Rayleigh [1], a necessary and sufficient condition for the stability of linearized system is stated as

$\frac{\text{d}}{\text{d}r}{\left(r^2\Omega \right)}^2>0,$ (1)

where $\Omega \left(r\right)$ is the angular velocity of the fluid at the distance r ( $r_1\le r\le r_2$ ) from the axis. More precisely, when the cylinders rotate in the same direction, only one mode of instability is present, which corresponds to the convection mode. When the cylinders rotate in opposite directions, two types or instability are presented. The second instability mode is of an oscillatory type. These results suggest that the “exchange of stability” may be valid when the cylinders rotate in the same direction, while it may not be valid when the cylinders rotate in opposite direction [16].

When the thermal varies in this area, the convection triggered by the thermal variation, also named as Bernard convection, will cause instability. There were many efforts about the stability or not on this case, see more details in [3] [17]. As to the coupled system with both rotation and the thermal convection, in [2] [18], the case has been studied in a rotation coordinate. However, in the work of [2] [18], in terms of the rotation and thermal convection, which one is the dominant role of the instability is not very obvious.

In this paper, we shall study the stability or instability criteria of the linearized incompressible inviscid fluid with both rotation and vertical background thermal variation. The vertical thermal variation is commonly observed in the atmosphere and in the oceans. Certainly, there also exists a horizontal variation of thermal across the latitudes due to differential heat radiation by the sun [19], however, in a small scale region, it is reasonable to consider only vertical thermal variation case. Namely, we shall study the following system:

$(\begin{array}{l}{\partial }_t\rho +\left(U\cdot \nabla \right)\rho =0,\\ {\partial }_{t}U+\left(U\cdot \nabla \right)U+\nabla P=\rho e_3,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}x\in {ℝ}^3,t>0,\\ \nabla \cdot U=0,\end{array}$ (2)

where $U={\left(U_1\mathrm{,}{U}_2\mathrm{,}{U}_3\right)}^{\text{T}}$ is the velocity. The scalers $\rho$ and P be the thermal and the pressure respectively. The system (2) also named as Boussinesq system, which is widely used to model the dynamics of the ocean or the atmosphere, see [20]. This system includes the weak nonlinear and dispersive effects, it can effectively interpret the dispersion wave in atmospheric dynamics [19] [21].

Figure 1. The geometry of two rotation coaxial cylinders.

The goal of present work is to understand the stability problem of system (2) connected with the rotation and the thermal variation. The main theorems state as following:

Theorem 1.1 The necessary and sufficient condition to be stable of the linearized system of (2) with axisymmetric perturbation is

$\frac{\text{d}\bar{\rho }}{\text{d}z}+k^2\frac{{\displaystyle {\int }_{r_1}^{r_2}}r{\xi }_r^2\Phi \left(r\right)\text{d}r}{{\displaystyle {\int }_{r_1}^{r_2}}r{\left(\frac{\text{d}}{\text{d}r}{\xi }_r+\frac{1}{r}{\xi }_r\right)}^2\text{d}r}>0.$ (3)

To be more precisly, our results state:

· The sufficient condition for the linear stable and unstable with axisymmetric perturbation:

1) When the rotation effect $\Phi \left(r\right)>0$ (defined in (34)) and $\frac{\text{d}\bar{\rho }}{\text{d}z}\ge 0$ in $r_1\le r\le r_2$, then the axisymmetric perturbations are always linear stable;

2) When the rotation effect $\Phi \left(r\right)>0$ and $\frac{\text{d}\bar{\rho }}{\text{d}z}<0$ in $r_1\le r\le r_2$, then write

$k_{\min }^2=\frac{{\displaystyle {\int }_{r_1}^{r_2}}r{\left(\frac{\text{d}{\xi }_r}{\text{d}r}+\frac{{\xi }_r}{r}\right)}^2\text{d}r}{{\displaystyle {\int }_{r_1}^{r_2}}r{\xi }_r^2\Phi \left(r\right)\text{d}r}\left(-\frac{\text{d}\bar{\rho }}{\text{d}z}\right),$ (4)

the modes with the wave number $k\ge k_{\min }$ are stable, and when $0\le k<k_{\min }$ the modes are unstable;

3) When the rotation effect always $\Phi \left(r\right)<0$ the axisymmetric perturbations are unstable for the wave number $k>k_{\min }$, and be stable when $0\le k\le k_{\min }$ ;

· As to the stationary state ( $\omega =0$, which named as a marginal state), the critical value for $\frac{\text{d}\bar{\rho }}{\text{d}z}$ is determined by a variation problem,

${\frac{\text{d}\bar{\rho }}{\text{d}z}|}_{critical}=\underset{{\xi }_r}{\min }\left\{-\frac{{\displaystyle {\int }_{r_1}^{r_2}}{k}^2\Phi \left(r\right)r{\left|{\xi }_r\right|}^2\text{d}r}{{\displaystyle {\int }_{r_1}^{r_2}}r{\left|D_{*}{\xi }_r\right|}^2\text{d}r}\right\};$

· On the contrary to the case without rotation, the principle of the exchange of stabilities is invalid.

Remark 1.1 For the case $\frac{\text{d}\bar{\rho }}{\text{d}z}\equiv 0$, then (3) is the well-known Rayleigh’s criteria [1]. In this theorem, we conclude that when $\frac{\text{d}\bar{\rho }}{\text{d}z}\ne +\infty$, the variation of thermal shall transfer the unstable modes. In other words, when $\frac{\text{d}\bar{\rho }}{\text{d}z}<0$, which implies hot fluid under the cold fluid, for the static fluid, the buoyancy tends to overturn the fluid. In terms of the rotation fluid, the thermal convection only affects the low-wave numbers modes, whether on the case $\Phi \left(r\right)>0$ or $\Phi \left(r\right)<0$.

Theorem 1.2 The necessary and sufficient condition to be stable of the linearized system of (2) with non-axisymmetric perturbation is

1) For the special situation: $\Phi \left(r\right)=\frac{\text{d}\bar{\rho }}{\text{d}z}=const.>0$, then system is linear stable. Meanwhile, when $\Phi \left(r\right)=\frac{\text{d}\bar{\rho }}{\text{d}z}=const.<0$, the linearized system is unstable;

2) When $\Phi \left(r\right)\ne \frac{\text{d}\bar{\rho }}{\text{d}z}$, the necessary and sufficient condition to be stable of the linearized system of (2) with non-axisymmetric perturbation is

$k^2=\frac{{\displaystyle {\int }_{r_1}^{r_2}}r\left[\left(\Phi -s^2\right){\xi }_r^2-\left(m^2{P}^2\right)/\left(r^2{s}^2\right)\right]\text{d}r}{{\displaystyle {\int }_{r_1}^{r_2}}r{P}^2/\left(s^2-\frac{\text{d}\bar{\rho }}{\text{d}z}\right)\text{d}r}>0.$ (5)

To be more precisly, our results state:

· The necessary condition for the linear stable:

1) When $\Phi \left(r\right)>s^2={\left(a+m\Omega \right)}^2$, then the necessary condition for the linear stable are $\frac{\text{d}\bar{\rho }}{\text{d}z}<{\left(a+m\Omega \right)}^2$ or $\frac{\text{d}\bar{\rho }}{\text{d}z}>\frac{k^2{r}^2+m^2}{m^2}{\left(a+m\Omega \right)}^2$ ;

2) When $\Phi \left(r\right)<s^2={\left(a+m\Omega \right)}^2$, $\frac{\text{d}\bar{\rho }}{\text{d}z}$ must satisfy ${\left(a+m\Omega \right)}^2<\frac{\text{d}\bar{\rho }}{\text{d}z}<\frac{k^2{r}^2+m^2}{m^2}{\left(a+m\Omega \right)}^2$.

· The sufficient condition for the linear stable and unstable:

1) When $\Phi \left(r\right)<s^2$ and $\frac{\text{d}\bar{\rho }}{\text{d}z}<s^2$, the linearized system is unstable;

2) When $\Phi \left(r\right)<s^2<\frac{\text{d}\bar{\rho }}{\text{d}z}$, the linearized system is stable.

The paper is organized as follows. In Section 3, we shall prove Theorem 1.1, and then theorem 1.2 shall be proved in Section 4.

2. The Perturbation Equations and the Basic State

We’re dealing with the system (2) in a coaxial cylinders area (see Figure 1), it’s more convenient to study the system (2) in cylindrical coordinates. Let

$e_r=\left(\frac{x}{r},\frac{y}{r},0\right),e_{\theta }=\left(-\frac{y}{r},\frac{x}{r},0\right),e_z=\left(0,0,1\right),$

with $r=\sqrt{x^2+y^2}$, and write

$U=u_r\left(r\mathrm{,}\theta \mathrm{,}z\mathrm{,}t\right)e_r+u_{\theta }\left(r\mathrm{,}\theta \mathrm{,}z\mathrm{,}t\right)e_{\theta }+u_z\left(r\mathrm{,}\theta \mathrm{,}z\mathrm{,}t\right)e_z\mathrm{,}$ (6)

we obtain

$(\begin{array}{l}{\partial }_t\rho +u_r{\partial }_r\rho +\frac{1}{r}{u}_{\theta }{\partial }_{\theta }\rho +u_z{\partial }_z\rho =0,\hfill \\ {\partial }_t{u}_r+u_r{\partial }_r{u}_r+\frac{1}{r}{u}_{\theta }{\partial }_{\theta }{u}_r+u_z{\partial }_z{u}_r-\frac{u_{\theta }^2}{r}=-{\partial }_{r}P,\hfill \\ {\partial }_t{u}_{\theta }+u_r{\partial }_r{u}_{\theta }+\frac{1}{r}{u}_{\theta }{\partial }_{\theta }{u}_{\theta }+u_z{\partial }_z{u}_{\theta }+\frac{u_{\theta }{u}_r}{r}=-\frac{1}{r}{\partial }_{\theta }P,\hfill \\ {\partial }_t{u}_z+u_r{\partial }_r{u}_z+\frac{1}{r}{u}_{\theta }{\partial }_{\theta }{u}_z+u_z{\partial }_z{u}_z=-{\partial }_{z}P+\rho ,\hfill \end{array}$ (7)

with

${\partial }_r{u}_r+\frac{u_r}{r}+\frac{1}{r}{\partial }_{\theta }{u}_{\theta }+{\partial }_z{u}_z=0.$ (8)

By writing $u\cdot \nabla =u_r\frac{\partial }{{\partial }_r}+\frac{u_{\theta }}{r}\frac{\partial }{{\partial }_{\theta }}+u_z\frac{\partial }{{\partial }_z}$, we investigate the following system:

$(\begin{array}{l}{\partial }_t\rho +\left(u\cdot \nabla \right)\rho =0,\hfill \\ {\partial }_t{u}_r+\left(u\cdot \nabla \right)u_r-\frac{u_{\theta }^2}{r}=-{\partial }_{r}P,\hfill \\ {\partial }_t{u}_{\theta }+\left(u\cdot \nabla \right)u_{\theta }+\frac{u_{\theta }{u}_r}{r}=-\frac{1}{r}{\partial }_{\theta }P,\hfill \\ {\partial }_t{u}_z+\left(u\cdot \nabla \right)u_z=-{\partial }_{z}P+\rho ,\hfill \\ \nabla \cdot u=0,U\cdot n=0,\hfill \end{array}$ (9)

where $n$ denotes the normal exterior vector.

Let $\left(u_r\mathrm{,}{u}_{\theta }\mathrm{,}{u}_z\right)$ are the velocity components in the cylindrical coordinates $\left(r\mathrm{,}\theta \mathrm{,}z\right)$ respectively, the velocity boundary conditions are

$(\begin{array}{l}{u_r|}_{r=r_1,r_2}=0,\\ {u_{\theta }|}_{r=r_1}={\Omega }_1{r}_1,{u_{\theta }|}_{r=r_2}={\Omega }_2{r}_2,\\ {u_z|}_{z=0,1}=0,\end{array}$ (10)

where $r=r_1$, $r=r_2$ are the radii of the two cylinders. ${\Omega }_1$ and ${\Omega }_2$ are the constant rotation velocities at $r=r_1$ and $r=r_2$ respectively.

In this paper, we shall study the stability of the following stationary Couette flow

${\bar{u}}_r={\bar{u}}_z=0,{\bar{u}}_{\theta }=V\left(r\right)=r\Omega \left(r\right),$ (11)

where $V\left(r\right)$ is an arbitrary smooth function of r.

In the experiments the basic temperature state is, a priori, a time-dependent state because the initial temperature gradient is not maintained by the boundary conditions. However, if we consider the time scale smaller than the typical diffusion time over the height of the cylinders, this gradient can be considered as constant in time. Therefore, we choose as time scale small enough, which allows us to write the basic state for the temperature as time-independent. In this paper, we consider the variation of background thermal as the following case:

$\bar{\rho }={\bar{\rho }}_{\left(z\right)}.$ (12)

According to the (11) and (12), the basic pressure distribution is determined by the system (9). It takes the form as

$\bar{P}={\displaystyle \int }\frac{V{\left(r\right)}^2}{r}\text{d}r+{\displaystyle \int {\bar{\rho }}_{\left(z\right)}\text{d}z}+const\mathrm{.}$ (13)

Consider an infinitesimal perturbation of the system (11)-(13), we write the perturbed state as

$\rho +\bar{\rho }\mathrm{,}{u}_r+{\bar{u}}_r\mathrm{,}{u}_{\theta }+{\bar{u}}_{\theta }\mathrm{,}{u}_z+{\bar{u}}_z\mathrm{,}P+\bar{P}\mathrm{.}$ (14)

Substitute (14) into (9), we get the system for the perturbations as

$(\begin{array}{l}{\partial }_t\rho +u_r{\partial }_r\rho +\frac{1}{r}\left(V+u_{\theta }\right){\partial }_{\theta }\rho +u_z{\partial }_z\rho +u_z{\partial }_z{\bar{\rho }}_{\left(z\right)}=0,\hfill \\ {\partial }_t{u}_r+u_r{\partial }_r{u}_r+\frac{1}{r}\left(V+u_{\theta }\right){\partial }_{\theta }{u}_r+u_z{\partial }_z{u}_r-\frac{{\left(V+u_{\theta }\right)}^2}{r}+{\partial }_{r}P+\frac{V{\left(r\right)}^2}{r}=0,\hfill \\ {\partial }_t{u}_{\theta }+u_r{\partial }_r\left(V+u_{\theta }\right)+\frac{1}{r}\left(V+u_{\theta }\right){\partial }_{\theta }{u}_{\theta }+u_z{\partial }_z{u}_{\theta }+\frac{\left(V+u_{\theta }\right)u_r}{r}+\frac{1}{r}{\partial }_{\theta }P=0,\hfill \\ {\partial }_t{u}_z+u_r{\partial }_r{u}_z+\frac{1}{r}\left(V+u_{\theta }\right){\partial }_{\theta }{u}_z+u_z{\partial }_z{u}_z+{\partial }_{z}P=\rho ,\hfill \\ {\partial }_r{u}_r+\frac{u_r}{r}+\frac{1}{r}{\partial }_{\theta }\left(V+u_{\theta }\right)+{\partial }_z{u}_z=0.\hfill \end{array}$ (15)

The linearized equations governing these perturbations system (15) is

$(\begin{array}{l}{\partial }_t\rho +\frac{V}{r}{\partial }_{\theta }\rho +u_z{\partial }_z{\bar{\rho }}_{\left(z\right)}=0,\hfill \\ {\partial }_t{u}_r+\frac{V}{r}{\partial }_{\theta }{u}_r-\frac{2V{u}_{\theta }}{r}+{\partial }_{r}P=0,\hfill \\ {\partial }_t{u}_{\theta }+\frac{V}{r}{\partial }_{\theta }{u}_{\theta }+\left(\frac{\text{d}V}{\text{d}r}+\frac{V}{r}\right)u_r+\frac{1}{r}{\partial }_{\theta }P=0,\hfill \\ {\partial }_t{u}_z+\frac{V}{r}{\partial }_{\theta }{u}_z+{\partial }_{z}P=\rho ,\hfill \\ {\partial }_r{u}_r+\frac{u_r}{r}+\frac{1}{r}{\partial }_{\theta }{u}_{\theta }+{\partial }_z{u}_z=0.\hfill \end{array}$ (16)

We analyze the disturbance by using the Normal-modes method. It is natural to write that the various quantities describing the perturbation have a $\left(t\mathrm{,}r\mathrm{,}\theta \mathrm{,}z\right)$ -dependence as

$(\begin{array}{l}\rho =\rho \left(r\right){\text{e}}^{i\left(kz+m\theta -\omega t\right)},\hfill \\ u_r=u_r\left(r\right){\text{e}}^{i\left(kz+m\theta -\omega t\right)},\hfill \\ u_{\theta }=u_{\theta }\left(r\right){\text{e}}^{i\left(kz+m\theta -\omega t\right)},\hfill \\ u_z=u_z\left(r\right){\text{e}}^{i\left(kz+m\theta -\omega t\right)},\hfill \\ P=P\left(r\right){\text{e}}^{i\left(kz+m\theta -\omega t\right)},\hfill \end{array}$ (17)

where $\omega$ is a constant (which may be comlex). m and k are positive integers, which are the oscillation frequency and the wave number of the disturbance in the z-direction respectively. $u_r\left(r\right)$, $u_{\theta }\left(r\right)$, $u_z\left(r\right)$ and $P\left(r\right)$ denote the amplitudes of the respective perturbations. From (16) and (17), then we get

$is\rho +u_z\frac{\text{d}\bar{\rho }}{\text{d}z}=\mathrm{0,}$ (18)

$is{u}_r-2\Omega u_{\theta }=-\frac{\text{d}P}{\text{d}r}\mathrm{,}$ (19)

$is{u}_{\theta }+\left(\Omega +\frac{\text{d}\left(r\Omega \right)}{\text{d}r}\right)u_r=-\frac{imP}{r}\mathrm{,}$ (20)

$is{u}_z=-ikP+\rho \mathrm{,}$ (21)

$\frac{\text{d}{u}_r}{\text{d}r}+\frac{u_r}{r}+\frac{im}{r}{u}_{\theta }+ik{u}_z=0.$ (22)

where we used the notation $s=-\omega +m\Omega$, which is named as the Doppler-shifted frequency. Physically, $\Omega =\frac{V\left(r\right)}{r}$ be the angular velocity and $\Omega +\text{d}\left(r\Omega \right)/\text{d}r=\left(1/r\right)\text{d}\left(r^2\Omega \right)/\text{d}r$ be the axial vorticity of the base flow.

To simplify the system (18)-(22), we consider the Lagrangian displacement ${\xi }_r$, ${\xi }_{\theta }$ and ${\xi }_z$ which describe the displacement of a fluid element in the perturbed flow relative the location at time t in the unperturbed flow. Thus we have variation of the velocity as

$\delta u=\frac{\partial \xi }{\partial t}+U\cdot \nabla \xi \mathrm{,}$ (23)

where U is the velocity field of the unperturbed flow. By using the Taylor’s formula, we also have

$\delta u=u+\xi \cdot \nabla u\mathrm{.}$ (24)

Combining (23) and (24), we get

$u=\frac{\partial \xi }{\partial t}+U\cdot \nabla \xi -\xi \cdot \nabla U\mathrm{.}$ (25)

Now, we proceed directly to the equations of the normal modes by letting $\xi =\left({\xi }_r\mathrm{,}{\xi }_{\theta }\mathrm{,}{\xi }_z\right){\text{e}}^{i\left(kz+m\theta -\omega t\right)}$, then there holds

$u_r=is{\xi }_r,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{u}_{\theta }=is{\xi }_{\theta }-r\frac{\text{d}\Omega }{\text{d}r}{\xi }_r,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{u}_z=is{\xi }_z.$ (26)

By using the solenoidal character of u and (26), we have

$s\rho +s\frac{\text{d}\bar{\rho }}{\text{d}z}{\xi }_z=\mathrm{0,}$ (27)

$s^2{\xi }_{\theta }-2i\Omega s{\xi }_r=\frac{im}{r}P\mathrm{,}$ (28)

$-s^2{\xi }_r-2\Omega \left(is{\xi }_{\theta }-r\frac{\text{d}\Omega }{\text{d}r}{\xi }_r\right)=-\frac{\text{d}P}{\text{d}r}\mathrm{,}$ (29)

$s^2{\xi }_z=ikP-\rho \mathrm{,}$ (30)

$\frac{\text{d}{\xi }_r}{\text{d}r}+\frac{{\xi }_r}{r}+\frac{im}{r}{\xi }_{\theta }+ik{\xi }_z=0.$ (31)

Eliminating ${\xi }_{\theta }$ between Equations (28) and (29), we have

$\left(s^2-2r\Omega \frac{\text{d}\Omega }{\text{d}r}\right){\xi }_r+\frac{2i\Omega }{s}\left(2i\Omega s{\xi }_r+\frac{im}{r}P\right)=\frac{\text{d}P}{\text{d}r}\mathrm{.}$ (32)

Rearranging the terms in this equation, we obtain

$\left[s^2-\Phi \left(r\right)\right]{\xi }_r=\frac{\text{d}P}{\text{d}r}+\frac{2m\Omega }{sr}P\mathrm{,}$ (33)

where

$\Phi \left(r\right)=2r\Omega \frac{\text{d}\Omega }{\text{d}r}+4{\Omega }^2=\frac{2\Omega }{r}\frac{\text{d}}{\text{d}r}\left(r^2\Omega \right)$ (34)

is the Rayleigh discriminant. Similarly, we also get

$\frac{\text{d}{\xi }_r}{\text{d}r}+\frac{{\xi }_r}{r}-\frac{2m\Omega }{rs}{\xi }_r=\left(\frac{m^2}{r^2{s}^2}+\frac{k^2}{s^2-\frac{\text{d}\bar{\rho }}{\text{d}z}}\right)P\mathrm{.}$ (35)

Due to the fluid is confined between two coaxial cylinders of radii $r_1$ and $r_2$, we must require that the radial components of the velocity vanish for these values of r. Thus, we consider the Equations (33) and (35) with the boundary conditions

$(\begin{array}{l}{{\xi }_r|}_{r=r_1,r_2}=0,\\ {\left(\frac{\text{d}P}{\text{d}r}+\frac{2m\Omega }{sr}P\right)|}_{r=r_1,r_2}=0.\end{array}$ (36)

3. Linear Stability for Axisymmetric Perturbations

In this section, we shall prove Theorem 1.1, the results for the axisymmetric perturbation case. Namely in (17), we take $m=0$ (which implies $s=-\omega +m\Omega =-\omega$ ), then Equations (33) and (35) read as

$\left({\omega }^2-\Phi \left(r\right)\right){\xi }_r=\frac{\text{d}P}{\text{d}r}\mathrm{,}$ (37)

$\frac{\text{d}{\xi }_r}{\text{d}r}+\frac{{\xi }_r}{r}=\frac{k^2}{{\omega }^2-\frac{\text{d}\bar{\rho }}{\text{d}z}}P.$ (38)

Eliminating P between these equations, we obtain

$\left({\omega }^2-\frac{\text{d}\bar{\rho }}{\text{d}z}\right)\left(D{D}_{*}-k^2\right){\xi }_r=-k^2\left(\Phi \left(r\right)-\frac{\text{d}\bar{\rho }}{\text{d}z}\right){\xi }_r.$ (39)

Write $\mathcal{M}={\omega }^2-\frac{\text{d}\bar{\rho }}{\text{d}z}$, we get

$\mathcal{M}\left(D{D}_{\mathrm{<mo>\; *\; </mo>}}-k^2\right){\xi }_r=-k^2\left(\Phi \left(r\right)-\frac{\text{d}\bar{\rho }}{\text{d}z}\right){\xi }_r$ (40)

where we used the notations

$D=\text{d}/\text{d}r\text{and}{D}_{\mathrm{<mo>\; *\; </mo>}}=\text{d}/\text{d}r+1/r\mathrm{.}$

Equation (40) together with the boundary conditions (36) constitute a characteristic value problem in terms of $\mathcal{M}$. Now let’s prove the Theorem 1.1.

We shall first show that the principle of the exchange of stabilities is invalid.

3.1. Exchange of Stability Is Invalid

Proof. Testing Equation (40) by $r{\xi }_r^{\mathrm{<mo>\; *\; </mo>}}$ (where ${\xi }_r^{\mathrm{<mo>\; *\; </mo>}}$ is the conjugate function of ${\xi }_r$ ), we obtain

$\mathcal{M}{\displaystyle {\int }_{r_1}^{r_2}}D{D}_{\mathrm{<mo>\; *\; </mo>}}{\xi }_r\left(r{\xi }_r^{\mathrm{<mo>\; *\; </mo>}}\right)\text{d}r-\mathcal{M}{\displaystyle {\int }_{r_1}^{r_2}}{k}^2{\xi }_r\left(r{\xi }_r^{\mathrm{<mo>\; *\; </mo>}}\right)\text{d}r=-k^2{\displaystyle {\int }_{r_1}^{r_2}}\left(\Phi \left(r\right)-\frac{\text{d}\bar{\rho }}{\text{d}z}\right){\xi }_r\left(r{\xi }_r^{\mathrm{<mo>\; *\; </mo>}}\right)\text{d}r$ (41)

Integrating by parts, from the left hand of side, we get

$\begin{array}{l}\mathcal{M}{\displaystyle {\int }_{r_1}^{r_2}}D{D}_{*}{\xi }_r\left(r{\xi }_r^{*}\right)\text{d}r-\mathcal{M}{\displaystyle {\int }_{r_1}^{r_2}}{k}^2{\xi }_r\left(r{\xi }_r^{*}\right)\text{d}r\\ =\mathcal{M}\left[D_{*}{\xi }_r\cdot {\left(r{\xi }_r^{*}\right)|}_{r=r_1,r_2}-{\displaystyle {\int }_{r_1}^{r_2}}{D}_{*}{\xi }_{r}D\left(r{\xi }_r^{*}\right)\text{d}r\right]-\mathcal{M}{\displaystyle {\int }_{r_1}^{r_2}}r{k}^2{\left|{\xi }_r\right|}^2\text{d}r\\ =-\mathcal{M}{\displaystyle {\int }_{r_1}^{r_2}}r\left({\left|D_{*}{\xi }_r\right|}^2+k^2{\left|{\xi }_r\right|}^2\right)\text{d}r.\end{array}$ (42)

Similarly, from the right hand side of the Equation (41) we get

$-k^2{\displaystyle {\int }_{r_1}^{r_2}}\left(\Phi \left(r\right)-\frac{\text{d}\bar{\rho }}{\text{d}z}\right){\xi }_r\left(r{\xi }_r^{*}\right)=-k^2{\displaystyle {\int }_{r_1}^{r_2}}\left[r\Phi \left(r\right){\left|{\xi }_r\right|}^2-\frac{\text{d}\bar{\rho }}{\text{d}z}r{\left|{\xi }_r\right|}^2\right]\text{d}r.$ (43)

So we have

$\mathcal{M}=\frac{k^2{\displaystyle {\int }_{r_1}^{r_2}}r\left(\Phi \left(r\right)-\frac{\text{d}\bar{\rho }}{\text{d}z}\right){\left|{\xi }_r\right|}^2\text{d}r}{{\displaystyle {\int }_{r_1}^{r_2}}r\left({\left|D_{\mathrm{<mo>\; *\; </mo>}}{\xi }_r\right|}^2+k^2{\left|{\xi }_r\right|}^2\right)\text{d}r}\mathrm{.}$ (44)

From (44), we conclude that the real and the imaginary parts must vanish separately. Noting that k, $\Phi \left(r\right)$ and $\frac{\text{d}\bar{\rho }}{\text{d}z}$ are all real, we get

$im\left(\mathcal{M}\right)=0.$ (45)

Recalling that $\mathcal{M}={\omega }^2-\frac{\text{d}\bar{\rho }}{\text{d}z}$, then (45) implies that $\omega$ be a real number or pure imaginary number. Let $\omega =a+bi$ then $a\cdot b=0$.

Now we prove the exchange of stability is invalid by using the proof of contradiction. Let $b=0$, by using (44), then

$a^2=\frac{k^2{\displaystyle {\int }_{r_1}^{r_2}}r\left(\Phi \left(r\right)-\frac{\text{d}\bar{\rho }}{\text{d}z}\right){\left|{\xi }_r\right|}^2\text{d}r}{{\displaystyle {\int }_{r_1}^{r_2}}r\left({\left|D_{\mathrm{<mo>\; *\; </mo>}}{\xi }_r\right|}^2+k^2{\left|{\xi }_r\right|}^2\right)\text{d}r}+\frac{\text{d}\bar{\rho }}{\text{d}z}\mathrm{,}k\in \mathbb{N}\mathrm{.}$ (46)

Since the Equation (46) holds for $k\in \mathbb{N}$, we conclude that $a\cancel{\equiv }0$, which implies the principle of exchange the stability is invalid.

3.2. The Critical Value of $\frac{d\bar{\rho }}{dz}$ at the Marginal State

The marginal state illustrates that the transition from stability to instability. In the following, we present the critical value for the $\frac{\text{d}\bar{\rho }}{\text{d}z}$ at the marginal state.

The equations governing the marginal state as (by setting $\omega =0$ )

$(\begin{array}{l}-\Phi \left(r\right){\xi }_r=\frac{\text{d}P}{\text{d}r},\\ D_{*}{\xi }_r=\frac{k^2}{-\frac{\text{d}\bar{\rho }}{\text{d}z}}P.\end{array}$ (47)

Eliminating P between these equations, and denoting $R=-\frac{1}{\frac{\text{d}\bar{\rho }}{\text{d}z}}$, we obtain

$(\begin{array}{l}D{D}_{*}{\xi }_r=-k^{2}R\Phi \left(r\right){\xi }_r,\\ {{\xi }_r|}_{r=r_1,r_2}=0.\end{array}$ (48)

We shall show that the critical value for R in (48) is a minimum for a variation problem. By testing $r{\xi }_r$ to the Equation (48), we get

${\displaystyle {\int }_{r_1}^{r_2}}D{D}_{*}{\xi }_r\left(r{\xi }_r\right)\text{d}r=-{\displaystyle {\int }_{r_1}^{r_2}}{k}^{2}R\Phi \left(r\right)r{\left|{\xi }_r\right|}^2\text{d}r.$ (49)

Integration by parts we have

$\begin{array}{c}{\displaystyle {\int }_{r_1}^{r_2}}D{D}_{*}{\xi }_r\left(r{\xi }_r\right)\text{d}r={D_{*}{\xi }_r\left(r{\xi }_r\right)|}_{r=r_1,r_2}-{\displaystyle {\int }_{r_1}^{r_2}}{D}_{*}{\xi }_{r}D\left(r{\xi }_r\right)\text{d}r\\ =-{\displaystyle {\int }_{r_1}^{r_2}}r{\left|D_{*}{\xi }_r\right|}^2\text{d}r.\end{array}$ (50)

From (49) and (50) we obtain

$R=\frac{{\displaystyle {\int }_{r_1}^{r_2}}r{\left|D_{*}{\xi }_r\right|}^2\text{d}r}{k^2{\displaystyle {\int }_{r_1}^{r_2}}\Phi \left(r\right)r{\left|{\xi }_r\right|}^2\text{d}r}:=\frac{I_1}{I_2}.$ (51)

In (51), the characteristic value of $R_{\min }$ will be a minimumin terms of the characteristic functions of a variation problem. To verify this fact, we denote $\delta R$ be the variation in R when ${\xi }_r$ is subjected to a small variation $\delta {\xi }_r$ which is also compatible with the boundary conditions on ${\xi }_r$.

According to the Equation (51), we obtain

$\delta R=\frac{1}{I_2}\left(\delta I_1-\frac{I_1}{I_2}\delta I_2\right)=\frac{1}{I_2}\left(\delta I_1-R\delta I_2\right)\mathrm{,}$ (52)

where

$\begin{array}{l}\delta I_1={\frac{\text{d}}{\text{d}\epsilon }\left({\displaystyle {\int }_{r_1}^{r_2}}r\left(\frac{\text{d}}{\text{d}r}+\frac{1}{r}\right)\left({\xi }_r+\epsilon \delta {\xi }_r\right)\left(\frac{\text{d}}{\text{d}r}+\frac{1}{r}\right)\left({\xi }_r+\epsilon \delta {\xi }_r\right)\right)|}_{\epsilon =0}\\ ={\displaystyle {\int }_{r_1}^{r_2}}r\left[\left(\frac{\text{d}}{\text{d}r}\delta {\xi }_r+\frac{1}{r}\delta {\xi }_r\right)\left(\frac{\text{d}}{\text{d}r}{\xi }_r+\frac{1}{r}{\xi }_r\right)+\left(\frac{\text{d}}{\text{d}r}{\xi }_r+\frac{1}{r}{\xi }_r\right)\left(\frac{\text{d}}{\text{d}r}\delta {\xi }_r+\frac{1}{r}\delta {\xi }_r\right)\right]\text{d}r\\ =2{\displaystyle {\int }_{r_1}^{r_2}}r\left(\frac{\text{d}}{\text{d}r}\delta {\xi }_r+\frac{1}{r}\delta {\xi }_r\right)\left(D_{*}{\xi }_r\right)\text{d}r,\end{array}$ (53)

and

$\delta I_2={\frac{\text{d}}{\text{d}\epsilon }\left({\displaystyle {\int }_{r_1}^{r_2}}\Phi \left(r\right)r\left({\xi }_r+\epsilon \delta {\xi }_r\right)\left({\xi }_r+\epsilon \delta {\xi }_r\right)\text{d}r\right)|}_{\epsilon =0}=2{\displaystyle {\int }_{r_1}^{r_2}}\Phi \left(r\right)r\delta {\xi }_r{\xi }_r\text{d}r.$ (54)

Integrations by parts, we further get

$\delta I_1={2\left(r\delta {\xi }_r{D}_{*}{\xi }_r\right)|}_{r=r_1,r_2}-2{\displaystyle {\int }_{r_1}^{r_2}}D{D}_{*}{\xi }_r\left(r\delta {\xi }_r\right)\text{d}r=-2{\displaystyle {\int }_{r_1}^{r_2}}D{D}_{*}{\xi }_r\left(r\delta {\xi }_r\right)\text{d}r.$ (55)

Thus

$\delta R=-\frac{2}{I_2}{\displaystyle {\int }_{r_1}^{r_2}}r\delta {\xi }_r\left[D{D}_{*}{\xi }_r+k^{2}R\Phi \left(r\right){\xi }_r\right]\text{d}r.$ (56)

From (56), it follows that $\delta R=0$ if and only if $D{D}_{\mathrm{<mo>\; *\; </mo>}}{\xi }_r+k^{2}R\Phi \left(r\right){\xi }_r=0$, which comes from the Equation (48).

3.3. The Linear Stability for Axsymmetric Perturbation

Proof. Recalling (45), for simplicity we write $Re\left(\mathcal{M}\right)=\mathcal{M}$. Next, we rewrite Equation (41) as

$\mathcal{M}{I}_1-k^2{I}_2+\frac{\text{d}\bar{\rho }}{\text{d}z}{k}^2{I}_3=\mathrm{0,}$ (57)

with

$\begin{array}{l}{I}_1={\displaystyle {\int }_{r_1}^{r_2}}r\left[{\left|\left(\frac{\text{d}}{\text{d}r}+\frac{1}{r}\right){\xi }_r\right|}^2+k^2{\left|{\xi }_r\right|}^2\right]\text{d}r\mathrm{,}\\ I_2={\displaystyle {\int }_{r_1}^{r_2}}r\Phi \left(r\right){\left|{\xi }_r\right|}^2\text{d}r\mathrm{,}\\ I_3={\displaystyle {\int }_{r_1}^{r_2}}r{\left|{\xi }_r\right|}^2\text{d}r\mathrm{.}\end{array}$ (58)

From (57), we get

$\mathcal{M}={\omega }^2-\frac{\text{d}\bar{\rho }}{\text{d}z}=\frac{k^2{I}_2-\frac{\text{d}\bar{\rho }}{\text{d}z}{k}^2{I}_3}{I_1}.$ (59)

To keep the linearized system stable, the necessary and sufficient condition is that $\omega$ be real. That is to say,

${\omega }^2=\frac{k^2{I}_2-\frac{\text{d}\bar{\rho }}{\text{d}z}{k}^2{I}_3}{I_1}+\frac{\text{d}\bar{\rho }}{\text{d}z}>0.$ (60)

It is equivalent to study

$k^2{I}_2-\frac{\text{d}\bar{\rho }}{\text{d}z}{k}^2{I}_3+\frac{\text{d}\bar{\rho }}{\text{d}z}{I}_1>0.$ (61)

Recalling (58), the inequality (61) is equivalent to the following one:

$\frac{\text{d}\bar{\rho }}{\text{d}z}{\displaystyle {\int }_{r_1}^{r_2}}r{\left|D_{*}{\xi }_r\right|}^2\text{d}r+{\displaystyle {\int }_{r_1}^{r_2}}{k}^{2}r{\xi }_r^2\Phi \left(r\right)\text{d}r>0.$ (62)

From (62), we get

Case 1. $\frac{\text{d}\bar{\rho }}{\text{d}z}=0$. In this case, the equilibrium state $\bar{\rho }$ is a constant. From the Equation (62), it is sufficient to study

${\displaystyle {\int }_{r_1}^{r_2}}r{\left|{\xi }_r\right|}^2\Phi \left(r\right)\text{d}r>0.$ (63)

According to the inequality (63), it is apparent that the linearized systems (33)-(36) are stable when $\Phi \left(r\right)$ is positive, and which is unstable when $\Phi \left(r\right)$ is negative. And if $\Phi \left(r\right)$ changes sign in the interval $\left(r_1\mathrm{,}{r}_2\right)$, Such case is a locally unstable situation.

This case is also the well-known Rayleigh’s criterion.

Case 2. $\frac{\text{d}\bar{\rho }}{\text{d}z}>0$. From the Equation (62), we obtain

1) If $\Phi \left(r\right)>0$, the inequality (62) always holds. In this case, the fluid is stable.

2) If $\Phi \left(r\right)<0$, we write

$k_{\min }^2=-\frac{{\displaystyle {\int }_{r_1}^{r_2}}r{\left|D_{*}{\xi }_r\right|}^2\text{d}r}{{\displaystyle {\int }_{r_1}^{r_2}}r{\xi }_r^2\Phi \left(r\right)\text{d}r}\left(\frac{\text{d}\bar{\rho }}{\text{d}z}\right),$ (64)

then for the modes $k\ge k_{\min }$, they all be unstable. For the modes $k<k_{\min }$, they are stable. Compare with the results of Rayleigh’s situation (see for example case 1), the condition $\frac{\text{d}\bar{\rho }}{\text{d}z}>0$, which means that cold fluid is under hot fluid, the buoyancy tends to stabilize the fluid although the rotation tends to turn over the fluid. In this case, the positive temperature gradient effect can stabilize the low modes, but the rotation effect is the dominant role to cause the instability.

Case 3. $\frac{\text{d}\bar{\rho }}{\text{d}z}<0$. Physically, under this condition, it means that hot fluid under the cold fluid.

1) If $\Phi \left(r\right)<0$, the fluid must be unstable and convection will occur. In this case, both the rotation and the buoyancy tend to overturn the flow.

2) If $\Phi \left(r\right)>0$, when $k>k_{\min }$ the fluid is stable and when $0\le k<k_{\min }$, it is unstable. Compare with Case 2-2), we conclude that the thermal affects only the small wave numbers (whatever the stabilize or destabilize the fluid).

4. Non-Axisymmetric Linear Stability Analysis

In this section, we study Theorem 1.2, the non-axisymmetric perturbation situation.

Proof. We consider non-axisymmetric perturbations, namely, $m\ne 0$. Noting that $s=-\omega +m\Omega$, we shall analyze the following cases.

Case 1. $s^2=0$. In this case, it is equivalent to

$(\begin{array}{l}\left(a+m\Omega \right)b=0,\\ {\left(a+m\Omega \right)}^2-b^2=0.\end{array}$ (65)

From which, we conclude that $b=0$, When $b=0$. In this situation, the perturbation is stable. Actually, from (26)-(31), this case means the trivial case $u_r=u_{\theta }=u_z=0$.

Case 2. $s^2=\frac{\text{d}\bar{\rho }}{\text{d}z}\ne 0$. In this case, it is equivalent to

$(\begin{array}{l}\left(a+m\Omega \right)b=0,\\ {\left(a+m\Omega \right)}^2-b^2=\frac{\text{d}\bar{\rho }}{\text{d}z}.\end{array}$ (66)

If $\frac{\text{d}\bar{\rho }}{\text{d}z}\ge 0$, then recalling (4.2)₁ we get $b=0$ which implies the stable situation. When $\frac{\text{d}\bar{\rho }}{\text{d}z}<0$, then recalling (4.2) we conclude that $b\ne 0$, which means unstable case.

Under this situation, by using (26)-(31) again, we conclude that $\Phi \left(r\right)\ne \frac{\text{d}\bar{\rho }}{\text{d}z}\left(z\right)$ or $\Phi \left(r\right)=\frac{\text{d}\bar{\rho }}{\text{d}z}\left(z\right)=const$. In the first case, it is corresponding to the trivial case $u_r=u_{\theta }=u_z=0$. For the second case, then $\Phi \left(r\right)=\frac{\text{d}\bar{\rho }}{\text{d}z}\left(z\right)=const>0$ it is stable and $\Phi \left(r\right)=\frac{\text{d}\bar{\rho }}{\text{d}z}\left(z\right)=const<0$ is unstable. In this situation, both the rotation and the thermal convection play the stabilization or destabilize role simutanously.

Case 3. $s^2\ne \frac{\text{d}\bar{\rho }}{\text{d}z}$ and $s^2\ne 0$. In this case, we write

$C\left(r\right)=\frac{m^2}{r^2{s}^2}+\frac{k^2}{s^2-\frac{\text{d}\bar{\rho }}{\text{d}z}}\mathrm{.}$ (67)

Since $m\ne 0$, which implies that $C\left(r\right)\ne 0$. From (35), we obtain

$P=\frac{1}{rC\left(r\right)}\left[\frac{\text{d}}{\text{d}r}\left(r{\xi }_r\right)-\frac{2m\Omega }{rs}r{\xi }_r\right]\mathrm{,}$ (68)

and

$\frac{\text{d}P}{\text{d}r}=\frac{\text{d}}{\text{d}r}\left[\frac{1}{rC\left(r\right)}\frac{\text{d}}{\text{d}r}\left(r{\xi }_r\right)\right]-\left[\frac{\text{d}}{\text{d}r}\left(\frac{2m\Omega }{rC\left(r\right)rs}r{\xi }_r\right)+\frac{2m\Omega }{rC\left(r\right)rs}\frac{\text{d}}{\text{d}r}\left(r{\xi }_r\right)\right].$ (69)

Substitute Equations (68) and (69) into Equation (33) and eliminat P, we obtain

$\begin{array}{l}\frac{\text{d}}{\text{d}r}\left[\frac{1}{rC\left(r\right)}\frac{\text{d}}{\text{d}r}\left(r{\xi }_r\right)\right]-\left[\frac{\text{d}}{\text{d}r}\left(\frac{2m\Omega }{rC\left(r\right)rs}\right)\right]\left(r{\xi }_r\right)-\frac{4m^2{\Omega }^2}{r^{3}C\left(r\right)s^2}\left(r{\xi }_r\right)\\ =\frac{s^2-\Phi \left(r\right)}{r}\left(r{\xi }_r\right).\end{array}$ (70)

By testing $-r{\xi }_r^{\mathrm{<mo>\; *\; </mo>}}$ (the conjugate function of $-r{\xi }_r$ ) to Equation (70), we obtain

$\begin{array}{l}-{\displaystyle {\int }_{r_1}^{r_2}}\frac{\text{d}}{\text{d}r}\left[\frac{1}{rC\left(r\right)}\frac{\text{d}}{\text{d}r}\left(r{\xi }_r\right)\right]\left(r{\xi }_r^{*}\right)\text{d}r+{\displaystyle {\int }_{r_1}^{r_2}}\frac{\text{d}}{\text{d}r}\left(\frac{2m\Omega }{r^{2}C\left(r\right)s}\right){\left|r{\xi }_r\right|}^2\text{d}r\\ +{\displaystyle {\int }_{r_1}^{r_2}}\frac{4m^2{\Omega }^2}{r^{3}C{s}^2}{\left|r{\xi }_r\right|}^2\text{d}r=-{\displaystyle {\int }_{r_1}^{r_2}}\frac{s^2-\Phi \left(r\right)}{r}{\left|r{\xi }_r\right|}^2\text{d}r.\end{array}$ (71)

Integrating by parts and using the boundary condition, we estimate the first term of (71) as

$\begin{array}{l}-{\displaystyle {\int }_{r_1}^{r_2}}\frac{\text{d}}{\text{d}r}\left[\frac{1}{rC\left(r\right)}\frac{\text{d}}{\text{d}r}\left(r{\xi }_r\right)\right]\left(r{\xi }_r^{*}\right)\text{d}r\\ =-{\left(\frac{1}{rC\left(r\right)}\frac{\text{d}}{\text{d}r}r{\xi }_r^{*}\right)|}_{r=r_1,r_2}+{\displaystyle {\int }_{r_1}^{r_2}}\frac{1}{rC\left(r\right)}\frac{\text{d}}{\text{d}r}\left(r{\xi }_r\right)\frac{\text{d}}{\text{d}r}\left(r{\xi }_r^{*}\right)\\ ={\displaystyle {\int }_{r_1}^{r_2}}\frac{1}{rC\left(r\right)}{\left|\frac{\text{d}}{\text{d}r}\left(r{\xi }_r\right)\right|}^2\text{d}r.\end{array}$ (72)

Similarly, there also holds

$\begin{array}{l}{\displaystyle {\int }_{r_1}^{r_2}}\frac{\text{d}}{\text{d}r}\left(\frac{2m\Omega }{r^{2}C\left(r\right)s}\right){\left|r{\xi }_r\right|}^2\text{d}r\\ ={\left(\frac{2m\Omega }{r^{2}C\left(r\right)s}{\left|r{\xi }_r\right|}^2\right)|}_{r=r_1\mathrm{,}{r}_2}-{\displaystyle {\int }_{r_1}^{r_2}}\frac{4m\Omega }{r^{2}C\left(r\right)s}\left|r{\xi }_r\right|\frac{\text{d}}{\text{d}r}\left|r{\xi }_r\right|\text{d}r\\ =-{\displaystyle {\int }_{r_1}^{r_2}}\frac{4m\Omega }{r^{2}C\left(r\right)s}\left|r{\xi }_r\right|\frac{\text{d}}{\text{d}r}\left|r{\xi }_r\right|\text{d}r\mathrm{.}\end{array}$ (73)

Combining (71)-(73), we get

${\displaystyle {\int }_{r_1}^{r_2}}\frac{1}{rC\left(r\right)}{\left(\frac{\text{d}}{\text{d}r}\left|r{\xi }_r\right|-\frac{2m\Omega }{rs}\left|r{\xi }_r\right|\right)}^2\text{d}r={\displaystyle {\int }_{r_1}^{r_2}}\frac{\Phi \left(r\right)-s^2}{r}{\left|r{\xi }_r\right|}^2\text{d}r.$ (74)

Noting that $s=\left(a+m\Omega \right)+bi$ and $s^2=\left[{\left(a+m\Omega \right)}^2-b^2\right]+2\left(a+m\Omega \right)bi$, we get

$\begin{array}{c}C\left(r\right)=\frac{m^2}{r^2{s}^2}+\frac{k^2}{s^2-\frac{\text{d}\bar{\rho }}{\text{d}z}}\\ =\frac{m^2}{r^2}\frac{\left[{\left(a+m\Omega \right)}^2-b^2\right]-2\left(a+m\Omega \right)bi}{{\left[{\left(a+m\Omega \right)}^2-b^2\right]}^2+4{\left(a+m\Omega \right)}^2{b}^2}\\ +k^2\frac{\left[{\left(a+m\Omega \right)}^2-b^2-\frac{\text{d}\bar{\rho }}{\text{d}z}\right]-2\left(a+m\Omega \right)bi}{{\left[{\left(a+m\Omega \right)}^2-b^2-\frac{\text{d}\bar{\rho }}{\text{d}z}\right]}^2+4{\left(a+m\Omega \right)}^2{b}^2}\\ :=A-Bi\end{array}$ (75)

where

$(\begin{array}{l}A=\frac{m^2\left[{\left(a+m\Omega \right)}^2-b^2\right]}{r^2{\left[{\left(a+m\Omega \right)}^2-b^2\right]}^2+4{\left(a+m\Omega \right)}^2{b}^2}+\frac{k^2\left[{\left(a+m\Omega \right)}^2-b^2-\frac{\text{d}\bar{\rho }}{\text{d}z}\right]}{{\left[{\left(a+m\Omega \right)}^2-b^2-\frac{\text{d}\bar{\rho }}{\text{d}z}\right]}^2+4{\left(a+m\Omega \right)}^2{b}^2},\\ B=\frac{2m^2\left(a+m\Omega \right)b}{r^2{\left[{\left(a+m\Omega \right)}^2-b^2\right]}^2+4{\left(a+m\Omega \right)}^2{b}^2}+\frac{2k^2\left(a+m\Omega \right)b}{{\left[{\left(a+m\Omega \right)}^2-b^2-\frac{\text{d}\bar{\rho }}{\text{d}z}\right]}^2+4{\left(a+m\Omega \right)}^2{b}^2}.\end{array}$ (76)

Substitute Equation (75) into Equation (74), we obtain

$\begin{array}{l}{\displaystyle {\int }_{r_1}^{r_2}}\frac{A+Bi}{r\left(A^2+B^2\right)}{\left[\frac{\text{d}\left(r{\xi }_r\right)}{\text{d}r}-\frac{2m\Omega }{r\left(a+m\Omega +bi\right)}\left(r{\xi }_r\right)\right]}^2\text{d}r\\ ={\displaystyle {\int }_{r_1}^{r_2}}\frac{\Phi \left(r\right)-{\left(a+m\Omega \right)}^2+b^2-2\left(a+m\Omega \right)bi}{r}{\left|r{\xi }_r\right|}^2\text{d}r.\end{array}$ (77)

The real and the imaginary parts must be considered separately. From (77) we obtain

$\begin{array}{l}{\displaystyle {\int }_{r_1}^{r_2}}\frac{A}{r\left(A^2+B^2\right)}{\left|\frac{\text{d}\left(r{\xi }_r\right)}{\text{d}r}\right|}^2\text{d}r-{\displaystyle {\int }_{r_1}^{r_2}}\frac{2m\Omega \left[A\left(a+m\Omega \right)+Bb\right]}{r^2\left(A^2+B^2\right)\left[{\left(a+m\Omega \right)}^2+b^2\right]}\frac{\text{d}}{\text{d}r}{\left(r{\xi }_r\right)}^2\\ +{\displaystyle {\int }_{r_1}^{r_2}}\left(4m^2{\Omega }^2\right)\frac{\left(A{\left[a+m\Omega \right]}^2-b^2\right)+2B\left(a+m\Omega \right)b}{r^3\left(A^2+B^2\right){\left[{\left(a+m\Omega \right)}^2-b^2\right]}^2+4{\left(a+m\Omega \right)}^2{b}^2}{\left|r{\xi }_r\right|}^2\text{d}r\\ ={\displaystyle {\int }_{r_1}^{r_2}}\frac{\Phi \left(r\right)-{\left(a+m\Omega \right)}^2+b^2}{r}{\left|r{\xi }_r\right|}^2\text{d}r.\end{array}$ (78)

and

$\begin{array}{l}{\displaystyle {\int }_{r_1}^{r_2}}\frac{B}{r\left(A^2+B^2\right)}{\left|\frac{\text{d}\left(r{\xi }_r\right)}{\text{d}r}\right|}^2\text{d}r+{\displaystyle {\int }_{r_1}^{r_2}}\frac{2m\Omega \left[Ab-\left(a+m\Omega \right)B\right]}{r^2\left(A^2+B^2\right)\left[{\left(a+m\Omega \right)}^2+b^2\right]}\frac{\text{d}}{\text{d}r}{\left(r{\xi }_r\right)}^2\\ -{\displaystyle {\int }_{r_1}^{r_2}}\left(4m^2{\Omega }^2\right)\frac{2A\left(a+m\Omega \right)b-B\left[{\left(a+m\Omega \right)}^2-b^2\right]}{r^3\left(A^2+B^2\right){\left[{\left(a+m\Omega \right)}^2-b^2\right]}^2+4{\left(a+m\Omega \right)}^2{b}^2}{\left|r{\xi }_r\right|}^2\text{d}r\\ ={\displaystyle {\int }_{r_1}^{r_2}}\frac{2\left(a+m\Omega \right)b}{r}{\left|r{\xi }_r\right|}^2\text{d}r.\end{array}$ (79)

It is obvious that $b=0$ (that is $s\in R$ ) is the necessary and sufficient condition for linearized stable to system (33)-(35).

At the beginning, we assume that $s\in ℝ$. Testing the Equation (33) by $r{\xi }_r^{\mathrm{<mo>\; *\; </mo>}}$, we obtain

$\begin{array}{c}{\displaystyle {\int }_{r_1}^{r_2}}r\left[s^2-\Phi \left(r\right)\right]{\left|{\xi }_r\right|}^2={\displaystyle {\int }_{r_1}^{r_2}}r{\xi }_r^{*}\frac{\text{d}P}{\text{d}r}+{\displaystyle {\int }_{r_1}^{r_2}}r{\xi }_r^{*}\frac{2m\Omega }{sr}P\text{d}r\\ ={\left(r{\xi }_r^{*}P\right)|}_{r=r_1,r_2}-{\displaystyle {\int }_{r_1}^{r_2}}P\text{d}\left(r{\xi }_r^{*}\right)+{\displaystyle {\int }_{r_1}^{r_2}}P\left({\xi }_r^{*}\frac{2m\Omega }{s}P\right)\text{d}r\\ =-{\displaystyle {\int }_{r_1}^{r_2}}P\left[\frac{\text{d}}{\text{d}r}\left(r{\xi }_r^{*}\right)-\frac{2m\Omega }{s}{\xi }_r^{*}\right]\text{d}r.\end{array}$ (80)

From (35), by taking the complex conjugate, we have

$\frac{\text{d}}{\text{d}r}\left(r{\xi }_r^{*}\right)-\frac{2m\Omega }{s}{\xi }_r^{*}=r\left(\frac{m^2}{r^2{s}^2}+\frac{k^2}{s^2-\frac{\text{d}\bar{\rho }}{\text{d}z}}\right)P^{*}.$ (81)

Substitute (81) into (80), we obtain

${\displaystyle {\int }_{r_1}^{r_2}}r\left(s^2-\Phi \left(r\right)\right){\left|{\xi }_r\right|}^2\text{d}r=-{\displaystyle {\int }_{r_1}^{r_2}}\left(\frac{m^2}{r^2{s}^2}+\frac{k^2}{s^2-\frac{\text{d}\bar{\rho }}{\text{d}z}}\right)r{\left|P\right|}^2\text{d}r,$ (82)

from which we get

$\begin{array}{c}{k}^2=\frac{{\displaystyle {\int }_{r_1}^{r_2}}r\left[\left(\Phi -s^2\right){\xi }_r^2-\left(m^2{P}^2\right)/\left(r^2{s}^2\right)\right]\text{d}r}{{\displaystyle {\int }_{r_1}^{r_2}}r{P}^2/\left(s^2-\frac{\text{d}\bar{\rho }}{\text{d}z}\right)\text{d}r}\\ =-\frac{{\displaystyle {\int }_{r_1}^{r_2}}\left[\frac{1}{s^2-\Phi }{\left(\frac{\text{d}P}{\text{d}r}+\frac{2m\Omega }{sr}P\right)}^2+\frac{m^2{P}^2}{r^2{s}^2}\right]r\text{d}r}{{\displaystyle {\int }_{r_1}^{r_2}}r{P}^2/\left(s^2-\frac{\text{d}\bar{\rho }}{\text{d}z}\right)\text{d}r}.\end{array}$ (83)

In the last step above, we use the Equation (33).

Noting that $k^2\ge 0$, and write

$I_1=-{\displaystyle {\int }_{r_1}^{r_2}}\left[\frac{1}{s^2-\Phi }{\left(\frac{\text{d}P}{\text{d}r}+\frac{2m\Omega }{sr}P\right)}^2+\frac{m^2{P}^2}{r^2{s}^2}\right]r\text{d}r,$ (84)

and

$I_2={\displaystyle {\int }_{r_1}^{r_2}}\frac{r{P}^2}{s^2-\frac{\text{d}\bar{\rho }}{\text{d}z}}\text{d}r\mathrm{,}$ (85)

we have the following results:

1) $\frac{\text{d}\bar{\rho }}{\text{d}z}<s^2$. From the Equation (85), we obtain $I_2>0$. When $\Phi \left(r\right)<s^2$ then $I_1<0$, combining with $I_2>0$, we obtain $k^2<0$. Therefore the assumption $s\in ℝ$ is not true. In other words, s must be complex, which means that the system is unstable.

2) $\frac{\text{d}\bar{\rho }}{\text{d}z}>{\left(-\omega +m\Omega \right)}^2=s^2$. From the Equation (85), we obtain $I_2<0$. When $\Phi \left(r\right)<s^2$ then $I_1<0$, combining with $I_2<0$, we obtain $k^2>0$. Therefore $k^2>0$ is equivalent to $s\in ℝ$. In this situation, the system is linear stable.

Moreover, we also give a necessary condition for the linear stable. To do this, let $b=0$, we get

$C\left(r\right)=\frac{m^2}{r^2{\left(a+m\Omega \right)}^2}+\frac{k^2}{{\left(a+m\Omega \right)}^2-\frac{\text{d}\bar{\rho }}{\text{d}z}}\mathrm{.}$ (86)

From the Equations (78)-(79) and (74), we get the following necessary condition for the linear stable:

3) When $\Phi \left(r\right)>s^2\ge 0$, namely $\Phi \left(r\right)>{\left(a+m\Omega \right)}^2\ge 0$. From the Equation (74), we know that $\frac{\text{d}\bar{\rho }}{\text{d}z}$ must satisfy $\frac{\text{d}\bar{\rho }}{\text{d}z}<{\left(a+m\Omega \right)}^2=s^2$ or $\frac{\text{d}\bar{\rho }}{\text{d}z}>\frac{k^2{r}^2+m^2}{m^2}{\left(a+m\Omega \right)}^2=\frac{k^2{r}^2+m^2}{m^2}{s}^2$ ;

4) When $\Phi \left(r\right)<s^2$, namely $\Phi \left(r\right)<{\left(a+m\Omega \right)}^2$. From the Equation (74), $\frac{\text{d}\bar{\rho }}{\text{d}z}$ must satisfy ${\left(a+m\Omega \right)}^2<\frac{\text{d}\bar{\rho }}{\text{d}z}<\frac{k^2{r}^2+m^2}{m^2}{\left(a+m\Omega \right)}^2$.

5. Further Discussion on the Axisymmetric Perturbations Case

To study the role of $\Phi \left(r\right)$ and $\frac{\text{d}\bar{\rho }\left(z\right)}{\text{d}z}$ in the stabilize or destabilize the flow in a more obvious way, in this section, we shall present some simulation by taking $\Phi \left(r\right)$ and $\frac{\text{d}\bar{\rho }\left(z\right)}{\text{d}z}$ with some special value. Moreover, for simplicity, we constrain for the case with axisymmetric perturbations on the domain $\pi \le r\le 2\pi$ and $0\le z\le \pi$. In this case, we get the system as

$(\begin{array}{l}{\partial }_t\rho +u_z{\partial }_z\bar{\rho }=0,(87\text{-}1)\hfill \\ {\partial }_t{u}_r-2\frac{V}{r}{u}_{\theta }+{\partial }_{r}P=0,(87\text{-2})\hfill \\ {\partial }_t{u}_{\theta }+\left(\frac{\text{d}V}{\text{d}r}+\frac{V}{r}\right)u_r=0,\pi \le r\le 2\pi ,0\le z\le \pi ,(87\text{-3})\hfill \\ {\partial }_t{u}_z+{\partial }_{z}P=\rho ,(87\text{-4})\hfill \\ {\partial }_r{u}_r+\frac{u_r}{r}+{\partial }_z{u}_z=0.(87\text{-5})\hfill \end{array}$

We first rewrite the linearized system (87) in a scalar equation.

Applying ${\partial }_z$ and ${\partial }_r$ to the Equations (87-2) and (87-4) respectively, eliminating P, we get

${\partial }_t\left({\partial }_r{u}_z-{\partial }_z{u}_r\right)+2\Omega {\partial }_z{u}_{\theta }={\partial }_r\rho \mathrm{.}$ (88)

Similarly, from (87-1) and (87-3), we obtain

${\partial }_t\left({\partial }_r\rho -2\Omega {\partial }_z{u}_{\theta }\right)+{\partial }_z\bar{\rho }{\partial }_r{u}_z-\Phi \left(r\right){\partial }_z{u}_r=0,$ (89)

where $\Phi \left(r\right)=2\Omega \left(\frac{\text{d}V\left(r\right)}{\text{d}r}+\frac{V\left(r\right)}{r}\right)$ is the Rayleigh discriminant. Combine the Equations (88) and (89), we have

${\partial }_{tt}\left({\partial }_r{u}_z-{\partial }_z{u}_r\right)=\Phi \left(r\right){\partial }_z{u}_r-{\partial }_z\bar{\rho }{\partial }_r{u}_z\mathrm{.}$ (90)

By using the incompressible conditions (87-5), as to the axisymmetric case, we write

$r{u}_r={\partial }_z\psi \text{and}r{u}_z=-{\partial }_r\psi ,\pi \le r\le 2\pi ,0\le z\le \pi ,$ (91)

with $\psi \left(t\mathrm{,}r\mathrm{,}z\right)$ be a scalar function and satisfying

${{\partial }_z\psi |}_{r=\pi ,2\pi }=0,{{\partial }_r\psi |}_{z=0,\pi }=0.$ (92)

Then we analyze the Equations (90)-(91), we get

$(\begin{array}{l}{\partial }_{tt}\left({\partial }_{rr}\psi +{\partial }_{zz}\psi -\frac{1}{r}{\partial }_r\psi \right)+\Phi \left(r\right){\partial }_{zz}\psi +{\partial }_z\bar{\rho }\left({\partial }_{rr}\psi -\frac{1}{r}{\partial }_r\psi \right)=0,\pi \le r\le 2\pi ,0\le z\le \pi ,(93\text{-}1)\\ {{\partial }_z\psi |}_{r=\pi ,2\pi }=0,{{\partial }_r\psi |}_{z=0,\pi }=0.(93\text{-}2)\end{array}$

Let $\psi \left(t\mathrm{,}r\mathrm{,}z\right)={\psi }_1\left(t\right){\gamma }_2\left(r\right){\gamma }_3\left(z\right)$, next by using the boundary condition (93-2), we take ${\gamma }_2\left(r\right)={\displaystyle {\sum }_{k_0}}{a}_m\sin \left(k_{0}r\right)$ and ${\gamma }_3\left(z\right)={\displaystyle {\sum }_k}{b}_k\sin \left(kz\right)$. From (93), we take

$\begin{array}{l}{C}_0\left(k_0,k,r\right){\partial }_{tt}{\psi }_1\left(t\right)\\ +\left[\left(k^2\Phi \left(r\right)+{\partial }_z\bar{\rho }{k}_0^2\right)\sin \left(k_{0}r\right)+{\partial }_z\bar{\rho }\frac{k_0}{r}\cos \left(k_{0}r\right)\right]{\psi }_1\left(t\right)=0,\end{array}$ (94)

where

$C_0\left(k_0,k,r\right)=\left(k_0^2+k^2\right)\sin \left(k_{0}r\right)+\frac{k_0}{r}\cos \left(k_{0}r\right).$

Case 1. $C_0\left(k_0\mathrm{,}k\mathrm{,}r\right)=0$. From Equation (94), to find a non-trivial solution, we have

$\left(k^2\Phi \left(r\right)+{\partial }_z\bar{\rho }{k}_0^2\right)\sin \left(k_{0}r\right)+{\partial }_z\bar{\rho }\left(z\right)\frac{k_0}{r}\cos \left(k_{0}r\right)=0.$ (95)

From which, we get $\Phi \left(r\right)={\partial }_z\bar{\rho }=\lambda$, with $\lambda$ be a uniform constant.

Write $\Gamma \left(t,r,z\right)={\partial }_{rr}\psi +{\partial }_{zz}\psi -\frac{1}{r}{\partial }_r\psi$, combing (93), we conclude that:

· When $\Phi \left(r\right)={\bar{\rho }}^{\prime }\left(z\right)=\lambda >0$, the general solution can be written as

$\Gamma =\Gamma \left(0,r,z\right)\cdot \left[A_1{\text{e}}^{\sqrt{\lambda }t}+A_2{\text{e}}^{-\sqrt{\lambda }t}\right],$ (96)

with $A_0$ and $A_1$ be constants. In this case, it is unstable for any non-trivial perturbation.

· When $\Phi \left(r\right)={\bar{\rho }}^{\prime }\left(z\right)=\lambda \le 0$, the solution is

$\Gamma \left(t,r,z\right)=\Gamma \left(0,r,z\right)\cdot \left[B_1\cos \left(\sqrt{-\lambda }t\right)+B_2\sin \left(\sqrt{-\lambda }t\right)\right],$ (97)

with $B_0$ and $B_1$ be constants. In this case, the linearized system is stable.

Case 2. $C_0\ne 0$. From the Equation (94), we get

${\partial }_{tt}{\psi }_1\left(t\right)+\frac{\left(k^2\Phi \left(r\right)+{\partial }_z\bar{\rho }{k}_0^2\right)\sin \left(k_{0}r\right)+{\partial }_z\bar{\rho }\frac{k_0}{r}\cos \left(k_{0}r\right)}{C_0\left(k_0,k,r\right)}{\psi }_1\left(t\right)=0.$ (98)

Write

$\Lambda =-\frac{\left(k^2\Phi \left(r\right)+{\partial }_z\bar{\rho }{k}_0^2\right)\sin \left(k_{0}r\right)+{\partial }_z\bar{\rho }\frac{k_0}{r}\cos \left(k_{0}r\right)}{C_0\left(k_0,k,r\right)}.$ (99)

To simplify this equation, we get

$\Lambda =-\frac{\left(k^2\Phi \left(r\right)+{\bar{\rho }}^{\prime }\left(z\right)k_0^2\right)r+k_0{\bar{\rho }}^{\prime }\left(z\right)\text{ctan}\left(k_{0}r\right)}{\left(k_0^2+k^2\right)r+m\cdot \text{ctan}\left(k_{0}r\right)}.$ (100)

It is obviously, the flow will be stable with $\Lambda \le 0$ and is unstable when $\Lambda >0$. In the following, we shall present some numerical simulation by taking $\Phi \left(r\right)$, ${\bar{\rho }}^{\prime }\left(z\right)$, and the wave number $k_0$, k with some special value.

· First, in Figure 2 and Figure 3, we take $\Phi \left(r\right)=-1$, ${\bar{\rho }}^{\prime }\left(z\right)=-2$ and $\Phi \left(r\right)=-2$, ${\bar{\rho }}^{\prime }\left(z\right)=-1$ respectively, In the graph, we also see that the ordinate values are all greater than zero on the interval. which are unstable cases.

· The next, in Figure 4 and Figure 5, we take $\Phi \left(r\right)=1$, ${\bar{\rho }}^{\prime }\left(z\right)=-1$ and $\Phi \left(r\right)=-1$, ${\bar{\rho }}^{\prime }\left(z\right)=1$ respectively, In that graph, we see that the ordinate values are less than zero on the interval. which are unstable cases also. In Figure 4, the instability comes from the thermal convection, more precisely, the higher frequency m with the variable z is more sensible. In Figure 5, the instability comes from the rotation. In this case, the higher frequency k with the variable r is more sensible, which implies the thermal buoyancy smoothed the lower frequency k with the variable r.

Remark 5.1 The effect from the boundary. During the numerical simulation above, whatever the sign for both $\Phi$ and ${\bar{\rho }}^{\prime }\left(z\right)$, there is rapid oscillation near $r=2\pi$. Since $\text{ctan}\left(r\right)\to \infty$ as $r\to 2\pi$, we guess this oscillation is from the function $\text{ctan}\left(r\right)$ at the beginning. However, from the results of numerical simulation, the oscillation is more likely an effect from the boundary. To clarify the situation, we shall take the following simulations:

· In Figure 6, we do the simulation with $r\in \left[\pi \mathrm{,}\frac{19}{10}\pi \right]$, oscillation vanished. Subsequently, In Figure 7, we do the simulation with $r\in \left[\pi \mathrm{,}\frac{199}{100}\pi \right]$, oscillation appeared. Apparently, there are some extra influence from the boundary in the interval $r\in \left(\frac{19}{10}\pi \mathrm{,}\frac{199}{100}\pi \right)$.

Figure 2. $\Phi =-1,{\bar{\rho }}^{\prime }=-2$.

Figure 3. $\Phi =-2,{\bar{\rho }}^{\prime }=-1$.

Figure 4. $\Phi =1,{\bar{\rho }}^{\prime }=-1$.

Figure 5. $\Phi =-1,{\bar{\rho }}^{\prime }=1$.

Figure 6. $\Phi =2,{\bar{\rho }}^{\prime }=1,\pi \le r\le \frac{19}{10}\pi$.

Figure 7. $\Phi =1,{\bar{\rho }}^{\prime }=2,\pi \le r\le \frac{199}{100}\pi$.

6. Conclusion

In this paper, we analyze the stability and instability criteria for the coupled thermal effects of fluids between coaxial rotating cylinders. The perturbation equation is analyzed by normal-modes method. We extend Rayleigh stability criterion by the analysis of axisymmetric perturbation in some cases, and we also analyze the case of non-axisymmetric perturbation, the results are presented in Theorem 1 and Theorem 2. Finally, through numerical simulation experiments, the results obtained by our experiments are consistent with the results obtained by our analysis on the specific cases of axisymmetric perturbations under certain given special values. For the fluid instability near the cylinder boundary, we also found some new problems waiting for us to further deal with.

Acknowledgements

This work is partially supported by the NSFC (No. 11971199) and Guangdong Basic and Applied Basic Research Foundation (No. 2020B1515310012).

NOTES

*First author.

^#Corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Rayleigh, F.R.S. (1887) On the Stability or Instability of Certain Fluid Motions II. Proceedings of the London Mathematical Society, s1-19, 67-75. https://doi.org/10.1112/plms/s1-19.1.67
[2]	Chandrasekhar, S. (1981) Hydrodynamic and Hydromagnetic Stability. Dover Publications, Mineola.
[3]	Drazin, P.G. and Reid, W.H. (2004) Hydrodynamic Stability. Vol. 2, Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9780511616938
[4]	Taylor, G.I. (1923) Stability of a Viscous Liquid Contained between Two Rotating Cylinders. Philosophical Transactions of the Royal Society of London, 223, 289-343. https://doi.org/10.1098/rsta.1923.0008
[5]	Rayleigh, L. (1917) On the Dynamics of Revolving Fluids. Proceedings of the Royal Society of London, 93, 148-154. https://doi.org/10.1098/rspa.1917.0010
[6]	Di Prima, R.C. and Swinney, H.L. (1981) Instabilities and Transition in Flow between Concentric Rotating Cylinders. Springer, Berlin. https://doi.org/10.1007/3-540-13319-4_16
[7]	Withjack, E.M. (1974) An Experimental Study of Couette Instability of Stratified Fluids. Journal of Fluid Mechanics, 66, 725-737. https://doi.org/10.1017/S0022112074000474
[8]	Andereck, C.D., Liu, S.S. and Swinney, H.L. (1986) Flow Regimes in a Circular Couette System with Independently Rotating Cylinders. Journal of Fluid Mechanics, 164, 155-183. https://doi.org/10.1017/S0022112086002513
[9]	Dubrulle, B., Dauchot, O., Daviaud, F., Longaretti, P.Y. and Zahn, J.P. (2005) Stability and Turbulent Transport in Taylor-Couette Flow from Analysis of Experimental Data. Physics of Fluids, 17, 531-578. https://doi.org/10.1063/1.2008999
[10]	Hua, B.L., Gentil, S.L. and Orlandi, P. (1997) First Transitions in Circular Couette Flow with Axial Stratification. Physics of Fluids, 9, 365-375. https://doi.org/10.1063/1.869237
[11]	Caton, F., Janiaud, B. and Hopfinger, E.J. (2000) Stability and Bifurcations in Stratified Taylor-Couette Flow. Journal of Fluid Mechanics, 419, 93-124. https://doi.org/10.1017/S0022112000001348
[12]	Boubnov, B.M., Gledzer, E.B. and Hopfinger, E.J. (1995) Stratified Circular Couette Flow: Instability and Flow Regimes. Journal of Fluid Mechanics, 292, 333-358. https://doi.org/10.1017/S0022112095001558
[13]	Hua, B.L., Moore, D.W. and Gentil, S.L. (1997) Inertial Nonlinear Equilibration of Equatorial Flows. Journal of Fluid Mechanics, 331, 345-371. https://doi.org/10.1017/S0022112096004016
[14]	Dubrulle, B., Mari, L., Normand, C., Richard, D., Hersant, F. and Zahn, J.-P. (2005) A Hydrodynamic Shear Instability in Stratified Disks. Astronomy and Astrophysics, 429, 1-13. https://doi.org/10.1051/0004-6361:200400065
[15]	Reid, W.H. (1960) Inviscid Modes of Instability in Couette Flow. Journal de Mathématiques Pures et Appliquées, 1, 411-422. https://doi.org/10.1016/0022-247X(60)90013-5
[16]	Park, J. and Billant, P. (2013) The Stably Stratified Taylor-Couette Flow Is Always Unstable Except for Solid-Body Rotation. Journal of Fluid Mechanics, 725, 262-280. https://doi.org/10.1017/jfm.2013.186
[17]	Yavneh, I., Mcwilliams, J.C. and Molemaker, M.J. (2001) Non-Axisymmetric Instability of Centrifugally Stable Stratified Taylor-Couette Flow. Communications on Pure & Applied Analysis, 448, 1-21. https://doi.org/10.1017/S0022112001005109
[18]	Hoffmann, C., Altmeyer, S., Heise, M., Abshagen, J. and Pfister, G. (2013) Axisymmetric Propagating Vortices in Centrifugally Stable Taylor-Couette Flow. Journal of Fluid Mechanics, 728, 458-470. https://doi.org/10.1017/jfm.2013.283
[19]	Pedlosky, J. (1987) Geophysical Fluid Dynamics. Second Edition, Springer, Berlin. https://doi.org/10.1007/978-1-4612-4650-3
[20]	Majda, A. (2003) Introduction to PDEs and Waves for the Atmosphere and Ocean. American Mathematical Society, Providence. https://doi.org/10.1090/cln/009
[21]	Pedlowsky, J. and Miles, J. (2004) Waves in the Ocean and Atmosphere: Introduction to Wave Dynamics. Applied Mechanics Reviews, 57, B20. https://doi.org/10.1115/1.1786585