Second-order curved boundary treatments of the lattice Boltzmann method for convection–diffusion equations

Accepted Manuscript Second-order curved boundary treatments of the lattice Boltzmann method for convection–diffusion equations

Juntao Huang, Zexi Hu, Wen-An Yong

PII: DOI: Reference:

S0021-9991(16)00009-7 http://dx.doi.org/10.1016/j.jcp.2016.01.008 YJCPH 6347

To appear in:

Journal of Computational Physics

Received date: Revised date: Accepted date:

2 November 2015 3 January 2016 8 January 2016

Please cite this article in press as: J. Huang et al., Second-order curved boundary treatments of the lattice Boltzmann method for convection–diffusion equations, J. Comput. Phys. (2016), http://dx.doi.org/10.1016/j.jcp.2016.01.008

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.

Second-order curved boundary treatments of the lattice Boltzmann method for convection-diffusion equations Juntao Huang Zhou Pei-Yuan Center for Applied Mathematics, Tsinghua University, Beijing 100084, China

Zexi Hu AML, Department of Engineering Mechanics, Tsinghua University, Beijing 100084, China

Wen-An Yong∗ Zhou Pei-Yuan Center for Applied Mathematics, Tsinghua University, Beijing 100084, China

Abstract In this paper, we present a kind of second-order curved boundary treatments for the lattice Boltzmann method solving two-dimensional convectiondiffusion equations with general nonlinear Robin boundary conditions. The key idea is to derive approximate boundary values or normal derivatives on computational boundaries, with second-order accuracy, by using the prescribed boundary condition. Once the approximate information are known, the second-order bounce-back schemes can be perfectly adopted. Our boundary treatments are validated with a number of numerical examples. The results show the utility of our boundary treatments and very well support our theoretical predications on the second-order accuracy thereof. The idea is quite universal. It can be directly generalized to 3-dimensional problems, multiple-relaxation-time models, and the Navier-Stokes equations. Keywords: Lattice Boltzmann method, Convection-diffusion equations, ∗

Corresponding author Email addresses: [email protected] (Juntao Huang), [email protected] (Zexi Hu), [email protected] (Wen-An Yong)

Preprint submitted to Journal of Computational Physics

January 11, 2016

Curved boundaries, Approximate boundary conditions, Second-order accuracy 1. Introduction The lattice Boltzmann method (LBM) is an effective and viable tool for simulating complex fluid flows. Historically, the method was evolved from the lattice gas automata [1, 2]. It is now well understood as a special discretization of the Boltzmann equation [3, 4]. The kinetic origin enables the method to naturally accommodate a variety of boundary conditions for flows with complex geometry. Besides its effectiveness and simplicity, the method has a clear and solid physical interpretation [5, 6]. These advantages make it have vast applications. Indeed, the method has been successfully used in a wide spectrum of areas including turbulent flows, microflows, multi-phase and multi-component flows, and particulate suspensions. It is becoming a serious alternative to traditional computational methods in fluid dynamics. We refer to [5, 7] for reviews of the method and its applications. In recent years, the LBM has been widely used to solve convectiondiffusion problems and most of the existing works treat Dirichlet and/or Neumann boundary conditions on straight boundaries. The interested reader is referred to [8, 9, 10] for detailed reviews thereabout. However, general Robin boundary conditions on complex boundaries do arise in many applications, such as chemical engineering and geology, and therefore have attracted a lot of attentions [8, 9, 10, 11, 12, 13]. Among these works, Refs. [8, 11, 12, 13] treated straight boundaries only. For curved boundaries, we only know Refs. [9, 10, 14, 15], which, except Ref. [15], are only first-order accurate for Neumann boundary conditions. In Ref. [15], the authors dealt with zero normal flux boundary conditions only and did not investigate the errors of temperature gradients, boundary values or fluxes. In addition, there are few works about nonlinear Robin boundary conditions. On the other hand, the standard boundary treatment for solid walls in the LBM is the so-called half-way bounce-back scheme [16, 17], which can be easily implemented especially for complex boundaries in porous media. With this boundary treatment, curved boundaries are approximated by stair steps (or called the zigzag boundary). However, as reported in Ref. [18], when applying this bounce-back idea to convection-diffusion equations with Neumann-type boundary conditions on curved boundaries, the numerical 2

solution does not converge to the exact solution as the mesh size tends to zero. Therefore, the authors of Ref. [18] presented a modified scheme with the intersection area accounted for to alleviate the discrepancy. Considerable improvements were achieved in their numerical results but the motivation was not explained clearly. In Ref. [19], we constructed a single-node secondorder boundary scheme which can deal with straight boundaries locating at any distance from the lattice nodes. When the boundary locates at half of the lattice grid, the scheme reduces to the half-way bounce-back scheme for general Robin boundary conditions (including Dirichlet and Neumann boundary conditions). In this work, we propose a kind of second-order curved boundary treatments for the lattice Boltzmann method solving two-dimensional convectiondiffusion equations with general nonlinear Robin boundary conditions. The treatments are based on the single-node bounce-back scheme constructed in Ref. [19]. The idea is to determine computational boundaries and derive approximate boundary information (values or normal derivatives), with second-order accuracy, by using the prescribed boundary condition. The computational boundaries are zigzag and the approximate boundary information is obtained only at the middle points of two adjacent lattice nodes with one of them in the computational domain another not (see Fig. 1). Once the approximate information are known, the second-order bounce-back schemes can be perfectly adopted. Our boundary treatments are numerically tested with a number of examples. They are (i) one-D transient convection-diffusion with Robin boundary conditions, (ii) Helmholtz equation in a square domain, (iii) steady heat conduction inside a circle with four types of boundary conditions (Dirichlet, Neumann, Robin and nonlinear Robin boundary conditions), and (iv) twoD time-dependent convection-diffusion problem in a quite irregular domain. We compute the relative errors between the numerical solutions and exact solutions as well as the gradient, boundary values and boundary fluxes. The numerical results show the utility of our boundary treatments and very well support our theoretical predications on the second-order accuracy thereof. Let us remark that our idea is quite universal. It can be directly generalized to 3-dimensional problems, multiple-relaxation-time (MRT) models, and other equations such as nonlinear convection-diffusion equations [20] and the Navier-Stokes equations. We will report our progresses in forthcoming publications. The paper is organized as follows. In Section 2 we introduce convection3

diffusion equations together with nonlinear Robin boundary conditions, present the lattice Boltzmann method, and review the half-way bounce-back scheme proposed in [19] for Robin boundary conditions. Section 3 contains our main ideas and boundary treatments together with some remarks. Numerical experiments are reported in Section 4. Some concluding remarks are given in Section 5 where we mention the possibilities to apply our idea to some important problems including the Navier-Stokes equations with moving boundaries. 2. Preliminaries Consider a convection-diffusion equation (CDE) for a physical quantity (concentration, temperature, etc.) C = C(t, x) defined for x ∈ Ω: ∂C + ∇ · (uC) = D∇2 C + S(C, t, x), ∂t

0 < t < T, x ∈ Ω,

(1)

with initial condition C(0, x) = C0 (x),

x ∈ Ω.

Here u = u(t, x) is a given velocity field, D is a diffusion coefficient, S = S(C, t, x) is a source term depending on the unknown C, and T is a prespecified time. On the boundary ∂Ω, a nonlinear Robin boundary condition is prescribed: ∂C α = f (C, t, x), 0 < t < T, x ∈ ∂Ω. (2) ∂n Here α is a constant, n is the outward unit normal vector to the boundary ∂Ω, and f is a given function of C, t and x. If f linearly depends on C, the Robin boundary condition is recovered: α 1 C + α2

∂C = α3 , ∂n

0 < t < T, x ∈ ∂Ω.

(3)

Here αk (k = 1, 2, 3) are given functions of t and x. 2.1. Lattice Boltzmann method For simplicity, we work only with the two-dimensional five-velocity (D2Q5) BGK model (see e.g. [18]). However, our ideas do not depend on the specific

4

model and can be easily implemented for Multiple-Relaxation-Time (MRT) models [18]. In the D2Q5 model, the five discrete velocities are  (0, 0), i = 1, ci = (cos[(i − 2)π/2], sin[(i − 2)π/2]), i = 2, 3, 4, 5. The evolution equation for the i-th distribution function gi = gi (t, x) reads as 1 (eq) gi (t + δt , x + cδt ci ) − gi (t, x) = − (gi (t, x) − gi (t, x)) + ωi δt S(C, t, x) (4) τ for i = 1, 2, · · · , 5. Here δt is a time step, c = h/δt is the lattice speed with h the lattice size, τ is the dimensionless relaxation time given by τ=

3δt 1 D+ 2 h 2

with D the diffusion coefficient in the convection-diffusion equation (1), (eq) (eq) gi = gi (t, x) is a local equilibrium distribution function (eq)

gi

= C(ωi +

u · ci ) 2c

(5)

with the macroscopic quantity C = C(t, x) defined as C=

5 

gi ,

i=1

and the weight coefficients ωi are  1/3, ωi = 1/6,

i = 1, i = 2, 3, 4, 5.

2.2. Bounce-back boundary conditions In this part, we review the modified bounce-back boundary scheme proposed in Ref. [19] for general Robin boundary conditions (3) with straight boundaries. It is of second-order accuracy. For the boundary intersecting with a lattice link at its middle point, the scheme is gα (t, x) =

−ph + q 2ωα α3 h gα¯ (t, x + hcα¯ ) + ph + q ph + q 5

(6)

with parameters p and q given by p = α1 −

u · cα α2 , D

q=

α 2 h2 . 3δt D

Here α1 , α2 and α3 are parameters in the Robin boundary condition (3). For Dirichlet boundary conditions, we have α1 = 1 and α2 = 0. Then the scheme (6) reduces to [18, 8] gα (t, x) = −gα¯ (t, x + hcα¯ ) + 2ωα α3 .

(7)

On the other hand, we have α1 = 0 and α2 = 1 for Neumann boundary conditions. Then the scheme (6) becomes gα (t, x) =

3δt (u · cα ) + h 6ωα Dδt gα¯ (t, x + hcα¯ ) + α3 . −3δt (u · cα ) + h −3δt (u · cα ) + h

(8)

3. Boundary treatments In this section, we present our treatments of the general nonlinear Robin boundary condition (2). As shown in Fig.1, the dot-dashed line (i.e. the boundary ∂Ω) is called the physical boundary. The black and white circles denote interior and exterior lattice nodes, respectively. The dashed line is parallel to the lattice grid piecewise and locates at half of the lattice link, which is called the computational boundary. Our main idea can be stated as: If we could get the boundary information (values or fluxes) on the computational boundary (white squares) precisely or with second-order accuracy, then we can achieve second-order accuracy by using scheme (6) (or (7),(8)). We will show how to obtain approximate boundary information on the computational boundary by using the prescribed boundary conditions. For definiteness, we take Point 2 (P2 ) as an example. Let Point 5 (P5 ) be the nearest point on the physical boundary ∂Ω to P2 . Such nearest points may not be unique but we will take one of them. If the physical boundary is smooth at P5 , it can be verify easily that the straight line linking P5 to P2 is orthogonal to the boundary. Namely, the unit outer normal vector n at P5 is equal to P 2 − P5 n = χΩ (P2 ) |P2 − P5 | with |P2 − P5 | the distance between P2 and P5 and χΩ (P2 ) = −1 if P2 ∈ Ω, otherwise χΩ (P2 ) = 1. Thus, the coordinate of P5 can be easily determined 6

Ϯ ϲ

ϵ

ϱ ϭ ϳ

ϯ

ϴ

ĞdžƚĞƌŝŽƌ ϰ

ŝŶƚĞƌŝŽƌ

Figure 1: Schematic map of the lattice and the curved boundary. Black circles: interior nodes; white circles: exterior nodes; dot-dashed lines: physical boundary; dashed lines: computational boundary; white squares: half-way nodes. The straight line linking P5 to P2 is orthogonal to the boundary at Point 5.

7

for a circle or ring domain and for a domain with straight boundaries. For a general (physical) boundary represented by f (x, y) = 0, the coordinate can be determined approximately with the formula (x5 , y5 ) ≈ (x2 , y2 ) −

f (x2 , y2 )∇f (x2 , y2 ) . |∇f (x2 , y2 )|2

(9)

It is not difficult to show with the Taylor expansion that this formula can approximate the exact coordinate by second-order accuracy. Next we obtain the boundary information for the following three different cases. We start with the Dirichlet boundary condition. 3.1. Dirichlet boundary condition In this case, the value of the unknown C is given on the physical boundary (the dot-dashed line). Then the value at P2 can be computed with C(tn , P2 ) = C(tn , P5 ) + (P2 − P5 ) · ∇C(tn , P5 ) + O(h2 ) ∂C (tn−1 , P1 )|P2 − P5 |χΩ (P2 ) + O(h2 ). = C(tn , P5 ) + ∂n Here C(tn , P5 ) is the exact value given by the Dirichlet boundary condition at P5 in the n-th time step, n is the unit outer normal vector at P5 , and we have used the facts that tn −tn−1 = O(h2 ), |P2 −P5 | = O(h) and |P1 −P5 | = O(h). Note that ∂C (t , P1 ) = n · ∇C(tn−1 , P1 ) and the gradient at lattice nodes ∂n n−1 can be obtained with the following formula [18]:   5 3 uδt 1 ∇C = C− ci g i . (10) τ h2 h i=1 Remark 1. The value of C at P2 can also be obtained by using the information of P9 as ∂C (tn , P9 )|P2 − P9 | + O(h2 ), ∂y ∂C (tn−1 , P1 )|P2 − P9 | + O(h2 ). = C(tn , P9 ) + ∂y

C(tn , P2 ) = C(tn , P9 ) +

As shown in Fig.1, P9 is the intersection of the physics boundary and the mesh grid. 8

3.2. Neumann boundary condition In this case, the normal derivative is given on the physical boundary and our goal is to obtain ∂C at P2 (see Fig.1). Notice that the identity ∂y ∂ ∂ ∂ = ny + nx , ∂y ∂n ∂t where nx (resp. ny ) is the x− (resp. y−) component of the outer normal vector n at P5 . Here and below, the tangential vector t on the boundary is defined as the outer normal vector n turning π/2 anticlockwise. That is, t = (tx , ty ) = (−ny , nx ). Thus, we only need to obtain the values of ∂C ∂n and ∂C at P . These two quantities can be computed by using the Taylor 2 ∂t expansion as follows. 1. Calculation of

∂C : ∂n

∂C ∂C ∂ 2C (tn , x5 )|P2 − P5 | + O(h2 ), (tn , P2 ) = (tn , P5 ) + ∂n ∂n ∂n2 2

2

2

2

∂ C where ∂∂nC2 := (n · ∇)2 C = n2x ∂∂xC2 + 2nx ny ∂x∂y + n2y ∂∂yC2 . The secondorder derivatives can be approximated to first-order accuracy by the difference quotations of first-order derivatives of the last time step:   ∂ 2C 1 ∂C ∂C (tn−1 , P3 ) − (tn−1 , P1 ) + O(h), (tn , P5 ) = ∂x2 h ∂x ∂x   ∂ 2C 1 ∂C ∂C (tn , P5 ) = (tn−1 , P1 ) − (tn−1 , P4 ) + O(h), ∂y 2 h ∂x ∂x   1 ∂C ∂C ∂ 2C (tn , P5 ) = (tn−1 , P1 ) − (tn−1 , P4 ) + O(h). ∂x∂y h ∂x ∂x

2. Calculation of

∂C : ∂t

∂C ∂C ∂ 2C (tn , P2 ) = (tn , P7 ) + (tn , P7 )|P2 − P7 | + O(h2 ), ∂t ∂t ∂n∂t ∂C ∂ 2C (tn−1 , P7 ) + (tn , P5 )|P2 − P7 | + O(h2 ), = ∂t ∂n∂t while ∂C (t , P7 ) is calculated by a linear interpolation of ∂C at P1 ∂t n−1 ∂t and P3 of the last time step, which guarantees a second-order accuracy 9

2

∂ C due to the diffusing scaling. For ∂n∂t (tn , P5 ), we recall that the normal derivatives along the boundary are known. Taking derivative on ∂C ∂n along the boundary, we get ∂ ∂ ∂C ( ) = (n · ∇C), ∂t ∂n ∂t ∂(∇C) ∂n + · ∇C, =n· ∂t ∂t ∂ 2C = + κt · ∇C, ∂n∂t ∂ 2C ∂C = +κ . ∂n∂t ∂t Here κ denotes the signed curvature at P5 [21]. For a general boundary represented by f (x, y) = 0, we have

κ=∇·( Hence

∂C ∂t

∇f (P5 ) ). |∇f (P5 )|

can be computed by

∂C ∂ ∂C ∂C ∂C (tn , P2 ) = (tn−1 , P7 )+( ( (tn , P5 ))−κ (tn−1 , P7 ))|P2 −P7 |. ∂t ∂t ∂t ∂n ∂t Remark 2. 1. When the intersecting point P7 (of the straight line passing P2 and P5 and the mesh grid) is outside the segment P1 P3 , say, P7 is between P3 and P8 , we calculate ∂C (t , P7 ) by a linear interpolation ∂t n at P3 and P8 . If the intersection point is outside P8 , then we directly approximate ∂C (t , P7 ) by ∂C (t , P8 ). Here extrapolation is always ∂t n ∂t n−1 avoided for the sake of numerical stability. From our numerical experiments, such a treatment can produce second-order accuracy results well. 2. For a general boundary, it may be not so easy to take the derivative ∂ ∂C along the curve. In this case, we may numerically calculate ∂t ( ∂n (tn , P5 )) as follows   ∂ ∂C 1 ∂C ∂C ( (tn , P5 )) = (tn , P5 + δht) − (tn , P5 ) ∂t ∂n δh ∂n ∂n (t , P5 + δht) with δ a fixed small number. Note that the vector n in ∂C ∂n n in the above formula is the unit outward normal vector at point (P5 + δht). In our last numerical example in Section 4, we will deal with an irregularly curved boundary with δ = 10−3 . 10

3.3. Nonlinear Robin boundary condition In this case, we have the boundary condition (2) with α = 0 and the goal is to obtain the normal derivative of C at P5 . To do this, we firstly calculate the value of C at P5 and the current time step by using the information at P1 in the last time step: C(tn , P5 ) = C(tn , P1 ) + (P5 − P1 ) · ∇C(tn , P1 ) + O(h2 ), = C(tn−1 , P1 ) + (P5 − P1 ) · ∇C(tn−1 , P1 ) + O(h2 ). Then the normal derivative at P5 is obtained by using the nonlinear Robin condition (2), ∂C 1 (tn , P5 ) = f (C(tn , P5 ), tn , P5 ). ∂n α Having this, we repeat the procedure in Subsection 3.2 dealing with the Neumann boundary condition and compute the value of ∂C at P2 . ∂y Remark 3. At this point, one may ask whether or not it is possible to compute C(tn , P2 ) by using the accurately obtained C(tn , P5 ) with the procedure in Subsection 3.1 dealing with the Dirichlet boundary condition. The answer is NO. In fact, our numerical experiments show that so doing is unstable. On the other hand, if it works, it means that the solutions could be obtained without using any boundary conditions. This is obviously false. 4. Numerical experiments In this section, we present some numerical results to validate our boundary treatments described in the previous section. In particular, we test the treatment of Robin boundary conditions by solving a time-dependent convection-diffusion problem defined in a quite irregular domain made up in [19]. 4.1. One-D transient convection-diffusion with Robin boundary conditions We start with a one-dimensional unsteady convection-diffusion equation ∂C ∂C ∂ 2C +u =D 2, ∂t ∂x ∂x

0 < x < L, t > 0,

with a blank initial condition C(t, x)|t=0 = 0, 11

0 < x < L,

y

γh

γh

1

O

1

x

Figure 2: Schematic depiction of the lattice for the one-D transient convection-diffusion with different γ. The spatial mesh size is h. The dots represent the lattice nodes in the interior domain. The Robin boundary condition and the Neumann boundary condition are applied at the left and right boundaries respectively. The periodic boundary conditions are imposed at the up and bottom boundaries.

a Robin boundary condition at the left endpoint: (uC − D

∂C )|x=0 = uCf , ∂x

t > 0,

and a Neumann boundary condition at the right endpoint: ∂C |x=L = 0, ∂x

t > 0.

Here u, D, Cf and L are all constants. A non-dimensional parameter is the Peclet number, which is defined as Pe = uL/D. In our numerical experiments reported below, we take D = 0.01, Cf = 50 and L = 1. To test our boundary treatments for two-dimensional problems, we extend this one-dimensional problem to a square domain with periodic boundary conditions imposed at the up and bottom boundaries of the domain and the specified boundary conditions applied at the left and right boundaries. As shown in Fig.2, the distance between the outermost lattice nodes and the left and right boundary is γh with h the spatial mesh size. 12

Relative error

10

10

10

10

10

-1

10

γ =0.01 γ =0.25 γ =0.5 γ =0.75 γ =0.99

-2

10

Relative error

10

Slope = 2

-3

-4

-6

10

-3

Mesh size

10

10

-2

(a)

γ =0.01 γ =0.25 γ =0.5 γ =0.75 γ =0.99

-3

Slope = 2

10

10

Pe = 1

-5

-2

-4

Pe = 10

-5

-6

10

-3

Mesh size

10

-2

(b)

Figure 3: The relative 2 -error vs. the mesh size with different γ and Pe numbers for 1D convection-diffusion at t = 1 with τ = 1.5. Left: Pe = 1; Right: Pe = 10.

When t is sufficiently small, the above problem has an asymptotic solution [22] 1 x − ut (x − ut)2 u2 t erfc( √ ) + exp(− ) C(t, x) =Cf 2 πD 4Dt 2 Dt (11) 1 ux u2 t ux x + ut

− (1 + + ) exp( ) erfc ( √ ) . 2 D D D 2 Dt In Ref.[19], this asymptotic solution was examined to be valid for t < 5 under the same parameters D, Cf and L. To see the accuracy of our boundary treatments, we examine the relative 2  -error in the whole computational domain:  (C num − C exa )2  . (12) E= (C exa )2 Here the summation is over all the interior lattice nodes, and C num and C exa denote the numerical and asymptotic solutions, respectively. Fig.3 shows the relative errors versus the mesh size at t = 1 in bilogarithm coordinates with different γ values and Pe numbers. For each γ we take 5 different meshes with node numbers 64, 128, 256, 512, 1024. The results show that the convergence orders of our new numerical scheme are all around 2 for Pe = 1 and Pe = 10 with γ = 0.01, 0.25, 0.50, 0.75 and 0.99. This numerical experiment shows that our boundary treatments produce second-order accurate results for a wide range of γ values. 13

4.2. Helmholtz equation in a square domain Next we use our boundary treatments to solve the 2D Helmholtz equation in a square domain: ∇2 C = (λ2 + μ2 )C,

0 < x < 1, 0 < y < 1,

with boundary conditions C = 0, 0 < x < 1, y = 1, C = exp(−λx), 0 < x < 1, y = 0, C = sinh(μ(1 − y))/sinh(μ), 0 < y < 1, x = 0, ∂C = 0, 0 < y < 1, x = 1. λC + ∂x Here λ and μ are constants. The exact solution to the Helmholtz equation is C exact (x, y) = exp(−λx)sinh(μ(1 − y))/sinh(μ).

(13)

In this numerical experiment, we take λ = μ = 1. The mesh grid is similar to that of the first numerical test, except that the distance between the outermost lattice nodes and the up/bottom boundaries is also γh. The reader is also referred to Ref. [19]. We set C = 0 at the beginning of the iteration. If the relative 2 -error between every 1000 time steps is smaller than a given tolerance (it is 10−10 in this numerical experiment), the current numerical solution will be output. To gain further insights into the boundary treatments on the accuracy, we define the relative errors at the boundary for C and its normal derivative, and the errors for the gradient ∇C as

num − C exact |2 boundary nodes |C

Eboundary = , (14) exact |2 |C boundary nodes

num − n · ∇C exact |2 boundary nodes |n · ∇C

, (15) Eboundary flux = 2 exact | boundary nodes |n · ∇C



and Egrad =

interior nodes



|∇C num − ∇C exact |2

exact |2 interior nodes |∇C

14

.

(16)

-2

10

γ =0.01 γ =0.25 γ =0.5

10

10

10

10

γ =0.75

10

γ =0.99 Slope = 2

-4

10

τ =0.75

-5

-6

10

-2

Mesh size

10

-2

10

γ =0.25

Relative error

10

10

10

10

10

-1

γ =0.01

10

10

Relative error of interior gradient

Relative error

10

-3

Relative error of interior gradient

10

-3

γ =0.5 γ =0.75

10

γ =0.99 -4

Slope = 2

10

-5

10

τ =1.5

-6

-7

10

-2

Mesh size

10

10

-1

(a)

-2

γ =0.01 γ =0.25 γ =0.5 -3

γ =0.75 γ =0.99 Slope = 2

-4

-5

τ =0.75

-6

10

-2

Mesh size

10

-1

-2

γ =0.01 γ =0.25 γ =0.5 -3

γ =0.75 γ =0.99 Slope = 2

-4

-5

τ =1.5

-6

10

-2

Mesh size

10

-1

(b)

Figure 4: Relative 2 -errors E and Egrad versus the mesh size h for the Helmholtz equation in a square domain with different parameters γ and τ . Left column: relative 2 -errors E with τ = 0.75 and τ = 1.5. Right column: relative 2 -errors Egrad with τ = 0.75 and τ = 1.5.

Here C exact denotes the analytic solution given by Eq.(13); n in Eq.(15) is the unit outer normal vector at the boundary nodes; C num and its flux at boundary nodes are computed by linear extrapolations of interior nodes (or by Eqs. (19)(20) in Ref.[19] for straight boundaries); and the gradient of numerical solutions in interior lattice nodes is computed by Eq.(10). Presented in Fig.4 are the relative 2 -errors E (defined by Eq.(12)) and Egrad (defined by Eq.(16)) between the exact solution and the output versus the mesh size for different γ and τ . It can be observed that the relative errors are about order 10−6 to 10−4 when the mesh size is around 10−2 . The slopes and thus the convergence orders are all around two for a wide range of γ values (γ = 0.01, 0.25, 0.5, 0.75, 0.99). Fig.5 is for the relative errors at the boundary for C and its normal derivative. It gives the same conclusion that the convergence orders are all around two. Moreover, the half-way case (i.e. 15

10

10

10

10

Relative error of boundary value

10

10

10

10

10

10

10

γ =0.01 γ =0.25 -3

γ =0.5 γ =0.75

10

γ =0.99 -4

Relative error of boundary flux

10

-2

Slope = 2

10

-5

τ =0.75

-6

-7

10

-2

Mesh size

10

10

10

-1

-2

10

γ =0.01

Relative error of boundary flux

Relative error of boundary value

10

γ =0.25 -3

γ =0.5 γ =0.75

10

γ =0.99 -4

Slope = 2

10

-5

10

τ =1.5

-6

-7

10

-2

Mesh size

10

10

-1

(a)

-2

γ =0.01 γ =0.25 γ =0.5 -3

γ =0.75 γ =0.99 Slope = 2

-4

-5

τ =0.75 -6

10

-2

Mesh size

10

-1

-2

γ =0.01 γ =0.25 γ =0.5 -3

γ =0.75 γ =0.99 Slope = 2

-4

-5

τ =1.5 -6

10

-2

Mesh size

10

-1

(b)

Figure 5: Relative 2 -errors Eboundary and Eboundary flux versus the mesh size h for the Helmholtz equation in a square domain with different parameters γ and τ . Left column: relative 2 -errors Eboundary with τ = 0.75 and τ = 1.5. Right column: relative 2 -errors Eboundary flux with τ = 0.75 and τ = 1.5.

16

γ = 0.5) produces much better results for boundary values. These numerical tests again validate our boundary scheme for straight boundaries. 4.3. Steady heat condution inside a circle In this part, we study steady heat conduction inside a circle: ∇2 C = 0,

0 ≤ r < r0 , 0 ≤ θ < 2π,

with four different types of boundary conditions. They are a Dirichlet boundary condition C = cos(kθ) at r = r0 , a Neumann boundary condition ∂C k = cos(kθ) at r = r0 , ∂n r0 a Robin boundary condition C+

∂C k = (1 + ) cos(kθ) at r = r0 , ∂n r0

and a nonlinear Robin boundary condition C2 +

∂C k = cos2 (kθ) + cos(kθ) at r = r0 . ∂n r0

Here r0 and k are constants. These four problems have the same exact solution [10] C(r, θ) = (r/r0 )k cos(kθ). In our numerical experiments, we set r0 = 0.4 and k = 4. 4.3.1. Classical bounce-back schemes Before solving this problem with our new boundary treatments, we use the classical bounce-back schemes (7)/(8) dealing with straight boundaries to solve this problem with a curved boundary. The results are presented in Fig.6. It is observed that the numerical solution converges to the exact solution with order around 1 to 1.5 for the Dirichlet boundary condition but does not converge for the Neumann boundary condition. To explain this non-convergence phenomenon, we notice that, when using the bounceback schemes, we essentially approximate the curved boundary with a zigzag 17

100

Dirichlet boundary condition Neumann boundary condition Slope = 1 Slope = 1.5

Relative error

10-1

10-2

10-3

10-4

10-2

Mesh size

Figure 6: Classical bounce-back schemes. Relative 2 -errors versus the mesh size h for the steady heat conduction inside a circle with Dirichlet/Neumann boundary conditions. Parameter: τ = 1.5.

boundary. For this heat conduction problem with the flux boundary condition, the energy which flows into the domain through its boundary depends on the length of the boundary. However, no matter how small the mesh size is, the length of the computational boundary (a zigzag boundary) is 8r0 which is different from that of the physical boundary 2πr0 . Similar numerical phenomena are also reported in Ref.[18]. These indicate that the classical bounce-back schemes are not reliable at least in solving Neumann boundary problems. 4.3.2. Dirichlet boundary condition Fig.7 illustrates the relative errors E, Egrad and Eboundary flux versus the mesh size h with different τ . For the relative errors at the curved boundary in Eqs.(14) and (15), we approximate C and its normal derivative at the boundary nodes by linear extrapolations of the corresponding values at interior nodes. When the mesh resolution is about 10−2 with τ = 0.75, the relative errors E, Egrad and Eboundary flux are around 1 × 10−3 , 2 × 10−3 and 6 × 10−3 , respectively. The relative error E has a convergence order around 2. For larger τ values, the interior gradient and boundary flux both converge with order 2, while for smaller τ values (τ = 0.75), they only converge with order 1.5 and 1, respectively. Also, the errors all get smaller when τ get smaller with the same mesh size. In Ref.[10], this problem with the same parameters was simulated and 18

10

0

τ =0.75 τ =1.5

Relative error

10

10

10

10

10

-1

Slope = 2

-2

-3

-4

-5

10

Relative error of interior gradient

10

10

10

10

10

Relative error of boundary flux

10

10

10

Mesh size

-2

0

τ =0.75 τ =1.5 τ =3 -1

Slope = 1.5 Slope = 2

-2

-3

-4

10

10

10

τ =3

Mesh size

-2

0

τ =0.75 τ =1.5 τ =3 -1

Slope = 1 Slope = 2

-2

-3

-4

10

Mesh size

-2

Figure 7: Relative 2 -errors E, Egrad and Eboundary flux versus the mesh size h for the steady heat conduction in a circle domain with a Dirichlet boundary condition and with different τ .

19

101 100

Relative error

10-1

Present scheme Scheme 2 in Li et al Slope = 1 Slope = 2

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

Relative error of interior gradient

101 100

10-2

Mesh size

10-1

Present scheme Scheme 2 in Li et al Slope = 1 Slope = 2

10-1 10-2 10-3 10-4 10-5

Relative error of boundary flux

10-1

10-2

Mesh size

10-1

Present scheme Scheme 2 in Li et al Slope = 1

10-2

10-3

10-2

Mesh size

10-1

Figure 8: Comparison in relative 2 -errors E, Egrad and Eboundary flux with schemes in Li et al. [10] for the steady heat conduction inside a circle with a Dirichlet boundary condition. Blue squares: our scheme; green triangles: Scheme 2 in Li et al. [10] (data from Fig.19, Fig.20 and Fig.21 in Ref.[10]). Parameter: τ = 0.75.

20

a second-order scheme was used to treat the Dirichlet boundary condition on the curved boundary. The authors also observed the degradation of the convergence orders for the derivative to physical quantities (See Fig.20 and Fig.21 in Ref.[10]). This is an interesting numerical phenomenon but we do not know the reason either. Fig.8 gives the relative 2 -errors E, Egrad and Eboundary flux of our scheme and the scheme in Ref.[10]. The data are extracted from Figs.19–21 of Ref.[10]. The convergence orders of our present scheme and the scheme in [10] are both around 2 for C, and around 1 for the gradient and boundary flux. 4.3.3. Neumann boundary condition Presented in Fig.9 are the relative 2 -errors E, Egrad and Eboundary versus the mesh size h. It is observed that the convergence orders are all around 2. Moreover, the relaxation time τ has an impact on the errors. A comparison with Li et. al.’s schemes is given in Fig.10. One can see that our scheme converges much faster (convergence order 2) than theirs (convergence order 1). Furthermore, our schemes have no oscillations. For this specific problem, it seems that theirs is better than ours in accuracy with low resolutions but ours is better with high resolutions. 4.3.4. Robin boundary condition Fig.11 illustrates the relative 2 -errors E, Egrad , Eboundary and Eboundary flux versus the mesh size h with different τ . They all have second-order accuracy. Note that the errors are all of order 10−3 to 10−2 when τ = 0.75 and h = 1/64. 4.3.5. Nonlinear Robin boundary condition The numerical results with different meshes h and τ are plotted in Fig.12. It shows that our scheme behaves very well for different τ and the convergence orders are all around 2. The rightmost point of the red line (τ = 3 and h = 1/16) in Fig.12 is missing as it is unstable with the corresponding parameter. This phenomenon may be caused by the nonlinearity of the boundary condition. 4.4. Time-dependent convection-diffusions in an irregular domain To examine the applicability of our boundary treatments in dealing with general curved boundaries, we solve the problem made up in Ref.[19] with an irregular domain. The boundary consists of zeros of the polynomial P (x, y) = x4 − 5x2 − 3x + 2y 4 − 6y 3 − y − 1 21

10

1

τ =0.75 τ =1.5 10

0

τ =3

Relative error

Slope = 2 10

10

10

10

10

Relative error of interior gradient

10

10

10

10

-1

-2

-3

-4

-5

10

-1

τ =0.75 τ =1.5 τ =3 Slope = 2 -2

-3

-4

10

Relative error of boundary value

10

10

10

10

Mesh size

-2

1

τ =0.75 τ =1.5

10

10

Mesh size

-2

0

τ =3 Slope = 2

-1

-2

-3

-4

10

Mesh size

-2

Figure 9: Relative 2 -errors E, Egrad and Eboundary versus the mesh size h for the steady heat conduction in a circle domain with a Neumann boundary condition and with different τ.

22

100

Relative error

10-1

Present scheme Coupled Scheme 2 in Li et al Slope = 1 Slope = 2

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

Relative error of interior gadient

100

10-1

10-2

Mesh size

10-1

Present scheme Coupled Scheme 2 in Li et al Slope = 1 Slope = 2

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

Relative error of boundary value

100

10-1

10-2

Mesh size

10-1

Present scheme Coupled Scheme 2 in Li et al Slope = 1 Slope = 2

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

10-2

Mesh size

10-1

Figure 10: Comparison in relative 2 -errors E, Egrad and Eboundary with the schemes in Li et al. [10] for the steady heat conduction inside a circle with a Neumann boundary condition. Blue squares: our scheme; green triangles: Coupled Scheme 2 in Li et al. [10] (data from Fig.22, Fig.23 and Fig.24 in Ref.[10]). Parameter: τ = 0.75.

23

1

10

Relative error of boundary value

10

τ =0.75 τ =1.5 10

0

τ =3

10

10

10

10

-1

10

-2

10

-3

-5

10

Relative error of interior gradient

10

10

10

10

10

10

-4

10

Mesh size

-2

0

τ =1.5 τ =3 -1

τ =1.5

Slope = 2

10

-2

10

-3

10

-4

10

10

Mesh size

-2

(a)

0

τ =3 Slope = 2

-1

-2

-3

-4

10

10

τ =0.75

Relative error of boundary flux

Relative error

10

τ =0.75

10

Slope = 2

1

Mesh size

-2

0

τ =0.75 τ =1.5 τ =3 -1

Slope = 2

-2

-3

-4

10

Mesh size

-2

(b)

Figure 11: Relative 2 -errors E, Egrad , Eboundary and Eboundary flux versus the mesh size h for the steady heat conduction in a circle domain with a Robin boundary condition and with different τ .

24

1

10

Relative error of boundary value

10

τ =0.75 τ =1.5 10

0

τ =3

Relative error

Slope = 2 10

10

10

10

-1

-2

10

-3

10

10

-4

-5

10

Relative error of interior gradient

10

10

10

10

10

10

Mesh size

-2

0

τ =1.5 τ =3 -1

τ =0.75 τ =1.5 τ =3 -1

Slope = 2

10

-2

10

-3

10

-4

10

10

Mesh size

-2

(a)

Slope = 2

-2

-3

-4

10

10

τ =0.75

Relative error of boundary flux

10

0

Mesh size

-2

0

τ =0.75 τ =1.5 τ =3 -1

Slope = 2

-2

-3

-4

10

Mesh size

-2

(b)

Figure 12: Relative 2 -errors E, Egrad , Eboundary and Eboundary flux versus the mesh size h for the steady heat conduction in a circle domain with a nonlinear Robin boundary condition and with different τ . The rightmost point of the red line (τ = 3 and h = 1/16) is missing as it is unstable with the corresponding parameter.

25

4 3

y

2 1 0 −1 −2 −3

−2

−1

0

1

2

3

x Figure 13: Schematic depiction of the computational domain with an irregular boundary.

and the computational domain is Ω = {(x, y) ∈ R2 |P (x, y) < 0} (see Fig.13). The time-dependent convection-diffusion equation is ∂t C + ∇ · (uC) = D∇2 C + S,

(x, y) ∈ Ω,

t > 0,

with a blank initial condition C = 0,

(x, y) ∈ Ω,

t=0

and a Robin boundary condition C+

∂C = α3 (t, x, y), ∂n

(x, y) ∈ ∂Ω,

t > 0.

Here the velocity field u = (1, 0)T , D is a positive constant, the source term is S = S(t, x, y) = (x3 + y 2 )ω cos(ωt) + (3x2 − D(6x + 2)) sin(ωt), and α3 = (x3 + y 2 + 3nx x2 + 2ny y) sin(ωt) with nx and ny the x- and y-components of the outward unit normal vector n, and ω a constant. This problem has an analytic solution C(t, x, y) = (x3 + y 2 ) sin(ωt). 26

0

10

Relative error of boundary value

10

τ =0.75 τ =1.5

10

10

10

10

-1

Relative error of interior gradient

10

10

10

10

10

Slope = 2

-2

10

-3

10

-4

10

-5

10

Mesh size 10

10

10

τ =3

-1

0

τ =1.5 -1

τ =3

10

Slope = 2

-2

10

-3

10

-4

-5

τ =0.75 τ =1.5 -1

10

10

Mesh size 10

-1

(a)

τ =3 Slope = 2

-2

-3

-4

-5

Mesh size 10

10

τ =0.75

Relative error of boundary flux

Relative error

10

0

-1

0

τ =0.75 τ =1.5 -1

τ =3 Slope = 2

-2

-3

-4

-5

Mesh size 10

-1

(b)

Figure 14: Relative 2 -errors E, Egrad , Eboundary and Eboundary flux versus the mesh size h at time T = 1 for the time-dependent convection-diffusion in an irregular domain with different τ .

In our numerical experiment, we take D = ω = 1 as in Ref.[19], place the computational domain in a square [−3, 3] × [−2, 4], and set a uniform N × N mesh on the square. The numerical results with different meshes are plotted in Fig.14. It shows that our scheme behaves very well for different τ . The convergence orders are all around 2. To figure out the effects of the relaxation time τ on the relative 2 -error, we plot the error versus τ at T = 1 with the lattice number N = 64 in Fig.15. From Fig.15, it is observed that the relative errors get larger when τ gets larger (although there will be less calculation costs). Thus, we recommend that the relaxation time τ should not be too large (τ < 3).

27

0.040

Physical quantity 0.035

Interior gradient Boundary value Boundary flux

Relative error

0.030

0.025

0.020

0.015

0.010

0.005

0.000 0.5

1.0

1.5

2.0

2.5

3.0

3.5

τ

Figure 15: The relative 2 -error E, Egrad , Eboundary and Eboundary flux versus the relaxation time τ at time T = 1 for time-dependent convection-diffusions in an irregular domain. The lattice number N = 64.

5. Conclusions and discussions In this paper, we propose a kind of second-order curved boundary treatments for the lattice Boltzmann method solving two-dimensional convectiondiffusion equations with general nonlinear Robin boundary conditions. The idea is to determine computational boundaries and derive approximate boundary information, with second-order accuracy, by using the prescribed boundary condition. The computational boundaries are zigzag and the approximate boundary information is obtained only at the middle points of two adjacent lattice nodes with one of them in the computational domain another not (see Fig. 1). Once the approximate information are known, the second-order bounce-back schemes can be perfectly adopted. The boundary treatments are validated with a number of examples. The numerical results confirm our theoretical predictions on the second-order accuracy of the treatments. In addition, we find that the classical bounce-back schemes are not reliable at least in solving Neumann boundary problems. In the end, we would like to point out that our idea to treat curved boundary conditions is quite universal. It is independent of the two-dimensional convection-diffusion equations. In particular, the treatment in Section 3.1 can be easily extended to the Navier-Stokes equations with non-slip boundary conditions [16, 17]. The accuracy and stability of this treatment are under our current investigation. For the Navier-Stokes equations with other 28

types of boundary conditions involving normal derivatives of velocities or pressures [23, 24, 25, 26], the calculations in Section 3.2 are promising and can be simplified when the boundary conditions as artificial ones are imposed at outflow and thus on straight boundaries. Furthermore, we are also considering to apply our idea to moving boundary problems of the Navier-Stokes equations [27, 28, 29] or convection-diffusion equations [30, 31]. In this case, the location of the computational boundary (see Fig. 1) may vary with time. But the cost of computing the coordinate of foot points on the boundary (P5 in Fig. 1) may not increase too much if we use the approximation formula (9). We expect to report our progresses on these important issues, together with comparisons with the immersed boundary-lattice Boltzmann method [32, 33, 34, 35, 36], in forthcoming publications. Acknowledgements The authors would like to thank Mr. Yue Gao (Tsinghua University) for valuable discussions. They also acknowledge the Tsinghua National Laboratory for Information Science and Technology for providing computing time. This work was financially supported by the National Natural Science Foundation of China (NSFC 11471185). References [1] U. Frisch, B. Hasslacher, Y. Pomeau, Lattice-gas automata for the Navier-Stokes equation, Physical review letters 56 (14) (1986) 1505. [2] S. Wolfram, Cellular automaton fluids 1: Basic theory, Journal of Statistical Physics 45 (3) (1986) 471–526. [3] X. He, L.-S. Luo, A priori derivation of the lattice Boltzmann equation, Physical Review E 55 (6) (1997) R6333. [4] X. He, L.-S. Luo, Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation, Physical Review E 56 (6) (1997) 6811. [5] S. Chen, G. D. Doolen, Lattice Boltzmann method for fluid flows, Annual Review of Fluid Mechanics 30 (1) (1998) 329–364.

29

[6] S. Succi, I. V. Karlin, H. Chen, Colloquium: Role of the H theorem in lattice Boltzmann hydrodynamic simulations, Reviews of Modern Physics 74 (4) (2002) 1203. [7] D. Yu, R. Mei, L.-S. Luo, W. Shyy, Viscous flow computations with the method of lattice Boltzmann equation, Progress in Aerospace Sciences 39 (5) (2003) 329–367. [8] T. Zhang, B. Shi, Z. Guo, Z. Chai, J. Lu, General bounce-back scheme for concentration boundary condition in the lattice-Boltzmann method, Physical Review E 85 (1) (2012) 016701. [9] Q. Chen, X. Zhang, J. Zhang, Improved treatments for general boundary conditions in the lattice Boltzmann method for convection-diffusion and heat transfer processes, Physical Review E 88 (3) (2013) 033304. [10] L. Li, R. Mei, J. F. Klausner, Boundary conditions for thermal lattice Boltzmann equation method, Journal of Computational Physics 237 (2013) 366–395. [11] Q. Kang, D. Zhang, S. Chen, X. He, Lattice Boltzmann simulation of chemical dissolution in porous media, Physical Review E 65 (3) (2002) 036318. [12] Q. Kang, P. C. Lichtner, D. Zhang, An improved lattice Boltzmann model for multicomponent reactive transport in porous media at the pore scale, Water Resources Research 43 (12) (2007) W12S14–. [13] L. Chen, Q. Kang, B. A. Robinson, Y.-L. He, W.-Q. Tao, Pore-scale modeling of multiphase reactive transport with phase transitions and dissolution-precipitation processes in closed systems, Physical Review E 87 (4) (2013) 043306. [14] I. Ginzburg, Generic boundary conditions for lattice Boltzmann models and their application to advection and anisotropic dispersion equations, Advances in Water Resources 28 (11) (2005) 1196–1216. [15] T. Geb¨ack, A. Heintz, A lattice Boltzmann method for the advectiondiffusion equation with Neumann boundary conditions, Communications in Computational Physics 15 (2) (2014) 487–505.

30

[16] A. J. Ladd, Numerical simulations of particulate suspensions via a discretized Boltzmann equation. part 1. theoretical foundation, Journal of Fluid Mechanics 271 (1994) 285–309. [17] A. J. Ladd, Numerical simulations of particulate suspensions via a discretized Boltzmann equation. part 2. numerical results, Journal of Fluid Mechanics 271 (1994) 311–339. [18] H. Yoshida, M. Nagaoka, Multiple-relaxation-time lattice Boltzmann model for the convection and anisotropic diffusion equation, Journal of Computational Physics 229 (20) (2010) 7774–7795. [19] J. Huang, W.-A. Yong, Boundary conditions of the lattice Boltzmann method for convection–diffusion equations, Journal of Computational Physics 300 (2015) 70–91. [20] B. Shi, Z. Guo, Lattice Boltzmann model for nonlinear convectiondiffusion equations, Physical Review E 79 (1) (2009) 016701. [21] M. P. Do Carmo, Differential geometry of curves and surfaces, Vol. 2, Prentice-hall Englewood Cliffs, 1976. [22] H. Brenner, The diffusion model of longitudinal mixing in beds of finite length. numerical values, Chemical Engineering Science 17 (4) (1962) 229–243. [23] I. Ginzburg, K. Steiner, Lattice boltzmann model for free-surface flow and its application to filling process in casting, Journal of Computational Physics 185 (1) (2003) 61–99. [24] M. Junk, Z. Yang, Outflow boundary conditions for the lattice Boltzmann method, Progress in Computational Fluid Dynamics, an International Journal 8 (1-4) (2008) 38–48. [25] M. Junk, Z. Yang, Pressure boundary condition for the lattice Boltzmann method, Computers & Mathematics with Applications 58 (5) (2009) 922–929. [26] Z. Yang, Lattice Boltzmann outflow treatments: Convective conditions and others, Computers & Mathematics with Applications 65 (2) (2013) 160–171. 31

[27] M. Bouzidi, M. Firdaouss, P. Lallemand, Momentum transfer of a Boltzmann-lattice fluid with boundaries, Physics of Fluids (1994present) 13 (11) (2001) 3452–3459. [28] P. Lallemand, L.-S. Luo, Lattice Boltzmann method for moving boundaries, Journal of Computational Physics 184 (2) (2003) 406–421. [29] P.-H. Kao, R.-J. Yang, An investigation into curved and moving boundary treatments in the lattice Boltzmann method, Journal of Computational Physics 227 (11) (2008) 5671–5690. [30] E. Attar, C. K¨orner, Lattice Boltzmann model for thermal free surface flows with liquid–solid phase transition, International Journal of Heat and Fluid Flow 32 (1) (2011) 156–163. [31] M. Markl, C. K¨orner, Free surface Neumann boundary condition for the advection–diffusion lattice Boltzmann method, Journal of Computational Physics 301 (2015) 230–246. [32] Z.-G. Feng, E. E. Michaelides, The immersed boundary-lattice Boltzmann method for solving fluid–particles interaction problems, Journal of Computational Physics 195 (2) (2004) 602–628. [33] X. Niu, C. Shu, Y. Chew, Y. Peng, A momentum exchange-based immersed boundary-lattice Boltzmann method for simulating incompressible viscous flows, Physics Letters A 354 (3) (2006) 173–182. [34] J. Zhang, P. C. Johnson, A. S. Popel, An immersed boundary lattice Boltzmann approach to simulate deformable liquid capsules and its application to microscopic blood flows, Physical biology 4 (4) (2007) 285. [35] C. Shu, N. Liu, Y.-T. Chew, A novel immersed boundary velocity correction–lattice Boltzmann method and its application to simulate flow past a circular cylinder, Journal of Computational Physics 226 (2) (2007) 1607–1622. [36] J. Wu, C. Shu, Implicit velocity correction-based immersed boundarylattice Boltzmann method and its applications, Journal of Computational Physics 228 (6) (2009) 1963–1979.

32