Global stability analysis of the natural convection between two horizontal concentric cylinders

International Journal of Heat and Mass Transfer 172 (2021) 121151 Contents lists available at ScienceDirect International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/hmt Global stability analysis of the natural convection between two horizontal concentric cylinders J.J. Serrano-Aguilera a , ∗, Francisco J. Blanco-Rodríguez b , L. Parras a a Universidad de Málaga, Escuela de Ingenierías Industriales, Campus de Teatinos s/n, 29071, Málaga, Spain b Área de Mecánica de Fluidos, Departamento de Ingeniería Aeroespacial y Mecánica de Fluidos, Universidad de Sevilla, Avenida de los Descubrimientos s/n, 41092, Sevilla, Spain a r t i c l e i n f o Article history: Received 16 September 2020 Revised 19 February 2021 Accepted 20 February 2021 Keywords: Natural convection Spectral methods Boussinesq approximation Buoyancy–driven flow Solar thermal energy a b s t r a c t In this investigation, the 2D flow between two horizontally positioned concentric cylinders (gravity perpendicular to the axis of the cylinders), where the inner cylinder is kept at constant temperature T i higher than the outer border temperature T o , is analyzed. Buoyancy forces initiate the movement of the fluid and the generated flow is studied in a fixed geometry for values of Prandtl numbers ( Pr) between 0.01 and 1, and Rayleigh numbers ( Ra ) between 10 2 and 5 ·10 6 . To solve the problem, a Chebyshev-Fourier spectral code is developed in polar coordinates (r, θ) respectively, and a complete map of steady-state solutions is obtained where regions with multiple solutions are identified. Later, a global stability study of the obtained stationary solutions is carried out, providing a transition curve to unstable areas as a function of the control parameters of the problem (Pr, Ra ) . Finally, to check the stability results, temporal evolution simulations are accomplished for several cases where dual solutions are presented, finding intermediate almost stationary solutions, and demonstrating that there exist typically single oscillating plume or double oscillating plume solutions (depending on the parameter space), where some of them have higher heat transfer coefficients, which may be interesting alternatives to improve heat exchange systems by means of passive control techniques. ©2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) 1. Introduction Convection is a major feature of the dynamics of the oceans, the atmosphere, and the interior of stars and planets [1,2] . Natural convection systems often exhibit surprisingly regular large scale patterns. Examples are the solar granulation and mesoscale convection in the Earth’s atmosphere as visualized by the cloud patterns seen in satellite images. Evidently, the presence of large-scale regular structures at extremely large Reynolds and Rayleigh numbers, often presumed in these systems to be asymptotically large, is compatible with a broad spectrum of highly chaotic motions at intermediate scales. Unfortunately, this aspect of turbulent convection cannot be easily investigated in laboratory experiments. Fluid motion driven by thermal gradients (thermal convection) is a common and important phenomenon in nature. Due to its theoretical interest and its various engineering applications such as ∗Corresponding author. E-mail addresses: jj.serr[email protected] (J.J. Serrano-Aguilera), [email protected] (F.J. Blanco-Rodríguez), [email protected] (L. Parras). heat exchangers, thermal energy storage systems, cooling of electronic components, transmission cables and solar thermal energy [3] , natural convection heat transfer in a horizontal cylindrical annulus kept at constant surface temperature has been the subject of interest of many researchers. Its simple geometry and well-defined boundary conditions allow to characterize the phenomenon by two dimensionless parameters: the Prandtl number ( P r) which is the ratio of viscous diffusivity to thermal diffusivity and the Rayleigh number ( Ra ) which can be seen as the ratio of the gravitational potential energy to the energy due to viscous dissipation and thermal diffusion. The aspect ratio A ≡D i /L which is the ratio of inner cylinder diameter D i to the gap width L constitutes the last parameter of the problem. Many fundamental aspects of patterns and their instabilities have been studied intensively over the past five decades in the context of Rayleigh-Bénard convection (RBC). In a traditional RBC experiment, a horizontal fluid layer of height dis confined between two thermally well conducting, parallel plates. When the difference T = T b −T t between the bottom-plate temperature T b and the top-plate temperature T t exceeds a critical value T crit , the conductive motionless state is unstable and convection sets in. The https://doi.org/10.1016/j.ijheatmasstransfer.2021.121151 0017-9310/© 2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) J.J. Serrano-Aguilera, F.J. Blanco-Rodríguez and L. Parras International Journal of Heat and Mass Transfer 172 (2021) 121151 Nomenclature A aspect ratio D diameter (m) JJacobian L gap width between cylinders (m) P rPrandtl number Ra Rayleigh number Nu Nusselt number T temperature (K) A differential operator B differential operator f generic function in the Newton-Raphson solver f m mass force (N ·kg −1 ) q solution column vector u velocity field (m ·s −1 ) R cylinder radius ratio dlayer height (m) ggravitational acceleration (m ·s −2 ) k thermal conductivity (W ·m −1 ·K −1 ) Nnumber of discretization nodes ppressure (Pa) rradial coordinate (m) ttime (s) u radial component of velocity field (m ·s −1 ) v azimuthal component of velocity field (m ·s −1 ) Abbreviations CVFEM control volume based finite element method FDM finite difference method FEM finite element method FVM finite volume method KG Kuehn & Goldstein LBM Lattice Boltzmann method MPS moving particle semi-implicit method RBC Rayleigh-Bénard convection S stable solution SPH smoothed particle hydrodynamics method SU dual stable-unstable solution U unstable solution UU dual unstable solution Greek Symbols αthermal diffusivity (m 2 ·s −1 ) βthermal volumetric expansion coefficient (K −1 ) T temperature difference (K) dimensionless distance λcomplex eigenvalue in global stability analysis λw wavelength (m) ∇differential operator nabla νkinematic viscosity (m 2 ·s −1 ) ωvorticity (s −1 ) ξChebyshev discretization variable (radial coordinate) temperature field with homogeneous boundary conditions ( K) stream function (m 2 ·s −1 ) ρdensity (kg ·m −3 ) θazimuthal coordinate Subscripts B Boussinesq BC boundary conditions b bottom c characteristic crit critical eq equivalent H hydrostatic i inner, imaginary part j iteration index m marginal max maximum o outer r radial s steady t top z z-axis (perpendicular to r −θdomain) θazimuthal Superscripts  perturbation ˆ dimensionless form average T transpose simplest pattern which can occur is that of straight, parallel convection rolls with a horizontal wavelength λw ≈2 d(wave number κ≈π/d). Such rolls can be found near onset; however, as the dimensionless distance ≡T / T crit −1 increases, the patterns often become progressively more complicated, and thus also more interesting. In their attempts to reach high Rayleigh numbers Ra, experimental researches have focused on convection layers with small aspect ratios (horizontal width to height) which do not permit the realizations of spontaneous large-scale patterns in highly turbulent convection. In fact, typical high Rayleigh number experiments are carried out with an aspect ratio of the order unity or less [4,5] . RBC is perhaps the most thoroughly investigated and understood pattern-forming system. The experimental setup is simple in principle and the basic physical mechanism (buoyancy vs dissipation), which can be described in terms of the Oberbeck- Boussinesq equations, is well understood. The basic flow field from the natural convection heat transfer in a horizontal cylindrical annulus for low value of Ra consists of two symmetric crescent-shaped eddies in which fluid rises near the inner hotter cylinder and sinks near the outer colder one. At high Ra, however, several kinds of convective flows which are dependent on P rand aspect ratio A can be developed. According to Kuehn and Goldstein [6] , the flow in a horizontal cylindrical annulus is steady over the range of Rayleigh number from 10 2 to 10 5 . Kuehn and Goldstein [7] also experimentally studied the flow patterns at Rayleigh numbers from 2 . 2 ·10 2 to 7 . 7 ·10 7 , and they found that the plume above the inner cylinder began to oscillate when Rayleigh number was near 2 ·10 5 and the entire plume was turbulent at Ra ≈2 ·10 6 . Their results also showed that the local heat transfer coefficient depended significantly on the eccentricity while the overall heat transfer coefficients changed by less than 10% with the change of eccentricity. Yoo [8] in his outstanding work, reported the occurrences of dual solutions at Ra larger than a critical value. Using a vorticity-streamfunction formulation, he observed dual steady solutions at Ra > Ra crit ≈3800 for wide gap annuli ( A = 2 ). Similar results were provided later by Mizushima and Hayashi [9] , Mizushima et al. [10] , Xin and Le Quéré [11] and Mercader et al. [12] . After the emergence of the nanofluids that have been shown to improve heat transfer process effectively, numerous research works have been also devoted to study heat transfer improvement in annuli systems using nanofluids [13,14] . In the last decades, a considerable amount of effort has been dedicated to study natural convection instabilities and bifurcations in horizontal concentric annuli. Powe et al. [15] experimentally investigated the bifurcation of natural convection of air ( P r = 0 . 7 ) by visualizing flow patterns, and categorized the flow patterns ob- 2 J.J. Serrano-Aguilera, F.J. Blanco-Rodríguez and L. Parras International Journal of Heat and Mass Transfer 172 (2021) 121151 Fig. 1. Sketch (left) and computational grid (right) of the problem. tained by their experiments and accumulated results by other researchers in a parameter space of (Ra, R ) where R ≡R o /R i = 1 + ( 2 /A ) is the cylinder radius ratio. This experiment was repeated and the classification was confirmed by Dyko et al. [16] . Busse [1] discussed the non-linear properties such as the dependence of the heat transport on Rayleigh and Prandtl numbers and the stability properties of thermal convection. Janssen and Henkes [17] analyzed the instabilities in three dimensional differentially-heated cavities with adiabatic horizontal walls. Petrone et al. [18] performed a stability analysis of numerical steady-state solutions, and provided a detail of the bifurcation diagram near the imperfect bifurcation for different radius ratio R = 1 . 2 , 1.4 and 2 at P r = 0 . 7 . Angeli et al. [19] provided a critical review of buoyancy-induced flow transitions in horizontal annuli. Soucasse et al. [20] studied the transitional regimes of natural convection in a differentially heated cubical cavity under the effects of wall and molecular gas radiation. Other well-known instability of the convective flows is the flickering instability [21,22] . In the present research work an extensive numerical analysis has been implemented in order to obtain a complete steady-state solution P r −Ra map for a fixed aspect ratio of A = 1 . 25 . Furthermore, a global stability analysis has been also performed and the marginal stability curve Ra crit = f(P r) has been obtained. Different regions with dual solutions have been identified and it is shown that in the unstable region the solution always reaches an attractor that appears as a vertical oscillating plume which is anchored to the inner cylinder and whose spatial and temporal structure can be explained through the flickering instability [21] , although its development is very limited by the upper surface. 2. Governing equations The problem of buoyancy (gravitational) induced motion of a fluid in a two-dimensional annulus is formulated using the Boussinesq approximation to solve the governing equations of mass, momentum and energy. Note that in the Boussinesq approximation, density variations are assumed to have a fixed part and another part that has a linear dependence on temperature: ρB (T ) = ρ[ 1 −β( T −T o ) ] , (1) where T is the temperature at a point within the fluid, βthe thermal volumetric expansion coefficient and T o is the temperature of the outer cylinder. Fig. 1 shows the geometry of the problem along with the polar coordinate system adopted in which the radial coordinate ris measured from the center of the system and θis measured clockwise from the upward vertical line. Moreover, the radial velocity u is positive radially outwards and the angular velocity v is positive in the clockwise direction. The governing equations for the natural convection in the region between horizontal concentric cylinders for an incompressible fluid can be written as ∇ ·u = 0 , (2) D u Dt = −1 ρ∇ p + ν∇ 2 u + ρB (T ) ρf m , (3) DT Dt = −α∇ 2 T . (4) where νis the kinematic viscosity, αis the thermal diffusivity, u = u ( r, θ) e r + v ( r, θ) e θ, D Dt ≡∂ ∂t + u ·∇is the material derivative and f m is the body force per unit mass which corresponds with gravity field for the case the annulus is not rotating. For the situation in which the fluid is quiescent ( u = 0 ) and at uniform temperature the above equations simplify as −∇p H + ρf m = 0 . (5) Eq. (5) reveals that we can redefine the pressure as p ∗= p −p H and the linear momentum equation can be written as D u Dt = −1 ρ∇ p ∗+ ν∇ 2 u + g β( T −T o ) e y , (6) where f m = −g e y . According to the Fig. 1 , the linear momentum equation can be written in polar coordinates as D u Dt = −1 ρ∇ p ∗+ ν∇ 2 u + g β( T −T o ) ( cos θe r −sin θe θ) , (7) The equations governing the phenomenon should also be deduced from the stream function-vorticity ( −ω) formulation. The velocity field is given as u = 1 r ∂ ∂θ e r −∂ ∂r e θ, (8) 3 J.J. Serrano-Aguilera, F.J. Blanco-Rodríguez and L. Parras International Journal of Heat and Mass Transfer 172 (2021) 121151 which automatically fulfils the incompressible constraint. The vorticity field for the resulting 2D problem has only one component perpendicular to the problem sketch as ω = ω e z which is given by ω = ( ∇ ×u ) ·e z = ∇ 2 , (9) Taking the curl of Eq. (7) and projecting onto the z-axis, the vorticity equation is obtained as Dω Dt = ν∇ 2 ω + g β1 r cos θ∂ ∂θ ( T −T o ) + sin θ∂ ∂r ( T −T o ) . (10) Furthermore, by setting the characteristic time and length of the problem as t c = L 2 αand l c = L = R o −R i , respectively D Dt ≡α L 2 D D ˆ t , ∇ 2 ≡1 L 2 ˆ ∇ 2 , ˆ T = T −T o T i −T o , r = L ˆ r , ω = α L 2 ˆ ω , = αˆ (11) the non-dimensional equations that govern the problem are: ˆ ω = ˆ ∇ 2 ˆ , (12) 1 P r D ˆ ω D ˆ t = ˆ ∇ 2 ˆ ω + Ra 1 ˆ r cos θ∂ ˆ T ∂θ + sin θ∂ ˆ T ∂ ˆ r , (13) D ˆ T D ˆ t = ˆ ∇ 2 ˆ T . (14) where P r = ν αand Ra = gβL 3 ( T i −T o ) να are the Prandtl and Rayleigh numbers, respectively. If we change the convention of sign of the vorticity ( ω → −ω) and use the free-fall velocity V f =  g L β( T i −T o ) instead of α/L as a characteristic velocity for scaling the velocity then the governing equations read ˆ ω = −ˆ ∇ 2 ˆ , (15) D ˆ ω D ˆ t =  P r Ra ˆ ∇ 2 ˆ ω −1 ˆ r cos θ∂ ˆ T ∂θ + sin θ∂ ˆ T ∂ ˆ r , (16) D ˆ T D ˆ t = 1 √ P r Ra ˆ ∇ 2 ˆ T , (17) which are the equations we have solved in the present work and coincide with those appearing in [23] . The bi-harmonic equation for the streamfunction ˆ is given by ∂ ˆ ∇ 2 ˆ  ∂ ˆ t + J ˆ ∇ 2 ˆ , ˆ =  P r Ra ˆ ∇ 4 ˆ + 1 ˆ r cos θ∂ ˆ T ∂θ + sin θ∂ ˆ T ∂ ˆ r , (18) ∂ ˆ T ∂ ˆ t + J( ˆ T , ˆ ) = 1 √ P r Ra ˆ ∇ 2 ˆ T . (19) where J(η, ˆ ) ≡1 ˆ r ∂ ˆ  ∂θ ∂η ∂ ˆ r −∂ ˆ  ∂ ˆ r ∂η ∂θ represents the convective term. The boundary conditions on the two impermeable isothermal walls are given by ˆ = ∂ ˆ  ∂ ˆ r = 0 , ˆ T = 1 , (20) on the inner cylinder ( ˆ r = R i /L ) and ˆ = ∂ ˆ  ∂ ˆ r = 0 , ˆ T = 0 , (21) on the outer cylinder ( ˆ r = R o /L ). In case symmetry with respect to the y axis is imposed, where only half of the annulus is taken as the computational domain, the following symmetric condition is applied along two vertical lines of symmetry at θ= 0 and θ= π: ˆ = ∂ ˆ  ∂θ = ∂ ˆ T ∂θ = 0 . (22) The solution of the heat equation with the same boundary conditions (which is only function of the radial coordinate) is an important reference to measure heat transfer efficiency of natural convection process. Therefore, in order to measure the convergence of the system and compare our solutions with those provided by other authors [6,24] , the local and average equivalent conductivity [23,25] are defined as follows: k eq ( ˆ r , θ) ≡Nu Nu c = −ˆ r ln (R ) ∂ ˆ T ∂n , ¯ k eq ( ˆ r ) ≡ k eq ( ˆ r , θ) ds 2 πˆ r = −ln (R ) 2 π ∂ ˆ T ∂n ds, (23) where Nu c ≡1 ln (R ) is the Nusselt number corresponding to the pure conductive flow and R = R o R i = 1 + 2 L D i = 1 + 2 A is the radius ratio. According to the geometric configuration we are dealing with ( ∂ ˆ T ∂n = ∂ ˆ T ∂r and ds = r dθ), the above equation can be further reduced to ¯ k eq,i = −ln (R ) 2 π( R −1 )  2π 0 ∂ ˆ T ∂ ˆ r dθ. (24) for the inner cylinder ( ˆ r = R i L ), and ¯ k eq,o = −R ln (R ) 2 π( R −1 )  2π 0 ∂ ˆ T ∂ ˆ r dθ. (25) for the outer cylinder ( ˆ r = R o L ). Fig. 2 shows the computational grid used for solving the problem where periodic boundary condition is employed in the azimuthal direction, θ. Additionally, Chebyshev discretization allows to radially concentrate the nodes near the walls ( ξ= ±1 ) in order to accurately solving both the thermal and velocity boundary layers. 3. Numerical method In recent years, a number of numerical methods have been applied to study natural convection in a horizontal cylindrical annulus, such as finite element method (FEM) [26,27] , finite difference method (FDM) [28] , finite volume method (FVM) [29,30] , lattice Boltzmann method (LBM) [31,32] , differential quadrature method [33] , control volume based finite element method (CVFEM) [34] , moving particle semi-implicit method (MPS) [35] and a smoothed particle hydrodynamics (SPH) method [36] . However, most of the previous studies did not consider the effect of the Prandtl number on the natural convection heat transfer in a horizontal cylindrical annulus. The objective of the present work is to extend our knowledge of flow patterns in a horizontal concentric annulus for a wide range of Rayleigh numbers ( 10 2 to 5 ·10 6 ) and Prandtl numbers (0.01 to 1). We perform a spectral numerical method for the spatial variables using Chebyshev and Fourier discretizations for the radial and azimuthal variables, respectively. Boundary conditions are directly implemented in the derivation matrices using the DM suite toolbox [37] . The dimensionless radial domain ( ˆ r ∈ [ R i /L, R o /L ] ) is mapped (see Fig. 2 ) in the Chebyshev domain ( ξ∈ [ −1 , 1] ) using the following path ˆ r = 1 2 ( ξ+ 1 ) + R i L = 1 2 ( ξ+ 1 ) + A 2 → ξ= 2 ˆ r −( A + 1 ) . (26) 4 J.J. Serrano-Aguilera, F.J. Blanco-Rodríguez and L. Parras International Journal of Heat and Mass Transfer 172 (2021) 121151 Fig. 2. Grid of the numerical domain. Since the derivation matrices implemented in the DM suite require Dirichlet boundary conditions for the temperature, we have decomposed the temperature field as ˆ T ( ˆ r , θ, ˆ t ) = ˆ ( ˆ r , θ, ˆ t ) + ˆ T BC ( ˆ r ) = ˆ ( ˆ r , θ, ˆ t ) + R o L −ˆ r  = ˆ ( ˆ r , θ, ˆ t ) + 1 + A 2 −ˆ r . (27) Substituting the above equation in the vorticity and temperature equations ( (15) –(17) ), we obtain ˆ ω = −ˆ ∇ 2 ˆ , (28) D ˆ ω D ˆ t =  P r Ra ˆ ∇ 2 ˆ ω −1 ˆ r cos θ∂ ˆ  ∂θ + sin θ∂ ˆ  ∂ ˆ r + sin θ, (29) D ˆ  D ˆ t −1 ˆ r ∂ ˆ  ∂θ = 1 √ P r Ra ˆ ∇ 2 ˆ −1 ˆ r . (30) The main advantage of the above formulation is that the deviatoric temperature field ˆ , has homogeneous Dirichlet boundary conditions at both the inner ( ˆ r = R i /L or ξ= −1 ) and outer ( R o /L or ξ= 1 ) walls [38] . To reduce the size of the system to solve, the flow variables used are ˆ and ˆ , and once defined the vector q = [ ˆ , ˆ ] T , the system of equations to solve can be written as A ∂q ∂ ˆ t = f ( ˆ t , q ) , (31) with A = ˆ ∇ 2 0 0 1 , (32) and f ( ˆ t , q ) = ⎛ ⎜ ⎜ ⎜ ⎝ −1 ˆ r ∂ ˆ  ∂θ ∂ ˆ ∇ 2 ˆ  ∂ ˆ r −∂ ˆ  ∂ ˆ r ∂ ˆ ∇ 2 ˆ  ∂θ +  Pr Ra ˆ ∇ 4 ˆ  +  1 ˆ r cos θ∂ ˆ  ∂θ + sin θ∂ ˆ  ∂ ˆ r −1  −1 ˆ r ∂ ˆ  ∂θ ∂ ˆ  ∂ ˆ r −∂ ˆ  ∂ ˆ r ∂ ˆ  ∂θ + 1 √ Pr Ra ˆ ∇ 2 ˆ −1 ˆ r + 1 ˆ r ∂ ˆ  ∂θ ⎞ ⎟ ⎟ ⎟ ⎠ . (33) 3.1. Steady solver To obtain the steady solution [ ˆ s , ˆ s ] T of Eq. (31) , a Newton- Raphson solver has been implemented as q j+1 = q j −∇ q f (q j ) −1 f j = q j −B (q j ) −1 f j , (34) where f j = f ( ˆ t , q j ) and B ≡∇ q f = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ −1 ˆ r ∂ ˆ j ∂θ ∂ ˆ ∇ 2 ∂ ˆ r + ∂ ˆ ∇ 2 ˆ j ∂ ˆ r ∂ ∂θ −∂ ˆ j ∂ ˆ r ∂ ˆ ∇ 2 ∂θ −∂ ˆ ∇ 2 ˆ j ∂θ ∂ ∂ ˆ r  +  Pr Ra ˆ ∇ 4 1 ˆ r cos θ∂ ∂θ + sin θ∂ ∂ ˆ r −1 ˆ r ∂ ˆ j ∂ ˆ r ∂ ∂θ −∂ ˆ j ∂θ ∂ ∂ ˆ r −∂ ∂θ  −1 ˆ r ∂ ˆ j ∂θ ∂ ∂ ˆ r −∂ ˆ j ∂ ˆ r ∂ ∂θ + 1 √ Pr Ra ˆ ∇ 2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , (35) and jmarks the index for each iteration of the algorithm. It is worth mentioning that the steady solutions have been sought taking advantage of their symmetry (anti-symmetric ˆ s and symmetric ˆ s ) with respect to the vertical axis ( θ= 0 ). It has been implemented by means of folding the differentiation matrices, which has allowed a significant reduction of the required computing resources by solving only within the θ∈ [0 , π] range. A steady solution is considered converged if max (| f j | ) < 10 −8 . 3.2. Unsteady solver To solve the unsteady problem a 4th order Runge-Kutta method has been implemented in the full domain (i.e. symmetry has not been imposed) as it is described in [39] . Defining the time step as ˆ t , the time evolution Eq. (31) are solved following the rule (where the superscript n defines the different time instants), q n +1 = q n + A −1 1 6 ˆ t ( k 1 + 2 k 2 + 2 k 3 + k 4 ) , (36) with k 1 = f ( ˆ t n , q n ) , (37) k 2 = f ( ˆ t n + ˆ t / 2 , q n + k 1 ˆ t / 2) , (38) k 3 = f ( ˆ t n + ˆ t / 2 , q n + k 2 ˆ t / 2) , (39) k 4 = f ( ˆ t n + ˆ t , q n + k 3 ˆ t ) . (40) 3.3. Global stability analysis Once the steady-state solutions are obtained, [ ˆ s , ˆ s ] T , we can perform a global stability analysis of the solution just defining a small variation of the solution around the steady-state solution, q = q s + q  . (41) 5 J.J. Serrano-Aguilera, F.J. Blanco-Rodríguez and L. Parras International Journal of Heat and Mass Transfer 172 (2021) 121151 Fig. 3. Comparison between the experimental local equivalent conductivity ¯ k eq by Kuehn and Goldstein (KG) [6] with the current numerical results (model) at Ra = 4 . 7 ·10 4 and Pr = 0 . 706 . If we define the perturbations q  = [ ˆ  , ˆ  ] T to be of the form, q  ( ˆ r , θ, ˆ t ) = ˜ q ( ˆ r , θ) e λˆ t , (42) the stability equations reduce to solving a generalized eigenvalue problem of the form, B ˜ q = λA ˜ q , (43) with A and B differential operators defined previously and λ= λr + i λi and ˜ q denoting the complex eigenvalue and eigenvector, respectively. Unlike in the steady solver, symmetry has not been forced in the stability analysis. It means that steady solution q s has been previously unfolded so that the stability analysis is applied over the full angular domain θ∈ [ 0 , 2 π] . 3.4. Validation For validation purposes the developed model has been solved using a mesh of N r ×N θ= 30 ×91 nodes (where 91 nodes span half of the angular domain). Its numerical solutions have been compared with the experimental results presented by Kuehn and Goldstein [6] . In Fig. 3 we can see the value of k eq at the inner and outer surfaces for Ra = 4 . 7 ·10 4 and P r = 0 . 706 obtained from the simulations made by the spectral code for different angular values ( θ), compared to the results obtained experimentally. It is noted that the numerical predictions are consistent with the experimental results measured on the surfaces of the concentric cylinders in [6] . Fig. 4 shows the results of the dimensionless temperature field as a function of the dimensionless radial coordinate inside the fluid, for different values of θcompared to the experimental results obtained for the same case in [6] . It presents an excellent agreement. Finally, Table 1 summarizes the results of k eq measured varying the values of Ra and P rfor A = 1 . 25 , and as it can be seen, values are obtained that vary in the third decimal place when they are compared to those of [24] . It is worth noting that due to conservation of energy any steady-state solution should verify that k eq,o = k eq,i and this is observed for our numerical results up to the fifth decimal place. 4. Numerical results This section details the procedure to obtain the complete map of steady-state solutions and shows the regions in the parameter Table 1 Comparison of the average equivalent conductivity ¯ k eq obtained by the current model with the experimental results provided by Kuehn & Goldstein [6] where inner and outer values are displayed separately. Current results are likewise compared with other numerical solutions provided in literature [24] . ¯ k eq Ra Pr A ≡D i /L ¯ k eq,i [6] ¯ k eq,o [6] ¯ k eq,io [24] ¯ k eq,io (current) 10 2 0.7 1.25 1.000 1.002 1.001 1.001 10 3 0.7 1.25 1.081 1.084 1.082 1.082 3 ·10 3 0.7 1.25 1.404 1.402 1.397 1.397 6 ·10 3 0.7 1.25 1.736 1.735 1.715 1.715 10 4 0.7 1.25 2.010 2.005 1.979 1.978 2 ·10 4 0.7 1.25 2.405 2.394 2.374 3 ·10 4 0.7 1.25 2.661 2.643 2.624 5 ·10 4 0.7 1.25 3.024 2.973 2.958 2.958 7 ·10 4 0.7 1.25 3.308 3.226 3.193 10 4 0.5 1.25 1.950 1.955 1.924 10 4 1.0 1.25 2.038 2.039 2.009 10 4 5.0 1.25 2.069 2.066 2.039 3 ·10 4 5.0 1.25 2.741 2.768 2.655 5 ·10 4 5.0 1.25 3.036 3.088 2.976 10 5 5.0 1.25 3.756 3.471 3.486 10 4 10.0 1.25 2.070 2.067 2.040 10 4 100.0 1.25 2.071 2.067 2.041 space where they are stable, which have been verified by a global stability study, as well as by solving the temporal evolution of the equations (28) –(30) . All the results presented in this research correspond to A ≡D i /L = 1 . 25 . 4.1. Map of steady state solutions The first problem to be addressed is the determination of steady-state solutions in a P r −Ra map. This is shown in Fig. 5 . To obtain this map, we have visited the different position of Ra ∈ [10 2 , 5 ·10 6 ] and P r ∈ [0 . 01 , 1] dividing the map in 75 and 149 points respectively according to a logarithmic scale. Therefore, this study aims to provide the steady state solution for 11175 cases in the range of Ra and P rused. Firstly, P rnumber is fixed and we increase the Ra number, solving the steady-state solution by means of Newton-Raphson iteration using as initial guess the result obtained for the previous Ra . After that sweep of the map, we repeat the process fixing Ra and increasing/decreasing P rnumbers from 10 −2 /1, respectively. This new sweep allows to obtain regions of dual steady-state solutions [40] . In Fig. 5 the different solutions are 6 J.J. Serrano-Aguilera, F.J. Blanco-Rodríguez and L. Parras International Journal of Heat and Mass Transfer 172 (2021) 121151 Fig. 4. Comparison between the experimental dimensionless radial temperature (circles) by Kuehn and Goldstein [6] with the current numerical model (continuous line) at Ra = 4 . 7 ·10 4 and Pr = 0 . 706 . Fig. 5. Map of steady solutions for A = 1 . 25 in the (Pr, Ra ) space. Some of the solutions analysed later are marked with name S for stable solutions and U for unstable solutions. marked, together with the special cases presented later. The solutions in the region hatched with lines sloping up to the right (blue lines) present a steady solution with one vertical plume emanating from the upper part of the inner cylinder. In the region hatched with lines sloping up to the left (red lines) the problem presents multiple plumes. The solutions used to discuss the structures created are marked by points and named using letters, S means that the solution is stable and U if unstable. There are regions where two steady-state solution coexist marked by two different letters, UU or SU depending whether both solutions are unstable or one stable and another unstable. Results of stability are discussed later. This map complements and fill the gaps of the work of [36] , where they presented only 24 different points in the Ra −P rdomain (see Fig. 10 ). Fig. 6 shows the stream function isocontours and isotherms of all the steady-state solutions where only a single solution has been found (single hatch pattern). The solutions for low Ra present a structure of two counter-rotating symmetrical vortices with respect to the y axis. For a given low Ra number ( Ra = 1 . 3397 ·10 2 ), increasing the Prandtl number we find the transition from solution S1 to S2, which basically increases the intensity of the rotation of the vortices and move their centers upwards. For P r = 0 . 8043 , increasing the Ra number, the solutions present the same kind of structure, but the intensity of the vortices increase and the center climbs up over the vertical height of the center of the annulus (see solutions S4 and S6). It can be also observed the appearance of an intense vertical plume where the temperature is maximum, placed on the y axis. For intermediate P rnumbers ( P r = 0 . 0733 , the solution passes from a stable solution in which the typical two counter-rotating vortices are deformed into four counter-rotating ones, to a solution where the flow is dominated by 6 counter rotating vortices (see transition from S3 to U2, U4). If the Rayleigh 7 J.J. Serrano-Aguilera, F.J. Blanco-Rodríguez and L. Parras International Journal of Heat and Mass Transfer 172 (2021) 121151 Fig. 6. Steady solutions corresponding to the region of the map ( Fig. 5 ) where only single steady solutions have been found. S and U stand for stable and unstable steady solutions respectively. Dimensionless temperature ˆ T have been plotted through colormap whereas isocontours describe the dimensionless stream function ˆ . 8 J.J. Serrano-Aguilera, F.J. Blanco-Rodríguez and L. Parras International Journal of Heat and Mass Transfer 172 (2021) 121151 number is further increased, the vortices at the upper part scale up in intensity and dominate (U5), creating a vertical plume that complicates further the structure of the flow in that region. It is worth noting that the temperature field near the inner surface shows an evolution from 5 plumes (U1, U2, U3) to 2 (U4 and U5). Finally, when Ra number is dominant (solution U6, for Ra = 3 . 2246 ·10 6 and P r = 0 . 1595 , the solution becomes a very intense plume flowing from the upper part of the inner cylinder. This kind of solution can be referred to as upwards vertical plume, since at the vertical axis ( y axis) the plume flows in the vertical direction. Regarding the two regions of dual steady-state solutions (regions with double hatch pattern), we show in Fig. 7 the two steady-state solutions found for the cases UU1, UU2 y UU3. The first case, UU1 corresponds to the region of dual solution for low values of P r( P r = 0 . 0128 ). Both solutions present basically the same structure with five plumes. Solutions UU2 to UU6 present two steady-state solutions, one with a very intense upward vertical plume, and the other one with 3 plumes, one vertical and the other one placed at θcrit ≈31 . 5 ◦. As we will discuss later, the solution with the upward vertical plume is the real attractor for the unsteady solution. Fig. 8 shows the flow structure in the dual solution region SU. There is a steady solution, that is the solution with one upwards vertical plume, and there also are another unstable solution presenting two plumes at angles ≈31 . 5 ◦as before. As we have shown in this section, the structure of the steady solutions obtained can be of a simple vertical plume or multiples plumes. The case of a single plume is easy to understand comparing with the results shown in [21] for an axisymmetric case, [41] for a 2D case of a variable density jet and [42] for a confined 2D jet/wake. Buoyancy forces are so intense that, due to the geometry of the problem, the boundary layer on the upper side of the inner cylinder starts to heat up the fluid creating a single planar jet. The force that generates the 2D jet is the buoyancy force that tries to initiate the movement of the fluid upwards, and once the movement has started, the flow becomes a planar variable density jet similarly to the one described in [41] . For Prandtl numbers of the order of unity, thermal and viscous diffusion becomes of the same order, and the described flow is created. However, for low P r numbers, the viscous diffusion is much lower than the thermal diffusion, allowing the appearance of more complicated structures. In these cases, two plumes are created at approximate angles of 31.5 ◦ on both sides of the symmetry plane. This allows to create a descending flow at the mid plane that moves against the buoyancy forces. The structure created can be described as two very interrelated counter-rotating vortices near the midplane. The stability characteristics of these more easily destabilized situations are analyzed in the next section. 4.2. Global stability Once the map of steady-state solutions to the problem of natural convection between two concentric cylinders has been obtained and described, it would be interesting to carry out a global stability analysis in order to elucidate the stability properties of each of the steady solutions found. To perform this task, the global stability equations described in Section 2 are used. Fig. 9 (a) shows the most unstable mode for the case of P r = 0 . 0733 as a function of the Rayleigh number. It is detected that the change of sign of the exponent λoccurs approximately for the critical value Ra crit ≈2438 . The eigenvector structure (real part) of the most unstable mode is also represented in Fig. 9 (b), (c) and (d) for the cases labelled S3, U2 and U4, respectively in the stability map. It can be seen that the structure is that of a sustained oscillation around the two 31.5 ◦ plumes. It is worth mentioning that while S3 and U2 solutions oscillate with a given frequency λmax,i  = 0 , the last solution U4 only has a non-null (positive) real part in the eigenvalue. Finally, repeating the process for all the P rnumber studied we have obtained the global stability map for the flow in this problem configuration ( A = 1 . 25 ). This is shown in Fig. 10 . The stability border (black continuous line) delimits the region of the P r −Ra map in which there is at least one stable steady solutions. As we will see below, this stable solution turns out to be the physical solution of the steady problem. It is worth mentioning that below the stability curve, the stable solutions are defined by a unique upwards vertical plume (hatch pattern with lines sloping up in blue color), but for low values of P r( P r ≤0 . 2 approximately) there is a small region of steady stable solutions with multiple plumes (hatch pattern with lines sloping down in red color). This numerical stability curve has been tested by means of unsteady numerical simulations in different points, obtaining the same values of Ra crit for the transition, where the whole angular domain θ∈ [ 0 , 2 π] has been regarded (no symmetry imposed). As it will be shown in the following section by means of transient simulations, whenever there is a solution with a vertical plume structure, it will act as a flow attractor. Therefore, the stability in the area of P r  0 . 25 is dominated by the stability of that solution and that is what is finally represented in Fig. 10 . Fig. 10 also displays some cases of the stability test by Yang et al. [36] accomplished by means of transient analysis. According to them, four states are regarded (see legend of Fig. 10 ): UP1 (unstable with one plume), UPM (unstable with multiple plumes), SP1 (stable with on plume) and SPM (stable with multiple plumes). It can be observed that the region with stable solutions resembles with the stable region found in the current research work. The only significant discrepancy is based on the observations made in the case where P r = 0 . 1 and Ra = 10 5 . Unlike the conclusion drawn by the aforementioned authors, as far as we have found, it seems to behave as an unstable solution, where the physical solution is a sustained oscillation around the steady state solution. As we discussed previously, the typical flow structures may contain a simple or multiple plumes. In the case of the simple plume, the buoyancy force creates the planar variable density jet. Bharawaj and Das [41] discussed that the only global instability that appear in this flow is a ‘puffing’ plume (varicose instability) but the one presented here is of sinuous type (see supplementary material). This instability is created due to the competition between the inertia and viscous forces in the shear layers of the planar jet. The flow can be seen as an oscillation over the steady state solution (equilibrium solution) as it can be observed in Fig. 10 by Tammisola et al. [42] . The case of multiple plumes have a more complex equilibrium. The structure that appears on the upper part of the cylinder is formed by two ascending plumes placed at approximately 31.5 ◦on both sides of the symmetry plane and a descending flow at the mid plane creating two counter-rotating vortices. As we have discussed previously, these cases are mainly controlled by thermal diffusion, so when they are in equilibrium, the heat transfer created on both cylinders and the convection due to the counter-rotating vortices equilibrate. If Ra is increased over a threshold, this heat transfer equilibrium is lost so the entire flow is destabilized, moving the solution to a single-point attractor as it is explained in the following subsection. The cause of the instability is created outside of the ascending jets, as it can be observed in figures 9 (b)-(c) where the structure of the most unstable mode for the same low Prandtl number ( P r = 0 . 0733 ) and three different Rayleigh numbers are plotted. The eigenvector corresponding to the most unstable mode has their maxima at approximate angles of 45 ◦in the lower region of the ascending jets so those parts will be the first destabilized to a new solution instead of oscillating around the steady state. 9 J.J. Serrano-Aguilera, F.J. Blanco-Rodríguez and L. Parras International Journal of Heat and Mass Transfer 172 (2021) 121151 of wall and molecular gas radiation, Phys. Fluids 26 (2) (2014) 024105, doi: 10. 1063/1.4864265 . [21] J.M. Lopez, F. Marques, Instability of plumes driven by localized heating, J. Fluid Mech. 736 (2013) 616–640, doi: 10.1017/jfm.2013.537 . [22] W. Coenen, L. Lesshafft, X. Garnaud, A. Sevilla, Global instability of low-density jets, J. Fluid Mech. 820 (2017) 187–207, doi: 10.1017/jfm.2017.203 . [23] D. Ho-Minh , D. Mai , T. Cong Cong , A Galerkin-RBF approach for the streamfunction-vorticity-temperature formulation of natural convection in 2D enclosured domains, CMES 44 (2009) . [24] C. Shu, Application of differential quadrature method to simulate natural convection in a concentric annulus, Int. J. Numer. MethodsFluids 30 (8) (1999) 977–993, doi: 10.1002/(SICI)1097-0363(19990830)30:8 < 977:: AIDFLD873 > 3.0.CO;2- J . [25] J.D. Chung, C.-J. Kim, H. Yoo, J.S. Lee, Numerical investigation on the bifurcative natural convection in a horizontal concentric annulus, Numer. Heat Transf. Part A 36 (3) (1999) 291–307, doi: 10.1080/104077899274778 . [26] F. Selimefendigil, H.F. Öztop, Conjugate natural convection in a nanofluid filled partitioned horizontal annulus formed by two isothermal cylinder surfaces under magnetic field, Int. J.Heat Mass Transf. 108 (2017) 156–171, doi: 10.1016/j. ijheatmasstransfer.2016.11.080 . [27] K. Kahveci , Stability of unsteady mixed convection in a horizontal concentric annulus, J. Appl. Fluid Mech. 9 (5) (2016) 2141–2147 . [28] J. Wu, P. Traoré, M. Zhang, A.T. Pérez, P.A. Váquez, Charge injection enhanced natural convection heat transfer in horizontal concentric annuli filled with a dielectric liquid, Int. J.Heat Mass Transf. 92 (2016) 139–148, doi: 10.1016/j. ijheatmasstransfer.2015.08.088 . [29] Y. Zhao , B. Zhang , A high-order characteristics upwind FV method for incompressible flow and heat transfer simulation on unstructured grids, Comput. Methods in Appl. Mech.Eng. 190 (5–7) (20 0 0) 733–756 . [30] E. Abu-Nada, Effects of variable viscosity and thermal conductivity of Al 2 O 3 - water nanofluid on heat transfer enhancement in natural convection, Int. J. Heat Fluid Flow 30 (4) (2009) 679–690, doi: 10.1016/j.ijheatfluidflow.2009.02. 003 . [31] K. Luo, H.-L. Yi, H.-P. Tan, Eccentricity effect on bifurcation and dual solutions in transient natural convection in a horizontal annulus, Int. J. Therm. Sci. 89 (2015) 283–293, doi: 10.1016/j.ijthermalsci.2014.11.020 . [32] Y. Cao, Y. Zhang, Investigation on the natural convection in horizontal concentric annulus using the variable property-based lattice Boltzmann flux solver, Int. J.Heat Mass Transf. 111 (2017) 1260–1271, doi: 10.1016/j.ijheatmasstransfer. 2017.04.071 . [33] C. Shu, H. Xue, Y. Zhu, Numerical study of natural convection in an eccentric annulus between a square outer cylinder and a circular inner cylinder using DQ method, Int. J.Heat Mass Transf. 44 (17) (2001) 3321–3333, doi: 10.1016/ S0 017-9310(0 0)0 0357-4 . [34] M. Sheikholeslami, M. Gorji-Bandpy, D. Ganji, S. Soleimani, Effect of a magnetic field on natural convection in an inclined half-annulus enclosure filled with Cu-water nanofluid using CVFEM, Adv. Powder Technol. 24 (6) (2013) 980–991, doi: 10.1016/j.apt.2013.01.012 . [35] Y. Liang, Z. Sun, G. Xi, L. Liu, Numerical models for heat conduction and natural convection with symmetry boundary condition based on particle method, Int. J.Heat Mass Transf. 88 (2015) 433–4 4 4, doi: 10.1016/j.ijheatmasstransfer.2015. 04.105 . [36] X. Yang, S.-C. Kong, Numerical study of natural convection in a horizontal concentric annulus using smoothed particle hydrodynamics, Eng. Anal. Boundary Ele. 102 (2019) 11–20, doi: 10.1016/j.enganabound.2019.02.007 . [37] J.A. Weideman, S.C. Reddy, A matlab differentiation matrix suite, ACM Trans. Math. Softw. 26 (4) (20 0 0) 465–519, doi: 10.1145/365723.365727 . [38] M. Napolitano, G. Pascazio, L. Quartapelle, A review of vorticity conditions in the numerical solution of the ζ−ψequations, Comput. Fluids 28 (2) (1999) 139–185, doi: 10.1016/S0 045-7930(98)0 0 024-3 . [39] U.M. Ascher , L.R. Petzold , Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations, first ed., Society for Industrial and Applied Mathematics, USA, 1998 . [40] J.-S. Yoo, Dual steady solutions in natural convection between horizontal concentric cylinders, Int. J. Heat Fluid Flow 17 (6) (1996) 587–593, doi: 10.1016/ S0142-727X(96)0 0 064-1 . [41] K.K. Bharadwaj, D. Das, Puffing in planar buoyant plumes: BiGlobal instability analysis and experiments, J. Fluid Mech. 863 (2019) 817–849, doi: 10.1017/jfm. 2018.1022 . [42] O. Tammisola, F. Lundell, P. Schlatter, A. Wehrfritz, L.D. Södeberg, Global linear and nonlinear stability of viscous confined plane wakes with co-flow, J. Fluid Mech. 675 (2011) 397–434, doi: 10.1017/jfm.2011.24 . 16