Immersed boundary lattice Boltzmann model ... - APS link manager

PHYSICAL REVIEW E 85, 016711 (2012)

Immersed boundary lattice Boltzmann model based on multiple relaxation times Jianhua Lu*,† School of Naval Architecture and Ocean Engineering, Huazhong University of Science and Technology, Wuhan 430074, People’s Republic of China

Haifeng Han State Key Laboratory of Coal Combustion, Huazhong University of Science and Technology, Wuhan 430074, People’s Republic of China

Baochang Shi‡ School of Mathematics and Statistics, Huazhong University of Science and Technology, Wuhan 430074, People’s Republic of China

Zhaoli Guo State Key Laboratory of Coal Combustion, Huazhong University of Science and Technology, Wuhan 430074, People’s Republic of China and State Key Laboratory of Enhanced Oil Recovery, Research Institute of Petroleum Exploration and Development, Beijing 100083, People’s Republic of China (Received 21 January 2011; published 27 January 2012) As an alterative version of the lattice Boltzmann models, the multiple relaxation time (MRT) lattice Boltzmann model introduces much less numerical boundary slip than the single relaxation time (SRT) lattice Boltzmann model if some special relationship between the relaxation time parameters is chosen. On the other hand, most current versions of the immersed boundary lattice Boltzmann method, which was first introduced by Feng and improved by many other authors, suffer from numerical boundary slip as has been investigated by Le and Zhang. To reduce such a numerical boundary slip, an immerse boundary lattice Boltzmann model based on multiple relaxation times is proposed in this paper. A special formula is given between two relaxation time parameters in the model. A rigorous analysis and the numerical experiments carried out show that the numerical boundary slip reduces dramatically by using the present model compared to the single-relaxation-time-based model. DOI: 10.1103/PhysRevE.85.016711

PACS number(s): 02.70.−c, 02.60.Lj, 47.11.−j

I. INTRODUCTION

Owing to its simplicity and efficiency, a relatively newly developed numerical scheme, the lattice Boltzmann method (LBM) has attained wide popularity in simulating fluid flow, especially with complex boundaries [1]. In the present lattice Boltzmann (LB) models, two main categories of LB models, the single relaxation time (SRT) and the multiple relaxation time (MRT) versions are the most widely used. As an advanced version of the SRT models, the MRT models possess more flexibility in choosing the relaxation times than the SRT ones [2,3]. As a result, the MRT models are superior to the SRT models in many aspects. Among these aspects, the numerical boundary slip may be one of the most important issues when we simulate fluid flow with complex boundaries such as pore scale fluid flow in porous media. If the SRT model is used, the model will inevitably suffer from the numerical boundary slip unless the relaxation time is chosen as a fixed value [4]; the numerical boundary slip will then introduce an unphysical phenomenon, the viscosity dependence of permeability. However, when the MRT model is employed [2,5], providing that a special combination for two relaxation time parameters is fixed [6,7] as a certain value, the unphysical dependence can be completely

* Present address: Department of Naval Architecture, Dalian University of Technology, Dalian City, 116024, People’s Republic of China. † [email protected] ‡ Corresponding author: [email protected]

1539-3755/2012/85(1)/016711(12)

eliminated. The works by Guo et al. [8] and Chai et al. [9] also showed the advantage of the MRT models compared with the SRT models in such simulations. As for the fluid flow problems with complex boundaries, those with moving boundaries are among the most important ones. Thus in the field of computational fluid dynamics, one of the main issues is handling moving boundaries. Among the existing methods of attacking the problem, the immersed boundary method (IBM), which was first developed by Peskin [10], may be the most popular. In this method, the fluid system is simulated on a fixed Eulerian grid and the moving boundaries are represented by a list of Lagrangian points. The basic idea is to introduce an extra force term acting on the fluid around the boundaries using the Delta function to represent the boundaries’ effect where the boundaries look “immersed” or “embedded.” The primary advantage is that the grid for fluid flow can always be the simple Cartesian one and should not be changed, regardless of the shapes of the moving boundaries. Due to such advantages compared to the body conforming methods, the proposed IBM and its variants have now been widely used in simulating the fluid-structure interaction systems. For more information see Refs. [11–13]. Owing to many common features, such as both using the Cartesian grid, the LBM and IBM were naturally integrated as a new method, which is referred to as the immersed boundary lattice Boltzmann method (IB-LBM). Several authors [14–25] devoted their efforts to developing the IBLBMs and improving their efficiency and accuracy since Feng and Michaelides [14] first combined the two methods.

016711-1

©2012 American Physical Society

JIANHUA LU, HAIFENG HAN, BAOCHANG SHI, AND ZHAOLI GUO

Generally speaking, the existing IB-LBM models differ in their approaches to calculate the required boundary force to enhance the nonslip boundary condition. For example, Feng developed both feedback forcing IB-LBM [14] and direct forcing IB-LBM [15]. In Ref. [14], the boundary force was calculated by Hooke’s law and the force was the function of the boundary surface deformation. In Ref. [15], the force was obtained using the Navier-Stockes equation. In Ref. [16], Niu used the moment exchange approach to get the boundary force. In Ref. [17], Dupuis proposed two other direct forcing approaches. As shown by Peng [21], the IB-LBM method is easier to implement compared to the reflection-type schemes, such as the interpolated bounce-back method. However, the method was found to have less accuracy than the usual boundary schemes. As Le and Zhang [18] pointed out, most existing IB-LBMs suffer from numerical boundary slip. Le and Zhang [18] indicated that the numerical boundary slip can be large enough to destroy the simulation’s validity and usefulness. In their work, a symmetric shear flow was simulated by the IB-LBM proposed by Dupuis [17] as an example. Both the numerical results and analytical solutions showed that the numerical boundary slip is strongly dependent on the relaxation time parameter. In their paper, an open problem was left to reduce the numerical boundary slip. As shown above, by simply setting a special relationship between the relaxation time parameters, the MRT model can be proved to be very effective to eliminate the similar numerical slip for pore scale fluid flow in porous media and then achieve the viscosity-permeability independence. Since the model can be used to deal with such a numerical boundary slip and then achieve certain accurate boundary conditions, it is possible to reduce the numerical slip in the IB-LBM for moving boundary problems by using a similar idea. Motivated by the possibility,


an MRT-based IB-LBM is proposed. Some numerical tests and a rigorous analysis for a simple flow are then carried out to validate the proposed method. It is found that the numerical boundary slip can be dramatically reduced by using the present method compared to the SRT-based IB-LBM.

II. CONVENTIONAL IMMERSED BOUNDARY LATTICE BOLTZMANN METHOD

For simplicity, a direct forcing version of IB-lBM proposed by Dupuis [17] is presented here. For the viscous incompressible flows in a two-dimensional domain containing moving boundaries, the governing equations for the flows can be described as follows: ∇ · u = 0, (1) ∂u (2) + (u · ∇)u = −∇p/ρ + ν∇ 2 u + g, ∂t where u is the velocity, p the pressure, and ρ the density of the fluid. The forcing term g accounts for the nonslip boundary condition and it is introduced to enable the use of the immersed boundary technique. To solve the above governing equations, the standard D2Q9 model is applied fi (x + ci δt,t + δt) − fi (x,t) ωi ρ 1 (eq) = − [fi (x,t) − fi (x,t)] + δt 2 ci · g, τ cs

where fi is the single particle density distribution function in the i direction, δt and ci δt are the time and space increments, respectively, and τ is the relaxation time due to collision for the velocity field. The discrete velocities ci ’s are defined as

⎧ ⎪ ⎨(0,0), ci = (cos[(i − 1)π/2], sin[(i − 1)π/2])), ⎪ √ ⎩ (cos[(2i − 9)π/4], sin[(2i − 9)π/4]) 2c,

(eq)

velocity c is δx/δt. fi is the local equilibrium distribution function as ci · u (ci · u)2 u2 (eq) fi = ωi ρ 1 + 2 + − cs 2cs4 2cs2 with the weight coefficient ⎧ i = 0, ⎨4/9, i = 1–4 ωi = 1/9, ⎩1/36, i = 5–8, √ and cs = c/ 3. The relationship between the density ρ, fluid velocity u, and fi can be written as ρ=

8 i=0

fi , ρu =

8

ci fi .

(3)

i = 0, i = 1,2,3,4, i = 5,6,7,8,

tion of the low Mach number, where ν = (τ − 1/2)cs2 δt.

(4)

To evaluate the force term g, Dupuis [17] followed the technique originally proposed by Mohd-Yusof [13]. If u∗ and ud denote the could-be boundary velocity without the boundary force and the desired boundary velocity with boundary force g, the governing equation (2) can be rewritten as u∗ (x,t + δt) − u(x,t) ∂u = ∂t δt = −(u · ∇)u − ∇p/ρ + ν∇ 2 u =

(5)

without the force term g. Including the force term leads to

i=0

From Eq. (3), we can obtain the governing equations (1) and (2) by the Chapman-Enskog expansion under the assump016711-2

ud (x,t + δt) − u(x,t) = + g(x,t). δt

(6)

IMMERSED BOUNDARY LATTICE BOLTZMANN MODEL . . .


By subtracting Eq. (5) from Eq. (6), we can get an expression of the boundary force as follows:

function is used. First, the could-be velocity U ∗b on a certain boundary point is determined by U ∗ (x b ) = D(x f − x b )u∗ (x f ), (8)

ud (x,t + δt) − u∗ (x,t + δt) . (7) δt To enable the IB-LBM model, the last step is to enable the communication of the body force between the Eulerian grids and the surface of the boundary. To achieve this, the Delta g(x,t) =

f

where the u∗ (x f )’s are the could-be velocity of the fluid points around the boundary point, D(x) is given as

1 + cos πx 1 + cos πy , |x| d and d d D(x) = 0 otherwise, 1 4d 2

|y| d,

(9)

where d is one-half of the immersed boundary layer thickness and has often been given as d = 2δx in the typical IBM simulations [18]. From Eq. (9), we can easily get that ⎧1 x = y = 0, 2, ⎪ ⎪ ⎪ d1 ⎨ , x = ±δx, y = 0 or x = 0,y = ±δx, 2 D(x) = 2d1 (10) ⎪ , x = y = ±δx or x = −y = ±δx, ⎪ 4d 2 ⎪ ⎩ 0, x = mδx, y = nδx if |m| 2 or |n| 2, if d = 2δx.

Similarly, the body force g(x f ) of the fluid point distributed from the nearby boundary point can be written as g(x f ) = D(x f − x b )g(x b ). (11) b

where |fi (x,t) = [f0 (x,t),f1 (x,t),f2 (x,t),f3 (x,t),f4 (x,t),

f5 (x,t),f6 (x,t),f7 (x,t),f8 (x,t)] ,

(13)

As shown by Le and Zhang [18], although the version of IB-LBM was validated by the simulation of flow past an impulsively started cylinder at a low and moderate Reynolds number [17], the model still suffers heavily from numerical boundary slip even for the simple symmetric Couette flow. It was concluded that a reasonable accuracy can be achieved only when the relaxation time τ is less than 2 for the system they studied.

and denotes the transpose operator, M is a matrix that projects f onto the moment space, S is a diagonal relaxation matrix S = diag(s0 ,s1 ,s2 ,s3 ,s4 ,s5 ,s6 ,s7 ,s8 )

(14)

and is here defined by S = diag(0,se ,s ,0,sq ,0,sq ,sν ,sν ),

(15)

| i (x,t) represents the moment vector

III. MULTIPLE-RELAXATION-TIME-BASED IMMERSED BOUNDARY LATTICE BOLTZMANN METHOD

As was shown in Ref. [18], the analytical solution of symmetric Couette flow by the IB-LBM version is strongly dependent on the relaxation time τ and then the numerical boundary slip gets born. Thus it is possible to reduce the numerical boundary slip in IB-LBM models if the dependence is destroyed. As we know, unlike the SRT-LB model on which most of the IB-IBM versions are based, the MRT-LB model has the maximum number of adjustable parameters allowed by the discrete velocity set [3]. So the MRT-based IB-LBM may make some difference. In this section, a MRT-based IB-LBM model is presented to this end. In this model, the standard D2Q9 MRT scheme proposed in Ref. [3] with the extra force model developed by Guo [8] is adapted for fluid flow |fi (x + ci δt,t + δt) − |fi (x,t) (eq)

= −M −1 S | i (x,t) − i (x,t) + (I − S/2)M|gi (x,t)δt},

(12)

| i (x,t) = (ρ,e,ε,jx ,qx ,jy ,qy ,pxx ,pxy ) , (eq)

the equilibrium value of which, | i (x,t), is written as (eq) (x,t) = ρ[1,−2 + 3(u2 + v 2 ),1 − 3(u2 + v 2 ), i

× u,−u,v,−v,u2 − v 2 ,uv] ,

(16)

where u = (u,v). There obviously exists a transformation matrix M between | i (x,t) and |fi (x,t) such that | i (x,t) = M|fi (x,t), (eq) (eq) | i (x,t) = M|fi (x,t), where ⎫ ⎧ 1 1 1 1 1 1 1 1 1⎪ ⎪ ⎪ ⎪ ⎪ ⎪ −4 −1 −1 −1 −1 2 ⎪ 2 2 2⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 4 −2 −2 −2 −2 1 ⎪ ⎪ 1 1 1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 0 1 0 −1 0 1 −1 −1 1⎬ ⎨ 0 −2 0 2 0 1 −1 −1 1 . M= ⎪ ⎪ ⎪ 0 0 1 0 −1 1 1 −1 −1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 0 0 −2 0 2 1 1 −1 −1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 0 ⎪ 1 −1 1 −1 0 0 0 0⎪ ⎪ ⎪ ⎭ ⎩ 0 0 0 0 0 1 −1 1 −1

016711-3


In Eq. (12), I is the unity matrix, |gi (x,t) represents the force term |gi (x,t) = (g0 ,g1 ,g2 ,g3 ,g4 ,g5 ,g6 ,g7 ,g8 ) with

gi = ωi ρ

(ug + gu)/2 : ci ci − cs2 I ci · g + 2 cs cs4

.

8 − sv . (20) 8(2 − sv ) As a result, the viscosity-permeability independence can be achieved when the pore scale simulation for fluid flow in porous media is performed by the MRT model plus halfway bounce back at the solid walls with such a relationship between sv and sq . Similarly, to reduce the numerical boundary slip on the immersed boundaries, sv and sq in S should satisfy a certain relationship if the MRT-based IBM is adapted. Through theoretical analysis for simple flow cases with moving boundaries, if the following formula holds, sq =

4(2 − sv ) , (21) 4 + 7sv it is possible to reduce the numerical slip dramatically. Similarly, if we substitute Eq. (21) for Eq. (19), here also keeps a constant 98 . Besides, no special restrictions need to be imposed on se and s . A detailed analysis will be presented in Sec.V. For the evaluation of the boundary force, the same procedure is applied as it was done in Ref. [18]. As shown in Sec. II, the procedure is carried out as the following steps. (1) Compute u∗ (x f ) using Eq. (12) without the boundary force term. (2) Compute U ∗ (x b ) using Eq. (8). (3) Compute g(x b ) using Eq. (48) rather than Eq. (7). [As shown by Kang et al. [24], a better accuracy can be obtained if Eq. (48) rather than Eq. (7) is adapted.] (4) Compute g(x f ) using Eq. (11). Unlike the SRT model with extra force term shown in Sec. II by Eq. (3), the relationship between the density ρ, fluid velocity u, and fi here can be written as

i=0

8

1 τ

for all i = 0–8, the MRT

IV. NUMERICAL TESTS

For the halfway bounce back rule, if keeps as a constant 3 value, such as 16 , an exact solution can be obtained for a straight channel, which means that sv and sq satisfy

fi , ρu =

1 . τ

(17)

However, if we want to eliminate or reduce the numerical slip coming from the boundaries, a special combination of certain relaxation time parameters should keep a constant value. As shown by Ginzburg et al. [6], the combination is 1 1 1 1 = . (19) − − sv 2 sq 2

8

sv =

If we set S as Eq. (14) and si = LBM collapses to the SRT LBM.

Similar to the general MRT model, the relaxation time parameter sv here also satisfies 1 1 . (18) − ν = cs2 sv 2

ρ=

By comparing Eq. (4) with Eq. (18), we have

sq =


In this section, four simple flows are simulated to test the efficiency of the present MRT-based IB-LBM model. For comparisons, we perform simulations for all cases by the SRT LBM at the same time. Here the SRT LBM is obtained by setting S as Eq. (14) and si = τ1 for all i = 0–8 in the MRT LBM. In the present simulations, se and s are both set as sv in the MRT version. As for the force scheme, we follow Kang et al. [24] to achieve better accuracy by using Eq. (48) rather than Eq. (7). Unless otherwise stated, all simulations are conducted over a grid of 200 × 200 which represents a computational domain [0.0,1.0] × [0.0,1.0]. The test problems are the two-dimensional symmetric Couette flow, asymmetric Couette flow, cylindrical Couette flow, and symmetric Poiseuille flow. For all the problems the theoretical solutions can be obtained. (a) Symmetric Couette flow: The test problems are set as the same as in Ref. [18]. As show in Fig. 1, two plain surfaces are also placed at y = 50δx and y = 150δx and move with opposite horizontal velocities U d and −U d . (b) Asymmetric Couette flow: The test problems are also set as the same as in Ref. [18]. Two plain surfaces are also placed at y = 50δx and y = 100δx and move with opposite horizontal velocities U d and −U d . (c) Cylindrical Couette flow: The test problems are still set as the same as in Ref. [18]. Two rings with radii of 45δx and 70δx are placed in the center of the computational domain. The outter ring is fixed while the inner one is rotating at the velocity of U d . (d) Symmetric Poiseuille flow: Two plain surfaces are placed at y = 50δx and y = 150δx; the same as test problem (a). The difference is that both the surfaces are still. To drive the fluid, a constant body force denoted by F is imposed on the horizontal direction x. For convenience, the x component of F is denoted by F . As for the boundary condition implementation, the generical periodic boundary conditions are applied in both the x and y directions. A wide range of parameters such as U d and F are investigated only if the simulation is stable.

ci fi + ρ gδt/2.

moving plate u = − U d y x moving plate u = U d

FIG. 1. The layout of symmetric Couette flow.

i=0

016711-4



0.8

1.5 s_v = 0.05 0.10 0.20 0.50 1.00 1.50

0.6 0.4 0.2

u −uT Ud 0

1

0.7 0.5

−0.2

0

−0.5

−0.4 −1

−0.6 0.2

0.2 0.2

u Ud

−0.8 0

1.2

0.4

y

0.6

0.8

−1.5 0

1

0.2

0.4

y

(a) 1.5 s_v = 0.05 0.10 0.20 0.50 1.00 1.50

0.005 u −uT Ud

0.9 0.23 0

−0.5

−0.01

−1

y

1

0.6

0.8

1

u

−0.005

0.4

0.8

0.5 Ud

0.2

0.6

1.1

1

0

−0.015 0

0.3

(b)

0.015 0.01

0.25

theory s_v = 0.05 0.10 0.20 0.50 1.00 1.50

1

−1.5 0

(c)

0.25

0.27

theory s_v = 0.05 0.10 0.20 0.50 1.00 1.50 0.2

0.4

y

0.6

0.8

1

(d)

T FIG. 2. (Color online) Profiles of the horizontal velocity and its errors (denoted by Uud and u−u ) normalized with respect to the theoretical Ud maximum velocity U d along the central vertical plain x = 100δx of the symmetric Couette flow under a different relaxation time parameter sv . Here uT represents the theoretical solution, U d is the maximum velocity for the theoretical solution. sv = 0.05,0.10,0.20,0.50,1.00,1.50. (a), (b) The results of the SRT IB-LBM model are given. (c), (d) The results of the present model are shown.

T The profiles of the horizontal velocity Uud and its errors u−u Ud along the central vertical plain x = 100δx of test problem (a) under different relaxation time parameters are plotted in Fig. 2. As shown in Figs. 2(a) and 2(b), the velocities computed from the SRT version deviate from the theoretical results greatly in almost all cases. As the relaxation time parameter sv decreases (τ increases), the deviation gets larger and the computed results then cannot be acceptable. The phenomena show a perfect agreement with those shown in Ref. [18]. However, the computed velocities by the present MRT-based IB-LBM make a difference indeed [Figs. 2(c) and 2(d)]. In fact, they show a better performance in all cases. The numerical boundary slip is now reduced to a rather low magnitude. Generally, the slips can be accepted whatever value sv is given, unlike those in the SRT version. However, a careful observation focused on the immersed boundary layer [Fig. 2(d)] shows that the numerical boundary slip is not eliminated completely. When sv is larger, the slip becomes more apparent. As a more general situation than test problem (a), we can find similar results for test problem (b). Figure 3 shows us T the profiles of the horizontal velocity Uud and its errors u−u Ud

along the central vertical plain x = 100δx of test problem (b) under different relaxation time parameters. As can be seen, the velocity profile obtained by SRT-based IB-LBM is strongly dependent on the the choice of the relaxation time parameter. For most of the sv tested, the model cannot give satisfactory results. However, when the present MRT-based IB-LBM is adapted, the dependence is almost invisible. Similar to that of test problem (a), the obtained velocity profiles with different values of the relaxation time parameter almost overlap with each other. However, the numerical slip still exists. The larger sv also corresponds to the larger slip. In test problems (a) and (b), the moving boundaries are straight lines. To test the model for more complex boundaries, we conduct the simulation on test problem (c). In the problem, the curved immersed boundaries are included. As a result, not only do the linear velocity profiles appear, but the nonlinear ones appear also. Figure 4 shows us the profiles of the T horizontal velocities Uud and their relative errors u−u along the Ud central horizontal plain y = 100δx under different relaxation time parameters. Similarly, the velocity profile obtained by SRT-based IB-LBM again shows a strong dependence on the

016711-5



1

1.5 s_v = 0.05 0.10 0.20 0.50 1.00 1.50

0.5 u −uT Ud

2 1 0.5

0

0

−0.5 −0.5 −1

0.2

0 0.23

u Ud

−1 0

1

0.4

y

0.6

0.8

−1.5 0

1

0.2

0.4

0.6

y

0.8

1

(b)

0.02

1.5 s_v = 0.05 0.10 0.20 0.50 1.00 1.50

0.015 0.01 0.005

u −uT

1

1

0.5

0.8

Ud

−0.005

0

−0.5

−0.01 −1

−0.015 0.2

0.4

y

0.6

0.8

0.24

0.26

u

0

−0.02 0

0.27

theory s_v = 0.05 0.10 0.20 0.50 1.00 1.50

(a)

Ud

0.25

1

−1.5 0

(c)

theory s_v = 0.05 0.10 0.20 0.50 1.00 1.50 0.2

0.4

y

0.6

0.8

1

(d)

T FIG. 3. (Color online) Profiles of the horizontal velocity and its errors (denoted by Uud and u−u ) normalized with respect to the theoretical Ud d maximum velocity U along the central vertical plain x = 100δx of the asymmetric Couette flow under a different relaxation time parameter sv . Here uT represents the theoretical solution, U d is the maximum velocity for the theoretical solution. sv = 0.05,0.10,0.20,0.50,1.00,1.50. (a), (b) The results of the SRT IB-LBM model are given. (c), (d) The results of the present model are shown.

relaxation time parameter. The numerical results for most values of sv tested are out of usefulness. Similar to those of test problems (a) and (b), when the present MRT-based IB-LBM is used, the obtained velocity profiles with different values of the relaxation time parameter again show almost no difference with each other. The numerical accuracy of the present model seems to be not affected although the curved immersed boundaries are introduced. However, the numerical slip again exists. The larger sv again corresponds to the larger slip. The test problems (a) through (c) are Couette-like ones. To test the model in more general cases, test problem (d), the symmetric Poiseuille flow is chosen. In the problem, only nonlinear velocity profiles are presented. Figure 5 shows T the profiles of the horizontal velocity errors u−u along the Uc central vertical plain x = 100δx under different relaxation time parameters. All values are normalized with respect to the 2 theoretical maximum velocity Uc = L8νF similar to those done in test problem (a). As shown in Fig. 5(a), the results of SRTLBM cannot be accepted in most cases except the cases where sv approaches 1.0. On the contrary, the present model always

gives us satisfactory results compared with those from the SRT-based IB-LBM model. Still the numerical boundary slip is not eliminated totally. Similar to test problems (a) through (c), the slip here is also more apparent when sv gets larger. For all test problems, the present MRT-based IBM model shows much better accuracy than the SRT-based one. However, 3 unlike the MRT model with Eq. (20) (means = 16 ) for flow in porous media, the numerical slip here still cannot be eliminated totally for the above test problems by the present model. As was shown in Ref. [26], if ≈ > 12 , the exact solution cannot be obtained even for Poiseuille flow. As we know, in the present MRT-based IBM, here is 98 . Is the value too high to eliminate the numerical slip for the test problems? To closely investigate this issue, the comparison study for test problems (a), (c), and (d) using a large range of is carried out. The corresponding relative errors (ER ) T| for velocity like |u−u are compared, where uT represents |uT | the analytical solution. The convergence rates for different values are also presented. Figures 6–8 show the relative errors as a function of . The maximum range of is [0.005,6.625]. Two values of sv ,

016711-6



0.4

1.5 s_v = 0.05 0.10 0.20 0.50 1.00 1.50

0.3 0.2 0.1

u −u T Ud

1.2 1 0.5 Ud

−0.1

0

−0.5

−0.2 −1

−0.3 0.2

0.4

x

0.6

0.8

1

−1.5 0

0.2

0.4

x

0.6

0.8

1

(b)

0.03

1.5 s_v = 0.05 0.10 0.20 0.50 1.00 1.50

0.02 0.01 u −u T

1.05 1 0.5 Ud

0

−0.01

−0.5

−0.02

−1

0.2

0.4

x

0.6

1 0.95 0.27

u

0

−0.03 0

0.3

theory s_v = 0.05 0.10 0.20 0.50 1.00 1.50

(a)

Ud

0.2 0.25

u

0

−0.4 0

0.7

0.8

−1.5 0

1

(c)

0.275

0.28

theory s_v = 0.05 0.10 0.20 0.50 1.00 1.50 0.2

0.4

x

0.6

0.8

1

(d)

T FIG. 4. (Color online) Profiles of the horizontal velocity and its errors (denoted by Uud and u−u ) normalized with respect to the theoretical Ud d maximum velocity U along the central vertical plain y = 100δx of the cylindrical Couette flow under a different relaxation time parameter sv . Here uT represents the theoretical solution, U d is the maximum velocity for the theoretical solution. sv = 0.05,0.10,0.20,0.50,1.00,1.50. (a), (b)The results of the SRT IB-LBM model are given. (c), (d) The results of the present model are shown.

0.05,1.00 are tested. The velocities are the horizontal ones of the central vertical plain for test problems (a) and (d) while they are the central horizontal plain for test problem (c). As shown by these figures, unlike the MRT model for porous flow, 12 here seems not to achieve better results. On the contrary, less accuracy loss appears for 1.0. At the same time, the errors seems always to decrease with and then increase with the parameter in the considered value range. As expected, = 98 in general obtains relative small errors for different values of sv . However, the minimum errors for all these test problems do not correspond to = 98 . In fact, the optimal values of for minimum errors always seems larger than 98 = 1.125. For test problem (c). which includes the most complex geometry, the optimal value is even large than 2.625. For a larger value of sv , the optimal value seems larger. By simple algebra, we can predict the optimal value theoretically for test problem (a), that is, = 1.47575,1.88564 for sv = 0.05,1.00. Figure 6 shows that the minimum errors are indeed obtained at these values. Figures 9–11 present the convergence rates for the three test problems under different values of , 0.005,1.125,2.625. Also two values of sv , 0.05,1.00 are tested. Four grid sizes 100 × 100, 200 × 200, 400 × 400, and 800 × 800 are used.

All these figures show that the approximate linear convergence rates can be achieved. In general, relatively smaller errors can be achieved for = 98 no matter what value sv takes. For smaller sv , the errors seem smaller also. As the results have shown in Fig. 7, Fig. 10 demonstrates that, whatever the grid size takes, = 2.625 gives the smallest errors for cylindrical Couette flow in which the curved immersed boundaries are introduced. In summary, the numerical results show that the numerical boundary slip can be reduced dramatically by using the present MRT-based IBM. However, the slip cannot be eliminated completely whatever value takes. = 1.125 always gives rather less accuracy loss. It seems that there exists an optimal value for if sv fixed, under the condition the minimum error can be achieved. In general, the minimum values seem always larger than 1.125. V. THEORETICAL ANALYSIS

To clarify how the present MRT-based IBM reduces the numerical boundary slip, we perform a rigorous analysis on the present model for the symmetric Couette flow following a

016711-7



8

0.005 s_v = 0.05 0.10 0.20 0.50 1.00 1.50

6

u −uT Uc

4

0

u −uT Uc

2

−0.005 −0.01

0

−0.015

−2

−0.02

−4 0

0.2

0.4

y

0.6

0.8

−0.025 0

1

s_v = 0.05 0.10 0.20 0.50 1.00 1.50 0.2

0.4

(a)

0.6

y

0.8

1

(b)

T FIG. 5. (Color online) Profiles of the horizontal velocity errors (denoted by u−u ) normalized with respect to the theoretical maximum Uc velocity Uc along the central vertical plain x = 100δx of the symmetric Poiseuille flow under a different relaxation time parameter sv . Here uT represents the theoretical solution, Uc is the maximum velocity for the theoretical solution. sv = 0.05,0.10,0.20,0.50,1.00,1.50. (a) The results of the SRT IB-LBM model are given. (b) The results of the present model are shown.

similar process to He et al. [4] and Le and Zhang [18]. For the simple two-dimensional, steady, and incompressible Couette flow in the x direction, the following equations are satisfied ρ = const, v = 0,

∂u = 0. ∂t

∂u = 0, ∂x

where

⎧ ⎨j, j = j + 1, ⎩j − 1.

(22)

If we denote the y position of the lattice node by the superscripts of mass distribution function fi , extra force term gi , and subscripts of x velocity u and force G, according to Eq. (12), at a steady state, we get j j

f¯ = f − M −1 SM f j − f j,(eq) i i i i j + (I − S/2)M gi δt , (23) j

where |f¯i represents the postcollision distribution function and satisfies j f (t + δt) = f¯j (t) , (24) i i

i = 0,1,3, i = 2,5,6, i = 4,7,8.

In the meantime, the momentum density in the x direction is j j

j Gj j j j ρuj = c f1 − f3 − f5 − f6 + f8 − f7 + δtρ . 2 (25) Thus, using Eqs. (23)–(25) and considering the momentum exchange due to applied force at a lattice node, we can prove that 1 + sv 2 − sv 3sv − 8 sv ρuj = ρuj + ρ(uj −1 + uj +1 ) + 3 6 12 × δtρ(Gj −1 + Gj +1 − 2Gj ) + sv δtρGj (26) by using sq =

0.02

2(2−sv ) , 4sv +2−sv

which is derived from Eq. (19).

0.04

s_v = 0.05 1.00

s_v = 0.05 1.00

0.035

0.015

ER

ER

0.03

0.01

0.025 0.02

0.005 0.015

0 0

0.5

1

1.5

2

2.5

0.01

3

0

2

4

6

8

Λ

Λ

FIG. 6. (Color online) Relative error of horizontal velocity u along the central vertical plain as a function of for symmetric Couette flow.

FIG. 7. (Color online) Relative error of horizontal velocity u along the central horizontal plain as a function of for cylindrical Couette flow.

016711-8



0.05

−3.5

s_v = 0.05 1.00

0.04

−4.5

ER

log2 ER

0.03

0.02

−5.5 −6.5 −7.5

0.01

−8.5 0

0

0.5

1

1.5

2

2.5

6.5

3

7

7.5

FIG. 8. (Color online) Relative error of horizontal velocity u along the central vertical plain as a function of for symmetric Poiseuille flow.

By a simple algebraic manipulation of Eq. (26) we can get 3 − 8 (Gj −1 + Gj +1 − 2Gj ) 0 = Gj + 12 uj −1 + uj +1 − 2uj , +ν δx 2

8.5

9

9.5

FIG. 10. (Color online) Comparison of convergence rates under different and sv for cylindrical Couette flow. The solid, dashed, and dotted lines correspond to = 1.125, 0.005, 2.625 while lines with asterisks and squares are for sv = 0.05,1.00.

know that

⎧ ⎨GB /2, Gj = GB /4, ⎩0,

(27)

which can be considered as a central finite different form of the simplified Navier-Stokes equation 0 = G(y) + ν

8

-log2 δx

Λ

∂ 2u . ∂y 2

(28)

j = j0 , j = j0 ± 1, |j − j0 | > 1,

(29)

where GB is the total force at the immersed boundary node. Owing to both the forcing symmetry and geometry symmetry of the present system, we have uj0 −1 = uj0 +1 . By substituting Eq. (29) into Eq. (27) we can get the velocity differences between two adjacent lattice layers δx 2 GB 9 + 8 · , 2ν 24 δx 2 GB 21 + 8 uj0 +1 − uj0 +2 = · , 2ν 24 δx 2 GB uj0 +2 − uj0 +3 = . 2ν uj0 − uj0 +1 =

Now let us focus on the velocity profile of the symmetric Couette flow obtained from the MRT-based IB-LBM. We start from the bottom plain at y = j0 δx (here j0 = 50 in our simulations) and move upward. From Eqs. (9) and (11) we

(30) (31) (32)

−4 −3

−5

−4 −5

−7

log2 ER

log2 ER

−6

−8

−6 −7

−9

−8

−10

−9

−11 6.5

7

7.5

8

8.5

9

−10 6.5

9.5

-log2 δx

7

7.5

8

8.5

9

9.5

- log2 δx

FIG. 9. (Color online) Comparison of convergence rates under different and sv for symmetric Couette flow. The solid, dashed, and dotted lines correspond to = 1.125, 0.005, and 2.625 while lines with asterisks and squares are for sv = 0.05,1.00.

FIG. 11. (Color online) Comparison of convergence rates under different and sv for symmetric Poiseuille flow. The solid, dashed, and dotted lines correspond to = 1.125, 0.005, and 2.625 while lines with asterisks and squares are for sv = 0.05,1.00.

016711-9


Beyond the immersed boundary layer, the boundary force is 0. By substituting the forces to Eq. (27), the velocity gradient is found to remain constant δx 2 GB , for j > j0 + 2. uj −1 − uj = (33) 2ν Combining Eqs. (30) and (31), we have δx 2 GB 15 + 8 · , (34) 2ν 12 whereas the theoretical prediction velocity gradient (here the assumption of the linear velocity profile for Couette flow is partly used) should be


thus u∗j can be written as u∗j =

uj0 − uj0 +2 =

2

δx GB (35) ν then the difference between Eqs. (34) and (35) can be written as δx 2 GB −9 + 8 · . (36) usj0 = uj0 − u¯ j0 = 2ν 12 Thus, to nullify the difference, or eliminate the numerical velocity slip from the MRT-based IB-LBM, we only need to impose usj0 ≡ 0, which yields Eq. (21) and leads to the present MRT-based IB-LBM. Assume that Eq. (21) is satisfied. Equations (30) and (31) can be rewritten as 3δx 2 GB uj0 − uj0 +1 = , (37) 8ν 5δx 2 GB . (38) uj0 +1 − uj0 +2 = 8ν Since the velocities at the center plain between the two moving plains are always zero, and the velocity gradients are known everywhere, the velocities on the horizontal plains j0 and j0 ± 1 can be easily deduced u¯ j0 − uj0 +2 =

uj0 ±1

Using Eq. (44) with Eq. (8) we can get the could-be velocity at the boundary node UB∗ = (u∗j0 + u∗j1 )/2, that is, UB∗ =

UB∗ =

If Eq. (7) has the form g(x,t) = 2

Using Eq. (42) and j ∗ f (t + δt) = f¯j (t) , i i

h 8 − 242sv − 367sv2 δx 2 GB · + U = . 2ν 2 144sv (4 + 7sv )

(43)

where the symbol is the mass distribution function corresponding to the could-be velocity u∗j , and

i=0

j∗

8 i=0

j∗

ci fi

(48)

(49)

As shown by Kang et al. [24], a better accuracy can be obtained if Eq. (48) rather than Eq. (7) is adapted. In the present work we follow Kang’s choice. With the relations as Eq. (49), we can deduce the analytical , computed boundary solutions for the bulk velocity gradient du dy velocity uj0 , and numerical boundary slip usj0 as follows: du dy bulk 2U d hδx

=

uj0 = Ud usj0 Ud

∗

fi , ρj∗ u∗j =

ud (x,t + δt) − u∗ (x,t + δt) , δt

then U d is written as

(39) (40)

(46)

then the relationship between the velocity of boundary node U d and the resulting boundary force GB can be obtained from Eqs. (7) and (46) h 200 − 2sv − 535sv2 δx 2 GB · + . (47) Ud = 2ν 2 144sv (4 + 7sv )

Next let us evaluate the boundary force term. As shown in Sec. III, u∗j should be computed first. To get u∗j , we drop the force term from the governing equation (23) j j

f¯ = f − M −1 SM f j − f j,(eq) . (42) i i i i

8

h 184 + 482sv + 199sv2 δx 2 GB · − 2ν 2 144sv (4 + 7sv )

d

where h (here h = 100 in the present simulations) is the gap between two moving plains. Therefore, the boundary velocity calculated from uj0 and uj0 +1 with Eq. (8) is δx 2 GB h 3 − . (41) UB = 2ν 2 8

ρj∗ =

δx 2 GB 2ν h 46 + (89 + 28 )sv − 8(7 − 32 + 18 2 )sv2 . · − 2 72sv [2 + (−1 + 4 )sv ] (45)

If = 98 ,

2

δx GB h · , 2ν 2 δx 2 GB h 3 − , = 2ν 2 4

uj0 =

1 + sv 2 − sv uj + (uj −1 + uj +1 ) 3 6 8 + 2(−5 + 12 )sv + (3 − 20 + 32 2 )sv2 − 12[2 + (−1 + 4 )sv ] 3 − 2sv δtGj . (44) × δt(Gj −1 + Gj +1 − 2Gj ) − 2

=

h 4

+

h 4

+

h 4 8−242sv −367sv2 288sv (4+7sv ) h 4 8−242sv −367sv2 288sv (4+7sv )

+

8−242sv −367sv2 288sv (4+7sv )

0 h 4

,

(50)

,

(51)

.

(52)

Following a similar procedure, the present MRT-based IBLBM can be derived for the symmetric Poiseuille flow. For details please refer to the Appendix. As observed in the numerical experiments presented in Sec. IV, the theoretical results Eqs. (50) and (51) also show that the numerical boundary slip in the present model still exist. The

016711-10


which is derived from Eq. (19), we can get U d as follows:

Ud =

1.008 LBM solution analytical solution fluid mechanics theory

Ud

1.006

uj0

exact solution of Couette flow cannot be obtained. As shown by Eqs. (37) and (38), the linear velocity gradient between the neighborhood grid in the boundary layer cannot be met. However, the total velocity gradient in the immersed boundary layer uj0 − uj0 +2 is linear, thus the leading numerical error is nullified. As a result, the boundary slip can be approximately eliminated when h is large enough and sv is in a proper value range. We can see it in Fig. 12 more clearly. In the figure, the numerical results and analytical solutions are both plotted under different sv for the symmetric Couette flow. From the figure, we can find that whatever sv changes, the numerical boundary slip always exists in the present model, but the slip can be always small enough generally. The figure also proves the correctness of both the analysis and the implementation in the present work. In addition, we can predict the optimal values for which correspond to the minimum numerical slips or errors. v) If we do not adapt Eq. (21) and only use sq = 4s2(2−s , v +2−sv


1.004

1.002

1 0.2

0.4

0.6

0.8

sv

1

1.2

1.4

FIG. 12. (Color online) Boundary slip velocity from the present MRT-based IB-LBM changing with the relaxation time parameter sv for the symmetric Couette flow, where sv changes in the range of [0.05,1.50].

h 2 − 137sv + 68 sv + 68sv2 − 304 sv2 + 144 2 sv2 δx 2 GB · + . 2ν 2 144sv [2 + (−1 + 4 )sv ]

It should be noted that Eq. (53) reduces to Eq. (49) if = 98 . √ −17+76sv + 217+2348sv +3328sv2 In fact, by setting = , Eq. (53) 72sv can also reduce to uj0 . For example, if sv = 0.05,1.00, = 1.47575,1.88564. Thus, the optimal values for are predicted as a result. Figure 6 shows the validation of the prediction. In the end, unlike the MRT model combined with the bounce back scheme with constant for pore scale flow fluid in porous media, loses its control of the numerical error in the present MRT model combined with the immersed boundary scheme for moving boundary problems. As shown in Eq. (53), the error also changes with sv . The control is first lost when u∗j is obtained and then UB∗ begins to show extra dependence on sv expect in the form of ν. Obviously, the main reason why the control is lost lies in the adoption of the procedure obtaining u∗j . So for further improvement of the present model the procedure should be carefully investigated. In addition, the proper choice of D(x) and the thickness of the immersed boundary layer may be also helpful for the method in the control of the numerical error by because the procedure to get u∗j is not the final one.

VI. CONCLUSION

In this work, an MRT-based IB-LBM is proposed to try to reduce the numerical boundary slip. Numerical experiments and rigorous analysis show that the present model has a much better performance than the corresponding SRT version. For the test problems, as the leading numerical error is eliminated, the numerical slip can be acceptable generally as has been shown. However, even for the symmetric Couette flow, the exact solution cannot be obtained and the numerical boundary

(53)

slip still cannot be eliminated completely in the present framework. Further investigations on the value of show that there seems always to exist optimal values of 1.125 which correspond to the minimum numerical slip or errors. The main reason may lie in the procedure in which u∗j is obtained. The more proper choice of D(x) and the thickness of the immersed boundary layer may also be helpful. As we know, the TRT model not only possesses the same advantage to reduce the numerical boundary slip as general MRT, but also has a simpler form [27]. It may be a better choice to develop similar TRT-based IB-LBM. Thus, a more careful investigation on these issues should be done in the future to improve the present model. ACKNOWLEDGMENTS

J. H. Lu acknowledges the financial support from the National Science Fund of China (Grant No. 60773195, 51006039, 51006040), the National Basic Research Program of China (Grant No. 2011CB707305), the Fundamental Research Funds for the Central Universities of China (DUT11RC(3)87) and a start-up grant of the postdoctoral research from Huazhong University of Science and Technology. The work is made possible by the High Performance Computing Center experimental testbed in SCTS/CGCL. Thanks should also be given to Professor J.F. Zhang, Dr. Z.H. Chai, and two anonymous referees for their helpful suggestions. APPENDIX: DERIVATION FOR THE SYMMETRIC POISEUILLE FLOW

We start from Eq. (26). Unlike the derivation for the symmetric Couette flow, the force term Gj has the following

016711-11


expression owing to the constant body force F imposed on the fluid: ⎧ ⎨GB /2 + F, j = j0 , Gj = GB /4 + F, j = j0 ± 1, (A1) ⎩F, |j − j0 | > 1, where F denotes the body force F in the horizontal direction. Following the same procedure, the velocity differences between two adjacent lattice layers can also be obtained δx 2 F δx 2 GB 9 + 8 · +1· , (A2) uj0 − uj0 +1 = 2ν 24 2ν δx 2 F δx 2 GB 21 + 8 uj0 +1 − uj0 +2 = · +3· , (A3) 2ν 24 2ν δx 2 GB δx 2 F uj0 +2 − uj0 +3 = +5· . (A4) 2ν 2ν Beyond the immersed boundary layer, the boundary force is 0 and the velocity gradient is written as follows:


gradient adjacent to the central plain must be met uj0 +h/2−1 − uj0 +h/2 = −(uj0 +h/2 − uj0 +h/2+1 ).

(A6)

With this condition, the relationship between GB and F can be obtained GB = −F h.

(A7)

From the theoretical prediction (here the assumption of the parabolic velocity profile for Poiseuille flow is partly used), we know that δx 2 F · 2(h − 2). (A8) 2ν On the other hand, from Eqs. (A2) and (A3) with Eq. (A7) we can get h(15 + 8 ) δx 2 F uj0 − uj0 +2 = 4− (A9) 2ν 12 u¯ j0 − uj0 +2 = −

then the numerical velocity slip here is

As we know, the velocity at the central plain j = j0 + h/2 always has the maximum value, then the opposite velocity

δx 2 F h 9 − 8 · . (A10) 2ν 12 Setting the numerical boundary slip usj0 = 0 also yields Eq. (21) and leads to the present MRT-based IB-LBM. Similar to the symmetric Couette flow case, the solution here is also not the optimal one and the exact solution cannot be obtained.

[1] C. K. Aidun and J. R. Clausen, Annu. Rev. Fluid Mech. 42, 439 (2010). [2] D. d’Humières, Prog. Astronaut. Aeronaut. 159, 450 (1992). [3] P. Lallemand and L.-S. Luo, Phys. Rev. E 61, 6546 (2000). [4] X. Y. He, Q. S. Zou, L. S. Luo, and M. Dembo, J. Stat. Phys. 87, 115 (1997). [5] C. X. Pan, L. S. Luo, and C. T. Miller, Comput. Fluids 35, 898 (2006). [6] I. Ginzburg and D. d’Humières, Phys. Rev. E 68, 066614 (2003). [7] D. d’Humières and I. Ginzburg, Comp. Math. Appl. 58, 823 (2009). [8] Z. L. Guo and C. G. Zheng, Int. J. Comput. Fluid Dyna. 22, 465 (2008). [9] Z. H. Chai, B. C. Shi, Z. L. Guo, and J. H. Lu, Commun. Comput. Phys. 8, 1052 (2010). [10] C. S. Peskin, J. Comput. Phys. 25, 220 (1977). [11] C. S. Peskin, Acta Numer. 11, 479 (2002). [12] R. Mittal and G. Iaccarino, Annu. Rev. Fluid Mech. 37, 239 (2005). [13] J. Mohd-Yusof, CTR Annu. Res. Briefs, 317 (1997). [14] Z. G. Feng and E. E. Michaelides, J. Comput. Phys. 195, 602 (2004).

[15] Z. G. Feng and E. E. Michaelides, J. Comput. Phys. 202, 20 (2005). [16] X. D. Niu, C. Shu, Y. T. Chew, and Y. Peng, Phys. Lett. A 354, 173 (2006). [17] A. Dupuis, P. Chatelain, and P. Koumoutsakos, J. Comput. Phys. 227, 4486 (2008). [18] G. G. Le and J. F. Zhang, Phys. Rev. E 79, 026701 (2009). [19] C. Shu, N. Y. Liu, and Y. T. Chew, J. Comput. Phys. 226, 1607 (2007). [20] Y. Sui, Y. T. Chew, P. Roy, and H. T. Hong-Tong Low, Int. J. Numer. Meth. Fluids 53, 1727 (2007). [21] Y. Peng and L. S. Luo, Prog. Comput. Fluid Dyna. 8, 156 (2008). [22] J. Wu and C. Shu, J. Comput. Phys. 228, 1963 (2009). [23] X. L. Yang, X. Zhang, Z. L. Li, and G. W. He, J. Comput. Phys. 228, 7821 (2009). [24] S. K. Kang and Y. A. Hassan, Int. J. Numer. Meth. Fluids 66, 1132 (2011). [25] H. K. Jeong, H. S. Yoon, M. Y. Ha, and M. Tsutahara, J. Comput. Phys. 229, 2526 (2010). [26] I. Ginzburg, F. Verhaeghe, and D. d’Humières, Commun. Comp. Phys. 3, 427 (2008). [27] I. Ginzburg, F. Verhaeghe, D. d’Humières, Commun. Comp. Phys. 3, 519 (2008).

δx 2 GB δx 2 F + [2(j − j0 ) − 1)] · , uj −1 − uj = 2ν 2ν for j > j0 + 2. (A5)

usj0 = uj0 − u¯ j0 =

016711-12