Study on Bifurcation and Dual Solutions in Natural Convection in a Horizontal Annulus with Rotating Inner Cylinder Using Thermal Immersed Boundary-Lattice Boltzmann Method

A numerical investigation has been carried out to understand the mechanism of the rotation effect on bifurcation and dual solutions in natural convection within a horizontal annulus. A thermal immersed boundary-lattice Boltzmann method was used to resolve the annular flow domain covered by a Cartesian mesh. The Rayleigh number based on the gap width is fixed at 104. The rotation effect on the natural convection is analyzed by streamlines, isotherms, phase portrait and bifurcation diagram. Our results manifest the existence of three convection patterns in a horizontal annulus with rotating inner cylinder which affect the heat transfer in different ways, and the linear speed (Ui*) determines the proportion of each convection. Comparison of average Nusselt number versus linear speed for the inner cylinder indicates the existence of the three different mechanisms which drive the convection in a rotation system. The convection pattern caused by rotation reduces the heat transfer efficiency. Our results in phase portraits also reveal the differences among different convection patterns.


Introduction
Natural convection in a horizontal annulus is a useful subject due to its importance in theoretical study and in many engineering applications, such as thermal energy storage systems, cooling systems in electronic components, brake assembly and bearing fittings. Considering its simple geometry, this system has become a classical model and has been well studied by many researchers. Early in 1966, Grigull and Hauf [1] investigated natural convection in a horizontal annulus experimentally. In 1969, Powe et al. [2] designed an experiment and delineated the convection with different types of flow depending on the Grashof number and radius ratio. Kuehn and Goldstein [3] studied the natural convection in concentric cylinders numerically and experimentally. Nguyen et al. [4] presented the flow patterns and heat transfer rates in terms of the radius ratio, the Rayleigh number and inversion parameter. Dyko et al. [5] presented a three-dimensional numerical and experimental study for flow between two horizontal coaxial cylinders at Rayleigh numbers approaching and exceeding the critical values. The instability and bifurcation phenomena in the natural convection in a horizontal annulus have also attracted the interest of many researchers. Yoo et al. [6][7][8] studied the effect of the Prandtl number on the bifurcation phenomenon of natural convection in an annulus with a narrow or wide gap. Results show that natural convection with high Prandtl number in a narrow annulus has bifurcation solutions. Petrone et al. [9,10] analyzed the bifurcation phenomena theoretically through a bifurcation diagram by numerical simulation. Luo et al. [11] reported the eccentricity effect on the dual solutions. Hu et al. [12] demonstrated three solutions in the wide gap annulus with a constant heat flux at wall. Zhang et al. [13] investigated the natural convection in a circular enclosure partitioned by a thin flat plate. Natural convection in other types of enclosures has also been studied by many other researchers [14,15].
It has been confirmed that the rotation effect plays a significant role in natural convection in an annulus in applications such as the heat transfer in bearing system, i.e., in the work of Lopez et al. [16] which studied the rotation effect on convection in a vertical annulus. For a square enclosure within which there is a rotating circular cylinder, Zhang et al. [17] found that the heat transfer can be conduction-dominated or mixed convection-dominated depending on the rotating speed. In this paper, we focus on the rotation effect on the natural convection in a horizontal annulus. The mechanism of the rotation effect is quantified by linear speed of rotational inner cylinder, and is presented and analyzed by the streamlines and isotherms at different dimensionless linear speeds. We also use the bifurcation diagram of rotation system to discuss the bifurcation solutions of natural convection in a horizontal annulus with stationary inner and outer cylinders. The present study of the rotation effect on natural convection provides a thorough insight in understanding the heat transfer in the bearing system problem.
The remainder of this paper is organized as follows. Section 2 describes the scope of simulations and briefly introduces the numerical method used in this paper. In Section 3, numerical results are compared with experiments. In Section 4, numerical results and the effect of rotation on the natural convection are presented. We summarize our findings and conclude in Section 5.

Numerical Model and Method
In this section, the configuration of natural convection in a horizontal annulus and a treatment for thermal flow with curved boundaries are introduced. We follow the idea of the coupled lattice BGK model (CLBGK) [21] to treat thermal flow field. The IB-LBM is used to treat the boundary conditions.

Numerical Model
The configuration of a two-dimensional natural convection is shown in Figure 1. R i and R o denote the radius of inner and outer cylinders, respectively. The flow is driven by gravity G α . The governing equations of the convection are [23]: where ρ is the fluid density, u α is velocity, G α is the buoyancy, F α is the body force, µ is the dynamic viscosity, T is the temperature, κ is thermal diffusivity and Φ is the heat source term. The inner and outer cylinder surfaces are maintained at different uniform temperatures T h and T c (T h > T c ), respectively. With the Boussinesq approximation [21], the fluid density is assumed to be a linear function of the temperature and the buoyancy can be written as: G α = −g α β(T − T 0 ) and β is the coefficient of thermal expansion and T 0 is the average temperature of fluid.
where ρ is the fluid density, α u is velocity, α G is the buoyancy, α F is the body force, µ is the dynamic viscosity, T is the temperature, κ is thermal diffusivity and Φ is the heat source term. The inner and outer cylinder surfaces are maintained at different uniform temperatures h T and c T ( c h T T > ), respectively. With the Boussinesq approximation [21], the fluid density is assumed to be a linear function of the temperature and the buoyancy can be written as: and β is the coefficient of thermal expansion and 0 T is the average temperature of fluid.
where c h is temperature difference between hot wall and cold wall, L is the characteristic length across the temperature field which is where Pr is the Prandtl number defined as where | | α q is the quantity of local heat flux and c Q is the average heat flux resulting from pure conduction.
where ∆T = T h − T c is temperature difference between hot wall and cold wall, L is the characteristic length across the temperature field which is R o − R i for the convection in a horizontal annulus. The Mach number used in present work is fixed at 0.01 to ensure that the flow is incompressible. The Rayleigh number (Ra): where Pr is the Prandtl number defined as Pr = ν/κ. The Nusselt number (Nu): where |q α | is the quantity of local heat flux and Q c is the average heat flux resulting from pure conduction.

Thermal IB-LBM
To solve the governing equations, the CLBGK model [23] is used in this paper. The general philosophy is to solve the velocity field and the temperature field using two independent lattice BGK equations, respectively, and then to combine them into one coupled model for the Boussinesq incompressible flow.
The lattice BGK equation with external force term reads: where f i is the density distribution function, f (eq) i is the equilibrium distribution function, c iα is the ith discretized velocity vector and ω f is the relaxation parameter. For the D2Q9 model [18], the equilibrium function is: where sound speed c s and the weight coefficient w i are: Using Chapman-Enskog expansion till the second order of approximation and considering the conserved quantity ρ = ∑ f i and ρu α = ∑ f i c iα + F α δ t /2, the relation between the relaxation parameter ω f and kinematic viscosity ν can be obtained as: The force scheme proposed by Buick [27] is applied in this work: where G α is the external force.
Coupling temperature field to flow field, the well-known Boussinesq approximation is often used in the study of natural convection. With the Boussinesq approximation, the gravity can be rewritten as: To simplify this formula, T 0 is chosen as −1/β and the gravity is: For the flow field, non-slip boundary condition is considered in all simulations. The direct force IB-LBM proposed by Feng [22] is applied to realize the dynamical boundary condition. The discretion of the derivative of time in above equation is related to the desired boundary condition (Equation (12)). This method is also termed as direct forcing model by Wu [25] penalty method by Feng et al. [21] and direct forcing method by Feng et al. 2005 [22]. They proposed a direct force term to realize the no-slip boundary condition as formed below.
Entropy 2018, 20, 733 and u D α is the desired value of velocity on the boundary. Based on their method, the force term used in a thermal IB-LBM has the following formula: where e yα is the unit vector in the y direction.
To simulate the thermal field, we follow the idea of CLBGK proposed by Guo [23] and use a 9-discrete-velocity model with the distribution function of temperature g i .
where Φ i is discretized heat source term which is related to heat source Φ. The thermal equilibrium function is: To compare with Equation (1) (c), the Chapman-Enskog expansion is also used to decide the relation between ω g and κ, which is: Following the idea of previous direct force model and considering u α ∂ α T = 0 at boundary, the isothermal boundary condition is realized by the following modified heat source: and T D is the desired value of temperature on the boundary. Since LBM deals with Eulerian flow points and boundaries are represented by Lagrangian points, Delta function is used to exchange the message between Eulerian points and Lagrangian points and has the following formula: The messages from both velocity field and thermal field are exchanged by Delta function mentioned in the above section.

Verification
The accuracy of the present method in natural convection in an annulus is verified. To qualify the results, the average Nusselt number (Nu)is defined by the following equation: where S is the area of the flow filed. For pure conduction, the exact solution of the temperature field in an annulus is: . Considering the definition of the average Nusselt number, its formula in an annulus yields: where the modified parameter γ is: The relation between γ and R β is shown in Figure 2.
field in an annulus is: (20) where radius ratio  (21) where the modified parameter γ is: The relation between γ and β R is shown in Figure 2. to verify the grid independence and the accuracy of the boundary treatment by comparing with the results from experiments [3]. Table 1 presents the results of average Nusselt number. It is shown that numerical results of the present method are in excellent agreement with experimental results. The relative errors between grid size of . and extrapolated results are only 1.39%, 1.02%, 1.29% and 1.17% The grid independence tests of present thermal IB-LBM for natural convection in horizontal annuli are conducted. The grid sizes of 128 × 128, 256 × 256 and 512 × 512(l.u.) are tested for 4 different Rayleigh numbers, namely, Ra = 10 2 , 10 3 , 10 4 and 2 × 10 4 to verify the grid independence and the accuracy of the boundary treatment by comparing with the results from experiments [3]. Table 1 presents the results of average Nusselt number. It is shown that numerical results of the present method are in excellent agreement with experimental results. The relative errors between grid size of 256 × 256(l.u.) and extrapolated results are only 1.39%, 1.02%, 1.29% and 1.17% for Ra = 10 2 , 10 3 , 10 4 and 5 × 10 4 , respectively. Results show that the grid size of 256 × 256(l.u.) which is used in following studies is fine enough.

Results and Discussion
In this section, the rotation effect on bifurcations of natural convection is analyzed with Pr = 0.7 and Ra = 10 4 in an annulus with its radius ratio R β = 1.5. The streamlines, isotherms and phase portrait of average Nusselt number will be further analyzed. In order to study the rotation effect of the inner cylinder, the linear speed of the inner cylinder is non-dimensionalized by the characteristic speed U 0 = |g α |β∆TL . As a matter of convenience, the Nusselt number Nu is used to represent the average Nusselt number in the flow field Nu in the following section.

Bifurcation Solutions in Stationary Concentric Cylinders
It has been found by many researchers that the bifurcation phenomenon appears in the natural convection in a horizontal annulus [7,11,12,31]. Different initial conditions will lead to different structures of convection. The schematic of different initial conditions used in our paper is shown in Figure 3. With different temperature distributions at Γ u , we obtained three different structure of convection. From isotherms, different numbers of hot peaks can be found in different states. In Figure 4c, for instance, one can find three hot peaks at (0.2, 0.65), (0.5, 0.8) and (0.8, 0.65), respectively. A number of hot peaks are used to distinguish these different states in this paper. In this subsection, the bifurcation phenomena are simulated, which mainly verify the accuracy of the thermal IB-LBM treatment and obtain the initial conditions for following rotating cases. consists of three hot peaks branches. 3-peak-state is presented in Figure 4c. A new pair of vortices appears at the top region. It is further observed that new vortex pair stretches and compresses the previous vortex pair and the previous vortex pair turns into two small vortices which beside the newly appeared large vortex pair. One also sees that 3-peak-state is very unstable from Figure 5c. This also suggests that with a small perturbation (inner cylinder rotate clockwise with its dimensionless linear speed exceeds 0.03), the flow state turns from 3-peak-state to 2-peak-state. It is further found that new pairs of vortices constantly appear at the top region of the annulus and the   The bifurcation solutions of natural convection in a horizontal annulus with stationary inner and outer cylinders are presented. We first use zero velocity and low temperature for all flow fields as zero initial condition and obtained 1-peak state. Then 2-peak state is obtained by using the final state of 1-peak state by removing internal energy step as its initial condition. Finally, 3-peak state is obtained by using the final state of 2-peak state with an injecting internal energy step as its initial condition. The treatment of removing/injecting internal energy step in dealing with the temperature field at top region Γ u is shown in Figure 3. When implementing removing internal energy step, temperature at top region Γ u is set as the lowest temperature T c , compulsively. The injecting internal energy step is to set the temperature at top region the highest temperature T h , compulsively. Figure 4 shows the results of streamlines and isotherms at Ra = 10 4 , Pr = 0.7, R β = 1.5. For 1-peak-state which is shown in Figure 4a, it is observed with zero initialization. It is seen that a pair of symmetrical vortices are formed in the annulus because of the buoyancy and one hot peak in the isotherms appears at the top region. According to the formula of local Nusselt number (Equation (4)), it is larger at top region than the remained region. After removing internal energy step, another solution with two hot peaks of this system is observed in Figure 4b. It is seen that a new pair of vortices appear at the top region of the annulus and velocity at middle of the top region turns from upwards to downwards. Injecting internal energy to Γ u region at 2-peak state, another bifurcation consists of three hot peaks branches. 3-peak-state is presented in Figure 4c. A new pair of vortices appears at the top region. It is further observed that new vortex pair stretches and compresses the previous vortex pair and the previous vortex pair turns into two small vortices which beside the newly appeared large vortex pair. One also sees that 3-peak-state is very unstable from Figure 5c. This also suggests that with a small perturbation (inner cylinder rotate clockwise with its dimensionless linear speed exceeds 0.03), the flow state turns from 3-peak-state to 2-peak-state. It is further found that new pairs of vortices constantly appear at the top region of the annulus and the Nusselt number gradually increases with increasing number of hot peaks at the same Ra, Pr and R β , which are 1.579, 1.791 and 1.888, respectively.    The phase portrait of Nusselt number is shown in Figure 5. In order to show the detail in the phase portrait, final phase trajectory of Nusselt number are zoomed in and presented in the inset of Figure 5. From Figure 5a,b, the phase portraits indicate that the Nusselt number increases with increasing hot peaks. Since the configuration is symmetric, the phase portraits in both Figure 5a,b indicate that the 2-peak-state and 3-peak-state are both stable. However, as discussed above, 3-peak-state will be unstable when the symmetry of the system breaks. Figure 5c shows the phase trajectory of the progress when a small perturbation is introduced into 3-peak-state. It is seen that any perturbation which breaks the symmetry of the system makes one of two previous small vortices compressed, the average Nusselt number decreases, and the compressed vortex finally disappears. Two large vortices beside the disappeared small vortex then merged into one vortex and the system turns to 2-peak-state. One also sees that the average Nusselt number of the newly formed 2-peak-state is smaller than that of the equilibrium system, vortices then resize under the buoyancy and the average Nusselt number increases till the system reaches equilibrium state. (a) Flow state turns from 1-peak-state to 2-peak-state after removing energy step; (b) Flow state turns from 2-peak-state to 3-peak-state after injecting energy step; (c) Flow state turns from 3-peak-state back to 2-peak-state after a small perturbation introduced into the system.

(c) 3-Peak-State
The phase portrait of Nusselt number is shown in Figure 5. In order to show the detail in the phase portrait, final phase trajectory of Nusselt number are zoomed in and presented in the inset of Figure 5. From Figure 5a,b, the phase portraits indicate that the Nusselt number increases with increasing hot peaks. Since the configuration is symmetric, the phase portraits in both Figure 5a,b indicate that the 2-peak-state and 3-peak-state are both stable. However, as discussed above, 3-peakstate will be unstable when the symmetry of the system breaks. Figure 5c shows the phase trajectory of the progress when a small perturbation is introduced into 3-peak-state. It is seen that any perturbation which breaks the symmetry of the system makes one of two previous small vortices compressed, the average Nusselt number decreases, and the compressed vortex finally disappears. Two large vortices beside the disappeared small vortex then merged into one vortex and the system Figure 5. The phase portrait of Nusselt number. (a) Flow state turns from 1-peak-state to 2-peak-state after removing energy step; (b) Flow state turns from 2-peak-state to 3-peak-state after injecting energy step; (c) Flow state turns from 3-peak-state back to 2-peak-state after a small perturbation introduced into the system.

Rotation Effect on Natural Convection
The effect of linear speed of rotational inner cylinder is further analyzed in this subsection. Since 3-peak-state is an unstable system for rotating inner cylinder, 1-peak-state and 2-peak-state are the only two bifurcation systems discussed in this subsection. The inner cylinder rotates clockwise and its dimensionless linear speed (U * i ) is ranging from 0 to 4. Figures 6a and 7a show the results obtained from stationary inner cylinder of 1-peak-state and 2-peak-state, respectively. The linear speed is increasingly considered after that.
3-peak-state is an unstable system for rotating inner cylinder, 1-peak-state and 2-peak-state are the only two bifurcation systems discussed in this subsection. The inner cylinder rotates clockwise and its dimensionless linear speed ( * i U ) is ranging from 0 to 4. Figures 6a and 7a show the results obtained from stationary inner cylinder of 1-peak-state and 2-peak-state, respectively. The linear speed is increasingly considered after that.     Figure 6. The streamlines and isotherms of rotation effect on 1-peak-state natural convection for It is seen that for 1-peak-state, two symmetrical vortices are formed in a stationary annulus as shown in Figure 6a. With increasing linear speed, left vortex is stretched and the right one is compressed, the hot peak moves clockwise and the width of hot peak increases as shown in Figure 6b. One also sees that the separated part of convection at top region deviates from the middle line more It is seen that for 1-peak-state, two symmetrical vortices are formed in a stationary annulus as shown in Figure 6a. With increasing linear speed, left vortex is stretched and the right one is compressed, the hot peak moves clockwise and the width of hot peak increases as shown in Figure 6b. One also sees that the separated part of convection at top region deviates from the middle line more obviously than that at bottom region. Left vortex wedges into the right one at top region and compresses the right vortex. The right vortex then approaches the outer cylinder. From Figure 6c, three main convection appear in the annulus. It is obtained that the first convection is the left vortex and the second one is the right vortex. This suggests that the third convection starts from right part of bottom region near the right vortex, clings to the hot inner cylinder, streams around the right vortex and finally back to the right part of bottom region. From Figure 6d, one also sees that the first convection, the left vortex is compressed and is located in the top region, the second convection, and the right main vortex disappears, the third convection is developed into a main convection in the annulus when the dimensionless linear speed is faster than U * i = 4.0. This mainly implies that the disappearance of the second convection leads to a rapid decrease of the Nusselt number after the critical number of dimensionless linear speed, which is between 3 and 4.
For 2-peak-state, a stable range of 0 ≤ U * i ≤ 0.28 is considered, in which 2-peak-state is stable. The streamlines and isotherms of rotation effect on 2-peak-state are shown in Figure 7, where four dimensionless speed, namely, U * i =0, 0.1, 0.2, 0.28 are considered. As shown in Figure 7a, one can see that a pair of main vortices and a pair of small vortices at top region emerge symmetrically in the annulus at stationary condition. From Figure 7b,c, it is seen that as linear speed of the inner cylinder increases, the symmetry of the flow are broken, the left vortices of both pairs are stretched and the right ones are compressed. It is also obtained that two hot peaks move clockwise and the width of left hot peak becomes larger than the right one. Figure 7d shows the result from highest linear speed before critical speed. It is seen that the compressed vortex of the small pair approaches the vertical position of the annulus. Amazingly, it is further found that the compressed small vortex disappears, and the stretched small vortex and the compressed main vortex merge with each other; the system then turns to 1-peak-state as the dimensionless linear speed is larger than a critical number. Figure 8 illustrates the phase portrait of the progress of 2-peak-state turning to 1-peak-state. In order to show the detail in the phase portrait, final phase trajectory of Nusselt number are zoomed and presented in the inset of Figure 8. The steady 2-peak-state flow field with U * i = 0.28 is used as the initial condition of this simulation, the speed of inner cylinder is increased to 0.29. After the non-equilibrium progress of small vortex disappearing and two vortex merging with each other, as analyzed before. This also suggests that the Nusselt number becomes that of steady 1-peak-state with inner cylinder speed equals 0.29 gradually. It is further found that the critical number of dimensionless linear speed is between 0.28 and 0.29, the bifurcation phenomenon disappears and the steady state is 1-peak-state when the linear speed is higher than the critical speed.
As shown in Figure 9, both duel solutions can be obtained till U * i = 0.28. With the increment of linear speed, the rotation effect leads to the decrement of the average Nusselt number, and it is more obvious in 2-peak-state than that in 1-peak-state, which mainly illustrates that 2-peak-state is more sensitive than 1-peak-state. It is obtained that the average Nusselt number of 2-peak-state is larger than that of 1-peak-state, i.e., heat convection is more important in 2-peak-state than that in 1-peak-state. It is implied that when U * i is greater than 0.29, 2-peak-state is unstable. A local minimum value of average Nusselt number which appears at U * i = 0.5, indicating two different main mechanisms driving the convection. For U * i ≤ 0.5, the system is over-damped for U * i ≥ 0.6, the system is under-damped. Figure 9 shows the phase trajectories of cases at U * i = 0.3 and U * i = 0.9 as examples of these two systems. This also suggests that the third convection seems leading to the behavior of under-damped trajectories while the first and the second convection seems leading to the behavior of over-damped trajectories. It is further found that a small increase occurs for the Nusselt number, when the third convection is taken control, in which this increment results in the local minimum value of average Nusselt number at around U * i = 0.5. Entropy 2018, 19, x 13 of 16 Figure 8. The phase portrait of 2-peak-state turns to 1-peak-state.
As shown in Figure 9, both duel solutions can be obtained till With the increment of linear speed, the rotation effect leads to the decrement of the average Nusselt number, and it is more obvious in 2-peak-state than that in 1-peak-state, which mainly illustrates that 2-peak-state is more sensitive than 1-peak-state. It is obtained that the average Nusselt number of 2-peak-state is larger than that of 1-peak-state, i.e., heat convection is more important in 2-peak-state than that in 1peak-state. It is implied that when   As shown in Figure 9, both duel solutions can be obtained till With the increment of linear speed, the rotation effect leads to the decrement of the average Nusselt number, and it is more obvious in 2-peak-state than that in 1-peak-state, which mainly illustrates that 2-peak-state is more sensitive than 1-peak-state. It is obtained that the average Nusselt number of 2-peak-state is larger than that of 1-peak-state, i.e., heat convection is more important in 2-peak-state than that in 1peak-state. It is implied that when

Conclusions
In this paper, numerical simulations have been performed to demonstrate the existence of bifurcation phenomena in natural convection in a horizontal annulus to investigate the rotation effect on it by using the thermal IB-LBM. Several conclusions can be summarized as follows.
In the first place, it is validated that numerical results of the present method are well consistent with experimental results. Our results manifest that 1-peak-state is the most stable system, 2-peak-state is more sensitive than 1-peak-state and 3-peak-state is the most sensitive system which will breakdown by a very small perturbation.
After that it is found that the rotation effect on 2-peak-state is more obvious than 1-peak-state, since 2-peak-state is more sensitive than 1-peak-state. It is noted that in the range of the present parameters, 2-peak-state is stable. When the linear speed of inner cylinder exceeds the critical speed, 2-peak-state breakdown and turns to 1-peak-state. From the flow pattern shown by streamlines, one may indicate that the right small vortex is compressed and finally disappears, then the left small vortex and the right main vortex merge with each other and the flow turns to 1-peak-state. Finally, three convections in this system which are affected by different linear speed of inner cylinder were discussed in detail. It is implied that the first and the second convections are contributed by the left and right main vortices, respectively. The third convection is a new vortex under the effect of rotation. With increasing linear speed, left main vortex is stretched and period of the first convection increases. Right main vortex is compressed and the scope of the second convection decreases. The third convection grows with the acceleration of rotation. It is further found that 2-peak-state is stable in the range of 0 ≤ U * i ≤ 0.28, and when the linear speed of inner cylinder exceeds the critical speed, 2-peak-state breakdown and turns to 1-peak-state.
Author Contributions: Y.W. contributed significantly to analysis and manuscript preparation. Z.W. performed the data analyses and wrote the manuscript. Y.Q. contributed the concept of the study. W.G. helped perform the analysis with constructive discussions.