1 Introduction

Natural convection in porous cavities has attracted many researchers due to its occurrence in various engineering applications, e.g., geothermal energy extraction, chemical catalytic reactors, electrical cooling, nuclear reactor systems, and building and thermal insulation systems^[1-4]. It has been shown in the literature that the natural convection may be affected by magnetic field, and the hydrodynamic flow and heat transfer in different porous and non-porous media have been investigated^[5-15].

The natural convection in porous media has been studied by many researches with both experimental and numerical methods. Pangrle et al.^[16] carried out an experimental study on the magnetic resonance imaging of the laminar fluid flow in porous tubes and shell systems. Sathiyamoorthy et al.^[17] numerically studied the natural convection flow in a square cavity filled with porous media. Sheikhzadeh et al.^[18] examined the steady state natural convection of the nanofluid inside a cavity filled with a porous medium numerically with the finite volume approach.

The lattice Boltzmann method (LBM) evolving from the lattice-gas automata method^[19] is a powerful approach for the simulation of complex fluid flow and complex physics^[20-24]. The existing lattice Boltzmann (LB) models for the flow in porous media (with or without magnetic field) can be generally classified into two divisions, i.e., the pore scale method^{[20, 25-27]} and the representative elementary volume (REV) scale method^[27-29]. In the pore scale method, the standard LB model is used to simulate the fluid flow in the pores, and the interaction between the solid phase and the fluid phase is accomplished with the no-slip bounce-back rule. With this method, one can obtain the detailed flow information of the pores. However, in this method, the detailed geometric information of the pores are needed, and each pore needs several lattice nodes in the simulation. Therefore, this method is very demanding and time-consuming for the computer performance, and is usually used for small domains^[28]. In the REV scale LB method^[27-30], an additional term is added to the standard LB equation in order to consider the effect of the porous media based on some semi-empirical models, e.g., the Darcy model, the Brinkman-extended Darcy model, and the Forchheimer-extended Darcy model. The main advantage of the REV-scale LB model is that it can be used in practical engineering applications with extensive computations^[31]. Seta et al.^[32] used the REV scale LB method to simulate the natural convection in a porous media. Mehrizi et al.^[33] simulated the force convection in a vented cavity filled with porous media by use of the LBM. Ashorynejad et al.^[34] studied the natural convection in a porous cavity in the presence of a vertical magnetic field with the LBM, and found that the magnetic field would reduce the heat transfer and fluid circulation within the cavity. Liu and He^[35] developed a non-orthogonal multiple-relaxation-time LB (MRT-LB) method to simulate the natural convection heat transfer in porous media at the REV scale based on the generalized non-Darcy model.

However, based on authors' knowledge, there is still no research on the effect of uniform magnetic field on the porous cavity with sinusoidal temperature distribution. Therefore, in this paper, we aim to use the LBM to simulate the MHD natural convection in a porous cavity with sinusoidal temperature distribution. The mathematical formulation for porous media is based on the Brinkman-Forchheimer equation model^[28]. The numerical results are presented in the form of streamlines, isotherms, and the Nusselt number, and the effects of different parameters, e.g., the Hartmann number, the porosity, the Darcy number, and the phase deviation, are studied.

2 Physical problem and macroscopic governing equations

Consider a two-dimensional (2D) natural convection flow in a square cavity with porous media filled with an electrically conducting fluid (water). The geometry of the problem is presented in Fig. 1.

The top and bottom walls of the cavity are thermally insulated, and the left and right walls are heated sinusoidally. A uniform magnetic field is exerted in the horizontal direction. In order to obtain the governing equations, a continuous model based on the REV scale is created for the porous media. A Cartesian reference is defined, and an element with an appropriate size (very larger than the volume of the porous cells) is examined to obtain the valid averaging. In addition, the following assumptions are considered:

(ⅰ) The fluid flow is Newtonian, steady, and laminar.

(ⅱ) The physical properties of water are considered to be constant except the variation of the density, which is approximated by the Boussinesq approximation.

(ⅲ) The magnetic Reynolds number is assumed to be very small, so that the included magnetic field and Hall effect can be negligible^[36].

(ⅳ) The porous media are assumed to be in the local thermal equilibrium state with the solid matrix and also be homogeneous and isotropic.

(ⅴ) The source term and radiation in the energy equation are negligible.

Therefore, the continuity, momentum, and energy equations for water can be written as follows^[37-39]:

(1)

(2)

(3)

where u is the velocity vector, T is the temperature, p is the pressure, ρ is the density, μ is the dynamic viscosity, α is the thermal diffusivity, and ε is the porosity. f is defined by^[39]

(4)

where j and v are the normal unit vector and velocity component in the y-direction, respectively. υ is the kinematic viscosity, g is the gravitational acceleration, β is the volume expansion coefficient, and T_m is the mean temperature. The Hartman number Ha is defined by Ha=, where σ is the electrical conductivity. K is the permeability, and is defined by

(5)

where Da is the Darcy number, and L is the characteristic length.

Introduce the following dimensionless variables:

(6)

The governing equations (1)-(4) can be written in dimensionless form as follows:

(7)

(8)

(9)

(10)

In the above equations, T_h and T_c are the temperatures of the hot and cold wall, respectively, and

3 LB model

In the present contribution, the standard D2Q9 model (see Fig. 2) is used for the flow and the temperature. The nine discrete velocities are^[40]

(11)

where c=Δx/Δt is the lattice velocity, Δx is the lattice spacing, Δt is the time step, ω_i (i=0, 1, 2, ..., 8) are the D2Q9 lattice constants, and

The Boltzmann equation used to solve the flow field is given by^[40]

(12)

where f_i is the density distribution function, and F_i is the external force term. It is demonstrated that^[32] the suitable choice for the force term is

(13)

In the above equation, the body force F is defined by Eq. (4).

τ_v is the relaxation time for the velocity field, and is given by

(14)

where c_s is the lattice speed of sound, and .

The equilibrium density distribution function f_i^eq for the D2Q9 model is defined by

(15)

The velocity vector u in the above equation is calculated by using a temporal velocity v as follows^[45]:

(16)

(17)

(18)

where K is given by Eq. (5).

The Boltzmann equation for solving the temperature field is

(19)

where g_i is the temperature distribution function, and τ_c is the relaxation time for the temperature field defined by

(20)

g_i^eq is the equilibrium distribution function for the temperature field, and is defined by

(21)

u and T can be calculated by the mentioned variables with

(22)

3.1 Boundary conditions

The bounce-back boundary conditions are applied on all solid boundaries, which means that the incoming boundary populations are equal to the out-going populations after the collision. Therefore, we impose the following conditions:

(ⅰ) The north boundary

(23)

(ⅱ) The south boundary

(24)

(ⅲ) The east boundary

(25)

(ⅳ) The west boundary

(26)

The adiabatic boundary condition is used on the north and south boundaries. For the north boundary, the following conditions are required:

(27)

For the south boundary,

(28)

For the east boundary, the east wall temperature is defined by , and

(29)

For the west boundary, the west wall temperature is defined by , and

(30)

4 Non-dimensional parameters

To ensure that the code works in near incompressible regime, the Mach number Ma should be less than 0.3. In this work, the Mach number is selected to be 0.1. By fixing the Rayleigh number, the Prandtl number, and the Mach number, the viscosity can be calculated by

(31)

where n is the lattice number in the y-direction, and the speed of sound is constant, i.e., . The Rayleigh number and the Prandtl number are defined by

(32)

The local Nusselt number on the left and right walls and the average Nusselt number are calculated as follows:

(33)

(34)

5 Grid testing and code validation

For the independency of grids, the variations of the average Nusselt number for the case, where

are studied for five different lattice sizes (see Fig. 3). It can be seen in Fig. 3 that the differences among the results of the meshes of 80× 80, 100× 100, and 120× 120 are very small. Therefore, the grid size of 100× 100 is chosen in this work.

To check the accuracy of the present numerical simulation, the present code is validated with the results of Ashorynejad et al.^[34]. As can be seen in Fig. 4, there is good agreement between the numerical results of the present study and the results of Ashorynejad et al.^[34].

Table 1 presents the results of the current study in comparison with those of Mejri et al.^[41]. It is seen that the two results agree well with each other.

6 Results and discussion

Figures 5 and 6 show the effects of Ha and the phase derivation γ on the streamlines and isotherms, respectively, where ε =0.6, and Da=10^-2. Figure 5 shows that, when Ha increases, the fluid flow weakens, leading to the symmetrization of the streamline cells at γ =0°. When γ increases, the four created cells transform to two cells, which will form identical shape at γ =180°. At γ =90°, a small vortex is created in the right corner of the cavity. When γ increases, this vortex disappears completely. Figure 6 shows that, for all γ and Ha, the isotherms are retained along the left sidewall. The heat transfer is constant on the left sidewall but various on the right sidewall. When γ increases, the heat transfer behavior of the left sidewall will change (see Fig. 6). It also shows that the thickness of the boundary layers increases when Ha increases, and thus the compactness of the isotherms near the walls decreases.

Figures 7 and 8 illustrate the effects of Da on the streamlines and isotherms for different Ha, where ε =0.9, and γ =90°. As can be seen in the streamlines, for all Ha, when Da increases, the shape of the streamlines changes. When Da=10^-2, there are two main cells in the streamlines. When Da increases, each of the two main cells tends to convert to two cells. This is because that decreasing Da will reduce the permeability and thus weaken the flow strength. It can also be seen from the isotherms that, the thickness of the boundary layer increases when Da decreases. This is due to the reduction of the permeability. When Ha increases, the thickness of boundary layer increases, similar to that when Da decreases.

Figure 9 presents the effect of the porosity on the sidewalls when Ha=0, Da=10^-2, and γ=90°. The curves drawn for the Nusselt number against y/L are approximately of sinusoidal shape like the thermal boundary. It is seen that, decreasing the porosity will reduce the local Nusselt number on the sidewalls. This is because that, when the porosity decreases, the porous media will be converted to the solid media. Therefore, it is expected that the heat transfer mechanism will change to conduction and thus the Nusselt number will decrease. From the other point of view, increasing the porosity will make the volume of the holes in the porous media larger, and thus the heat transfer mechanism will be convection and the Nusselt number will increase.

Figure 10 illustrates the local Nusselt number on the sidewalls for different γ, where Ha=30, Da=10^-2, and ε =0.6. The local Nusselt number curves are approximately of sinusoidal shape like the thermal boundary along the sidewalls. This clearly shows that the local heat transfer is directly affected by the temperature distribution on the surface. The maximum value of the local Nusselt number on the left wall is obtained at γ =180°. But the local Nusselt number along the right side wall is greatly affected by the change of the phase deviation. The effects of Ha on the local Nusselt number on both sidewalls are presented in Fig. 11, where γ =90°, Da=10^-2, and ε =0.6. Since the magnetic field has the tendency to slow down the movement of the fluid in the cavity, increasing Ha will weaken the natural convection and thus reduce the local Nusselt number on both sidewalls.

Figure 12 depicts the effects of Da on the local Nusselt number of both sidewalls, where γ =90°, Ha=0, and ε =0.6. It can be seen that, decreasing Da will obviously reduce the local Nusselt number on both sidewalls. This is due to the fact that, reducing Da will slow down the fluid flow and thus weaken the natural convection.

Table 2 presents the effects of the porosity on the average Nusselt number for different Da and γ, where Ha=0. At Da=10^-2, when the porosity increases, the average Nusselt number increases for all of the phase derivation. However, at Da=10^-4, the permeability is very small, and the effect of increasing the porosity on the average Nusselt number is not noticeable. The maximum and minimum average Nusselt numbers are obtained at γ =90° and γ =0°, respectively. Moreover, decreasing Da weakens the natural convection and reduces the average Nusselt number.

Figure 13 shows the effects of the porosity and Ha on the average Nusselt number, where Da=10^-2, and γ =90°. It can be seen that, increasing the porosity and Ha will decrease the average Nusselt number. When Ha increases, the natural convection will be weakened and thus the local Nusselt number will decrease. However, increasing the porosity will enhance the average Nusselt number.

7 Conclusions

In the present work, the MHD natural convection in a porous cavity with sinusoidal temperature distribution is studied with the LBM. The main conclusions are as follows:

(ⅰ) The local heat transfer is directly affected by the temperature distribution on both sidewalls, i.e., larger heat transfer occurs when the temperature is higher.

(ⅱ) At γ=90°, heat transfer and fluid flow will decline with an increase in Ha while enhance with an increase in the porosity.

(ⅲ) When Ha=0 and Da=10^-2, the average Nusselt number will be creased by increasing the porosity for all γ. However, it has a fixed value when Da=10^-4. The maximum and minimum average Nusselt numbers are obtained at γ=90° and γ =0°, respectively.

(ⅳ) Decreasing Da has a noticeable effect on the local Nusselt number on both sidewalls. When Da decreases, the local Nusselt number on both sidewalls declines.