- Email: [email protected]
Multiple-relaxation-time lattice Boltzmann modeling of incompressible flows in porous media
Physica A xx (xxxx) xxx–xxx
Contents lists available at ScienceDirect
Physica A journal homepage: www.elsevier.com/locate/physa
Multiple-relaxation-time lattice Boltzmann modeling of incompressible flows in porous media Q1
Qing Liu, Ya-Ling He ∗ Key Laboratory of Thermo-Fluid Science and Engineering of MOE, School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi, 710049, PR China
highlights • • • •
A D2Q8 MRT-LB model is proposed for simulating incompressible porous flows at the REV scale. The generalized non-Darcy model is employed to describe the momentum transfer in porous media. The generalized Navier–Stokes equations can be recovered through the Chapman–Enskog analysis in the moment space. The MRT-LB model is demonstrated by numerical simulations of several typical two-dimensional porous flows.
article
info
Article history: Received 21 October 2014 Received in revised form 12 January 2015 Available online xxxx Keywords: Lattice Boltzmann model Multiple-relaxation-time Porous media Incompressible flows REV scale
abstract In this paper, a two-dimensional eight-velocity multiple-relaxation-time (MRT) lattice Boltzmann (LB) model is proposed for incompressible porous flows at the representative elementary volume scale based on the Brinkman–Forchheimer-extended Darcy model. In the model, the porosity is included into the pressure-based equilibrium moments, and the linear and nonlinear drag forces of the porous matrix are incorporated into the model by adding a forcing term to the MRT-LB equation in the moment space. Through the Chapman–Enskog analysis, the incompressible generalized Navier–Stokes equations can be recovered. Numerical simulations of several typical porous flows are carried out to validate the present MRT-LB model. It is found that the present numerical results agree well with the analytical solutions and/or other numerical results reported in the literature. © 2015 Published by Elsevier B.V.
1. Introduction Fluid flow and related transport phenomena in porous media have gained significant research interest due to the importance of related technological and industrial applications, which include contaminant transport in groundwater, crude oil exploration and extraction, radioactive waste management, hydrogeology and so on [1–3]. For incompressible flows in porous media at the representative elementary volume (REV) scale, the Darcy model, the Brinkman-extended Darcy model and the Forchheimer-extended Darcy model have been widely employed. However, the Darcy model and the two extended (Brinkman and Forchheimer) models have some intrinsic limitations in simulating porous flows [4,5]. In order to overcome the deficiencies of the above mentioned models, the Brinkman–Forchheimer-extended Darcy model (also called the generalized model) has been developed by several research groups [4–7]. In the generalized model, the viscous and inertial forces are incorporated into the momentum equation by using the local volume-averaging technique at the REV scale. The Darcy model and the two extended models can be regarded as the limiting cases of the generalized model. Due to the
∗
Corresponding author. E-mail address: [email protected] (Y.-L. He).
http://dx.doi.org/10.1016/j.physa.2015.01.067 0378-4371/© 2015 Published by Elsevier B.V.
1
2 3 4 5 6 7 8 9 10
2
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
similarity with the Navier–Stokes equations, the generalized model can be used to simulate transient flows through porous media. Moreover, as reported by Vafai and Kim [8], numerical results based on the Brinkman–Forchheimer-extended Darcy formulation have been shown to be in good agreement with the experimental predictions. In the past several decades, various traditional numerical methods, such as the finite volume (FV) method, the finite difference (FD) method, and the finite element (FE) method, have been employed to study porous flows based on the generalized non-Darcy model. The lattice Boltzmann (LB) method, as a bran-new mesoscopic numerical technique originates from the lattice-gas automata (LGA) method [9], has achieved significant success in modeling complex fluid flows and simulating complex physics in fluids [10–16]. Owing to its kinetic background, the LB method has some distinctive advantages over the traditional numerical methods (e.g., see Ref. [17]). Recently, the LB method has been successfully applied to simulate fluid flows in porous media at the REV scale [18–25]. In the REV scale method, the detailed geometric structure of the media is ignored, and the standard LB equation is modified by adding an additional term to account for the presence of the porous media. Therefore, LB method at the REV scale can be used for systems with large computational domain. It is worth mentioning that the REV LB method is very effective for simulating fluid flows in the region which is partially filled with a porous medium. As reported in Ref. [21], the discontinuity of the velocity-gradient at the porous medium/free-fluid interface can be well captured by the LB method without including the stress boundary condition into the simulation model. To the best of our knowledge, most of the existing REV LB models [18–25] for incompressible porous flows employ the Q2 Bhatnagar–Gross–Krook (BGK) model (also called the single-relaxation-time model) [26] to represent the collision process. Although the BGK-LB model has become the most popular one in the LB community because of its extreme simplicity, there are several well-known criticisms on this model, such as the numerical instability at low values of viscosity [27–29] and the inaccuracy in treating boundary conditions [30]. Fortunately, these shortcomings of the BGK-LB model can be easily addressed by employing the multiple-relaxation-time (MRT) model proposed by d’Humières [31]. Hence, the aim of this paper is to develop a new MRT-LB model for incompressible porous flows at the REV scale based on the generalized model. In the model, a pressure-based MRT-LB equation with the eight-by-eight collision matrix [32] is constructed in the framework of the standard MRT-LB method. The remainder of this paper is organized as follows. In Section 2, the MRT-LB model for incompressible porous flows at the REV scale is presented in detail. In Section 3, numerical tests of the MRT-LB model are performed for the porous Poiseuille flow, porous Couette flow, lid-driven flow in a square porous cavity, and natural convection flow in a square porous cavity. Finally, a brief conclusion is made in Section 4.
28
2. MRT-LB model for incompressible flows in porous media
29
2.1. Macroscopic governing equations
30 31 32 33
The fluid flow is assumed to be two-dimensional, laminar and incompressible. For incompressible flows through porous media at the REV scale, the generalized model proposed by Nithiarasu et al. [4] is employed in the present study. The dimensional governing equations (generalized Navier–Stokes equations) of the generalized model can be written as follows:
∇ · u = 0,
(1)
36
u 1 ∂u + (u · ∇) = − ∇ (φ p) + υe ∇ 2 u + F, (2) ∂t φ ρ0 fluid velocity and pressure, respectively, φ is the where ρ0 is the average fluid density, u and p are the volume-averaged porosity, and υe is the effective kinetic viscosity. F = Fx , Fy denotes the total body force induced by the porous matrix and
37
other external forces, which can be expressed as [6,22]
34
35
38
39 40 41 42
43
44 45 46
47
48
F=−
φυ
φ Fφ u − √ |u|u + φ a, K K
(3)
where υ is the kinetic viscosity of the fluid, K is the permeability, Fφ is the geometric function, a is the body force due to an external force, and |u| =
u2x + u2y , in which ux and uy are the x- and y-components of the macroscopic velocity u,
respectively. Based on Ergun’s relation [33], the geometric function Fφ and the permeability K of the porous media can be expressed as [34] 1.75 Fφ = , 150φ 3
K =
φ 3 d2p 150 (1 − φ)2
,
(4)
where dp is the diameter of the solid particle. The flow governed by the generalized Navier–Stokes equations (1) and (2) are characterized by the porosity φ and several dimensionless parameters: the Darcy number Da, the viscosity ratio J, and the Reynolds number Re, which are defined as K υe LU , J = , Re = , L2 υ υ where L and U are the characteristic length and velocity of the system, respectively. Da =
(5)
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
3
Fig. 1. Discrete velocities of the D2Q8 model.
2.2. MRT-LB model
1
In this subsection, a two-dimensional eight-velocity (D2Q8) MRT-LB model is presented to study incompressible porous flows. According to Refs. [35,36], the MRT-LB equation with an explicit treatment of the forcing term can be written as f (x + eδt , t + δt ) − f (x, t ) = −M−1 3 m (x, t ) − m(eq) (x, t ) + δt M−1
I−
3 2
S,
(6)
where M is the 8 × 8 orthogonal transformation matrix, 3 is the non-negative 8 × 8 diagonal relaxation matrix, and I is the 8 × 8 identity matrix. The boldface symbols, f, m, m(eq) , and S represent 8-dimensional (column) vectors: f (x, t ) = (f1 (x, t ) , f2 (x, t ) , . . . , f8 (x, t )) , T
2 3
4
5 6 7
f (x + eδt , t + δt ) = (f1 (x + e1 δt , t + δt ) , . . . , f8 (x + e8 δt , t + δt )) ,
8
m (x, t ) = (m1 (x, t ) , m2 (x, t ) , . . . , m8 (x, t ))T ,
9
T
m
(eq)
(eq)
(x, t ) = m1
(eq)
(x, t ) , m2
(eq)
(x, t ) , . . . , m8
T (x, t ) ,
10
S = (S1 , S2 , . . . , S8 )T ,
11
where T is the transpose operator, fi (x, t ) is the discrete distribution function with velocity ei at position x and time t (at the REV scale), m (x, t ) and m(eq) (x, t ) are the velocity moments and the corresponding equilibrium moments, respectively, and {Si |i = 1, 2, . . . , 8} are components of the forcing term S. In the D2Q8 model, the eight discrete velocities {ei |i = 1, 2, . . . , 8} (see Fig. 1) are given by [32,37]
ei =
(cos [(i − 1) π /2] , sin [(i − 1) π /2]) c , √ (cos [(2i − 9) π /4] , sin [(2i − 9) π /4]) 2c ,
i=1∼4
(7)
i = 5 ∼ 8,
where c = δx /δt is the lattice speed,√ in which δt and δx are the discrete time step and lattice spacing, respectively. The sound speed of the D2Q8 model is cs = c / 3. In the MRT model, the lattice speed c is usually set to be 1, which leads to δx = δt . The transformation matrix M linearly transforms the discrete distribution functions f ∈ V = R8 (velocity space) to their velocity moments m ∈ M = R8 (moment space): m = Mf,
f = M−1 m.
(8)
For the D2Q8 model, the eight velocity moments {mi |i = 1, 2, . . . , 8} are: m = (m1 , m2 , m3 , m4 , m5 , m6 , m7 , m8 )
= p, e, jx , qx , jy , qy , pxx , pxy
T
1 −1 0 0 1 −2 −1 0
1 −1 −1 2 0 0 1 0
1 −1 0 0 −1 2 −1 0
1 1 1 1 1 1 0 1
,
(9)
1 1 −1 −1 1 1 0 −1
1 1 −1 −1 −1 −1 0 1
14 15
16
17 18 19 20
21
23
where m1 = p is the pressure, m2 = e is related to energy, m3,5 = jx,y are components of the momentum J = jx , jy , m4,6 = qx,y are related to energy flux, and m7,8 = pxx,xy are related to the diagonal and off-diagonal components of the stress tensor, respectively [27]. With the ordering of the velocity moments {mi }, the transformation matrix M can be constructed via the Gram–Schmidt orthogonalization procedure (c = 1) [32]: 1 −1 1 −2 M= 0 0 1 0
13
22
T
12
1 1 1 1 . −1 −1 0 −1
24
25 26 27 28
(10)
29
4
1
2
3
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
The equilibrium moments m(eq) corresponding to the velocity moments m are defined by m(eq) =
5 3
φp +
e(eq) = α2 φ p +
5
q(xeq) =
6
q(yeq) =
7
p(xxeq) =
8
p(xyeq) =
9
11 12 13
16
17
18 19 20 21 22
1 2 1 2 3 2 3 2
3
γ2
4
ρ0 u2x + u2y φ
24 25
26
28
29
31
,
γ1 γ3
ρ0 u2x − u2y φ
,
ρ0 ux uy . φ
(12)
3 = diag (s1 , s2 , s3 , s4 , s5 , s6 , s7 , s8 )
= diag 1, se , 1, sq , 1, sq , sυ , sυ .
(13)
The evolution process of the MRT-LB equation (6) consists of two steps: the collision step and streaming step [27]. Usually, the collision step is implemented in the moment space: m
+
=m−3 m−m
(eq)
3
+ δt I −
2
S,
(14)
where m+ represent the velocity moments after collision, while the streaming step is still carried out in the velocity space: fi (x + ei δt , t + δt ) = fi + (x, t ) ,
(15)
where f+ = M−1 m+ . In order to derive the correct equations of hydrodynamics, the forcing term S in the moment space should be chosen appropriately. For the D2Q8 MRT-LB model, the forcing term S in the moment space can be chosen as 0 4u · F / (5φ) F x − Fx , S = ρ0 Fy −Fy 2 u F − u F /φ y y x x ux Fy + uy Fx /φ
(16)
where F = Fx , Fy is given by Eq. (3). The fluid velocity u is defined as
ρ0 u =
8
ei fi +
δt 2
ρ0 F.
(17)
Note that the total body force F also contains the velocity u, so Eq. (17) is a nonlinear equation for u. According to Eqs. (3) and (17), the velocity u can be calculated explicitly by v
u= l0 +
30
(11)
c2 ρ0 uy ,
i =1 27
,
c1 ρ0 ux ,
23
T
To get the correct generalized Navier–Stokes equations (1) and (2), the parameters are chosen as follows: α2 = −1, c1 = c2 = −2, γ2 = 0, and γ1 = γ3 = 2/3. In the above equilibrium moments, the incompressibility approximation have been employed i.e., the fluid density ρ = ρ0 + δρ ≈ ρ0 (δρ is the density fluctuation), and J = jx , jy ≈ ρ0 u. The relaxation matrix 3 is given by
14 15
φ
3
, e(eq) , jx , q(xeq) , jy , q(yeq) , p(xxeq) , p(xyeq)
where
4
10
2 ρ0 |u|2
l20 + l1 |v|
,
(18)
where v, l0 and l1 are given by v=
8 i =1
ei fi /ρ0 +
δt 2
φ a,
l0 =
1 2
1+φ
δt υ 2 K
,
l1 = φ
δt Fφ √ . 2
K
(19)
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
5
According to Ref. [38], the fluid pressure p can be determined by p=
8
cs2
φ (1 − ω0 )
1
fi + ρ0 s0 (u) ,
(20)
where ω0 = 4/9 and s0 (u) = −2ρ0 |u|2 / (3φ). Here, si (u) = ωi 3ei · u + 4.5 (ei · u)2 /φ − 1.5|u|2 /φ with ω0 = 4/9, ω1∼4 = 1/9 and ω5∼8 = 1/36. The effective kinetic viscosity υe in the model is defined as
υe = cs2 τυ −
1 2
δt
with s7 = s8 = sυ = 1/τυ . The equilibrium distribution function fi velocity space is given by fi
(eq)
= ρ0 ωi
2
i =1
ei · u |u|2 φp (ei · u)2 + 2 + − , 2 4 ρ0 cs cs 2φ cs 2φ cs2
(eq)
i = 1 ∼ 8.
3 4 5
(21)
6
corresponding to the distribution function fi in the
7 8
(22)
The above expression is useful when implementing the boundary conditions. Through the Chapman–Enskog analysis [32,35,39] of the MRT-LB equation (6) in the moment space, the generalized Navier–Stokes equations (1) and (2) can be recovered exactly in the incompressible limit (see Appendix A for details). It should be noted that, as φ → 1 and Da → ∞, the present MRT-LB model reduces to the incompressible MRT-LB model [32] for free-fluid flows without a porous medium. When Fφ = 0, the Brinkman-extended Darcy equation can be obtained from the present MRT-LB model. 3. Numerical tests
9
10 11 12 13 14 15
16
In this section, numerical simulations of several typical two-dimensional porous flows are carried out to validate the proposed MRT-LB model. The testing problems include the Poiseuille flow in a channel filled with porous media, the Couette flow between two parallel plates filled with porous media, the lid-driven flow in a square porous cavity, and the natural convection flow in a square porous cavity. The present numerical results are compared with the analytical solutions and/or other numerical results in previous studies. In simulations, we set ρ0 = 1, c = 1, and δt = δx = δy = 1. The relaxation rates {si |1 ≤ i ≤ 8} are chosen as follows: s1 = s3 = s5 = 1, s2 = 1.1, s4 = s6 = 1.2, and s7 = s8 = sυ = 1/τυ . Unless otherwise stated, the non-equilibrium extrapolation scheme [40] is employed to treat the velocity and temperature boundary conditions in all computations. 3.1. Poiseuille flow in a channel filled with porous media
18 19 20 21 22 23 24
25
The first test problem is the Poiseuille flow in a 2D channel filled with porous media. The height of the channel is L, and the flow is driven by an external force a = (ax , 0) along the channel direction. When the flow is fully developed in the x-direction (along the channel direction), the governing equation of the flow can be expressed as
υe ∂ ux υ Fφ − ux − √ u2x + ax = 0 2 φ ∂y K K
17
26 27 28
2
(23)
with ux (x, 0) = ux (x, L) = 0, and the y-direction velocity component uy is zero everywhere. The Reynolds number is defined by Re = Lu0 /υ , where u0 is the maximum velocity of the porous Poiseuille flow (without the nonlinear drag force, i.e., Fφ = 0) along the centerline of the channel [22]: ax K
1 − cosh−1 υ √ where ϑ = φυ/ (K υe ). u0 =
ϑL 2
,
(24)
In computations, the porosity φ is set to be 0.1, the viscosity ratio J is set to be 1, the Darcy number Da changes from 10−6 to 10−3 , and the Reynolds number Re changes from 0.01 to 100. The relaxation rate sυ is set to be 5/3 (τυ = 0.6) with a Nx × Ny = 80 × 80 square mesh. The nonequilibrium extrapolation scheme is employed to the top and bottom plates for no-slip velocity boundary condition, and the periodic boundary conditions are employed to the inlet and outlet of the channel. At t = 0, the velocity moments are set to be their equilibrium values, and the velocity field is initialized to be 0. The velocity profiles predicted by the present MRT-LB model for different Re and Da are shown in Fig. 2. The numerical solutions given by Guo and Zhao [22] using the FD method are also included in the figure for comparison. From Fig. 2 it is clearly seen that the present numerical results agree well with the FD solutions reported in the literature.
29
30 31 32
33
34 35 36 37 38 39 40 41 42
6
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
Fig. 2. Velocity profiles of the porous Poiseuille flow for different Re and Da with φ = 0.1 and J = 1. Solid lines represent the FD solutions [22] and symbols represent the present MRT-LB results.
1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
17
3.2. Couette flow between two parallel plates filled with porous media We now apply the MRT-LB model to simulate the Couette flow in a 2D channel filled with fluid-saturated porous media. The flow is driven by the top plate moving along the x-direction with a uniform velocity u0 , while the bottom plate is fixed. The Reynolds number of this flow is defined by Re = Lu0 /υ . At the steady state, the flow still follows Eq. (24) (ax = 0) with ux (x, 0) = 0 and ux (x, L) = u0 . In simulations, the porosity φ is set to be 0.1, the viscosity ratio J is set to be 1, and the boundary and initial conditions are the same as those in the porous Poiseuille flow. For, 10, 50, and 1000, the kinetic viscosity υ is set to be 0.4, 0.08, 0.032, and 0.0008 based on a Nx × Ny = 80 × 80 square mesh, respectively. In Fig. 3, the velocity profiles predicted by the present MRT-LB model are presented for various Re and Da. The FD solutions in Ref. [22] are also included in Fig. 3 for comparison. Excellent agreement can be found between the present MRT-LB results and the FD solutions. The MRT-LB model is also applied to simulate the modified Couette flow in a 2D channel partially filled with a fluidsaturated porous medium, which has been investigated in Refs. [21,22,41]. A porous layer is positioned in the lower part of the channel such that there is a gap between the top plate and the medium. The porosity φ is set to be 0.1 for 0 ≤ y/L ≤ 1/2 (porous region) and 1 for 1/2 < y/L ≤ 1 (free-fluid region). For small Re and Da, the nonlinear drag force can be neglected, and the Brinkman-extended Darcy model is applied in the porous region. At the steady state, the analytical solution of the flow can be expressed as [22,41]: ux (y) =
√
18 19 20 21
d0 + d1 (y/L − 1/2) d0 exp [ϑ (y/L − 1/2)]
1/2 < y/L ≤ 1 0 ≤ y/L ≤ 1/2,
(25)
where ϑ = φυ/ (K υe ), d0 = 2ϑ Ku0 / (2ϑ K + φ), and d1 = 2φ u0 / (2ϑ K + φ). In simulations, Fφ is set to be zero, and the fluid kinetic viscosity υ is set to be 0.16 with a Nx × Ny = 80 × 80 lattice. The relaxation rate sυ is determined by sυ = 1/ 0.5 + J υ/ cs2 δt (for free-fluid region, J = 1). The velocity profiles predicted by the present MRT-LB model for J = 1 and 4 at Da = 0.001 and Re = 0.01 are plotted in Fig. 4. As shown, the present numerical results agree well with the
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
7
Fig. 3. Velocity profiles of the porous Couette flow for different Re and Da with φ = 0.1 and J = 1. Solid lines represent the FD solutions [22] and symbols represent the present MRT-LB results.
Fig. 4. Velocity profiles of the porous Couette flow for different viscosity ratio J with Re = 0.01 and Da = 0.001. Solid lines represent the analytical solutions from Eq. (25) and symbols represent the present MRT-LB results.
analytical ones. The discontinuity of the velocity-gradient at the porous medium/free-fluid interface is well captured by the present MRT-LB model without any special treatment for the boundary condition at the interface in simulations.
1 2
8
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
(a) Da = 10−2 .
(b) Da = 10−3 .
(c) Da = 10−4 .
Fig. 5. Streamlines of the lid-driven flow in a porous cavity for different Da with φ = 0.1 and Re = 10: (a) Da = 10−2 ; (b) Da = 10−3 ; (c) Da = 10−4 .
1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
3.3. Lid-driven flow in a square porous cavity In this subsection, we apply the MRT-LB model to simulate the incompressible flow in a square cavity filled with a fluidsaturated porous medium. The top wall of the cavity moves from left to right with a uniform velocity u0 , while the bottom, right, and left walls of the cavity are fixed. The Reynolds number of this flow is defined by Re = Lu0 /υ , where L is the height of the cavity. In simulations, we set φ = 0.1, Re = 10, J = 1, and the kinetic viscosity υ is set to be 0.032 based on a Nx × Ny = 128 × 128 square mesh. The velocity moments are set to be their equilibrium values, and the velocity field is initialized to be 0 at t = 0. The streamlines predicted by the present MRT-LB model for different Darcy numbers with φ = 0.1 and Re = 10 are shown in Fig. 5. From the figure we can see that, as Da decreases, the vortex inside the cavity becomes weaker, and the boundary layer near the top wall becomes thinner. This can be attributed to the lower permeability (K = DaL2 ) of the porous medium which results in lower fluid velocity inside the cavity. To be more informative, the horizontal velocity component ux in the vertical midplane x = L/2 and the vertical velocity component uy in the horizontal midplane y = L/2 are plotted in Fig. 6. The FD solutions in Ref. [22] are also plotted in Fig. 6 for comparison. Obviously, the present numerical results are in good agreement with the FD solutions. As φ → 1 and Da → ∞, the present MRT-LB model reduces to the incompressible MRT-LB model for free-fluid flows without a porous medium. We now present the numerical results for the cases in which φ → 1 and Da → ∞ with 100 ≤ Re ≤ 1000. In simulations, we set φ = 0.9999, Da = 108 , and υ = 0.0128. In Fig. 7, the velocity profiles through the center of the cavity are plotted for different Reynolds numbers. The benchmark solutions of Ghia et al. [42] are also included in Fig. 7 for comparison. Excellent agreement can be found between the present MRT-LB results and the benchmark solutions.
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
9
Fig. 6. Velocity profiles through the center of the cavity: (a) horizontal velocity component ux in the vertical midplane x = L/2; (b) vertical velocity component uy in the horizontal midplane y = L/2. Solid lines represent the FD solutions [22] and symbols represent the present MRT-LB results.
In what follows, the comparison between the present MRT-LB model and the BGK-LB model [22] is made through studying the lid-driven flow in a porous cavity at low viscosity based on a Nx × Ny = 128 × 128 lattice. In simulations, we set Re = 10, Da = 0.0001, φ = 0.1, and J = 1. The computation results are presented in Fig. 8. It can be seen from Fig. 8(c) that the Contour lines of the vertical velocity component obtained by the BGK model [22] show severe oscillations at υ = 5.12 × 10−4 (u0 = 4 × 10−5 ), but the MRT model can give stable results under the same condition (see Fig. 8(b)). The comparison indicates that the MRT model is more stable than the BGK model when the viscosity approaches zero. The MRT model can significantly improve the numerical stability of the LB schemes by separating the relaxation rates of the conserved (hydrodynamic) and nonconserved (kinetic) moments. 3.4. Natural convection flow in a square porous cavity Natural convection flow in a square cavity filled with a porous medium (see Fig. 9) has been studied extensively by many researchers [4,7,43,44] based on the generalized model. For this problem, the velocity field is solved by the present MRTLB model, while the temperature field is solved by the D2Q5 MRT-LB model [44] (see Appendix B for details). The top and bottom walls of the cavity are adiabatic, while the left and right vertical walls are kept at constant but different temperatures Th and Tc (Th > Tc ), respectively. The Prandtl number Pr and Rayleigh number Ra of this flow are defined as Pr = υ/α and Ra = g β ∆TL3 Pr /υ 2 , respectively, where ∆T = Th − Tc is the temperature difference, α is the thermal diffusivity of the fluid, g is the gravitational acceleration, β is the thermal expansion coefficient, and L is the distance between the walls. Based on the Boussinesq approximation, the body force (buoyancy force) a is given by a = g β (T − T0 ) j, where T0 = (Th + Tc ) /2 is the reference temperature, and j is the unit vector in the y-direction. In simulations, we set Pr = 1, J = 1, σ = 1 (thermal capacity ratio), γ = αe /α = 1 (αe is the effective thermal diffusivity of the porous media), ζ0 = 1, ζ1 = ζ2 = 1/τT , ζ3 = ζ4 = 1.1, Th = 21, and Tc = 1. According to Refs. [13,45],
1 2 3 4 5 6 7 8
9
10 11 12 13 14 15 16 17 18 19 20
10
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
Fig. 7. Velocity profiles through the center of the cavity: (a) horizontal velocity component ux in the vertical midplane x = L/2; (b) vertical velocity component uy in the horizontal midplane y = L/2. Solid lines represent the benchmark solutions [42] and symbols represent the present MRT-LB results. 1
the dimensionless relaxation times of the velocity and temperature fields can be determined by
√ 2
τυ =
1 2
+
MaJL 3 Pr c δt
√
Ra
,
τT = √
1 2
+
γ cs2 (τυ − 0.5) , 2 Pr Jσ csT
(26)
√
12
2 respectively, where Ma = uc /cs = 3uc is the Mach number (uc = β g ∆TL is the characteristic velocity), and csT = 2c 2 /5 (csT is the sound speed of the D2Q5 model). For incompressible fluids considered in this study, Ma is set to be 0.15 in all computations. The grid sizes of 120 × 120, 200 × 200, and 250 × 250 are employed for Da = 10−2 , 10−4 , and 10−6 , respectively. Fig. 10 illustrates the streamlines and isotherms for different Rayleigh numbers and Darcy numbers (Darcy–Rayleigh number Ra∗ = RaDa = 1000) with Pr = 1 and φ = 0.6. From the figure we can observe that, as Da decreases, the velocity and thermal boundary layers become thinner near the hot and cold vertical walls. As Da increases to 10−2 , the streamlines and isotherms are less crowded near the vertical walls and more convective mixing occurs inside the cavity. To quantify the results, the average Nusselt numbers of the hot wall are calculated and listed in Table 1. The numerical results obtained by Nithiarasu et al. [4,7] using the FE method are also included in Table 1 for comparison. To sum up, the numerical results of the present MRT-LB model agree well with those results reported in previous studies.
13
4. Conclusion
3 4 5 6 7 8 9 10 11
14 15 16 17 18
In this paper, an MRT-LB model with the eight-by-eight collision matrix has been presented for simulating incompressible flows in porous media at the REV scale. The key point of the MRT-LB model is to include the porosity into the pressure-based equilibrium moments, and add a forcing term to the MRT-LB equation in the moment space to account for the linear and nonlinear drag forces of the solid matrix based on the generalized model. Through the Chapman–Enskog analysis in the moment space, the generalized Navier–Stokes equations can be recovered exactly. Numerical simulations of the porous Poiseuille
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
(a) MRT, υ = 0.256.
11
(b) MRT, υ = 5.12 × 10−4 .
(c) BGK, υ = 5.12 × 10−4 . Fig. 8. Contour lines of the vertical velocity component uy for the lid-driven flow in a porous cavity: (a) MRT, υ = 0.256; (b) MRT, υ = 5.12 × 10−4 ; (c) BGK, υ = 5.12 × 10−4 .
Fig. 9. Natural convection flow in a square cavity filled with a porous medium.
flow, porous Couette flow, lid-driven flow in square porous cavity, and natural convection flow in a square porous cavity have been carried out to demonstrate the present MRT-LB model. The present numerical results agree well with the analytical solutions and/or other numerical solutions reported in previous studies. Moreover, the comparison between the MRT
1 2 3
12
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
(a) Da = 10−2 , Ra = 105 .
(b) Da = 10−4 , Ra = 107 .
(c) Da = 10−6 , Ra = 109 . Fig. 10. Streamlines and isotherms for different Ra and Da with Pr = 1 and φ = 0.6: (a) Da = 10−2 , Ra = 105 ; (b) Da = 10−4 , Ra = 107 ; (c) Da = 10−6 , Ra = 109 . Table 1 Comparisons of the average Nusselt numbers for different Ra, Da and φ with Pr = 1.0 (grid sizes: 120 × 120 for Da = 10−2 , 200 × 200 for Da = 10−4 , and 250 × 250 for Da = 10−6 ). Da
10
−2
10−4
10−6
1 2
φ = 0.4
Ra
3
10 104 105 105 106 107 107 108 109
φ = 0.6
φ = 0.9
Ref. [4]
Ref. [7]
Present
Ref. [4]
Ref. [7]
Present
Ref. [4]
Ref. [7]
Present
1.010 1.408 2.983 1.067 2.550 7.810 1.079 2.970 11.460
1.008 1.359 2.986 1.064 2.580 7.677 1.074 2.969 11.699
1.008 1.364 3.005 1.065 2.605 7.770 1.076 3.020 11.476
1.015 1.530 3.555 1.071 2.725 8.183 1.079 2.997 11.790
1.012 1.489 3.430 1.066 2.686 8.452 1.074 2.982 12.098
1.012 1.495 3.451 1.068 2.714 8.540 1.076 3.037 11.878
1.023 1.640 3.910 1.072 2.740 9.202 1.08 3.00 12.01
– – – – – – – – –
1.018 1.637 3.930 1.070 2.798 9.244 1.077 3.041 12.142
and BGK models has also been made. It is found that the stability of the MRT-LB model is much better than the BGK-LB model when the viscosity approaches zero. This feature makes the present MRT-LB model more useful in practical applications.
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
13
Acknowledgment
1
This work was supported by the National Key Basic Research Program of China (973 Program) (2013CB228304).
2
Appendix A. Chapman–Enskog analysis of the D2Q8 MRT-LB model
3
The Chapman–Enskog expansion method [32,35,39] is adopted to derive the generalized Navier–Stokes equations (1) and (2) from the D2Q8 MRT-LB model. To this end, the following expansions in time and space are introduced: fi (x + ei δt , t + δt ) =
∞
ϵn (∂t + ei · ∇)n fi (x, t ) , n!
4 5
(A.1a)
6
(A.1b)
7
(A.1c)
8
where ϵ = δt is a small expansion parameter, S1 = (S11 , . . . , S81 )T , F1 = Fx1 , Fy1 . With the above expansions, we can derive the following equations from Eq. (6) as consecutive orders of the parameter ϵ in the moment space as
10
i =0
fi = fi
(0)
+ ϵ fi(1) + ϵ 2 fi(2) + · · · ,
∂t = ϵ∂t1 + ϵ ∂t2 , 2
∂j = ϵ∂j1 ,
S = ϵ S1 ,
F = ϵ F1 ,
ϵ 0 : m(0) = m(eq) ,
3
(A.2a)
11
(A.2b)
12
(A.2c)
13
ϵ 1 : D˜ 1 m(0) = −3′ m(1) + I − S1 , 2 δt ˜ 3 3 m(1) + D I − S1 = −3′ m(2) , ϵ 2 : ∂t2 m(0) + D˜ 1 I − 1 2
2
2
˜ 1 = MD1 M−1 = ∂t I + Cj ∂j1 (j = x, y), D1 = ∂t I + ∂j1 diag e1j , e2j , . . . , e8j , Cj = M eij I M−1 , 3′ = 3/δt , and where D 1 1
(1) T
m(1) = 0, e(1) , −δt ρ0 Fx1 /2, q(x1) , −δt ρ0 Fy1 /2, q(y1) , p(xx1) , pxy
.
(A.3)
Cx and Cy can be given explicitly by
0
0
0 3 4 0 Cx = 0 0 0
0
0
0
1 4 1 0 0 0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0 0 0
0 0 0
−1
0 0 1
0 0
0 0 0 1
0
0
0
2
1
0 0 0 3 Cy = 4 0 0
3
3
0
0
0 2
3
3
0 0 0 0 1 3 0
−
3
0
1 2 0 0
0
0 , 1 1 0 0
0
4 1 0
1 1
0 2
0 0
3 0 0
3 0 0
0
0
0
0
0
0
0 1
0 1
3
3
0
0
2
1
3
3
−
0
0
0
0
0
0
0 0 1
1 1
∂ t1
5 3
φp +
2 ρ0 |u| 3
φ
2
15
−
2 1 0 0
0 .
0 0
(A.4)
+ ∂x1 (ρ0 ux ) + ∂y1 ρ0 uy = 0,
1 1 s2 ∂t1 (−φ p) − ∂x1 (ρ0 ux ) − ∂y1 ρ0 uy = −s′2 e(1) + 1 − S21 , 3 3 2 ρ0 u2x ρ0 ux uy δt s3 ∂t1 (ρ0 ux ) + ∂x1 φ p + + ∂y1 = s′3 ρ0 Fx1 + 1 − S31 , φ φ 2 2 ρ0 u2x − u2y ρ0 ux uy s4 ∂t1 (−ρ0 ux ) + ∂x1 −φ p − + ∂y1 = −s′4 q(x1) + 1 − S41 , φ φ 2 ρ0 u2y ρ 0 ux uy δt s5 ∂t1 ρ0 uy + ∂x1 + ∂y1 φ p + = s′5 ρ0 Fy1 + 1 − S51 , φ φ 2 2 ρ0 u2x − u2y ρ 0 ux uy s6 ∂t1 −ρ0 uy + ∂x1 + ∂y1 −φ p + = −s′6 q(y1) + 1 − S61 , φ φ 2
17
0
From Eq. (A.2b), the following equations at the t1 time scale can be obtained:
14
16
0
1 1
9
18
(A.5a)
19
(A.5b)
20
(A.5c)
21
(A.5d)
22
(A.5e)
23
(A.5f)
24
14
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
1
∂t1
2
∂t1
3
4
5
∂t2
10
11 12
13 14
15
16
17
18 19
20
22
23
5 3
+
3
3
(A.5g)
2
1 1 s8 S81 . + ∂x1 ρ0 uy + ∂y1 (ρ0 ux ) = −s′8 p(xy1) + 1 − 3
3
(A.5h)
2
φp +
δt
2
2 ρ0 |u|2
φ
3
δt
∂ t1
2
∂t1
= 0, s3
1−
1−
2
s3 2
(A.6a)
ρ0 Fx1 + ∂x1
S31 + ∂x1
1 4
1
1−
4
1−
s2 2
s2 (1) 1 s7 (1) s8 (1) e + 1− pxx + ∂y1 1 − pxy 2 2 2 2
S21 +
s8 1 s7 = 0, 1− S71 + ∂y1 1 − S81 2 2 2
(A.6b)
δt s5 s2 (1) 1 s7 (1) s8 (1) 1 ∂t2 ρ0 uy − ∂t1 1 − ρ0 Fy1 + ∂x1 1 − pxy + ∂y1 1− e − 1− pxx 2
+
δt
2
∂t1
2
1−
s5 2
2
S51 + ∂x1
1−
s8 2
4
S81 + ∂y1
1 4
1−
s2 2
2
S21 −
2
1
s7
2
2
1−
2
S71
= 0,
(A.6c)
Following the approach in Ref. [32], we can obtain:
∂t1
5 3
φp +
2 ρ0 |u|2
φ
3
= 0,
(A.7)
Combining Eqs. (A.5a) and (A.6a) with Eq. (A.7) leads to the following incompressible continuity equation
∂x ux + ∂y uy = 0,
(A.8)
(1)
(1)
Note that e(1) , pxx and pxy in Eqs. (A.6b) and (A.6c) are unknowns to be determined. According to Eqs. (A.5b), (A.5g) and (A.5h), we can get:
1 1 s2 −s′2 e(1) = ∂t1 (−φ p) − ∂x1 (ρ0 ux ) − ∂y1 ρ0 uy − 1 − S21 , 3 3 2 ρ0 u2x − u2y 2 2 s7 ′ (1) −s7 pxx = ∂t1 + ∂x1 (ρ0 ux ) − ∂y1 ρ0 uy − 1 − S71 , φ 3 3 2 1 1 s8 ρ0 ux uy ′ (1) + ∂x1 ρ0 uy + ∂y1 (ρ0 ux ) − 1 − −s8 pxy = ∂t1 S81 . φ 3 3 2
(A.9a)
(A.9b)
(A.9c)
Neglecting the terms of order O(|u|3 ) and higher-order terms of the form uj ∂k (uk uj ), using Eqs. (A.5c) and (A.5e), we can obtain:
∂t1
21
ρ 0 ux uy φ
∂t2 (ρ0 ux ) −
8
9
φ
2 2 s7 S71 , + ∂x1 (ρ0 ux ) − ∂y1 ρ0 uy = −s′7 p(xx1) + 1 −
From Eq. (A.2c), the following equations at the t2 time scale can be obtained:
6
7
ρ0 u2x − u2y
∂t1 ∂t1
ρ0 u2x φ
ρ0 u2y
=
=
φ
ρ 0 ux uy φ
2ρ0 ux Fx1
,
(A.10a)
,
(A.10b)
ρ0 ux Fy1 + uy Fx1 = . φ
(A.10c)
φ 2ρ0 uy Fy1
φ
With the above equations, we can obtain: ′ (1)
24
−s 2 e
=
25
−s′7 p(xx1) =
26
−s′8 p(xy1)
4ρ0 15
2s2 ρ0 ux Fx1 + uy Fy1
∂x1 ux + ∂y1 uy +
5
φ
ρ0 ux Fx1 − uy Fy1 ∂x1 ux − ∂y1 uy + s7 , 3 φ s8 ρ0 ux Fy1 + uy Fx1 ρ0 ∂x1 uy + ∂y1 ux + = . 3 2 φ 2ρ0
,
(A.11a)
(A.11b)
(A.11c)
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
15
Substituting Eq. (A.11) into Eq. (A.6), the following equations at the t2 time scale can be derived:
∂ t 2 ux −
δt
15
δt
1 s2
1
1
−
2
1
δt
∂x1 ∂x1 ux + ∂y1 uy −
3
1 s7
−
1 2
2
∂y1 ∂x1 uy + ∂y1 ux = 0, 3 s8 2 δt 1 δt 1 1 1 ∂ t 2 uy − − ∂x1 ∂x1 uy + ∂y1 ux − − ∂y1 ∂x1 ux + ∂y1 uy −
3
+
δt
−
s8
1
2
1
−
1
∂x1 ∂x1 ux − ∂y1 uy
15
s2
(A.12a)
4
2
∂y1 ∂x1 ux − ∂y1 uy = 0.
(A.12b) 3 s7 2 Combining Eq. (A.12) (t2 time scale) with Eqs. (A.5c) and (A.5e) (t1 time scale), the following equations can be derived (∂t = ϵ∂t1 + ϵ 2 ∂t2 ):
∂ t ux + ∂ x
2
ux
φ
+ ∂y
ux uy
δt 1 1 ∂x ∂x ux + ∂y uy ∂x (φ p) + − ρ0 2 15 s2 δt 1 δt 1 1 1 + − ∂x ∂x ux − ∂y uy + − ∂y ∂x uy + ∂y ux + Fx ,
= −
φ
1
∂t uy + ∂x
ux uy
φ
+ ∂y
u2y
6 7
8
δt 1 1 ∂y (φ p) + − ∂x ∂x uy + ∂y ux ρ0 3 s8 2 δt 1 δt 1 1 1 + − ∂y ∂x ux + ∂y uy − − ∂y ∂x ux − ∂y uy + Fy . (A.13b)
9
s7
2
3
s8
2
1
= −
φ
5
(A.13a)
3
3
15
s2
2
3
s7
2
With the aid of the incompressible continuity equation ∂x ux + ∂y uy = 0, the incompressible generalized Navier–Stokes equations (1) and (2) can be obtained from Eq. (A.13). The effective kinetic viscosity is defined as υe = cs2 (τυ − 0.5) δt with s7 = s8 = sυ = 1/τυ . Appendix B. D2Q5 MRT-LB model
10 11 12
13
The macroscopic temperature equation of natural convection flow in porous media can be written as [4,7,43]
14
∂T + u · ∇ T = ∇ · (αe ∇ T ) , (B.1) ∂t where σ is the thermal capacity ratio, and αe is the effective thermal diffusivity of the porous media. For the temperature σ
field, the D2Q5 MRT-LB equation is given by g (x + eδt , t + δt ) − g (x, t ) = −N
16 17
(eq)
2 n (x, t ) − n (B.2) where N is a 5 × 5 orthogonal transformation matrix, and 2 is a diagonal relaxation matrix. The boldface symbols, g, n, and −1
15
(x, t ) ,
n(eq) are 5-dimensional (column) vectors:
18 19 20
g (x, t ) = (g0 (x, t ) , g1 (x, t ) , . . . , g4 (x, t )) ,
21
g (x + eδt , t + δt ) = (g0 (x + e0 δt , t + δt ) , . . . , g4 (x + e4 δt , t + δt ))T ,
22
n (x, t ) = (n0 (x, t ) , n1 (x, t ) , . . . , n4 (x, t ))T ,
23
T
(eq)
n(eq) (x, t ) = n0
(x, t ) , n(1
eq)
(x, t ) , . . . , n(4
eq)
(x, t )
T
,
24
where gi (x, t ) is the temperature distribution function, n (x, t ) and n(eq) (x, t ) are the velocity moments and the corresponding equilibrium values, respectively. The five discrete velocities {ei |i = 0, 1, . . . , 4} of the D2Q5 model are given by
ei =
(0, 0) , (cos [(i − 1) π/2] , sin [(i − 1) π /2]) c ,
i=0 i = 1 ∼ 4.
(B.3)
The transformation matrix N is given by
1
0
1 1 0 1 1
1 0 1 1 −1
1
−1 0 1 1
1 0 −1 , 1 −1
where |1⟩ = (1, 1, 1, 1, 1)T , |ex ⟩ = (0, 1, 0, −1, 0)T , and |ey ⟩ = (0, 0, 1, 0, −1)T .
27
28
T N = |1⟩, |ex ⟩, |ey ⟩, |5e2 − 4⟩, |e2x − e2y ⟩ 0 N= 0 −4
25 26
29
(B.4)
30
31
16
1 2 3 4 5
6
7 8 9
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
The transformation matrix N linearly transforms the discrete distribution functions g ∈ V = R5 to their velocity moments n ∈ M = R5 : n = Ng,
g = N −1 n .
(B.5)
In the system of {gi }, only the temperature T ≡ n0 =
4
i =0
(eq)
gi is conserved quantity. The equilibrium moments ni
|
i = 0, 1, . . . , 4 for the velocity moments {n|i = 0, 1, . . . , 4} are defined as follows [44]:
(eq)
(eq)
ux T
uy T
(eq) (eq) (eq) , n2 = , n3 = ϖ T , n4 = 0, (B.6) σ σ where ϖ is a constant (−4 < ϖ < 1). In the present work, ϖ is set to be 0. The diagonal relaxation matrix 2 is given by: 2 = diag (1, ζ1 , ζ2 , ζ3 , ζ4 ) . (B.7) 2 2 2 The effective thermal diffusivity αe is given by αe = σ csT (τT − 0.5) δt with ζ1 = ζ2 = 1/τT and csT = (4 + ϖ ) c /10.
n0
= T,
n1
=
References [1] K. Vafai, Handbook of Porous Media, Marcel Dekker, New York, 2000. [2] I. Pop, D.B. Ingham, Convective Heat Transfer: Mathematical and Computational Modeling of Viscous Fluids and Porous Media, Pergamon, Oxford, 2001. [3] D.A. Nield, A. Bejan, Convection in Porous Media, third ed., Springer, New York, 2006. [4] P. Nithiarasu, K.N. Seetharamu, T. Sundararajan, Natural convective heat transfer in a fluid saturated variable porosity medium, Int. J. Heat Mass Transfer 40 (1997) 3955–3967. [5] K. Vafai, C.L. Tien, Boundary and inertia effects on flow and heat transfer in porous media, Int. J. Heat Mass Transfer 24 (1981) 195–203. [6] C.T. Hsu, P. Cheng, Thermal dispersion in a porous medium, Int. J. Heat Mass Transfer 33 (1990) 1587–1597. [7] P. Nithiarasu, K. Ravindran, A new semi-implicit time stepping procedure for buoyancy driven flow in a fluid saturated porous medium, Comput. Methods Appl. Mech. Engrg. 165 (1998) 147–154. [8] K. Vafai, S.J. Kim, On the limitations of Brinkman–Forchheimer-extended Darcy equation, Int. J. Heat Fluid Flow 16 (1995) 11–15. [9] U. Frisch, B. Hasslacher, Y. Pomeau, Lattice-gas automata for the Navier–Stokes equation, Phys. Rev. Lett. 56 (1986) 1505. [10] R. Benzi, S. Succi, M. Vergassola, The lattice Boltzmann equation: theory and applications, Phys. Rep. 222 (1992) 145–197. [11] S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329–364. [12] Y. Wang, Y.L. He, T.S. Zhao, G.H. Tang, W.Q. Tao, Implicit-explicit finite-difference lattice Boltzmann method for compressible flows, Internat. J. Modern Phys. C 18 (2007) 1961–1983. [13] Y.L. He, Y. Wang, Q. Li, Lattice Boltzmann Method: Theory and Applications, Science Press, Beijing, 2009. [14] Q. Li, K.H. Luo, X.J. Li, Lattice Boltzmann method for relativistic hydrodynamics: issues on conservation law of particle number and discontinuities, Phys. Rev. D 86 (2012) 085044. [15] Y.L. He, Q. Liu, Q. Li, Three-dimensional finite-difference lattice Boltzmann model and its application to inviscid compressible flows with shock waves, Physica A 392 (2013) 4884–4896. [16] Y.X. Xia, Y.H. Qian, Lattice Boltzmann simulation for forced two-dimensional turbulence, Phys. Rev. E 90 (2) (2014) 023004. [17] S. Succi, Lattice Boltzmann across scales: from turbulence to DNA translocation, Eur. Phys. J. B 64 (2008) 471–479. [18] M.A.A. Spaid, F.R. Phelan Jr., Lattice Boltzmann methods for modeling microscale flow in fibrous porous media, Phys. Fluids 9 (1997) 2468–2474. [19] O. Dardis, J. McCloskey, Lattice Boltzmann scheme with real numbered solid density for the simulation of flow in porous media, Phys. Rev. E 57 (1998) 4834–4837. [20] D.M. Freed, Lattice-Boltzmann method for macroscopic porous media modeling, Internat. J. Modern Phys. C 9 (1998) 1491–1503. [21] N.S. Martys, Improved approximation of the Brinkman equation using a lattice Boltzmann method, Phys. Fluids 6 (2001) 1807–1810. [22] Z. Guo, T.S. Zhao, Lattice Boltzmann model for incompressible flows through porous media, Phys. Rev. E 66 (2002) 036304. [23] Q. Kang, D. Zhang, S. Chen, Unified lattice Boltzmann method for flow in multiscale porous media, Phys. Rev. E 66 (2002) 056307. [24] S. Anwar, M.C. Sukop, Regional scale transient groundwater flow modeling using lattice Boltzmann methods, Comput. Math. Appl. 58 (2009) 1015–1023. [25] S.D.C. Walsh, H. Burwinkle, M.O. Saar, A new partial-bounceback lattice-Boltzmann method for fluid flow through heterogeneous media, Comput. Geosci. 35 (2009) 1186–1193. [26] P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases. I. small amplitude processes in charged and neutral one-component systems, Phys. Rev. 94 (1954) 511–525. [27] P. Lallemand, L.-S. Luo, Theory of the lattice Boltzmann method: dispersion, dissipation, isotropy, Galilean invariance, and stability, Phys. Rev. E 61 (2000) 6546. [28] D. d’Humières, I. Ginzburg, M. Krafczyk, P. Lallemand, L.-S. Luo, Multiple-relaxation-time lattice Boltzmann models in three dimensions, Philos. Trans. R. Soc. Lond. Ser. A 360 (2002) 437–451. [29] Q. Li, Y.L. He, G.H. Tang, W.Q. Tao, Improved axisymmetric lattice Boltzmann scheme, Phys. Rev. E 81 (2010) 056707. [30] I. Ginzburg, D. d’Humières, Multireflection boundary conditions for lattice Boltzmann models, Phys. Rev. E 68 (2003) 066614. [31] D. d’Humières, Generalized lattice-Boltzmann equations, in: B.D. Shizgal, D.P. Weaver (Eds.), Rarefied Gas Dynamics: Theory and Simulations, in: Prog. Astronaut. Aeronaut., vol. 159, AIAA, Washington, DC, 1992, pp. 450–458. [32] R. Du, B. Shi, Incompressible MRT lattice Boltzmann model with eight velocities in 2D space, Internat. J. Modern Phys. C 20 (2009) 1023–1037. [33] S. Ergun, Fluid flow through packed columns, Chem. Eng. Prog. 48 (1952) 89–94. [34] K. Vafai, Convective flow and heat transfer in variable-porosity media, J. Fluid Mech. 147 (1984) 233–259. [35] M.E. MaCracken, J. Abraham, Multiple-relaxation-time lattice-Boltzmann model for multiphase flow, Phys. Rev. E 71 (2005) 036701. [36] Q. Li, K.H. Luo, X.J. Li, Lattice Boltzmann modeling of multiphase flows at large density ratio with an improved pseudopotential model, Phys. Rev. E 87 (2013) 053301. [37] Y.H. Qian, D. d’Humeier, P. Lallemand, Lattice BGK models for Navier–Stokes equation, Europhys. Lett. 17 (1992) 479–484. [38] Z. Guo, B. Shi, N. Wang, Lattice BGK model for incompressible Navier–Stokes equation, J. Comput. Phys. 165 (2000) 288–306. [39] Z. Chai, T.S. Zhao, Effect of the forcing term in the multiple-relaxation-time lattice Boltzmann equation on the shear stress or the strain rate tensor, Phys. Rev. E 86 (2012) 016705. [40] Z.L. Guo, C.G. Zheng, B.C. Shi, Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method, Chin. Phys. 11 (2002) 366. [41] N. Martys, D.P. Bentz, E.J. Garboczi, Computer simulation study of the effective viscosity in Brinkman’s equation, Phys. Fluids 6 (1994) 1434–1439. [42] U. Ghia, K.N. Ghia, C.T. Shin, High-Re solutions for incompressible flow using the Navier–Stokes equations and a multigrid method, J. Comput. Phys. 48 (1982) 387–411. [43] Z. Guo, T.S. Zhao, A lattice Boltzmann model for convection heat transfer in porous media, Numer. Heat Transfer B 47 (2005) 157–177. [44] Q. Liu, Y.L. He, Q. Li, W.Q. Tao, A multiple-relaxation-time lattice Boltzmann model for convection heat transfer in porous media, Int. J. Heat Mass Transfer 73 (2014) 761–775. [45] Q. Li, Y.L. He, Y. Wang, G.H. Tang, An improved thermal lattice Boltzmann model for flows without viscous heat dissipation and compression work, Internat. J. Modern Phys. C 19 (2008) 125–150.
Contents lists available at ScienceDirect
Physica A journal homepage: www.elsevier.com/locate/physa
Multiple-relaxation-time lattice Boltzmann modeling of incompressible flows in porous media Q1
Qing Liu, Ya-Ling He ∗ Key Laboratory of Thermo-Fluid Science and Engineering of MOE, School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi, 710049, PR China
highlights • • • •
A D2Q8 MRT-LB model is proposed for simulating incompressible porous flows at the REV scale. The generalized non-Darcy model is employed to describe the momentum transfer in porous media. The generalized Navier–Stokes equations can be recovered through the Chapman–Enskog analysis in the moment space. The MRT-LB model is demonstrated by numerical simulations of several typical two-dimensional porous flows.
article
info
Article history: Received 21 October 2014 Received in revised form 12 January 2015 Available online xxxx Keywords: Lattice Boltzmann model Multiple-relaxation-time Porous media Incompressible flows REV scale
abstract In this paper, a two-dimensional eight-velocity multiple-relaxation-time (MRT) lattice Boltzmann (LB) model is proposed for incompressible porous flows at the representative elementary volume scale based on the Brinkman–Forchheimer-extended Darcy model. In the model, the porosity is included into the pressure-based equilibrium moments, and the linear and nonlinear drag forces of the porous matrix are incorporated into the model by adding a forcing term to the MRT-LB equation in the moment space. Through the Chapman–Enskog analysis, the incompressible generalized Navier–Stokes equations can be recovered. Numerical simulations of several typical porous flows are carried out to validate the present MRT-LB model. It is found that the present numerical results agree well with the analytical solutions and/or other numerical results reported in the literature. © 2015 Published by Elsevier B.V.
1. Introduction Fluid flow and related transport phenomena in porous media have gained significant research interest due to the importance of related technological and industrial applications, which include contaminant transport in groundwater, crude oil exploration and extraction, radioactive waste management, hydrogeology and so on [1–3]. For incompressible flows in porous media at the representative elementary volume (REV) scale, the Darcy model, the Brinkman-extended Darcy model and the Forchheimer-extended Darcy model have been widely employed. However, the Darcy model and the two extended (Brinkman and Forchheimer) models have some intrinsic limitations in simulating porous flows [4,5]. In order to overcome the deficiencies of the above mentioned models, the Brinkman–Forchheimer-extended Darcy model (also called the generalized model) has been developed by several research groups [4–7]. In the generalized model, the viscous and inertial forces are incorporated into the momentum equation by using the local volume-averaging technique at the REV scale. The Darcy model and the two extended models can be regarded as the limiting cases of the generalized model. Due to the
∗
Corresponding author. E-mail address: [email protected] (Y.-L. He).
http://dx.doi.org/10.1016/j.physa.2015.01.067 0378-4371/© 2015 Published by Elsevier B.V.
1
2 3 4 5 6 7 8 9 10
2
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
similarity with the Navier–Stokes equations, the generalized model can be used to simulate transient flows through porous media. Moreover, as reported by Vafai and Kim [8], numerical results based on the Brinkman–Forchheimer-extended Darcy formulation have been shown to be in good agreement with the experimental predictions. In the past several decades, various traditional numerical methods, such as the finite volume (FV) method, the finite difference (FD) method, and the finite element (FE) method, have been employed to study porous flows based on the generalized non-Darcy model. The lattice Boltzmann (LB) method, as a bran-new mesoscopic numerical technique originates from the lattice-gas automata (LGA) method [9], has achieved significant success in modeling complex fluid flows and simulating complex physics in fluids [10–16]. Owing to its kinetic background, the LB method has some distinctive advantages over the traditional numerical methods (e.g., see Ref. [17]). Recently, the LB method has been successfully applied to simulate fluid flows in porous media at the REV scale [18–25]. In the REV scale method, the detailed geometric structure of the media is ignored, and the standard LB equation is modified by adding an additional term to account for the presence of the porous media. Therefore, LB method at the REV scale can be used for systems with large computational domain. It is worth mentioning that the REV LB method is very effective for simulating fluid flows in the region which is partially filled with a porous medium. As reported in Ref. [21], the discontinuity of the velocity-gradient at the porous medium/free-fluid interface can be well captured by the LB method without including the stress boundary condition into the simulation model. To the best of our knowledge, most of the existing REV LB models [18–25] for incompressible porous flows employ the Q2 Bhatnagar–Gross–Krook (BGK) model (also called the single-relaxation-time model) [26] to represent the collision process. Although the BGK-LB model has become the most popular one in the LB community because of its extreme simplicity, there are several well-known criticisms on this model, such as the numerical instability at low values of viscosity [27–29] and the inaccuracy in treating boundary conditions [30]. Fortunately, these shortcomings of the BGK-LB model can be easily addressed by employing the multiple-relaxation-time (MRT) model proposed by d’Humières [31]. Hence, the aim of this paper is to develop a new MRT-LB model for incompressible porous flows at the REV scale based on the generalized model. In the model, a pressure-based MRT-LB equation with the eight-by-eight collision matrix [32] is constructed in the framework of the standard MRT-LB method. The remainder of this paper is organized as follows. In Section 2, the MRT-LB model for incompressible porous flows at the REV scale is presented in detail. In Section 3, numerical tests of the MRT-LB model are performed for the porous Poiseuille flow, porous Couette flow, lid-driven flow in a square porous cavity, and natural convection flow in a square porous cavity. Finally, a brief conclusion is made in Section 4.
28
2. MRT-LB model for incompressible flows in porous media
29
2.1. Macroscopic governing equations
30 31 32 33
The fluid flow is assumed to be two-dimensional, laminar and incompressible. For incompressible flows through porous media at the REV scale, the generalized model proposed by Nithiarasu et al. [4] is employed in the present study. The dimensional governing equations (generalized Navier–Stokes equations) of the generalized model can be written as follows:
∇ · u = 0,
(1)
36
u 1 ∂u + (u · ∇) = − ∇ (φ p) + υe ∇ 2 u + F, (2) ∂t φ ρ0 fluid velocity and pressure, respectively, φ is the where ρ0 is the average fluid density, u and p are the volume-averaged porosity, and υe is the effective kinetic viscosity. F = Fx , Fy denotes the total body force induced by the porous matrix and
37
other external forces, which can be expressed as [6,22]
34
35
38
39 40 41 42
43
44 45 46
47
48
F=−
φυ
φ Fφ u − √ |u|u + φ a, K K
(3)
where υ is the kinetic viscosity of the fluid, K is the permeability, Fφ is the geometric function, a is the body force due to an external force, and |u| =
u2x + u2y , in which ux and uy are the x- and y-components of the macroscopic velocity u,
respectively. Based on Ergun’s relation [33], the geometric function Fφ and the permeability K of the porous media can be expressed as [34] 1.75 Fφ = , 150φ 3
K =
φ 3 d2p 150 (1 − φ)2
,
(4)
where dp is the diameter of the solid particle. The flow governed by the generalized Navier–Stokes equations (1) and (2) are characterized by the porosity φ and several dimensionless parameters: the Darcy number Da, the viscosity ratio J, and the Reynolds number Re, which are defined as K υe LU , J = , Re = , L2 υ υ where L and U are the characteristic length and velocity of the system, respectively. Da =
(5)
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
3
Fig. 1. Discrete velocities of the D2Q8 model.
2.2. MRT-LB model
1
In this subsection, a two-dimensional eight-velocity (D2Q8) MRT-LB model is presented to study incompressible porous flows. According to Refs. [35,36], the MRT-LB equation with an explicit treatment of the forcing term can be written as f (x + eδt , t + δt ) − f (x, t ) = −M−1 3 m (x, t ) − m(eq) (x, t ) + δt M−1
I−
3 2
S,
(6)
where M is the 8 × 8 orthogonal transformation matrix, 3 is the non-negative 8 × 8 diagonal relaxation matrix, and I is the 8 × 8 identity matrix. The boldface symbols, f, m, m(eq) , and S represent 8-dimensional (column) vectors: f (x, t ) = (f1 (x, t ) , f2 (x, t ) , . . . , f8 (x, t )) , T
2 3
4
5 6 7
f (x + eδt , t + δt ) = (f1 (x + e1 δt , t + δt ) , . . . , f8 (x + e8 δt , t + δt )) ,
8
m (x, t ) = (m1 (x, t ) , m2 (x, t ) , . . . , m8 (x, t ))T ,
9
T
m
(eq)
(eq)
(x, t ) = m1
(eq)
(x, t ) , m2
(eq)
(x, t ) , . . . , m8
T (x, t ) ,
10
S = (S1 , S2 , . . . , S8 )T ,
11
where T is the transpose operator, fi (x, t ) is the discrete distribution function with velocity ei at position x and time t (at the REV scale), m (x, t ) and m(eq) (x, t ) are the velocity moments and the corresponding equilibrium moments, respectively, and {Si |i = 1, 2, . . . , 8} are components of the forcing term S. In the D2Q8 model, the eight discrete velocities {ei |i = 1, 2, . . . , 8} (see Fig. 1) are given by [32,37]
ei =
(cos [(i − 1) π /2] , sin [(i − 1) π /2]) c , √ (cos [(2i − 9) π /4] , sin [(2i − 9) π /4]) 2c ,
i=1∼4
(7)
i = 5 ∼ 8,
where c = δx /δt is the lattice speed,√ in which δt and δx are the discrete time step and lattice spacing, respectively. The sound speed of the D2Q8 model is cs = c / 3. In the MRT model, the lattice speed c is usually set to be 1, which leads to δx = δt . The transformation matrix M linearly transforms the discrete distribution functions f ∈ V = R8 (velocity space) to their velocity moments m ∈ M = R8 (moment space): m = Mf,
f = M−1 m.
(8)
For the D2Q8 model, the eight velocity moments {mi |i = 1, 2, . . . , 8} are: m = (m1 , m2 , m3 , m4 , m5 , m6 , m7 , m8 )
= p, e, jx , qx , jy , qy , pxx , pxy
T
1 −1 0 0 1 −2 −1 0
1 −1 −1 2 0 0 1 0
1 −1 0 0 −1 2 −1 0
1 1 1 1 1 1 0 1
,
(9)
1 1 −1 −1 1 1 0 −1
1 1 −1 −1 −1 −1 0 1
14 15
16
17 18 19 20
21
23
where m1 = p is the pressure, m2 = e is related to energy, m3,5 = jx,y are components of the momentum J = jx , jy , m4,6 = qx,y are related to energy flux, and m7,8 = pxx,xy are related to the diagonal and off-diagonal components of the stress tensor, respectively [27]. With the ordering of the velocity moments {mi }, the transformation matrix M can be constructed via the Gram–Schmidt orthogonalization procedure (c = 1) [32]: 1 −1 1 −2 M= 0 0 1 0
13
22
T
12
1 1 1 1 . −1 −1 0 −1
24
25 26 27 28
(10)
29
4
1
2
3
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
The equilibrium moments m(eq) corresponding to the velocity moments m are defined by m(eq) =
5 3
φp +
e(eq) = α2 φ p +
5
q(xeq) =
6
q(yeq) =
7
p(xxeq) =
8
p(xyeq) =
9
11 12 13
16
17
18 19 20 21 22
1 2 1 2 3 2 3 2
3
γ2
4
ρ0 u2x + u2y φ
24 25
26
28
29
31
,
γ1 γ3
ρ0 u2x − u2y φ
,
ρ0 ux uy . φ
(12)
3 = diag (s1 , s2 , s3 , s4 , s5 , s6 , s7 , s8 )
= diag 1, se , 1, sq , 1, sq , sυ , sυ .
(13)
The evolution process of the MRT-LB equation (6) consists of two steps: the collision step and streaming step [27]. Usually, the collision step is implemented in the moment space: m
+
=m−3 m−m
(eq)
3
+ δt I −
2
S,
(14)
where m+ represent the velocity moments after collision, while the streaming step is still carried out in the velocity space: fi (x + ei δt , t + δt ) = fi + (x, t ) ,
(15)
where f+ = M−1 m+ . In order to derive the correct equations of hydrodynamics, the forcing term S in the moment space should be chosen appropriately. For the D2Q8 MRT-LB model, the forcing term S in the moment space can be chosen as 0 4u · F / (5φ) F x − Fx , S = ρ0 Fy −Fy 2 u F − u F /φ y y x x ux Fy + uy Fx /φ
(16)
where F = Fx , Fy is given by Eq. (3). The fluid velocity u is defined as
ρ0 u =
8
ei fi +
δt 2
ρ0 F.
(17)
Note that the total body force F also contains the velocity u, so Eq. (17) is a nonlinear equation for u. According to Eqs. (3) and (17), the velocity u can be calculated explicitly by v
u= l0 +
30
(11)
c2 ρ0 uy ,
i =1 27
,
c1 ρ0 ux ,
23
T
To get the correct generalized Navier–Stokes equations (1) and (2), the parameters are chosen as follows: α2 = −1, c1 = c2 = −2, γ2 = 0, and γ1 = γ3 = 2/3. In the above equilibrium moments, the incompressibility approximation have been employed i.e., the fluid density ρ = ρ0 + δρ ≈ ρ0 (δρ is the density fluctuation), and J = jx , jy ≈ ρ0 u. The relaxation matrix 3 is given by
14 15
φ
3
, e(eq) , jx , q(xeq) , jy , q(yeq) , p(xxeq) , p(xyeq)
where
4
10
2 ρ0 |u|2
l20 + l1 |v|
,
(18)
where v, l0 and l1 are given by v=
8 i =1
ei fi /ρ0 +
δt 2
φ a,
l0 =
1 2
1+φ
δt υ 2 K
,
l1 = φ
δt Fφ √ . 2
K
(19)
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
5
According to Ref. [38], the fluid pressure p can be determined by p=
8
cs2
φ (1 − ω0 )
1
fi + ρ0 s0 (u) ,
(20)
where ω0 = 4/9 and s0 (u) = −2ρ0 |u|2 / (3φ). Here, si (u) = ωi 3ei · u + 4.5 (ei · u)2 /φ − 1.5|u|2 /φ with ω0 = 4/9, ω1∼4 = 1/9 and ω5∼8 = 1/36. The effective kinetic viscosity υe in the model is defined as
υe = cs2 τυ −
1 2
δt
with s7 = s8 = sυ = 1/τυ . The equilibrium distribution function fi velocity space is given by fi
(eq)
= ρ0 ωi
2
i =1
ei · u |u|2 φp (ei · u)2 + 2 + − , 2 4 ρ0 cs cs 2φ cs 2φ cs2
(eq)
i = 1 ∼ 8.
3 4 5
(21)
6
corresponding to the distribution function fi in the
7 8
(22)
The above expression is useful when implementing the boundary conditions. Through the Chapman–Enskog analysis [32,35,39] of the MRT-LB equation (6) in the moment space, the generalized Navier–Stokes equations (1) and (2) can be recovered exactly in the incompressible limit (see Appendix A for details). It should be noted that, as φ → 1 and Da → ∞, the present MRT-LB model reduces to the incompressible MRT-LB model [32] for free-fluid flows without a porous medium. When Fφ = 0, the Brinkman-extended Darcy equation can be obtained from the present MRT-LB model. 3. Numerical tests
9
10 11 12 13 14 15
16
In this section, numerical simulations of several typical two-dimensional porous flows are carried out to validate the proposed MRT-LB model. The testing problems include the Poiseuille flow in a channel filled with porous media, the Couette flow between two parallel plates filled with porous media, the lid-driven flow in a square porous cavity, and the natural convection flow in a square porous cavity. The present numerical results are compared with the analytical solutions and/or other numerical results in previous studies. In simulations, we set ρ0 = 1, c = 1, and δt = δx = δy = 1. The relaxation rates {si |1 ≤ i ≤ 8} are chosen as follows: s1 = s3 = s5 = 1, s2 = 1.1, s4 = s6 = 1.2, and s7 = s8 = sυ = 1/τυ . Unless otherwise stated, the non-equilibrium extrapolation scheme [40] is employed to treat the velocity and temperature boundary conditions in all computations. 3.1. Poiseuille flow in a channel filled with porous media
18 19 20 21 22 23 24
25
The first test problem is the Poiseuille flow in a 2D channel filled with porous media. The height of the channel is L, and the flow is driven by an external force a = (ax , 0) along the channel direction. When the flow is fully developed in the x-direction (along the channel direction), the governing equation of the flow can be expressed as
υe ∂ ux υ Fφ − ux − √ u2x + ax = 0 2 φ ∂y K K
17
26 27 28
2
(23)
with ux (x, 0) = ux (x, L) = 0, and the y-direction velocity component uy is zero everywhere. The Reynolds number is defined by Re = Lu0 /υ , where u0 is the maximum velocity of the porous Poiseuille flow (without the nonlinear drag force, i.e., Fφ = 0) along the centerline of the channel [22]: ax K
1 − cosh−1 υ √ where ϑ = φυ/ (K υe ). u0 =
ϑL 2
,
(24)
In computations, the porosity φ is set to be 0.1, the viscosity ratio J is set to be 1, the Darcy number Da changes from 10−6 to 10−3 , and the Reynolds number Re changes from 0.01 to 100. The relaxation rate sυ is set to be 5/3 (τυ = 0.6) with a Nx × Ny = 80 × 80 square mesh. The nonequilibrium extrapolation scheme is employed to the top and bottom plates for no-slip velocity boundary condition, and the periodic boundary conditions are employed to the inlet and outlet of the channel. At t = 0, the velocity moments are set to be their equilibrium values, and the velocity field is initialized to be 0. The velocity profiles predicted by the present MRT-LB model for different Re and Da are shown in Fig. 2. The numerical solutions given by Guo and Zhao [22] using the FD method are also included in the figure for comparison. From Fig. 2 it is clearly seen that the present numerical results agree well with the FD solutions reported in the literature.
29
30 31 32
33
34 35 36 37 38 39 40 41 42
6
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
Fig. 2. Velocity profiles of the porous Poiseuille flow for different Re and Da with φ = 0.1 and J = 1. Solid lines represent the FD solutions [22] and symbols represent the present MRT-LB results.
1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
17
3.2. Couette flow between two parallel plates filled with porous media We now apply the MRT-LB model to simulate the Couette flow in a 2D channel filled with fluid-saturated porous media. The flow is driven by the top plate moving along the x-direction with a uniform velocity u0 , while the bottom plate is fixed. The Reynolds number of this flow is defined by Re = Lu0 /υ . At the steady state, the flow still follows Eq. (24) (ax = 0) with ux (x, 0) = 0 and ux (x, L) = u0 . In simulations, the porosity φ is set to be 0.1, the viscosity ratio J is set to be 1, and the boundary and initial conditions are the same as those in the porous Poiseuille flow. For, 10, 50, and 1000, the kinetic viscosity υ is set to be 0.4, 0.08, 0.032, and 0.0008 based on a Nx × Ny = 80 × 80 square mesh, respectively. In Fig. 3, the velocity profiles predicted by the present MRT-LB model are presented for various Re and Da. The FD solutions in Ref. [22] are also included in Fig. 3 for comparison. Excellent agreement can be found between the present MRT-LB results and the FD solutions. The MRT-LB model is also applied to simulate the modified Couette flow in a 2D channel partially filled with a fluidsaturated porous medium, which has been investigated in Refs. [21,22,41]. A porous layer is positioned in the lower part of the channel such that there is a gap between the top plate and the medium. The porosity φ is set to be 0.1 for 0 ≤ y/L ≤ 1/2 (porous region) and 1 for 1/2 < y/L ≤ 1 (free-fluid region). For small Re and Da, the nonlinear drag force can be neglected, and the Brinkman-extended Darcy model is applied in the porous region. At the steady state, the analytical solution of the flow can be expressed as [22,41]: ux (y) =
√
18 19 20 21
d0 + d1 (y/L − 1/2) d0 exp [ϑ (y/L − 1/2)]
1/2 < y/L ≤ 1 0 ≤ y/L ≤ 1/2,
(25)
where ϑ = φυ/ (K υe ), d0 = 2ϑ Ku0 / (2ϑ K + φ), and d1 = 2φ u0 / (2ϑ K + φ). In simulations, Fφ is set to be zero, and the fluid kinetic viscosity υ is set to be 0.16 with a Nx × Ny = 80 × 80 lattice. The relaxation rate sυ is determined by sυ = 1/ 0.5 + J υ/ cs2 δt (for free-fluid region, J = 1). The velocity profiles predicted by the present MRT-LB model for J = 1 and 4 at Da = 0.001 and Re = 0.01 are plotted in Fig. 4. As shown, the present numerical results agree well with the
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
7
Fig. 3. Velocity profiles of the porous Couette flow for different Re and Da with φ = 0.1 and J = 1. Solid lines represent the FD solutions [22] and symbols represent the present MRT-LB results.
Fig. 4. Velocity profiles of the porous Couette flow for different viscosity ratio J with Re = 0.01 and Da = 0.001. Solid lines represent the analytical solutions from Eq. (25) and symbols represent the present MRT-LB results.
analytical ones. The discontinuity of the velocity-gradient at the porous medium/free-fluid interface is well captured by the present MRT-LB model without any special treatment for the boundary condition at the interface in simulations.
1 2
8
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
(a) Da = 10−2 .
(b) Da = 10−3 .
(c) Da = 10−4 .
Fig. 5. Streamlines of the lid-driven flow in a porous cavity for different Da with φ = 0.1 and Re = 10: (a) Da = 10−2 ; (b) Da = 10−3 ; (c) Da = 10−4 .
1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
3.3. Lid-driven flow in a square porous cavity In this subsection, we apply the MRT-LB model to simulate the incompressible flow in a square cavity filled with a fluidsaturated porous medium. The top wall of the cavity moves from left to right with a uniform velocity u0 , while the bottom, right, and left walls of the cavity are fixed. The Reynolds number of this flow is defined by Re = Lu0 /υ , where L is the height of the cavity. In simulations, we set φ = 0.1, Re = 10, J = 1, and the kinetic viscosity υ is set to be 0.032 based on a Nx × Ny = 128 × 128 square mesh. The velocity moments are set to be their equilibrium values, and the velocity field is initialized to be 0 at t = 0. The streamlines predicted by the present MRT-LB model for different Darcy numbers with φ = 0.1 and Re = 10 are shown in Fig. 5. From the figure we can see that, as Da decreases, the vortex inside the cavity becomes weaker, and the boundary layer near the top wall becomes thinner. This can be attributed to the lower permeability (K = DaL2 ) of the porous medium which results in lower fluid velocity inside the cavity. To be more informative, the horizontal velocity component ux in the vertical midplane x = L/2 and the vertical velocity component uy in the horizontal midplane y = L/2 are plotted in Fig. 6. The FD solutions in Ref. [22] are also plotted in Fig. 6 for comparison. Obviously, the present numerical results are in good agreement with the FD solutions. As φ → 1 and Da → ∞, the present MRT-LB model reduces to the incompressible MRT-LB model for free-fluid flows without a porous medium. We now present the numerical results for the cases in which φ → 1 and Da → ∞ with 100 ≤ Re ≤ 1000. In simulations, we set φ = 0.9999, Da = 108 , and υ = 0.0128. In Fig. 7, the velocity profiles through the center of the cavity are plotted for different Reynolds numbers. The benchmark solutions of Ghia et al. [42] are also included in Fig. 7 for comparison. Excellent agreement can be found between the present MRT-LB results and the benchmark solutions.
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
9
Fig. 6. Velocity profiles through the center of the cavity: (a) horizontal velocity component ux in the vertical midplane x = L/2; (b) vertical velocity component uy in the horizontal midplane y = L/2. Solid lines represent the FD solutions [22] and symbols represent the present MRT-LB results.
In what follows, the comparison between the present MRT-LB model and the BGK-LB model [22] is made through studying the lid-driven flow in a porous cavity at low viscosity based on a Nx × Ny = 128 × 128 lattice. In simulations, we set Re = 10, Da = 0.0001, φ = 0.1, and J = 1. The computation results are presented in Fig. 8. It can be seen from Fig. 8(c) that the Contour lines of the vertical velocity component obtained by the BGK model [22] show severe oscillations at υ = 5.12 × 10−4 (u0 = 4 × 10−5 ), but the MRT model can give stable results under the same condition (see Fig. 8(b)). The comparison indicates that the MRT model is more stable than the BGK model when the viscosity approaches zero. The MRT model can significantly improve the numerical stability of the LB schemes by separating the relaxation rates of the conserved (hydrodynamic) and nonconserved (kinetic) moments. 3.4. Natural convection flow in a square porous cavity Natural convection flow in a square cavity filled with a porous medium (see Fig. 9) has been studied extensively by many researchers [4,7,43,44] based on the generalized model. For this problem, the velocity field is solved by the present MRTLB model, while the temperature field is solved by the D2Q5 MRT-LB model [44] (see Appendix B for details). The top and bottom walls of the cavity are adiabatic, while the left and right vertical walls are kept at constant but different temperatures Th and Tc (Th > Tc ), respectively. The Prandtl number Pr and Rayleigh number Ra of this flow are defined as Pr = υ/α and Ra = g β ∆TL3 Pr /υ 2 , respectively, where ∆T = Th − Tc is the temperature difference, α is the thermal diffusivity of the fluid, g is the gravitational acceleration, β is the thermal expansion coefficient, and L is the distance between the walls. Based on the Boussinesq approximation, the body force (buoyancy force) a is given by a = g β (T − T0 ) j, where T0 = (Th + Tc ) /2 is the reference temperature, and j is the unit vector in the y-direction. In simulations, we set Pr = 1, J = 1, σ = 1 (thermal capacity ratio), γ = αe /α = 1 (αe is the effective thermal diffusivity of the porous media), ζ0 = 1, ζ1 = ζ2 = 1/τT , ζ3 = ζ4 = 1.1, Th = 21, and Tc = 1. According to Refs. [13,45],
1 2 3 4 5 6 7 8
9
10 11 12 13 14 15 16 17 18 19 20
10
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
Fig. 7. Velocity profiles through the center of the cavity: (a) horizontal velocity component ux in the vertical midplane x = L/2; (b) vertical velocity component uy in the horizontal midplane y = L/2. Solid lines represent the benchmark solutions [42] and symbols represent the present MRT-LB results. 1
the dimensionless relaxation times of the velocity and temperature fields can be determined by
√ 2
τυ =
1 2
+
MaJL 3 Pr c δt
√
Ra
,
τT = √
1 2
+
γ cs2 (τυ − 0.5) , 2 Pr Jσ csT
(26)
√
12
2 respectively, where Ma = uc /cs = 3uc is the Mach number (uc = β g ∆TL is the characteristic velocity), and csT = 2c 2 /5 (csT is the sound speed of the D2Q5 model). For incompressible fluids considered in this study, Ma is set to be 0.15 in all computations. The grid sizes of 120 × 120, 200 × 200, and 250 × 250 are employed for Da = 10−2 , 10−4 , and 10−6 , respectively. Fig. 10 illustrates the streamlines and isotherms for different Rayleigh numbers and Darcy numbers (Darcy–Rayleigh number Ra∗ = RaDa = 1000) with Pr = 1 and φ = 0.6. From the figure we can observe that, as Da decreases, the velocity and thermal boundary layers become thinner near the hot and cold vertical walls. As Da increases to 10−2 , the streamlines and isotherms are less crowded near the vertical walls and more convective mixing occurs inside the cavity. To quantify the results, the average Nusselt numbers of the hot wall are calculated and listed in Table 1. The numerical results obtained by Nithiarasu et al. [4,7] using the FE method are also included in Table 1 for comparison. To sum up, the numerical results of the present MRT-LB model agree well with those results reported in previous studies.
13
4. Conclusion
3 4 5 6 7 8 9 10 11
14 15 16 17 18
In this paper, an MRT-LB model with the eight-by-eight collision matrix has been presented for simulating incompressible flows in porous media at the REV scale. The key point of the MRT-LB model is to include the porosity into the pressure-based equilibrium moments, and add a forcing term to the MRT-LB equation in the moment space to account for the linear and nonlinear drag forces of the solid matrix based on the generalized model. Through the Chapman–Enskog analysis in the moment space, the generalized Navier–Stokes equations can be recovered exactly. Numerical simulations of the porous Poiseuille
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
(a) MRT, υ = 0.256.
11
(b) MRT, υ = 5.12 × 10−4 .
(c) BGK, υ = 5.12 × 10−4 . Fig. 8. Contour lines of the vertical velocity component uy for the lid-driven flow in a porous cavity: (a) MRT, υ = 0.256; (b) MRT, υ = 5.12 × 10−4 ; (c) BGK, υ = 5.12 × 10−4 .
Fig. 9. Natural convection flow in a square cavity filled with a porous medium.
flow, porous Couette flow, lid-driven flow in square porous cavity, and natural convection flow in a square porous cavity have been carried out to demonstrate the present MRT-LB model. The present numerical results agree well with the analytical solutions and/or other numerical solutions reported in previous studies. Moreover, the comparison between the MRT
1 2 3
12
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
(a) Da = 10−2 , Ra = 105 .
(b) Da = 10−4 , Ra = 107 .
(c) Da = 10−6 , Ra = 109 . Fig. 10. Streamlines and isotherms for different Ra and Da with Pr = 1 and φ = 0.6: (a) Da = 10−2 , Ra = 105 ; (b) Da = 10−4 , Ra = 107 ; (c) Da = 10−6 , Ra = 109 . Table 1 Comparisons of the average Nusselt numbers for different Ra, Da and φ with Pr = 1.0 (grid sizes: 120 × 120 for Da = 10−2 , 200 × 200 for Da = 10−4 , and 250 × 250 for Da = 10−6 ). Da
10
−2
10−4
10−6
1 2
φ = 0.4
Ra
3
10 104 105 105 106 107 107 108 109
φ = 0.6
φ = 0.9
Ref. [4]
Ref. [7]
Present
Ref. [4]
Ref. [7]
Present
Ref. [4]
Ref. [7]
Present
1.010 1.408 2.983 1.067 2.550 7.810 1.079 2.970 11.460
1.008 1.359 2.986 1.064 2.580 7.677 1.074 2.969 11.699
1.008 1.364 3.005 1.065 2.605 7.770 1.076 3.020 11.476
1.015 1.530 3.555 1.071 2.725 8.183 1.079 2.997 11.790
1.012 1.489 3.430 1.066 2.686 8.452 1.074 2.982 12.098
1.012 1.495 3.451 1.068 2.714 8.540 1.076 3.037 11.878
1.023 1.640 3.910 1.072 2.740 9.202 1.08 3.00 12.01
– – – – – – – – –
1.018 1.637 3.930 1.070 2.798 9.244 1.077 3.041 12.142
and BGK models has also been made. It is found that the stability of the MRT-LB model is much better than the BGK-LB model when the viscosity approaches zero. This feature makes the present MRT-LB model more useful in practical applications.
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
13
Acknowledgment
1
This work was supported by the National Key Basic Research Program of China (973 Program) (2013CB228304).
2
Appendix A. Chapman–Enskog analysis of the D2Q8 MRT-LB model
3
The Chapman–Enskog expansion method [32,35,39] is adopted to derive the generalized Navier–Stokes equations (1) and (2) from the D2Q8 MRT-LB model. To this end, the following expansions in time and space are introduced: fi (x + ei δt , t + δt ) =
∞
ϵn (∂t + ei · ∇)n fi (x, t ) , n!
4 5
(A.1a)
6
(A.1b)
7
(A.1c)
8
where ϵ = δt is a small expansion parameter, S1 = (S11 , . . . , S81 )T , F1 = Fx1 , Fy1 . With the above expansions, we can derive the following equations from Eq. (6) as consecutive orders of the parameter ϵ in the moment space as
10
i =0
fi = fi
(0)
+ ϵ fi(1) + ϵ 2 fi(2) + · · · ,
∂t = ϵ∂t1 + ϵ ∂t2 , 2
∂j = ϵ∂j1 ,
S = ϵ S1 ,
F = ϵ F1 ,
ϵ 0 : m(0) = m(eq) ,
3
(A.2a)
11
(A.2b)
12
(A.2c)
13
ϵ 1 : D˜ 1 m(0) = −3′ m(1) + I − S1 , 2 δt ˜ 3 3 m(1) + D I − S1 = −3′ m(2) , ϵ 2 : ∂t2 m(0) + D˜ 1 I − 1 2
2
2
˜ 1 = MD1 M−1 = ∂t I + Cj ∂j1 (j = x, y), D1 = ∂t I + ∂j1 diag e1j , e2j , . . . , e8j , Cj = M eij I M−1 , 3′ = 3/δt , and where D 1 1
(1) T
m(1) = 0, e(1) , −δt ρ0 Fx1 /2, q(x1) , −δt ρ0 Fy1 /2, q(y1) , p(xx1) , pxy
.
(A.3)
Cx and Cy can be given explicitly by
0
0
0 3 4 0 Cx = 0 0 0
0
0
0
1 4 1 0 0 0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0 0 0
0 0 0
−1
0 0 1
0 0
0 0 0 1
0
0
0
2
1
0 0 0 3 Cy = 4 0 0
3
3
0
0
0 2
3
3
0 0 0 0 1 3 0
−
3
0
1 2 0 0
0
0 , 1 1 0 0
0
4 1 0
1 1
0 2
0 0
3 0 0
3 0 0
0
0
0
0
0
0
0 1
0 1
3
3
0
0
2
1
3
3
−
0
0
0
0
0
0
0 0 1
1 1
∂ t1
5 3
φp +
2 ρ0 |u| 3
φ
2
15
−
2 1 0 0
0 .
0 0
(A.4)
+ ∂x1 (ρ0 ux ) + ∂y1 ρ0 uy = 0,
1 1 s2 ∂t1 (−φ p) − ∂x1 (ρ0 ux ) − ∂y1 ρ0 uy = −s′2 e(1) + 1 − S21 , 3 3 2 ρ0 u2x ρ0 ux uy δt s3 ∂t1 (ρ0 ux ) + ∂x1 φ p + + ∂y1 = s′3 ρ0 Fx1 + 1 − S31 , φ φ 2 2 ρ0 u2x − u2y ρ0 ux uy s4 ∂t1 (−ρ0 ux ) + ∂x1 −φ p − + ∂y1 = −s′4 q(x1) + 1 − S41 , φ φ 2 ρ0 u2y ρ 0 ux uy δt s5 ∂t1 ρ0 uy + ∂x1 + ∂y1 φ p + = s′5 ρ0 Fy1 + 1 − S51 , φ φ 2 2 ρ0 u2x − u2y ρ 0 ux uy s6 ∂t1 −ρ0 uy + ∂x1 + ∂y1 −φ p + = −s′6 q(y1) + 1 − S61 , φ φ 2
17
0
From Eq. (A.2b), the following equations at the t1 time scale can be obtained:
14
16
0
1 1
9
18
(A.5a)
19
(A.5b)
20
(A.5c)
21
(A.5d)
22
(A.5e)
23
(A.5f)
24
14
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
1
∂t1
2
∂t1
3
4
5
∂t2
10
11 12
13 14
15
16
17
18 19
20
22
23
5 3
+
3
3
(A.5g)
2
1 1 s8 S81 . + ∂x1 ρ0 uy + ∂y1 (ρ0 ux ) = −s′8 p(xy1) + 1 − 3
3
(A.5h)
2
φp +
δt
2
2 ρ0 |u|2
φ
3
δt
∂ t1
2
∂t1
= 0, s3
1−
1−
2
s3 2
(A.6a)
ρ0 Fx1 + ∂x1
S31 + ∂x1
1 4
1
1−
4
1−
s2 2
s2 (1) 1 s7 (1) s8 (1) e + 1− pxx + ∂y1 1 − pxy 2 2 2 2
S21 +
s8 1 s7 = 0, 1− S71 + ∂y1 1 − S81 2 2 2
(A.6b)
δt s5 s2 (1) 1 s7 (1) s8 (1) 1 ∂t2 ρ0 uy − ∂t1 1 − ρ0 Fy1 + ∂x1 1 − pxy + ∂y1 1− e − 1− pxx 2
+
δt
2
∂t1
2
1−
s5 2
2
S51 + ∂x1
1−
s8 2
4
S81 + ∂y1
1 4
1−
s2 2
2
S21 −
2
1
s7
2
2
1−
2
S71
= 0,
(A.6c)
Following the approach in Ref. [32], we can obtain:
∂t1
5 3
φp +
2 ρ0 |u|2
φ
3
= 0,
(A.7)
Combining Eqs. (A.5a) and (A.6a) with Eq. (A.7) leads to the following incompressible continuity equation
∂x ux + ∂y uy = 0,
(A.8)
(1)
(1)
Note that e(1) , pxx and pxy in Eqs. (A.6b) and (A.6c) are unknowns to be determined. According to Eqs. (A.5b), (A.5g) and (A.5h), we can get:
1 1 s2 −s′2 e(1) = ∂t1 (−φ p) − ∂x1 (ρ0 ux ) − ∂y1 ρ0 uy − 1 − S21 , 3 3 2 ρ0 u2x − u2y 2 2 s7 ′ (1) −s7 pxx = ∂t1 + ∂x1 (ρ0 ux ) − ∂y1 ρ0 uy − 1 − S71 , φ 3 3 2 1 1 s8 ρ0 ux uy ′ (1) + ∂x1 ρ0 uy + ∂y1 (ρ0 ux ) − 1 − −s8 pxy = ∂t1 S81 . φ 3 3 2
(A.9a)
(A.9b)
(A.9c)
Neglecting the terms of order O(|u|3 ) and higher-order terms of the form uj ∂k (uk uj ), using Eqs. (A.5c) and (A.5e), we can obtain:
∂t1
21
ρ 0 ux uy φ
∂t2 (ρ0 ux ) −
8
9
φ
2 2 s7 S71 , + ∂x1 (ρ0 ux ) − ∂y1 ρ0 uy = −s′7 p(xx1) + 1 −
From Eq. (A.2c), the following equations at the t2 time scale can be obtained:
6
7
ρ0 u2x − u2y
∂t1 ∂t1
ρ0 u2x φ
ρ0 u2y
=
=
φ
ρ 0 ux uy φ
2ρ0 ux Fx1
,
(A.10a)
,
(A.10b)
ρ0 ux Fy1 + uy Fx1 = . φ
(A.10c)
φ 2ρ0 uy Fy1
φ
With the above equations, we can obtain: ′ (1)
24
−s 2 e
=
25
−s′7 p(xx1) =
26
−s′8 p(xy1)
4ρ0 15
2s2 ρ0 ux Fx1 + uy Fy1
∂x1 ux + ∂y1 uy +
5
φ
ρ0 ux Fx1 − uy Fy1 ∂x1 ux − ∂y1 uy + s7 , 3 φ s8 ρ0 ux Fy1 + uy Fx1 ρ0 ∂x1 uy + ∂y1 ux + = . 3 2 φ 2ρ0
,
(A.11a)
(A.11b)
(A.11c)
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
15
Substituting Eq. (A.11) into Eq. (A.6), the following equations at the t2 time scale can be derived:
∂ t 2 ux −
δt
15
δt
1 s2
1
1
−
2
1
δt
∂x1 ∂x1 ux + ∂y1 uy −
3
1 s7
−
1 2
2
∂y1 ∂x1 uy + ∂y1 ux = 0, 3 s8 2 δt 1 δt 1 1 1 ∂ t 2 uy − − ∂x1 ∂x1 uy + ∂y1 ux − − ∂y1 ∂x1 ux + ∂y1 uy −
3
+
δt
−
s8
1
2
1
−
1
∂x1 ∂x1 ux − ∂y1 uy
15
s2
(A.12a)
4
2
∂y1 ∂x1 ux − ∂y1 uy = 0.
(A.12b) 3 s7 2 Combining Eq. (A.12) (t2 time scale) with Eqs. (A.5c) and (A.5e) (t1 time scale), the following equations can be derived (∂t = ϵ∂t1 + ϵ 2 ∂t2 ):
∂ t ux + ∂ x
2
ux
φ
+ ∂y
ux uy
δt 1 1 ∂x ∂x ux + ∂y uy ∂x (φ p) + − ρ0 2 15 s2 δt 1 δt 1 1 1 + − ∂x ∂x ux − ∂y uy + − ∂y ∂x uy + ∂y ux + Fx ,
= −
φ
1
∂t uy + ∂x
ux uy
φ
+ ∂y
u2y
6 7
8
δt 1 1 ∂y (φ p) + − ∂x ∂x uy + ∂y ux ρ0 3 s8 2 δt 1 δt 1 1 1 + − ∂y ∂x ux + ∂y uy − − ∂y ∂x ux − ∂y uy + Fy . (A.13b)
9
s7
2
3
s8
2
1
= −
φ
5
(A.13a)
3
3
15
s2
2
3
s7
2
With the aid of the incompressible continuity equation ∂x ux + ∂y uy = 0, the incompressible generalized Navier–Stokes equations (1) and (2) can be obtained from Eq. (A.13). The effective kinetic viscosity is defined as υe = cs2 (τυ − 0.5) δt with s7 = s8 = sυ = 1/τυ . Appendix B. D2Q5 MRT-LB model
10 11 12
13
The macroscopic temperature equation of natural convection flow in porous media can be written as [4,7,43]
14
∂T + u · ∇ T = ∇ · (αe ∇ T ) , (B.1) ∂t where σ is the thermal capacity ratio, and αe is the effective thermal diffusivity of the porous media. For the temperature σ
field, the D2Q5 MRT-LB equation is given by g (x + eδt , t + δt ) − g (x, t ) = −N
16 17
(eq)
2 n (x, t ) − n (B.2) where N is a 5 × 5 orthogonal transformation matrix, and 2 is a diagonal relaxation matrix. The boldface symbols, g, n, and −1
15
(x, t ) ,
n(eq) are 5-dimensional (column) vectors:
18 19 20
g (x, t ) = (g0 (x, t ) , g1 (x, t ) , . . . , g4 (x, t )) ,
21
g (x + eδt , t + δt ) = (g0 (x + e0 δt , t + δt ) , . . . , g4 (x + e4 δt , t + δt ))T ,
22
n (x, t ) = (n0 (x, t ) , n1 (x, t ) , . . . , n4 (x, t ))T ,
23
T
(eq)
n(eq) (x, t ) = n0
(x, t ) , n(1
eq)
(x, t ) , . . . , n(4
eq)
(x, t )
T
,
24
where gi (x, t ) is the temperature distribution function, n (x, t ) and n(eq) (x, t ) are the velocity moments and the corresponding equilibrium values, respectively. The five discrete velocities {ei |i = 0, 1, . . . , 4} of the D2Q5 model are given by
ei =
(0, 0) , (cos [(i − 1) π/2] , sin [(i − 1) π /2]) c ,
i=0 i = 1 ∼ 4.
(B.3)
The transformation matrix N is given by
1
0
1 1 0 1 1
1 0 1 1 −1
1
−1 0 1 1
1 0 −1 , 1 −1
where |1⟩ = (1, 1, 1, 1, 1)T , |ex ⟩ = (0, 1, 0, −1, 0)T , and |ey ⟩ = (0, 0, 1, 0, −1)T .
27
28
T N = |1⟩, |ex ⟩, |ey ⟩, |5e2 − 4⟩, |e2x − e2y ⟩ 0 N= 0 −4
25 26
29
(B.4)
30
31
16
1 2 3 4 5
6
7 8 9
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
Q. Liu, Y.-L. He / Physica A xx (xxxx) xxx–xxx
The transformation matrix N linearly transforms the discrete distribution functions g ∈ V = R5 to their velocity moments n ∈ M = R5 : n = Ng,
g = N −1 n .
(B.5)
In the system of {gi }, only the temperature T ≡ n0 =
4
i =0
(eq)
gi is conserved quantity. The equilibrium moments ni
|
i = 0, 1, . . . , 4 for the velocity moments {n|i = 0, 1, . . . , 4} are defined as follows [44]:
(eq)
(eq)
ux T
uy T
(eq) (eq) (eq) , n2 = , n3 = ϖ T , n4 = 0, (B.6) σ σ where ϖ is a constant (−4 < ϖ < 1). In the present work, ϖ is set to be 0. The diagonal relaxation matrix 2 is given by: 2 = diag (1, ζ1 , ζ2 , ζ3 , ζ4 ) . (B.7) 2 2 2 The effective thermal diffusivity αe is given by αe = σ csT (τT − 0.5) δt with ζ1 = ζ2 = 1/τT and csT = (4 + ϖ ) c /10.
n0
= T,
n1
=
References [1] K. Vafai, Handbook of Porous Media, Marcel Dekker, New York, 2000. [2] I. Pop, D.B. Ingham, Convective Heat Transfer: Mathematical and Computational Modeling of Viscous Fluids and Porous Media, Pergamon, Oxford, 2001. [3] D.A. Nield, A. Bejan, Convection in Porous Media, third ed., Springer, New York, 2006. [4] P. Nithiarasu, K.N. Seetharamu, T. Sundararajan, Natural convective heat transfer in a fluid saturated variable porosity medium, Int. J. Heat Mass Transfer 40 (1997) 3955–3967. [5] K. Vafai, C.L. Tien, Boundary and inertia effects on flow and heat transfer in porous media, Int. J. Heat Mass Transfer 24 (1981) 195–203. [6] C.T. Hsu, P. Cheng, Thermal dispersion in a porous medium, Int. J. Heat Mass Transfer 33 (1990) 1587–1597. [7] P. Nithiarasu, K. Ravindran, A new semi-implicit time stepping procedure for buoyancy driven flow in a fluid saturated porous medium, Comput. Methods Appl. Mech. Engrg. 165 (1998) 147–154. [8] K. Vafai, S.J. Kim, On the limitations of Brinkman–Forchheimer-extended Darcy equation, Int. J. Heat Fluid Flow 16 (1995) 11–15. [9] U. Frisch, B. Hasslacher, Y. Pomeau, Lattice-gas automata for the Navier–Stokes equation, Phys. Rev. Lett. 56 (1986) 1505. [10] R. Benzi, S. Succi, M. Vergassola, The lattice Boltzmann equation: theory and applications, Phys. Rep. 222 (1992) 145–197. [11] S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329–364. [12] Y. Wang, Y.L. He, T.S. Zhao, G.H. Tang, W.Q. Tao, Implicit-explicit finite-difference lattice Boltzmann method for compressible flows, Internat. J. Modern Phys. C 18 (2007) 1961–1983. [13] Y.L. He, Y. Wang, Q. Li, Lattice Boltzmann Method: Theory and Applications, Science Press, Beijing, 2009. [14] Q. Li, K.H. Luo, X.J. Li, Lattice Boltzmann method for relativistic hydrodynamics: issues on conservation law of particle number and discontinuities, Phys. Rev. D 86 (2012) 085044. [15] Y.L. He, Q. Liu, Q. Li, Three-dimensional finite-difference lattice Boltzmann model and its application to inviscid compressible flows with shock waves, Physica A 392 (2013) 4884–4896. [16] Y.X. Xia, Y.H. Qian, Lattice Boltzmann simulation for forced two-dimensional turbulence, Phys. Rev. E 90 (2) (2014) 023004. [17] S. Succi, Lattice Boltzmann across scales: from turbulence to DNA translocation, Eur. Phys. J. B 64 (2008) 471–479. [18] M.A.A. Spaid, F.R. Phelan Jr., Lattice Boltzmann methods for modeling microscale flow in fibrous porous media, Phys. Fluids 9 (1997) 2468–2474. [19] O. Dardis, J. McCloskey, Lattice Boltzmann scheme with real numbered solid density for the simulation of flow in porous media, Phys. Rev. E 57 (1998) 4834–4837. [20] D.M. Freed, Lattice-Boltzmann method for macroscopic porous media modeling, Internat. J. Modern Phys. C 9 (1998) 1491–1503. [21] N.S. Martys, Improved approximation of the Brinkman equation using a lattice Boltzmann method, Phys. Fluids 6 (2001) 1807–1810. [22] Z. Guo, T.S. Zhao, Lattice Boltzmann model for incompressible flows through porous media, Phys. Rev. E 66 (2002) 036304. [23] Q. Kang, D. Zhang, S. Chen, Unified lattice Boltzmann method for flow in multiscale porous media, Phys. Rev. E 66 (2002) 056307. [24] S. Anwar, M.C. Sukop, Regional scale transient groundwater flow modeling using lattice Boltzmann methods, Comput. Math. Appl. 58 (2009) 1015–1023. [25] S.D.C. Walsh, H. Burwinkle, M.O. Saar, A new partial-bounceback lattice-Boltzmann method for fluid flow through heterogeneous media, Comput. Geosci. 35 (2009) 1186–1193. [26] P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases. I. small amplitude processes in charged and neutral one-component systems, Phys. Rev. 94 (1954) 511–525. [27] P. Lallemand, L.-S. Luo, Theory of the lattice Boltzmann method: dispersion, dissipation, isotropy, Galilean invariance, and stability, Phys. Rev. E 61 (2000) 6546. [28] D. d’Humières, I. Ginzburg, M. Krafczyk, P. Lallemand, L.-S. Luo, Multiple-relaxation-time lattice Boltzmann models in three dimensions, Philos. Trans. R. Soc. Lond. Ser. A 360 (2002) 437–451. [29] Q. Li, Y.L. He, G.H. Tang, W.Q. Tao, Improved axisymmetric lattice Boltzmann scheme, Phys. Rev. E 81 (2010) 056707. [30] I. Ginzburg, D. d’Humières, Multireflection boundary conditions for lattice Boltzmann models, Phys. Rev. E 68 (2003) 066614. [31] D. d’Humières, Generalized lattice-Boltzmann equations, in: B.D. Shizgal, D.P. Weaver (Eds.), Rarefied Gas Dynamics: Theory and Simulations, in: Prog. Astronaut. Aeronaut., vol. 159, AIAA, Washington, DC, 1992, pp. 450–458. [32] R. Du, B. Shi, Incompressible MRT lattice Boltzmann model with eight velocities in 2D space, Internat. J. Modern Phys. C 20 (2009) 1023–1037. [33] S. Ergun, Fluid flow through packed columns, Chem. Eng. Prog. 48 (1952) 89–94. [34] K. Vafai, Convective flow and heat transfer in variable-porosity media, J. Fluid Mech. 147 (1984) 233–259. [35] M.E. MaCracken, J. Abraham, Multiple-relaxation-time lattice-Boltzmann model for multiphase flow, Phys. Rev. E 71 (2005) 036701. [36] Q. Li, K.H. Luo, X.J. Li, Lattice Boltzmann modeling of multiphase flows at large density ratio with an improved pseudopotential model, Phys. Rev. E 87 (2013) 053301. [37] Y.H. Qian, D. d’Humeier, P. Lallemand, Lattice BGK models for Navier–Stokes equation, Europhys. Lett. 17 (1992) 479–484. [38] Z. Guo, B. Shi, N. Wang, Lattice BGK model for incompressible Navier–Stokes equation, J. Comput. Phys. 165 (2000) 288–306. [39] Z. Chai, T.S. Zhao, Effect of the forcing term in the multiple-relaxation-time lattice Boltzmann equation on the shear stress or the strain rate tensor, Phys. Rev. E 86 (2012) 016705. [40] Z.L. Guo, C.G. Zheng, B.C. Shi, Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method, Chin. Phys. 11 (2002) 366. [41] N. Martys, D.P. Bentz, E.J. Garboczi, Computer simulation study of the effective viscosity in Brinkman’s equation, Phys. Fluids 6 (1994) 1434–1439. [42] U. Ghia, K.N. Ghia, C.T. Shin, High-Re solutions for incompressible flow using the Navier–Stokes equations and a multigrid method, J. Comput. Phys. 48 (1982) 387–411. [43] Z. Guo, T.S. Zhao, A lattice Boltzmann model for convection heat transfer in porous media, Numer. Heat Transfer B 47 (2005) 157–177. [44] Q. Liu, Y.L. He, Q. Li, W.Q. Tao, A multiple-relaxation-time lattice Boltzmann model for convection heat transfer in porous media, Int. J. Heat Mass Transfer 73 (2014) 761–775. [45] Q. Li, Y.L. He, Y. Wang, G.H. Tang, An improved thermal lattice Boltzmann model for flows without viscous heat dissipation and compression work, Internat. J. Modern Phys. C 19 (2008) 125–150.