Natural convection between infinite concentric cylinders

M. Prud'homme and P. Le Quéré

Introduction

Many configurations of industrial interest are cylindrical containers isolated from their surroundings by an annular space filled with a low conductivity fluid such as air. The fluid in this annular space is submitted to a horizontal temperature difference resulting in a global motion whose net vertical mass flux is zero if the annular space is sealed at both ends. In many cases the container is set in an environment which is vertically thermally stratified. If the annular space is thin compared to its height one can assume that the flow takes place in an infinite annular slot whose wall temperatures increase linearly with height in such a way that at any elevation the temperature difference across the walls is constant. In the next section we will determine the base flow structure and the corresponding Nusselt number. One also wants to determine wether the base flow is stable, to determine the nature of the bifurcation, and to investigate the properties of the fully non-linear flows resulting from the instabilities.

Base flow

Consider the space between two vertical concentric cylinders of inner and outer radius $R_1$ and $R_2$. We assume that the wall temperature of the innner cylinder is $T_0 + \frac{\Delta T }{2}+ \alpha z$ and than that of the outer cylinder reads $T_0 - \frac{\Delta T }{2} + \alpha z$, where $T_0$ is a reference temperature, $\Delta T$ is the constant temperature difference across the slot and $\alpha$ the constant vertical temperature gradient. The slot is filled with a fluid of thermal diffusivity $\kappa$ and of kinematic diffusivity $\nu$. We assume that the fluid follows the Boussinesq approximation, ie $\rho = \rho_0 ( 1 - \beta (T- T_0))\,$ where $\beta = - \frac{1}{\rho_0} \frac{\partial \rho}{\partial T}$. Looking for a parallel solution such as $(U,W)= (0, w(r))$ the governing equations reduce to : $\begin{matrix} \frac{ d P}{ d z} - g \beta \alpha z & = &\nu D w + g \beta \Theta \\ w \alpha &= & \kappa D (\theta) \frac{ }{ } \end{matrix}$

where $D = \frac{1}{r} \frac{\partial}{\partial r} r \frac{\partial }{\partial r}$ and $\Theta(r) = T(r,z) - T_0 - \alpha z \,$. The left and right members of the first equation are functions of $z$ and $r$ respectively and thus equal to a constant which can be written $g \beta \Theta_r$, resulting in the system : $\begin{matrix} \nu D(w) + g \beta ( \Theta - \Theta_r) & = &0 \\ w \alpha &= & \kappa D (\Theta) \end{matrix}$

Upon introducing $\Delta R = R_2-R_1$ as reference length and $\Delta T$ as reference temperature difference, the elimination of $\Theta$ between the two equations shows that $w$ satisfies the fourth order differential equation : $D^2( w) + Ra S w = 0 \,$ where $Ra = \frac{g \beta \Delta T \Delta R^3 } { \nu \kappa}$, and $S =\frac{\alpha \Delta R}{\Delta T} \,$. Thus $w$ expresses as $w = C_1 J_0(\lambda \gamma r) + C_2 Y_0(\lambda \gamma r) + C_3 I_0(\lambda \gamma r) + C_4 K_0(\lambda \gamma r)\,$ where $J_0, Y_0$ are Bessel functions and $I_0, K_0$ modified Bessel functions of the first and second kind respectively with $\lambda =1+ i$, and $\gamma$ such as $4 \gamma^4 = Ra S \,$. As $w = \frac{1}{r}\frac{\partial \Psi}{\partial r}$, it follows that $\frac{d \Psi}{d r} = r (C_1 J_0(\lambda \gamma r) + C_2 Y_0(\lambda \gamma r) + C_3 I_0(\lambda \gamma r) + C_4 K_0(\lambda \gamma r)).$ Integrating this equation, making use of algebraic properties of Bessel functions \cite{abram} and imposing that the two walls are no-slip and iso-flowrates boundaries shows that the 4 coefficients satisfy the following linear system $\begin{bmatrix} J_0(\lambda \gamma r_1) & Y_0(\lambda \gamma r_1) & I_0(\lambda \gamma r_1)& K_0(\lambda \gamma r_1)\\ J_0(\lambda \gamma r_2) & Y_0(\lambda \gamma r_2) & I_0(\lambda \gamma r_2)& K_0(\lambda \gamma r_2) \\ J_1(\lambda \gamma r_1) & Y_1(\lambda \gamma r_1) & I_1(\lambda \gamma r_1)& -K_1(\lambda \gamma r_1) \\ J_1(\lambda \gamma r_2) & Y_1(\lambda \gamma r_2) & I_1(\lambda \gamma r_2)& -K_1(\lambda \gamma r_2) \\ \end{bmatrix} \begin{bmatrix} C_1\\ C_2 \\ C_3 \\ C_4 \\ \end{bmatrix} = \begin{bmatrix} 0\\ 0\\ \lambda \gamma / r_1 \\ \lambda \gamma / r_2\\ \end{bmatrix}$

Due to the asymptotic behavior of the Bessel functions, this system becomes ill-conditioned when the product $\gamma r$ becomes large (values of 10 or so), that is when the solution is composed of two distinct boundary layers.

This situation provides an alternative way to determine the base flow, through a matched asymptotic expansion of the inner and outer solutions. The procedure is described in Prud'homme and Le Quéré. As expected the expansion converges quickly when the two boundary layers are indeed distinct, ie when the stream function achieves a constant value in the central core region. The figure above gives a comparison of the base flow solution obtained when both methods work equally well. Both methods are thus complementary and provide access to the base flow solution for the entire range of governing parameters.

Linear instability

The stability of this base flow is first investigated in the framework of the linear stability theory, ie under the assumption of disturbances of infinitesimal amplitude expanded in normal modes. Perturbations are written as $\psi (r) \exp (i k x + \sigma t )$ where $k$ is the real wave umber and $\sigma = \sigma_r + i \sigma_i$ is the complex growth factor. $\sigma_i \le 0$ corresponds to instability and $- \sigma_r/k$ is the wave speed. The generalized eigenvalue problem is solved for $k$ fixed with the help of the library. A secant method is used to determine the value of $Ra$ for which $\sigma_i = 0$ . The linear stability was solved for a wide range of $\gamma$ and Prandtl number values.

A sample result is presented in figure 2, corresponding to a flow of air in the case of strong curvature effects and a strong stratification. The figure shows that the instability mode is concentrated along the inner heated cylinder, and its phase velocity indicates that it is an upward moving wave.

Weakly nonlinear instability

Additional computations were performed in order to investigate the nature of the bifurcation. These calculations involve projection of the nonlinear on the adjoint of the most unstable eigenmode in order to compute the Landau coefficients. Details are given in Prud'homme and Le Quéré . The results show that the Hopf bifurcation can be either subcritical or supercritical depending on the parameters.