# Double-Population Link-Wise Artificial Compressibility Method for Simulating Two-Dimensional Magnetoconvection

Authors: Mohamed Hamdi; Souheil Elalimi; Sassi Ben Nasrallah
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

The recently developed Link-Wise Artificial Compressibility Method (LW-ACM) with double population is presented, validated, and then applied to magnetoconvective flow. The validation step evolves in two stages. First of all, we used the problem of mixed convection over a vertical duct preceded by a sudden expansion in order to validate the multiple relaxation times (MRT) code. Then, we explore the accuracy and the performance of the LW-ACM versus the MRT model and the classic Bhatnagar-Gross-Krook (BGK) model by using a differentially heated cavity problem for Prandtl numberP
r
=
0.71. The numerical results for values of Nusselt numbers and process time clearly showed the accuracy and the efficiency of LW-ACM. The adequate accuracy of the proposed model encouraged us to apply it to a magnetoconvection problem. Dynamic and thermal fields are presented and analyzed at P
r
=
0.054 and Rayleigh number R
a
=
5
×
10
4 for Hartmann number H
a up to 60. Results revealed that the effects of magnetic fields on dynamic and thermal fields are important and that the heat transfer in the cavity is a decreasing function of magnetic field strength.


---

## Body

## 1. Introduction

The Lattice Boltzmann Method, commonly called LBM, was introduced as a new alternative for the numerical simulation of physical phenomena that can address many problems in physics [1]. Lattice Boltzmann Method sustains today a rapid evolution in terms of physical models, computer implementation, and engineering applications. This method gains an increasing attention and becomes a convenient alternative to classical Computational Fluid Dynamics (CFD) methods. In fact, authors of LBM exceeded the validation step of the method for a wide range of physical problems and many of them are oriented to make the method more efficient than other CFD methods. However, more studies are still needed to make the method useful in engineering practice. The evolutionary process and the formulation of LBM can be highlighted as the following: The Lattice Gas Cellular Automat (LGCA), the continuous Boltzmann-BGK equation, and then the Grad’s Hermite-quadrature expansion. Since its development, LBM demonstrated a successful progress not only in solving viscous flow problems but also in heat transfer area [2]. The first connection of the lattice Boltzmann equation (LBE) to the Boltzmann equation has been established by He and Luo [3]. These formulations can construct models that recover incompressible Navier-Stokes Equations. The Bhatnagar-Gross-Krook (BGK) approximation [4] is the most popular lattice Boltzmann model. This approach is derived from the Enskog equation. In this model the collision operator requires the same relaxation time for each physical quantity. In spite of its simplicity, BGK model suffers from numerical instability at high Reynolds (or Rayleigh number) number. To avoid this restriction, the direct way is to use a large number of grid points; however this will cost large computer resources and lower the computational efficiency. The instability problems of LBM have been extensively studied by many authors through theoretical and numerical investigations. Thus, care was taken to justify the techniques used to stabilize the method without affecting its accuracy.To remove numerical instability in LBM, some authors [2] used the multiple relaxation times model (MRT). The advantage of this model is that it has adjustable parameters with an additional degree of freedom. These parameters can be determined by optimizing the hydrodynamic properties of the model and linear stability analysis of the LBE evolution operator. Another way to remedy the stability problem in LBM is to return into the entropic formulation of LBM. Some authors [5] used entropic lattice Boltzmann schemes. The derivation of ELB can be performed in many ways and the most popular one is derived from the analog of the discrete Boltzmann H function of standard statistical mechanics.The artificial compressibility method (ACM) introduced by Chorin in 1967 [6] to solve the incompressible Navier-Stokes Equations (NSE) is recently developed with an easy formulation known as Link-Wise Artificial Compressibility Method (LW-ACM) that benefits from some similarities between its classic formulation and Lattice Boltzmann Methods (LBM). Asinari et al. [7] developed the LW-ACM by a finite set of discrete directions (links) on a regular lattice mesh similar to LBM and demonstrated the stability and accuracy of the proposed model. It should be mentioned that the link between LBM and ACM was observed earlier by He et al. [8] and they performed a comparison between the two methods. Their results showed that LBM and ACM are closely related to each other. The similarities between them appear in the continuous form of the macroscopic governing equations, while they differ from each other in their discrete forms. Also they found that difference between LBM and ACM does not affect the momentum results, but it has an impact on the pressure fields. Our recent paper [9] reports numerical results for attached and separated flows by using LW-ACM. The numerical results are relevant for exploiting the applicability of the method to these setups and the computational efficiency of such model is substantiated.In fact, LW-ACM can be used to solve convective flows. Obrecht and Kuznik [10] used a hybrid thermal LW-ACM to solve the problem of a differentially heated cubic cavity. Their results showed that this scheme remains stable until Rayleigh number R
a
=
10
8. In fact, LW-ACM can be used with double-population model and simulate convective flows and save the memory usage compared to MRT-LBM [11, 12]. Compared to the MRT implementation, the main advantage of the LW-ACM formulation is its simplicity and suitability for parallel implementation. The stability of LBM simulations continues to be the subject of current studies. Recent papers dealing with LBM are oriented towards its efficiency against classical CFD solvers. Thus, much attention has been accorded to the performance and the efficiency of the implementation. The present study deals with the accuracy and the efficiency of the LW-ACM model and then applying it to a magnetoconvective flow problem. The remainder of the present paper is organized as follows. In Section 2, we describe the problem under investigation. The LBM equations with double-population LW-ACM formulation are presented in Section 3. In Section 4, the prescribed LBM model is validated and compared to BGK and MRT results. Then the LW-ACM code is used to compute the magnetoconvection flow in a differentially heated cavity. In particular, the effects of magnetic strength are discussed.


## 2. Problem Formulation and Boundary Conditions

The problem under investigation is sketched in Figure1 and detailed as the following. A uniform and horizontal magnetic field B is applied to a square cavity filled with an electrically conducting, viscous, and incompressible fluid with negligible viscous dissipation and radiation effects. The cavity of height H is heated from below and cooled from above with constant temperatures T
h and T
c, respectively T
h
>
T
c. The upper and lower walls generate a vertical temperature gradient. The problem is considered as two-dimensional. Using the Boussinesq approximation, the fluid density varies as ρ
=
ρ
0
1
-
β
T
-
T
0 where β is the thermal expansion coefficient and ρ
0 is the fluid density at the reference temperature T
0
=
T
h
+
T
c
/
2.Figure 1
Cavity configuration and boundary conditions.The classic thermohydrodynamic equations governing the flow in the cavity with the Boussinesq approximation take the form(1a)
∇
.
u
=
0
(1b)
∂
t
u
+
u
.
∇
u
=
-
1
ρ
0
∇
p
+
ν
∇
2
u
-
β
g
T
-
T
c
e
y
+
1
ρ
0
J
×
B
(1c)
∂
t
T
+
u
.
∇
T
=
α
∇
2
T
u, p, and T are the velocity, pressure, and temperature of the fluid, respectively. ν, α, and β are the kinematic viscosity, the thermal diffusivity, and the thermal expansion of the fluid. e
y is the unit vector in the vertical direction. g is the gravity and J
=
σ
-
∇
φ
+
u
×
B is the current density obtained by using Ohm’s Law, and φ is the electric potential. The effect of the electric field E
=
-
∇
φ is neglected as it vanishes everywhere in the domain.The boundary conditions are as follows. For velocities, no-slip boundary conditions are applied on each side of the cavity. For thermal boundary conditions, two types of boundary conditions are used here. The top and bottom walls are considered isothermal and the side walls adiabatic.In order to express the governing equations in dimensionless form, we introduce the following scaling:(2)
x
X
,
Y
→
x
x
,
y
H
,
t
→
t
α
H
2
,
u
U
,
V
→
u
u
,
v
H
α
,
p
→
p
H
2
ρ
0
α
2By introducing this scaling and previous assumption, we obtain the dimensionless form of (1a)–(1c):(3a)
∇
.
u
=
0
(3b)
∂
t
u
+
u
.
∇
u
=
-
∇
p
+
P
r
∇
2
u
+
R
a
P
r
θ
-
H
a
2
P
r
V
e
y
(3c)
∂
t
θ
+
u
.
∇
θ
=
∇
2
θ
θ
=
T
-
T
c
/
T
h
-
T
c and P
r is the Prandtl number P
r
=
ν
/
α. R
a is the Rayleigh number R
a
=
β
g
T
h
-
T
c
H
3
/
ν
α and H
a
=
B
H
σ
/
μ is the ratio of Lorentz force to the viscous force known as Hartmann number. σ is the electrical conductivity and μ is the dynamic viscosity. Hence the problem the problem is fundamentally governed by the prescribed dimensionless parameters R
a
,
P
r
,
H
a and the boundary conditions can be written in dimensionless quantities as follows:(4)
u
X
=
0
,
Y
=
u
X
=
1
,
Y
=
u
X
,
Y
=
0
=
u
X
,
Y
=
1
=
0
,
θ
X
,
Y
=
0
=
1
,
θ
X
,
Y
=
1
=
0
,
∂
θ
∂
X
X
=
0
=
∂
θ
∂
X
X
=
1
=
0
.


## 3. Double-Population LW-ACM Method

The artificial compressibility method substitutes the equation of conservation of mass and momentum (1a)-(1b) into the alternative form of the momentum equation (5a) and the artificial compressibility equation (ACE) (5b) as follows:(5a)
∂
t
′
ρ
′
+
∇
.
u
′
=
0
(5b)
∂
t
′
u
′
+
R
e
×
u
′
.
∇
u
′
=
-
∇
p
′
+
∇
2
u
′
+
F
′
R
e is the Reynolds number based on some reference velocity u
0 and reference length H. With the dimensionless time t
′
=
t
×
ν
/
H
2, the velocity u
′
=
u
/
u
0, the pressure p
′
=
p
×
H
/
ρ
0
ν
u
0, and the dimensionless external force F
′
=
F
×
ν
u
0
/
H
2. F is equal to the sum of buoyant and magnetic force. u
0 is a reference velocity that can be taken equal to the characteristic velocity in the buoyant flow in our case. Note that by this choice Reynolds number is equal to the square root of Grashof number such that R
e
=
G
r
=
R
a
/
P
r. The artificial density ρ
′ can be determined easily by using the equation of state p
′
=
c
s
2
ρ
′. The artificial speed of sound is c
s
=
1
/
δ, where δ is a disposable parameter, analogous to a relaxation parameter, and known as the artificial compressibility. The auxiliary system with (5a)-(5b) converges much more rapidly to the same steady solution as (1a)-(1b). The introduction of the artificial compressibility δ into the equation of motion does not affect the results and acts in the same manner as the relaxation parameter in LBM.The explicit time-marching and working on a regular Cartesian grid make some analogies between ACM and LBM. This is the starting point of the link-wise formulation of the ACM. The LW-ACM is developed in the same manner as LBM models starting from the ACM and using analog integration strategy. The particle distributions defined for the finite set of the discrete particle velocity vectorsc
i at a site r at time t are denoted f
i
r
,
t, i
=
0
,
.
.
8. In two dimensions, the direction of a single particle probability distribution function is limited to nine directions as follows: (
0,0
) for i
=
0, ±
1,0 and0
,
±
1 for i
=
1
-
4 and ±
1
,
±
1 for i
=
5
-
8. In the LW-ACM, the evolution process of the distribution function f
i can be written as(6a)
f
i
r
+
c
i
,
t
+
1
-
f
i
e
q
r
,
t
=
2
1
-
τ
m
f
i
e
q
,
o
r
+
c
i
,
t
-
f
i
e
q
,
o
r
,
t
+
F
i
r
,
t
(6b)
f
i
e
q
,
o
r
,
t
=
1
2
f
i
e
q
ρ
,
u
-
f
i
e
q
ρ
,
-
uwhere c
i is the discrete speed, τ
m is the relaxation time linked to the kinematic viscosity in a lattice unit ν
l
b
m and given by τ
m
=
0.5
+
3
ν
l
b
m, and ρ and u are the density and velocity in the lattice unit, respectively. F
i is the projection of the body force onto the velocity space and will be detailed later. f
i
e
q is the local equilibrium distribution function that has an appropriately prescribed functional dependence on the local hydrodynamic properties as follows:(7)
f
i
e
q
ρ
,
u
=
ω
i
ρ
1
+
3
u
.
c
i
+
4.5
u
.
c
i
2
-
1.5
u
2The weighting factors ω
i are given by ω
0
=
4
/
9, ω
1
-
4
=
1
/
9 and ω
5
-
8
=
1
/
36. The pressure and the flow velocity in lattice unit are obtained through moment summations in the velocity space as follows:(8)
p
r
,
t
=
1
3
∑
i
f
i
r
,
t
,
u
=
∑
i
c
i
f
i
r
,
t
3
p
r
,
tThe speed of sound remains constant throughout this study and is equal to 1
/
3 in a lattice unit. The asymptotic analysis in [7] showed that the pressure p
p
h and the velocity u
p
h in the physical units can be easily computed by taking p
p
h
=
p
-
p
0
/
ϵ
2 and u
p
h
=
u
/
ϵ, where ϵ
=
1
/
N and N is the number of mesh points along the characteristic length of the problem. Quantities in a lattice units p and u are readily computed in the LW-ACM code.To solve the advection-diffusion equation for the temperature, we use the D2Q9 model. Note that using the same stencil as for dynamic fields (D2Q9) is more convenient in terms of implementation (see [11]). The thermal evolution equation can be written as follows:(9a)
g
i
r
+
c
i
,
t
+
1
-
g
i
e
q
r
,
t
=
2
1
-
τ
t
g
i
e
q
,
o
r
+
c
i
,
t
-
g
i
e
q
,
o
r
,
t
(9b)
g
i
e
q
,
o
r
,
t
=
1
2
g
i
e
q
θ
,
u
-
g
i
e
q
θ
,
-
u
τ
t is the relaxation time in lattice a unit and is related to the thermal diffusivity and g
i
e
q is the equilibrium distribution function defined as follows:(10)
g
i
e
q
θ
,
u
=
ω
i
θ
1
+
3
u
.
c
i
+
4.5
u
.
c
i
2
-
1.5
u
2Terms of second order in u are removed from (10) in order to maintain the linearity of the energy equation with regard to the velocity. The temperature in a lattice unit is computed by conserving only the first moment, i.e., θ
r
,
t
=
∑
i
g
i
r
,
t. The relation between the diffusion coefficient in lattice unit α
l
b
m and the relaxation time is taken as τ
t
=
0.5
+
3
α
l
b
m. With this formulation, the model solves the following energy equation:(11)
∂
t
θ
+
u
.
∇
θ
=
α
l
b
m
∇
2
θThe introduction of the force term in equation (6a) can be done in various ways. It consists of changing the equilibrium distribution functions by applying the contribution of the force in each direction. The body force F
i giving rise to the external force F is defined such that(12)
F
i
=
3
ω
i
c
i
.
F
,
F
=
ρ
0
β
g
θ
-
θ
0
-
H
a
2
μ
N
y
2
u
y
e
ywhere N
y is the number of nodes in the y direction. Note that this assumption is based on the low magnetic Reynolds number (or quasistatic) approximation. Otherwise the magnetic induction model should be taken into consideration to have reliable results [13]. In this case magnetic field should be computed separately using a third distribution function [14].As for LBM, the treatment of boundary condition in LW-ACM is not applied to the macroscopic quantities but directly applied to the distribution function at the boundary limit. Let us supposer
f is the fluid node at which we want to compute the unknown distribution function (close to the solid wall r
b ). i and i
- are the direction of the discrete velocity and its inverse direction, respectively. The no-slip boundary condition is ensured by reversing the velocity of the particle which wants to enter a solid area:(13)
f
i
r
f
,
t
+
1
=
f
i
e
q
r
f
,
t
+
2
1
-
τ
m
f
i
-
e
q
,
o
r
f
,
tNote that this assumption supposes that the wall is halfway between r
f and r
b. For thermal boundary conditions, isothermal boundaries are applied by using the normal equilibrium condition. To determine unknown densities for straight boundaries, the normal equilibrium condition gives(14)
g
i
-
r
f
,
t
+
1
=
-
g
i
r
f
,
t
+
g
i
e
q
r
f
,
t
+
g
i
-
e
q
r
f
,
tAdiabatic boundary conditions are insured by setting the heat transferred to zero:(15)
g
i
-
r
f
,
t
+
1
=
g
i
r
f
,
t
+
ϕ
nwhere the heat flux ϕ
n
=
0 for insulated surfaces.By this manner, LW-ACM satisfies the ACM previously introduced by Chorin [6]. For computer implementation, the kinematic viscosity ν and thermal diffusivity α can be easily computed using P
r and R
a and other parameters in simulations:(16)
ν
=
P
r
β
g
Δ
T
H
3
R
a
1
/
2
,
α
=
ν
P
r
=
β
g
Δ
T
H
3
P
r
R
a
1
/
2As can be seen from (3a)–(3c), the dynamic and thermal fields depend on the dimensionless parameters, namely, the Prandtl, Rayleigh, and Hartmann numbers. In the classic CFD solvers, this type of problem is fully defined by the prescribed nondimensional parameters P
r
,
R
a
,
H
a. However, for LW-ACM, as for LBM implementation, these three parameters are not enough to compute τ
m and τ
t since the thermal expansion β is known as a physical property but not known in lattice units. For this reason, it is better to add another parameter, namely, the Mach number. Hence, the Mach number is fixed within the incompressible limit (
M
a
≪
1
). This choice is motivated by the fact that the Mach number adjustment contributes to the instabilities over LBM simulations on the grid scale [15, 16] and it avoids possible velocity overshoots in the ACM calculations as suggested in [6]. In this way, the characteristic velocity in the buoyant flow u
b
=
β
g
∆
T
H is taken ​​ sufficiently small compared to the speed of sound c
s for all simulations.The compact form of (6a) and (9a) simplifies the adaptation of an existing BGK-LBM code to LW-ACM. Motivated by this fact, the algorithm of the double-population LW-ACM looks as follows:(i)
Setup of dimensionless parameters:R
a, P
r, H
a, and M
a.
(ii)
Initialization ofρ, θ, u, f
i
ρ
,
u, and g
i
θ
,
u using (7) and (10).
(iii)
Advection off
i and g
i.
(iv)
Compute the equilibrium functionsf
i
e
q
,
o and g
i
e
q
,
o at fluid nodes using (6b) and (9b).
(v)
Apply the dynamic and thermal boundary conditions.
(vi)
Computeθ
=
∑
i
g
i and the external force F.
(vii)
Adding the projection of the force termF
i to the distribution function f
i using (12).
(viii)
Compute the conserved quantitiesρ
=
∑
i
f
i, u
=
∑
i
c
i
f
i
/
∑
i
f
i.
(ix)
Ifu and θ reach the convergence criterion, computation terminated.
(x)
Otherwise, processes (iii)-(viii) are repeated.Results are assumed to be converged when we reach the following convergence criterion:(17a)
∑
i
,
j
u
i
,
j
,
t
+
1
-
u
i
,
j
,
t
2
∑
i
,
j
u
i
,
j
,
t
2
<
10
-
8
(17b)
m
a
x
θ
i
,
j
,
t
+
1
-
θ
i
,
j
,
t
<
10
-
6The subscript sequence i
,
j represents the space Cartesian coordinates x and y in the standard basis e
x
,
e
y with the summation is over the entire system and 2 denotes the L
2 norm. The volume average Nusselt number which gives an idea about the overall heat transfer in the whole cavity is computed as follows:(18)
N
u
=
1
+
u
y
θ
H
α
Δ
θ
=
1
+
1
N
x
α
l
b
m
Δ
θ
∑
i
,
j
u
y
i
,
j
θ
i
,
j
u
y is the y-component of the velocity and the summation is over the entire system.


## 4. Validation

### 4.1. Validation of the MRT-LBM Code

As it has been developed in the introduction, LW-ACM is a new formulation of ACM that operates on lattices similar to LBM. Some important points have to be discussed and analyzed before performing simulations using this model, such as its efficiency in terms of speed and accuracy against BGK which is the most widely used model and MRT which is the most stable and accurate one [17].The MRT-LBM code is validated by two-dimensional simulation of Rayleigh-Bénard convection at different Rayleigh numbers with a Prandtl number of 0.71 (air). The critical Rayleigh number for the onset of the Rayleigh-Bénard convection isR
a
c
=
1707.74 which agrees with the theoretical prediction [18]. As the Rayleigh number increased, the steady two-dimensional convection rolls become unstable. The wavy instability and periodic motion observed are in good agreement with the well-known experimental observations and theoretical predictions. Figure 2 shows a snapshot of the temperature distribution for R
a
=
10
4 and P
r
=
0.71 using MRT-LBM model. At the first level of the simulation, the generation of the mushroom-like isotherms (plumes) is observed. It is caused by the instability and the fluctuation of the instantaneous thickness of temperature boundary layer. This phenomenon can be observed in turbulent convection at high Rayleigh number. The isotherms show a higher level of convective activity and thin thermal boundary layers. It is important to mention here that similar results of snapshot contour plots are found experimentally by Sparrow et al. [19] using electrically heated horizontal copper surfaces situated in a water medium. They used an electrochemical technique to facilitate physical observations. As observed from Figure 2, one can affirm that the MRT-LBM code reproduces the mushroom-like appearance of the isotherms as observed by the reference.Figure 2
Snapshot contour plot of temperature forR
a
=
10
4 and P
r
=
0.71 using MRT-LBM model.Another validation of the MRT-LBM code was obtained by the mixed convection heat transfer of fluid over a vertical duct preceded with a double-step expansion. The fluid is injected into the channel at a cooler temperatureT
c and with a fully developed velocity profile. Channel expansion ratio A
=
L
/
h is taken equal to 30 and the overall length of the computational domain is L
+
35
h, where h is the step height. The walls downstream of the step are maintained at a constant temperature T
h, while the other walls are treated adiabatically. At the inlet, a fully developed parabolic profile v was enforced. With f
2 known from the nonequilibrium part of the population normal to the boundary, f
5 and f
6 can be found as f
5
r
i
n
=
f
7
r
i
n
+
1
/
2
f
3
r
i
n
-
f
1
r
i
n
+
1
/
6
ρ
v and f
6
r
i
n
=
f
8
r
i
n
+
1
/
2
f
1
r
i
n
-
f
3
r
i
n
+
1
/
6
ρ
v. At all walls, bounce-back boundary conditions were applied. The artificial open boundary condition is applied in the outlet of the channel by copying the fields from the penultimate column f
i
r
o
u
t
=
f
i
r
o
u
t
-
1. Both the ‘bounce-back’ and the ‘antibounce-back’ rules [2] were used to implement isothermal and adiabatic boundary conditions, respectively.The results for mixed convection flow with airP
r
=
0.71 have been compared with numerical results reported in using finite-volume method. Computations were performed on 702
×
62 grid size. Velocity profiles at the particular axial position Y=5 for fixed Reynolds number R
e
=
114 and for the Grashof number ranging from G
r
=
10
2 to G
r
=
10
4 have been compared with those obtained by Tsui and Shu [20] in Figure 3. Here the Reynolds number is defined as R
e
=
v
0
H
/
ν, where v
0 is the centerline fluid velocity at the inlet. The Grashof number is also based on the step height G
r
=
β
g
T
h
-
T
c
H
3
/
ν
2. It is shown that the numerical results of the MRT-LBM model compare well with the numerical data. Hence, results of MRT-LBM code can be used as a benchmark when testing the LW-ACM. More importantly, it is well known that the MRT-LBM scheme is more accurate, stable, and efficient than the BGK-LBM scheme [17].Figure 3
Mixed convection in vertical sudden expansion (Left) and velocity profiles at y/h=5 (Right) with the MRT-LBM model forR
e
=
114 and P
r
=
0.71.Both streamlines and steady state isotherms are shown in Figure4. As depicted, the asymmetric flow is detected at Grashof number G
r
=
10
2. For G
r
=
2
×
10
3, the flow becomes symmetric and the central region is nearly stagnant. However, for Grashof number G
r
=
10
4, flow reversal takes place near the centerline region. Similar results are reported by Tsui and Shu [20] as well.Figure 4
Streamlines (left) and isotherms (right) for (Pr=0.71) and Grashof numbers (a)G
r
=
10
2, (b) G
r
=
2
×
10
3, and (c) G
r
=
10
4.


### 4.2. Comparison of the LW-AC Method with Other LB Methods

Table1 shows the computed Nusselt numbers and CPU times for Pr=0.71 obtained by the BGK-LBM, MRT-LBM, and the double LW-ACM method for a square cavity with differentially heated side walls problem at Rayleigh number from 103 to 105. The flow configuration and boundary conditions are as follows. The right and left boundaries are maintained at constant temperatures θ
=
1and θ
=
0, respectively, while the other two walls are kept adiabatic. The no-slip boundary condition is applied at walls. The computing of averaged Nusselt number representing the heat transfer rate at the isothermal walls is defined as follows:(19)
N
u
=
1
H
Δ
θ
∫
0
H
∂
θ
∂
x
d
y
≃
1
N
y
Δ
θ
∑
1
N
y
3
θ
0
,
j
-
4
θ
1
,
j
+
θ
2
,
j
2 where H is the cavity height and N
y is the number of lattice nodes in the vertical direction. The results of Wang et al. [2] using MRT-LBM model are used as reference. As it can be seen from all chosen Rayleigh numbers, the LW-ACM model gives results closer to the reference than the BGK-LBM model. On the other hand, the MRT-LBM model gives analogous results as those reported in the reference even we use a coarse mesh. The relative deviation between LW-ACM and the reported results obtained by the coarse 642 mesh is about 0.044
% for R
a
=
10
3, 0.022
% for R
a
=
10
4, and 0.079
% for R
a
=
10
5. As expected, at low Rayleigh number, values of the Nusselt number are close to the reference data while the deviation is more pronounced at higher Rayleigh number. On the other hand, the maximum deviation of the finest grid solution is about 0.035
% at R
a
=
10
5. The double LW-ACM model gives useful results within acceptable errors and it may be used as an effective tool for thermal problems. The table also provides the CPU time required to reach the steady state when using BGK-LBM, MRT-LBM, and LW-ACM models. Clearly, we can see that the CPU time for BGK-LBM model is the shortest. This gives an explanation of why BGK-LBM is the largely used and the most referred model. Except for the first computation using the coarse mesh at R
a
=
10
3, the CPU time required for the MRT model is always high. The comparison between several models with respect to the CPU time shows that the slowness of the MRT model compared to the BGK model l
=
C
P
U
M
R
T
-
C
P
U
B
G
K
/
C
P
U
B
G
K is about 23.2
%, 19.7
%, and 32
% using the intermediate 1282 grid size and at R
a
=
10
3, R
a
=
10
4, and R
a
=
10
5, respectively. It is true that MRT-LBM model improves the accuracy of simulations, but it requires more time to reach steady state. The double LW-ACM in turn requires about 11
%, 12.7
%, and 14
% using the same grid size and Rayleigh number.Table 1
Computed Nu and CPU time with different models,P
r
=
0.71 and R
a
=
10
3
-
10
5.
Ra
Grid size
Wang et al. (2013)
BGK-LBM
Δ(%)
CPU
MRT-LBM
Δ(%)
CPU
LW-ACM
Δ(%)
CPU
103
642
1.1178
1.1170
0.071
30.11
1.1174
0.035
40.11
1.1173
0.044
44.65
1282
1.1171
0.062
63.31
1.1175
0.026
78.01
1.1174
0.035
70.23
2562
1.1172
0.053
81.22
1.1178
0.0
110.22
1.1178
0.0
95.368
104
642
2.2448
2.2430
0.080
45.21
2.2445
0.013
54.26
2.2443
0.022
48.05
1282
2.2436
0.053
76.02
2.2443
0.022
91.01
2.2440
0.035
85.67
2562
2.2444
0.017
140.21
2.2446
0.008
180.41
2.2444
0.017
156.51
105
642
4.5216
4.5172
0.097
65.25
4.5185
0.068
80.21
4.5180
0.079
72.38
1282
4.5185
0.068
150.21
4.5195
0.046
198.25
4.5193
0.050
171.32
2562
4.5110
0.013
320.00
4.5202
0.030
400.25
4.5200
0.035
361.49The second step in the validation process of the LW-ACM code is the comparison of the given results with those reported in the literature using stability analysis. The average Nusselt number for Rayleigh number ranging fromR
a
=
10
3 to R
a
=
10
5 has been compared with those obtained by Clever and Busse [21] in Figure 5. It is shown that our numerical results compare well with both numerical and analytical data obtained from the correlation N
u
=
1.56
R
a
/
R
a
c
0.296. At low Rayleigh number, the grid resolution is sufficient to solve the flow and heat transfer. However, for high Rayleigh number, the LW-ACM model slightly underestimates the heat transfer rate.Figure 5
Variation of Nusselt number for 2D Rayleigh-Bénard problem with Rayleigh number for Pr=0.71 (air). Results are obtained using LW-ACM.


## 4.1. Validation of the MRT-LBM Code

As it has been developed in the introduction, LW-ACM is a new formulation of ACM that operates on lattices similar to LBM. Some important points have to be discussed and analyzed before performing simulations using this model, such as its efficiency in terms of speed and accuracy against BGK which is the most widely used model and MRT which is the most stable and accurate one [17].The MRT-LBM code is validated by two-dimensional simulation of Rayleigh-Bénard convection at different Rayleigh numbers with a Prandtl number of 0.71 (air). The critical Rayleigh number for the onset of the Rayleigh-Bénard convection isR
a
c
=
1707.74 which agrees with the theoretical prediction [18]. As the Rayleigh number increased, the steady two-dimensional convection rolls become unstable. The wavy instability and periodic motion observed are in good agreement with the well-known experimental observations and theoretical predictions. Figure 2 shows a snapshot of the temperature distribution for R
a
=
10
4 and P
r
=
0.71 using MRT-LBM model. At the first level of the simulation, the generation of the mushroom-like isotherms (plumes) is observed. It is caused by the instability and the fluctuation of the instantaneous thickness of temperature boundary layer. This phenomenon can be observed in turbulent convection at high Rayleigh number. The isotherms show a higher level of convective activity and thin thermal boundary layers. It is important to mention here that similar results of snapshot contour plots are found experimentally by Sparrow et al. [19] using electrically heated horizontal copper surfaces situated in a water medium. They used an electrochemical technique to facilitate physical observations. As observed from Figure 2, one can affirm that the MRT-LBM code reproduces the mushroom-like appearance of the isotherms as observed by the reference.Figure 2
Snapshot contour plot of temperature forR
a
=
10
4 and P
r
=
0.71 using MRT-LBM model.Another validation of the MRT-LBM code was obtained by the mixed convection heat transfer of fluid over a vertical duct preceded with a double-step expansion. The fluid is injected into the channel at a cooler temperatureT
c and with a fully developed velocity profile. Channel expansion ratio A
=
L
/
h is taken equal to 30 and the overall length of the computational domain is L
+
35
h, where h is the step height. The walls downstream of the step are maintained at a constant temperature T
h, while the other walls are treated adiabatically. At the inlet, a fully developed parabolic profile v was enforced. With f
2 known from the nonequilibrium part of the population normal to the boundary, f
5 and f
6 can be found as f
5
r
i
n
=
f
7
r
i
n
+
1
/
2
f
3
r
i
n
-
f
1
r
i
n
+
1
/
6
ρ
v and f
6
r
i
n
=
f
8
r
i
n
+
1
/
2
f
1
r
i
n
-
f
3
r
i
n
+
1
/
6
ρ
v. At all walls, bounce-back boundary conditions were applied. The artificial open boundary condition is applied in the outlet of the channel by copying the fields from the penultimate column f
i
r
o
u
t
=
f
i
r
o
u
t
-
1. Both the ‘bounce-back’ and the ‘antibounce-back’ rules [2] were used to implement isothermal and adiabatic boundary conditions, respectively.The results for mixed convection flow with airP
r
=
0.71 have been compared with numerical results reported in using finite-volume method. Computations were performed on 702
×
62 grid size. Velocity profiles at the particular axial position Y=5 for fixed Reynolds number R
e
=
114 and for the Grashof number ranging from G
r
=
10
2 to G
r
=
10
4 have been compared with those obtained by Tsui and Shu [20] in Figure 3. Here the Reynolds number is defined as R
e
=
v
0
H
/
ν, where v
0 is the centerline fluid velocity at the inlet. The Grashof number is also based on the step height G
r
=
β
g
T
h
-
T
c
H
3
/
ν
2. It is shown that the numerical results of the MRT-LBM model compare well with the numerical data. Hence, results of MRT-LBM code can be used as a benchmark when testing the LW-ACM. More importantly, it is well known that the MRT-LBM scheme is more accurate, stable, and efficient than the BGK-LBM scheme [17].Figure 3
Mixed convection in vertical sudden expansion (Left) and velocity profiles at y/h=5 (Right) with the MRT-LBM model forR
e
=
114 and P
r
=
0.71.Both streamlines and steady state isotherms are shown in Figure4. As depicted, the asymmetric flow is detected at Grashof number G
r
=
10
2. For G
r
=
2
×
10
3, the flow becomes symmetric and the central region is nearly stagnant. However, for Grashof number G
r
=
10
4, flow reversal takes place near the centerline region. Similar results are reported by Tsui and Shu [20] as well.Figure 4
Streamlines (left) and isotherms (right) for (Pr=0.71) and Grashof numbers (a)G
r
=
10
2, (b) G
r
=
2
×
10
3, and (c) G
r
=
10
4.


## 4.2. Comparison of the LW-AC Method with Other LB Methods

Table1 shows the computed Nusselt numbers and CPU times for Pr=0.71 obtained by the BGK-LBM, MRT-LBM, and the double LW-ACM method for a square cavity with differentially heated side walls problem at Rayleigh number from 103 to 105. The flow configuration and boundary conditions are as follows. The right and left boundaries are maintained at constant temperatures θ
=
1and θ
=
0, respectively, while the other two walls are kept adiabatic. The no-slip boundary condition is applied at walls. The computing of averaged Nusselt number representing the heat transfer rate at the isothermal walls is defined as follows:(19)
N
u
=
1
H
Δ
θ
∫
0
H
∂
θ
∂
x
d
y
≃
1
N
y
Δ
θ
∑
1
N
y
3
θ
0
,
j
-
4
θ
1
,
j
+
θ
2
,
j
2 where H is the cavity height and N
y is the number of lattice nodes in the vertical direction. The results of Wang et al. [2] using MRT-LBM model are used as reference. As it can be seen from all chosen Rayleigh numbers, the LW-ACM model gives results closer to the reference than the BGK-LBM model. On the other hand, the MRT-LBM model gives analogous results as those reported in the reference even we use a coarse mesh. The relative deviation between LW-ACM and the reported results obtained by the coarse 642 mesh is about 0.044
% for R
a
=
10
3, 0.022
% for R
a
=
10
4, and 0.079
% for R
a
=
10
5. As expected, at low Rayleigh number, values of the Nusselt number are close to the reference data while the deviation is more pronounced at higher Rayleigh number. On the other hand, the maximum deviation of the finest grid solution is about 0.035
% at R
a
=
10
5. The double LW-ACM model gives useful results within acceptable errors and it may be used as an effective tool for thermal problems. The table also provides the CPU time required to reach the steady state when using BGK-LBM, MRT-LBM, and LW-ACM models. Clearly, we can see that the CPU time for BGK-LBM model is the shortest. This gives an explanation of why BGK-LBM is the largely used and the most referred model. Except for the first computation using the coarse mesh at R
a
=
10
3, the CPU time required for the MRT model is always high. The comparison between several models with respect to the CPU time shows that the slowness of the MRT model compared to the BGK model l
=
C
P
U
M
R
T
-
C
P
U
B
G
K
/
C
P
U
B
G
K is about 23.2
%, 19.7
%, and 32
% using the intermediate 1282 grid size and at R
a
=
10
3, R
a
=
10
4, and R
a
=
10
5, respectively. It is true that MRT-LBM model improves the accuracy of simulations, but it requires more time to reach steady state. The double LW-ACM in turn requires about 11
%, 12.7
%, and 14
% using the same grid size and Rayleigh number.Table 1
Computed Nu and CPU time with different models,P
r
=
0.71 and R
a
=
10
3
-
10
5.
Ra
Grid size
Wang et al. (2013)
BGK-LBM
Δ(%)
CPU
MRT-LBM
Δ(%)
CPU
LW-ACM
Δ(%)
CPU
103
642
1.1178
1.1170
0.071
30.11
1.1174
0.035
40.11
1.1173
0.044
44.65
1282
1.1171
0.062
63.31
1.1175
0.026
78.01
1.1174
0.035
70.23
2562
1.1172
0.053
81.22
1.1178
0.0
110.22
1.1178
0.0
95.368
104
642
2.2448
2.2430
0.080
45.21
2.2445
0.013
54.26
2.2443
0.022
48.05
1282
2.2436
0.053
76.02
2.2443
0.022
91.01
2.2440
0.035
85.67
2562
2.2444
0.017
140.21
2.2446
0.008
180.41
2.2444
0.017
156.51
105
642
4.5216
4.5172
0.097
65.25
4.5185
0.068
80.21
4.5180
0.079
72.38
1282
4.5185
0.068
150.21
4.5195
0.046
198.25
4.5193
0.050
171.32
2562
4.5110
0.013
320.00
4.5202
0.030
400.25
4.5200
0.035
361.49The second step in the validation process of the LW-ACM code is the comparison of the given results with those reported in the literature using stability analysis. The average Nusselt number for Rayleigh number ranging fromR
a
=
10
3 to R
a
=
10
5 has been compared with those obtained by Clever and Busse [21] in Figure 5. It is shown that our numerical results compare well with both numerical and analytical data obtained from the correlation N
u
=
1.56
R
a
/
R
a
c
0.296. At low Rayleigh number, the grid resolution is sufficient to solve the flow and heat transfer. However, for high Rayleigh number, the LW-ACM model slightly underestimates the heat transfer rate.Figure 5
Variation of Nusselt number for 2D Rayleigh-Bénard problem with Rayleigh number for Pr=0.71 (air). Results are obtained using LW-ACM.


## 5. Magnetic Field Effect

In this section, numerical simulations have been performed for square cavity heated from below and cooled from above in the presence of magnetic fields. This can evaluate the double-population LW-ACM in solving more complex flows. Recently, this problem sustains a considerable attention because of a number of its wide variety of applications in engineering and technology such as metallic alloys, magnetoelectrochemistry (MEC), and astrophysical and environmental systems. Thus, it has been shown both experimentally and numerically that magnetic fields can be used to control thermal convection which is important in problems involving metallurgy, microstructure devices, and protein crystals and under reduced gravity conditions. In order to evaluate the potential effects of magnetic fields in heat transfer area, Benos et al. [22] solved analytically and numerically two-dimensional MHD natural convection flow in an internally heated horizontal shallow cavity. Their comparison of the analytical and numerical results showed the validity and the correctness of their analysis. It is of interest to note that the recent progress in superconductivity and electric resistance of metals at very low temperature allowed researchers to reach larger magnetic fields up to 15 T by using super-conducting magnets. Effects of strong magnetic field on two-dimensional natural convection have been studied by Pirmohammadi and Ghassemi [23]. They found that when the Hartmann number is sufficiently large, the convection is suppressed. In fact, in spite of the existence of many papers dealing with the interaction between convection and magnetic fields, there is as yet no complete understanding of magnetoconvection and further research works still needed to provide certain important features. Motivated by this fact, we apply the LW-ACM to the two-dimensional magnetoconvection problem as described in Section 2.Figure6 shows the streamlines, the steady state isotherms and the pressure contours at Rayleigh number R
a
=
5
×
10
4 without magnetic fields (
H
a
=
0
) and with magnetic fields for Hartmann number 60. We observe that the thermal gradient gives rise to two convection rolls characterizing the instability of Rayleigh-Bénard. The number that governs this type of flow is the Rayleigh number. For all chosen Rayleigh numbers, an adjustment in the flow pattern is observed when magnetic fields are applied. The rotational flow with a single elliptic vortex is observed for high Hartmann number. Isotherms become more uniformly spaced throughout the cavity and the intensity of convection is considerably decreased by the drag induced by the magnetic field, as indicated by a weak distortion of the isothermal lines. The magnetic field reduces the Nusselt number and the fluid velocity. When the Hartmann number is sufficiently large, the mushroom-like isotherms (plumes) caused by the instability of the thermal boundary layer become thinner which tends to homogenize the temperature throughout the cavity. Without magnetic fields, the streamlines form two symmetric vortices with respect to the vertical centerline of the enclosure and when a magnetic field is applied with a sufficiently high value of Hartmann number, the two vortices move from the center of the cavity to the top. The bifurcation of the flow takes place for all cases. The two large vortices appear on the center of the enclosure, but with increasing in Hartmann number, the two vortices are reduced in size located near the top and bottom of the right horizontal top wall. The core of vortex depends on the magnitude of magnetic field. Also, the increase in the value of the magnetic field moves the vortex into the top wall where the temperature is cold. We note also a change in the shape of vortices from the circular shape to elliptical one. One of the application of such result is in melting processes of an electrically conducting fluid metals when the fluid motion presents turbulent convection flow. The applied magnetic field stabilizes or cancels unwanted oscillations in the flow.Figure 6
(a) Streamlines, (b) steady state isotherms, and (c) pressure contours atP
r
=
0.054 and R
a
=
5
×
10
4 for Hartmann number H
a from 0 to 60 (from top to bottom) using LW-ACM model.Figure7 presents the evolution of the averaged Nusselt number for different Hartmann numbers. The convergence history indicates that the volume average Nusselt number has a temporal damped oscillation at the beginning of the simulation. This results from incipient instabilities since all the Nusselt numbers have converged to a constant values after a number of iterations. The Nusselt number keeps decreasing with respect of Hartmann number; we can highlight that when H
a
=
60, the effect of magnetic fields on the Nusselt number (and so on the heat transferred through the walls) is significant.Figure 7
Convergence of the volume-averaged Nusselt number with different Hartmann number atR
a
=
5
×
10
4 using LW-ACM (for interpretation of the colors in the figure(s), the reader is referred to the web version of this article.).For the sake of completeness, computations of the volume-averaged Nusselt number for different magnetic field magnitude give4.261, 4.185, 2.850, and 2.768 for H
a
=
0, H
a
=
20, H
a
=
40, and H
a
=
60, respectively. In agreement with what has been concluded from the thermal and flow fields in Figure 6, increasing Hartmann number is associated with a decrease in Nusselt number and the heat transfer in the cavity. In fact, increasing magnetic strength leads to an increase the Lorenz force which prevents buoyant force and substantially suppresses the convection. More importantly, the decrease in the heat transfer reaches 1.78
% when H
a increases by 20 compared to the case without magnetic field. When H
a increases from 20 to 40, the decrease in heat transfer is about 31.89
%.Results of natural convective flow in the presence of magnetic fields make it possible to judge the ability of the double-population LW-ACM model to predict realistic flows. Thus, the method seems to be able to predict quantities of engineering interest like the volume-averaged Nusselt number accurately. However, as MRT-LBM and LW-ACM methods are comparable, the question is why should double-population LW-ACM be preferred to BGK-LBM with local grid refinement and MRT-LBM?First, the LW-ACM method is a formulation of ACM that operates on a lattice. Therefore, it does not involve an enforcement of mass conservation at each time step. This causes abrupt changes in particle distributions and breaks the mass conservation especially at the boundary for BGK-LBM. This generates errors near the boundary and destroys uniform convergence to the solution of LBM schemes.Second, the double-population LW-ACM has the advantage of lower memory usage and computational effort requirement. Thus, if we look carefully at (6a)-(6b) and (9a)-(9b), only the distribution functions f
i and g
i have to be stored whereas, for a MRT-LBM code, it is necessary to add much many quantities at each position r since the transformation matrix has eigenvectors connected with the velocity set.Finally, the implementation of MRT-LBM code is somehow sensitive to the choice of the relaxation parameters which need to be adjusted by a linear stability analysis of the LBE evolution, while only physical quantities governing the problem are needed in order to implement a double LW-ACM code.


## 6. Conclusions

In this work, a double-population Link-Wise Artificial Compressibility Method (LW-ACM) has been developed to solve convective flows in the presence of a magnetic field. Computation of mixed convective flow over a vertical duct preceded by a sudden expansion was performed in order to validate the MRT model. Then, results of natural convection flows in a differentially heated cavity problem using LW-ACM were compared and validated with those obtained using BGK-LBM and MRT-LBM models and the accuracy of the prescribed model was evaluated. The numerical results for values of Nusselt numbers and process time prove that LW-ACM is an efficient tool for fluid flow and heat transfer computations. Therefore, from a practical point of view, this eliminates the need for using the MRT model which is more expensive and more demanding in terms of memory. Then we applied the numerical model to study magnetoconvection in electrically conducting fluid withP
r
=
0.054, for R
a
=
5
×
10
4 and for Hartmann number up to 60. Results showed that the presence of magnetic fields affects directly the buoyant force. The transversal velocity changed considerably which has a direct manifestation on the volume average Nusselt number. This trend demonstrates that LW-ACM is competitive tool and a good alternative to LBM models for engineering applications.



---

*Source: 7834731-2019-03-17.xml*