Numerical simulation of natural convection in a square enclosure filled with nanofluid using the two-phase Lattice Boltzmann method

Considering interaction forces (gravity and buoyancy force, drag force, interaction potential force, and Brownian force) between nanoparticles and a base fluid, a two-phase Lattice Boltzmann model for natural convection of nanofluid is developed in this work. It is applied to investigate the natural convection in a square enclosure (the left wall is kept at a high constant temperature (TH), and the top wall is kept at a low constant temperature (TC)) filled with Al2O3/H2O nanofluid. This model is validated by comparing numerical results with published results, and a satisfactory agreement is shown between them. The effects of different nanoparticle fractions and Rayleigh numbers on natural convection heat transfer of nanofluid are investigated. It is found that the average Nusselt number of the enclosure increases with increasing nanoparticle volume fraction and increases more rapidly at a high Rayleigh number. Also, the effects of forces on nanoparticle volume fraction distribution in the square enclosure are studied in this paper. It is found that the driving force of the temperature difference has the biggest effect on nanoparticle volume fraction distribution. In addition, the effects of interaction forces on flow and heat transfer are investigated. It is found that Brownian force, interaction potential force, and gravity-buoyancy force have positive effects on the enhancement of natural convective heat transfer, while drag force has a negative effect.


Background
Compared with common fluids such as water, nanofluid, using nanoscale particles dispersed in a base fluid, has an effect of enhancing the performance of natural convection heat transfer due to its high heat conductivity coefficient. Many researchers investigated nanoparticles and nanofluid in recent years. Wang et al. [1] synthesized stimuli-responsive magnetic nanoparticles and investigated the effect of nanoparticle fraction on its cleavage efficiency. Bora and Deb [2] developed a novel bioconjugate of stearic acid-capped maghemite nanoparticle (γ-Fe 2 O 3 ) with bovine serum albumin. Guo et al. [3] produced magnetic nanofluids containing γ-Fe 2 O 3 nanoparticles using a two-step method, measured their thermal conductivities and viscosity, and tested their convective heat transfer coefficients. Pinilla et al. [4] investigated the growth of Cu nanoparticles in a plasmaenhanced sputtering gas aggregation-type growth region. Yang and Liu [5] produced a kind of stable nanofluid by surface functionalization of silica nanoparticles. Zhu et al. [6] developed a wet chemical method to produce stable CuO nanofluids. Nadeem and Lee [7] investigated the steady boundary layer flow of nanofluid over an exponential stretching surface. Wang and Fan [8] reviewed the nanofluid research in the last 10 years.
Natural convection is applied in many fields, and extensive researches have been performed. Oztop et al. [9] and Ho et al. [10] respectively investigated natural convection in partially heated rectangular enclosures and discussed the effects of viscosity and thermal conductivity of nanofluid on laminar natural convection heat transfer in a square enclosure by a finite-volume method. Saleh et al. [11] investigated heat transfer enhancement utilizing nanofluids in a trapezoidal enclosure by a finite difference approach. Ghasemi et al. [12], Santra et al. [13], and Aminossadati et al. [14] numerically simulated natural convection in a triangular enclosure and studied the behavior of natural convection heat transfer in a differentially heated square cavity, described a study on natural convection of a heat source embedded in the bottom wall of an enclosure, and used the SIMPLE algorithm to solve the governing equation. Kargar et al. [15] used computational fluid dynamics and an artificial neural network to investigate the cooling performance of two electronic components in an enclosure. Abu-Nada et al. [16] investigated the effect of variable properties on natural convection in enclosures filled with nanofluid, and the governing equations are solved by an efficient finite-volume method. Hwang et al. [17] investigated the thermal characteristics of natural convection in a rectangular cavity heated from below by Jang and Choi's model [18].
The Lattice Boltzmann method is a new way to investigate natural convection. Compared with the above traditional methods, the Lattice Boltzmann method has many merits including that boundary conditions can be conveniently dealt with, the transform between macroscopic and microscopic equations is easily achieved, the details of the fluid can be presented, and so on. In addition, nanofluid as the media can enhance heat transfer due to factors such as nanofluids having higher thermal conductivity and the nanoparticles in the fluid disturbing the laminar flow. Therefore, many researchers undertook investigations on the natural convection of nanofluids by the Lattice Boltzmann method. Barrios et al. [19] developed a Lattice Boltzmann model and applied it to investigate the natural convection of an enclosure with a partially heated left wall. Peng et al. [20] presented a simple a Lattice Boltzmann model without considering thermal diffusion, and this model is easily applied because it does not contain a gradient term. He et al. [21] proposed a new Lattice Boltzmann model which introduced an internal energy distribution function to simulate the temperature field, and the result has a good agreement with the benchmark solution. Nemati et al. [22] simulated the natural convection of a liddriven flow filled with Cu-water, CuO-water, and Al 2 O 3 -water nanofluids and discussed the effects of nanoparticle volume fraction and Reynolds number on the heat transfer. Wang et al. [23] presented a Lattice Boltzmann algorithm to simulate the heat transfer of a fluid-solid fluid, and the result has a satisfactory agreement with the published data. Dixit et al. [24] applied the Lattice Boltzmann method to investigate the natural convection of a square cavity at high Rayleigh numbers. Peng et al. [25] developed a 3D incompressible thermal Lattice Boltzmann model for natural convection in a cubic cavity. The above Lattice Boltzmann methods are all single-phase models, and the nanofluid was seen as a single-phase fluid without considering the interaction forces between nanoparticles and water. In addition, the effects of these interaction forces on heat transfer were disregarded.
There are few two-phase lattice Boltzmann models that consider the interaction forces between nanoparticles and a base fluid for natural convection in an enclosure. Xuan et al. [26] proposed a two-phase Lattice Boltzmann model to investigate sudden-start Couette flow and convection in parallel plate channels without researching the effect of forces on volume fraction distribution of nanoparticles. Because these forces were not investigated before our work, the effects of forces between water and nanoparticles on the fluid flow patterns were unknown. In addition, as we know, the nanoparticles in the fluid easily gather together and deposit, especially at high volume fraction. Hence, the nanoparticle distribution in the fluid flow is important for nanofluid application, which is another objective in our paper. However, the single-phase model cannot be used to investigate nanoparticle distribution. Furthermore, natural convection of a square enclosure (left wall kept at a high constant temperature (T H ), and top wall kept at a low constant temperature (T C )) filled with nanofluid is not investigated in the published literatures. In this paper, a two-phase Lattice Boltzmann model is proposed and applied to investigate the natural convection of a square enclosure (left wall kept at a high constant temperature (T H ), and top wall kept at a low constant temperature (T C )) filled with Al 2 O 3water nanofluid and the inhomogeneous distribution of nanoparticles in the square enclosure.

Lattice Boltzmann method
The density distribution function for a single-phase fluid is calculated as follows: where τ σ f is the dimensionless collision-relaxation time for the flow field, e α is the lattice velocity vector, the subscript α represents the lattice velocity direction, f σ α r; t ð Þ is the distribution function of the nanofluid with velocity e α (along the direction α) at lattice position r and time t, f σeq α r; t ð Þ is the local equilibrium distribution function, δ t is the time step, δ x is the lattice step, the order numbers α = 1,. . .,4 and α = 5,. . .,8, respectively represent the rectangular directions and the diagonal directions of the lattice, F σ 0 α is the external force term in the direction of the lattice velocity without interparticle interaction, G = − β(T nf − T 0 )g is the effective external force, where g is the gravity acceleration, β is the thermal expansion coefficient, T nf is the temperature of the nanofluid, and T 0 is the mean value of the high and low temperature of the walls.
A nanofluid is a two-phase fluid constituted by nanoparticles and a base fluid, and there are interaction forces (gravity and buoyancy force, drag force, interaction potential force, and Brownian force) between nanoparticles and the base fluid. Thus, the macroscopic density and velocity fields are simulated using the density distribution function by adding the forces term.
where F σ α represents the total interparticle interaction forces, and B α is one of the weight coefficients.
Because the total interparticle interaction forces cannot be optionally added in the lattice Boltzmann equation, we introduce an unknown coefficient in the total interparticle interaction forces. In order to enable the lattice Boltzmann equation including the total interparticle interaction forces to recover to the Navier-Stokes equation, based on the mass and momentum conservation, we used multi-scale technique to deduce the unknown coefficient which is equal to Due to the very long derivation process, we directly gave the final result in the paper. The weight coefficient B α is given as: For the two-dimensional nine-velocity LB model (D2Q9) considered herein, the discrete velocity set for each component α is: The density equilibrium distribution function is chosen as follows: where c 2 s ¼ c 2 3 is the lattice's sound velocity, and w α is the weight coefficient.
The macroscopic temperature field is simulated using the temperature distribution function.
where τ T is the dimensionless collision-relaxation time for the temperature field. The temperature equilibrium distribution function is chosen as follows: In the case of no internal forces and external forces, the macroscopic temperature, density and velocity are respectively calculated as follows: Considering the internal and external forces, the macroscopic velocities for nanoparticles and base fluid are modified to: where F p represents the total forces acting on the nanoparticles, F w represents the total forces acting on the base fluid, and L x L y represents the total number of lattices.
When the internal forces and external forces are considered, energy between nanoparticles and base fluid is exchanged, and the macroscopic temperature for nanoparticles and base fluid is then given as: where Φ αβ is the energy exchange between nanoparticles and base fluid, , and h αβ is the convective heat transfer coefficient of the nanofluid. The corresponding kinematic viscosity and thermal diffusion coefficients are respectively defined as follows: The dimensionless collision-relaxation times τ f and τ T are respectively given as follows: where Ma = 0.1, H = 1, c = 1, δt = 1, and the other parameters equations are given as follows: From Equations 18 and 19, the collision-relaxation time for the flow field and the temperature field can be calculated. For water phase, the τ f collision-relaxation times are respectively 0.51433 and 0.501433 at Ra = 10 3 and Ra = 10 5 , and the collision-relaxation time τ T is 0.5. For nanoparticle phase, the τ f collision-relaxation times are respectively 0.50096 and 0.500096 at Ra = 10 3 and Ra = 10 5 , and the collision-relaxation time τ T is 0.500025.

Interaction forces between base fluid and nanoparticles
As noted before, a nanofluid is, in reality, a kind of two-phase fluid. There are interaction forces between liquid and nanoparticles which affect the behavior of the nanofluid. The external forces include gravity and buoyancy forces F H , and the interparticle interaction forces include drag force (Stokes force) F D , interaction potential F A , and Brownian force F B . We introduce them as follows.    The gravity and buoyancy force is given as: where a is the radius of a nanoparticle, and Δρ ' is the mass density difference between the suspended nanoparticle and the base fluid. The drag force (Stokes force) is given as: where μ is the viscosity of the fluid, and Δu is the velocity difference between the nanoparticle and the base fluid.
The interaction potential is presented as [27]: where A is the Hamaker constant, and L cc is the centerto-center distance between particles.
The interaction potential force is shown as: where n i is the number of the particles within the adjacent lattice i, n i = ρ σ V/m σ , m σ is the mass of a single nanoparticle, and V is the volume of a single lattice. The Brownian force is calculated as [28]: where G i is a Gaussian random number with zero mean and unit variance, which is obtained from a program written by us, and C = 2γk B T = 2 × (6πηa) k B T, γ is the surface tension, k B is the Boltzmann constant, T is the absolute temperature, and η is the dynamic viscosity. The total per unit volume forces acting on nanoparticles of a single lattice is: where n is the number of the particles in the given lattice, and V is the lattice volume. In a nanofluid, the forces acting on the base fluid are mainly drag force and Brownian force. Thus the force acting on the base fluid in a given lattice is:

Results and discussion
The two-phase Lattice Boltzmann model is applied to simulate the natural convection heat transfer in a square cavity which is shown in Figure 1. The square cavity is filled with the Al 2 O 3 -water nanofluid. The thermophysical properties of water and Al 2 O 3 are given in Table 1. The height and the width of the enclosure are both H. The left wall is kept at a high constant temperature (T H ), and the top cold wall is kept at a low constant temperature (T C ). The boundary conditions of the other walls (right wall and bottom wall) are all adiabatic. The initial conditions for the four walls are given as follows: In the simulation, a non-equilibrium extrapolation scheme is adopted to deal with the boundary, and the criteria of the program convergence for the flow field and the temperature field are respectively given as follows: where ε is a small number, for example, for Ra = 1 × 10 3 , ε 1 = 10 −6 , and ε 2 = 10 −6 . About 2 weeks is needed to achieve the equilibrium state for the low Rayleigh number (Ra = 1 × 10 3 ), and about 1 month for the high Rayleigh number (Ra = 1 × 10 5 ). The Nusselt number can be expressed as: The heat transfer coefficient is computed from: The thermal conductivity of the nanofluid is defined by: Substituting Equations 33 and 34 into Equation 32, the local Nusselt number along the left wall can be written as: The average Nusselt number is determined from: In order to perform a grid independence test and validate the Lattice Boltzmann model proposed in this work, we used another square enclosure, because there are exact solutions for this square enclosure. The left wall is kept at a high constant temperature (T H ), and the right wall is kept at a low constant temperature (T C ). The boundary conditions of the other walls (top wall and bottom wall) are all adiabatic, and the other conditions are the same as those in Figure 1.
As shown in Table 2, the grid independence test is performed in a square enclosure using successively sized grids, 128 × 128, 192 × 192, 256 × 256, and 320 × 320 at Ra = 1 × 10 5 , Pr = 0.7. It can be seen from Table 2 that there is a bigger difference between the result obtained with grid sizes 128 × 128 and 192 × 192 and the result available from the literature [30] than when compared with the result obtained with grids 256 × 256 and 320 × 320. In addition, the result with grid 256 × 256 and the result with grid 320 × 320 are very close. In order to accelerate the numerical simulation, a grid size of 256 × 256 was chosen as a suitable one which can guarantee a grid-independent solution.
In order to validate the Lattice Boltzmann model proposed in this work, the temperature distribution at midsections of the enclosure at Ra = 1 × 10 5 , Pr = 0.7 is compared with the numerical results from Khanafer et al. [31] and experimental results from Krane et al. [32] in Figure 2. It can be seen that the results of this paper have a good agreement with those numerical [31]   and experimental [32] results. They are closer to the experimental [32] than the numerical [31] results. In addition, the Nusselt number results at different Rayleigh numbers of this paper are compared with other published literature listed in Table 3, and it can be seen that the results are in good agreement. Due to the temperature balance between water and nanoparticles, the temperature nephogram of water and nanoparticles for each of the nanoparticle fractions are identical. The temperature nephograms of nanofluid at Ra = 1 × 10 3 and Ra = 1 × 10 5 are presented in Figure 3. It can be seen that isotherms are more crooked with the higher Rayleigh number, which denotes that the heat transfer characteristic transforms from conduction to convection.
Because there are fewer nanoparticles than water molecules, and the drag force of nanoparticles on water is small, the velocity vectors of nanofluid with different nanoparticle fractions have such small differences that it is difficult to distinguish them. However, the differences can be observed in the Nusselt number distribution. For this reason, only the velocity vectors of nanofluid components with ϕ = 0.03 at different Rayleigh numbers are given as an example in Figure 4. Separating the nanofluid into its two constitutive components, it can be seen that the velocity vectors of the water component are larger than those of the nanoparticle component due to the law of conservation of momentum. The velocity difference between the water component and the nanoparticle component gives rise to the drag force. In addition, it can be seen that velocity increases with Rayleigh number, which can also explain that the heat transfer characteristic transforms from conduction to convection.
Driving force and interaction forces have a big effect on nanoparticle volume fraction distribution and the flow and heat transfer characteristics of the nanofluid. The main driving force in this work is the temperature difference. Interaction forces between nanoparticles and base fluid include gravity-buoyancy force, drag force, interaction potential force, and Brownian force. In order to compare the effects of these forces, the ranges of them are presented in Table 4. We used doubleprecision variables in our code. From Table 4, we can find that the temperature difference driving force F S is much bigger than the other forces (interaction forces between nanoparticles and base fluid). The driving force has the greatest effect on nanoparticle volume fraction distribution, and the effects of other forces on nanoparticle volume fraction distribution can be ignored in this case. However, these other forces play an important role in the flow and heat transfer of the nanofluid. Apart from the temperature difference driving force, the Brownian force is much larger than other forces, which is different from other two-phase fluids. For this reason, the Brownian force can enhance the heat transfer of the nanofluid by disturbing the flow boundary layer and the thermal boundary layer. Drag force comes about due to the velocity difference between nanoparticles and water molecules, and the nanoparticles in the water decrease the velocity of nanofluid in the enclosure, which in turn attenuates the natural convection of nanofluid in the enclosure. The interaction potential force prevents the nanoparticles from gathering together and keeps the nanoparticles dispersed in the water. In addition to the above forces, there is the gravity-buoyancy force, that is, the sum of gravity of the nanoparticles themselves and the buoyancy force of the water. The gravity-buoyancy force and temperature difference driving force together give rise to the velocity vectors of the nanofluid within the enclosure. In summary, Brownian force, interaction potential force, and gravity-buoyancy force contribute to the enhanced natural convective heat transfer, while drag force contributes to the attenuation of heat transfer.
The temperature difference driving force distribution in the square at different Rayleigh numbers is given in Figure 5. From Figure 5, we can see that the temperature difference driving force along the left wall (high temperature) and the top wall (low temperature) is high. Its direction along the high-temperature wall is upward, and that along the low-temperature wall is downward, while the temperature difference driving force in other regions far away from the two walls (left wall and top wall) is small. From Figure 3, it can be seen that the temperature gradient near the left wall and the top wall is higher than that in other regions, which causes a high temperature difference driving force near there. Similarly, the temperature gradient in other regions is small, causing only a low temperature difference driving force in that vicinity. In addition, it can be seen that the same driving force line at a high Rayleigh number becomes more crooked than that at a low Rayleigh number. This is because the driving force is caused by the temperature difference (temperature gradient); a bigger temperature gradient causes the same driving force line to become more crooked. It can be seen from Figure 3 that isotherms are more crooked at a higher Rayleigh number, and the isotherm changes correspond to the changes of temperature gradient. Thus, the conclusion that the same driving force line at a high Rayleigh number becomes more crooked than that that at a low Rayleigh number is obtained. Figures 6 and 7 give the density distribution of the water phase at Ra = 1 × 10 3 and Ra = 1 × 10 5 . For a low Rayleigh number (Ra = 1 × 10 3 ), when the water near the left wall is heated, its density decreases and flows upward, so the density of water near the top right corner also becomes smaller. Then when the water is cooled by the top wall, the density of the water becomes larger. Then the denser water flows downward to the lower right corner, and so, the density of water in the lower right corner is larger than that in other regions. Because the temperature gradient (corresponding to the temperature difference driving force) is small and the temperature is high in the lower left corner, the density of water in the lower left corner is thus low. For a high Rayleigh number (Ra = 1 × 10 5 ), the temperature gradient and the corresponding driving force become larger, then the lower-density water, including that in the lower left corner, rises to the top right corner. The denser water is cooled by the top wall and flows downward to the lower right corner, and the area where the denser water in the lower right corner becomes larger. Figures 8 and 9 respectively present the nanoparticle distribution of nanofluid with volume fractions at Ra = 1 × 10 3 and Ra = 1 × 10 5 . For a low Rayleigh number (Ra = 1 × 10 3 ), the driving force along the left wall is upward, and many nanoparticles are driven to the top right corner, which contributes to the high nanoparticle volume fraction in the top right corner. However, the temperature gradient in the lower left corner is small and causes a correspondingly small temperature difference driving force. Thus, many nanoparticles are left in the lower left corner, which contributes to the high nanoparticle volume fraction in the lower left corner. There is a large temperature gradient in the lower right corner, and the large driving force displaces the nanoparticles off the lower right corner, which contributes to the low nanoparticle volume fraction in the lower right corner. For a high Rayleigh number (Ra = 1 × 10 5 ), the convection heat transfer is enhanced and the velocity of the nanofluid becomes larger, and the temperature gradient and the corresponding driving force become bigger. Thus, many nanoparticles from the bottom are driven to the top by the driving force, which contributes to the low nanoparticle volume fraction at the bottom and a high nanoparticle volume fraction at the top. In addition, we can see that the nanoparticle volume fraction distribution is opposite to that of the water-phase density distribution. From Table 4, we can see that the temperature difference driving force is the biggest one, and the changes of the water-phase density and the inhomogeneous nanoparticle distribution are mainly due to the driving force. Through the above analysis, it is found that the nanoparticles migrate to locations where the water density is small, and thus, the conclusion that the nanoparticle volume fraction distribution is opposite to that of the waterphase density distribution is obtained.
It is also found that almost all the isolines behave with oscillations in Figures 6, 7, 8, 9, but smooth isolines are given in Figures 3 and 5. Due to the ruleless Brownian movement of nanoparticles, it is difficult for nanofluid to achieve a complete equilibrium state, which is the difference compared with other common two-phase fluids. In order to expediently judge the equilibrium state and save time, we choose the temperature equilibrium states of water phase and nanoparticle phase as the whole nanofluid equilibrium state in the computation. When the water-phase and nanoparticle-phase temperatures all achieve equilibrium state, the whole nanofluid (temperature distribution, velocity vectors, density distribution, and nanoparticle volume fraction distribution) is considered as being in an equilibrium state. Hence, the temperature isolines in Figures 3 and 5 look smooth due to a complete equilibrium state, and the density distribution in Figures 6  and 7 and nanoparticle volume fraction distribution in Figures 8 and 9 behave with oscillations due to an approximate equilibrium state. Although the interparticle interaction forces have little effect on heat transfer, they play an important role on the nanoparticle distribution. Figure 10 shows the Nusselt number distribution along the heated surface using Al 2 O 3 -water nanofluid at Ra = 10 3 . It can be seen that the Nusselt number along the heated surface increases with nanoparticle volume fraction at low Y (0 < Y < 0.58) and decreases with nanoparticle volume fraction at high Y (0.58 < Y < 1). Because  the heat transfer is more sensitive to thermal conductivity than viscosity at low Y, while it is more sensitive to viscosity than thermal conductivity at high Y. Figure 11 shows Nusselt number distribution along the heated surface using Al 2 O 3 -water nanofluid at Ra = 10 5 . It can be seen that the Nusselt number along the heated surface increases with nanoparticle volume fraction at low Y (0 < Y < 0.875) and decreases with nanoparticle volume fraction at high Y (0.875 < Y < 1). Compared with Figure 7, the Nusselt number becomes larger, and the enhanced heat transfer section also gets longer. The high Rayleigh number increases the velocity and then enhances the heat transfer. Figure 12 presents the average Nusselt numbers at different Rayleigh numbers. Although the Nusselt number distribution along the heated surface increases with nanoparticle volume fraction in one section and decreases in the other section, the average Nusselt numbers at Ra = 10 3 and Ra = 10 5 both increase with nanoparticle volume fraction. For this square enclosure (left wall is kept at a high constant temperature (T H ), and top cold wall is kept at a low constant temperature (T C )), adding nanoparticles can enhance the average heat transfer at both a low and a high Rayleigh number. In addition, the enhancement of the average Nusselt numbers is much more pronounced at a high Rayleigh number than at a low Rayleigh number.

Conclusion
A 2D two phase Lattice Boltzmann model has been developed for nanofluids and the simulation results of this two-phase Lattice Boltzmann model are in good agreement with published experimental results. This model is applied to investigate the natural convection of a square enclosure filled with Al 2 O 3 nanofluid. The effects of different nanoparticle fractions and Rayleigh numbers on natural convection heat transfer of nanofluid are investigated. In addition the effects of forces on the nanoparticles volume fraction distribution and the heat transfer are also investigated.
It is found that the Nusselt number distribution along the heated surface firstly increases, and then decreases with Y at both low and high Rayleigh numbers. Average Nusselt numbers of the whole square enclosure both increase with nanoparticles volume fraction at a low and a high Rayleigh number. In addition, the enhancement of the average Nusselt numbers is much more pronounced at a high Rayleigh number than at a low Rayleigh number.
It is found that the temperature difference driving force is the biggest force and has the greatest effect on nanoparticle volume fraction distribution. For a low Rayleigh number, the nanoparticle volume fraction is low in the lower right corner and high in the top right corner and lower left corner. For a high Rayleigh number, the nanoparticle volume fraction is low at the bottom and high at the top.
Apart from the temperature difference driving force, Brownian force, interaction potential force, and gravitybuoyancy force contribute to the enhanced natural convective heat transfer, while the drag force contributes to the attenuation of heat transfer.
Nomenclature a radius of nanoparticle (m) A Hamaker constant B a weight coefficient c reference lattice velocity c s lattice sound velocity c p specific heat capacity (J/kg K) e α lattice velocity vector f σ α density distribution function f σeq α local equilibrium density distribution function F σ 0 α dimensionless external force in direction of lattice velocity F σ α dimensionless total interparticle interaction forces F S dimensionless temperature difference driving forces F B dimensionless Brownian force F H dimensionless gravity and buoyancy force F D dimensionless drag force F A dimensionless interaction potential force g dimensionless gravitational acceleration G dimensionless effective external force G i Gaussian random number h aβ convective heat transfer coefficient (W/(m 2 K)) H dimensionless characteristic length of the square cavity k thermal conductivity coefficient (W/m/K) k B Boltzmann constant L cc center-to-center distance between particles (m) Ma Mach number m σ mass of a single nanoparticle (kg) n i number of the particles within the adjacent lattice i Nu Nusselt number Pr Prandtl number r position vector Ra Rayleigh number t time (s) T σ α temperature distribution function T σeq α local equilibrium temperature distribution function T dimensionless temperature T 0 dimensionless average temperature (T 0 = (T H + T C )/2) T H dimensionless hot temperature T C dimensionless cold temperature u σ dimensionless macro-velocity u c dimensionless characteristic velocity of natural convection