Transitions and bifurcations of transient natural convection in a horizontal annulus with radiatively participating medium are numerically investigated using the coupled lattice Boltzmann and direct collocation meshless (LB-DCM) method. As a hybrid approach based on a common multi-scale Boltzmann-type model, the LB-DCM scheme is easy to implement and has an excellent flexibility in dealing with the irregular geometries. Separate particle distribution functions in the LBM are used to calculate the density field, the velocity field and the thermal field. In the radiatively participating medium, the contribution of thermal radiation to natural convection must be taken into account, and it is considered as a radiative term in the energy equation that is solved by the meshless method with moving least-squares (MLS) approximation. The occurrence of various instabilities and bifurcative phenomena is analyzed for different Rayleigh number *Ra* and Prandtl number *Pr* with and without radiation. Then, bifurcation diagrams and dual solutions are presented for relevant radiative parameters, such as convection-radiation parameter *Rc* and optical thickness *τ*. Numerical results show that the presence of volumetric radiation changes the static temperature gradient of the fluid, and generally results in an increase in the flow critical value. Besides, the existence and development of dual solutions of transient convection in the presence of radiation are greatly affected by radiative parameters. Finally, the advantage of LB-DCM combination is discussed, and the potential benefits of applying the LB-DCM method to multi-field coupling problems are demonstrated.

## I. INTRODUCTION

Natural convection in a horizontal annulus is important in many applications such as thermal energy storage systems, cooling of electronic components and transmission cables, etc. Due to its simple geometry and well-defined boundary conditions, the system has been studied extensively by many researchers. Kuhen and Goldsein^{1} conducted an experimental and theoretical investigation of natural convection in concentric and eccentric horizontal cylindrical annuli. Mojtabi *et al.*^{2} discussed the energy stability of a natural convective flow in a horizontal annular space. Wang *et al.*^{3} numerically investigated the low Rayleigh number convection in a horizontal annulus. Guj and Stella^{4} reported the numerical and experimental buoyancy driven flow in horizontal eccentric annuli. They studied the effect of the horizontal eccentricity and found that the average Nusselt number is nearly independent of the horizontal eccentricity. Kuo *et al.*^{5} numerically studied the combined natural convection and radiation in a two-dimensional horizontal annulus filled with a radiatively participating gray medium using a control-volume-based finite difference method and a spectral collocation method coupled with an influence matrix technique. Dyko *et al.*^{6} presented a numerical study of three-dimensional supercritical states of natural convection in a narrow-gap horizontal cylindrical annulus. Wordsworth *et al.*^{7} experimentally studied the turbulent flow in a differentially heated rotating annulus. Carrasco-Teja *et al.*^{8} analyzed the effects of rotation and axial motion of the inner cylinder of an eccentric annular duct during the displacement flow between two Newtonian fluids of differing density and viscosity. Fattahi *et al.*^{9} simulated the natural and mixed convection heat transfer problems by lattice Boltzmann model based on double-population approach. Cao *et al.*^{10} investigated the natural convection in eccentric annulus over a wide range of Rayleigh numbers using the dissipative particle dynamics model with energy conservation (eDPD). Besides, natural convection with volumetric radiation has been investigated by different methods. Lari *et al.*^{11,12} analyzed the combined heat transfer of radiation and natural convection in a square cavity containing participating gases by the combined finite volume and discrete ordinates method. Moufekkir *et al.*^{13,14} investigated the heat transfer by natural convection and radiation in an enclosure filled with an isotropic scattering medium, using the coupled multiple relaxation time lattice Boltzmann method and finite difference method.

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 (*Pr* = 0.7) by visualizing flow patterns, and categorized the flow patterns obtained by their experiments and accumulated results by other researchers in a parameter space of (*Ra*, *R*) where *R* is radius ratio. This experiment was repeated and the classification was confirmed by Dyko *et al.*^{16} Busse^{17} 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 *et al.*^{18} analyzed the instabilities in three dimensional differentially heated cavities with adiabatic horizontal walls. Yoo,^{19} in his outstanding work, reported the occurrences of dual solutions at *Ra* larger than a critical value. He observed dual steady solutions at *Ra* > *Ra*_{c} ≈ 3800 for wide gap annuli (*R* = 2) by the vorticity-streamfunction method. Similar results were provided later by Mizushima *et al.*,^{20,21} xin *et al.*^{22} and Mercader *et al.*^{23} Petrone *et al.*^{24} 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 *Pr* = 0.7. Angeli *et al.*^{25} provided a critical review of buoyancy-induced flow transitions in horizontal annuli. Soucasse *et al.*^{26} studied the transitional regimes of natural convection in a differentially heated cubical cavity under the effects of wall and molecular gas radiation.

As mentioned above, much work has been done for the bifurcation phenomena of natural convection. However, few studies have been concerned about the effect of volumetric radiation on flow instability and bifurcation. One major reason for this lack of research is that the inclusion of thermal radiation increases the complexity of the analysis greatly. The original complex equations of motion for pure convection become even more complicated with the addition of the integro-differential radiative transfer equation. Besides, transient process must be considered to present the onset and development of bifurcation, which will obviously increase the computational cost. Thanks to the inherently transient property, the lattice Boltzmann method is an excellent approach for the investigation of thermal and hydrodynamic instabilities and bifurcations including unsteady or periodic process. In the present work, we develop the coupled lattice Boltzmann and direct collocation meshless method to simulate multi-field coupling problems.

The lattice Boltzmann method (LBM) is a powerful numerical technique based on kinetic theory for simulating fluid flows and modeling the physical process in fluids. Owing to its mesoscopic origin, in comparison with conventional fluid dynamics solvers, it offers such advantages as simple calculation procedure, simple and efficient implementation for parallel computation, and easy and robust handling of complex geometries. In 1998, Chen and Doolen^{27} first applied the LBM to the fluid flow problems. Then, it was developed to solve a wide range of heat transfer problems.^{28,29} Quite recently, natural convection problems have also been investigated by the LBM.^{30} Besides, LBM can be coupled to a conventional numerical scheme to provide a convenient way to simulate multi-field coupling problems such as the thermo-convective problems, thermo-magnetic systems, double diffusion models, coupled convection and radiation heat transfer problems, etc. In recent years, three different combinations have been proposed, namely, the combination of lattice-Boltzmann and finite difference method (LB-FDM),^{13,14,31} the hybrid technique of lattice-Boltzmann and finite volume method (LB-FVM)^{32–35} and the combination of lattice-Boltzmann and finite element method (LB-FEM).^{36} One of the main advantages of LBM is that it is easy and robust to deal with the natural convection problems in complex geometries. To simulate the combined heat transfer of radiation and natural convection in irregular geometries, in addition to the LBM for natural convection, we need to find another efficient approach to solve the radiative transfer equation (RTE).

Many numerical methods have been developed to solve radiative transfer in participating media in recent years, such as Monte Carlo method,^{37} discrete transfer method,^{38} discrete ordinates method (DOM),^{39} finite volume method (FVM),^{40} finite element method (FEM)^{41} and meshless method.^{42,43} Among them, the direct collocation meshless (DCM) method is an extremely good choice to couple with LBM, for the reason that the DCM method is simple and easy to implement without any additional grids. Besides, the DCM method is also easy to deal with irregular geometries.

The objective of the present work is to use the LB-DCM method to investigate the bifurcation and dual solutions in transient natural convection with consideration of volumetric radiation in a horizontal annulus, mainly focusing on the analyses of the effects of various parameters such as convection-radiation parameter and optical thickness on thermal and hydrodynamic instabilities and bifurcations. The remainder of the paper is organized as follows. In the next section, the lattice Boltzmann equations for natural convection is briefly outlined. Separate particle distribution functions in the LBM are used to calculate the density field, the velocity field and the thermal field, and the curved boundary treatment is detailed presented. In Section III, the direct collocation meshless method with MLS approximation is introduced for solving the RTE. The validation and application of the LB-DCM method are provided in Section V. Finally, the conclusions are drawn in Section VI.

## II. LATTICE BOLTZMANN EQUATIONS FOR NATURAL CONVECTION

### A. Physical model and computational domain

The considered configuration is shown in Fig. 1. Its inner and outer walls are at temperatures *T*_{h} and *T*_{c}, respectively. The radius ratio *R* is assumed to be 2 in this work. Thermophysical properties, except density, are assumed constant. Density is considered to vary in the Boussinesq sense. The homogeneous medium in the annulus is absorbing, emitting, and scattering, and its optical properties, such as the extinction coefficient *β* and scattering albedo *ω*, are constant. The boundaries of the annulus are considered to be diffuse and gray.

Fig. 2 shows one-quarter of the computational physical domain. For the LBM, the lattice cells are the basic computational units in which the fluid particles undergo collision and streaming steps. For the DCM, the lattice nodes are taken as collocation points and used for moving least-squares (MLS) approximation. Therefore, no additional nodes and grids are required for the hybrid LB-DCM approach.

### B. Lattice Boltzmann equation for density and velocity fields

In the LBM, separate particle distribution functions are used to compute the density (and velocity) and thermal fields. The distribution function is obtained by solving the lattice-Boltzmann equation which is a special discretization of the kinetic Boltzmann equation. The macroscopic quantities of the simulated fluid can then be derived by calculating the hydrodynamic moments of the distribution function. Unlike the second-order PDEs in the Navier-Stokes approach, the LBM uses only the first order PDEs, and the discretized form of LB equation can be expressed as^{9,27}

where *f*_{α} is the distribution function in the direction *α*, *e*_{α}is the microscopic velocity, τ_{ν} is the relaxation time, *v* is the kinematical viscosity, |$f_\alpha ^{eq} $|$f\alpha eq$is the equilibrium distribution function, and *F*_{α}is the external force in the direction of lattice velocity.

In this work the D2Q9 model is used for the discretization of velocity space for the two dimensional case. The nine velocities in the D2Q9 (Fig. 2) lattices are given by

in which, the streaming speed is defined as *c* = Δ*x*/Δ*t*, where Δ*x* and Δ*t* are the lattice cell size and the lattice time step, respectively.

For the D2Q9 lattices (Fig. 2) used in the present work, the relaxation time τ_{ν} is computed from^{9,27}

The equilibrium distribution function of Eq. (1) is expressed as

where *w*_{a} is the equilibrium distribution weight for direction *α*, given as

The fluid mass density *ρ* can be evaluated from Eq. (6a), whereas the velocity ** u** is contained in the momentum density of Eq. (6b)

In order to incorporate buoyancy force in the model, the Boussinesq approximation is applied. Therefore, the force term in Eq. (1) needs to be calculated as below in vertical direction (y):

where *g*_{y} is the acceleration of gravity acting in the y-direction of the lattice links; *β* is the thermal expansion coefficient; |${\bm e}_{\alpha y} $|$e\alpha y$ is the y-component of |${\bm e}_\alpha $|$e\alpha $; the terms ρ(**x**, *t*) and *T*(**x**, *t*) are the density and temperature, respectively.

### C. Lattice Boltzmann equation for thermal field

For the computation of thermal field, in the presence of volumetric radiation, the lattice Boltzmann equation governing the temperature field can be given by^{28,30}

where *g*_{α} is the temperature distribution function in the *α* direction, τ_{T} is the relaxation time, |$g_\alpha ^{eq} $|$g\alpha eq$ is the equilibrium distribution function, and *q*_{R} is the radiative heat flux.

In Eq. (8), τ_{T} is represented as

in which, *χ* is the thermal diffusion coefficient. Equilibrium distribution function for temperature field can be expressed as

where *T* is the fluid temperature and can be evaluated from

In the case of natural convection, the flow and the temperature fields depend on two non-dimensional parameters, namely, the Prandtl number and Rayleigh number presented in Eq. (12) and Eq. (13), respectively,

where Δ*T* is the temperature difference and *l* (= *R*_{o}-*R*_{i}) is the characteristic length of the computational domain. With the inclusion of volumetric radiation, another non-dimensional parameter, the convection-radiation parameter (*Rc*), is given by

where *k* is the thermal conductivity, *σ* is the Stefan Boltzmann constant, and *T*_{ref} (= (*T*_{h} + *T*_{c})/2) is the reference temperature.

### D. Curved boundary treatment in lattice Boltzmann method

In Fig. 3, a part of curved wall separates the solid region from the fluid region. The lattice node on the fluid side of the boundary is denoted as |${\bm x}_{f} $|$xf$ and that on the solid side is denoted as |${\bm x}_b $|$xb$. The hollow small circles on the boundary, |${\bm x}_w $|$xw$, denote the intersections of the wall with various lattice links. The boundary velocity at |${\bm x}_w $|$xw$ is |${\bm u}_w $|$uw$. The fraction of an intersected link in the fluid region, denoted by Δ, is defined as

The |$\tilde f_{\bar \alpha } ({\bm x}_b,t)$|$f\u0303\alpha \xaf(xb,t)$ and |$\tilde g_{\bar \alpha } ({\bm x}_b,t)$|$g\u0303\alpha \xaf(xb,t)$ that are post-collision distribution functions coming from a solid node |${\bm x}_b $|$xb$ to a fluid node |${\bm x}_{f} $|$xf$, are not known. To complete the subsequent streaming step, |$\tilde f_{\bar \alpha } ({\bm x}_b,t)$|$f\u0303\alpha \xaf(xb,t)$ and |$\tilde g_{\bar \alpha } ({\bm x}_b,t)$|$g\u0303\alpha \xaf(xb,t)$ are required to be obtained since they give |$\tilde f_{\bar \alpha } ({\bm x}_{f},t + \delta t)$|$f\u0303\alpha \xaf(xf,t+\delta t)$ and |$\tilde g_{\bar \alpha } ({\bm x}_{f},t + \delta t)$|$g\u0303\alpha \xaf(xf,t+\delta t)$ after streaming. Many methods are developed to obtain the |$\tilde f_{\bar \alpha } ({\bm x}_b,t)$|$f\u0303\alpha \xaf(xb,t)$ and |$\tilde g_{\bar \alpha } ({\bm x}_b,t)$|$g\u0303\alpha \xaf(xb,t)$, such as the bounce-back scheme,^{44,45} the interpolation scheme,^{46} the non-equilibrium extrapolation scheme,^{47} etc.

For getting the density and velocity fields in the curved boundary, the present work adopts the YMS scheme proposed by Yu *et al.*^{46} that satisfies the no-slip boundary condition and has second order accuracy. The post-collision distribution function can be expressed as

For obtaining the temperature field in the curved boundary we use the non-equilibrium extrapolation scheme introduced by Guo *et al.*^{47} Distribution function for temperature is divided into two parts, equilibrium and non-equilibrium

in which, equilibrium and non-equilibrium parts are define as

with

## III. DIRECT COLLOCATION MESHLESS FOR RTE

### A. Radiative transfer equation

The divergence of radiative heat flux |$\nabla \bullet {\bm q}_R $|$\u2207\u2022qR$ in Eq. (8) is given as

in which, *n* is the refractive index of the participating medium, and the incident radiation *G* distribution can be calculated as

The radiative intensity *I* is obtained by solving the RTE in each ray direction. The discrete-ordinates equation of radiative transfer can be written as

where |${\bm \Omega }^m$|$\Omega m$ is the vector of transport direction, and *S*^{m}is the source term, given by

in which, μ_{m}, η_{m} and ξ_{m} are direction cosines in direction *m*; β = κ_{a} + κ_{s}is the extinction coefficient where κ_{a} and κ_{s} are the absorption and scattering coefficients, respectively; Φ is the scattering phase function; |$w^{m^{\prime} } $|$wm\u2032$ is the weight corresponding to the direction *m*′.

Considering the opaque surface, the radiative boundary condition is given by

in which |$I_w^m $|$Iwm$ is the radiation intensity at the boundary at the direction *m*, ε_{w} is the wall emissivity, |${\bm n}_w$|$nw$ is the unit normal vector of the boundary surface, and |$s^{m^{\prime} } $|$sm\u2032$ is the unit vector in the direction *m*′.

In the application of DCM to solving the RTE, the moving least-squares (MLS) approximation is employed to construct the trial functions by using collocation points, and the DCM is adopted for constructing the linear equations as mentioned in Ref. 48.

### B. Moving least-squares approximation

As shown in Fig. 2, consider a spatial subdomain *V*_{x} in the neighborhood of a point **x** denoted as the domain of the definition of the MLS approximation for the trial function at **x**. To approximate the distribution of a generic trial function *I*^{m}(**x**) in *V*_{x} over a number of local nodes {**x**_{i}} (*i* = 1, 2, …, *n*), the MLS approximation |$\tilde I^m ({\bf x})$|$I\u0303m(x)$ of *I*^{m}(**x**), ∀**x** ∈ *V*_{x}, is given by

where *p*_{j} is the *j*th polynomal basis function and *a*_{j} is the corresponding expansion coefficient. To give the local projection, a local inner product defined in Ref. 48 is expressed as

where *w* is a weight function with compact support. In this article, the complete second order polynomial basis (**p**^{T}(**x**) = [1, *x*, *y*, *x*^{2}, *xy*, *y*^{2}]) is used for MLS approximation. The least square approach can then be formulated with the help of the local inner product as

A solution to the minimization problem under the extreme value condition of ∂*J*_{x}(**a**)/∂*a*_{j} = 0 yields

which can be written in a matrix form as

where

A nodal basis function (at node *i*) from the MLS approximation can then be written as

and the local polynomial approximation [Eq. (25)] can be formulated into the equivalent nodal basis expansion as

The first-order partial derivative of the radiative intensity is required for the solution to the RTE. Here, the evaluation of the partial derivative is by the ‘diffuse derivative’ approach,^{49} and yields

where ∇ is considered as a gradient operator. To define the weight function *w*(*r*) appeared in Eq. (31) and Eq. (34), the 4th order Wendland's function,^{50} given by Eq. (35), is chosen that is compactly supported and continuous up to the 4th order derivatives.

where *r*_{0} = 2*r*/*d*; *r* is the distance between **x**_{i} and **x**; *d* is the diameter of the support domain, represented as *d* = *k*(*p* + 1)*Δh* in which *p* is the order of polynomial basis, *Δh* is the local discrete length, and *k* = 1.5 in this paper.

### C. Collocation-based meshless approach

The basic idea of the DCM approach for solving the radiative transfer equation is to force the residual of the equations to be zeros on each collocation point. By using the MLS approximation mentioned above, and substituting Eq. (33) and its first-order partial derivatives [Eq. (34)] into Eq. (22), we get the discretization forms of the RTE and its boundary condition as

## IV. SOLUTION PROCEDURE AND DIMENSIONLESS PARAMETERS

As mentioned in section I, the simulation of transient combined radiation and convection is always difficult and time consuming. Among many numerical schemes, the LB-DCM method is a good choice with an excellent flexibility in dealing with multi-fields coupling problems and irregular geometries. Besides, the LB-DCM method is stable and efficient.

In the solution of the coupled radiation and convection heat transfer, the LB-DCM is applied as a transient solver with non-dimensional time step of *Δζ* = *χΔt*/*L*^{2} = 10^{−5}. In each time step, the LBM is used to calculate the velocity and temperature fields as mentioned in section II, and the radiative intensity is solved by DCM method as presented in section III. Through coupled iterating, the transient solution of whole field is obtained. Besides, the non-dimensional forms of all variables used in LB-DCM approach are given as follows.

The convective and radiative Nusselt numbers at the inner cylinder are calculated from the corresponding heat fluxes as:

and average Nusselt number *Nu*_{ave} at the wall is the linearly-averaged value of *Nu*. The total Nusselt number is

The maximum relative error of *Nu*_{total} at the inner cylinder is taken as stopping criterion for the global iteration.

## V. RESULTS AND DISCUSSIONS

This section contains three parts. First, the grid independence test and code validation are performed. Then, natural convection and bifurcation with and without the contribution of volumetric radiation are investigated. Finally, bifurcation diagrams are presented for different relevant radiative parameters, namely, convection-radiation parameters *Rc* and optical thickness *τ*. In all calculations, the entire annulus is considered without symmetry assumptions.

### A. Grid independence test and Code verifications

Grid independence of the results has been tested for two different Rayleigh numbers, namely, 10^{4} and 10^{5}. The distributions of local convective Nusselt number along the inner cylinder for different grid numbers are shown in Fig. 4. The radiative parameters are set as *Rc* = 1 and τ = 10. It is seen from Fig. 4(a) that the results suffer from numerical oscillation if grid number is less than 32 × 32. Besides, the maximum relative error between the results for two finest grid discretization of 128 × 128 and 192 × 192 is only about 0.5%, and consequently no further refinement is carried out. Therefore, for *Ra* = 10^{4}, the grid number of 128 × 128 can be considered for the further analysis with low computational cost and good accuracy. However, the required grid number increases to 192 × 192 for *Ra* = 10^{5} as shown in Fig. 4(b).

In addition, the variation of radial velocity at the centerline of the annular gap for different grid numbers is also provided in Fig. 5. It can be seen from Fig. 5 that the velocity profiles are less sensitive to the number of grids than the distributions of convective Nusselt number and give stable results for comparatively less number of grids. Besides, the symmetry of velocity distributions under irregular grids indicates that the LB-DCM algorithm has a good stability and flexibleness.

Because of lack of related research on coupled natural convection and radiation problems in a horizontal annulus, the codes for convection and radiation are separately verified. First, we test the correctness of LBM code for the natural convection without consideration of volumetric radiation, and the results are compared with the benchmark solutions reported by Yoo.^{19} The distributions of the radial velocity of upward flow at the centerline of the annular gap are presented in Fig. 6 for *Pr* = 0.7 and 0.3 at different *Ra*. When *Pr* = 0.7 (Fig. 6(a)), the maximum *U* of upward flow are located at the uppermost point (φ = 0) of the annulus for all *Ra*, i.e., the most strong fluid flow directed upward occurs at the top of the annulus. Accordingly, the upward flow at *Pr* = 0.7 is stable and the transition from upward to downward flow does not occur. However, the maximum *U* of the upward flow for *Pr* = 0.3, does not always increase with *Ra* as shown in Fig. 6(b). Besides, the position of maximum velocity moves upward with the increase of *Ra*, and finally, the flow field is not stable any longer and a transition from upward to downward flow occurs at a critical value of *Ra* = 6.5 × 10^{4}. As presented in Fig. 6, a good agreement is achieved between the numerical results of the present work and those obtained by Yoo,^{19} and the maximum error is within 2%. In addition, the validation of DCM method for radiation has been reported in our previous work.^{43}

### B. Radiation effects on natural convection and bifurcation

In this part, natural convection problems in the presence of volumetric radiation are numerically simulated. In accordance with Yoo's report,^{19} we name the crescent-shaped convection in which the fluid at the top of the annulus ascends as ‘upward flow’, and call the solution branch of the upward flow as ‘U_Branch’. The flow pattern in which the fluid at the top of the annulus descends by forming two counter-rotating eddies in a half annulus is named as ‘downward flow’, and the corresponding solution branch is called as ‘D_Branch’.

Figure 7 presents the temperature fields obtained by LBM-DCM approach for different *Pr* at *Ra* = 4 × 10^{4}, and the results with and without considering volumetric radiation have been compared. For the cases without radiation (Fig. 7(a)), if the simulation is started from the zero initialization, a steady-state downward flow is reached for *Pr* = 0.3, while upward flows are obtained for *Pr* = 0.7 and *Pr* = 1. Then from the Fig. 7(b) we can see that the introduction of radiation will affect the flow structure significantly. In the case of natural convection considering radiation, the fluid close to the lower half of the inner cylinder is still being heated and rises along the inner cylinder. While it is observed that the fluid detaches from the cylinder wall earlier than that in the pure convection case, which causes a downward shift of the centers of circulations. As shown in Fig. 7(b), the pair of thermal plumes becomes more lagging at *Pr* = 0.3 under the effect of radiation; a transition from upward to downward flow occurs at *Pr* = 0.7 and the sole (central) thermal plume becomes less distinct and the isotherms become less distorted at *Pr* = 1. In a word, volumetric radiation influences both the flow type and temperature distribution of natural convection.

Figure 8 plots the non-dimensional radial velocity at the top of center of the annulus (*U*_{0}), showing the effect of thermal radiation on bifurcation of natural convection. Simulations are conducted with and without the effects of volumetric radiation for *Pr* = 0.7, *Rc* = 1 and *τ* = 10. In the absence of radiation, the convection approaches to a steady state upward flow regardless of the initial condition for *Ra* < 2845. However, for *Ra* ≥ 2845 the steady state solution depends upon the initial condition. U_branch (solid line) and D_branch (dashed line) steady flow solutions can be obtained for zero initialization and special initialization, respectively. It is obviously seen that the bifurcation point of natural convection in this case is *Ra*_{S} = 2845. When considering radiation contribution, the thermal energy tends to be redistributed more evenly inside the annulus and consequently smaller *U*_{0} can be observed. In addition, the presence of radiation changes the temperature gradient of the fluid, and results in an increase in the flow critical value for bifurcation from *Ra*_{S} = 2845 to *Ra*_{S} = 6485.

As shown in Fig. 9, dual solutions and bifurcation of natural convection with and without consideration of radiation are visually illustrated by temperature contours, in which isotherms are equal difference with the dimensionless temperature values between 0 and 1. It can be seen from Fig. 9(a) that the U_branch starts at low *Ra* number with zero initialization of the thermal and flow fields. By gradually increasing the value of *Ra*, the isotherms become more curved and the upward flow gets stronger. However, D_branch must starts at high *Ra* number by initializing the flow field to the downward flow, with the decrease of *Ra*, the downward flow gets weaker and finally a transition from downward flow to upward flow occurs at the critical *Ra* value. With the effect of volumetric radiation, isotherms are less distorted and the critical *Ra* is increased as presented in Fig. 9(b). It is concluded that the bifurcation and the existence of dual solution of natural convection are greatly affected by the volumetric radiation.

### C. Bifurcation phenomena for different radiative parameters

In this section we further observe and analyze the instability and bifurcation phenomena of natural convection in radiatively fluid with different radiative parameters, namely, convection-radiation parameter *Rc* and optical thickness *τ*. The enclosure walls are assumed black and the *Pr* is set as 0.71.

Figure 10 shows the bifurcation diagrams in different *Rc* for *τ* = 10. With strong radiation effects (*Rc* = 0.5), the temperature distribution becomes more uniform and is dominated by the diffusive radiative transport when compared with pure natural convection, and the temperature gradients at the cylinder surfaces also become more uniform, resulting in a big increase in the critical Rayleigh number (*Ra*_{S}) as presented in Fig. 10(a). With the increase of *Rc*, radiation effects get weaker, and the *Ra*_{S} is getting close to the critical Rayleigh number of the pure natural convection. Besides, it is seen from Fig. 10(b) ∼ (d) that the center of circulation in flow patterns gets higher as the *Rc* increases.

Figure 11 provides the existence and development of dual solutions in detail. For *Ra* < *Ra*_{S}, same temperature distributions for U_branch and D_Branch can be found regardless of the initial condition. While for *Ra* ≥ *Ra*_{S}, both upward and downward flow can be observed with zero and special initialization, respectively.

Bifurcation diagrams are presented in Fig. 12 for the case of combined natural convection and radiation, with *τ* varying from 1 to 50 for *Rc* = 1. At a high value of *τ* such as *τ* = 50, the effect of thermal radiation is similar to that of conduction, and the change of *Ra*_{S} is small. With the decrease of *τ*, the influence of radiation increases, and the *Ra*_{S} reaches its maximum value at *τ* = 5. However, with further decreasing of *τ*, the radiative effects decrease as radiative energy is less absorbed at lower optical thickness *τ*. With a very small *τ*, the fluid can be regarded as the radiatively non-participating transparent medium.

It is seen from Fig. 13 that the changes of thermal fields for different optical thickness *τ* are quite significant; the thermal plume is weakest and the isotherms are least distorted for *τ* = 5.

## VI. CONCLUSIONS

Radiation effects on transitions and bifurcations of transient natural convection in a horizontal annulus with participating media were numerically investigated by using the hybrid approach of lattice Boltzmann and direct collocation meshless (LB-DCM) method which is inherently transient, easy to implement on lattice grids and has an excellent flexibility in dealing with irregular geometries.

The validation of LB-DCM code was conducted, and the results of the present work were compared with those obtained by other numerical methods reported in the literatures. From the comparisons, we find that this approach is stable and very accurate for the multi-field coupling problems. Then, natural convection without and with volumetric radiation was simulated, and the radiation effects on natural convection in a horizontal annulus were analyzed.

An extensive simulation was carried out to establish the bifurcation diagrams of coupled natural convection and radiation heat transfer. It was concluded that the radiation tends to redistribute the thermal energy more evenly inside the annulus, resulting in an increase of critical Rayleigh number (*Ra*_{S}) with which the bifurcation phenomena in natural convection can be observed. Particularly, the *Ra*_{S} in the case considering strong radiation effects (*Rc* = 0.5, *τ* = 10) is five times of that obtained in pure natural convection. Finally, a large number of bifurcation diagrams were presented with different convection-radiation parameters *Rc* and optical thickness *τ* for the further investigation into the effects of thermal radiation on the instabilities and bifurcation in the natural convection. With the increase of *Rc*, radiation effects become weaker, and the *Ra*_{S} gets closer to the value of pure natural convection. As *τ* increases, the contribution of radiation increases and when *τ* = 5 the *Ra*_{S} reaches its maximum value; while with the further increasing of *τ*, the radiation effect is weakened and consequently the *Ra*_{S} descends.

## ACKNOWLEDGMENTS

This work is supported by the Foundation for Innovative Research Groups of the National Natural Science Foundation of China (Grant No. 51121004) and the National Natural Science Foundation of China (Grant No. 51176040).