G Petrone
email: petrone@univ-mlv.fr

E Chénier
email: chenier@univ-mlv.fr

G Lauriat
email: lauriat@univ-mlv.fr

Stability of free convection in air-filled horizontal annuli: Influence of the radius ratio

Linear stability of two-dimensional natural convection in air-filled horizontal annuli is numerically investigated for radius ratios in the range 1.2 ≤ R ≤ 3 and for Rayleigh numbers less than 10 4 . Bifurcation diagrams are obtained for various radius ratios and the main thresholds are tracked as a function of R. A new instability mode has been highlighted which breaks the symmetry of the basic flow. This result demonstrates the need of modeling the annular gap without assuming flow symmetry. In addition to bifurcation maps drawn in the Rayleigh number-radius ratio plane, a map of possible flow-patterns is also established. This map allows to foresee the number of solutions and the corresponding flow structures.

Natural convection between horizontal isothermal concentric cylinders has been extensively studied because of its importance in many technological applications such as heat exchangers. From a theoretical point of view, this problem has been considered in a number of studies because a large variety of flow structures is encountered in this configuration according to the aspect ratio, Rayleigh and Prandtl numbers.

For low Ra, the basic flow field consists of two two-dimensional crescentshaped cells, symmetrical with respect to the vertical plane containing the axes of the cylinders. In each of the two annular half spaces, the fluid goes upwards and downwards, respectively along the hot inner and cold outer cylinders. The conduction is the major mode of heat transfer between the differentially heated boundaries of the annulus. As the Rayleigh number is increased, the center of rotation of the main cells moves upward and a thermal plume starts to form at the upper part of the annulus with an impingement region on the outer cylinder. The distribution of the thermal fluxes along the inner and outer cylinders shows that the largest part of the heat convected within the annulus is extracted from the lower part of the inner cylinder.

From their experimental work, Powe et al. [START_REF] Powe | Free convective flow patterns in cylindrical annuli[END_REF][START_REF] Bishop | Natural convective oscillatory flow in cylindrical annuli[END_REF] depicted flow regimes and spatial patterns for air-filled annuli as a function of the radius ratio and Rayleigh number. For wide gap annuli, i.e. for radius ratio larger than R = 1.71, transitions happen from two-dimensional steady to oscillatory flows. For moderate gap annuli (1.24 < R < 1.71) a three-dimensional spiral flow is observed past the transition, while for narrow gap annuli (R < 1.24) the basic one-cellular flow changes to multi-cellular flows. Later, experimental and numerical studies conducted by Rao et al. [START_REF] Rao | Flow patterns of natural convection in horizontal cylindrical annuli[END_REF] yielded results in qualitative agreement with the Powe experimental map for moderate and narrow gap annuli. However, they failed to detect the oscillatory two-dimensional flow régime for wide gap annuli. Vafai and Ettefagh [START_REF] Vafai | An investigation of transient three-dimensional buoyancy-driven flow and heat transfer in a closed horizontal annulus[END_REF] and more recently Cadiou [START_REF] Cadiou | Contribution à l'étude numérique des transferts de chaleur par convection naturelle ou par convection mixte dans une cavité annulaire horizontale[END_REF] and Dyko et al. [START_REF] Dyko | A numerical and experimental investigation of stability of natural convective flows within a horizontal annulus[END_REF][START_REF] Dyko | Three-dimensional natural convective states in a narrow-gap horizontal annulus[END_REF][START_REF] Dyko | On the presence of odd transverse convective rolls in narrow-gap horizontal annuli[END_REF] numerically studied the effect of the axial confinement for various radial and axial aspect ratios. The existence and validity of two-dimensional solutions were established either in the core region, the size of which decreases when increasing Ra, or for transient flows. The three-dimensional simulations performed for moderate and large gap annuli [START_REF] Dyko | A numerical and experimental investigation of stability of natural convective flows within a horizontal annulus[END_REF] present a good accordance with the prediction of the Powe chart. Some discrepancies are however observed for narrow annuli (R < 1.2) in which fully three-dimensional flows, due to the development of spiral vortex, were computed. Therefore, the flow exhibits in the upper region of the narrow-gap annulus either transverse rolls or a combination of longitudinal and transverse rolls with respect to the cylinder axis [START_REF] Dyko | Three-dimensional natural convective states in a narrow-gap horizontal annulus[END_REF]. For small enough radius ratio, an odd number of transverse rolls can occur so that a longitudinal flow is present at the mid axial plane [START_REF] Dyko | On the presence of odd transverse convective rolls in narrow-gap horizontal annuli[END_REF].

A three-dimensional stability analysis was performed by Choi and Kim [START_REF] Choi | Three-dimensional linear stability of natural convective flow between concentric horizontal cylinders[END_REF] using the linear theory and more recently by Dyko and al. [START_REF] Dyko | A numerical and experimental investigation of stability of natural convective flows within a horizontal annulus[END_REF] using both the linear theory and the energy method. The neutral stability curves as well as the critical wave numbers were computed as a function of the radius ratio.

The unstable mode superimposed on the basic flow gives rise to spiral motions similar to those observed in numerical or experimental simulations. As R tends to 1, Ra c approaches 1708, the critical Rayleigh number which corresponds to the onset of the natural convection in a horizontal layer.

Amongst the two-dimensional studies, a large variety of works was focused on multi-cellular flow patterns and multiplicity of solutions for various sets of parameters (Yoo [START_REF] Yoo | Dual steady solutions in natural convection between horizontal concentric cylinders[END_REF][START_REF] Yoo | Dual free-convective flows in a horizontal annulus with a constant heat flux wall[END_REF], Cadiou et al. [START_REF] Cadiou | Natural convection in a narrow horizontal annulus: the effects of thermal and hydrodynamic instabilities[END_REF][START_REF] Cadiou | Etude numérique de la stabilité des écoulements de convection naturelle dans des espaces annulaires horizontaux de faible épaisseur[END_REF], Desrayaud et al. [START_REF] Desrayaud | Thermoconvective instabilities in a narrow horizontal air-filled annulus[END_REF] and Mizushima et al. [START_REF] Mizushima | Exchange of instability modes for natural convection in a narrow horizontal annulus[END_REF]16]). These results reveal the existence of an imperfect bifurcation: the Rayleigh-Bénard pitchfork bifurcation is broken into a couple of one virtual transcritical and one saddle-node bifurcation from which two stable branches of multi-cellular solutions emerge. The mechanisms of formation or merging of the cells within the upper part of the annulus are indeed the same on the two branches and are due to thermal instabilities. The number of cells mainly depends on the radius ratio and Rayleigh number. For a fixed supercritical value of Ra, the number of cells increases as the radius ratio decreases. When R → 1 the number of cells tends toward infinity and the classical Rayleigh-Bénard problem characterized by a pitchfork bifurcation applies. The multi-cellular flows calculated in these studies was shown to undergo an unsteady secondary instability by increasing the Rayleigh number, provided the radius ratio is small enough (for example R = 1.14). The resulting periodic flow is composed of cellular patterns located into the basic flow in the lateral regions. This second type of instability is hydrodynamic in its origin as in air-filled vertical slots.

An overview of the published numerical studies (2D-as well as 3D-modeling) shows that very few considered a possible breaking of the vertical symmetry [START_REF] Petrone | Brisure de symétrie des écoulements de convection naturelle entre deux cylindres coaxiaux horizontaux[END_REF][START_REF] Petrone | Instabilities of two dimensional natural convection in a horizontal annulus[END_REF] because flow-symmetry was assumed in order to save CPU-time. As an example of such a symmetry breaking, Petrone et al. [START_REF] Petrone | Brisure de symétrie des écoulements de convection naturelle entre deux cylindres coaxiaux horizontaux[END_REF] computed steady and oscillatory dissymmetrical flows for low Prandtl number fluids. For air as working fluid and a fixed radial ratio R = 1.2 [START_REF] Petrone | Instabilities of two dimensional natural convection in a horizontal annulus[END_REF], the appearances of dissymmetrical disturbances give rise to unstable and non-symmetrical flows that had never been reported in the current literature.

Therefore, the present study is motivated by the need of revisiting bifurcation and flow structure maps in the (Ra, R) plane in order to unify the results published for two-dimensional natural convection in horizontal air-filled annuli. The paper is divided into two main parts. A set of numerical methods especially developped for linear stability analysis and bifurcation tracking are first presented and validated. The results are then discussed into two sections. The first section is devoted to the flow description and stability analysis for a fixed radius ratio R = 1.2. The effects of increases in R on the flow patterns and critical Rayleigh number are then investigated.

2 Equations and numerical methods

Governing equations

Consider two horizontal concentric cylinders of radii r * i and r * o = R × r * i > r * i held at uniform temperatures T * i and T * o < T * i , respectively (Fig. 1). The two-dimensional governing equations for laminar, incompressible flows of air with constant physical properties except the density in the buoyancy term are written in cylindrical coordinates and the Boussinesq approximation is invoked. By scaling the length by the gap width d and the time by α/d, by introducing the dimensionless temperature difference T = (T * -0.5

× (T * o + T * i ))/(T * i -T * o )
and the radial coordinate r = (r * -r * i )/d, the dimensionless governing equations in cylindrical coordinates (r, θ), where θ is measured from the downward vertical, read :

∇.u = 0 ∂ ∂t (f u) + ∇. (u.u) -ηv 2 = -f ∂p ∂r -RaP rf T cos θ + P r ∇ 2 u -2η 2 f ∂v ∂θ ∂ ∂t (f v) + ∇. (v.u) + ηuv = -f ∂p ∂θ + RaP rf T sin θ + P r ∇ 2 v + 2η 2 f ∂u ∂θ ∂ ∂t (f T ) + ∇. (T.u) = ∇ 2 T (1)
where η = R -1 and f = ηr + 1. The operators are defined as follows:

∇.(X.u) = (∂/∂r) (f X.u) + (∂/∂θ) (ηX.v) and ∇ 2 X = (∂/∂r)(f [∂X/∂r]) + (∂/∂θ)(η 2 /f × [∂X/∂θ])
, where X(r, θ) is a scalar quantity. The boundary conditions are periodic in the azimuthal direction and verify the following conditions in the radial direction:

0 ≤ θ ≤ 2π : r = 0      u = 0 T = +0.5 r = 1      u = 0 T = -0.5 (2)

Numerical methods

The governing equations (1) with boundary conditions (2) were solved numerically in the primitive variable formulation using a finite volume method on a non-uniform, staggered-structured grid. Each elementary mesh is defined by [r i , r i+1 ] × [θ j , θ j+1 ] with

r i = tanh c r 2 i-1 Nr -1 + tanh(c r ) 2 tanh(c r ) i = 1, N r + 1 θ j =          π exp 2πc θ j-1 N θ -1 exp(πc θ ) -1 j = 1, N θ 2 + 1 2π -θ N θ +2-j j = N θ 2 + 1, N θ + 1 (3)
where N r is an any integer and N θ an even integer, c r and c θ are two negative or zero real parameters used to modify the grid distribution in the upper part of the annulus (near θ = π) and close to the cylinder surfaces. The transport terms of the momentum and energy equations were discretized with a secondorder centered scheme. A time splitting method [START_REF] Karniadakis | High order splitting method for incompressible Navier-Stokes equations[END_REF] was used to uncouple the velocity and pressure fields.

Time-marching numerical schemes

To perform the linear stability analysis, several numerical methods have been implemented and are now described. In order to simplify the presentation, notations are introduced in cartesian coordinate system where, in contrast to cylindrical coordinate, no coupling of the radial and azimuthal velocity components is present in the diffusion terms of the momentum equations. The discrete velocity components and temperature are combined into one single vector S(t) = (u, T ) t . Let note δΩ the boundary of an elementary cell Ω and n the unit outwards vector, then the spatially discretized transport, diffusive and Boussinesq operators are noted T (u), D and B respectively so that:

T (u)S = Nr×N θ i=1    δΩ i    (u ⊗ u).n T u.n    ds    DS = Nr×N θ i=1    δΩ i    P r ∇u .n ∇T .n    ds    BS = Nr×N θ i=1    Ω i    P r (0, T ) t 0    dV   
By introducing the non-homogeneous part of the boundary conditions as vector bc, the set of the governing unsteady equations reads

∂ ∂t S = (T (u) + D + Ra B) S + bc (4)
The steady or unsteady flows are obtained by using a time integration of Eq. ( 4) with a first order scheme in which the diffusive terms are treated by an implicit method and the convective terms by an explicit method.

The linear stability analysis of a steady state S 0 consists in studying the growth rates of infinitesimal disturbances δS = (δu, δT ) t as a function of only one parameter : the bifurcation parameter. Therefore, the solution S 0 + δS satisfies Eq. ( 4):

∂ ∂t (S 0 + δS) = T (u 0 + δu) + D + Ra B (S 0 + δS) + bc that is equivalent to ∂ ∂t δS = (J[T (u 0 )] + D + Ra B) δS + O( δS 2 )
where J[T (u 0 )]δS is the Frechet derivative defined by J[T (u 0 )]δS = T (u 0 )δS + T (δu)S 0 . By neglecting the second-order terms, the solution of the previous equation reads:

δS(t) = exp (J[T (u 0 )] + D + Ra B) t δS(t = 0) (5)
When time tends towards infinity, the disturbance amplifies or is damped according to the sign of the eigenvalue of the matrix (J[T (u 0 )] + D + B) having the largest real part. The evaluation of the dominant spectrum with the largest real parts is achieved by an Arnoldi method [START_REF] Arnoldi | The principle of minimized iterations in the solution of the matrix eigenvalue problem[END_REF] as proposed by Mamun and Tuckerman [START_REF] Mamun | Asymmetry and Hopf bifurcation in spherical Couette flow[END_REF]: the evaluation of the right-hand side of Eq. ( 5) is performed by a first order temporal scheme in which the convective terms have been linearized in perturbations.

Non-linear systems

Although the time marching code allows to compute stable flows, its convergence rate tends towards zero when the parameters are close to the bifurcation thresholds. To evaluate efficiently steady flows near the transitions or the unstable steady solutions, the following non-linear system is solved:

(T (u 0 ) + D + Ra B) S 0 + bc = 0 (6)
using the Newton method. The k th step reads:

     J[T (u (k) 0 )] + D + Ra B δS (k+1) = T (u (k) 0 ) + D + Ra B S (k) 0 + bc S (k+1) 0 = S (k) 0 -δS (k+1) (7)
In order to compute an eigenvalue, a possible alternative to the Arnoldi method is to solve directly the non-linear eigen-problem [START_REF] Chénier | Stability of the axisymmetric buoyant-capillary flows in a laterally heated liquid bridge[END_REF]:

             (J[T (u 0 )] + D + Ra B) ψ -λψ = 0 c t .Re(ψ) = 1 c t .Im(ψ) = 0 (8)
where (λ, ψ) is the couple of complex eigenvalue with its corresponding eigenvector and c is an arbitrary real vector used to ensure the normalized condition for ψ: in this work, c ≡ S 0 . Therefore, the k th step of the Newton algorithm consists in computing:

                                 J[T (u (k) 0 )] + D + Ra B δψ (k+1) -λ (k) δψ (k+1) -δλ (k+1) ψ (k) = J[T (u (k) 0 )] + D + Ra B ψ (k) -λ (k) ψ (k) c t .Re(δψ (k+1) ) = c t .Re(ψ (k) ) -1 c t .Im(δψ (k+1) ) = c t .Im(ψ (k) ) ψ (k+1) = ψ (k) -δψ (k+1) λ (k+1) = λ (k) -δλ (k+1) (9)
To track the Hopf bifurcations as a function of a given parameter, Griewank and Reddien [START_REF] Griewank | The calculation of Hopf points by a direct method[END_REF] proposed to seek simultaneously the parameter threshold, the steady state S 0 , the complex eigenmode ψ and the corresponding eigenvalue λ = 0 + iω at the transition:

                   (T (u 0 ) + D + Ra B) S 0 + bc = 0 (J[T (u 0 )] + D + Ra B) ψ -iωψ = 0 c t .Re(ψ) = 1 c t .Im(ψ) = 0 (10)
In our problem, the bifurcation parameter is Ra so that the k th step of the Newton method reads:

                                                                       J[T (u (k) 0 )] + D + Ra B δS (k+1) + (δRa (k+1) B)S (k) 0 = T (u (k) 0 ) + D + Ra B S (k) 0 + bc J[T (u (k) 0 )] + D + Ra B δψ (k+1) + (δRa (k+1) B)ψ (k) -iω (k) δψ (k+1) -iδω (k+1) ψ (k) = J[T (u (k) 0 )] + D + Ra B ψ (k) -iω (k) ψ (k) c t .Re(δψ (k+1) ) = c t .Re(ψ (k) ) -1 c t .Im(δψ (k+1) ) = c t .Im(ψ (k) ) ψ (k+1) = ψ (k) -δψ (k+1) ω (k+1) = ω (k) -δω (k+1) Ra (k+1) = Ra (k) -δRa (k+1) (11)
A particular attention must be paid to the zero imaginary part of the eigenvalue. Indeed, equations ( 10) define a saddle-node bifurcation point [START_REF] Chénier | Stability of the axisymmetric buoyant-capillary flows in a laterally heated liquid bridge[END_REF]. The corresponding Newton algorithm is similar to Eq. ( 11) with Im(ψ

(k) ) = 0, ω (k)
= 0 as well as their corresponding increment Im(δψ (k+1) ) = 0 and

δω (k+1) = 0.
The linear systems arising from the Newton algorithm are iteratively solved either by the BiConjugate-Gradient Squared, as proposed in [START_REF] Joly | Preconditioned biconjugate gradient methods for numerical reservoir simulation[END_REF], or by the Generalized Minimal RESidual methods coming from the library [START_REF] Frayssé | A set of GMRES routines for real and complex arithmetics on high performance computers[END_REF]. The GMRES method is based on the minimization of the residual on a Krylov subspace. The main advantage and drawback of this method are highly linked: the residual decreases in a monotone manner as a function of the increase in the Krylov dimension that ensures convergence. However, the required space storage, proportional to the number of the base vectors defining the Krylov subspace, may increase dramatically. Hence, the BCGS method was preferred for solving the linear equations (Eqs. ( 7) and ( 9)) but failed for the Hopf system Eq. ( 11) which was solved by using the GMRES algorithm.

Implementation of the linear systems

The matrices of the linear systems (Eqs.7, 9, 11) are not explicitly built up because of their huge sizes. Each matrix vector product is obtained thanks to the method proposed by Mamun and Tuckerman [START_REF] Mamun | Asymmetry and Hopf bifurcation in spherical Couette flow[END_REF] for the steady state problem (Eq. ( 6)) and extended by Chénier et al. [START_REF] Chénier | Stability of the axisymmetric buoyant-capillary flows in a laterally heated liquid bridge[END_REF] for the eigenvalue problems (Eq. ( 8)) and the saddle-node and Hopf bifurcations (Eq. ( 10)). The principle of this method is based on the temporal code (Eq. ( 4)) discretized in time with a first order scheme as described in section 2.2.1. Let S (0) and S (1) be two consecutive iterations so that for a time step ∆t:

S (1) -S (0) ∆t = T (u (0) ) + Ra B S (0) + (D)S (1) + bc
The solution S (1) gives:

S (1) = 1 ∆t I -D -1 1 ∆t S (0) + T (u (0) ) + Ra B S (0) + bc therefore S (1) -S (0) = 1 ∆t I -D -1 T (u (0) ) + D + Ra B S (0) + bc ( 12
)
where I is the identity matrix. The idea of the method proposed by Mamun and Tuckerman [START_REF] Mamun | Asymmetry and Hopf bifurcation in spherical Couette flow[END_REF] is therefore to substitute Eq. ( 6) by its preconditioned expression:

-D -1 T (u (0) ) + D + Ra B S (0) + bc = 0 obtained by subtracting two consecutive iterations S (1) -S (0) for large enough time step ∆t (Eq. 12). All the matrix vector products in Eqs. [START_REF] Dyko | Three-dimensional natural convective states in a narrow-gap horizontal annulus[END_REF][START_REF] Choi | Three-dimensional linear stability of natural convective flow between concentric horizontal cylinders[END_REF][START_REF] Yoo | Dual free-convective flows in a horizontal annulus with a constant heat flux wall[END_REF] are then constructed on the same principle by adapting the initial temporal code (Eq. 4) to the problem to be resolved.

The use of the projection algorithm [START_REF] Karniadakis | High order splitting method for incompressible Navier-Stokes equations[END_REF] in the context of the Mamum and Tuckerman method [START_REF] Mamun | Asymmetry and Hopf bifurcation in spherical Couette flow[END_REF] requires some technical adaptations. Indeed, the choice of a large ∆t for the pressure calculation increases dramatically the number of iterations of the GMRES and BCGS methods [START_REF] Chénier | Etude de la stabilité linéaire des écoulements thermocapillaires et thermogravitationnels en croissance cristalline[END_REF]. To overcome this problem by keeping the Mamum and Tuckerman method, a second time step (or parameter) ∆t p is introduced by adding the

u (n) -u (n)
∆tp expression into the discretized Navier-Stokes equations. The projection step, which defines the pressure equation, is then constructed with ∆t p whereas the diffusion part uses the large time step ∆t. The degree of freedom coming from the parameter value ∆t p is therefore used to reduce as much as possible the number of iterations necessarily to converge with the GMRES and BCGS methods. Classicaly, the ∆t p value is in the range [10 -3 , 10 -1 ].

Additional numerical aspects

Although many algorithms allow to compute the same physical quantities (steady-states, eigenvalues or bifurcations), their respective use are often complementary. The Newton method converges by an iterative process and needs an appropriate initial guess. This initialization is therefore ensured either by a time marching code or by using previous results obtained with the Newton procedure. The Arnoldi method (5) allows to evaluate the leading eigenvalues without necessary using an initial guess. However, as the gap between the real parts of the interested and dominant eigenvalues increases, the computations are less accurate. Moreover, numerous iterations may be necessary to distinguish two modes having almost identical growth rates.

The following remarks can be brought:

• The evaluation of the critical Rayleigh number is of course sensitive to the mesh size used for computations. This dependency is linked to the slope s of the curve, representing the real part of the eigenvalue as a function of Ra. In the vicinity of the threshold value, Re(λ) evolves quite linearly as a function of Ra, and therefore an absolute error of order ε on the eigenvalue computation induces a threshold displacement of the order of ε/s. So, if the slope is very small, the error on the critical parameter may be highly amplified.

• Some convergence problems of the Newton algorithm may be observed when the Jacobian matrix of the non-linear application has too small eigenvalues or is not invertible. In that situation, theorical studies suggest a linear (instead of a quadratic) convergence of the Newton method. This behaviour should probably be observed if the iterative resolution of the linear system at each Newton step was exact. For the present computations, the BCGS and GMRES methods failed to converge for a small enough stopping criteria: a too large stopping criteria does not improve enough the Newton steps, so that the global convergence is finally not ensured.

Validation

Comparisons with the literature

In order to validate the numerical schemes, several usual computations have been performed and then compared with the literature. For Ra = 3000, P r = 0.7, R = 1.6 and a 60 × 240 mesh, the steady flow obtained with the time marching code (Eq. ( 4)) on the full cavity consists in two crescent-shape cells, symmetrically located within the annulus. As shown in table [START_REF] Powe | Free convective flow patterns in cylindrical annuli[END_REF], the average Nusselt numbers at the cylinder surfaces are in good agreement with previous works. No significant differences were obtained between the values computed by the Newton (Eq. ( 7)) or by the temporal (Eq. ( 4)) codes.

The growth rates of the dominant linear perturbation associated with the symmetrical two-cellular flow at R = 1.2 and P r = 0.7 and achieved with the Newton method (Eq. ( 9)) are in accordance with the data obtained by the spectral method used by Mizushima et al. [16] (Tab. 2). For a fixed 60 × 240 mesh, the relative difference between the results does not exceed 1.5 %, except for small growth rates where the error reaches 4 %.

Comparisons were also made about the location and nature of the first transition of the symmetrical two-cellular flow at R = 1.2 and P r = 0.7 (Tab. 3).

Although the transition is differently named according to the authors, it deals with the same physical process. The present Rayleigh number values at the thresholds are close to the previous works, with a maximum relative difference less than 1.5 %. The virtual transcritical bifurcation threshold, so-called by Mizushima et al. [16], has been obtained by interpolating the eigenvalues as a function of Ra so that the growth rate is zero whereas the saddle-node bifurcation is directly computed thanks to the Newton algorithm (Eq. ( 11)).

Comparisons between the numerical methods and convergence properties

Comparisons are discussed here about the use of the different methods presented in this paper. At first, the focus is on the growth rate calculations carried out by the Arnoldi and Newton schemes. The Arnoldi method is based on a first order time marching code that implies a linear convergence of solutions as a function of the time step (Fig. 2). The linear extrapolation based on the six most accurate eigenvalues obtained using the Arnoldi method for a zero time-step (lim ∆t→0 (λ) = -2.4391) presents a good agreement with the value achieved with the Newton procedure (λ = -2.4330), with a relative discrepancy less than 0.5 %. The evaluation of the Hopf bifurcation point for R = 1.3 was then compared. It can be performed, either by computing iteratively the steady (Eq. ( 6)) and eigenvalue (Eq. ( 8)) problems for successive Rayleigh numbers (Fig. 3), or by simultaneously solving the steady and eigen-problems at the bifurcation point (Eq. ( 10)). In both cases, the critical threshold and pulsation are almost identical, Ra ≈ 6091 and Im(λ) ≈ 12.72 when using a 60 × 240 mesh. The threshold Ra c (Fig. 4a) and the pulsation Im(λ) (Fig. 4b) of the Hopf bifurcation at R = 1.3 was also studied as a function of the average mesh size (h 2 r = 1/N 2 r , h 2 θ = (2π/N θ ) 2 ), and for several grid distributions parametrized by c r and c θ . The numerical results are plotted in these figures by symbols as well as by continuous lines which correspond to linear polynomial fits as a function of the average squared grid size. Each grid distribution c r and c θ is associated to a particular symbol. The lower and upper abscissa axis are labelled by the cell numbers N r and N θ respectively. The set of four (four filled symbols) and three curves (four open symbols), which converge respectively towards (Ra c ≈ 6034, Im(λ) ≈ 12.48) with h r → 0 and (Ra c ≈ 6188, Im(λ) ≈ 13.13) with h θ → 0, correspond to the influence of the average grid sizes h r and h θ , for fixed meshes in the second direction, namely (N θ = 240, c θ = -0.75) and (N r = 60, c r = 0). Therefore, the best radial mesh with N θ = 240 (filled symbols) is obtained for c r = 0, i.e. a regular mesh, and it is preferable to employ c θ = -1.00 when N r = 60 (open symbols). In this last case, decreasing c θ from -1.00 to -1.25 does not improve the convergence rate because the computed Ra c and Im(λ) are indeed quite identical. Whereas it is worthwhile to reduce the mesh size in the azimuthal direction in order to improve the spatial resolution in the upper part of the annulus, that does not seem useful close to the radial boundaries. The last two curves in Figs. 4(a) and 4(b) which tend to (Ra c ≈ 6152, Im(λ) ≈ 13.00) are achieved by increasing proportionally the cell numbers in both directions. By considering relative differences between the extrapolated and computed values less than 5 % and by taking into account the computation time and the memory storage requirements, a (N r , N θ ) = (60, 240) mesh with (c r , c θ ) = (0, -1) was retained to compute the Hopf bifurcations as a function of the radius ratio. With this particular choice at R = 1.3, the threshold and pulsation are equal to (Ra c = 6091, Im(λ) = 12.72) and the relative discrepancies are about 1 % and 2 %, respectively.

Results

It can be easily proved that if (u 0 , v 0 , T 0 ) is a solution of the coupled equations (1) and boundary conditions (2), then (u 1 , v 1 , T 1 ) = Γ(u 0 , v 0 , T 0 ) defined by Γ :

             Γ(u(r, θ)) = u(r, -θ) Γ(v(r, θ)) = -v(r, -θ) Γ(T (r, θ)) = T (r, -θ) (13)
is also a solution. A particular and interesting situation occurs when Γ(u 0 , v 0 , T 0 ) = (u 0 , v 0 , T 0 ). In that case, the flow is itself symmetric through the transformation Γ, noted Γ-symmetry in what follows, and the computation of such a solution requires only to solve equations (1) with the following symmetrical boundary conditions:

0 ≤ θ ≤ π : r = 0      u = 0 T = +0.5 r = 1      u = 0 T = -0.5 (14) and 0 ≤ r ≤ 1 : θ = 0, π                ∂u ∂θ = 0 v = 0 ∂T ∂θ = 0 (15)
This simplified mathematical model labelled π-annulus has been extensively employed, for instance in [START_REF] Yoo | Dual steady solutions in natural convection between horizontal concentric cylinders[END_REF][START_REF] Mizushima | Exchange of instability modes for natural convection in a narrow horizontal annulus[END_REF]16,[START_REF] Cheddadi | Twodimensional convective bifurcation in a horizontal annulus[END_REF][START_REF] Kim | Numerical investigation on bifurcative natural convection in an air-filled horizontal annulus[END_REF][START_REF] Chung | Numerical investigation on the bifurcative natural convection in a horizontal concentric annulus[END_REF], in order to decrease both the CPU time and memory space storage.

The results are organized into two sections. The flows and the bifurcation diagram are first presented for a narrow cavity with R = 1.2. Secondly, the bifurcations are tracked as a function of the radius ratio.

3.1 Radius ratio R=1.2

Bifurcation diagram

The results are presented in the form of a bifurcation diagram representing the steady radial velocity component u(0.5, π) as a function of Ra (Fig. 5). This component is useful in determining quickly the direction of the rotating cells in the upper region of the annulus when the flow is Γ-symmetric. For a given Rayleigh number, the steady states lying on the continuous and dashed lines are stable and unstable, respectively. Below Ra l s 1 = 1911, only the well known basic couple of crescent-shaped cells, symmetrically located on both side of the annulus are observed. Just above Ra l s 1 , two new steady solutions are found lying on the unstable branch as shown in the enlargement plotted in Fig. 5. For convenience, the branches starting from Ra = 0 and Ra l s 1 will be called in what follows as the basic-branch and lower-branch, respectively: the quantities expressed on those branches are labeled ". b " and ". l ". The sharp growth of u(0.5, π) as well as the appearance of the lower-branch of solutions are due to thermal instabilities similar to those observed in Rayleigh-Bnard experiments.

The effect of the curvature is to break the classical pitchfork bifurcation into a saddle-node bifurcation at Ra l s 1 = 1911 and into a transcritical-like bifurcation at Ra t = 1917, values in excellent agreement with those reported in [16] (see Tab. 3). However, the comparison of the enlargement with [16] shows differences on the stability solutions lying on the lower part of the lower-branch, below u(0.5, π) ≈ 0. Indeed, unlike to the Mizushima et al. results which leads to stable steady states, our solutions are unstable. An unstable disturbance is still present above Ra l s 1 (Fig. 6). By describing the curves through decreases in Ra, the eigenvalue λ l 2 which gives rise to the saddle-node bifurcation, decreases sharply from zero whereas the growth rate λ l 1 of a new disturbance slowly decreases and reaches zero at Ra l f 1 . It is worth noting that the slope of λ l 1 as a function of Ra is of order 10 -4 and therefore Ra l f 1 is rather sensitive to numerical errors (see section 2.2.4). The streamlines of these two eigenmodes, plotted in Fig. 7, are characterized by two different symmetries. The mode giving rise to the saddle-node bifurcation is associated to the Γ-symmetry whereas the other eigenvector satisfies the ∆-symmetry defined by ∆ :

             ∆(u(r, θ)) = -u(r, -θ) ∆(v(r, θ)) = v(r, -θ) ∆(T (r, θ)) = -T (r, -θ). (16)
At the bifurcation point Ra l f 1 , a new branch of solutions thus appears (Fig. 5). This branch is locally tangent to the central manifold: the steady states lying on that new branch are then qualitatively made of the Γ-symmetrical basic flow superimposed to the ∆-symmetrical perturbation. Consequently, the Γ-symmetry is broken. Since the governing equations always verify the Γ-symmetry property, another perturbation is also possible. It corresponds to the disturbance obtained through the Γ-symmetry applied to the perturbation drawn in Fig. 7a. It should be noted that the branch based on this second disturbance is not distinguishable on the bifurcation diagram (Fig. 5) because u(0.5, π) is unchanged through the Γ-symmetry. Thanks to all these considerations and to the stability of solutions lying on this new branch, it is obvious that the critical threshold Ra l f 1 corresponds to a sub-critical pitchfork bifurcation. Above Ra l f 1 , the steady states are stable up to Ra = 5000.

The dominant spectrum computed for flows lying on the basic-branch is shown in Fig. 8. For Ra < Ra t = 1917, the two real eigenvalues λ b 1 and λ b 2 grow up simultaneously to almost zero values. It should be noted that the part of the curve λ b 2 vs Ra for Ra < Ra t is in a good agreement with the results of Mizushima et al. [16]. On the other hand, the second part differs because it was computed in [16] for solutions lying on the upper part of the lower-branch. Therefore their eigenvalues λ 2 are positive: the branch jumping defines the virtual-transcritical bifurcation. Close to Ra t (Fig. 8) both eigenvalues sepa-rate so that λ b 2 decreases sharply while λ b 1 crosses the zero-axis, first from negative to positive values at Ra b f 1 = 2068, and then after in the reverse direction at Ra b f 2 = 2456. These transitions correspond to two pitchfork bifurcations (Fig. 5) and are associated to one ∆-symmetric disturbance looking like the one presented in Fig. 7a. As for the pitchfork bifurcation of the lower-branch of the bifurcation diagram, these two transitions give rise to new solutions breaking the Γ-symmetry. In particular, the branch starting from the Ra b f 1 threshold and made of unstable dissymmetrical solutions (Fig. 9), joins the pitchfork bifurcation defined at Ra l f 1 . For Ra > Ra b f 2 , the flow is found stable until Ra = 5000.

The solutions lying on the almost horizontal part of the lower-branch (see enlargement in Fig. 5) defined for Ra > Ra t are two times unstable. At Ra ≈ 2435 (Fig. 5), an imperfect pitchfork bifurcation occurs which seem to be made of both one transcritical (Ra ≈ 2435) and one saddle-node (Ra s 2 ≈ 2431) bifurcations (Fig. 10a). An appropriate scaling of Fig. 10a clearly shows two couples of saddle-node and pitchfork bifurcations close to Ra f s ≈ 2434.45 and Ra ≈ 2434.95 (Fig. 10b). By following the upper curve of Fig. 10a, another pitchfork appears which gives rise to dissymmetrical solutions on a new branch at Ra f 2 ≈ 2950 before reaching the basic-branch of the bifurcation diagram (Fig. 5) at Ra b f 2 through one saddle-node at Ra s 3 ≈ 3510. It is therefore possible to cross continuously from the lower to the basic-branch of the bifurcation diagram by following a set of curves. In what follows, each part of the continuation curve are labelled by "uns i " to denote that solutions are unstable to i disturbances with positive growth rates. Therefore, it exists two different paths to join the bifurcation points Ra l f and Ra u f 2 . The easiest one is constituted of three pitchfork bifurcations and described by the way Ra l

f 1 uns 1 → Ra b f 1 uns 1 → Ra b f 2 .
Both four saddle-nodes and pitchfork bifurcations take place on the much more complicated second path summed up by

Ra l f 1 uns 1 → Ra l s 1 uns 2 → Ra f s uns 4 → Ra s 2 uns 3 → Ra f 2 uns 2 → Ra s 3 uns 1 → Ra b f 2 .
Although the pitchfork bifurcations described in Fig. 10b occur in a very narrow gap of Ra, and could be considered as spurious, they are in fact linked to the appearance of the two other transitions at Ra s 2 and Ra f 2 .

Additional comments should also be added on the symmetry breaking of the perturbations. Indeed, this results could be surprisingly new since numerous works having dealt with this problem. In fact, the use of Γ-symmetrical models are probably only one part of the explanation. Several computations on the full annulus, with periodic conditions in the azimuthal direction, have already been performed without showing such instability. The main reason is probably given through examination of the growth rate values of the ∆-symmetrical disturbances represented in Figs. 6 and8. These values are so small that convergence with a temporal code may easily be achieved without observing the destabilizing effects of this perturbation. Because the steady solutions used to initialize the temporal codes are often characterized by the Γ-symmetry, the initial magnitude of the ∆-disturbance is of the order of the numerical errors. For instance, in the linear point of view and for a growth rate λ ≈ 0.025 (Fig. 8), a dimensionless time of 1/λ ≈ 40 is required to multiply the magnitude of the disturbance by only e 1 . This behavior is not new and explains the need of using suitable disturbances to evolve from conductive to convective flows in the classical Rayleigh-Bénard flow simulations.

Flow fields

For Ra < Ra t , the flow pattern labelled C + consists in two large cells on each part of the annulus due to the Γ-symmetry. Since only Γ-symmetrical flows are described, the flow patterns are plotted in a half region of the annulus. By following the basic-branch of the bifurcation diagram by increasing Ra (Fig. 5), the flow evolves continuously from one large crescent-cell C + to three cells, two cells being located at the top of the annulus and the other forms the crescent-shaped cell. In order to investigate the birth and growth of the upper cells, their angular sizes at the middle between the two cylinders are measured as a function of Ra (Fig. 11). For convenience, the cells rotating in the same direction as the large crescent-cell C + are also noted C + or c + , whether they are stretching out on the entire annulus gap or not. The counterrotating cells are labelled C -. Slightly above Ra t , a co-rotating cell c + grows very fast with Ra but stays confined in the large C + (Fig. 12a). The splitting of the daughter cell from the mother one occurs at Ra ≈ 2270 and gives rise to one counter-rotating cell C -located between them (Fig. 12b). From this point and until Ra = 5000, the steady flow remains three-cellular; this pattern is therefore denoted C + C -C + . It should be noticed that globally the angular size of the largest C + decreases in favor of the two upper eddies. The Rayleigh number value which corresponds to the onset of the C + C -C + pattern is in very good agreement with Kim and Ro [START_REF] Kim | Numerical investigation on bifurcative natural convection in an air-filled horizontal annulus[END_REF] and Chung et al. [START_REF] Chung | Numerical investigation on the bifurcative natural convection in a horizontal concentric annulus[END_REF] who proposed Ra = 2270 and Ra = 2330, respectively, as the threshold values for the flow transition between 1 and 3 cells. On the lower part of the continuation diagram (Fig. 5), a reverse flow appears just after the saddle-node bifurcation point at Ra l s 1 . Above this point, the negative radial velocity u(0.5, π) (see enlargement in Fig. 5) indicates the formation of a counter-rotating cell C - adjacent to the main crescent-shaped flow (Fig. 13). Like on the basic-branch of the bifurcation diagram, a co-rotating cell c + is set up slightly before Ra l f 1 but disappears at Ra ≈ 3500 (Fig. 14). The C -C + → C -c + C + → C -C + transition is therefore observed by increasing Ra on this branch. Some of the flow pattern modifications can be correlated with the unstable regions. The generation of the co-rotating nucleus c + on the basic-branch occurs when the growth rate of the disturbance is becoming positive, while the C + C -C + establishment precedes the flow stabilization. On the lower-branch, the onset of the re-circulation nucleus c + is detected at the location which characterises the re-stabilization of the steady flow.

A slow increase in Ra from zero shows that the steady flow pattern evolves from the crescent-shaped cell C + at first to c + C + and then to C -c + C + structures when the solutions reach the lower-branch as it is presented on the diagram in Fig. 15. The cell nucleus then disappears and the two-cell flow structure C -C + is the only possible. It is worth noting that the stable steady state tracking does not allow to reach the three-cell flows C + C -C + . To end, Figure 15 shows a hysteretical behavior between patterns lying on the lower and on the basic branches at Ra l f 1 and Ra b f 1 .

Effect of the radius ratio: bifurcation maps

The effects of the aspect ratio on the flow stability are discussed in this paragraph for 1.2 ≤ R ≤ 3 and Ra ≤ 10 4 . Increases in R lead to displacements of the bifurcation thresholds and also modify the shape of the bifurcation diagram, especially the basic-branch. Consequently, we will first focus on the fundamental changes occurring on the bifurcation diagram for small radius ratio R ∈]1.24, 1.28[. The critical values will be then presented in the (R, Ra)plane.

The basic and other unstable branches of the bifurcation diagrams for R = 1. u ", appears for Ra ≥ Ra u s (Fig. 17). The origin of the upper-branch is related to another saddle-node bifurcation. Indeed, by increasing R, the upper part of the branch b 2 (R = 1.25) and the basic branch b 1 (R = 1.25) approach themselves, merge together for R c ∈]1.25, 1.26[, and finally separate from each other after having exchanged parts of their branches. Consequently, no solution exists above R c in the vicinity of the previous contact point. Finally, the two saddle-node bifurcations belonging to the basic-branch at Ra s 1 (b 1 ) and Ra s 2 (b 1 ) collapse in the R-interval ]1.26, 1.28[. Increases in R upon a small interval has thus generated a new independent branch and modified substantially the solutions lying on the basic-branch.

Concerning the effects produced by increasing the radius ratio on the flow patterns corresponding to the different branches, no qualitative difference is found for 1.2 ≤ R ≤ 1.25: the flow structures are similar to those at R = 1.2. On the other hand, for 1.26 ≤ R < 1.35, the C + C -C + pattern is not observed on the basic-branch which is now characterized by c + C + structures: a c + nucleus appears, grows and then disappears. The C + C -C + is located on the upper-branch because this new branch takes its origin from a part of the previous basic-branch. For R ≥ 1.35 the fluid flows on the basic-branch are mono-cellular in the range of Ra investigated, while the upper-branch consists systematically in three-cellular patterns. Concerning the lower-branch of the bifurcation diagram, the C -c + C + flow is observed for R < 1.24 only. The c + nucleus size becomes smaller and smaller by increasing R and then the nucleus disappears. For larger radius ratios, the flow is two-cellular.

The bifurcation thresholds and the stability regions are presented in the following as two-dimensional diagrams plotted in the (Ra, R)-plane for the lower, basic and upper-branches. The lower-branch, starting from the saddle-node Ra l s 1 is plotted in Fig. 18. The variation of Ra l s 1 as a function of R is found in good agreement with the results of by Mizushima et al. [16] since the largest relative difference is only 0.25 % for R ∈ [1.2, 3]. For small radius ratios (R ≤ 1.5), the critical parameters of the pitchfork and saddle-node bifurcations are very close. Above R = 1.6, the pitchfork threshold Ra l f 1 abruptly increases to finally collapses at R ≈ 2 with another pitchfork bifurcation appearing on the same branch , but at Ra l f 2 . The stability region of solutions is therefore bounded by the two curves giving the pitchfork bifurcation thresholds. Indeed, for R ≥ 2.1 and Ra ≤ 10 4 , any stable flow occurs on this branch as it is checked by the study of the growth rate of the dominant perturbations. The small region of this branch characterized by the C -c + C + flow pattern is visible in the bottom left side of Fig. 18 bounded by dashed lines. On the basic-branch (Fig. 19) and for R ∈ merge between R = 1.25 R = 1.3 so that the associated unstable region disappears: the flow is therefore stable until the onset of the oscillatory disturbance at Ra b h . Finally the Ra-domain of instability due to unsteady perturbations is reduced when increasing R, and disappears past R ≈ 1.32. Consequently, the solutions may be continuously followed from Ra = 0 to Ra = 10 4 without meeting unstable regions. The upper-branch, starting from the saddle-node at Ra u s , consists in unstable solutions until the pitchfork bifurcation at Ra u f 2 occurs as it is shown in Fig. 20. It should be noted that superimposing Fig. 20 and Fig. 19, the two curves representing the bifurcation thresholds Ra b f 2 and Ra u f 2 match well and form a coherent continuous curve. As the radius ratio increases, the threshold Ra u f 2 moves quickly along the upper-branch up to Ra = 10 4 .

From the bifurcation maps shown in Figs. [START_REF] Petrone | Instabilities of two dimensional natural convection in a horizontal annulus[END_REF], 19 and 20 it can be concluded that various stable solutions and flow patterns set in for a given couple (R, Ra). The locations of the flow patterns are gathered in Fig. 21. As can be seen, three, two or only one solutions co-exist with a number of possible flows which decreases when the radius ratio increases.

Conclusion

Two-dimensional natural convection in air-filled horizontal annuli has been investigated numerically and linear stability analyses are presented for radius ratio in the range 1.2 ≤ R ≤ 3 and for Ra ≤ 10 4 . The governing equations were integrated by a classical time-marching procedure or directly solved in their non-linear steady form by using the Newton algorithm. The linear stability analysis was conducted by investigating the growth rates of infinitesimal disturbances about the steady states. The resulting eigen-problems were solved either by the Arnoldi method or by the Newton method. Solutions of large non-linear systems of algebraic equations were also computed in order to properly evaluate the Hopf and saddle-node bifurcations when using the Newton method. Validations of the numerical tools and comparisons with previously published results are presented and are shown in good agreement. From the bifurcation diagrams and bifurcation maps in (R, Ra)-planes, complex flow behaviors are pointed out, especially for small radius ratios where solutions are highly modified through small increases in Ra. New pitchfork bifurcation thresholds are determined thanks to the full annulus model adopted in the computations. These transitions are proved to be produced by perturbations breaking the symmetry about the vertical plane containing the axes of the cylinders. Flow patterns were discussed and a map of the possible flow structures is presented in the (R, Ra)-plane. For the first time, multiplicity of solutions and flow patterns have been clearly established and described. G. Fig. 6. Dominant eigenvalues for flows lying on the lower part of the lower-branch, below u(0.5, π) ≈ 0 (Fig. 5).

G. Petrone et al. G. Petrone et al. 10 Enlargements of the bifurcation diagram (Fig. 5). 32

✲ ❄ ✻ ❄ ✻ ✻ ❄ C + C -C + C -C + Ra C + c + C + C -c + C + Ra b f 1 Ra l f 1 Ra b f 2
11 Angular sizes of cells for flows lying on the basic-branch of the bifurcation diagram (Fig. 5). 33

12 Streamlines of steady states lying on the basic-branch of the bifurcation diagram (Fig. 5). 34

13 Angular sizes of cells for flows lying on the lower-branch of the bifurcation diagram (Fig. 5). 35

14 Streamlines of steady states lying on the lower-branch of the bifurcation diagram (Fig. 5). 36

15 Flow patterns as a function of Ra. 37 18 Bifurcation thresholds and stability regions in the plane (Ra, R) for the lower-branch of the bifurcation diagram (Fig. 5). 40

19 Bifurcation thresholds and stability regions in the plane (Ra, R) for the basic-branch of the bifurcation diagram (Fig. 5). 41 20 Bifurcation thresholds and stability regions in the plane (Ra, R) for the upper-branch of the bifurcation diagram (Fig. 5). 42

21 Possible flow patterns in the (R, Ra) plane. 43

  24 and R = 1.25 are plotted in Fig.16. Such a small increase in R brings closer the branches labelled b 2 (R = 1.24) and b 3 (R = 1.24) so that they come into contact, very close to or at the saddle-node bifurcation points defined by Ra s (b 2 , R = 1.24) and Ra s (b 3 , R = 1.24), exchanging parts of their curves each other, and forming two new branches labelled b 2 (R = 1.25) and b 3 (R = 1.25). By crossing the critical radius ratio R c ∈]1.24, 1.25[, the number of solutions in the interval Ra ∈]Ra s (b 2 , R), Ra s (b 3 , R)[ changes from zero to two. This finding characterizes a saddle-node transition for the bifurcation parameter R. If the contact points are exactly located at the thresholds Ra s (b 2 , R = 1.24) and Ra s (b 3 , R = 1.24), a bifurcation with a co-dimension two is therefore expected in the (Ra, R) plane at Ra s (b 2 , R c ) = Ra s (b 3 , R c ). The most important changes for the flow stability occur in the R-interval ]1.25, 1.26[ where a new independent branch of solutions, called as the upper-branch and labelled with ".

  [1.2, 1.25], the pitchfork bifurcation thresholds Ra b f 1 and Ra b f 2 move so that the instability region grows slowly. The changes in shape of the bifurcation diagram act significantly on the flow stability because they occur in the instability region bounded by the pitchfork transitions. It can be thus concluded that the flow stabilizes not on the basic-branch at Ra b f 2 but rather on the upper-branch at Ra u f 2 . From Ra b f 1 , the flows lying on the basic-branch are unstable, firstly because of the disturbance associated to the pitchfork bifurcation, and secondly because the onset of an oscillatory perturbation gives rise to a Hopf bifurcation located at Ra b h . This last transition moves towards by increasing Ra so that a pitchfork bifurcation at Ra b f 3 appears just above R = 1.28. Thus, for R ∈]1.28, 1.29], a region of stable flows is distinguishable if Ra ∈]Ra b f 3 , Ra b h [. At last, the two pitchfork transitions at Ra b f 3 and Ra b f 1

Fig. 1 .

 1 Fig. 1. Geometry of the annulus.

Fig. 3 .

 3 Fig. 3. Eigenvalue as a function of Ra. Real (left axis) and imaginary (right axis) parts of the eigenvalues for R = 1.3.

Fig. 4 .

 4 Fig. 4. Influence of the average mesh size, labelled by (N r , N θ ), and distribution parameters (c r , c θ ) on (a) the critical Rayleigh number Ra c and (b) the pulsation Im(λ) for R = 1.3.

Fig. 5 .

 5 Fig. 5. Bifurcation diagram and enlargement for Ra ∈ [1900 ; 1950]: u(0.5, π) versus Ra. Continuous (resp. dashed) lines stand for branches of stable (resp. unstable) solutions.

Fig. 7 . 29 - 2 Fig. 8 .Fig. 9 .

 729289 Fig. 7. Streamlines of perturbations associated to the bifurcations.

Fig. 10 .

 10 Fig.10. Enlargements of the bifurcation diagram (Fig.5).

Fig. 11 .

 11 Fig. 11. Angular sizes of cells for flows lying on the basic-branch of the bifurcation diagram (Fig. 5).

Fig. 12 .

 12 Fig.[START_REF] Cadiou | Natural convection in a narrow horizontal annulus: the effects of thermal and hydrodynamic instabilities[END_REF]. Streamlines of steady states lying on the basic-branch of the bifurcation diagram (Fig.5).

Fig. 13 .

 13 Fig.[START_REF] Cadiou | Etude numérique de la stabilité des écoulements de convection naturelle dans des espaces annulaires horizontaux de faible épaisseur[END_REF]. Angular sizes of cells for flows lying on the lower-branch of the bifurcation diagram (Fig.5).

Fig. 14 .

 14 Fig. 14. Streamlines of steady states lying on the lower-branch of the bifurcation diagram (Fig. 5).

❄

Fig. 15 .

 15 Fig. 15. Flow patterns as a function of Ra.

Fig. 16 .

 16 Fig. 16. Elements of the bifurcation diagrams at R = 1.24 and R = 1.25. The arrows indicate the part of the curves coming closer with the increasing radius ratio.

Fig. 17 .

 17 Fig. 17. Elements of the bifurcation diagrams at R = 1.25 and R = 1.26. The arrows indicate the part of the curves coming closer with the increasing radius ratio.

Fig. 18 .Fig. 19 .

 1819 Fig.[START_REF] Petrone | Instabilities of two dimensional natural convection in a horizontal annulus[END_REF]. Bifurcation thresholds and stability regions in the plane (Ra, R) for the lower-branch of the bifurcation diagram (Fig.5).

Fig. 20 .

 20 Fig.[START_REF] Arnoldi | The principle of minimized iterations in the solution of the matrix eigenvalue problem[END_REF]. Bifurcation thresholds and stability regions in the plane (Ra, R) for the upper-branch of the bifurcation diagram (Fig.5).

Fig. 21 . 24 3 25 4 26 5 27 6 28 7 29 8

 21242526272829 Fig. 21. Possible flow patterns in the (R, Ra) plane.

47 16

 16 Elements of the bifurcation diagrams at R = 1.24 and R = 1.25. The arrows indicate the part of the curves coming closer with the increasing radius ratio. 38 17 Elements of the bifurcation diagrams at R = 1.25 and R = 1.26. The arrows indicate the part of the curves coming closer with the increasing radius ratio. 39

  Petrone et al.

		-2.32 0.015
			Re(λ)
			12.8
		-2.34 0.01
			Im(λ)
		0.005
		-2.36
		0
	λ Re(λ)	-2.38	12.7	Im(λ)
		-0.005
		-2.4
		-0.01	Numerical results
		-2.42	Linear fit
		-0.015	12.6
		-2.44 -0.02 6060 0	0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 0.0045 6070 6080 6090 6100 6110 δt Ra

Fig. 2. Eigenvalues as a function of time step δt for Ra = 1500, P r = 0.7 and R = 1.2.

G. Petrone et al.

 [START_REF] Cheddadi | Twodimensional convective bifurcation in a horizontal annulus[END_REF][28] [START_REF] Chung | Numerical investigation on the bifurcative natural convection in a horizontal concentric annulus[END_REF] Present result N u (N r × N θ ) Comparisons of the average Nusselt number N u = 1/(π) π 0 N u dθ on the radial boundaries in a half cavity for Ra = 3000, P r = 0.7 and R = 1.6. 2 Comparisons of eigenvalues for R = 1.2 and P r = 0.7 as a function of Ra.

G. Petrone et al.

G. Petrone et al.

List of Tables

1

Comparisons of the average Nusselt number N u = 1/(π) π 0 N u dθ on the radial boundaries in a half cavity for Ra = 3000, P r = 0.7 and R = 1.6.