NUMERICAL ANALYSIS OF NATURAL CONVECTION IN INTERNALLY FINNED HORIZONTAL ANNULI

Detailed numerical analysis is presented for natural convection heat transfer in internally finned horizontal annuli. Governing equations are discretized using the finite volume method, and solved using SIMPLE algorithm with Quick scheme. The results show that the flow and heat transfer can reach steady state when the Rayleigh number is below 2×104. When the Rayleigh number is greater than 3×104, two different types of numerical solutions under the same parameters can be obtained for different initial conditions. The critical Rayleigh numbers with two different initial conditions are different from steady to unsteady solutions. The oscillatory flow undergoes several bifurcations and ultimately evolves to a chaotic flow as the Rayleigh numbers increase.


INTRODUCTION
Natural convection in concentric annuli has been studied extensively due to both their theoretical interests and engineering applications such as the storage system of thermal energy, cooling system of electronic and mechanical parts, cooling of electric transmission cable, etc. Researchers (Kuhen and Goldstein, 1978;Kumar, 1988;Cao and Zhang, 2017;Cao, 2017;Zhang and Cao, 2018) studied the natural convection in horizontal annuli, where the inner cylinder is heated by the application of a constant heat flux and the outer cylinder is isothermally cooled. The Rayleigh number and Prandtl number are constants in the studies (Cao and Zhang, 2017;Cao, 2017). Cao and Zhang (2017) analyzed the necessity of considering the total variation in fluid properties in terms of the convective heat transfer, Cao (2017) stressed the essential role of the temperature difference ratio. The flow instability behaviors of the thermal convection were discussed with Rayleigh number varying from 10 4 to 10 6 (Zhang and Cao, 2018). Other researchers (Mahmoudinezhad et al., 2018;Tsui et al., 1984;Yang et al., 1988) reported a number of experimental and numerical studies for a variety of configurations, boundary conditions and fluid properties. The flow and thermal fields have been well studied for the natural convection induced by internal heating in horizontal concentric annuli.
Many helpful results from the research of natural convection indicated that the heat transfer between horizontal annuli is limited by the surface area of the inner cylinder. As a result, fins attached to the inner cylinder are sometimes employed to increase the surface area of heat transfer, leading to an increase in the heat transfer between the cylinders. Chai and Patankar (1993) numerically investigated the effect of radial fins on laminar natural convection in horizontal annuli. Their results indicate that the fin orientation shows no significant effect on average Nusselt number, and the average Nusselt number increases with increasing Rayleigh number and decreases with increasing fin height. Farinas et al. (1997) studied numerically the effect of internal fins on the * Corresponding author. Email: zhangkun52015@163.com. flow pattern, temperature distribution and heat transfer between concentric horizontal cylinders. The results indicated that the heat transfer is about the same for three different fin configurations, but the best fin efficiency is associated with the fin with a round tip. Sheikhzadeh et al. (2013) investigated the effects of the Rayleigh numbers, number of fins, and length of the fins on the flow and heat transfer in the internally finned horizontal annuli. Ha and Kim (2004) studied the fin configuration and ratio of annulus gap width to the inner cylinder radius on the flow and heat transfer. Kiwan and Zeitoun (2008) studied the effect of fin conductivity ratio, Darcy number, and Rayleigh number on the average Nusselt number for fins made of porous material when attached to the inner cylinder of the annulus between two concentric cylinders. Rahnama et al. (1999) studied the effect of fin orientation for Rayleigh numbers less than 10 6 and the optimum length for the maximum heat transfer rate. The existing numerical investigation of natural convection in internally finned horizontal annuli has focused on unique and steady flow fields under the same boundary conditions. The recently studies about the horizontal annular with fin (Ryad et al., 2019;Nemati et al., 2020;Hadidi et al., 2020). The practical engineering applications such as Thermal Energy Storage Units (Ali Rabienataj Darzi et al. (2016); Hu et al., 2018).
However, the flow fields in natural convection under the same boundary conditions are always not unique and most of flows are unsteady in the world. The previous studies indicated that a laminar steady flow and heat transfer is observed for lower Rayleigh number, while unsteady self-sustained phenomenon and multiply solutions are obtained at high Rayleigh number (Petrone et al., 2006;Powe et al., 1969Powe et al., , 1971. Many numerical simulations have been carried out to study the natural convection instabilities and bifurcations in horizontal concentric annuli. Powe et al. (1969Powe et al. ( , 1971 visualized the flow patterns with smoke as a tracer and obtained the Rayleigh number at which the flow will change from a steady flow to an unsteady flow. Rao et al. (1985) investigated the transient oscillatory phenomena and the numerical determination of the critical Rayleigh number at which unsteady flow occurs. Cheddadi et al. (1992) studied the natural convection bifurcation in a horizontal annulus, and showed that flow pattern was not unique and depended on the initial conditions at high Rayleigh number. Yoo (1999Yoo ( , 2003 investigated the occurrences of dual solutions and the bifurcation sequences to the chaos for the natural convection in horizontal concentric annuli. Similar results were provided by Chung et al. (1999), and Desrayaud et al. (2000). Mizushima et al. (2001) investigated a rigorous determination of the behavior of the system near the bifurcation point, confirming that the dual solutions were due to the static bifurcation. Hu et al. (2015) studied three types of steady-state solutions and the effects of the Rayleigh number on three kinds of flow patterns.
In addition, the non-linear characteristics of natural convection in a more complex domain were discussed by many researchers. Lou et al. (2015) investigated the bifurcation phenomena and the existence of dual solutions in natural convection in eccentric annulus. Yu et al. (2010) studied the transient natural convective heat transfer of a low Prandtl number fluid inside a horizontal circular cylinder with an inner coaxial triangular cylinder. Their results indicated that temporal phases of the flow development are identified as initializing developing, transitioning, and steady/quasi-steady state or oscillating. Zhang et al. (2011Zhang et al. ( , 2014 numerically investigated nonlinear characteristics of natural convection heat transfer in a cylindrical envelope with an internal concentric cylinder with slots. Their results indicated that the oscillatory flow undergoes several bifurcations and ultimately evolves to a chaotic flow. Zhao et al. (2018) investigated the characteristics of evolution from laminar to chaotic nature convection inside the horizontal annulus with an internal slotted circle, they found that after the first bifurcation, the oscillatory flow underwent several bifurcations and finally transited into a chaotic flow. Park (2018, 2019) investigated bifurcation phenomena of natural convection in horizontal annulus under selfinduced circular magnetic fields. They found that magnetic field can boost the flow instability and induce flow bifurcation even in the Rayleigh number region where the bifurcation does not appear.
The objective of this paper is to numerically study natural convection heat transfer in internally finned horizontal annuli. The governing equations will be discretized using the SIMPLE algorithm with Quick scheme. The effect of initial conditions on the final flow structure and the multiplicity of solutions related to static bifurcation will be investigated. The space trajectory of velocity at the sample point, the temporal evolution and power spectrum of the average dimensionless equivalent thermal conductivity (Keq) will be studied to the effect of the increasing Rayleigh numbers on the routes to chaos.

PROBLEM FORMULTION
A schematic diagram of the geometries analyzed in the present study is shown in Fig. 1. The inner and outer cylinders, of radii ri and ro, are kept at uniform but different temperatures Ti and To, respectively. The fins of height of Lf and thickness of d are placed on the inner cylinder. As a result of the temperature difference between the two cylinders, density gradient occurs which leads to natural convection.
The problem is considered to be two-dimensional and unsteady convection in annular space between long horizontal concentric cylinders. The fluid properties are also assumed to be constant, except for the density in the buoyancy term, which follows the Boussinesq approximation. The non-dimensional governing equation for the problem can be written as where the source terms are The dimensionless parameters in Eqs.
(1)-(6) are defined as follows The Rayleigh number Ra, Prandtl number Pr, reference velocity UR , and relative thermal conductivity K, are defined as where β, α , ν are thermal expansion coefficient, thermal diffusivity, kinematic viscosity of the fluid, respectively. K is the relative thermal conductivity, which takes the value of one in the fluid region and ks/kf in the solid region. L is the gap width of annulus, and g is the gravitational acceleration. The boundary conditions for Eqs. (1)-(4) are as following: The following periodical boundary conditions are given at θ = 0 and θ = 2π To observe the total heat transfer effect, the average dimensionless equivalent thermal conductivity based on the whole outer circle is defined as where Q is the total heat transfer rate.

NUMERICAL PROCEDURES
The governing equations were discretized on a staggered mesh, and the central difference scheme was adopted for the discretization of the diffusion terms. The SIMPLE algorithm with Quick scheme was used for handling the pressure-velocity coupling. To verify the numerical code used in the current investigation, the results of the present code for clear annulus are tested and compared with the experimental results obtained by Kuehn and Goldstein (1978). It is noted that Ram is a modified Rayleigh number introduced by Alshahrani and Zeitoun: The average dimensionless equivalent conductivity Keq of present code are in good agreement with other experimental and numerical results in the range of Ram = 1 to 4, the maximum error is less than 6%. Thus, the present computational method can be used for investigating the flow and heat transfer characteristics in the internally finned annuli.

Steady-state solutions
When the parameters are Lf/L = 0.5, Nf = 6, d/L = 0.05, Pr = 0.701, Ra = 10 5 , the grid independence of numerical results is studied. Five different grid sizes in Table 1 are used, which are 80×50, 120×60, 160×70, 200×80, and 240×90, respectively. Numerical experiments showed that the relative error of Keq at Fo = 10 corresponding to grid number of 200×80 and 240×90 is less than 1%, 200×80 grid points are adequate to yield accurate results. Table 2 shows that the Keq at Fo = 10 of different time steps for 200×80 grid points. As is observed from this table, the difference between Keq obtained using the time steps of 0.01 and 0.001 is negligible. Thus, in the following simulations, the time step of 0.01 is adopted. The numerical simulation revealed that the steady and symmetric solution at lower Rayleigh numbers can be reached from different initial conditions. The isotherms and velocity vector fields at Ra = 10 4 obtained by current and previous methods (Chai and Patankar, 1993) were shown in Fig. 2. Both the two numerical solutions have the same results for these cases with low Rayleigh numbers. The fluid is heated by the internally finned cylinder with higher temperature. The fluid moves upwards and downwards along the surfaces of inner and outer cylinders, respectively. When the fluid meets the fins located at the center of the annuli, it is forced to bend and moves along the fins. There are six vortices formed in the enclosure of annuli. Two in the right side of the enclosure rotate in the clockwise and another two in the left side rotate in the counterclockwise direction. At the upper part of the enclosure, the fluid moves upwards along the fins attached to the inner cylinder. When the ascending fluid meets the cold outer cylinder, it changes its direction and moves downwards along the outer cold wall. The two vortices in this domain rotate in the opposite directions. At the bottom part of the enclosure, the outer cylinder at lower temperature is located at the lower part relatively in the gravitational direction than the inner cylinder at higher temperature. Therefore, the natural convection can be limited by the surface of outer cylinder, and the fluid between the inner and outer cylinders is in a stable condition and is almost stagnant.

Static branching solution
The effects of the initial conditions on the final solution are calculated and analyzed numerically in this section. When the parameters are Ro = 2, Ri = 1, Lf / L = 0.5, Nf = 6, d = 0.05, Pr = 0.701, two different types of flow patterns at Ra = 5×10 4 can be reached from different initial fields even though other conditions are the same, as is shown in Fig.3. The reason for this phenomenon is the nonlinear characteristics of this system, and the initial conditions play a determining role in problems involving instability. The final bifurcation results are regarded to be sensitive to the initial conditions, and the initial conditions are expected to determine which of the branches the solutions will finally be led to. When the initial condition is given zero velocity and uniform field, the five-roll flow pattern is obtained and a big roll cell at the top of the enclosure rotates in the counterclockwise direction in Fig. 3(a). Therefore, the asymmetric temperature field can be reached with axis-symmetric boundary condition. When the initial condition is given six-roll flow field in Fig. 2, the number of rolls in flow field remains unchanged and two roll cells rotate in the clockwise and counterclockwise directions at the top of the enclosure in Fig. 3(b). For the second type of six roll flow pattern, the streamlines and temperature field are in good agreement with the results obtained by Chai (1993). However, the first type of five roll flow pattern in this paper can't be obtained by the numerical method in Chai (1993). Obviously, the numerical method in Chai (1993) that the right halves of the annuli were chosen as the whole solution domains is not fit for some certain flow and heat transfer, which corresponding to the non-linear phenomenon. The maximum difference of the two types of flow pattern occurs at the top part of enclosure. The outer cylinder at a constant low temperature is located at the upper part relatively in the gravitational direction than the inner cylinder at a constant high temperature. The fluid between the two fins at the top of enclosure is in an unstable condition and shows the pattern similar to the flow and heat transfer for the case of Rayleigh Benard convection in the horizontal rectangular enclosure. From these results, the authors deduce that the third type of flow pattern should be obtained if the initial condition is feasible to be selected. In the third flow pattern, the number of roll cells in the enclosure is five and the big vortex at the top of the enclosure rotates in the clockwise direction, which is opposite with the direction of rotation in the second type. The numerical results indicate that the values of average dimensionless equivalent conductivity Keq for six-roll and five-roll flow fields are 3.68 and 3.24, respectively. The relative difference in Keq is higher than 10%. This is because the heat transfer rate depends on the flow pattern besides flow velocity and surface area. In the six-roll flow pattern, two opposite rolls occur above the internal cylinder and the isotherm curved more. obviously.The local heat transfer rate is high at the region of impingement.  Fig. 4 shows that the value of Keq for different Rayleigh numbers. For the case of Ra lower than 2×10 4 , the only solution of flow and heat transfer can be reached even from the different initial conditions. As the Rayleigh number increases, a static bifurcation occurred at a critical Rayleigh number between 2×10 4 and 3×10 4 . The multiplicity of solutions under the same parameters occurs in a narrow band of Rayleigh numbers from 3×10 4 to 10 5 . In this scope, the unstable change of flow patterns is usually caused by small initial perturbations. The heat transfer rate in the six-roll flow patterns is enhanced more obviously than that in the fiveroll flow pattern

Self-sustained oscillation
Numerical simulations of the unsteady natural convection are performed at different Rayleigh numbers. When time infinities are sufficiently long, the oscillatory phenomenon of natural convection from the zero initial velocity and uniform temperature fields can be observed for Ra = 10 6 . Attention is especially paid to the non-periodic oscillated solutions in the duration of dimensionless time from Fo = 2000 to Fo = 2080. The first type of the oscillatory flow field reached from zero-roll pattern initial field is shown in Fig. 5. It can be seen that the asymmetric oscillation phenomenon exists above the inner cylinder, although the boundary conditions are symmetric and stable. There are two vortices between the two fins at the top of enclosure. The top vortex on the right is significantly larger than the left vortex in Fig.5 (a). The top vortex on the left moves up and shrinks with time gradually. Almost at the same time, the top vortex on the right becomes larger and occupies all of the space between the two fins, as is shown in Fig. 5 (b). Then the big vortex becomes elliptic and moves up, providing space for a third vortex to develop in Fig.5 (c). The third vortex grows larger and another small vortex at the top disappears in Fig.5 (d). The vortex tends to moves backwards and forwards between the two fins, as is shown in Fig.5 (d-h). Finally, the vortex moves up along the wall of the left fin and the two vortices between the two fins occur on the left and right at the upper part of the enclosure, respectively. The temperature field is also altered obviously, as depicted by Fig. 6. It should be noted that the change of heat transfer rate are closely related to these flow pattern. The largest temperature gradient and local heat transfer rate occur at the right top of outer cylindrical envelope. The isotherm lines spacing on the right side in the annulus space are smaller than the lines on the left side, so the overall heat transfer rate on the right side is larger than that on the left side. The maximum of heat transfer rate appear at dimensionless time Fo = 2000. The values of average dimensionless equivalent conductivity Keq at Fo = 2010, 2020 and 2080 are greater than those values for other cases. The two vortices at the upper part of the enclosure occur along the surface of the outer cylinder at Fo = 2010, 2020 and 2080, as are shown in Fig.5 (a-c) and Fig. 5 (h), those flow patterns are more favorable for the enhancement of heat transfer. When the initial condition is given a six-roll pattern, the second type of oscillatory flow field in the duration of the dimensionless time from 2000 to 2080 can also be obtained for Ra = 10 6 , as are shown in Fig. 7. There are four vortices in the annuli: two between the fins at the top and two at each side of inner cylinder. The two vortices between the two fins at the upper part of the enclosure have nearly the same sizes. Correspondingly, the temperature field at different time is shown in Fig.  8. The isotherms are closely packed at the top and the local heat transfer in this region is highest while that at the bottom is lowest. Compared with the above oscillatory solution, the oscillation in the second type of flow field is inconspicuous and the number of the big vortices in the annuli remains unchanged. The two big vortices between the two fins at the top of inner cylinder are favorable for increasing the heat transfer rate at the top surface of outer cylinder. Therefore, the heat transfer rate can be enhanced by the second type of flow pattern more obviously.
The time signal of average thermal conductivity Keq for different initial conditions at Ra = 10 6 is shown in Fig. 9. When time infinities are sufficiently long, the mean values of Keq with the two different initial conditions for zero-roll and six-roll flow pattern are 6.87 and 7.88, respectively. Obviously, the second type of oscillatory flow field is more favorable for the enhancement of heat transfer. The maximum value of the amplitude of Keq is lower than 0.23 for the case of the initial condition with six-roll flow pattern and it is equal to 1.8 for the other case. The oscillatory phenomenon is more obvious for the first type of flow field. More numerical simulation indicated that there exists a critical Rayleigh number Rac, which can dramatically change the flow structure. For the case of Ra<Rac, the natural convection is steady. The flow turns to oscillated solution after the first Hopf bifurcation for Ra>Rac. The initial conditions play a determining role in problems involving instability and multiplicity of solutions. With the increasing Rayleigh numbers, the flow and heat transfer with different initial conditions develop in the different ways. When the zero velocity and uniform temperature fields are selected as the initial conditions, the flow and heat transfer transit from steady to unsteady at Rac between 10 5 and 2×10 5 . When the six-roll flow fields obtained for the previously lower Rayleigh number are selected as the initial conditions, the steady state solutions can still persist for the Rayleigh number beyond 2×10 5 up to 6×10 5 . For the case of Ra>7×10 5 , the non-periodic oscillated solution can be reached from the initial fields with the six roll flow pattern. Therefore, the critical Rayleigh numbers with two different initial conditions are different from steady to unsteady solutions. The first type of flow pattern reached from zero-roll flow pattern tends to lose stability more easily because of the lower value of the critical Rayleigh number.

On the route to chaos
There are several ways to recognize the route to chaos: analyzing their time history, phase space trajectories, power spectral method and so on. The effects of the increasing Rayleigh numbers on the route to chaos were analyzed by using phase space of velocity at the sample point (27,27), as is shown in Fig. 10. The sample point is selected in order to analyze the nonlinear characteristics of the complicated system further.
There is no special reason to select a particular point since if the system is periodic it will be so at every point in the whole computational domain. However, some points near boundary can't be selected due to small amplitude. A power spectrum of Hamming window was obtained to identify the route from steady to chaotic. The time history and power spectrum analysis of average thermal conductivity Keq are shown in Fig.  11 and Fig. 12, respectively. The data of dimensionless time less than 1000 is not shown in those figures for the convenience of presentation. The results indicated that the time signal of average thermal conductivity Keq is a constant in Fig. 11(a), and the phase trajectory of velocity in Fig.10. (a) is a limit point for Ra = 10 4 . The limit point loses its stability and becomes a steady cycle after the first Hopf bifurcation, as is shown  in Fig. 10(b). It can be seen from Fig. 11(a) and Fig. 11(b) that the periodic solutions of average thermal conductivity Keq appear at Ra = 2×10 5 , and correspondingly only one harmonic peak occurs in Fig. 12(b). There is one more bifurcation at successive critical values of Rayleigh numbers from 2×10 5 to 10 7 . For a truly chaotic system, there should be a continuous broadband spectrum in Fig. 12(c). It is ascertained in Fig.  11(c) and Fig. 12(c) that periodicity is lost at Ra = 10 7 . The results show that the flow undergoes several bifurcations and ultimately evolves to a chaotic flow after the first bifurcation.

CONCLUSIONS
Numerical analyses of natural convection in internally finned horizontal annuli are carried out. Nonlinear characteristics of natural convection with increasing Rayleigh numbers are discussed and analyzed. Flow and heat transfer characteristics are investigated for a Rayleigh number range of 10 4 to 10 7 , while the Pr number is taken to be 0.71.
The results indicated that the static bifurcation and self-sustained oscillation of the natural convection can be obtained at different Rayleigh numbers. When the Rayleigh numbers are small enough, the flow and heat transfer is unique and steady state. As the increasing Rayleigh numbers, two different types of flow with six or five-roll pattern occurs in a narrow band of Rayleigh numbers from 3×10 4 to 10 5 . The initial conditions play a determining role in problems involving instability and determine which of the branches the solutions will finally be led to. When the Rayleigh number is greater than a critical value, a time-dependent behavior is observed. The five-roll flow pattern tends to lose stability more easily because of the lower value of the critical Rayleigh number. The flow undergoes several bifurcations and ultimately evolves to a chaotic flow after the first bifurcation. There is one more bifurcation at successive critical values of Rayleigh numbers from 2×10 5 to 10 7 . The natural convection undergoes many bifurcations and transits from the steady state to the determinist chaos with the increasing Rayleigh numbers.