# Anisotropic thermal lattice Boltzmann simulation of 2D natural convection in a square cavity

François Dubois, Chao-An Lin, Mohamed Mahdi Tekitek

► To cite this version:

François Dubois, Chao-An Lin, Mohamed Mahdi Tekitek. Anisotropic thermal lattice Boltzmann simulation of 2D natural convection in a square cavity. Computers and Fluids, 2016, 124, pp.278-287. 10.1016/j.compfluid.2015.10.015. hal-03177125

HAL Id: hal-03177125

https://hal.science/hal-03177125v1

Submitted on 15 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Anisotropic thermal lattice Boltzmann simulation of 2D natural convection in a square cavity

François Dubois<sup>a,b</sup>, Chao-An Lin<sup>c</sup>, Mohamed Mahdi Tekitek<sup>c,d,*</sup>

Natural convection in a square cavity is simulated by multiple relaxation time (MRT) lattice Boltzmann method (LBM) with a separate distribution function to solve for the temperature distribution. The Rayleigh numbers examined range from $Ra = 10^3$ to $Ra = 10^6$. The simulations are performed for anisotropic thermal case and compared to isotropic thermal case.

## 1. Introduction

Lattice Boltzmann method (LBM) [2,26,37,39] has been successfully applied to various hydrodynamic problems and the major advantage of the LBM is explicit formulation. However, its application to non-isothermal problem is limited because of the numerical instability for thermal models [32]. In general, there are three thermal lattice Boltzmann methods (TLBM) named the multi speed approach [1], the passive scalar approach and the double population approach.

The multi-speed approach adopts a single distribution function in order to obtain the macroscopic dynamic and thermal equations [1]. However, the Prandtl number is fixed and this approach suffers from lack of numerical stabilities.

The passive scalar approach also called hybrid method, consists of solving velocity field using LBM and the macroscopic temperature is solved by different numerical methods (e.g. finite difference or finite volume) [24,30]. This approach is more stable than the multi speed approach. But it has two disadvantages: first the viscous heat dissipation and compression work done by pressure cannot be incorporated, and second the simplicity of LBM is lost.

The double population method, is first used by He et al. [19]. This approach can be regarded as another version of the passive scalar method. In fact to solve macroscopic temperature another LBM distribution is used. This model has a better numerical stability than the multi speed approach, and the viscous heat dissipation and compression work done by the pressure can be solved implicitly. Peng et al. [35] proposed a simplified thermal energy distribution model where the compression work done by the pressure and the viscous heat dissipation are neglected. By introducing a forcing function, Guo et al. [17] proposed a thermal lattice BGK equation with viscous heat dissipation in the incompressible limit.

The thermally driven cavity with adiabatic top and bottom walls (also called natural convection in a square cavity) is a classical benchmark to examine the accuracy of the scheme. The solution is given for 4 values of the Rayleigh number ($Ra$), ($Ra = 10^3, 10^4, 10^5$ and $10^6$). The value of the Prandtl number ($Pr$) is equal to 0.71, which corresponds to a cavity filled by air. The reference solution of this problem is given by De Vahl Davis [4].

To validate double population LBM method few researchers [10, 16, 18, 21, 29] have carried out the above problem. We note here that most of this works are using simple relaxation time (SRT), also called Lattice Boltzmann Bhatngar-Gross-Krook (LBGK). This is due to extreme simplicity of this method. Even that LBGK suffers from lack of numerical stability and inaccuracy in implementing boundary conditions.

<footnote><sup>a</sup> Department of Mathematics, Paris-Sud University, Orsay, France</footnote>

<footnote><sup>b</sup> Conservatoire National des Arts et Métiers, LMSSC, Paris, France</footnote>

<footnote><sup>c</sup> Department of Power Mechanical Engineering, National Tsing Hua University, Hsinchu 30013, Taiwan</footnote>

<footnote><sup>d</sup> Department of Mathematics, Faculty of Sciences of Tunis, University of Tunis El Manar, 2092 El Manar, Tunisia</footnote>

<footnote>* Corresponding author at: Department of Mathematics, Faculty of Sciences of Tunis, University of Tunis El Manar, 2092 El Manar, Tunisia. Tel: +216 98901903.

E-mail address: mahdi.tekitek@math.u-psud.fr, memotek@gmail.com (M.M. Tekitek).</footnote>
In this paper we present a double population approach using multiple relaxation time lattice Boltzmann method (MRT-LBM) [38] with D2Q9 lattice model for solving velocity field and another D2Q9 for solving macroscopic temperature. The choice of D2Q9 model for thermal is to be able to model anisotropic thermal diffusion. We note that anisotropic thermal diffusion have many applications in diffusivity study of gas diffusion layer (see [34] for measurements and [11] for computations), advection and anisotropic-dispersion equation, porous media (see e. g. [7,13,15]).

First we consider natural convection in a square cavity when the flow is laminar (i.e. Rayleigh number is less than $10^6$.). To validate our model, we choose isotropic thermal diffusivity (i.e. diffusivity in x direction $\kappa_x$ is equal to diffusivity in y direction $\kappa_y$). Then we consider anisotropic thermal diffusion. In fact we consider two cases $\kappa_x = \kappa_y/2$ and $\kappa_x = 2\kappa_y$, and we compare the solution to the isotropic one.

This paper is organized as follows. In Section 2, a brief overview of the MRT D2Q9 for advection-diffusion and the MRT D2Q9 for fluid problem. After in Section 3 we introduce the thermal LBM for the simulation of a Boussinesq fluid in a square cavity. In Section 4, results are presented and discussed. Finally, in Section 5 we conclude.

## 2. Multi relaxation time lattice Boltzmann method

### 2.1. Dynamic field

The multi relaxation time (MRT) lattice Boltzmann method [23] can be expressed as:

$$m_l^*(\vec{x}, t) = m_l(\vec{x}, t) - S_{lj}[m_j(\vec{x}, t) - m_j^{eq}(\vec{x}, t)] \quad (1)$$

$$f_i(\vec{x} + \vec{e}_i \Delta t, t + \Delta t) = M_{il}^{-1} m_i^*(\vec{x}, t) \quad (2)$$

$m_i^*$, Eq. (1), is the collision at the moment space and Eq. (2) represents the streaming operation. Here, **M** is a matrix that transforms the distribution function **f** to the velocity moment, **m**=**M****f**, and **S** is the relaxation matrix. These will be defined later.

Based on the particle distribution functions, the macroscopic density and velocity are defined as:

$$\sum_i f_i = \rho, \quad \sum_i f_i \vec{e}_i = \rho \vec{u}. \quad (3)$$

For the present 2D applications, D2Q9 model are adopted to model fluid problems and the particle speed $\vec{e}_i$ are defined as,

$$\left\{ \begin{array}{l} \vec{e}_0 = 0, \\ \vec{e}_i = (\cos[\pi(i-1)/2], \sin[\pi(i-1)/2])c, \\ \quad \text{for } i = 1, 2, 3, 4, \\ \vec{e}_i = (\cos[\pi(i-4-1/2)/2], \sin[\pi(i-4-1/2)/2])\sqrt{2}c, \\ \quad \text{for } i = 5, 6, 7, 8, \end{array} \right.$$

where $c = dx/dt$ is the lattice speed, and $dx$ and $dt$ are the lattice width and time step, respectively. Here, $dt$ is chosen to be equal to $dx$, thus $c = 1$. Moreover, the speed of sound is $C_s = c/\sqrt{3}$.

The transformation matrix **M** and the velocity moment vector **m** are defined as,

$$\underbrace{\begin{bmatrix} m_0 (\rho) \\ m_1 (\epsilon) \\ m_2 (\epsilon) \\ m_3 (j_x) \\ m_4 (q_x) \\ m_5 (j_y) \\ m_6 (q_y) \\ m_7 (p_{xx}) \\ m_8 (p_{xy}) \end{bmatrix}}_{\mathbf{m}} = \underbrace{\begin{bmatrix} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ -4 & -1 & -1 & -1 & -1 & 2 & 2 & 2 \\ 4 & -2 & -2 & -2 & -2 & 1 & 1 & 1 \\ 0 & 1 & 0 & -1 & 0 & 1 & -1 & 1 \\ 0 & -2 & 0 & 2 & 0 & 1 & -1 & 1 \\ 0 & 0 & 1 & 0 & -1 & 1 & 1 & -1 \\ 0 & 0 & -2 & 0 & 2 & 1 & 1 & -1 \\ 0 & 1 & -1 & 1 & -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & -1 & 1 \end{bmatrix}}_{\mathbf{M}} \underbrace{\begin{bmatrix} f_0 \\ f_1 \\ f_2 \\ f_3 \\ f_4 \\ f_5 \\ f_6 \\ f_7 \\ f_8 \end{bmatrix}}_{\mathbf{f}}$$

and the equilibria of the velocity moments $\mathbf{m}^{eq}$ are,

$$\left\{ \begin{aligned} \rho^{eq} &= \rho, & e^{eq} &= -2\rho + \frac{3}{\rho}(j_x^2 + j_y^2), & \epsilon^{eq} &= \rho - \frac{3}{\rho}(j_x^2 + j_y^2), \\ j_x^{eq} &= j_x, & q_x^{eq} &= -j_x, & j_y^{eq} &= j_y, & q_y^{eq} &= -j_y, \\ p_{xx}^{eq} &= \frac{1}{\rho}(j_x^2 - j_y^2), & p_{xy}^{eq} &= \frac{j_x j_y}{\rho}. \end{aligned} \right. \quad (5)$$

The relaxation matrix **S** is a diagonal matrix, i.e.,

$$\mathbf{S} = \mathrm{diag}[s_0, s_1, s_2, s_3, s_4, s_5, s_6, s_7, s_8] \quad (6)$$

where $s_0 = s_3 = s_5 = 0$ enforces mass and momentum conservation before and after collision [23]. Here, $s_4 = s_6 = s_q$, $s_7 = s_8 = s_v$, thus the viscosity formulation is the same as that by the SRT model as shown in [23], i.e. the corresponding kinematic viscosity is $\nu = (\frac{1}{s_v} - \frac{1}{2})C_s^2 dt$ in the simulation. More specifically, we choose to use the following relationship between two relaxation rates: $s_q = 8 \frac{(2-s_v)}{(8-s_v)}$ where $s_q = s_4 = s_6$. See [8,14] for more details.

We note that, the MRT model can recover to SRT model if $s_1 = s_2 = s_4 = s_6 = s_7 = s_8 = s_v$.

### 2.2. Thermal field

The thermal field is modeled using the passive scalar approach to enhance the numerical stability, where a separate distribution function is used to solve for the temperature distribution [27,35,36]. The D2Q9 model introduced in the above section is adopted. The evolution of the scalar MRT LB scheme is given as:

$$\tilde{m}_l^*(\vec{x}, t) = \tilde{m}_l(\vec{x}, t) - \sigma_{lj} [\tilde{m}_j(\vec{x}, t) - \tilde{m}_j^{eq}(\vec{x}, t)] \quad (7)$$

$$g_i(\vec{x} + \vec{e}_i \Delta t, t + \Delta t) = M_{il}^{-1} \tilde{m}_i^*(\vec{x}, t) \quad (8)$$

Here, $f_i$ is replaced by $g_i$ in Eq. (2), because $g_i$ is now the energy distribution function. The transformation matrix is the same as in equation (4), thus $\tilde{m} = \mathbf{M}g$. Again, $\tilde{m}_i^*$ is the scalar collision at the moment space. $\sigma$ is the diagonal relaxation matrix, i.e.

$$\sigma = \mathrm{diag}[\sigma_0, \sigma_1, \sigma_2, \sigma_3, \sigma_4, \sigma_5, \sigma_6, \sigma_7, \sigma_8] \quad (9)$$

where $\sigma_0 = 0$ enforces energy conservation, $\tilde{m}_0 = \sum_i g_i = T$, before and after collision.

The equilibrium values $\tilde{m}_i^{eq}$ of the nonconserved moments are given by (see [22] for more details):

$$\tilde{m}_1^{eq} = \tilde{\alpha}T + 3T(u^2 + v^2),$$

$$\tilde{m}_2^{eq} = \tilde{\beta}T,$$

$$\tilde{m}_3^{eq} = uT,$$

$$\tilde{m}_4^{eq} = uT(-1 + 3(u^2 + v^2)),$$

$$\tilde{m}_5^{eq} = vT,$$

$$\tilde{m}_6^{eq} = vT(-1 + 3(u^2 + v^2))$$

$$\tilde{m}_7^{eq} = a_x T + T(u^2 - v^2),$$

$$\tilde{m}_8^{eq} = a_y T + T(uv),$$

where $\mathbf{V} = (u, v)$ is the dynamic field.

Using Taylor expansion [5,6] or Chapman–Enskog procedure [12], the advection diffusion [15] equation with an-isotropic coefficient up to order two in $\Delta t$ can be expressed as:

$$\frac{\partial T}{\partial t} + U \frac{\partial T}{\partial t} + V \frac{\partial T}{\partial y} = \frac{C_s^2 dt}{2} \left( \frac{1}{2} - \frac{1}{\sigma_3} \right) (\tilde{\alpha} + 3a_x + 4) \frac{\partial^2 T}{\partial x^2} \\ + \frac{C_s^2 dt}{2} \left( \frac{1}{2} - \frac{1}{\sigma_5} \right) (\tilde{\alpha} + 3a_x + 4) \frac{\partial^2 T}{\partial y^2} \\ + \frac{3C_s^2 dt}{2} a_y \left( \frac{1}{\sigma_3} + \frac{1}{\sigma_5} - 1 \right) \frac{\partial^2 T}{\partial xy},$$
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Configuration of natural convection in a square cavity.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Left: Boundary node $x_b$ in the bottom of the domain $\Omega$. Right: Boundary node $x_b$ in the right of the domain $\Omega$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Isobars (P) of flow fields for Pr=0.71. From top to bottom $Ra = 10^3, 10^4, 10^5$, and $10^6$ respectively for mesh size $105^2, 155^2, 205^2$ and $255^2$. From left to right : (left) $\kappa_x = \kappa_y/2$, (center) isotropic case $\kappa_x = \kappa_y$, (right) $\kappa_x = 2\kappa_y$.</figcaption></figure>
<table><caption>**Table 1**

Grid dependence and order of accuracy for numerical results of simulating 2-D natural convection of air in a square cavity for the case Ra = 10⁶.</caption><thead><tr><th>Mesh size</th><th>N&#773;u</th><th>u<sub>max</sub></th><th>y</th><th>v<sub>max</sub></th><th>x</th></tr></thead><tbody><tr><td>47<sup>2</sup></td><td>9.5729</td><td>64.6766</td><td>0.8414</td><td>222.8704</td><td>3.6585.10<sup>-2</sup></td></tr><tr><td>57<sup>2</sup></td><td>9.2079</td><td>64.8453</td><td>0.8508</td><td>218.5256</td><td>4.3859.10<sup>-2</sup></td></tr><tr><td>81<sup>2</sup></td><td>8.8913</td><td>64.8181</td><td>0.8456</td><td>218.6303</td><td>4.3209.10<sup>-2</sup></td></tr><tr><td>161<sup>2</sup></td><td>8.7828</td><td>64.8375</td><td>0.8478</td><td>220.1784</td><td>4.0372.10<sup>-2</sup></td></tr><tr><td>225<sup>2</sup></td><td>8.8006</td><td>64.8393</td><td>0.8511</td><td>220.7288</td><td>3.7777.10<sup>-2</sup></td></tr><tr><td>321<sup>2</sup></td><td>8.8139</td><td>64.8403</td><td>0.8489</td><td>220.5115</td><td>3.8940.10<sup>-2</sup></td></tr><tr><td>615<sup>2</sup></td><td>8.8226</td><td>64.8391</td><td>0.8495</td><td>220.5741</td><td>3.8211.10<sup>-2</sup></td></tr><tr><td>Reference solution</td><td>8.8241</td><td>64.8372</td><td>0.8495</td><td>220.5739</td><td>3.7993.10<sup>-2</sup></td></tr><tr><td><strong>Accuracy</strong></td><td>2.109</td><td>2.0154</td><td>1.9261</td><td>2.0592</td><td>1.7571</td></tr><tr><td>Davis [4]</td><td>8.817</td><td>64.630</td><td>0.850</td><td>219.360</td><td>3.8.10<sup>-2</sup></td></tr><tr><td>Le Quéré [25]</td><td>8.825</td><td>64.483</td><td>0.850</td><td>220.559</td><td>3.8.10<sup>-2</sup></td></tr></tbody></table>

<table><caption>**Table 2**

Grid dependence and order of accuracy for numerical results of simulating 2-D natural convection of air in a square cavity for the case Ra = 10⁸.</caption><thead><tr><th>Mesh size</th><th>Nu</th><th>U<sub>max</sub></th><th>Y</th><th>V<sub>max</sub></th><th>X</th></tr></thead><tbody><tr><td>495<sup>2</sup></td><td>29.97</td><td>313.97</td><td>0.926</td><td>2219.3</td><td>0.011</td></tr><tr><td>987<sup>2</sup></td><td>30.20</td><td>319.81</td><td>0.927</td><td>2222.8</td><td>0.011</td></tr><tr><td>2016<sup>2</sup></td><td>30.22</td><td>321.59</td><td>0.928</td><td>2222.7</td><td>0.012</td></tr><tr><td>Le Quéré [25]</td><td>30.22</td><td>321.88</td><td>0.928</td><td>2222.2</td><td>0.012</td></tr></tbody></table>

<table><caption>**Table 3**

Comparison of predicted numerical results. Davis [4], Mayne et al. [31], Liu et al. [28], Dixit et al. [10], Kuznik et al. [21], Mezrhab et al. [33] and Wang et al. [38].</caption><thead><tr><th>Ra</th><th>[4]</th><th>[31]</th><th>[28]</th><th>[21]</th><th>[33]</th><th>[38]</th><th>Present</th></tr></thead><tbody><tr><td rowspan="5">10<sup>3</sup></td><td>u<sub>max</sub></td><td>3.649</td><td>3.649</td><td>3.649</td><td>3.636</td><td>3.667</td><td>3.649</td><td>3.649</td></tr><tr><td>y</td><td>0.813</td><td>0.812</td><td>0.810</td><td>0.809</td><td>-</td><td>0.813</td><td>0.814</td></tr><tr><td>v<sub>max</sub></td><td>3.697</td><td>3.696</td><td>3.698</td><td>3.686</td><td>3.714</td><td>3.697</td><td>3.697</td></tr><tr><td>x</td><td>0.178</td><td>0.179</td><td>0.180</td><td>0.174</td><td>-</td><td>0.178</td><td>0.176</td></tr><tr><td>N&#773;u</td><td>1.117</td><td>1.114</td><td>1.115</td><td>1.117</td><td>1.112</td><td>1.117</td><td>1.117</td></tr><tr><td rowspan="5">10<sup>4</sup></td><td>u<sub>max</sub></td><td>16.178</td><td>16.179</td><td>16.154</td><td>16.167</td><td>16.202</td><td>16.183</td><td>16.188</td></tr><tr><td>y</td><td>0.823</td><td>0.823</td><td>0.820</td><td>0.821</td><td>-</td><td>0.823</td><td>0.822</td></tr><tr><td>v<sub>max</sub></td><td>19.617</td><td>19.617</td><td>19.614</td><td>19.597</td><td>19.644</td><td>19.627</td><td>19.632</td></tr><tr><td>x</td><td>0.119</td><td>0.119</td><td>0.120</td><td>0.120</td><td>-</td><td>0.118</td><td>0.119</td></tr><tr><td>N&#773;u</td><td>2.243</td><td>2.259</td><td>2.229</td><td>2.246</td><td>2.241</td><td>2.244</td><td>2.243</td></tr><tr><td rowspan="5">10<sup>5</sup></td><td>u<sub>max</sub></td><td>34.730</td><td>34.774</td><td>34.508</td><td>34.962</td><td>34.805</td><td>34.743</td><td>34.748</td></tr><tr><td>y</td><td>0.855</td><td>0.853</td><td>0.855</td><td>0.854</td><td>-</td><td>0.854</td><td>0.856</td></tr><tr><td>v<sub>max</sub></td><td>68.590</td><td>68.692</td><td>68.595</td><td>68.578</td><td>68.630</td><td>68.631</td><td>68.652</td></tr><tr><td>x</td><td>0.066</td><td>0.066</td><td>0.065</td><td>0.067</td><td>-</td><td>0.065</td><td>0.065</td></tr><tr><td>N&#773;u</td><td>4.519</td><td>4.483</td><td>4.489</td><td>4.518</td><td>4.519</td><td>4.521</td><td>4.517</td></tr><tr><td rowspan="5">10<sup>6</sup></td><td>u<sub>max</sub></td><td>64.630</td><td>64.691</td><td>63.456</td><td>64.133</td><td>64.793</td><td>64.827</td><td>64.842</td></tr><tr><td>y</td><td>0.850</td><td>0.846</td><td>0.848</td><td>0.860</td><td>-</td><td>0.849</td><td>0.849</td></tr><tr><td>v<sub>max</sub></td><td>219.360</td><td>220.833</td><td>219.788</td><td>220.537</td><td>219.663</td><td>220.550</td><td>220.669</td></tr><tr><td>x</td><td>0.037</td><td>0.038</td><td>0.036</td><td>0.038</td><td>-</td><td>0.037</td><td>0.037</td></tr><tr><td>N&#773;u</td><td>8.799</td><td>8.881</td><td>8.750</td><td>8.792</td><td>8.817</td><td>8.819</td><td>8.806</td></tr></tbody></table>

Let $a_x = a_y = 0$, the above advection diffusion equation is reduced to equation anisotropic diffusion coefficient and is expressed as:

$$ \frac{\partial T}{\partial t} + \mathbf{V} \cdot \nabla T = \kappa_x \frac{\partial^2 T}{\partial x^2} + \kappa_y \frac{\partial^2 T}{\partial y^2}, \qquad (10) $$

where the values of x-diffusivity $κ_x$ and y-diffusivity $κ_y$ are:

$$ \kappa_x = C_s^2 dt \frac{\tilde{\alpha} + 4}{2} \left( \frac{1}{\sigma_3} - \frac{1}{2} \right), \quad \kappa_y = C_s^2 dt \frac{\tilde{\alpha} + 4}{2} \left( \frac{1}{\sigma_5} - \frac{1}{2} \right). \qquad (11) $$

Here, $\tilde{\alpha}$ and $\tilde{\beta}$ are -2 and 1, respectively. The present D2Q9 advection diffusion equation can accommodate thermal problem with isotropic (i. e. $κ_x = κ_y = κ$) diffusivity [7]. When $U = V = 0$, Eq. (10) reduces to diffusion equation.

**Remark** Note here that using D2Q5 for thermal problem (see [3]) and isotropic diffusivity is sufficient, faster and requires less memory than using D2Q9. The advantages of using Thermal D2Q9 is the ability

<table><caption>**Table 4**

Convergence for the case Ra$_y$ = 10⁶ and $κ_x = \frac{κ_y}{2}$.</caption><thead><tr><th>Mesh size</th><th>N&#773;u</th><th>u<sub>max</sub></th><th>y</th><th>v<sub>max</sub></th><th>x</th></tr></thead><tbody><tr><td>47<sup>2</sup></td><td>23.1332</td><td>50.5699</td><td>0.8919</td><td>168.2357</td><td>3.2749.10<sup>-2</sup></td></tr><tr><td>57<sup>2</sup></td><td>22.2512</td><td>50.7018</td><td>0.9019</td><td>164.9560</td><td>3.9260.10<sup>-2</sup></td></tr><tr><td>81<sup>2</sup></td><td>21.4861</td><td>50.6805</td><td>0.8963</td><td>165.0350</td><td>4.3209.10<sup>-2</sup></td></tr><tr><td>161<sup>2</sup></td><td>21.2239</td><td>50.6957</td><td>0.8987</td><td>166.2036</td><td>3.6139.10<sup>-2</sup></td></tr><tr><td>225<sup>2</sup></td><td>21.2805</td><td>50.6999</td><td>0.9000</td><td>166.5744</td><td>3.3333.10<sup>-2</sup></td></tr><tr><td>321<sup>2</sup></td><td>21.2991</td><td>50.6979</td><td>0.8998</td><td>166.4551</td><td>3.4857.10<sup>-2</sup></td></tr><tr><td>615<sup>2</sup></td><td>21.3201</td><td>50.6970</td><td>0.9005</td><td>166.5023</td><td>3.4205.10<sup>-2</sup></td></tr><tr><td>Reference solution</td><td>21.3237</td><td>50.6955</td><td>0.9005</td><td>166.5022</td><td>3.4009.10<sup>-2</sup></td></tr><tr><td><strong>Accuracy</strong></td><td>2.1009</td><td>2.0582</td><td>1.9368</td><td>2.0855</td><td>1.7009</td></tr></tbody></table>

<table><caption>**Table 5**

Convergence for the case Ra$_y$ = 10⁶ and $κ_x = 2κ_y$.</caption><thead><tr><th>Mesh size</th><th>N&#773;u</th><th>u<sub>max</sub></th><th>y</th><th>v<sub>max</sub></th><th>x</th></tr></thead><tbody><tr><td>47<sup>2</sup></td><td>3.8261</td><td>119.7447</td><td>0.8569</td><td>285.6379</td><td>4.0420.10<sup>-2</sup></td></tr><tr><td>57<sup>2</sup></td><td>3.6802</td><td>120.0571</td><td>0.8665</td><td>280.0695</td><td>4.8457.10<sup>-2</sup></td></tr><tr><td>81<sup>2</sup></td><td>3.5537</td><td>120.0067</td><td>0.8612</td><td>280.2036</td><td>4.7738.10<sup>-2</sup></td></tr><tr><td>161<sup>2</sup></td><td>3.5103</td><td>120.0426</td><td>0.8634</td><td>282.1877</td><td>4.4604.10<sup>-2</sup></td></tr><tr><td>225<sup>2</sup></td><td>3.5174</td><td>120.0460</td><td>0.8668</td><td>282.8932</td><td>4.1737.10<sup>-2</sup></td></tr><tr><td>321<sup>2</sup></td><td>3.5227</td><td>120.0478</td><td>0.8645</td><td>282.6147</td><td>4.3022.10<sup>-2</sup></td></tr><tr><td>615<sup>2</sup></td><td>3.5262</td><td>120.0456</td><td>0.8652</td><td>282.6949</td><td>4.2216.10<sup>-2</sup></td></tr><tr><td>Reference solution</td><td>3.5268</td><td>120.0421</td><td>0.8652</td><td>282.6946</td><td>4.1976.10<sup>-2</sup></td></tr><tr><td><strong>Accuracy</strong></td><td>2.1747</td><td>2.0581</td><td>1.8891</td><td>2.3270</td><td>1.7507</td></tr></tbody></table>

<table><caption>**Table 6**

Predicted numerical results for different cases: $κ_x = \frac{κ_y}{2}$, $κ_x = κ_y$ and for $κ_x = 2κ_y$.</caption><thead><tr><th>Ra</th><th>κ<sub>x</sub> = <sup>κ<sub>y</sub></sup>&frasl;<sub>2</sub></th><th>κ<sub>x</sub> = κ<sub>y</sub></th><th>κ<sub>x</sub> = 2κ<sub>y</sub></th></tr></thead><tbody><tr><td rowspan="3">10<sup>3</sup></td><td>u<sub>max</sub></td><td>3.3705</td><td>3.6496</td><td>3.8185</td></tr><tr><td>y</td><td>0.8142</td><td>0.8142</td><td>0.8142</td></tr><tr><td>v<sub>max</sub></td><td>3.4515</td><td>3.6973</td><td>3.8428</td></tr><tr><td rowspan="5">10<sup>4</sup></td><td>x</td><td>0.1761</td><td>0.1761</td><td>0.1857</td></tr><tr><td>N&#773;u</td><td>2.4957</td><td>1.1179</td><td>0.5226</td></tr><tr><td>u<sub>max</sub></td><td>12.3628</td><td>16.1881</td><td>21.1512</td></tr><tr><td>y</td><td>0.8290</td><td>0.8225</td><td>0.8225</td></tr><tr><td>v<sub>max</sub></td><td>16.0159</td><td>19.6323</td><td>24.0035</td></tr><tr><td rowspan="5">10<sup>5</sup></td><td>x</td><td>0.1064</td><td>0.1193</td><td>0.1322</td></tr><tr><td>N&#773;u</td><td>5.3711</td><td>2.2438</td><td>0.9261</td></tr><tr><td>u<sub>max</sub></td><td>23.5783</td><td>34.7486</td><td>56.0032</td></tr><tr><td>y</td><td>0.8560</td><td>0.8560</td><td>0.8609</td></tr><tr><td>v<sub>max</sub></td><td>53.5863</td><td>68.6527</td><td>86.0559</td></tr><tr><td rowspan="5">10<sup>6</sup></td><td>x</td><td>0.0609</td><td>0.0658</td><td>0.0707</td></tr><tr><td>N&#773;u</td><td>10.9325</td><td>4.5177</td><td>1.8262</td></tr><tr><td>u<sub>max</sub></td><td>50.6999</td><td>64.8428</td><td>120.0525</td></tr><tr><td>y</td><td>0.9000</td><td>0.8490</td><td>0.8647</td></tr><tr><td>v<sub>max</sub></td><td>166.5744</td><td>220.6695</td><td>282.8172</td></tr><tr><td rowspan="2">10<sup>7</sup></td><td>x</td><td>0.0333</td><td>0.0372</td><td>0.0411</td></tr><tr><td>N&#773;u</td><td>21.2805</td><td>8.8062</td><td>3.5197</td></tr></tbody></table>

to model non isotropic thermal problem and the possibility to cancel the dependence of thermal diffusivity on the advection velocity (for more details see [22]).

## 2.3. Coupling of dynamic and thermal fields

With the Boussinesq approximation, the buoyancy term is assumed to depend linearly on the temperature as,

$$ F_y = \beta g_y (T - T_{ref}) \qquad (12) $$

where $β$ is the thermal expansion coefficient, $g_y$ is the acceleration due to gravity, and $T_{ref}$ is the reference temperature.

To perform the coupling, the buoyancy force $F_y$ is added in moment space before and after the collision process of the LB scheme as described by Eq. (1). The procedure goes like this [9]:
<figure><img src="image_5.png" /><figcaption>Fig. 4. First component of velocity (U) of flow for Pr=0.71. From top to bottom Ra = 10³, 10⁴, 10⁵, and 10⁶ respectively for mesh size 105², 155², 205² and 255². From left to right : (left) κx = κy/2, (center) isotropic case κx = κy, (right) κx = 2κy.</figcaption></figure>

* The y direction momentum ($\vec{j}_y$) and energy flux ($q_y$) are modified by adding half of the external force $F_y$, i.e.

$$ \vec{j}_y = j_y + \frac{\Delta t}{2} F_y, \quad \bar{q}_y = q_y - \frac{\Delta t}{2} F_y, $$

* Compute the equilibrium moments in Eq. (5) using $\vec{j}_y$ and $\bar{q}_y$ to replace $j_y$ and $q_y$.

* Perform collision in Eq. (1).

* Post collision y direction momentum and energy flux are modified by adding another half of the external force, i.e.,

$$ \vec{j}_y^* = j_y^* + \frac{\Delta t}{2} F_y, \quad \bar{q}_y^* = q_y^* - \frac{\Delta t}{2} F_y. $$

* Perform streaming in Eq. (2) using $\vec{j}_y^*$ and $\bar{q}_y^*$ to replace $j_y^*$ and $q_y^*$.

Other forms of forcing term accounted for the discrete effect could also be adopted [17]. It is noted that the compressibility may influence the results, and this can be eliminated by incompressible model

[20]. However, since the present Mach number is low, therefore this influence could be neglected [35].

## 2.4. Geometry and boundary conditions

Natural convection in a square cavity $\Omega = ]0, H[ \times ]0, H[$ (see Fig. 1) is considered, where the flow is bounded by a stationary square enclosure with sidewalls maintained at different temperatures and driven by the buoyancy force. For laminar convection in this flow configuration, the viscous heat dissipation is assumed to be negligible. The temperature difference between the walls introduces a temperature gradient in the fluid, and the consequent density difference induces a convective fluid motion. The left wall is at the higher uniform temperature $T_l$ and the right wall is at the lower uniform temperature $T_r$. Both the top and bottom walls are adiabatic, i.e. $\partial T / \partial y = 0$. The summary of the boundary conditions is shown below.

$$ u = v = 0 \text{ on } \partial \Omega \qquad (13) $$
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Second component of velocity (V) of flow for Pr=0.71. From top to bottom Ra = 10³, 10⁴, 10⁵, and 10⁶ respectively for mesh size 105², 155², 205² and 255². From left to right: (left) $κ_x = κ_y/2$, (center) isotropic case $κ_x = κ_y$, (right) $κ_x = 2κ_y$.</figcaption></figure>

$$T = T_l \text{ on } \{0\} \times [0, H] \qquad (14)$$

$$T = T_r \text{ on } \{H\} \times [0, H] \qquad (15)$$

$$\frac{\partial T}{\partial y} = 0 \text{ on } [0, H] \times \{0\} \text{ and } [0, H] \times \{H\} \qquad (16)$$

For Dirichlet boundary condition for the velocity (13) at the walls of the cavity, the classical half way bounce-back boundary condition is adopted. So, for example, consider the bottom wall for a boundary node $x_b$ (see left figure of Fig. 2), the following bounce-back boundary condition is applied.

$$f_2(x_b, t + \Delta t) = f_4(x_e, t + \Delta t) = f_4^*(x_b, t),$$

$$f_5(x_b, t + \Delta t) = f_7(x_c, t + \Delta t) = f_7^*(x_b, t),$$

$$f_6(x_b, t + \Delta t) = f_8(x_d, t + \Delta t) = f_8^*(x_b, t).$$

For the thermal boundary condition, the Dirichlet boundary conditions given by Eqs. (15) and (14) on the left and right wall of the domain $Ω$ are introduced. For a given constant temperature $T$, this

can be archived using the following scheme in boundary node $x_b$ on the right wall (see right figure of Fig. 2):

$$g_3(x_b, t + \Delta t) = -g_1(x_e, t + \Delta t) + \frac{1}{36}(4 - \tilde{\alpha} - 2\tilde{\beta})T,$$

$$g_7(x_b, t + \Delta t) = -g_5(x_c, t + \Delta t) + \frac{1}{36}(4 + 2\tilde{\alpha} + \tilde{\beta})T,$$

$$g_6(x_b, t + \Delta t) = -g_8(x_d, t + \Delta t) + \frac{1}{36}(4 + 2\tilde{\alpha} + \tilde{\beta})T,$$

For the Neumann boundary condition on the top and bottom wall of the domain $Ω$ given by Eq. (16), the classical “bounce back” scheme is adopted. Consider a boundary node $x_b$ in the bottom wall (see right figure of Fig. 2), the following scheme is used.

$$g_2(x_b, t + \Delta t) = g_4(x_e, t + \Delta t) = g_4^*(x_b, t),$$

$$g_5(x_b, t + \Delta t) = g_7(x_c, t + \Delta t) = g_7^*(x_b, t),$$

$$g_6(x_b, t + \Delta t) = g_8(x_d, t + \Delta t) = g_8^*(x_b, t).$$

For more detail about how to reconstruct the above boundary condition for thermal problem see [7].
<figure><img src="image_7.png" /><figcaption>Fig. 6. Streamlines of flow fields (Ψ) for Pr=0.71. From top to bottom Ra = 10³, 10⁴, 10⁵, and 10⁶ respectively for mesh size 105², 155², 205² and 255². From left to right : (left) κx = κy / 2, (center) isotropic case κx = κy, (right) κx = 2κy.</figcaption></figure>

In the present parallel implementation, the single program multiple data (SPMD) environment is employed. Message-Passing-Interface (MPI) is adopted for the communication between the processors. The domain decomposition is done on direction of the computational domain, where the ghost cells are adopted along the inter-processor boundary.

The Nusselt number is also an important dimensionless parameter in describing the convective heat transport. Its average in the whole flow domain is defined as,

$$ \overline{Nu} = \frac{1}{\kappa \Delta T} \int_{0}^{H} q_{x}(x, y) dy \qquad (17) $$

where $q_x(x, y) = uT(x, y) - \kappa \partial T(x, y) / \partial x$ is the local heat flux in the horizontal direction.

To compare with previous results, the main quantities to compute are : $u_{max}$, y, $v_{max}$, x and $\overline{Nu}$. Where $u_{max}$ and its location y, the maximum vertical velocity on the horizontal mid-plane of the cavity, $v_{max}$ and its location x, the maximum horizontal velocity on the vertical mid-plane of the cavity and the average Nusselt number $\overline{Nu}$.

We compute, for some cases, the maximum stream function $\psi_{max}$ on the whole domain. Where the stream function is determined from:

$$ \nabla \cdot (\nabla \times \psi) = \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y}. $$

## 3. Numerical results and discussion

### 3.1. Isotropic case

Let consider the isotropic thermal case, i. e. the x- thermal diffusivity $κx$ is equal to y- thermal diffusivity $κy$ equal to a given thermal diffusivity $κ$. In this case we fix $σ_3 = σ_5$ to have $κ_x = κ_y$, described by Eq. (11), equal to the given thermal diffusivity $κ$. For the present natural convection within the square cavity as shown in Fig. 1, the major control parameter is the Rayleigh number $Ra = \beta g \Delta T H^3 Pr / \nu^2$ associated with the heat transfer within the fluid, where H is the height or width of the cavity.
<figure><img src="image_8.png" /><figcaption>Fig. 7. Isotherms of flow fields (T) for Pr=0.71. From top to bottom Ra = 10³, 10⁴, 10⁵, and 10⁶ respectively for mesh size 105², 155², 205² and 255². From left to right : (left) κx = κy/2, (center) isotropic case κx = κy, (right) κx = 2κy.</figcaption></figure>

Note here that for the computation of $\overline{Nu}$ described by Eq. (17), the temperature gradient $\frac{\partial T}{\partial x}$ is needed. To calculate this gradient we do not do any additional interpolation method. In fact this quantity can be evaluated by using Taylor expansion [8] up to second-order of the non-conserved moment $\tilde{m}_1$:

$$ \tilde{m}_1 = -\lambda^2 \Delta t \frac{1}{\sigma_{\kappa}} \left[ \frac{4 + \tilde{\alpha}}{6} \frac{\partial T}{\partial x} \right] + O(\Delta t^2). $$

All the velocities are normalized using the diffusion velocity $\kappa/H$. The temperature are dimensionless, locations x and y are normalized using H.

## 3.2. Stability and admissible grid

Let define the Mach number as follows:

$$ Ma = \frac{U}{C_s}, \qquad (18) $$

where the quantity $U = \sqrt{\alpha g \Delta T H} = \sqrt{\frac{Ra}{Pr}} \frac{\nu}{H}$ is the characteristic velocity in thermal convective flows. So, to keep stability of the LBM scheme related dynamic field, the Mach number should be smaller than critical value $\frac{3}{10}$ (see [38] for more details). As in numerical simulation the parameters are fixed as $Pr = 0.71$ and $\nu = 0.01$, Eq. (18) gives a constrain on mesh size $H$. In fact, the mesh size must verify

$$ \frac{10}{3} \frac{\nu}{C_s} \sqrt{\frac{Ra}{Pr}} < H. \qquad (19) $$

Example for $Ra = 10^7$ and $Ra = 10^8$ the mesh size must satisfy $H > 216$ and $H > 685$, respectively.

## 3.3. Grid dependence

We begin by the study of the grid dependence and the accuracy of the scheme. In fact, Table 1 gives the results for the simulation for Rayleigh number $Ra = 10^6$ by using some different mesh sizes. We note that the calculated values approach the values given by the
benchmark of de Vahl Davis [4] and the benchmark of Le Quéré [25]. In other hand table shows also second order accuracy for every measured quantities. Where the accuracy is calculated by least squares method using relative error between the solution obtained by mesh size $N^2$ and the reference one obtained by fine mesh ($N_{ref}^2 = 1001^2$).

We have also see the grid dependance in the case of Ra = $10^8$ for the following mesh size : 495 × 495, 987 × 987 and 2016 × 2016. The Table 2 shows that the quantities calculated quickly approach the values given by the benchmark of the Le Quéré [25]. We note here that the first grid size do not resolve the stability condition Eq. (19). But this size still gives a good solution.

### 3.4. Isotropic test case

Many papers study the square heated cavity, for the following 4 values of the Ra numbers : Ra = $10^3$, $10^4$, $10^5$ and $10^6$. So we compare our results to the following results: de Vahl Davis [4] benchmark solution where second order finite differences scheme and a Richardson extrapolation scheme are used, Mayne et al. [31] using h-adaptive finite element method, Kuznik et al. [21] and Liu et al. [28] using TLBM based on the BGK and Mezrhab et al. [33] and Wang et al. [38] using TLBM based on MRT D2Q9 for flow and MRT D2Q5 for temperature.

Table 3 shows the numerical results where the domain is covered by a lattice sizes of 105 × 105, 155 × 155, 205 × 205 and 255 × 255, respectively for Ra = $10^3$, $10^4$, $10^5$ and $10^6$ compared to the results obtained by the methods listed above. The simulated results are contrasted with the benchmark solutions of De Vahl Davis [4] and the agreements are satisfactory. It is also noted that differences of the predicted velocities and average Nusselt number are less than 0.1%. The middle column of Figs. 3–7 show the solution predicted by the present double D2Q9 MRT LBE method for Rayleigh numbers Ra = $10^3$, for Ra = $10^4$, Ra = $10^5$ and Ra = $10^6$ for isotropic case.

### 3.5. Anisotropic test case

In this section the effect of the anisotropy is performed. Let $κ_x$ the x thermal diffusivity and $κ_y$ the y thermal diffusivity. So we define the x Rayleigh number $Ra_x = \beta g \Delta T^3 \kappa_x / \nu$ associated with the x thermal diffusivity and $Ra_y = \beta g \Delta T^3 \kappa_y / \nu$ associated with the y thermal diffusivity. The choice of the anisotropy will be for different Rayleigh numbers fixed $Ra_y = 10^3, \dots, 10^6$ as follows:

* x thermal diffusivity given by $κ_x = \frac{κ_y}{2}$.

* x thermal diffusivity given by $κ_x = \frac{κ_y}{2}$.

The average of Nusselt number in the whole flow domain is defined now as,

$$ \overline{Nu} = \frac{1}{\kappa_x \Delta T} \int_0^H q_x(x, y) dy \quad (20) $$

where $q_x(x, y) = uT(x, y) - \kappa_x \partial T(x, y) / \partial x$ is the local heat flux in the horizontal direction.

To compare with isotropic results, the same main quantities to compute are : $u_{max}$, $y_{max}$, $x_{max}$ and $\overline{Nu}$. We remark that all the velocities are normalized using the y diffusion velocity $κ_y/H$.

First we study the convergence and the accuracy of the scheme for anisotropic case. In fact, Tables 4 and 5 give the results for the simulation in case of Rayleigh number $Ra_y = 10^6$, for $κ_x = \frac{κ_y}{2}$ and $κ_x = 2κ_y$, respectively, by using some different mesh sizes. Tables 4 and 5 show also second order accuracy for every measured quantities. Here the accuracy is calculated by least squares method using relative error between the solution obtained by mesh size $N^2$ and the reference one obtained by fine mesh ($N_{ref}^2 = 1001^2$).

Table 6 shows the numerical results for three different cases: $κ_x = \frac{κ_y}{2}$, $κ_x = κ_y$ (isotropic case) and $κ_x = 2κ_y$. The domain is covered by a lattice sizes of 105 × 105, 155 × 155, 205 × 205 and 255 × 255, respectively for $Ra_y = 10^3$, $10^4$, $10^5$ and $10^6$.

We note here that the $u_{max}$ and $v_{max}$ increase when the $κ_x$ increase. This is due to the fact of the imposed hot wall (at $x = 0$) and cold wall (at $x = 1$) is in x direction.

Figs. 3–7 show the solution predicted by the present double D2Q9 MRT LBE method for Rayleigh numbers $Ra_y = 10^3$, $10^4$, $10^5$ and $10^6$ for anisotropic case. We refind here that the effect of x thermal diffusivity $κ_x$ is more important than y thermal diffusivity $κ_y$. In fact when $κ_x = 2κ_y$ the velocity of the fluid is bigger.

## 4. Conclusion

In this paper, a multi-relaxation time thermal lattice Boltzmann scheme has been applied to compute natural convection flow within differential heated square cavity. For Rayleigh number under $10^6$ the present results compare favorably with previous benchmark solutions. Then anisotropic thermal diffusion is investigated ($κ_x = \frac{κ_y}{2}$ and $κ_x = 2κ_y$). The solution is compared to the isotropic case. We note the ability of double D2Q9 population to resolve anisotropic problem. Finally we remark when the x thermal diffusivity $κ_x$ increase the velocity of the fluid increase and the convergence of the scheme to the steady state is faster.

## Acknowledgments

The authors thank Pierre Lallemand (Beijing Computational Science Research Center, Beijing, China) for helpful discussion during the elaboration of this work. Many thanks also to the support by the Taiwan National Science Council (Grant 94-2212-E-007-059) and the computational facilities provided by the Taiwan National Center for High-Performance Computing.

## References

[1] Alexander FJ, Chen S, Sterling JD. Lattice Boltzmann thermodynamics. Phys Rev E 1993;47:2249.

[2] Bespalko D, Pollard A, Uddin M. Analysis of the pressure fluctuations from an LBM simulation of turbulent channel flow. Comput Fluids 2012;54:143.

[3] Contrino D, Lallemand P, Asinaria P, Luo L. Lattice-Boltzmann simulations of the thermally driven 2D square cavity at high Rayleigh numbers. J Comput Phys 2014;275:257-72.

[4] De Vahl Davis G. Natural convection of air in a square cavity: a bench mark numerical solution. Int J Numer Methods Fluids 1983;3:249.

[5] Dubois F. Equivalent partial differential equations of a lattice Boltzmann scheme. Comput Math Appl 2008;55:1141-9.

[6] Dubois F, Lallemand P. Towards higher order lattice Boltzmann schemes. J Stat Mech: Theory Exp 2009;2009:06006.

[7] Dubois F, Lallemand P, Tekitek MM. Using lattice Boltzmann scheme for anisotropic diffusion. Finite Volumes for Complex Applications V, 795. London: ISTE, Wiley; 2008.

[8] Dubois F, Lallemand P, Tekitek MM. On a superconvergent lattice Boltzmann boundary scheme. Comput Math Appl 2010;59:2141.

[9] Dubois F, Lallemand P, Tekitek MM. Taylor Expansion Method for Linear Lattice Boltzmann Schemes with an External Force: Application to Boundary Conditions. Lecture Notes in Computational Science and Engineering, 99. European Workshop HONOM 2013. Bordeaux, France, March 18-22, Springer: 2014. p. 89-107.

[10] Dixit HN, Babu V. Simulation of high Rayleigh number natural convection in a square cavity using the lattice Boltzmann method. Int J Heat Mass Transf 2006;49:727.

[11] Flückiger R, Freunberger SA, Kramer D, Wokaun A, Scherer GG, Büchi FN. Anisotropic, effective diffusivity of porous gas diffusion layer materials for PEFC. Electrochim Acta 2008;54:551.

[12] Frisch U, d'Humières D, Hasslacher B, Lallemand P, Pomeau Y, Rivet J-P. Lattice gas hydrodynamics in two and three dimensions. Complex Syst 1987;1:649-707.

[13] Ginzburg I. Equilibrium-type and link-type lattice Boltzmann models for generic advection and anisotropic-dispersion equation. Adv Water Resour 2005;28:1171-95.

[14] Ginzburg I, Adler P. Boundary flow condition analysis for the three-dimensional lattice Boltzmann model. J Phys II 1994;4:191-214.

[15] Ginzburg I, d'Humières D. Lattice Boltzmann and analytical modeling of flow processes in anisotropic and heterogeneous stratified aquifers. Adv Water Resour 2007;30:2202-34.

[16] Guo Z, Shi B, Zheng C. A coupled lattice BGK model for the Boussinesq equations. Int J Numer Methods Fluids 2002;39:325.

[17] Guo ZL, Zheng CG, Shi BC. Discrete lattice effects on the forcing term in the lattice Boltzmann method. Phys Rev E 2002;65:1. Article Number: 046308.
[18] Guo Z, Zheng C, Shi B. Thermal lattice Boltzmann equation for low Mach number flows: decoupling model. Phys Rev E 2007;75:15.

[19] He X, Chen S, Doolen GD. A novel thermal model for the lattice Boltzmann method in incompressible limit. J Comput Phys 1998;146:282.

[20] He XY, Luo LS. Lattice Boltzmann model for the incompressible Navier-Stokes equation. J Stat Phys 1997;88:927.

[21] Kuznik F, Vareilles J, Rusaouen G, G Krauss G. A double-population lattice Boltzmann method with non-uniform mesh for the simulation of natural convection in a square cavity. Int J Heat Fluid Flow 2007;28:862.

[22] Lallemand P. Improve LBE Model, Summer School on the Lattice Boltzmann Method. Beijing Computational Science Research Center; 2012.

[23] Lallemand P, Luo L. Theory of the lattice Boltzmann method: dispersion, dissipation, isotropy, Galilean invariance, and stability. Phys Rev E 2000;61:6546-62.

[24] Lallemand P, Luo L-S. Hybrid finite-difference thermal lattice Boltzmann equation. Int J Mod Phys B 2003;17:41.

[25] Quéré PL. Accurate solutions to the square differentially heated cavity at high Rayleigh number. Comput Fluids 1991;20:19.

[26] Lin L-S, Chen Y-C, Lin C-A. Multi relaxation time lattice Boltzmann simulations of deep lid driven cavity flows at different aspect ratios. Comput Fluids 2011;45:233.

[27] Lin K-H, Liao C-C, Lien S-Y, Lin C-A. Thermal lattice Boltzmann simulations of natural convection with complex geometry. Comput Fluids 2012;69:35.

[28] Liu CH, Lin KH, Mai HC, Lin CA. Thermal boundary conditions for thermal lattice Boltzmann simulations. Comput Math Appl 2010;59:2178.

[29] Mai H-C, Lin K-H, Yang C-H, Lin C-A. A thermal lattice Boltzmann model for flows with viscous heat dissipation. Comput Model Eng Sci 2010;61:45.

[30] Massaioli F, Benzi R, Succi S. Exponential tails in two-dimensional Rayleigh-Benard convection. Europhys Lett 1993;21:305.

[31] Mayne DA, Usmani A, Crapper M. h-Adaptive finite element solution of unsteady thermally driven cavity problem. Int J Numer Methods Heat Fluid Flow 2000;10:598.

[32] McNamara G, Garcia AL, Alder BJ. Stabilization of thermal lattice Boltzmann models. J Stat Phys 1995;81:395.

[33] Mezrhab A, Moussaoui MA, Jami M, Naji H, Bouzidi M. Double MRT thermal lattice Boltzmann method for simulation convective flows. Phys Lett A 2010;374:3499.

[34] Pasaogullari U, Mukherjee PP, Wang C-Y, Chen KS. Anisotropic heat and water transport in a PEFC cathode gas diffusion layer. J Electrochem Soc 2007;154:823.

[35] Peng Y, Shu C, Chew YT. Simplified thermal lattice Boltzmann model for incompressible thermal flows. Phys Rev E 2003;68:026701.

[36] Shi Y, Zhao TS, Guo ZL. Thermal lattice Bhatnagar-Gross-Krook model for flows with viscous heat dissipation in the incompressible limit. Phys Rev E 2004;70:066310.

[37] Shih C-H, Wu C-L, Chang L-C, Lin C-A. Lattice Boltzmann simulations of incompressible liquid-gas systems on partial wetting surfaces. Philos Trans R Soc A - Math Phys Eng Sci 2011;369:2510.

[38] Wang J, Wang DH, Lallemand P, Luo L-S. Lattice Boltzmann simulations of thermal convective flows in two dimensions. Comput Math Appl 2013;65:262-86.

[39] Yu DZ, Mei RW, Luo L-S, Shyy W. Viscous flow computations with the method of lattice Boltzmann equation. Prog Aeosp Sci 2003;39:329.