Lattice model effects on the accuracy of the boundary condition implementations for the convection–diffusion lattice Boltzmann method

Accepted Manuscript

Lattice Model Effects on the Accuracy of the Boundary Condition Implementations for the Convection-Diffusion Lattice Boltzmann Method Liangqi Zhang , Shiliang Yang , Zhong Zeng , Jia Wei Chew PII: DOI: Reference:

S0045-7930(18)30577-2 https://doi.org/10.1016/j.compfluid.2018.08.029 CAF 3994

To appear in:

Computers and Fluids

Received date: Accepted date:

19 June 2018 30 August 2018

Please cite this article as: Liangqi Zhang , Shiliang Yang , Zhong Zeng , Jia Wei Chew , Lattice Model Effects on the Accuracy of the Boundary Condition Implementations for the Convection-Diffusion Lattice Boltzmann Method, Computers and Fluids (2018), doi: https://doi.org/10.1016/j.compfluid.2018.08.029

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights

The second-order boundary scheme is extended to the D2Q9 lattice model Consistent implementations of Dirichlet, Neumann, Robin conditions are derived The lattice model effects on the accuracy of the boundary scheme are investigated Second-order accuracy of macroscopic scalar distribution preserved by both models D2Q9 model more accurate than D2Q5 model for curved boundaries

AC

CE

PT

ED

M

AN US

CR IP T

    

1

ACCEPTED MANUSCRIPT

Lattice Model Effects on the Accuracy of the Boundary Condition Implementations for the Convection-Diffusion Lattice Boltzmann Method Liangqi Zhanga, Shiliang Yanga, Zhong Zengb,c, Jia Wei Chewa,d, 1 a

School of Chemical and Biomedical Engineering, Nanyang Technological University, Singapore 637459, Singapore

b

CR IP T

Department of Engineering Mechanics, Colledge of Aerospace Engineering, Chongqing University, Chongqing, 400044, PR China

c

State Key Laboratory of Coal Mine Disaster Dynamics and Control, Chongqing University, Chongqing, 400044, PR China

d

Singapore Membrane Technology Center, Nanyang Environment and Water Research Institute, Nanyang Technological University, Singapore 637141, Singapore

AN US

Abstract

An alternative second-order boundary scheme for convection-diffusion lattice Boltzmann method was proposed in our recent work with advantages of locality and consistency for various boundary conditions. In this study, the boundary scheme, previously derived for the D2Q5 lattice model, is firstly extended to the D2Q9 model with more lattice directions, and then the accuracy

M

of the boundary condition implementations based on different lattice models is compared based on various boundary conditions, namely, Dirichlet condition, Neumann condition and linear

ED

Robin condition. Conclusions are drawn from the numerical validations that: (i) for the Neumann condition and Robin condition at straight boundaries, the D2Q5 model is proved to be more

PT

accurate than the D2Q9 model, and the inaccurate treatments for the corner points by the D2Q9 model would further affect its accuracy for the scalar gradients distribution; (ii) for curved

CE

boundaries, the D2Q9 model is demonstrated to be more accurate than the D2Q5 model that has less lattice directions; and (iii) the second-order accuracy of the resulting macroscopic scalar

AC

distribution is preserved by both lattice models for all the numerical tests. Furthermore, a modified scheme for the Dirichlet condition treatments at straight boundaries is proposed for the D2Q5 lattice model, and leads to obvious accuracy improvements. Keywords: lattice Boltzmann method; convection-diffusion equation; Neumann condition; Robin condition

1

Author for correspondence. Tel: +65 6316 8916 (J.W. Chew); E-mail: [email protected] (J.W. Chew)

2

ACCEPTED MANUSCRIPT

1. Introduction The lattice Boltzmann (LB) method, derived from the kinetic theory, has developed into a convenient alternative for the traditional computational fluid dynamics (CFD) method, and

CR IP T

achieved wide applications in various hydrodynamic systems as well as the coupled heat and mass transfer processes [1-10]. In addition to the advances with respect to the LB model

developments for recovering accurate macroscopic equations (namely, the Navier-Stokes (N-S) equations and the coupled convection-diffusion equation (CDE)) [11-20], accurate boundary condition implementations are also essential to the overall integrity of the LB simulations and

AN US

boundary schemes with second-order accuracy are highly desired [21-23].

It should be noted that a large portion of the existing boundary schemes for the CDE LB method is derived as a direct extension from their hydrodynamic schemes counterpart. For example, Dirichlet condition treatments for the CDE LB method was proposed by Zhang et al. [24, 25] based on the inherent bounce-back rule for velocity boundary conditions, and then modified by

M

Chen et al. [22] through a midpoint bounce-back scheme [26]. When extending these two bounce-back rule based schemes to Neumann and Robin conditions, first-order finite difference

ED

approximation of the spatial derivative was necessary to reduce the Neumann and Robin conditions to the Dirichlet condition [22, 24]. Alternatively, another attempt for the extended application of the bounce-back rule to the CDE LB method was carried out by Yoshida et al.

PT

[19]. Subsequently, to improve the accuracy of the bounce-back rule based scheme by Yoshida et al. [19], the interpolation of the distribution functions from the neighboring fluid nodes was

CE

introduced by Li et al. [21] and the coefficients for the involved distribution functions were theoretically determined to satisfy the Dirichlet and the Neumann conditions with second-order

AC

accuracy. Numerical validations confirmed the second-order accuracy of the Li et al. scheme [21] for straight boundaries, but only first-order accuracy was demonstrated for the Neumann condition at curved boundaries. Moreover, a second-order bounce-back scheme for the CDE LB method was recently proposed by Huang et al. [27] based on their previous single-node bounceback scheme [28], and second-order Taylor expansion and complicated interpolation were applied to preserve the overall accuracy. Despite the complicated practical implementations, the

3

ACCEPTED MANUSCRIPT

Huang et al. scheme [27] provides second-order treatments for the Neumann condition with nonzero flux and the Robin conditions at curved boundaries for the first time. Besides, some other boundary schemes are directly derived for the CDE LB method without a hydrodynamic counterpart, such as the local scheme based on the moments of the distribution functions for the Robin condition [23], the mass conservation scheme by Gebäck et al. [29], and

CR IP T

the direct extrapolation scheme for the Neumann condition with zero flux by Huang et al. [30]. In our recent work [31], a second-order boundary scheme was proposed based on a second-order expansion of the distribution functions at the boundary nodes, and then consistent

implementations of the Dirichlet, Neumann and Robin conditions were developed for both

straight and curved geometries. Other than the verified second-order accuracy, the proposed

AN US

boundary scheme has the advantage of its locality and consistency, namely, no information from the neighboring fluid nodes is required and all the three kinds of boundary conditions are directly implemented without any additional transformation.

Although the above boundary schemes as well as the LB models for CDE are mainly developed

M

for the D2Q5 lattice models, the D2Q9 model with more lattice directions is reported to be helpful in recovering accurate macroscopic equations especially for some complicated convection-diffusion systems [13, 15, 32, 33], and recent numerical validations without

ED

boundary scheme effects demonstrated that the D2Q9 model was more accurate for systems with higher convection strength [17]. Therefore, the extension of the boundary schemes to the D2Q9

PT

lattice model without affecting its second-order accuracy is helpful to expand the applications of the CDE LB method further. However, it has been reported that a complicated Cartesian

CE

Decomposition Method (CDM) is necessary for extending the Li et al. scheme to the D2Q9 model for Neumann condition treatments, and the resulting derived scheme with more lattice

AC

directions leads to reduced accuracy [17]. Moreover, the second-order bounce-back scheme by Huang et al. [27] as well as its predecessor (i.e. the single-node scheme [28]) are both built for D2Q5 lattice models and their extension to the D2Q9 are still not determined. In this work, our newly proposed boundary schemes for CDE LB method is firstly extended to the D2Q9 lattice model, and then the lattice model effects on the accuracy of the boundary implementations are investigated with respect to the theoretical derivations and numerical validations. It is well demonstrated that, in addition to its locality and consistency, our proposed boundary scheme also 4

ACCEPTED MANUSCRIPT

has the advantage of flexibility for various lattice models. Moreover, a modified treatment for the Dirichlet condition at straight boundaries is proposed for the D2Q5 model and the resulting accuracy improvements are confirmed by numerical tests. Lastly, we make our recommendations for practical CDE LB simulations based on the present theoretical and numerical investigations of the lattice model effects on the accuracy of the boundary implementations. Even though the

CR IP T

present discussions are conducted for two-dimensional (2D) systems, the implementation on 3D cases can be made straightforwardly.

2. Formulations

AN US

In this part, the newly proposed boundary scheme is firstly extended to the D2Q9 lattice model, and the original boundary scheme for the D2Q5 model is also provided as its reduced version. Thus, the comparisons for the two sets of developed boundary condition implementations are performed based on the corresponding theoretical derivations. Furthermore, a modified scheme for the Dirichlet condition treatments at straight boundaries is introduced for the D2Q5 lattice

M

model.

ED

2.1 CDE LB model

The LB equation with Bhatnagar-Gross-Krook (BGK) collision operator [34] is considered

PT

f  x  ξ  t , t   t   f  x, t    eq 

 x, t 



 f eq 



(1)

are respectively the distribution function and its equilibrium form

CE

where f  x, t  and f

f  1

at the position x and time t for the macroscopic scalar  (temperature T, concentration C, etc.).

AC

ξ is the discrete particle velocity vector defined by the adopted lattice model

 0 (0, 0) ξ   , for the D2Q5 model c(cos[(  1) / 2],sin[(  1) / 2])   1, 2,3, 4

(2a)

(0, 0)  0  ξ  c(cos[(  1) / 2],sin[(  1) / 2])   1, 2,3, 4 , for the D2Q9 model (2b)   2c(cos[(2  1) / 4],sin[(2  1) / 4])   5, 6, 7,8 5

ACCEPTED MANUSCRIPT

where c   x  t is the lattice velocity, and  x and  t are respectively the discrete spatial and time steps, which are both prescribed to be dimensionless constant 1 in the following LB simulations. Besides, excluding the effects from the velocity square terms, the following equilibrium distribution function expression is applied for both lattice models eq 

 u  ξ    w 1   RT   

(3)

CR IP T

f

where u is the fluid velocity, R and T  are respectively the ideal gas constant and reference temperature, with RT   1 3 , and w is the weight coefficient given by the lattice models

4 9,  w  1 9, 1 36, 

 0 , for the D2Q5 model   1, 2,3, 4

AN US

1 3, w   1 6,

 0   1, 2,3, 4 , for the D2Q9 model   5, 6, 7,8

(4a)

(4b)

M

Additionally, the single relaxation parameter in Eq. (1),  , is related to the diffusion coefficient

ED

D in the recovered macroscopic equations as

PT

 1 D     RT  t  2

(5)

CE

Moreover, the macroscopic scalar  is calculated as

   f

(6)



AC

and the following computing formula is applied for the scalar gradients calculations in the present numerical validations

3  t 1     2 u   f ξ   x x  

2.2 Expression of the distribution function f

6

(7)

ACCEPTED MANUSCRIPT

The Chapman-Enskog (CE) analysis is firstly implemented for the above CDE LB model with the following multi-scale expansions     2 t t0 t1

(8a)

CR IP T

   x x0 f  f(0)   f(1)   2 f(2)

(8b)

(8c)

where  is the expansion parameter. Applying Eq. (8) to the second-order Taylor expansion of

AN US

Eq. (1) 2

   t2     1  eq   t    i f       i  f   f  f  t  x 2  t  x  i  i   

we obtain

    0 : f  f

eq 

M

0



    (0) 1 (1)   i f  f   xi   t   t0

 (0)  1     (1) 1 (2) f  1      i f  f   t1 x0i   t   2   t0

PT

2 :

(9)

(10a) (10b)

ED

1 : 



(10c)

AC

CE

Moments of the Eqs. (10b) and (10c) are determined as

  f(0)   t0  x0 j

  f(0) i   t0  x0 j

f  

 f

 j  0

(0)

 i j  

1

(0)

 1    f(0)  1    t1   2  x0 j

7

f  

 t

(11a)

f  

 j  0

(1)

 i

(1)

(11b)

(11c)

ACCEPTED MANUSCRIPT

Substituting the moments of the equilibrium distribution function in Eq. (3) f  

(0)

 ,

f  

 j   u ,

(0)

f  

i j   RT  ij

(0)

(12)

and neglecting the terms higher than O(Ma 2 ) , with Ma being the Mach number, we obtain

 i   tRT 

(1)

 x0i

CR IP T

 f

(13)

Then, the macroscopic equation recovered from the CDE LB model is derived as    D  x j

  

AN US

    u j    t x j x j

(14)

Introducing the Hermite expansion technique [35, 36], the dimensionless distribution function

RT  fˆ  f (with 0  1 being the constant density) is expressed as

0

M

 1 fˆ (ξˆ, x, t )   (ξˆ ) ai( n ) ( x, t ) H i( n ) ( ξˆ) n 0 n !

(15)

RT 

, and the leading order Hermite polynomials are defined as

CE

PT

velocity ξˆ  ξ

ED

where H i( n ) (ξˆ ) is the tensorial Hermite polynomials with respect to the dimensionless particle

H (0) (ξˆ )  1

(16a)

H i(1) (ξˆ )  ˆi

(16b)

AC

Besides, the expansion coefficient in Eq. (15) is determined as ˆ ( n ) (ξˆ )d ξˆ  1 ai( n ) ( x, t )   fH i

0

 fH

(n) i

(ξˆ )dξ

(17)

which is naturally related to the moments of the distribution function. Therefore, excluding the neq second order term  2 f(2) in Eq. (8c), the non-equilibrium part of the distribution f  is

expanded in accordance with its moments 8

ACCEPTED MANUSCRIPT

f  

neq 

neq 

 0,

 f  f

f  

eq 

  f  1

(18)

 i    f1 i   tRT 

neq 



f neq    w t

 xi

  i xi

(19)

(20)

CR IP T

f

Thus, the second-order expression for the distribution function f , derived from the CE analysis and the Hermite expansion technique, is obtained as

(21)

AN US

   u  f  f eq   f neq   w      i i   t  i  RT  xi   2.3 Present boundary scheme for the D2Q9 lattice model

The distribution functions at the computational boundary nodes are firstly divided into two subsets, namely, I in and I loc respectively for the unknown distribution functions fin and the

M

remaining known ones floc after the streaming step. Particularly, N in and N loc denote the corresponding two subsets for the indices of the distributions functions  . The essence of the

ED

boundary implementation in the LB method is the determination of the unknown distribution functions fin at the boundary nodes in accordance with the macroscopic physical constraints at

PT

the boundary surface. However, it is challenging to build the closure relations at the boundary nodes since the number of the unknown distribution functions fin is always larger than the

CE

number of the physical constraints for the macroscopic fluid variables, especially for the D2Q9 lattice model with more lattice directions. Therefore, the Robin condition treatments proposed by

AC

Meng et al. [23] are hard to be applied to the CDE LB model with D2Q9 lattice structure, and the extension of the Li et al. scheme from the D2Q5 model [21] to the D2Q9 model requires more complicated interpolations [17]. However, in the present boundary scheme, in addition to the macroscopic physical boundary constraints, the local known distribution functions floc are introduced by applying the derived expression in Eq. (21) to complete the closure at the same boundary nodes, and thus the unknown distribution functions fin are determined locally.

9

ACCEPTED MANUSCRIPT

2.3.1 Straight boundaries The boundary nodes are directly placed on the straight boundary surface, as demonstrated in Fig.

CR IP T

1a, and the division for the local distribution functions at boundary node P is carried out as

AN US

Fig. 1 Schematic description of a boundary node at a straight boundary (a) and around a curved boundary (b) for the D2Q9 lattice model, with the red dashed lines denoting the directions of the unknown distribution functions.

M

I in   f 2 , f5 , f6  , I loc   f0 , f1 , f3 , f 4 , f7 , f8

ED

N in  2,5,6 , N loc  0,1,3, 4, 7,8

(22a) (22b)

S ,

 S  S , n  

(23)

 S  S and are the spatial derivatives of the scalar  in the normal (n) and tangential n  

CE

where

PT

At the boundary surface, the physical constraints are provided for the macroscopic variables

(   ) directions, respectively. Besides, the distribution functions at the boundary nodes are

AC

related to the macroscopic variables in Eq. (23) by   u f  x, t   w 1   i i RT  

  S   S  S   w  t        n      n 

 S  S  A  S  B  C n  

with 10

(24)

ACCEPTED MANUSCRIPT

  u A  w 1   i i RT  

  , B   w t n , C  w t  

(25)

where  n and   are the discrete particle velocity components in respectively the normal and tangential directions. Hereinafter, the unknown distribution functions fin are derived from the

CR IP T

known ones floc in accordance with various boundary conditions. Practically, the unknown macroscopic variables at the boundary surface are firstly determined by the known distribution functions floc , then the unknown distribution functions fin are readily obtained through Eq. (24).

AN US

(i) Dirichlet condition

S With the known macroscopic scalar  at the boundary surface given by the Dirichlet condition,

 S the tangential derivative can also be directly derived and here regarded as an ‘implicit  

constraint’ for the Dirichlet condition. Therefore, the normal derivative is the only unknown

M

macroscopic variable in Eq. (24), and is determined by a selected known distribution function

ED

from I loc

PT

 S  loc  S  loc   f  A  S  C  B , with   N and B  0    n 

(26)

(ii) Neumann condition

CE

The normal gradient of the macroscopic scalar

 S is given by the Neumann condition, and the n

AC

remaining macroscopic variables are obtained as

  S  loc  S   floc  B  A with   N and ξ // n  n  

(27a)

 S  loc  S  S   f  A   B C with   N loc and C  0  n    

(27b)

(iii) Linear Robin condition 11

ACCEPTED MANUSCRIPT

With the following definition of the Robin condition  S 1   2  3 n S

(28)

all the unknown macroscopic variables at the boundary surface are determined directly by

CR IP T

substituting the Robin condition into the expression for floc in Eq. (24)

 S   floc  B 3 2   A  B 1 2  , with   N loc and ξ // n  S   3  1 S   2 n

 S is then derived by Eq. (27b).  

(30)

AN US

and the tangential derivative

(29)

2.3.2 Curved boundaries

The curved geometries cause great difficulties in the boundary implementations for the CDE LB

M

method especially for the boundary conditions with respect to the scalar gradients, namely, the Neumann and Robin conditions. Although the second-order bounce-back rule by Huang et al. is

ED

claimed to be the first second-order curved boundary scheme applicable for the Dirichlet, Neumann and Robin conditions, its practical implementations introduce complicated interpolations and second-order derivative approximations. Moreover, the extension of the

PT

Huang et al. scheme to the D2Q9 lattice model is still not available currently. To address these issues, in the present scheme, the curved geometries are firstly represented by a local curvilinear

CE

coordinate system, and the macroscopic information at the physical boundary surface is transferred to the computational boundary nodes by appropriate Taylor expansion. Thus, the

AC

boundary implementations for straight boundaries are extended to curved geometries with second-order accuracy, and the advantages of its locality and consistency are preserved. The schematic for the boundary nodes near the curved boundary surface is given in Fig. 1b, and the distribution functions at the boundary node P are divided as

I in   f3 , f 4 , f7 , f8  , I loc   f0 , f1 , f2 , f5 , f6

12

(31a)

ACCEPTED MANUSCRIPT

N in  3, 4,7,8 , N loc  0,1, 2,5, 6

(31b)

Besides, the normal projection point P S of the boundary node P at the boundary surface is introduced with  n , the normal distance of P to the boundary surface. Furthermore, the local polar coordinate system is applied, with r and  respectively the local radius and azimuth, for an

coordinate system is defined as

gi 

r xi

AN US

g1  cos ex  sin  e y

CR IP T

accurate representation of the curved geometries. The covariant local basis gi for the curvilinear

g2  r sin  ex  r cos e y

(32) (33a) (33b)

where r  r cos  ex  r sin  e y is the arbitrary position vector, and e x and e y are unit basis vectors

M

i in the Cartesian coordinate system. Then, the contravariant basis g is derived from the

orthogonality rule g i  g j   ij , with  ij being the Kronecker delta

PT

ED

g1  cos  ex  sin  e y

g2 

 sin  cos  ex  ey r r

(34a)

(34b)

CE

Therefore, the macroscopic scalar gradients  as well as the distribution function expression

AC

f in Eq. (21) are rewritten in their tensorial form as  

 i g xi

(35a)

      ξ   i g i    j g j   i  i x  x 

(35b)

   ui  f  w      i   t i  i  RT  x  

(35c)

13

ACCEPTED MANUSCRIPT

with  i  ξ  g i . Furthermore, the distribution functions at the boundary node P are related to the macroscopic variables at the boundary surface P S by applying the following Taylor expansions  S  P   n  O  x 2  r

  P     P S   O  x 

CR IP T

  P     PS  

(36a) (36b)

where the normal distance of the boundary node P to the boundary surface  n is directly determined as  n  r  P   r  P s  . Applying the approximation in Eq. (36), we obtain   S i     t    i x   

AN US

  S    i u i  f  P   w   S   n  1   r RT       S  S S     A   B  C r 

M

with

ED

  ui  A  w 1   i  RT   

C   tw  

(38a)

(38b)

(38c)

CE

PT

  u i   B  w 1   i   n   t r  RT    

(37)

It should be noted that, in Eq. (36b), only the first-order approximation of the scalar gradients

AC

 is applied, but the resulting error terms in the distribution function expression in Eq. (37) are expected to be second-order, since  is directly related to the non-equilibrium part of the distribution function that itself is of the order O  t  . In the following, determinations of the in unknown distribution functions f are carried out using Eq. (37) with the prepared macroscopic

information at the boundary surface. (i) Dirichlet condition 14

ACCEPTED MANUSCRIPT

With the direct macroscopic constraint for  S and the ‘implicit constraint’ for

 S , the 

unknown distribution functions can be obtained by Eq. (37) once the remaining unknown macroscopic variable

 S is determined r

CR IP T

 S  loc  S  loc   f  A  S  C  B , with   N and B  0   r  (ii) Neumann condition

 S given by the Neumann condition, the macroscopic r

AN US

With the known information

(39)

information at the boundary surface is computed as

f 0loc  S  n w0 r

(40)

 S  loc  S  loc   f  A  S  B  C , with   N and C  0 r   

(41)

M

S 

ED

(iii) Linear Robin condition

 S  3 (i.e., the definition of the Robin condition), the r

PT

With the relation 1 S   2

AC

CE

macroscopic information at the boundary surface are derived as

 f loc         S   0  n 3  1  n 1  2   2   w0  S   3  1 S   2 r

 S In addition, Eq. (41) is then employed for the tangential derivative . 

2.4 Present boundary scheme for the D2Q5 lattice model

15

(42)

(43)

CR IP T

ACCEPTED MANUSCRIPT

Fig. 2 Schematic description of a boundary node at a straight boundary (a) and around a curved boundary (b) for the D2Q5 lattice model, with the red dashed lines denoting the directions of the unknown distribution functions.

AN US

The D2Q5 model defined in Eq. (4a) is here regarded as a simplified lattice model from the D2Q9 model by neglecting the lattice directions in the diagonal directions, and thus the computational efficiency is highly improved in terms of both the basic evolutions and the boundary implementations. As schematically demonstrated in Fig. 2 for a computational boundary node at straight (Fig. 2a) and curved (Fig. 2b) boundaries, the number of the unknown

M

distribution functions is obviously reduced for both kinds of geometries by the D2Q5 model. For

ED

the straight boundary, the only one unknown distribution function is along the inner normal direction and thus the tangential derivative

 S is unnecessary in the practical boundary  

PT

implementations, which simplifies the boundary implementations. On the other hand, for the curved boundary, the above implementations for the D2Q9 lattice model are consistently applied

CE

to the D2Q5 lattice model, and thus detailed discussions are not included here. Therefore, boundary implementations based on the D2Q5 model are provided only for the straight boundary as follows. Firstly, the division of the distribution functions at the boundary node P in Fig. 2a is

AC

introduced as

I in   f 2  , I loc   f0 , f1 , f3 , f4

(44a)

N in  2 , N loc  0,1,3, 4

(44b)

16

ACCEPTED MANUSCRIPT

Since only the distribution function along the inner normal direction is unknown, the boundary conditions are efficiently implemented by Eq. (24) after the determination of  S and

 S . n

With the given  S , the normal derivative

 S is determined as n

CR IP T

(i) Dirichlet condition

 S   floc  A  S  B , with   N loc and ξ // n n

AN US

(ii) Modified scheme for Dirchlet condition

As discussed above, an ‘implicit constraint’ on the tangential derivative

(45)

 S is derived from the  

Dirichlet condition, and can be naturally introduced to the D2Q9 model based implementations but ignored by the D2Q5 model. Therefore, a modified scheme is introduced by including the

M

‘implicit constraint’, namely, the distribution functions along the tangential directions (i.e., f1 , f3 for the present boundary node P in Fig. 2a) are also updated by Eq. (24) with the macroscopic

(iii) Neumann condition

 S S , the macroscopic scalar  is obtained as n

PT

With the given

ED

information at the boundary nodes.

CE

  S  loc  S   floc  B  A , with   N and ξ // n  n  

(46)

AC

(iv) Linear Robin condition Applying the relation 1 S   2

 S  3 , the necessary macroscopic information is derived by n

Eqs. (29) and (30). It is demonstrated from above theoretical derivations that the present boundary scheme, originally proposed for the D2Q5 lattice model, can be readily extended to the D2Q9 model 17

ACCEPTED MANUSCRIPT

without affecting the theoretical accuracy, and its advantages over the other existing secondorder boundary scheme for the CDE LB method in terms of its locality and consistency are also preserved. Moreover, the original boundary treatments within the framework of the D2Q5 model are here regarded as the reduced implementations of the boundary scheme based on the D2Q9 model. Thus, the resulting D2Q5 lattice model has the advantage of reducing the computational

CR IP T

cost due to lesser unknown distribution functions at the boundary nodes. Particularly, for the straight boundary treatments, only the distribution function along the inner normal direction is unknown, thus the scalar gradient at the tangential direction becomes unnecessary in the practical boundary implementations and the computational efficiency is further improved. However, the ‘implicit constraint’ on the tangential derivative, derived from the Dirichlet condition at straight

AN US

boundaries, is naturally included by the D2Q9 model based boundary scheme but ignored by the implementations with the D2Q5 model. Therefore, a modified Dirichlet condition treatment based on the D2Q5 model for straight boundaries is provided by updating the distribution functions along the tangential direction, and the modified scheme is expected to be helpful for the accuracy improvements. Therefore, the D2Q9 lattice model with the additional diagonal

M

lattice directions would inevitably lead to higher computational cost than the D2Q5 model, but the ‘implicit constraint’ for the Dirichlet condition at a straight boundary is naturally included by

ED

the D2Q9 model based implementations, and more lattice directions are also expected to be helpful for curved boundary treatments when representing the curved geometries. The present theoretical investigations on the lattice model effects are further validated by the following

PT

numerical tests, and the corresponding analytical solutions are introduced for comparing the

CE

order of the numerical accuracy.

AC

3. Numerical results

In this part, four numerical tests with available analytical solutions for various macroscopic constraints and boundary geometries are firstly adopted here for numerical investigations of the lattice model effects on the accuracy of the boundary implementations for the CDE LB method. Besides, a more practical test is also included to validate the applicability of the boundary schemes based on two lattice models. The unified LB model with BGK collision operator and the equilibrium distribution function excluding the velocity square terms (i.e. Eq. (3)) is applied for 18

ACCEPTED MANUSCRIPT

all the simulations to remove the effects from the CDE LB models. Moreover, the relative errors for the macroscopic scalar  and its gradients  with the following definitions are employed to evaluate the numerical accuracy of the boundary implementations.



E2   

numerical

  analytical

2

i, j

CR IP T



(47a) analytical 2

i, j

E2I    

   

numerical

   

analytical 2

interior nodes

   

analytical 2

AN US

interior nodes

(47b)

where macroscopic variables with the superscript ‘numerical’ and ‘analytical’ denote the numerical results and the corresponding exact solution, respectively. Besides, the relative error for the scalar  is determined based on all the lattice nodes in the entire computational domain,

M

while E2I    is defined for the interior fluid nodes by excluding the boundary nodes. Moreover, for the accuracy evaluation of the curved boundary treatments, the additional relative

ED

errors are defined based on the computational boundary nodes for the boundary value and

PT

boundary flux

AC

CE

E2B   

where n 

E2B n  



 numerical   analytical

2

boundary nodes



 analytical

(48a)

2

boundary nodes



n numerical  n analytical

2

boundary nodes



n analytical

2

(48b)

boundary nodes

 is the normal derivative. Particularly, we note that, due to the non-uniqueness of n

the solutions for the Neumann condition tests, additional treatments are introduced for the

19

ACCEPTED MANUSCRIPT

boundary nodes to restrict the numerical results to the analytical solutions in the practical simulations. Therefore, the calculation of E2   for the Neumann condition tests is also based on the interior fluid nodes. 3.1 One-dimensional (1D) transient convection-diffusion D2Q5, t=1 D2Q9, t=1

10

Analytical Solution t=5

4

Analytical Solution t=1 D2Q5, t=2 D2Q9, t=2

40

C(t,x)

Analytical Solution t=2 D2Q5, t=5 D2Q9, t=5

6

D2Q5, t=1 D2Q9, t=1

50

Analytical Solution t=1 D2Q5, t=2 D2Q9, t=2

8

C(t,x)

(b)

CR IP T

(a) 12

Analytical Solution t=2 D2Q5, t=5 D2Q9, t=5

30

Analytical Solution t=5

Pe=10

20

Pe=1 2 0 0.0

0.2

0.4

0.6

0.8

AN US

10

0 0.0

1.0

0.2

0.4

0.6

0.8

1.0

x/L

x/L

Fig. 3 Concentration distributions at various times obtained from the boundary implementations based on both

M

lattice models for cases of: (a) Pe=1; and (b) Pe=10.

ED 10-2

AC

CE

PT

Relative Error for Concentration

10-1

10-3

10-4

10-5

D2Q5, Pe=1 D2Q9, Pe=1 D2Q5, Pe=10 D2Q9, Pe=10 Slope=-2 20

40

60

80

100 120

L

Fig. 4 Relative error for concentration versus mesh resolution L   N x  1  x obtained from the boundary implementations based on both lattice models at Pe values of 1 and 10. The reference dashed pink line has a slope of -2.

20

ACCEPTED MANUSCRIPT

The present 1D unsteady-state test is defined as C C  2C u  D 2 , 0  x  L, t  0 t x x

(49a)

C  x, t  t 0  0,

(49b)

C    uC  D   uC f , x  x 0 

C x

 0,

CR IP T

0 xL

t 0

t 0

xL

(49c)

(49d)

AN US

with the diffusion coefficient D=0.01 and C f  50 for all the test cases. By the dimensionless Peclet number Pe  uL D , the convective velocity u is related to L   N x  1  x for each test case, and N x denotes the number of lattice nodes in the x direction. In addition to the Robin condition (left) and the Neumann condition (right) treatments based on the present boundary

M

scheme, the distribution functions at the bottom and top boundaries are subjected to the following conditions for all the lattice directions

(50a)

f  i, N y   f  i, N y  1

(50b)

PT

ED

f  i,1  f  i, 2 

Thus, the present 1D test can be represented by the adopted 2D CDE LB model. With mesh

CE

resolution N x  41 , the concentration profiles along the x direction at various times are obtained by the boundary implementations based on both lattice models, and demonstrated in Fig. 3 in

AC

comparison with the analytical solution

1  x  ut C  x, t   C f  erfc  2  2 Dt 

  x  ut 2  u 2t   exp      D 4 Dt   

1  ux u t   x  ut  ux   1    exp   erfc  2 D D  D  2 Dt 2

21

  

for 0  L  1, t  5

(51)

ACCEPTED MANUSCRIPT

which is valid for the prescribed D and C f . For both test cases of Pe=1 (Fig. 3a) and Pe=10 (Fig. 3b), the numerical results from the boundary implementations based on the two lattice models agree well with each other and also with the analytical solutions. Moreover, the numerical accuracy of the boundary scheme, demonstrated by the convergence order of the relative error (Eq. (47a)) versus various mesh resolutions, is compared in Fig. 4 for the two lattice models

CR IP T

considered here. Obviously, the convergence trends are parallel to the reference line with a slope of -2, thus confirming the second-order numerical accuracy of the boundary implementations based on both lattice models. Besides, the relative error versus mesh resolution slopes for both lattice models are similar to each other, which again demonstrates that the present boundary

3.2 Heat transfer in a 2D channel flow

AN US

schemes based on both lattice models are equivalently accurate for the present 1D test.

To compare the accuracy of the boundary implementations for the 2D Dirichlet and Neumann conditions, the 2D heat transfer test is considered here with the following governing equation

  2T  2T  T u  D 2  2  x y   x

M

(52)

where the constant velocity u relates to the diffusion coefficient D by the Peclet number

ED

Pe  uH D , with H   N y  1  x being the height of the channel and Pe =20 is adopted here.

PT

With the periodic condition applied in the x direction, two sets of boundary conditions, i.e., the Dirichlet and Neumann conditions, are prescribed for the bottom and the top boundaries. The

AC

CE

Dirichlet condition and the resulting analytical solution are defined as

with   k 1 

T  x,0   T  x, H   cos  kx  , with k  2 L

(53a)

  exp   y   exp    H  y     TD  x, y   Real exp  ikx      exp   H   1   

(53b)

iU . On the other hand, the Neumann condition and its derived exact solution are Dk

given as

22

ACCEPTED MANUSCRIPT

T y

 y 0

T y

 cos  kx  H

(54a)

yH

  exp   y   exp    H  y     TN  x, y   Real exp  ikx       H  exp   H   1    Relative Error for Interior Gradients

10-1

10-2

10-3

10-4

D2Q5, =0.75 D2Q9, =0.75 D2Q5, =1.5 D2Q9, =1.5 D2Q5, =3.0 D2Q9, =3.0 Slope=-2

20

40

60

80

H

CR IP T

(b)

Dirichlet Condition

Dirichlet Condition

10-1

10-2

10-3

D2Q5, =0.75 D2Q9, =0.75 D2Q5, =1.5 D2Q9, =1.5 D2Q5, =3.0 D2Q9, =3.0 Slope=-2

AN US

Relative Error for Temperature

(a)

10-4

20

100 120

(54b)

40

60

80

100 120

H

Fig. 5 Relative errors at various relaxation parameters (  ) for temperature (a) and the interior temperature

M

gradient (b) versus mesh resolution H   N y  1  x for heat transfer in a 2D channel flow with the Dirichlet

pink line has a slope of -2.

(b)

PT

10-1

Relative Error for Interior Gradients

Dirichlet Condition

10-2

D2Q5, =0.75 Modified Scheme, =0.75 D2Q5, =1.5 Modified Scheme, =1.5 D2Q5, =3.0 Modified Scheme, =3.0 Slope=-2

CE

10-3

10-4

AC

Relative Error for Temperature

(a)

ED

condition, obtained from the boundary implementations based on both lattice models. The reference dashed

20

40

60

80

Dirichlet Condition

10-1

10-2

10-3

10-4

D2Q5, =0.75 Modified Scheme, =0.75 D2Q5, =1.5 Modified Scheme, =1.5 D2Q5, =3.0 Modified Scheme, =3.0 Slope=-2

20

100 120

40

60

80

100 120

H

H

Fig. 6 Relative errors at various relaxation parameters (  ) for temperature (a) and the interior temperature gradient (b) versus mesh resolution H   N y  1  x for heat transfer in a 2D channel flow with the Dirichlet condition, obtained from the original boundary implementations based on the D2Q5 lattice model and its modified implementations. The reference dashed pink line has a slope of -2.

23

ACCEPTED MANUSCRIPT

For the Dirichlet condition test, the relative errors for the entire temperature distribution (Eq. (47a)) and the interior temperature gradients (Eq. (47b)) from the two lattice models are obtained with various relaxation parameters, and their convergence trends with respect to the mesh resolution are demonstrated in Fig. 5. Firstly, the second-order accuracy of the resulting temperature field as well as the interior gradients are confirmed for both lattice models. The

CR IP T

relative difference between the two models is greater for temperature (Fig. 5a) than for the temperature gradients (Fig. 5b), and the D2Q9 lattice model is more accurate than the D2Q5 model in both cases. As discussed in Sec. 2.4, the ‘implicit constraint’ derived from the Dirichlet condition cannot be included by the boundary implementations based on the D2Q5 model, and thus a modified scheme is developed by updating the distribution functions at the tangential

AN US

directions. As demonstrated in Figs. 6a and 6b, the relative errors for the temperature and the interior temperature gradients are both effectively reduced by the modified scheme, and the

Neumann Condition

10-2

10-4

D2Q5, =0.75 D2Q9, =0.75 D2Q5, =1.5 D2Q9, =1.5 D2Q5, =3.0 D2Q9, =3.0 Slope=-2

20

ED

10-3

(b)

40

60

80

Neumann Condition

10-1

M

10-1

Relative Error for Interior Gradients

(a)

PT

Relative Error for Interior Temperature

accuracy improvements confirm our previous theoretical analysis.

10-2

10-3

10-4

D2Q5, =0.75 D2Q9, =0.75 D2Q5, =1.5 D2Q9, =1.5 D2Q5, =3.0 D2Q9, =3.0 Slope=-2

20

100 120

40

60

80

100 120

H

H

CE

Fig. 7 Relative errors at various relaxation parameters (  ) for the interior temperature (a) and the interior temperature gradient (b) versus mesh resolution H   N y  1  x for heat transfer in a 2D channel flow with

AC

the Neumann condition, obtained from the boundary implementations based on both lattice models. The reference dashed pink line has a slope of -2.

The relative errors for the Neumann condition test are all determined based on the interior fluid nodes, and the obtained results are demonstrated in Fig. 7, which shows that both lattice models can preserve the second-order accuracy of the boundary implementation for the Neumann 24

ACCEPTED MANUSCRIPT

condition at a straight boundary. For the cases of  =0.75 and 1.5, the errors from both lattice models for the interior temperature field (Fig. 7a) agree well with each other, but for a larger relaxation parameter  =3.0, the D2Q5 model is more accurate. As for the interior temperature gradients field in Fig 7b, equivalent accuracy is demonstrated for the  =1.5 case, but the D2Q5 model is again more accurate for other cases. Additionally, convergent solutions cannot be

CR IP T

obtained by the D2Q9 lattice model with the mesh resolution N y =21 and relaxation parameter  =3.0. This is because the large diffusion coefficient D caused by a relatively larger relaxation parameter  and a coarse mesh resolution, i.e., a small H, lead to a large convective velocity u for the fixed Pe, and the simulation diverge due to the resulting increased compressibility errors.

10-2

10-3

10-4

D2Q5, Pe=0 D2Q9, Pe=0 D2Q5, Pe=20 D2Q9, Pe=20 D2Q5, Pe=40 D2Q9, Pe=40 D2Q5, Pe=80 D2Q9, Pe=80 Slope=-2

20

10

-1

40

60

80

10-2

D2Q5, Pe=0 Modified Scheme, Pe=0 D2Q5, Pe=20 Modified Scheme, Pe=20 D2Q5, Pe=40 Modified Scheme, Pe=40 D2Q5, Pe=80 Modified Scheme, Pe=80 Slope=-2

10-3

10-4

20

100 120

40

80

100 120

ED

Neumann Condition =0.75

10-1 10-2

PT

Relative Error for Temperature

CE

60

H

H

(c)

AC

Dirichlet Condition =0.75 & D2Q5

AN US

-1

Relative Error for Temperature

10

(b)

Dirichlet Condition =0.75

M

Relative Error for Temperature

(a)

10-3 10-4 10-5 10-6

D2Q5, Pe=0 D2Q9, Pe=0 D2Q5, Pe=20 D2Q9, Pe=20 D2Q5, Pe=40 D2Q9, Pe=40 D2Q5, Pe=80 D2Q9, Pe=80 Slope=-2

20

40

60

80

100 120

H

Fig. 8 The convection strength effects on the accuracy of the temperature distribution for heat transfer in a 2D channel flow: (a) Dirichlet condition test by boundary implementations based on both lattice modes; (b) Dirichlet condition test by the original D2Q5 based boundary treatment and its modified implementation; (c) Neumann condition test by boundary implementations based on both lattice modes. The reference dashed pink line has a slope of -2.

25

ACCEPTED MANUSCRIPT

Moreover, the effects of the convection strength (i.e., magnitude of Pe) on the accuracy of the present boundary implementations are investigated with the present two tests, and the convergence of the temperature errors for various Pe values are exhibited in Fig. 8. The following are clear: (i) for all the test cases, the temperature errors increase with the increase of the convection strength; (ii) the D2Q9 lattice model is more accurate than the D2Q5 model for

CR IP T

the Dirichlet condition test at various Pe values; (iii) the accuracy improvements by the modified Dirichlet boundary treatments over the original D2Q5 implementations are further confirmed by tests with various Pe; and (iv) the accuracy of the D2Q9 model is similar with that of the D2Q5 model for the Neumann test cases for Pe=20, 40, 80, but the D2Q5 model is more accurate when

3.3 Helmholtz equation in a square domain -1

10-2

10-4 10-5 10-6

D2Q5, =0.75 D2Q9, =0.75 D2Q5, =1.5 D2Q9, =1.5 D2Q5, =3.0 D2Q9, =3.0 Slope=-2

20

M

10-3

40

60

80

Relative Error for Interior Gradients

(b) 10

ED

Relative Error for Concentration

(a) 10

AN US

Pe=0, i.e., convection effects are absent.

100 120

-1

10-2 10-3 10-4 10-5 10-6

D2Q5, =0.75 D2Q9, =0.75 D2Q5, =1.5 D2Q9, =1.5 D2Q5, =3.0 D2Q9, =3.0 Slope=-2

20

L

40

60

80

100 120

PT

L

Fig. 9 Relative errors at various relaxation parameters (  ) for concentration (a) and the interior

CE

concentration gradient (b) versus mesh resolution

L   N x  1  x

for the Helmholtz equation in a

square domain, obtained from the boundary implementations based on both lattice models. The

AC

reference dashed pink line has a slope of -2.

26

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

ED

Fig. 10 The distribution of the relative errors for the Helmholtz equation in a square domain obtained from boundary implementations based on both lattice models: (a) concentration errors for the D2Q5 model; (b) concentration gradients errors for the D2Q5 model; (c) concentration errors for the D2Q9 model; (d)

CE

PT

concentration gradients errors for the D2Q9 model.

The lattice model effects on the Robin condition implementations at the straight boundaries are

AC

investigated by the following test  2C    2   2  C ,

0  x  1, 0  y  1

(55a)

C  x, 0  exp   x 

(55b)

C  x,1  0

(55c)

C  0, y   sinh   1  y   sinh   

(55d)

27

ACCEPTED MANUSCRIPT

C 1, y  

C 1, y  0 x

(55e)

with     1 , and its analytical solution is derived as C  x, y   exp   x  sinh   1  y   sinh   

(56)

CR IP T

The relative errors for the concentration and the interior concentration gradients fields resulting from the present boundary implementations are demonstrated in Fig. 9. The concentration errors from the D2Q5 model converge with a second-order rate with respect to the mesh resolution, while the D2Q9 model make the relative errors converge at a lower rate for the cases of   0.75 and 1.5 as in Fig 9a. Besides, the D2Q5 model is proved to be more accurate than the D2Q9 for

AN US

the present boundary implementation, since the relative error values from the D2Q5 for the both the concentration distribution in Fig. 9a and the interior gradients in Fig. 9b are smaller than those of the D2Q9 model. Also, as demonstrated in Fig. 9b, the interior gradient errors from the D2Q5 model converge faster than those of the D2Q9 model. Moreover, the error distributions for

M

the case of N x  21 and   0.75 are given in Fig. 10, which indicates that: (i) both lattice models lead to high concentration errors at the bottom left corner of the square (Figs. 10a and 10c), but the D2Q9 model results in additional errors at the right Robin condition boundary (Fig.

ED

10c); (ii) comparison of the concentration gradient error distribution in Figs. 10b and 10d demonstrate that the maximum error value for the D2Q9 model is 10 times larger than that of the

PT

D2Q5 model, and significant errors are caused by the D2Q9 model at the boundary lattice nodes near the two bottom corner points. Additionally, a schematic implementation for the bottom left

CE

corner is given in Fig. 11, and it is demonstrated that the distribution function f 5 at the corner node would enter the fluid region at the next step and can be either determined by the bottom

AC

boundary condition or by the condition at the left boundary, whereas such defects are avoided by the D2Q5 lattice model. Therefore, we conclude that the inappropriate treatments at the corner point by the D2Q9 model is the cause underlying the increase of concentration gradient errors as in Figs. 9b and 10d, which would further lead to larger relative error for the concentration distribution in Fig. 9a.

28

CR IP T

ACCEPTED MANUSCRIPT

Fig. 11 Schematic description of bottom left corner point treatments for the D2Q5 lattice model (a) and the

AN US

D2Q9 lattice model (b).

3.4 Steady-state heat conduction inside a circle

The curved boundary implementations based on the two lattice models are compared with the present test with the following governing equation

0  r  r0 , 0    2

M

2T  0,

(57)

ED

where r0   Nr  1  x is the radius of the circular computational domain. Three kinds of boundary conditions are considered here:

CE

PT

(i) Dirichlet condition

T  cos  k 

(58)

T k  cos  k  r r0

(59)

T  k   1   cos  k  r  r0 

(60)

AC

(ii) Neumann condition

(iii) linear Robin condition

T

29

ACCEPTED MANUSCRIPT

with the constant k  4 . And the following exact solution is used for all the three tests with various boundary conditions k

r T  r ,     cos  k   r0 

10-2

10-4

10

(c)

10

ED

10-5

PT

D2Q5, =0.75 D2Q9, =0.75 D2Q5, =1.5 D2Q9, =1.5 D2Q5, =3.0 D2Q9, =3.0 Slope=-1.5

10-4

10-4

10-5

100

r0

(d)

M

-1

10-2

D2Q5, =0.75 D2Q9, =0.75 D2Q5, =1.5 D2Q9, =1.5 D2Q5, =3.0 D2Q9, =3.0 Slope=-2

10-3

10

Dirichlet Condition

10-3

10-2

100

r0

10

-1

AN US

D2Q5, =0.75 D2Q9, =0.75 D2Q5, =1.5 D2Q9, =1.5 D2Q5, =3.0 D2Q9, =3.0 Slope=-2

10-3

10

Dirichlet Condition

CR IP T

Relative Error for Boundary Value

10

-1

10-5

Relative Error for Interior Gradient

(b)

Dirichlet Condition

Relative Error for Boundary Flux

Relative Error for Temperature

(a)

(61)

100

Dirichlet Condition

10-1

10-2

10-3

10-4

D2Q5, =0.75 D2Q9, =0.75 D2Q5, =1.5 D2Q9, =1.5 D2Q5, =3.0 D2Q9, =3.0 Slope=-1.5 Slope=-1

10

r0

100

CE

r0

Fig. 12 Relative error at various relaxation parameters (  ) for temperature (a), boundary value (b), interior temperature gradient (c), and boundary flux (d) versus mesh resolution r0   Nr  1  x for the heat conduction

AC

inside a circle with the Dirichlet condition, obtained from the curved boundary implementations for both lattice models.

Similarly, the lattice model effects on the curved boundary treatments are investigated by comparing the resulting relative errors. Firstly, the convergence of the errors from the two lattice models are compared in Fig. 12 for the Dirichlet condition test. Since only first-order Taylor 30

ACCEPTED MANUSCRIPT

expansion in Eq. (36b) is applied by the present curved boundary scheme, the convergence rate of the temperature gradients errors is slower than the second-order rate as in Fig. 12c for the interior gradients and Fig. 12d for the boundary flux, but the second-order convergence trend is preserved by both models for the errors of the entire temperature distribution in Fig. 12a and the boundary temperature values in Fig. 12b. Secondly, equivalent accuracy of the two lattice

CR IP T

models is demonstrated for the temperature gradients distributions in Figs. 12c and 12d, and the errors for entire temperature distribution obtained from the two models also agree with each other for the relaxation parameters  =1.5 and 3.0 in Fig. 12a. On the other hand, for the errors of the boundary temperature values in Fig. 12b, the D2Q9 model is slightly more accurate than the D2Q5 model, and the entire temperature distribution for the case of  =0.75 from the D2Q9

AN US (b)

Neumann Condition

10-1

D2Q5, =0.75 D2Q9, =0.75 D2Q5, =1.5 D2Q9, =1.5 D2Q5, =3.0 D2Q9, =3.0 Slope=-2

10-3 10-4 10-5 10

M

10-2

Relative Error for Interior Gradient

(a)

ED

Relative Error for Interior Temperature

model is obviously more accurate than the D2Q5 model in Fig. 12a.

100

10

10-2

D2Q5, =0.75 D2Q9, =0.75 D2Q5, =1.5 D2Q9, =1.5 D2Q5, =3.0 D2Q9, =3.0 Slope=-1.5

10-3

10-4

10-5 10

100

r0

r0

(c)

Relative Error for Boundary Flux

AC

CE

PT

Neumann Condition

10-1

D2Q5, =0.75 D2Q9, =0.75 D2Q5, =1.5 D2Q9, =1.5 D2Q5, =3.0 D2Q9, =3.0 Slope=-1

10-2

10-3

Neumann Condition

-1

10

100

r0

31

ACCEPTED MANUSCRIPT

Fig. 13 Relative error at various relaxation parameters (  ) for the temperature (a), interior temperature gradients (b), and boundary flux (c) versus mesh resolution r0   Nr  1  x for the heat conduction inside a circle with Neumann condition, obtained from the curved boundary implementations for both lattice models.

10-2 D2Q5, =0.75 D2Q9, =0.75 D2Q5, =1.5 D2Q9, =1.5 D2Q5, =3.0 D2Q9, =3.0 Slope=-2

10-3

10-4

100

r0

D2Q5, =0.75 D2Q9, =0.75 D2Q5, =1.5 D2Q9, =1.5 D2Q5, =3.0 D2Q9, =3.0 Slope=-2

10-3

10-4

10-5 10

100

r0

(d)

(c)

Robin Condition 10

10-2

AN US

10

10

Robin Condition -1

CR IP T

-1

10-5

-1

D2Q5, =0.75 D2Q9, =0.75 D2Q5, =1.5 D2Q9, =1.5 D2Q5, =3.0 D2Q9, =3.0 Slope=-1.5

10-3

10-4

10-5 10

M

10-2

ED

Relative Error for Interior Gradient

Relative Error for Boundary Value

10

(b)

Robin Condition

100

Relative Error for Boundary Flux

Relative Error for Temperature

(a)

Robin Condition

10-1

10-2

10-3

10-4

D2Q5, =0.75 D2Q9, =0.75 D2Q5, =1.5 D2Q9, =1.5 D2Q5, =3.0 D2Q9, =3.0 Slope=-1.5 Slope=-1

10

r0

100

PT

r0

Fig. 14 Relative error at various relaxation parameters (  ) for temperature (a), boundary value (b), interior

CE

temperature gradient (c), boundary flux (d) versus mesh resolution r0   Nr  1  x for the heat conduction inside a circle with linear Robin condition, obtained from the curved boundary implementations for both lattice

AC

models.

Comparisons based on the Neumann condition test is given in Fig. 13, with Figs. 13a, 13b and 13c showing respectively the errors of the interior temperature, interior temperature gradient and the boundary flux. Consistent with the Dirichlet condition test, the convergence trend for the temperature distribution in Fig. 13a is second-order, while the convergence order for the errors of the interior gradients in Fig. 13b and the boundary flux in Fig. 13c are lower than second-order. 32

ACCEPTED MANUSCRIPT

Besides, as demonstrated in Fig. 13a, the interior temperature distribution from the D2Q9 lattice model is more accurate than that of the D2Q5 model for all the test cases. As for the temperature gradients, the accuracy of the two lattice models for the interior gradients are similar with each other in Fig. 13b, while the errors for the boundary flux from the D2Q9 lattice model are smaller than those of the D2Q5 model for cases of  =0.75 and 1.5 in Fig. 13c.

CR IP T

Lastly, the lattice model effects on the Robin condition treatments at curved boundaries are

evaluated by the linear Robin condition test defined by Eqs. (57) and (60), as depicted in Fig. 14. Again, the second-order accuracy of the obtained temperature distribution and the boundary temperature values from both lattice models are confirmed by the convergence trends of the corresponding relative errors, as respectively shown in Figs. 14a and 14b. Besides, the

AN US

degradation of the accuracy for the solutions of the temperature gradients, due to the first-order approximation Eq. (36b), is also demonstrated in Fig. 14c for the interior gradients and Fig. 14d for the boundary flux. Moreover, although the accuracy of the interior temperature gradients and the boundary flux resulted from both lattice models agree well with each other in Figs. 14c and 14d, respectively, the D2Q9 model is generally more accurate than the D2Q5 model for the

M

entire temperature distribution and the boundary temperature values.

Therefore, it is concluded from the comparisons in Figs. 12a, 13a and 14a for the temperature

ED

distribution and Figs. 12b and 14b for the boundary temperature values that the D2Q9 lattice model is generally more accurate than the D2Q5 model in terms of determining the temperature

CE

PT

field for the curved boundary treatments due to the additional diagonal directions.

3.5 2D Natural Convection in a square Enclosure with a circular cylinder in the center

AC

Lastly, the lattice model effects of the present boundary scheme are further evaluated by a more practical benchmarking test on the coupling between the macroscopic velocity and the temperature fields. Boundary conditions are defined as follows: the stationary non-slip condition is applied at both the four outer straight walls and the inner circular surface, while a higher temperature Th and lower temperature Tc are prescribed at the inner cylinder surface and the outer straight walls, respectively. A body force term  az    g T  Tm  is added to the

33

ACCEPTED MANUSCRIPT

momentum equation by invoking the Boussinesq approximation, where Tm  Th  Tc  2 is the reference temperature, g is the gravitational acceleration and  is the thermal-expansion coefficient. The dimensionless parameters in the present test are the Prandtl Number Pr   D and the Rayleigh number Ra  g  Th  Tc  L3 D , with  representing the kinematic viscosity

CR IP T

of the flow, L  5r0 the length of the outer boundary and r0 the radius of the inner circular cylinder. The present boundary schemes based on the two lattice models are applied to

implement the Dirichlet temperature conditions, while our previous second-order curved

boundary scheme [37] and the non-equilibrium extrapolation scheme by Guo [38] are introduced

AC

CE

PT

ED

M

AN US

for the hydrodynamic conditions at the inner surface and the outer straight walls, respectively.

Fig. 15 Isotherms in the region between the circular cylinder and the cavity walls obtained from the present boundary scheme based on the D2Q9 lattice model for Ra values of: (a) 103; (b) 104; (c) 105; and (d) 106.

34

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

PT

Fig. 16 Streamlines in the region between the circular cylinder and the cavity walls obtained from the present

CE

boundary scheme based on the D2Q9 lattice model for Ra values of: (a) 103; (b) 104; (c) 105; and (d) 106.

For a fixed Pr  0.71 , test cases of Ra  103 , 104 , 105 and 106 are considered here to first verify

AC

the validity of the boundary scheme based on the D2Q9 model. With the mesh resolution of

N x  201 , convergent results obtained from the boundary schemes based on both lattice models

agree well with each other. In particular, the temperature distributions and the flow structures resulting from the D2Q9 model based boundary scheme are demonstrated in Figs. 15 and 16, and accord well with the reference results in the literature[39-41]. Furthermore, the surface-averaged

35

ACCEPTED MANUSCRIPT

Nusselt number Nu obtained from the boundary schemes based on the two lattice models are compared in Table 1 using the following definition Nu 

1 L T Nu ds , with Nu   L 0 n

(62) wall

CR IP T

Quantitative comparisons in Table 1 demonstrate that: (i) the surface-averaged Nusselt number

Nu obtained from the boundary schemes based on both lattice models accords well with the reference results[39-41] for all the test cases; and (ii) the differences caused by the two lattice models are negligible, and thus the present boundary schemes based on the two lattice models are equivalently accurate for the test considered here.

Li et al.[39]

M

103 104 105 106

Present Boundary Scheme D2Q5 D2Q9 Model Model 3.168 3.168 3.225 3.226 4.904 4.904 8.856 8.861

3.167 3.226 4.912 8.966

Reference Results Peng et al. [40] 3.24 4.84 8.75

Kim et al.[41] 3.396 3.414 5.138 9.390

ED

Ra

AN US

Table 1 The surface-averaged Nusselt number Nu for various Ra values.

4. Conclusions

PT

In this work, the lattice model effects on the boundary implementations for the CDE LB method are investigated by firstly extending our previously proposed second-order boundary scheme to

CE

the D2Q9 lattice model. Theoretical derivations demonstrate that, the present boundary scheme is flexible for both lattice models, and its extension to the D2Q9 model is straightforward while

AC

preserving the second-order accuracy and its advantages over other existing second-order boundary schemes [21, 27]. Moreover, the original D2Q5 boundary treatments are viewed as directly reduced implementations of the boundary scheme based on the D2Q9 model, and the straight boundary implementations are greatly simplified by the D2Q5 lattice model, since the only one unknown distribution function is along the inner normal direction and the tangential derivative at the boundary nodes is free of the practical boundary implementations. However, for the Dirichlet condition at a straight boundary, the constraint on the tangential derivative 36

 is  

ACCEPTED MANUSCRIPT

directly derived as an ‘implicit constraint’, and can be naturally included by the D2Q9 lattice model but ignored by the D2Q5 model. Therefore, a modified Dirichlet condition treatment at a straight boundary is developed for the D2Q5 model by additionally updating the distribution functions in the tangential direction. Furthermore, the lattice model effects on the boundary implementations are numerically

CR IP T

evaluated by directly comparing the errors from the two lattice models with respect to the corresponding analytical solutions. It is demonstrated that: (i) the accuracy of boundary

implementations based on both lattice models agrees well with each other for the 1D unsteadystate test; (ii) the D2Q9 model is slightly more accurate than the D2Q5 model for the Dirichlet condition treatments at straight boundaries, and the implementation for the D2Q5 model is

AN US

effectively improved by the modified scheme; (iii) for the treatments of the Neumann and Robin condition at the straight boundaries, the D2Q5 model is demonstrated to be generally more accurate than the D2Q9 model; and (iv) the D2Q9 model is helpful to improve the accuracy of the obtained temperature distribution for curved boundary treatments. Moreover, degradation of the accuracy is demonstrated in 2D straight boundary treatments based on the D2Q9 model for

M

the third test, and we attribute such loss of accuracy to the improper treatments of the corner points.. Therefore, the D2Q5 lattice model, with the modified treatments for the Dirichlet

ED

condition, is highly advantageous for straight boundary implementations. Besides, for curved boundary treatments, the D2Q5 model can also be viewed as a compromised option based on the

PT

tradeoff between computational efficiency and numerical accuracy.

CE

Acknowledgements

We would like to thank the financial support from the National Research Foundation (NRF),

AC

Prime Minister’s Office, Singapore under its Campus for Research Excellence and Technological Enterprise (CREATE) program, and the Joint Singapore-Germany Research Project Fund (SGPPROG3-019). This work is also financially supported by the National Natural Science Foundation of China (No. 11572062 and No. 41672292) and the Program for Changjiang Scholars and Innovative Research Team in University (No. IRT13043).

37

ACCEPTED MANUSCRIPT

Reference

AC

CE

PT

ED

M

AN US

CR IP T

[1] Aidun CK, Clausen JR. Lattice-Boltzmann Method for Complex Flows. Annual Review of Fluid Mechanics. 2010;42:439-72. [2] Li Q, Luo KH, Kang QJ, He YL, Chen Q, Liu Q. Lattice Boltzmann methods for multiphase flow and phase-change heat transfer. Progress in Energy and Combustion Science. 2016;52:62-105. [3] Wang J, Chen L, Kang Q, Rahman SS. The lattice Boltzmann method for isothermal micro-gaseous flow and its application in shale gas flow: A review. International Journal of Heat and Mass Transfer. 2016;95:94-108. [4] Chen L, Kang Q, Mu Y, He Y-L, Tao W-Q. A critical review of the pseudopotential multiphase lattice Boltzmann model: Methods and applications. International Journal of Heat and Mass Transfer. 2014;76:210-36. [5] Chen L, Kang Q, Tang Q, Robinson BA, He Y-L, Tao W-Q. Pore-scale simulation of multicomponent multiphase reactive transport with dissolution and precipitation. International Journal of Heat and Mass Transfer. 2015;85:935-49. [6] Kang J, Prasianakis NI, Mantzaras J. Thermal multicomponent lattice Boltzmann model for catalytic reactive flows. Physical Review E. 2014;89:063310. [7] Chen L, Kang Q, Carey B, Tao W-Q. Pore-scale study of diffusion–reaction processes involving dissolution and precipitation using the lattice Boltzmann method. International Journal of Heat and Mass Transfer. 2014;75:483-96. [8] Gillissen JJJ, Looije N. Boundary conditions for surface reactions in lattice Boltzmann simulations. Physical Review E. 2014;89:063307. [9] Pedersen J, Jettestuen E, Vinningland JL, Hiorth A. Improved Lattice Boltzmann Models for Precipitation and Dissolution. Transport in Porous Media. 2014;104:593-605. [10] Li X, Chen J, Xu M, Huai X, Xin F, Cai J. Lattice Boltzmann simulation of catalytic reaction in porous media with buoyancy. Applied Thermal Engineering. 2014;70:586-92. [11] Aursjø O, Jettestuen E, Vinningland JL, Hiorth A. An improved lattice Boltzmann method for simulating advective–diffusive processes in fluids. Journal of Computational Physics. 2017;332:363-75. [12] Chai Z, Shi B. A novel lattice Boltzmann model for the Poisson equation. Appl Math Model. 2008;32:2050-8. [13] Chai Z, Zhao TS. Lattice Boltzmann model for the convection-diffusion equation. Physical Review E. 2013;87:063309. [14] Guo Z, Zheng C, Shi B, Zhao TS. Thermal lattice Boltzmann equation for low Mach number flows: Decoupling model. Physical Review E. 2007;75:036704. [15] Huang R, Wu H. A modified multiple-relaxation-time lattice Boltzmann model for convection– diffusion equation. Journal of Computational Physics. 2014;274:50-63. [16] Li L, Mei R, Klausner JF. Multiple-relaxation-time lattice Boltzmann model for the axisymmetric convection diffusion equation. International Journal of Heat and Mass Transfer. 2013;67:338-51. [17] Li L, Mei R, Klausner JF. Lattice Boltzmann models for the convection-diffusion equation: D2Q5 vs D2Q9. International Journal of Heat and Mass Transfer. 2017;108:41-62. [18] Yi H-L, Yao F-J, Tan H-P. Lattice Boltzmann model for a steady radiative transfer equation. Physical Review E. 2016;94:023312. [19] Yoshida H, Nagaoka M. Multiple-relaxation-time lattice Boltzmann model for the convection and anisotropic diffusion equation. Journal of Computational Physics. 2010;229:7774-95. [20] Zheng L, Shi B, Guo Z. Multiple-relaxation-time model for the correct thermohydrodynamic equations. Physical Review E. 2008;78:026705. [21] Li L, Mei R, Klausner JF. Boundary conditions for thermal lattice Boltzmann equation method. Journal of Computational Physics. 2013;237:366-95.

38

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

[22] Chen Q, Zhang X, Zhang J. Improved treatments for general boundary conditions in the lattice Boltzmann method for convection-diffusion and heat transfer processes. Physical Review E. 2013;88:033304. [23] Meng X, Guo Z. Boundary scheme for linear heterogeneous surface reactions in the lattice Boltzmann method. Physical Review E. 2016;94:053307. [24] Zhang T, Shi B, Guo Z, Chai Z, Lu J. General bounce-back scheme for concentration boundary condition in the lattice-Boltzmann method. Physical Review E. 2012;85:016701. [25] Zhang T, Shi B, Guo Z, Chai Z, Lu J. Erratum: General bounce-back scheme for concentration boundary condition in the lattice-Boltzmann method [Phys. Rev. E 85, 016701 (2012)]. Physical Review E. 2013;88:029903. [26] Yin X, Zhang J. An improved bounce-back scheme for complex boundary conditions in lattice Boltzmann method. Journal of Computational Physics. 2012;231:4295-303. [27] Huang J, Hu Z, Yong W-A. Second-order curved boundary treatments of the lattice Boltzmann method for convection–diffusion equations. Journal of Computational Physics. 2016;310:26-44. [28] Huang J, Yong W-A. Boundary conditions of the lattice Boltzmann method for convection–diffusion equations. Journal of Computational Physics. 2015;300:70-91. [29] Gebäck T, Heintz A. A Lattice Boltzmann Method for the Advection-Diffusion Equation with Neumann Boundary Conditions. Communications in Computational Physics. 2015;15:487-505. [30] Huang HB, Lu XY, Sukop MC. Numerical study of lattice Boltzmann methods for a convection– diffusion equation coupled with Navier–Stokes equations. Journal of Physics A: Mathematical and Theoretical. 2011;44:055001. [31] Zhang L, Yang S, Zeng Z, Chew JW. Consistent Second-order Boundary Implementations for Convection-Diffusion Lattice Boltzmann Method. Physical Review E. 2017, submitted. [32] Shi B, Guo Z. Lattice Boltzmann simulation of some nonlinear convection–diffusion equations. Computers & Mathematics with Applications. 2011;61:3443-52. [33] Shi B, Guo Z. Lattice Boltzmann model for nonlinear convection-diffusion equations. Physical Review E. 2009;79:016701. [34] Bhatnagar PL, Gross EP, Krook M. A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Physical review. 1954;94:511. [35] Grad H. Note on N‐dimensional hermite polynomials. Communications on Pure and Applied Mathematics. 1949;2:325-30. [36] Grad H. On the kinetic theory of rarefied gases. Communications on pure and applied mathematics. 1949;2:331-407. [37] Zhang L, Zeng Z, Xie H, Tao X, Zhang Y, Lu Y, et al. An alternative second order scheme for curved boundary condition in lattice Boltzmann method. Computers and Fluids. 2015:193. [38] Zhao-Li G, Chu-Guang Z, Bao-Chang S. Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method. Chinese Physics. 2002;11:366. [39] Li L, Mei R, Klausner JF. Heat Transfer Evaluation on Curved Boundaries in Thermal Lattice Boltzmann Equation Method. Journal of Heat Transfer. 2013;136:012403--14. [40] Peng Y, Chew YT, Shu C. Numerical simulation of natural convection in a concentric annulus between a square outer cylinder and a circular inner cylinder using the Taylor-series-expansion and leastsquares-based lattice Boltzmann method. Physical Review E. 2003;67:026701. [41] Kim BS, Lee DS, Ha MY, Yoon HS. A numerical study of natural convection in a square enclosure with a circular cylinder at different vertical locations. International Journal of Heat and Mass Transfer. 2008;51:1888-906.

39