Thermal Lattice Boltzmann Flux Solver for Natural Convection of Nanofluid in a Square Enclosure

Abstract

In the present study, mathematical modeling was performed to simulate natural convection of a nanofluid in a square enclosure using the thermal lattice Boltzmann flux solver (TLBFS). Firstly, natural convection in a square enclosure, filled with pure fluid (air and water), was investigated to validate the accuracy and performance of the method. Then, influences of the Rayleigh number, of nanoparticle volume fraction on streamlines, isotherms and average Nusselt number were studied. The numerical results illustrated that heat transfer was enhanced with the augmentation of Rayleigh number and nanoparticle volume fraction. There was a linear relationship between the average Nusselt number and solid volume fraction. and there was an exponential relationship between the average Nusselt number and Ra. In view of the Cartesian grid used by the immersed boundary method and lattice model, the immersed boundary method was chosen to treat the no-slip boundary condition of the flow field, and the Dirichlet boundary condition of the temperature field, to facilitate natural convection around a bluff body in a square enclosure. The presented numerical algorithm and code implementation were validated by means of numerical examples of natural convection between a concentric circular cylinder and a square enclosure at different aspect ratios. Numerical simulations were conducted for natural convection around a cylinder and square in an enclosure. The results illustrated that nanoparticles enhance heat transfer in higher Rayleigh number, and the heat transfer of the inner cylinder is stronger than that of the square at the same perimeter.

1. Introduction

Natural convection has received widespread attention by many researchers because it is relevant to many engineering applications, such as heat exchangers, solar energy and nuclear reactors. Conventional fluids, such as water and ethylene glycol mixture, are not effective heat transfer medias, due to low thermal conductivity. Therefore, nanofluids have gained attention as an alternative and effective heat transfer medium, due to having higher thermal conductivities [1]. There are two main research approaches for studying nanofluids: experiments and numerical simulations. In view of experiments, Song et al. [2] measured the thermal performance of SiC nanofluid in a water pool boiling experiment, and investigated the enhancement for critical heat flux. Nikhah et al. [3] carried out an experimental investigation on the convective boiling of dilute CuO-water nanofluids in an upward flow inside a conventional heat exchanger. Alkasmoul et al. [4] investigated the turbulent flow of Al₂O₃-water, TiO₂-water and CuO-water nanofluids in a heated, horizontal tube with a constant heat flux. The results showed that the efficiency of nanofluids in enhancing heat transfer was not high for turbulent flows. Qi et al. [5] carried out an experimental study on boiling heat transfer of an α-Al₂O₃-water nanofluid.

More researchers have applied numerical methods to study the performance of nanofluids. Khanafer et al. [6] directly solved the macroscopic governing equations to investigate heat transfer enhancement in a two-dimensional enclosure utilizing nanofluids for various pertinent parameters, including Grashof numbers and volume fractions. The results indicated that heat transfer increased with the volumetric fraction of the copper nanoparticles in water at any given Grashof number. Fattahi et al. [7] carried out a study on water-based nanofluid, containing Al₂O₃ or Cu nanoparticles, in a square cavity for Rayleigh number 10³–10⁶ and solid volume fraction 0–0.05, by means of the lattice Boltzmann method. The results indicated that the average Nusselt number increased by increasing the solid volume fraction and the effects of solid volume fraction on Cu were stronger than on Al₂O₃. He et al. [8] applied the single-phase lattice model to simulate convection heat transfer utilizing Al₂O₃-water nanofluid in a square cavity. Qi et al. [9] applied the two-phase lattice Boltzmann model for natural convection of nanofluid. From the above analysis, the lattice Boltzmann method (LBM) has obtained remarkable achievements in simulating incompressible viscous laminar nanoflow. Saadat et al. [10] developed a compressible LB model on standard lattices to solve supersonic flows involving shock waves, based on the consistent D2Q9 LB model, and with the help of appropriate correction terms introduced into the kinetic equations to compensate for deviations in the hydrodynamic limit. Huang et al. [11] improved the lattice Boltzmann model with a self-tuning equation of state to simulate the thermal flows beyond the Boussinesq and ideal-gas approximations. Hosseini et al. [12] derived the appropriate form of the correction term for the space- and time-discretized LB equations, through a Chapman–Enskog analysis for different orders of the equilibrium distribution function. As a mesoscopic approach, LBM can easily solve the macroscopic variables used by distribution functions and the linear streaming and collision processes can effectively simulate the nonlinear convection and diffusion effects in the macroscopic state. With the development of Lattice models in recent years, LBM can solve various flow problems successfully, including incompressible, compressible and thermal flows, by introducing a variety of applicable models. However, the solutions of flow for High Mach number and turbulence problems of complex shape are limited because the standard LBM is strictly limited to using the uniform Cartesian mesh due to the lattice uniformity for flow.

Recently, the idea of coupling the LBM and conventional methods (including finite difference method and finite volume method) has been proposed for computational fluid dynamics. It effectively combines the merits of macroscopic and mesoscopic methods. The coupling algorithm can be divided into the whole region coupling algorithm and the partition coupling algorithm. The whole region coupling algorithm solves the different variables used by different numerical algorithms. Nie et al. [13] and Mezrhab et al. [14] used the LBM-FDM coupling method to solve natural convection problems, in which LBM solved flow problems and FDM analyzed heat transfer. Chen et al. [15] used the LBM-FDM coupling method to solve the two-phase interface convection problem, in which LBM solved the velocity field and FDM solved the concentration field. Mishra et al. [16] used LBM-FVM to solve heat conduction and radiation problems. Sun and zhang [17] used LBM-FVM for conduction and radiation in irregular geometry. The partition coupling algorithm divides the whole region into several sub-regions and realizes the coupling function through information transfer between the sub-regions. Luan et al. [18,19,20] simulated complex flows in porous media using LBM-FVM. Chen et al. [21,22,23] used LBM-FVM to study the multi-scale flow, multi-component mass transfer, proton conduction and electrochemical reaction processes. Li et al. [24,25] used LBM-FVM to study natural convection and the solid–liquid variation problem. Feng et al. [26] developed a thermal lattice Boltzmann model with a hybrid recursive regularization collision operator on standard lattices for simulation of subsonic and sonic compressible flows without shock by LBM-FVM. Essentially, the main advantage of the above two coupling methods is to improve the calculation efficiency of LBM and expand the applications of macroscopic computational fluid dynamics.

A new coupling idea gas been proposed in the past five years. This coupling method adopts the finite volume method to discretize macroscopic governing equations and uses local lattice Boltzmann equation solutions to calculate interface flux, on the basis of considering migration and collision processes. This method realizes the coupling of the macroscopic method and the mesoscopic model and is named the lattice Boltzmann flux solver (LBFS). Yang et al. [27,28] proposed LBFS based on compressible models, which is suitable for calculating viscous and compressible multi-component flows. Shu et al. [29] and Wang et al. [30,31,32] developed LBFS for incompressible viscous flow problems. This method integrates the advantages of the macroscopic method and the mesoscopic model, to not only realize the unified solution of non-viscous flux and viscous flux, but also to improve calculation efficiency without using a uniform grid in the whole calculation domain. Based on the above development, Wang et al. [33] developed the thermal lattice Boltzmann flux solver (TLBFS) and successfully used it to simulate the natural convection problem. Cao [34] proposed a variable property-based lattice Boltzmann flux solver (VPLBFS) for thermal flows with partial or total variation in fluid properties in the low Mach number limit.

In this paper, we attempted to build mathematical modeling to simulate the natural convection of Al₂O₃/water nanofluid in a square enclosure using the thermal lattice Boltzmann flux solver (TLBFS), which is a coupling method combining the finite volume method to discretize the macroscopic governing equations in space, and reconstructed flux solutions at the interface between two adjacent cell centers by using the single-relaxation-time Lattice Boltzmann model. The top mpotivating priority of this paper was to establish a simple and effective numerical calculation method to solve natural convection problems. Therefore, it was necessary to introduce the boundary treatment technique in the solver. Tong et al. [35] applied the multiblock lattice Boltzmann method with a fixed Eulerian mesh, and the fouling layer was represented by an immersed boundary with Lagrangian points. The shape change of the fouling layer could be carried out by deforming the immersed boundary, while keeping the mesh of flow simulation unchanged. Suzuki et al. [36] simulated lift and thrust generation by a butterfly-like flapping wing body model by means of immersed boundary lattice Boltzmann simulations. The immersed boundary method is an effective and simple method to treat solid surface boundary conditions and the numerical method based on a non-body-fitted grid can avoid the abundant work involved in grid generation. Therefore, the immersed boundary method was applied to implement the no-slip boundary condition and Dirichlet boundary condition was applied for natural convection around a bluff body in a square enclosure with the purpose of effective treatment of surface boundaries. Natural convection problems were investigated at different Rayleigh numbers and nanoparticle volume fractions. Influences of the Rayleigh number and nanoparticle volume fraction on the streamlines, isotherms and average Nusselt number were studied.

2. Governing Equations and Numerical Method

2.1. The Macroscopic Governing Equations

For incompressible thermal nanofluid, in consideration of single phase and constant properties flow conditions, the macroscopic governing equations of natural convection in a two-dimensional enclosure can be written as follows:

Continuity equation;

$$\frac{\partial {\rho }_{nf}}{\partial t}+\nabla \cdot {\rho }_{nf}\mathit{u}=0$$

(1)

Momentum equation;

$$\frac{\partial }{\partial t}\left({\rho }_{nf}\mathit{u}\right)+\nabla \left({\rho }_{nf}\mathit{u}\mathit{u}\right)=-\nabla p+{\mu }_{nf}\nabla \left[\left(\nabla \mathit{u}+{\left(\nabla \mathit{u}\right)}^T\right)\right]+{\mathit{F}}_{nf}$$

(2)

Energy equation;

$$\frac{\partial }{\partial t}\left({\rho }_{nf}e\right)+\nabla \left({\rho }_{nf}\mathit{u}e\right)={\chi }_{nf}{\nabla }^2\left({\rho }_{nf}e\right)$$

(3)

where ρ, u, p and m represent fluid density, velocity, pressure, dynamic viscosity coefficient, respectively; e stands for internal energy defined as e = DRT/2, where D is the dimension, R is the gas constant and T represents the temperature; χ is the thermal diffusivity. The subscript nf denotes the nanofluid.

Natural convection heat transfer in nanofluids is studied in a two-dimensional enclosure. Nanoparticles considered to be spherical and frictional forces are neglected. The flow is assumed as laminar with a single-phase homogeneous mixture. The buoyancy force always plays an essential role as an external force. Using the Boussinesq approximation, the force source term can be defined as:

$${\mathit{F}}_{nf}={\rho }_{nf}{\beta }_{nf}g\left(T-T_m\right)\mathit{j}$$

(4)

where g represents the gravity acceleration, β is the thermal expansion coefficient and T_m is the average temperature.

According to Chapman-Enskog analysis, the relationships can be established between the fluxes and the distribution functions of the lattice Boltzmann model. Based on the thermal lattice Boltzmann flux solver (TLBFS), the governing Equations (1)–(3) can be rewritten as:

$$\frac{\partial {\rho }_{nf}}{\partial t}+\nabla \cdot \left({\displaystyle \sum _{\alpha }{e}_{\alpha }{f}_{\alpha }^{eq}}\right)=0$$

(5)

$$\frac{\partial {\rho }_{nf}\mathit{u}}{\partial t}+\nabla \cdot {\prod }_1={\mathit{F}}_{nf}$$

(6)

$$\frac{\partial {\rho }_{nf}e}{\partial t}+\nabla \cdot {\prod }_2=0$$

(7)

where

$${\prod }_1={\displaystyle \sum _{\alpha =0}^N{(e_{\alpha })}_{\beta }{(e_{\alpha })}_{\gamma }}\left[f_{\alpha }^{eq}+(I-\frac{1}{2{\tau }_v})f_{\alpha }^{neq}\right]$$

(8)

$${\prod }_2={\displaystyle \sum _{\alpha =0}^N{e}_{\alpha }}\left[g_{\alpha }^{eq}+(I-\frac{1}{2{\tau }_c})g_{\alpha }^{neq}\right]$$

(9)

$${\tau }_v={\mu }_{nf}/({\rho }_{nf}{c}_s^2{\delta }_t)+0.5$$

(10)

$${\tau }_c={\chi }_{nf}/\left(2c_s{}^2{\delta }_t\right)+0.5$$

(11)

From the above process, the macroscopic flow variables and fluxes can be computed by equilibrium and non-equilibrium distribution functions of the lattice model for the governing equations of nanofluid. Equations (8) and (9) are used to solve the macroscopic flow variables, and fluxes can be evaluated by the thermal lattice Boltzmann flux solver, which is introduced in detail in the next section. The force source term is added at the cell center during the calculation process.

2.2. Thermal Lattice Boltzmann Flux Solver

The discrete term of the governing Equations (5)–(7) by finite volume method:

$$\frac{d{\mathit{W}}_i}{dt}=\frac{1}{\Delta V_i}{\displaystyle \sum _k{\mathit{R}}_k}d{S}_k+\mathit{F}$$

(12)

where $\mathit{W}={[{\rho }_{nf},{\rho }_{nf}u,{\rho }_{nf}v,{\rho }_{nf}e]}^{\mathrm{T}}$; dV_i and dS_k are the volume of ith control volume and the area of the kth interface. For the 2D case, the D2Q9 lattice velocity model [37] is used for momentum and energy fluxes. The expression of the fluxes R_k at the cell interfaces is as followed:

$${\mathit{R}}_k=(\begin{array}{c}{n}_x\left(f_1^{eq}-f_3^{eq}+f_5^{eq}-f_6^{eq}-f_7^{eq}+f_8^{eq}\right)+n_y\left(f_2^{eq}-f_4^{eq}+f_5^{eq}+f_6^{eq}-f_7^{eq}-f_8^{eq}\right)\\ n_x\left(f_1^{^}+f_3^{^}+f_5^{^}+f_6^{^}+f_7^{^}+f_8^{^}\right)+n_y\left(f_5^{^}-f_6^{^}+f_7^{^}-f_8^{^}\right)\\ n_x\left(f_5^{^}-f_6^{^}+f_7^{^}-f_8^{^}\right)+n_y\left(f_2^{^}+f_4^{^}+f_5^{^}+f_6^{^}+f_7^{^}+f_8^{^}\right)\\ n_x\left(g_1^{^}-g_3^{^}+g_5^{^}-g_6^{^}-g_7^{^}+g_8^{^}\right)+n_y\left(g_2^{^}-g_4^{^}+g_5^{^}+g_6^{^}-g_7^{^}-g_8^{^}\right)\end{array})$$

(13)

$$f_{\alpha }^{^}=f_{\alpha }^{eq}+(1-\frac{1}{2{\tau }_v})f_{\alpha }^{neq}$$

(14)

$$g_{\alpha }^{^}=g_{\alpha }^{eq}+(1-\frac{1}{2{\tau }_c})g_{\alpha }^{neq}$$

(15)

From Equations (13)–(15), it can be seen that the important segment to solve fluxes is to accurately evaluate the $f_{\alpha }^{eq}$, $f_{\alpha }^{^}$ and $g_{\alpha }^{^}$ terms.

The simplified thermal lattice Boltzmann model with BGK approximation can be written as:

$$f_{\alpha }\left(\mathit{r}+{\mathit{e}}_{\alpha }{\delta }_t,t+{\delta }_t\right)-f_{\alpha }\left(\mathit{r},t\right)=-\frac{1}{{\tau }_v}\left[f_{\alpha }\left(\mathit{r},t\right)-f_{\alpha }^{eq}\left(\mathit{r},t\right)\right]$$

(16)

$$g_{\alpha }\left(\mathit{r}+{\mathit{e}}_{\alpha }{\delta }_t,t+{\delta }_t\right)-g_{\alpha }\left(\mathit{r},t\right)=-\frac{1}{{\tau }_c}\left[g_{\alpha }\left(\mathit{r},t\right)-g_{\alpha }^{eq}\left(\mathit{r},t\right)\right]$$

(17)

In which equilibrium density distribution function and equilibrium internal energy distribution function is given as:

$$f_{\alpha }^{eq}\left(\mathit{r},t\right)=\rho w_{\alpha }\left[1+\frac{{\mathit{e}}_{\alpha }\cdot \mathit{u}}{c_s^2}+\frac{{\left({\mathit{e}}_{\alpha }\cdot \mathit{u}\right)}^2-{\left(c_s\left|\mathit{u}\right|\right)}^2}{2c_s^4}\right]$$

(18)

$$g_{\alpha }^{eq}\left(\mathit{r},t\right)=\left\{\begin{array}{cc}-\frac{2\rho }{3}\frac{{\left|\mathit{u}\right|}^2}{c^2},& \alpha =0\\ \frac{\rho e}{9}\left[\frac{3}{2}+\frac{3}{2}\cdot \frac{{\mathit{e}}_{\alpha }\cdot \mathit{u}}{c^2}+\frac{9}{2}\cdot \frac{{\left({\mathit{e}}_{\alpha }\cdot \mathit{u}\right)}^2}{c^4}-\frac{3}{2}\cdot \frac{{\left|\mathit{u}\right|}^2}{c^2}\right],& \alpha =1,2,3,4\\ \frac{\rho e}{36}\left[3+6\cdot \frac{{\mathit{e}}_{\alpha }\cdot \mathit{u}}{c^2}+\frac{9}{2}\cdot \frac{{\left({\mathit{e}}_{\alpha }\cdot \mathit{u}\right)}^2}{c^4}-\frac{3}{2}\cdot \frac{{\left|\mathit{u}\right|}^2}{c^2}\right],& \alpha =5,6,7,8\end{array}\right.$$

(19)

Using the second-order Taylor series expansion, Equations (16) and (17) can be transformed as below:

$${\delta }_t\left(\frac{\partial }{\partial t}+{\mathit{e}}_{\alpha }\cdot \nabla \right)f_{\alpha }+\frac{{\delta }_t^2}{2}{\left(\frac{\partial }{\partial t}+{\mathit{e}}_{\alpha }\cdot \nabla \right)}^2{f}_{\alpha }+\frac{1}{\tau }\left(f_{\alpha }-f_{\alpha }^{eq}\right)+O\left({\delta }_t^3\right)=0$$

(20)

$${\delta }_t\left(\frac{\partial }{\partial t}+{\mathit{e}}_{\alpha }\cdot \nabla \right)g_{\alpha }+\frac{{\delta }_t^2}{2}{\left(\frac{\partial }{\partial t}+{\mathit{e}}_{\alpha }\cdot \nabla \right)}^2{g}_{\alpha }+\frac{1}{\tau }\left(g_{\alpha }-g_{\alpha }^{eq}\right)+O\left({\delta }_t^3\right)=0$$

(21)

By the multi-scale Chapman-Enskog expansion, the distribution function, the temporal and spatial derivatives, the non-equilibrium density and energy distribution functions can be transformed into an expression only related to the equilibrium distribution functions and can be derived from:

$$f_{\alpha }^{neq}\left(\mathit{r},t\right)=-{\tau }_v\left[f_{\alpha }^{eq}\left(\mathit{r},t\right)-f_{\alpha }^{eq}\left(\mathit{r}-{\mathit{e}}_{\alpha }{\delta }_t,t-{\delta }_t\right)\right]$$

(22)

$$g_{\alpha }^{neq}\left(\mathit{r},t\right)=-{\tau }_c\left[g_{\alpha }^{eq}\left(\mathit{r},t\right)-g_{\alpha }^{eq}\left(\mathit{r}-{\mathit{e}}_{\alpha }{\delta }_t,t-{\delta }_t\right)\right]$$

(23)

From Figure 1, the flow properties of eight vertices of the D2Q9 model can be evaluated by interpolation with the given flow properties at the cell centers of two adjacent control volumes. The values r_i, r_i+1 and r are defined as the physical positions of the two cell centers and their interfaces, respectively. The interpolation formulation can be given as:

$$\psi \left(\mathit{r}-{\mathit{e}}_{\alpha }{\delta }_t,t-{\delta }_t\right)=\left\{\begin{array}{c}\psi \left({\mathit{r}}_i\right)+\left(\mathit{r}-{\mathit{e}}_{\alpha }{\delta }_t-{\mathit{r}}_i\right)\cdot \nabla \psi \left({\mathit{r}}_i\right)\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\mathit{r}-{\mathit{e}}_{\alpha }{\delta }_t in {\mathsf{\Omega }}_i\\ \psi \left({\mathit{r}}_{i+1}\right)+\left(\mathit{r}-{\mathit{e}}_{\alpha }{\delta }_t-{\mathit{r}}_{i+1}\right)\cdot \nabla \psi \left({\mathit{r}}_{i+1}\right)\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\mathit{r}-{\mathit{e}}_{\alpha }{\delta }_t in {\mathsf{\Omega }}_{i+1}\end{array}\right.$$

(24)

where $\psi$ stands for the flow properties, including ρ, u, v and e. $f_{\alpha }^{eq}\left(\mathit{r}-{\mathit{e}}_{\alpha }{\delta }_t,t-{\delta }_t\right)$ and $g_{\alpha }^{eq}\left(\mathit{r}-{\mathit{e}}_{\alpha }{\delta }_t,t-{\delta }_t\right)$ can be obtained by the corresponding equilibrium density distribution function and energy distribution function. Then, the flow properties of the cell interface can be written as:

$$\rho \left(\mathit{r},t\right)={\displaystyle \sum _{\alpha =0}{f}_{\alpha }^{eq}\left(\mathit{r}-{\mathit{e}}_{\alpha }{\delta }_t,t-{\delta }_t\right)}$$

(25)

$$\rho \left(\mathit{r},t\right)u\left(\mathit{r},t\right)={\displaystyle \sum _{\alpha =0}{\mathit{e}}_{\alpha }{f}_{\alpha }^{eq}\left(\mathit{r}-{\mathit{e}}_{\alpha }{\delta }_t,t-{\delta }_t\right)}$$

(26)

$$\rho \left(\mathit{r},t\right)e\left(\mathit{r},t\right)={\displaystyle \sum _{\alpha =0}{g}_{\alpha }^{eq}\left(\mathit{r}-{\mathit{e}}_{\alpha }{\delta }_t,t-{\delta }_t\right)}$$

(27)

Next, $f_{\alpha }^{eq}\left(\mathit{r},t\right)$ and $g_{\alpha }^{eq}\left(\mathit{r},t\right)$ can also be easily solved by distribution functions. After obtaining the equilibrium distribution functions, the fluxes can be evaluated according to Equation (13).

2.3. Computational Sequence

The complete numerical simulation procedures for each time step of the proposed method are summarized below.

• According to the fluid properties of the nanofluid, determine initial velocity and temperature field;
• Based on the grid size, identify a streaming time step at each interface and then the single relaxation parameters, including dynamic viscosity and the thermal diffusivity;
• Apply the D2Q9 model to compute the density and energy equilibrium distribution functions $f_{\alpha }^{eq}\left(\mathit{r}-{\mathit{e}}_{\alpha }{\delta }_t,t-{\delta }_t\right)$ and $g_{\alpha }^{eq}\left(\mathit{r}-{\mathit{e}}_{\alpha }{\delta }_t,t-{\delta }_t\right)$ around the middle point r of each interface;
• Compute the macroscopic flow properties of nanofluid at the cell interface and then compute $f_{\alpha }^{eq}\left(\mathit{r},t\right)$ and $g_{\alpha }^{eq}\left(\mathit{r},t\right)$ by the equilibrium distribution functions of the D2Q9 model;
• Compute $f_{\alpha }^{^}$ and $g_{\alpha }^{^}$ terms, then the fluxes at the cell interface can be solved by Equation (13);
• Calculate the force source term and add this term to the fluxes;
• Solve Equations (5)–(7) to obtain the macroscopic flow properties of the nanofluid;
• Repeat steps (3)–(7) until the following convergence criterion is satisfied.

3. Numerical Examples of Natural Convection in a Square Enclosure

3.1. Problem Description

The computational domain and boundary conditions are shown in Figure 2. From this figure, it can be seen that the no slip boundary condition was applied on four walls. The adiabatic condition was set on the top and bottom walls and temperatures of 1 and 0 were applied on the left and right walls, respectively. The non-dimensional parameters, Prandtl number Pr and Rayleigh number Ra, were applied to determine the dynamic similarity as follows:

$$Pr=\nu /\chi$$

(28)

$$Ra=\frac{V_c^2\cdot L^2}{\nu \cdot \chi }$$

(29)

where L = 1 is the characteristic length of the square cavity and V_c is the characteristic thermal velocity which is constrained by the low Mach number limit. In the present simulations, V_c = 0.1 was set in order to ensure incompressible viscous flow.

In the present study, Al₂O₃/water nanofluid was used. The thermophysical properties of the water and nanoparticles are listed in

Table 1. Thermophysical properties of fluid and nanoparticles.

Properties	Fluid Phase (Water)	Solid Phase (Al2O3)
ρ (kg/m3)	997	3880
cp (J/kg·K)	4179	765
β (1/K)	0.00021	0.0000085
k (W/m·K)	0.613	40
μ (kg/m·s)	0.000855	-

. The homogeneous model for nanofluid was adopted. Physical properties of the nanofluids, including density, specific heat and thermal expansion coefficient, were obtained using the classical formula developed for conventional solid–liquid mixtures as follows:

$${\rho }_{nf}=\left(1-\varphi \right){\rho }_f+\varphi {\rho }_s$$

(30)

$${\left(\rho c_p\right)}_{nf}=\left(1-\varphi \right){\left(\rho c_p\right)}_f+\varphi {\left(\rho c_p\right)}_s$$

(31)

$${\beta }_{nf}=\left(1-\varphi \right){\beta }_f+\varphi {\beta }_s$$

(32)

where ϕ refers to the volume concentration of nanoparticles and the subscripts s, f denote the particle and base fluids.

The effective viscosity and thermal conductivity of the nanofluid strongly affect the heat transfer rate and flow characteristics of nanofluids. The effective viscosity could be estimated by experimental correlation for 47 nm Al₂O₃/water nanofluid by Angue Mintsa et al. [38] and thermal conductivity was given by Gherasim et al. [39] as follows:

$${\mu }_{nf}=0.904e^{14.8\varphi }\mu$$

(33)

$${\kappa }_{nf}=\left(1.72\varphi +1.0\right){\kappa }_f$$

(34)

In the present simulations, the convergence criterion for flow field and temperature field were respectively given as follows:

$$Error1=\frac{{\displaystyle {\sum }_{ij}\left|\sqrt{{\left(u\left(i,j,t+{\delta }_t\right)\right)}^2+{\left(v\left(i,j,t+{\delta }_t\right)\right)}^2}-\sqrt{{\left(u\left(i,j,t\right)\right)}^2+{\left(v\left(i,j,t\right)\right)}^2}\right|}}{{\displaystyle {\sum }_{ij}\sqrt{{\left(u\left(i,j,t+{\delta }_t\right)\right)}^2+{\left(v\left(i,j,t+{\delta }_t\right)\right)}^2}}}\le 1\times {10}^{-7}$$

(35)

$$Error2=\frac{{\displaystyle {\sum }_{ij}\left|T\left(i,j,t+{\delta }_t\right)-T\left(i,j,t\right)\right|}}{{\displaystyle {\sum }_{ij}T\left(i,j,t+{\delta }_t\right)}}\le 1\times {10}^{-7}$$

(36)

3.2. Natural Convection of Pure Fluid in a Square Enclosure

To testify as to the accuracy and performance of the lattice Boltzmann flux solver based on the population model, the classical natural convection in a square enclosure filled with air and water was studied at Ra = 10³, 10⁴, 10⁵ and 10⁶.

Firstly, a grid independent study was conducted on five different uniform grids of 101 × 101, 151 × 151, 201 × 201, 251 × 251 and 301 × 301 for the natural convection problem at Ra = 10⁶ and Pr = 0.7. As shown in

Table 2. Grid independent study on uniform of natural convection at Ra = 10⁶.

Method	Grids	Nuavg	y	umax	x	vmax
Present	101 × 101	8.788	0.855	64.22	0.0350	217.60
	151 × 151	8.809	0.850	64.66	0.0367	219.32
	201 × 201	8.816	0.853	64.99	0.0375	220.05
	251 × 251	8.819	0.854	66.18	0.0380	220.17
	301 × 301	8.819	0.852	66.99	0.0383	220.14
De Vahl Davis [40]		8.800	0.850	64.63	0.0397	219.36
Hortmann et al. [41]		8.825	-	-	-	-

, when the mesh size was 201 × 201, or even larger, the average Nusselt number did not change much and the value was between the benchmark solutions of Davis [40] and Hortmann et al. [41]. When the mesh size was larger than 151 × 151, the maximum horizontal velocity on the vertical mid-plane, the maximum vertical velocity on the horizontal mid-plane and their locations were in agreement with the benchmark solutions of Davis [40]. The above results illustrated grid independence on uniform grids of 201 × 201, for the case of Ra = 10⁶.

Based on the above results, the grid independent study was conducted on non-uniform grids by using the size of less than 201 × 201.

Table 3. Grid independent study on non-uniform of natural convection at Ra = 10⁶.

Grids	Nuavg	y	umax	x	vmax
81 × 81	8.803	0.856	64.77	0.0358	218.98
101 × 101	8.811	0.855	64.99	0.0349	219.16
121 × 121	8.816	0.854	65.04	0.0400	219.37
141 × 141	8.818	0.853	65.03	0.0388	220.15
161 × 161	8.819	0.853	65.03	0.0379	220.14
181 × 181	8.820	0.853	65.06	0.0372	220.22

shows the numerical results of six different non-uniform grids of natural convection at Ra = 10⁶. From this table, the results were close to the data of uniform grids of 201 × 201 when the non-uniform mesh was more than 121 × 121. In order to ensure the accuracy and efficiency of numerical simulations, the non-uniform grid of 141 × 141 was chosen to simulate natural convection in a square enclosure.

The average Nusselt number results at different Rayleigh numbers are listed in

Table 4. Comparison of average Nusselt numbers at different Rayleigh numbers.

Ra	Air				Water
	Present	Davis [40]	Khanafer et al. [6]	Qi [9]	Present	Kahveci [42]	Lai and Yang [1]
103	1.118	1.118	1.118	1.118	1.119	-	1.128
104	2.246	2.243	2.245	2.247	2.278	2.274	2.286
105	4.522	4.519	4.522	4.522	4.725	4.722	4.729
106	8.818	8.800	8.826	8.808	9.204	9.230	9.173

, and it can be seen that the numerical simulation results were in good agreement with previous literature results at different Rayleigh numbers. This illustrated the accuracy of the present method for natural convection.

Figure 3 shows the temperature distribution at horizontal midsections of the enclosure. For the enclosure filled with air, the results of Ra = 10⁵ were compared with the numerical results of Khanafer et al. [6] and the experimental results of Krane and Jessee [43]. For the enclosure filled with water, the results were compared with numerical results of Lai and Yang [1]. It was noted from the comparisons that the solutions were in excellent agreement. This illustrated that the method in this paper could capture the temperature field very well.

The streamlines and isotherms of air and water at various Rayleigh numbers are shown in Figure 4 and Figure 5, respectively. It can be seen that the natural convection and heat transfer between the wall and fluid were enhanced as Ra increased. For Ra ≤ 10⁴, the flow characteristic was to appear as a central vortex. For Ra > 10⁴, the central vortex became more expanded and finally broke up into two vortices so that temperature boundary layers were formed. The above phenomenon agreed well with previous studies.

3.3. Natural Convection of Nanofluid in a Square Enclosure

After validating the numerical method for natural convection in a square enclosure filled with pure fluid, the natural convection in a square enclosure filled with Al₂O₃-water nanofluid of nanoparticles having volume fraction ϕ = 1–4% at Ra = 10³–10⁶ was simulated to validate the present numerical algorithm. The presented averaged Nusselt numbers were compared with the numerical results of Lai and Yang [1] and listed in

Table 5. Comparison of average Nusselt numbers with the previous studies.

Ra	ϕ	Present	Ref [1]	Relative Error (%)
103	0.01	1.139	1.147	0.697
	0.02	1.158	1.167	0.771
	0.03	1.177	1.186	0.756
	0.04	1.196	1.206	0.829
104	0.01	2.317	2.326	0.387
	0.02	2.357	2.366	0.380
	0.03	2.396	2.406	0.416
	0.04	2.435	2.445	0.409
105	0.01	4.807	4.811	0.008
	0.02	4.890	4.894	0.008
	0.03	4.972	4.977	0.010
	0.04	5.054	5.059	0.010
106	0.01	9.366	9.331	0.375
	0.02	9.528	9.492	0.386
	0.03	9.688	9.653	0.363
	0.04	9.849	9.813	0.367

. It shows that there was a good agreement and the relative errors were less than 0.8%, which further illustrated that the present numerical method could simulate the natural convection of nanofluid at different Rayleigh numbers and nanoparticle volume fractions.

In the present numerical simulations, the effect of nanoparticle suspensions (Al₂O₃-water) on flow and temperature characteristics for Ra = 10³–10⁶ and nanoparticles volume fraction ϕ = 0–10% were studied. The variation of average Nusselt number against solid volume fraction for different Rayleigh numbers is shown in Figure 6a and the variation of average Nusselt number against Rayleigh number for different solid volume fractions is shown in Figure 6b.

Numerical results indicated that average Nusselt number increased with the increase of Ra and ϕ. This illustrated that the function of heat transfer was enhanced with the augmentation of nanofluid thermal conductivity, which indicated that the major mechanism of heat transfer in flowing fluid was thermal dispersion. At the same Ra, the relationship of the average Nusselt number and solid volume fraction was almost linear. At the same solid volume fraction, the relationship of the average Nusselt number and Ra presented an exponential form. At higher Rayleigh number, the greater the heat transfer rate that could be obtained.

Figure 7 and Figure 8 indicate the isotherms and streamlines of nanofluid (Al₂O₃-water) at Ra = 10³–10⁶ and ϕ = 0%, 5% and 10%, which show the effect of volume fraction and Ra on flow field and temperature field very well. From Figure 7, it can be seen that heat transfer between the wall and fluid were enhanced as Ra increased. As the volume fraction of nanoparticles increased, the isotherm changed slightly. That was because the mixture flow became more viscous, due to the nanoparticles. The velocity of flow fluid reduced and then natural convection weakened. However, the function of heat transfer in total computational domain was enhanced, which was attributed to the augmentation of nanofluid thermal conductivity.

From Figure 8, it can be observed that the flow appeared as a central vortex for lower Ra. As Ra increased, the central vortex became more expanded and finally broke up into two vortices, so that temperature boundary layers were formed. For pure fluid, the vortex formed in the enclosure as a result of the buoyancy effect. By increasing the volume fraction of nanoparticles, the intensity of streamlines increased, due to the high energy transport through the flow as a result of irregular motion of the ultra-fine particles.

4. Numerical Examples of Natural Convection around Bluff Body in a Square Enclosure

4.1. Problem Description

The boundary condition-enforced immersed boundary method was chosen for treatment of the solid boundary conditions. Based on the immersed boundary method and thermal lattice Boltzmann flux solver (IB-TLBFS), the macroscopic governing equations can be rewritten as:

$$\frac{\partial {\rho }_{nf}}{\partial t}+\nabla \cdot \left({\displaystyle \sum _{\alpha }{e}_{\alpha }{f}_{\alpha }^{eq}}\right)=0$$

(37)

$$\frac{\partial {\rho }_{nf}\mathit{u}}{\partial t}+\nabla \cdot {\prod }_1={\mathit{F}}_{nf}+{\mathit{f}}_b$$

(38)

$$\frac{\partial {\rho }_{nf}e}{\partial t}+\nabla \cdot {\prod }_2=q_b$$

(39)

where the force source term ${\mathit{f}}_b$ and the heat source term $q_b$ are both generated by the immersed boundary. To solve the governing equations, the calculation process is divided into two steps: the first step predicts the state variables without taking account of the boundary function and the second step corrects velocity and temperature by the immersed boundary method.

In this work, the implicit velocity correction scheme proposed by Wang et al. [44] was be applied in view of satisfaction of the no slip boundary. The implicit heat source scheme proposed by Ren et al. [45] was applied for the Dirichlet boundary conditions of the temperature field.

Natural convection of a heated bluff body in a square enclosure was studied. The physical models, computational domain and boundary conditions are presented in Figure 9. All boundaries were no-slip and isothermal boundary conditions. The flow was assumed to be laminar and driven by the temperature difference.

Numerical investigations were carried on two types of bluff bodies, a circular cylinder and a square. The four side walls of the outer square enclosure were cooled isothermally at T_C and the side length was L. The wall of the inner bluff body was heated isothermally at T_H and D and a represent the diameter of the circular cylinder and the side length of the square, respectively. For fixed Rayleigh number, numerical simulation cases were designed to have a fixed perimeter for different bluff bodies and the influences of geometry on the heat transfer is discussed in detail.

4.2. Natural Convection in the Annulus between Concentric Circular Cylinder and Square Enclosure

After validating the numerical algorithm of the thermal lattice Boltzmann flux solver, natural convection in the annulus between concentric circular cylinder and square enclosure at Ra = 10⁴, 10⁵ and 10⁶ were simulated to validate the immersed boundary method and code implementation. Numerical simulations were conducted for three different aspect ratios (Ar = 1.67, 2.5 and 5.0). The average Nusselt number was also computed and compared with reference data in the literature.

The computed average Nusselt numbers are compared in

Table 6. Comparison of average Nusselt numbers at different Ra and Ar.

Ra	Ar	Present	Ren et al. [45]	Shu and Zhu [46]	Moukalled and Acharya [47]
104	1.67	5.425	5.303	5.40	5.826
	2.50	3.256	3.161	3.24	3.331
	5.00	2.090	2.051	2.08	2.071
105	1.67	6.285	6.171	6.21	6.212
	2.50	4.954	4.836	4.86	5.080
	5.00	3.809	3.704	3.79	3.825
106	1.67	11.943	11.857	12.00	11.620
	2.50	9.002	8.546	8.90	9.374
	5.00	6.110	5.944	6.11	6.107

with those of Ren et al. [45], Shu et al. [46] and Moukalled et al. [47]. From this table, it can be seen that the present results of the method combining IBM and TLBFS agreed very well with reference data. Besides this, the results revealed that the average Nusselt number greatly depended on Rayleigh number and aspect ratio. Due to buoyancy-induced convection, the average Nusselt number increased with increase of Ra, while it decreased with increase of Ar, due to the effect of annulus gap space.

The streamlines and isotherms in the annulus at various Rayleigh numbers and aspect ratios are shown in Figure 10. Conduction dominated the flow field and a relatively weak convective flow could be observed in the annulus at the lower Ra. As the Rayleigh number increased, the strength of the convective flow grew and the center of the recirculation eddy changed its position. When Ra = 10⁶, a relatively stronger convective flow dominated the fluid field and a higher temperature gradient could be observed. In contrast, stronger convective flow and higher temperature gradient could be observed in the case of lower values of Ar.

4.3. Natural Convection of Nanofluid between Bluff Body and Square Enclosure

In the present study, numerical investigations of natural convection between heated bluff body and square enclosure were conducted for nanoparticles having volume fractions of ϕ = 0%, 2% and 4% and Rayleigh numbers of Ra = 10⁴, 10⁵ and 10⁶.

The averaged Nusselt numbers are listed in

Table 7. Comparison of average Nusselt numbers at various Ra and ϕ.

	Cylinder			Square
ϕ	Ra = 104	Ra = 105	Ra = 106	Ra = 104	Ra = 105	Ra = 106
0.00	3.131	5.080	9.144	2.9432	4.8675	8.6541
0.02	3.132	5.098	9.200	2.9447	4.8884	8.7186
0.04	3.134	5.117	9.263	2.9463	4.9103	8.7883

. The numerical simulation results indicated that average Nusselt number increased with increase of Ra and ϕ, which was the same as occurred in natural convection of square enclosure. By comparison, the averaged Nusselt number of the natural convection around a circular cylinder in an enclosure was greater than that of the square at the same calculation conditions. This illustrated that a smooth geometrical shape was beneficial to heat transfer.

Figure 11 and Figure 12 present the distribution of the isotherms for different Rayleigh numbers (Ra = 10⁵ and 10⁶) and values of nanoparticle volume fractions (ϕ = 0 and 0.04). An overview of this figure indicated that the thermal fields strongly depended on Rayleigh number. When Ra = 10⁵ or even lower, the isotherms of ϕ = 0 were almost close to that of ϕ = 0.04, which illustrated that nanoparticle volume fraction played a smaller role in heat transfer and flow pattern. When Ra = 10⁶, there were significant differences between the isotherms of ϕ = 0 and ϕ = 0.04, which illustrated that nanoparticle volume fraction played a role in heat transfer and flow pattern for high Ra. The thickness of the thermal boundary layer decreased as the volume fraction increased, which was due to the increasing conduction heat transfer by adding nanoparticle volume fraction.

Figure 13 shows the streamlines for natural convection around a circular cylinder and square at nanoparticle volume fractions ϕ = 0.04 and Ra = 10⁶. From Table 7, it can be seen that the preferable heat transfer effect could be acquired by the cylinder in comparison with the square at the same perimeter. That was because the velocity and temperature gradients around the sharp corners of the square dramatically changed, which prevented the heat transfer effect.

5. Conclusions

The thermal lattice Boltzmann flux solver (TLBFS) was applied to simulate natural convection of nanofluid in a square enclosure. This method couples the finite volume method and lattice Boltzmann models to realize the solution of incompressible thermal flow. To validate the accuracy and performance of this method, natural convection in a square enclosure filled with pure fluid (air and water) was first studied. There were good agreements with previous literature. Numerical investigations of fluid flow and convective heat transfer were performed. The effects of some parameters, such as the Rayleigh number (Ra), and volume fraction of nanoparticles (ϕ), on natural convection were analyzed. With increase in the Rayleigh number and nanoparticle volume fraction, the heat transfer rate increased and the nanofluid flow became more viscous and this led to a decrease in nanofluid motion velocity. The average Nusselt number was an increasing exponential function of the Rayleigh number and an increasing linear function of the nanoparticle volume fraction. Then, natural convection around a bluff body in a square enclosure was studied by a method combining TLBFS and immersed boundary method. Natural convection problems in the annulus between concentric circular cylinder and square enclosure without nanofluid were simulated, which validated the feasibility of the numerical algorithm and code implement. Numerical investigations of natural convection between heated bluff body (cylinder and square) and square enclosure were conducted for different nanoparticle volume fractions and Rayleigh numbers. The numerical results illustrated that heat transfer effect increased with increase of Ra and ϕ. At lower Ra, the function of heat transfer with the augmentation of nanofluid thermal conductivity was counteracted by the more viscous flow. Nevertheless, nanoparticles played a better role in enhancing natural convection at higher Ra. The above results declare that the TLBFS is a promising method for heat transfer of nanofluids of the future.