Convection-Diffusion Lattice Boltzmann Scheme for Irregular Lattices

Journal of Computational Physics 160, 766–782 (2000) doi:10.1006/jcph.2000.6491, available online at http://www.idealibrary.com on

Convection-Diffusion Lattice Boltzmann Scheme for Irregular Lattices R. G. M. van der Sman∗ and M. H. Ernst† ∗ Agrotechnological Research Institute, P.O. Box 17, Wageningen, NL-6700 AA, The Netherlands; †Institute for Theoretical Physics, University of Utrecht, Utrecht, The Netherlands E-mail: [email protected], [email protected] Received March 26, 1999; revised November 24, 1999

In this paper, a lattice Boltzmann (LB) scheme for convection diffusion on irregular lattices is presented, which is free of any interpolation or coarse graining step. The scheme is derived using the axioma that the velocity moments of the equilibrium distribution equal those of the Maxwell–Boltzmann distribution. The axioma holds for both Bravais and irregular lattices, implying a single framework for LB schemes for all lattice types. By solving benchmark problems we have shown that the scheme is indeed consistent with convection diffusion. Furthermore, we have compared the performance of the LB schemes with that of finite difference and finite element schemes. The comparison shows that the LB scheme has a similar performance as the one-step second-order Lax–Wendroff scheme: it has little numerical diffusion, but has a slight dispersion error. By changing the relaxation parameter ω the dispersion error can be balanced by a small increase of the numerical diffusion. °c 2000 Academic Press Key Words: lattice Boltzmann; convection diffusion; grid refinements.

1. INTRODUCTION

In the last decade lattice Boltzmann schemes have been successfully applied to the analysis of a variety of complex physical phenomena, such as turbulent flow, natural convection, and multi-phase flow [1–4]. However, less complicated phenomena like (convection) diffusion have hardly been studied [5–8]. This is probably due to the vast reservoir of alternative finite element and finite difference schemes, for solving the convection-diffusion problem. But these simple phenomena are ideal for investigating ways to improve the LB methodology [8]. Hence, in this paper we investigate whether the LB methodology can be extended to irregular grids in a natural way. The problem of irregular lattices has previously been addressed in a few papers [9–12], using either coarse-graining or interpolation techniques. These techniques imply a significant departure from the traditional framework of the lattice Boltzmann scheme, thereby 766 0021-9991/00 $35.00 c 2000 by Academic Press Copyright ° All rights of reproduction in any form reserved.

CONVECTION-DIFFUSION LATTICE BOLTZMANN SCHEMES

767

losing its attractive properties. Moreover, these techniques exhibit significant numerical diffusion and do not satisfy conservation laws [12]. In this paper we will derive the LB scheme for irregular lattices by applying the same framework as for Bravais lattices, which is developed in Refs. [8, 13]. The key element of this framework is that the velocity moments of the equilibrium particle distribution function must equal those of the classical Maxwell–Boltzmann distribution. For convection diffusion it is sufficient that the velocity moments up to second order are satisfied [4]. In that case the scheme is Galilean invariant and will show little numerical diffusion. It has been shown that the constraints for the velocity moments are satisfied for highly symmetric lattices, such as the hexagonal lattice and the nine-velocity square lattice [1, 4, 14]. For simplicity sake, we only consider lattices having lattice (Wigner–Seitz) cells with only twofold rotational symmetry. Hereby, we extend our previous studies on convectiondiffusion problems on orthorhombic lattices [6, 7]. In this paper we include rest particles, which are shown to give a major improvement to the accuracy of the scheme [4]. Before deriving the scheme for irregular lattices, we first apply the framework to derive the convection-diffusion scheme for the orthorhombic lattice with rest particles. Subsequently, the same framework is applied to derive the scheme for convection diffusion on 2-D irregular lattices with rectangular lattice cells. Finally, the LB-schemes are compared with a number of the traditional finite difference and finite element schemes using benchmark problems, in order to assess the merits and shortcomings of the LB schemes. 2. LB SCHEME FOR ORTHORHOMIC LATTICES

In this section, a LB scheme is derived for convection diffusion on orthorhombic lattices. For the convection-difussion problems considered in this paper we assume: (1) isotropic diffusion and (2) an externally imposed velocity field, which is uniform and time independent. Under these assumptions the convection-diffusion equation reads ∂t ρg + u · ∇ρg = D∇ 2 ρg .

(1)

Here ρg is the convected physical quantity, which can be a mass density of a tracer or an energy density (i.e., temperature), u is the velocity field, and D is the (thermal) diffusivity. Lattice Boltzmann schemes describe convection diffusion by the time evolving particle distribution function gi (x, t). This function states the number density of particles at lattice site x and time t moving with velocity ci = 1xi /1t along the lattice link connecting the sites x − 1xi and x. The dynamics on the macroscopic scale is then obtained by summing P the particle distribution over all states; i.e., the density is ρg (x, t) = i gi (x, t). At each time step, the lattice gas particles propagate to neighbouring lattice sites, where they collide with other particles. Furthermore, the particles can change their momentum by interaction with externally imposed fields, such as the velocity field in the convectiondiffusion problem. The propagation and collisions of lattice gas particles are described by the so-called lattice Boltzmann equation, which is a discretisation of the classical Boltzmann equation, having a linearised collision integral. In its most general formulation the lattice Boltzmann equation reads X ¤ £ eq Äi j g j (x, t) − g j (x, t) . (2) gi (x + 1xi , t + 1t) − gi (x, t) = j

768

VAN DER SMAN AND ERNST

eq

Here, gi (x, t) is the local equilibrium distribution, which is invariant under collisions. The operator Äi j controls the collisions between the lattice gas particles. In the more simplified case of the lattice-BGK scheme [14], the collision operator reduces to Äi j = −ωδi j . Because of its computational simplicity this paper is restricted to the lattice-BGK scheme. eq The local equilibrium distribution, gi (x, t), follows from the requirement that its velocity moments must equal the moments of the classical Maxwell–Boltzmann distribution [13]. For second-order accurate solutions the following constraints have to be satisfied at each lattice site for each cartesian component α, β, cf. Ref. [4], X X X

eq

(3)

eq

(4)

gi = ρg

i

ci,α gi = jαeq = ρg u α ,

i eq

eq

ci,α ci,β gi = 5αβ = ρg cs2 δαβ + ρg u α u β .

(5)

i eq

Here jα is a component of the equilibrium mass flux density, 5αβ is a component of the equilibrium momentum flux density tensor, and cs is a free model parameter. Because of our restriction to orthorhombic lattices, constraint Eq. (5) is not satisfied in the case of flow fields u, which are not parallel to one of the principal axes of the lattice. But for uniform flow fields, parallel to one of the principles axes of the lattice, the scheme will be Gallilean invariant. The form of the equilibrium distribution, compatible with the two-fold rotational symmetry of the orthorhombic lattice, and satisfying constraints Eqs. (3)–(5) for α = β, is given by · ¸ ci · u (ci · u)2 eq gi = wi ρg 1 + 2 + cs ci2 cs2 X eq g0 = ρg − gi

if i 6= 0

(6) (7)

i6=0

with weight factors, satisfying

P

i wi

wi =

= 1, given by cs2 2ci2

w0 = 1 −

if i 6= 0 X i6=0

wi = 1 −

(8) 2 cs0 . cs2

(9)

The index i = 0 denotes rest particles. Because w0 must be positive, the thermal velocity cs can only be set to a value in the range 0 ≤ cs ≤ cs0 . Here, cs0 is the thermal velocity for a lattice gas without rest particles. By applying the Chapman–Enskog procedure, cf. [8], one obtains the expression for the diffusion coefficient ¶ µ 1 1 − 1t (10) D = cs2 ω 2

CONVECTION-DIFFUSION LATTICE BOLTZMANN SCHEMES

769

which is an identical expression as derived for diffusion [8], and hydrodynamics [1]. The range of the relaxation parameter is 0 ≤ ω ≤ 2. However, in the range of ω < 1 the scheme is consistent with diffusion for a limited range of large time-scale and long wave-lengths [8]. Hence, the practical range of the relaxation parameter is 1 ≤ ω ≤ 2. 3. LB SCHEME FOR IRREGULAR LATTICES

Starting from the assumption, as observed by Koelman [1], that the constraints for the equilibrium fluxes also hold for irregular grids, a convection-diffusion LB scheme for these grids is derived. Recall that we only consider irregular grids with twofold rotational symmetry and 2d + 1 particle velocities, as sketched in Fig. 1. Furthermore, isotropic diffusion and uniform time-independent flow fields are assumed. The lattice gas particles associated with the particle density distribution gi (x, t) at lattice site x, are thought to be located within the Wigner–Seitz cell. This control volume is defined as the lattice cell with boundaries at the midpoint of the lattice links 1xi and normal to these links. In the case of lattice links with opposite direction having unequal lengths, the lattice site is not in the centre of the Wigner–Seitz cell. This is similar to the so-called cell vertex finite volume method [18]. As in LB schemes on Bravais lattices, the particle velocities are defined by ci = 1xi /1t, so the particles always move to neighbouring lattice sites at subsequent time steps. For irregular grids the velocities may vary with the location of the lattice site. Hence, we denote the velocities as a function of their location: ci (x) = 1xi (x)/1t. Note that particles

FIG. 1. Irregular lattice with rectangular Wigner–Seitz cells (indicated with dashed lines). Also shown are the pre-collision velocity vectors, ci (x) (i = 1, 2, 3, 4), of particles populating lattice site x, located between regions with different lattice spacing, i.e., ci 6= ci∗ .

770

VAN DER SMAN AND ERNST

propagating from the same lattice site but in opposite direction may have different velocities, i.e., ci (x) 6= ci∗ (x), with i∗ indicating the opposite direction of i. It proves to be more convenient to work with the particle number distribution Ni (x, t) = gi (x, t)1V (x), instead of the traditionally used particle number density distribution gi (x, t). Here, 1V (x) is the volume of the Wigner–Seitz cell surrounding lattice site x. Another important notice to make is that in normal physical practice, fluxes are really defined only for the surfaces enclosing the control volume considered. Fortunately, for regular lattices the definition of the equilibrium fluxes at the lattice site does not introduce an error. But for irregular lattices this definition does produce an inconsistent scheme. Hence, the constraints Eqs. (3)–(5) must be applied to the fluxes at the boundaries of the Wigner–Seitz cell, surrounding the lattices site considered. The constraints for an irregular lattice read as X eq Ni (x) = ρg 1V (x) (11) M(x) = i eq

eq

Ni (x) − Ni∗ (x − 1xi ) 1t = ρg (x)(ei · u)1Si (x)

eq

0i (x) =

eq ci (x)Ni (x)

+ ci∗ (x − − 1xi ) 1t ¤ £ = ρg (x)cs2 + ρg (x)(ei · u)2 1Si (x).

eq

Fi (x) =

(12)

eq 1xi )Ni∗ (x

(13) eq

Here, M(x) is the total mass of the particles in the Wigner–Seitz cell located at x. 0i (x) is the net equilibrium mass flux arriving at x through the surface area of the Wigner–Seitz cell, 1Si (x), located between lattice sites x and x − 1xi . Notice that lattice spacing 1xi depends on the location x of the lattice cell. The unit vector ei = ci /ci indicates the direction eq of the particle velocity ci . Fi (x) is the force the lattice gas exerts on the boundary of the eq Wigner–Seitz cell, midway between the lattice sites x and x − 1xi . Or in other words Fi (x) is the the momentum flux arriving at x through the surface area 1Si (x). eq eq By dividing the mass flux 0i (x) and the momentum flux Fi (x) by the surface area of the Wigner–Seitz cell, one obtains respectively the equilibrium mass flux density jαeq (x) and the equilibrium momentum flux density 5eq αα (x), crossing the particular boundary of the Wigner–Seitz cell. Starting from Eqs. (11)–(13) and the velocity set {ci (x)}, one obtains after some straightforward algebra, the expression for the equilibrium distribution, · ¸ ci · u ci · u eq Ni (x) = wi (x)ρg (x)1V (x) 1 + 2 + 2 2 if i 6= 0 (14) cs cs ci X eq eq N0 (x) = ρg (x)1V (x) − Ni (x) (15) i6=0

with the local weight factor wi (x) defined as wi (x) = eq

eq

cs2 . ci (x)[ci (x) + ci∗ (x)]

(16)

After recalling that Ni = gi 1V , one observes that the expressions Eqs. (6)–(7) for the

CONVECTION-DIFFUSION LATTICE BOLTZMANN SCHEMES

771

eq

equilibrium particle density distribution gi also hold for irregular grids. The effects of variable lattice spacing are absorbed in the weight factors wi (x). After having constructed the equilibrium distribution for irregular grids, we arrive at the problem of how to define the lattice Boltzmann equation. This problem is not trivial, as for irregular grids the traditional lattice Boltzmann equation for Bravais lattices, Eq. (2), does not hold for irregular grids. This becomes evident when considering the case of global eq equilibrium (ρg (x) = Ni (x)/1V (x) = ρ0 ). Clearly, a proper lattice Boltzmann equation should leave the global equilibrium distribution invariant. This is clearly not the case for eq eq Eq. (2) which requires that Ni (x − 1xi , t + 1t) = Ni (x, t), which is generally not true for irregular grids. However, because of Eq. (12) invariance of the global equilibrium distribution is obtained by eq

eq

eq

Ni∗ (x − 1xi , t + 1t) = Ni (x, t) − 0i (x, t)1t.

(17)

We assume that also for the general case, Eq. (17) describes the evolution of the equilibrium part of the distribution function. If the description of the evolution of the non-equilibrium part of the particle distribution, neq eq Ni = Ni − Ni , is also known, then we are able to construct the lattice Boltzmann equation for irregular grids. This description can be obtained by observing the behaviour of the nonequilibrium part of the particle distribution in the case of zero flow field u = 0 and a constant density gradient, ∇ρg (x) = constant. For a regular lattice the non-equilibrium particle distribution of a density field with a neq constant gradient is given by Ni = −1V (x)wi ci · ∇ρg 1t/ω, cf. Ref. [8]. This expression is independent of the lattice spacing. Analysis shows that this expression is valid for irregular grids as well. Hence, the non-equilibrium part of the particle distribution should evolve in the same way for both regular and irregular grids, as is described by neq

neq

Ni∗ (x − 1xi , t + 1t) = (1 − ω)Ni∗ (x, t).

(18)

By adding Eqs. (17)–(18), one finally arrives at the complete lattice Boltzmann equation for irregular grids eq

eq

neq

Ni∗ (x − 1xi , t + 1t) = Ni (x, t) − 0i (x, t)1t + (1 − ω)Ni∗ (x, t).

(19)

Here, the equilibrium distribution is given by Eqs. (14)–(15), the equilibrium mass flow neq is given by Eq. (12), and the non-equilibrium distribution function is given by Ni (x) = eq Ni (x) − Ni (x). The expression for the diffusion coefficient, Eq. (10), also holds for irregular lattices. eq eq eq Notice that for regular grids, Ni (x, t) − 0i (x, t)1t = Ni∗ (x, t), and consequently Eq. (19) becomes identical to the lattice-BGK scheme, which is obtained by setting Äi j = ωδi j in Eq. (2). 4. NUMERICAL ANALYSIS

The consistency and accuracy of the LB scheme for both regular and irregular grids are analysed numerically. The analysis is done by comparing numerical solutions to the analytical solutions of some benchmark problems. These benchmarks are:

772

VAN DER SMAN AND ERNST

• The 1-D steady state problem with inhomogeneous boundary conditions and uniform flow field. • The 1-D transient problem with a Gaussian hill as an initial density field in a uniform flow field. 4.1. Gradient resolution. In order to observe the behaviour of the lattice Boltzmann scheme at boundaries with steep gradients, we solve a 1-D benchmark problem considering inhomogeneous boundary conditions and a uniform flow field (u). With the boundary conditions denoted as ρ(0) = ρl and ρ(L) = ρr , the solution of the benchmark reads as ρ(x) = ρl + (ρr − ρl )

[1 − exp(Pe · x/L)] . [1 − exp(Pe)]

(20)

Here Pe = u L/D is the Peclet number. For Peclet numbers Pe ≥ 0.1 the density field has a steep gradient near the right boundary. These gradients are known to induce numerical oscillations in various numerical schemes. Hence, this benchmark is a good test for the ability of the LB scheme to resolve steep gradients. The benchmark is solved for the case of L = 20, 1x = 1, 1t = 1, ρg (0) = 200, and ρg (L) = 100. The velocity field is uniform and positive u(x) = u > 0. The boundary conditions are imposed by the following constraints at the exterior lattice sites, X X gi (x = 0) = 200, and gi (x = L) = 100. (21) i

i

The implementation of these boundary conditions is straightforward: after collision and propagation, new lattice gas particles have to be injected into the computational domain. The amount of particles to be injected is determined by the constraints of Eq. (21). Furthermore, we have set the thermal velocity cs2 = 12 , the grid Fourier number Fo∗ = D1t/1x 2 = 0.1, and the grid Peclet number Pe∗ = u1x/D = 0.1, 0.5, 1, 5. The numerical solutions together with the analytical solutions are shown in Fig. 2a. Observing these results

FIG. 2. Comparison of numerical solution (symbols) with analytical solution (lines) for boundary value problem ρg (0) = 200 and ρg (L) = 100, with L = 20. (a) The solution for a Bravais lattice with lattice spacing 1x = 1 and grid Peclet number Pe∗ = 0.1, 0.5, 1, and 5. (b) The solution for an irregular lattice, which is graded from 1x = 2 (near x = 0) to 1x = 0.1 (near x = L). Here the grid Peclet number assumes the values Pe∗ = 0.1, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0. Note the numerical oscillations in the upper right corner of Fig.2a, induced by the steep gradient near the boundary.

CONVECTION-DIFFUSION LATTICE BOLTZMANN SCHEMES

773

one sees that the LB scheme quite accurately resolves the analytical solution for low grid Peclet numbers, Pe∗ ≤ 1. At higher grid Peclet numbers spurious oscillations occur, which are induced by the sharp gradient at the right boundary. The magnitude of the spurious oscillations can be decreased a little by lowering cs or increasing the Courant number Cr = u1t/1x, but the effect cannot be eliminated. The benchmark is also solved using a lattice-BGK scheme for an irregular lattice, which is refined at the right boundary. The lattice spacing is varied from 1x = 2.0 (x = 0) to 1x = 0.1 (x = L). Because the lattice spacing varies, the grid Peclet and Fourier numbers vary with the location x, i.e., Pe∗ = Pe∗ (x) and Fo∗ = Fo∗ (x). Hence, if in the text below we refer to the grid Peclet and Fourier numbers, it is implied it is with reference to the unit lattice spacing 1x = 1 (if not stated otherwise). The parameter settings of the previous calculations are maintained: cs2 = 12 and Fo∗ = 0.1. The Peclet number Pe∗ is varied from 0.1 to 10. The numerical solutions are depicted in Fig. 2b, together with the analytical solutions. Figure 2b shows that by using grid refinements numerical oscillations are eliminated, even at high Peclet numbers Pe∗ À 1. By decreasing the lattice spacing 1x at locations with steep gradients one lowers the local grid Peclet number Pe∗ (x) = u1x(x)/D to the regime of Pe∗ (x) ≤ 2, where no oscillations occur. If the gradient is small, the local grid Peclet number can be large. 4.2. Transient behaviour. The transient behaviour of the lattice Boltzmann scheme is investigated by solving a 1-D benchmark problem of a Gaussian hill in a uniform velocity field. The numerical solution is calculated by the LB scheme on a regular lattice. The consistency and stability of the scheme is analysed with the method of moments, by which errors in phase velocity, diffusion, and symmetry of the Gaussian hill are calculated. The solution of this benchmark is given by ρ(x) = ρ0

exp[−(x − x0 )2 /2σ 2 (t)] √ . 2π σ (t)

(22)

Here, ρ0 is the initial height of the Gaussian hill, x0 is the initial position, and σ0 the initial width of the hill. Hence, the initial Gaussian profile is described by £ ± ¤ ρg (x, 0) = ρ0 exp −(x − x0 )2 2σ02 .

(23)

Method of moments. The density profiles after a travel time t are analysed using the method of moments [15]. The moments of the density field are defined as Z M0 (t) =

ρg (x, t) d x ≈

X

Z M1 (t) =

xρg (x, t) d x ≈ Z

M2 (t) =

ρg (n1x, t)1x,

X

xρg (n1x, t)1x,

(25)

n

(x − µ)2 ρg (x, t) d x ≈

X (x − µ)2 ρg (n1x, t)1x,

(26)

n

Z M3 (t) =

(24)

n

(x − µ)3 ρg (x, t) d x ≈

X (x − µ)3 ρg (n1x, t)1x. n

(27)

774

VAN DER SMAN AND ERNST

The moments of a Gaussian distribution, travelling with velocity u, are equal to M0 (t) = M0 ,

(28)

M1 (t)/M0 = µ(t) = x0 + ut,

(29)

M2 (t)/M0 = σ 2 (t) = σ02 + 2Dt,

(30)

M3 (t)/M0 = S(t) = 0.

(31)

From the change in time of the mean value, one can obtain the error in the average flow ˆ i.e., the numerical diffusivity, velocity δ uˆ = [µ(te ) − x0 ]/te . The error in the diffusivity δ D, is obtained from ¤ £ 2 σ (te ) − σ02 ˆ − D. (32) δD = 2te The third-order moment is related to the skewness of the distribution. The error in skewness is [15] δ Sˆ =

S(t) . 6te M0

(33)

The benchmark is solved on a Bravais lattice with a lattice spacing 1x = 1 and 128 lattice spacings long. The timestep is set to 1t = 1. Initially, the center of the Gaussian profile is located at x0 = 32. Other parameters are set equal to σ02 = 8 and ρ˜ 0 = 100. The initial particle distribution corresponding with the initial Gaussian hill is set equal to the first-order perturbation distribution, as derived in Ref. [8], eq

gi (x, t = 0) = gi (ρg ) − wi

1t (ci − u) · ∇ρg (x, t = 0), ω

(34)

where ω is the collision frequency defined in Eq. (19). At the boundaries of the lattice, periodic boundary conditions are applied. Two sets of calculations are performed: (1) at a moderate grid Peclet number Pe∗ = 10, and (2) at a high grid Peclet number Pe∗ = 1000. There are two free parameters in the LB scheme, which are varied in our analysis: the Courant number in the range of Cr = u1t/1x between 0.01 ≤ Cr ≤ 1, and the relaxation parameter between the range of 1 ≤ ω < 2. The range of ω < 1 is not investigated, as it is known from our previous study [8] that in this regime inconsistency with diffusion can occur already at moderate gradients. The density profiles at time te = 401x/Cr have been analysed using the method of ˆ found with the method moments. The errors in the diffusivity and the skewness, δ Dˆ and δ S, of moments, are plotted in Fig. 3 as a function of the relaxation parameter ω and the Courant number Cr. The error in velocity is not shown, as the velocity calculated from the first moment is found to be equal to the pre-set value u up to machine accuracy for all simulations. Observing Fig. 3, one sees that the numerical diffusion is small for low values of ω and Cr. For the case ω = 1 the error is zero (up to machine accuracy) for both cases of Pe∗ = 10 and Pe∗ = 1000. The error δ Dˆ increases for ω → 2 or Cr → 1. The increase in numerical ˆ which is preferable for diffusion is accompanied with an decrease of the skewness δ S, calculations with high grid Peclet numbers.

CONVECTION-DIFFUSION LATTICE BOLTZMANN SCHEMES

775

FIG. 3. Contour plots of the errors in diffusivity and skewness of the lattice Boltzmann scheme using the method of moments for the cases Pe∗ = 10 and Pe∗ = 1000. The relaxation parameter is varied in the range 1 ≤ ω < 2 and the Courant number in the range 0.01 < Cr ≤ 1. The diffusivity error is in %, and the plotted value of the ˆ The regions of instability are also indicated. skewness is 100 times δ S.

There are some limits for the values of ω and Cr due to instability of the LB scheme (the spurious oscillations grow exponentially). At high values of ω two regions of stability remain, i.e., Cr → 0 and Cr → 1. At Pe∗ = 10 the region of Cr → 1 is extremely small and at Pe∗ = 1000 the region of Cr → 0 is extremely small. Therefore, they are not shown in Fig. 3. For any value of the Courant number, there is a range of values for ω which gives stable results. By a suitable choice of ω, any combination of the grid Peclet number Pe∗ and Courant number (satisfying the CFL-condition |Cr| ≤ 1) can be reached with the lattice Boltzmann scheme.

5. COMPARISON WITH TRADITIONAL SCHEMES

In order to assess the merits and shortcomings of the LB scheme, its solution of the above investigated benchmark is compared to those of traditional numerical schemes. The following schemes are used: standard Galerkin (which is equivalent with the finite difference scheme with central differencing in space and forward differencing in time), denoted as CDFD; finite difference scheme with “optimal” upwinding and forward differencing in time (FD+) [16]; finite element scheme with first-order streamline upwinding and implicit time integration (SUPG); and a Galerkin finite element scheme with quadratic elements and Adams–Bashfort (semi-implicit) time integration (ABG). It must be noted that the FD+ scheme is equivalent to a second-order one-step Lax– Wendroff scheme [18, 19]. Furthermore, it is also equivalent to the lattice Boltzmann scheme with the relaxation parameter set to ω = 1, as we show in the Appendix. The solutions of the FD+ scheme are actually calculated with our lattice Boltzmann code. The solutions of the other traditional numerical schemes are obtained by using the general purpose finite element package FIDAP [17].

776

VAN DER SMAN AND ERNST

TABLE I Diffusivity and Skewness Errors and L2 -Norm of Various Schemes Pe∗

Scheme

ˆ δ D/D (%)

δ Sˆ (%)

L 2 norm

10 10 10 10 1000 1000 1000 1000

CDFD FD+ LB ABG SUPG FD+ LB ABG

−100.00 0.00 0.44 1.50 45835.00 −0.03 0.43 4.64

1.60 1.40 0.25 0.02 6.28 5.58 0.84 6.64

0.528 0.168 0.033 0.006 0.968 0.449 0.159 0.066

All above-mentioned numerical schemes, except ABG, are low-order schemes (first or second order in space and/or time). ABG is a high-order scheme (fourth order in space and second order in time) and is similar to the more familiar Crank–Nicolson Galerkin scheme. The benchmark problem is solved for the case Pe∗ = 10 and Cr = 0.2, and the case ∗ Pe = 1000 and Cr = 0.4. The solutions of the LB scheme are calculated using the values ω = 1.8 in the case of Pe∗ = 10 and ω = 1.4 in the case of Pe∗ = 1000. The values of the relaxation parameter are well within the region of stability and show little skewness (see Fig. 3). The simulation results are shown in Fig. 4 and Table I. Figure 4 shows the Gaussian profile at time t = 40/Cr. The errors in diffusivity and skewness are calculated with the method of moments and are listed in Table I. Furthermore, we have calculated the L 2 -norm as a measure of the overall accuracy. We have defined the L 2 -norm as sà N ! X 2 L2 = |ρg (xi ) − ρˆ g (xi )| /N . (35) i=1

Here ρg (xi ) is the exact solution at grid point xi and ρˆ g (xi ) is the numerical solution. As indicated in Fig. 4 and Table I, the standard low-order finite difference and finite element schemes (CDFD and SUPG) clearly show the problems that can arise when solving transient convection-diffusion problems: the numerical scheme either shows large spurious oscillations, or has large numerical diffusion, resulting in a poor overall accuracy. Better performance is obtained by the FD+ scheme and the lattice Boltzmann scheme. Both schemes have little or zero numerical diffusion, but have numerical oscillations due to dispersion errors. This behaviour is typical for a second-order scheme. The accuracy of the lattice Boltzmann scheme can be improved over that of the FD+ scheme by choosing an appropriate value for the relaxation parameter ω. Hereby, the numerical oscillations are reduced at the expense of a little numerical diffusion. The overall accuracy of the high-order scheme (ABG) is one order better than all other lower-order schemes. The ABG scheme has little numerical diffusion and dispersion errors (which is an effect of the odd-order truncation error, and hence, is absent in the fourth-order ABG scheme). Irregular lattice. The transient 1-D benchmark is also solved on an irregular grid, for the lattice Boltzmann scheme, the Lax–Wendroff scheme (FD+), and the ABG scheme.

CONVECTION-DIFFUSION LATTICE BOLTZMANN SCHEMES

777

FIG. 4. Comparison of numerical solution (symbols and dashed lines) with analytical solution (solid lines) for the transient 1-D benchmark problem. The LB scheme is compared with various other schemes for grid Peclet numbers and Courant numders: (1) Pe∗ = 10 and Cr = 0.2 (left side of figure), and (2) Pe∗ = 1000 and Cr = 0.4 (right side of figure). For these two cases the values of the relaxation parameter of the LB scheme are respectively ω = 1.8 and ω = 1.4.

The locations of the nodes of the irregular grid xn are specified by xn = n

for n = 0, 1, . . . , 64,

xn+1 − xn = 0.75 + 0.5 cos[π(n − 64)/128]

(36) for n = 65, . . . , 196.

(37)

778

VAN DER SMAN AND ERNST

FIG. 5. Numerical solution (symbols and dashed lines) and analytical solution (solid lines) for the transient 1-D benchmark problem using an irregular grid. Solutions are shown for times t = 40/Cr (left peak) and t = 80/Cr (right peak). The ABG scheme is solved with the grid Peclet number Pe∗ = 100 and the Courant number of Cr = 0.2. The LB scheme is solved with the grid Peclet number Pe∗ = 100, and various Courant numbers and relaxation parameter values, as indicated in the legend of the graphs.

The solution is computed for the grid Peclet number Pe∗ = 100. In Fig. 5 the solutions of the various schemes are depicted for time steps t = 40/Cr and t = 80/Cr. From Fig. 5 and Table II, one can see that the behaviour of the various schemes is similar for regular lattices and irregular grids. The second-order schemes (Lax–Wendroff and LB) show numerical oscilations and the high-order scheme (ABG) shows improved accuracy with little dispersion errors. Again the LB scheme performs better than the Lax–Wendroff

TABLE II L2 -Norm of Benchmark Solution on Irregular Grid for Various Schemes t

Scheme

L 2 norm

40/Cr

LB(ω = 1.0) LB(ω = 1.4) LB(ω = 1.8) ABG

0.295 0.073 0.129 0.017

80/Cr

LB(ω = 1.0) LB(ω = 1.4) LB(ω = 1.8) ABG

0.631 0.184 0.211 0.014

CONVECTION-DIFFUSION LATTICE BOLTZMANN SCHEMES

779

scheme. The results concerning the overall accuracy indicate there are some optimal values for the relaxation parameter and Courant number. It is left to the reader to investigate how to find these values for his type of problem.

6. PERFORMANCE OF THE LB SCHEME IN 2-D

The performance of the LB scheme is also investigated for two-dimensional lattices. Hence, the benchmark of the propagation of a Gaussian profile in a uniform velocity field is solved for both regular and irregular 2-D lattices. The benchmark is solved for the following parameter setting: Pe∗ = 100, Cr = 0.1, ω = 1.4, and σ02 = 8. The direction of the uniform velocity field is taken at an angle with the principle axes of the lattice, i.e., u = u(ˆex + 12 eˆ y ). In this way the occurrence of crosswind diffusion can be investigated. Surface plots of the Gaussian profile at time t = 20/Cr are shown in Figs. 6 and 7 for the Bravais lattice and the irregular grid, respectively. By calculating the moments from the simulation results, using M2,x x =

XX x

M2,yy =

(38)

(y − µ y )2 ρ(x, y),

(39)

y

XX x

(x − µx )2 ρ(x, y)

y

the diffusion coefficients Dx x and D yy in x- and y-direction are estimated. For both the Bravais lattice and the irregular grid it is found that Dx x = D yy = 0.001, implying that diffusion is isotropic and numerical (crosswind) diffusion is insignificant.

FIG. 6. Gaussian profile at t = 20/Cr on a Bravais lattice, solved with the LB scheme. The centre of the profile is located at (16, 16) at t = 0 and has the velocity u = (1, 12 ). Other parameter settings are Pe∗ = 100 and Cr = 0.1.

780

VAN DER SMAN AND ERNST

FIG. 7. Gaussian profile at t = 20/Cr on an irregular 2-D grid, solved with the LB-scheme. The centre of the profile is located at (16, 16) at t = 0 and has the velocity u = (1, 12 ). Other parametersettings are Pe∗ = 100 and Cr = 0.1.

7. CONCLUSIONS

In this paper, a convection-diffusion lattice Boltzmann scheme for irregular grids is presented. Like the LB schemes for Bravais lattices, our scheme follows the traditional two-step mechanism of (1) the collision of lattice gas at the lattices sites, and subsequently (2) the propagation to adjacent lattice sites. Hence, there is no need of any interpolation step or coarse graining step, which is needed in previous attempts to map the LB scheme to irregular grids [12]. The scheme for irregular grids is derived using the same framework, which is traditionally used to derive LB schemes for Bravais lattices [1, 4, 8, 13]. The keypoint of this framework is that velocity moments of the equilibrium distribution equal those of the classical Maxwell– Boltzmann distribution. By solving benchmark problems we have shown that our schemes are consistent with convection diffusion for both regular and irregular grids. It must be noted that we have tested the consistency for only uniform flows. In subsequent papers we will investigate these schemes for the more general case of non-uniform flows. The behaviour of the LB schemes is similar to second-order finite volume/difference schemes, like the Lax–Wendroff scheme. (Even more, we have shown that the LB scheme is identical to the Lax–Wendroff scheme in the case of the relaxation parameter set to ω = 1.) The convection-diffusion LB scheme has little numerical diffusion, but has some numerical dispersion. This behaviour is typical for a second-order scheme. The numerical dispersion, induced by sharp gradients, can be reduced by an appropriate choice of the relaxation parameter or by lattice refinements. Its accuracy has shown to be less than that of higherorder finite element schemes, but is obtained at much lower computational costs. Hence, for (3-D) problems with a large computational domain or with a complicated geometry the simple and straightforward lattice Boltzmann scheme can be a valuable choice. The equivalence of the lattice Boltzmann scheme with the Lax–Wendroff scheme is a theme worthy of further investigation. This might reveal the relationship of the LB scheme

CONVECTION-DIFFUSION LATTICE BOLTZMANN SCHEMES

781

with other numerical schemes and may eventually lead to the crossover of ideas and concepts between the two research areas, which still evolve independently. Finally, we point out the value of studying simple physical phenomenon like (convection) diffusion with the lattice Boltzmann method. These problems are ideal for investigating fundamental aspects of the LB methodology and further extensions to it. As such, we have been able to extend the traditional methodology to irregular grids. Likewise, other extensions of the methodology like variable timestepping, boundary conditions, and adaptive grids, etc., can be studied. Once established for simple physical phenomena, it can take little effort to apply these extensions to more complicated phenomena like hydrodynamics. Hence, the LB schemes presented in this paper can be the starting point for more complex schemes like (1) third-order accurate convection-diffusion schemes (for regular and irregular grids), and (2) hydrodynamics schemes for irregular grids. Both these schemes require that third-order velocity moments are satisfied. These constraints can be satisfied by taking a larger particle velocity set ci (for a 2-D lattice a 13-velocity set with rest particles will probably suffice). If these schemes can be established, they will mean a major step in the further development of the lattice Boltzmann method. APPENDIX: RELATION WITH FINITE DIFFERENCE SCHEMES

The lattice-BGK scheme is identical to a Lax–Wendroff finite difference scheme, if the relaxation parameter is set equal to ω = 1. Below, this relation is shown for a 1-D Bravais lattice. In this case the diffusion coefficient is given by D = 12 cs2 1t, and the grid Fourier number is equal to Fo∗ = 12 cs2 /ci2 . The Courant number is defined as Cr = u/ci . The change of density in time at lattice site x is described by X X eq gi (x, t + 1t) = gi (x − ci 1, t) ρg (x, t + 1t) = =

i eq g1 (x

− 1x, t) +

i eq g0 (x, t)

eq

+ g2 (x + 1x, t).

(40)

2 = ci2 , the equilibrium distribution functions are With cs0

g0 (x) = ρg (x)[1 − 2Fo∗ − Cr2 ] · 1 eq g1 (x − 1x) = ρg (x − 1x) Fo∗ + Cr + 2 · 1 eq g2 (x + 1x) = ρg (x + 1x) Fo∗ − Cr + 2 eq

(41) 1 2 Cr 2

¸

¸ 1 2 Cr . 2

(42) (43)

After substitution in Eq. (40) it follows that ρg (x, t + 1t) − ρg (x, t) µ ¶ 1 2 ∗ = Fo + Cr [ρg (x + 1x, t) + ρg (x − 1x, t) − 2ρg (x, t)] 2 1 − Cr[ρg (x + 1x, t) − ρg (x − 1x, t)]. 2

(44)

This expression is identical to the finite difference scheme with “so-called” optimal upwinding [16], which is actually a one-step second-order Lax–Wendroff scheme [19].

782

VAN DER SMAN AND ERNST

In the limit of low Courant numbers (Cr → 0), the Lax–Wendroff scheme is identical with the forward time central space differencing. Whereas, in the opposite limit of high Courant numbers (Cr → 1), this scheme is identical with forward time full upwind differencing. REFERENCES 1. J. M. V. A. Koelman, A simple lattice Boltzmann scheme for Navier–Stokes fluid flow, Europhys. Lett. 15(6), 603 (1991). 2. J. G. M. Eggels and J. A. Somers, Numerical simulation of free convective flow using the lattice Boltzmann scheme, Int. J. Heat Fluid Flow 15, 357 (1996). 3. E. G. Flekkoy, Lattice BGK model for miscible fluids, Phys. Rev. E 47, 4247 (1993). 4. M. R. Swift et al., Lattice Boltzmann simulations of liquid gas and binary fluid systems, Phys. Rev. E 54(5), 5041 (1996). 5. D. Wolf-Gladrow, A lattice Boltzmann equation for diffusion, J. Statist. Phys. 79(5/6), 1023 (1995). 6. R. G. M. Van der Sman, Lattice Boltzmann scheme for natural convection in porous media, Int. J. Mod. Phys. C 8(4), 879 (1997). 7. R. G. M. Van der Sman, Solving the vent hole design problem with a convection-diffusion lattice Boltzmann scheme, Int. J. Comp. Fluid Dynam. 11(3–4), 237–248 (1999). 8. R. G. M. Van der Sman and M. H. Ernst, Diffusion lattice Boltzmann scheme on an orthorhombic lattice, J. Statist. Phys. 94(1/2), (1999). 9. F. Nannelli and S. Succi, The lattice Boltzmann equation on irregular lattices, J. Statist. Phys. 68, 401 (1992) 10. X. He et al., Some progress in lattice Boltzmann method. Part I. Non-uniform mesh grids, J. Comput. Phys. 129, 357 (1996). 11. N. Cao et al., Physical symmetry and lattice symmetry in the lattice Boltzmann method, Phys. Rev. E 55(1), R21 (1997). 12. H. Chen, Volumetric formulation of the lattice Boltzmann method for fluid dynamics: Basic concept, Phys. Rev. E 58(3), 3955 (1998). 13. G. McNamara and B. Alder, Analysis of the lattice Boltzmann treatment of hydrodynamics, Phys. A 194, 218 (1993). 14. Y. H. Qian, D. d’Humieres, and P. Lallemand, Lattice BGK models for Navier–Stokes equation, Europhys. Lett. 17(6), 479 (1992). 15. P. Muldrup et al., An accurate and numerically stable model for one-dimensional solute transport in soils, Soil Sci. 153(4), 261 (1992). 16. P. Muldrup et al., Removing numerically induced dispersion from finite difference models for solute and water transport in unsaturated soils, Soil Sci. 157(3), 153 (1994). 17. FIDAP v7.0, Fluid Dynamics International (Evaston, IL,1993). 18. J. C. Tannehill, D. A. Anderson, and R. H. Pletcher, Computational Fluid Mechanics and Heat Transfer (Taylor & Francis, London, 1994). 19. S. Succi et al., An integer realization of a Lax scheme for transport processes in multiple component fluid flows, J. Comput. Phys. 152, 493 (1999).