Comparison of passive scalar transport models coupled with the Lattice Boltzmann method

Computers and Mathematics with Applications (

)

–

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

Comparison of passive scalar transport models coupled with the Lattice Boltzmann method V.E. Küng a, *, F. Osmanlic b , M. Markl a , C. Körner a a

Chair of Materials Science and Engineering for Metals, Friedrich-Alexander-Universität Erlangen-Nürnberg, Martensstr. 5, 91058 Erlangen, Germany b Joint Institute of Advanced Materials and Processes, Dr.-Mack-Str. 81, 90762 Fürth, Germany

article

info

Article history: Available online xxxx Keywords: Passive scalar transport Lattice Boltzmann method Advection

a b s t r a c t A complex fluid flow often involves passive scalars that are transported with the flow field but can be neglected regarding the fluid motion, such as concentration distributions, markers, or temperature fields. Using the Lattice Boltzmann (LB) method for fluid dynamic simulations, different approaches for coupling advection solvers to the fluid motion exist. The aim of this work is to give an overview of four different methods used for solving the advection equation. Two of these solvers are based on finite differences, namely the Lax–Wendroff method (Lax and Wendroff, 1960) and an extension from Succi et al. (1999) that artificially reduces numerical diffusion. The other two methods Onishi et al. (2005) and Osmanlic and Körner (2016) take advantage of the already available LB distribution functions, instead of depending on macroscopic quantities only. In the ansatz of Onishi et al. (2005) the local concentration flux is calculated from distribution functions directly, while in a similar approach used by Osmanlic and Körner (2016) an additional case distinction is made between incoming and outgoing flux. Besides investigating the order of accuracy and grid dependency, numerical diffusion is examined closely with the objective to simulate fluid flow that is solely determined by advection and contains no physical diffusion. © 2018 Elsevier Ltd. All rights reserved.

1. Introduction Passive scalar transport means the advection of a scalar field with fluid motion and has to be solved in many simulations of complex fluids ranging from ocean physics [1] to chemical applications [2]. Typical examples for these scalars are tracer particles [3,1], concentration distributions of dilute multicomponent systems [2,4] or solutes [5,6], as well as temperature fields [7]. The Lattice Boltzmann model provides an effective and parallel numerical simulation method for complex fluid systems. Instead of solving the Navier–Stokes equation directly, the Boltzmann equation is solved [8]. Therefore, it is well suited to handle multicomponent [9] and multiphase [10,11] fluids as well as complex boundary conditions, as needed for the simulation of free surfaces [12–14] or porous media [15]. Coupling transport to fluid flow means to solve the convection–diffusion equation for a given velocity field. Convection treatment with finite differences, however, is difficult for low diffusivity and gives rise to wriggles at high Péclet numbers [16]. In order to reduce this behaviour, higher order schemes were developed, such as higher order upwind

*

Corresponding author. E-mail address: [email protected] (V.E. Küng).

https://doi.org/10.1016/j.camwa.2018.01.017 0898-1221/© 2018 Elsevier Ltd. All rights reserved.

Please cite this article in press as: V.E. Küng, et al., Comparison of passive scalar transport models coupled with the Lattice Boltzmann method, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.01.017.

2

V.E. Küng et al. / Computers and Mathematics with Applications (

)

–

schemes [17,18] and Galerkin methods [19]. An extensive review on coupling transport and fluid flow was done by Dawson et al. [20], however not including the LB method. In the context of LB, a multi distribution function method was applied to handle temperature as a passive scalar by Shan [21] and later further developed [22,23]. Merks et al. [24] introduced and validated the Moment Propagation method to introduce advection–diffusion of tracers which was further improved for higher Péclet numbers at low Reynolds numbers [25]. An extension of LB models for anisotropic convection–diffusion equations and development of corresponding boundary conditions was done by Ginzburg [26,27]. In their equilibrium approach antinumerical diffusion can be introduced to counteract numerical dispersion terms generating a method that is especially suited for advection driven flows. A more straightforward approach is the coupling to finite differences [4] or Lax–Wendroff schemes [28]. The design of the LBM allows direct calculation of transport from distribution functions, which was for example utilized by Onishi et al. [29] and Osmanlic et al. [30]. The focus of this work lies on pure advective transport without physical diffusion at low Reynolds numbers, as needed to trace fluid motion or for advection of concentration fields in Lattice Boltzmann models. For this purpose four different solvers, namely (i) the Lax–Wendroff scheme, (ii) an extension to it from Succi et al. [28], (iii) the construction used by Onishi et al. [29] and (iv) a similar construction by Osmanlic et al. [30], are compared. The paper is structured as follows: At first the four methods, regarded in this comparison, are described briefly. The numerical dispersion and the corresponding diffusion constants are investigated subsequently. In Section 3 the order of accuracy is compared for the four methods and their grid dependency is studied. Concluding, the advantages and disadvantages of each solver are discussed. 2. Transport solvers With a given velocity field u(x, t), the scalar field c(x, t) at location x and time t, is calculated from the advection–diffusion equation

∂ c(x, t) + u(x, t) · ∇ c(x, t) − D ∇ 2 c(x, t) = 0 ∂t

(1)

with diffusion constant D. As the focus of this work lies on advective transport the diffusion constant is D = 0 and Eq. (1) is reduced to the advection equation

∂ c(x, t) + u(x, t) · ∇ c(x, t) = 0 . ∂t

(2)

The fluid flow field is defined by the continuity equation

∇ ·u=0

(3)

for constant density ρ , and the Navier–Stokes equation

ρ

(

) ∂u + u · ∇ u = −∇ p + ρν∇ 2 u + Fext ∂t

(4)

with viscosity ν and external force Fext . The fluid movement itself is calculated with the LB method [8,31] using the single relaxation time collision operator (BGK) [32]. The LBM solves the microscopic Boltzmann equation on a discrete Cartesian grid with a discrete set of Q velocities {ci }. In this work simulations were carried out in 2 dimensions with 9 velocities (D2Q9). The corresponding distribution functions are [8] fi (x + ei , t + 1) = fi (x, t) −

) 1( eq fi (x, t) − fi (x, t) ,

τ

(5)

eq

with the relaxation time constant τ . fi is the equilibrium distribution function [33] that corresponds to the velocity direction ei . The macroscopic quantities can be recovered as moments of the continuous distribution function, therefore they are given as sums over the discrete set of functions [8]:

ρ (x, t) =

Q −1 ∑

fi (x, t)

(6)

i=0

u(x, t) =

1

Q −1 ∑

ρ (x, t)

i=0

fi (x, t)ci .

(7)

With the Chapman–Enskog expansion the continuity and Navier–Stokes equations can be recovered, and the relation

τ = 3ν + 0.5

(8)

between relaxation constant and fluid viscosity ν is obtained. Please cite this article in press as: V.E. Küng, et al., Comparison of passive scalar transport models coupled with the Lattice Boltzmann method, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.01.017.

V.E. Küng et al. / Computers and Mathematics with Applications (

)

–

3

The transport solvers regarded in this comparison can be divided into two groups, depending on their basic concept: The first group contains two finite difference based methods that depend on macroscopic quantities, while the second group of solvers takes advantage of the Lattice Boltzmann distribution functions. 2.1. Finite difference schemes Lax–Wendroff scheme. The Lax–Wendroff scheme is a widely used standard method [34,16,35,36] for solving advection problems with high stability. The scalar c at time t + 1 at position x is given by [37,28] c(x, t + 1) =

d 1∑

2

[(

1 d

α=1

) ( ) ] 1 c(x − eα , t) + c(x + eα , t) , + u− − u+ α α

(9)

d

with eα being the unit vectors in Cartesian coordinates and d the dimension. The velocities u± α = uα (x ± eα , t) are the velocities of the neighbouring cells to position x. However, Eq. (9) is actually the finite difference representation of Eq. (1) with the numerical diffusion constant Dlax α =

1 2d

(1 − d u2α )

(10)

instead of the pure advection equation (2) [34]. This leads to numerical dispersion of initial distributions which is especially large for small velocities and reaches its maximum for u = 0. Extension from Succi. In order to eliminate this numerical diffusion, Succi et al. [28] extended the Lax–Wendroff method to explicitly include the diffusion term of Eq. (1) with an ‘anti-diffusion’ constant. It is defined as the negative value of the Lax–Wendroff diffusion constant: Dα = −Dlax α . Solving this term with finite differences leads to c(x, t + ∆t) =

+

1∑ 2

[(

d

α

∑

1

+

u− α

)

c(x − eα , t) +

(

1 d

−

u+ α

)

]

c(x + eα , t)

[c(x + eα , t) − 2c(x, t) + c(x − eα , t)] Dα

α

−

−

1 ∑ 8

[

u− α uβ (x − eα , t) c(x − eα + eβ , t) − c(x − eα − eβ , t)

[

]

α,β̸=α

u+ α uβ (x

+ eα , t) c(x + eα + eβ , t) − c(x + eα − eβ , t)

[

]

]

.

(11)

The last sum is required in order to cancel out cross terms that arise for dimensions larger than 1 [28]. Both the Lax–Wendroff scheme as well as its extension from Succi are based on finite differences and can be coupled to any fluid dynamic solver, since they only depend on macroscopic quantities. If coupling to the LB method, however, also the local distribution functions fi are known for each time step. Instead of depending on macroscopic quantities only, the second group of solvers takes advantage of these distribution functions. 2.2. Lattice Boltzmann schemes The distribution functions fi (x, t) contain information about the fluid flow in velocity direction ei , when i indicates one of the discrete velocities. Both methods in this section regard the advancing of the scalar field c in time, as a quantity which is transported from one cell to its neighbours by the fluid flux given by the distribution functions fi . Instead of the Cartesian unit-vectors eα now the velocity directions ei are used; both are illustrated in Fig. 1. Onishi’s method. In Onishi’s method [29], the flux of the scalar c out of one cell is the local concentration multiplied by a weighting factor: c(x, t)Wi (x). The flux into the cell is correspondingly defined by the distribution function in the neighbouring cell c(x + ei , t) multiplied by a weighting factor W¯i (x + ei ) that corresponds to the fluid movement in opposite direction of i, indicated by the bar. For the scalar quantity at position x and time t + 1 we therefore obtain [29]: c(x, t + 1) = c(x, t) +

Q −1 ∑ i=1

⎡

⎤

⎣−c(x, t)Wi (x) + c(x + ei , t)W¯i (x + ei )⎦       flux out

(12)

flux in

with the weighting factors fi (x) Wi (x) = ∑ − C ωi . j fj (x)

(13)

Please cite this article in press as: V.E. Küng, et al., Comparison of passive scalar transport models coupled with the Lattice Boltzmann method, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.01.017.

4

V.E. Küng et al. / Computers and Mathematics with Applications (

)

–

Fig. 1. Velocity set of D2Q9 method. The black numbers are the fluid directions i, while the Cartesian unit vectors are indicated in blue and numbered with α . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

C ωi again is an ‘anti diffusion’ term, containing the Lattice Boltzmann weight-factors ωi and an adjustable control parameter C that can assume values between 0 and 1. If C = 0 no correction is applied and if C = 1 the numerical diffusion is fully corrected, similar to Succi’s method.

Distinction of cases. Following a similar philosophy as in Onishi’s ansatz we regard the mass flux between cells, as implemented by Osmanlic et al. [30] and Klassen et al. [38]. The mass exchange between a cell at position x and its neighbour in direction ei is [12]

∆mi = f¯i (x + ei , t) − fi (x, t)

(14)

where ¯i again indicates the opposite direction to i. If the mass flux ∆mi is negative, concentration is transported from the cell considered to its neighbour, with ∆ci = ∆mi c(x, t). If the mass change ∆mi is positive, the concentration is transported from the neighbouring cell at x + ei into the cell at position x and therefore: ∆ci = ∆mi c(x + ei , t). Hence, the concentration at time t + 1 is given by [30] c(x, t + 1) = c(x, t) +

∑ i

∆mi m(x, t + ∆t)

{

c(x, t), if ∆mi ≤ 0

c(x + ei , t), if ∆mi > 0.

(15)

Thus, a case distinction is introduced into the transport method. Hereafter the different methods will be named: ‘Lax–Wendroff’ for the Lax–Wendroff scheme [37], ‘Succi’ for Succi’s extension of the Lax–Wendroff scheme [28], ‘Onishi’ for Onishi’s solver [29], and ‘Osmanlic’ if the case distinction is used [30]. 3. Numerical dispersion and diffusion constant In order to investigate dispersion and numerical accuracy (see Section 4), a simple 1D system is chosen: a scalar Gaussian distribution advected by a fluid with constant velocity. The standard deviation of the Gaussian is σ = 50 and the fluid velocity u = 0.01. The initial distribution is located at x0 = 10σ = 500 and the box length is N0 = 1500. 3.1. Numerical dispersion and oscillations Regarding this set-up it becomes apparent that the results are divided into two cases: In the first case we observe numerical dispersion as shown in Fig. 2, while in the other case the distribution functions are not broadened but oscillations appear over time (Fig. 3). In Fig. 2 the initial Gaussian distribution is shown, as well as the distribution after 50 000 time steps, as it is obtained from different solvers. In the graph broadening of the initial distribution and reduction of the maximal value is visible. The Lax–Wendroff method shows the largest numerical diffusion. Onishi’s method without correction terms (C = 0) has less diffusion than Lax–Wendroff, but the Gaussian still spreads severely. However, this will be reduced when the correction term is considered. Finally, the Osmanlic method shows far less diffusion than the other two. Interpreting the spreading of the Gaussian as diffusion, it is possible to deduce the corresponding diffusion constant, as discussed in Section 3.2. The solvers with a correction term hindering numerical diffusion (Fig. 3), maintain the original Gaussian profile. However, oscillations around zero occur in succession of the Gauss peak due to overshooting of the correction term. In order to portray these oscillations in Fig. 3 a slightly altered set-up is chosen: the Gaussian distribution is σ = 5 and the fluid velocity u = 0.1. For comparison the graph also shows Osmanlic’s distribution. Although in many cases small oscillations in a scalar field are bearable, a lot of scalars are physically restricted to positive values and negative values would be un-physical or lead to stability issues. Moreover, these oscillations become quite high Please cite this article in press as: V.E. Küng, et al., Comparison of passive scalar transport models coupled with the Lattice Boltzmann method, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.01.017.

V.E. Küng et al. / Computers and Mathematics with Applications (

)

–

5

Fig. 2. The figure displays the dispersion of a Gaussian distribution with standard deviation σ = 50. The distribution was advected by a fluid with velocity u = 0.01. The graph shows the initial distribution at time t = 0 and after t = 50 000 time steps solved with the methods of Lax–Wendroff (black), Onishi without correction (C = 0, blue), Succi (dashed green) and Osmanlic (red). The dotted lines indicate the decay of the maximal value over time for each solver. The grey line shows the analytical solution. The full length of the system is N0 = 1500 and the initial position of the Gauss-peak is x0 = 10σ = 500. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3. Gaussian distribution with standard deviation σ = 5 advected by a fluid with velocity u = 0.1. The graph shows the initial distribution at time t = 0 and after t = 2000 time steps solved with the methods of Succi (dashed, green) and Onishi with full correction (C = 1, dotted, blue) as well as Osmanlic (solid, red). The analytical solution is shown by the grey solid line. The initial Gauss-peak is positioned at x0 = 50 and the full box length is N0 = 1000. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

for non-continuous initial distributions, for example theta functions, as illustrated in Fig. 4. Both graphs also demonstrate that Succi’s extension of the Lax–Wendroff scheme is equivalent to Onishi’s method, when the numerical diffusion is fully corrected with C = 1. When reducing the correction term in the latter case, the numerical diffusion quickly begins to outweigh oscillation. 3.2. Diffusion constant Due to numerical diffusion the maximum value of a Gaussian distribution decays with

√ A(t) = A0

σ2 , σ + 2Dt 2

(16)

with A being the amplitude of the Gauss distribution, σ the deviation and D the diffusion constant. In Fig. 5 the dimensionless numerical diffusion constant, obtained with Eq. (16), is plotted against the dimensionless fluid velocity. For the Lax–Wendroff scheme we clearly see that the numerical results match the known diffusion constant of Eq. (10). It is zero for velocity u = 1 and has it highest value at zero velocity with DLax = 0.5. A similar equation can be derived for the method used by Osmanlic et al. for the one dimensional case. For fluid flow in +x direction Eq. (15) becomes the forward in time backwards in space (FTBS) finite difference representation of the advection equation. Accordingly the numerical diffusion constant is DOs =

1 2

(1 − u)u .

(17)

Please cite this article in press as: V.E. Küng, et al., Comparison of passive scalar transport models coupled with the Lattice Boltzmann method, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.01.017.

6

V.E. Küng et al. / Computers and Mathematics with Applications (

)

–

Fig. 4. Advection of the theta distribution c(x) = Θ (x) in a fluid with velocity u = 0.01. The distribution at time t = 10 000 is shown for the methods from Succi(dashed), Onishi with C = 1 (dotted), and Osmanlic (solid, red), as well as the theoretical distribution (solid, grey). The full box length is N0 = 500. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5. Numerical diffusion constants of different methods plotted against fluid velocity. The data points are obtained from numerical simulations of a Gaussian function advected with constant fluid velocity (see Fig. 2). The lines are the theoretical curves and not interpolated.

The full derivation is shown in the Appendix. Thus, the diffusion is zero for u = 0 as well as u = 1 and the highest value D = 0.25 is observed at u = 0.5. The diffusion constant of Onishi’s solver is approximately constant over different velocities and can be reduced by increasing C . Since the Lattice Boltzmann method is only correct for small fluid velocities the Péclet number is regarded in that regime (see Fig. 6). It is defined as Pe = u/D and relates advective and diffusive transport. For the Lax–Wendroff method we get PeLax =

2u 1−u

−→ 0, for u → 0

(18)

which means that diffusion outweighs the advection when the transport is slow enough. The same holds for Onishi’s method since the diffusion constant is not velocity dependent. In Osmanlic’s case, however, we get PeOs =

2 1−u

−→ 2, for u → 0 .

(19)

This means that the advection always outweighs the diffusion even for very small velocities. Please cite this article in press as: V.E. Küng, et al., Comparison of passive scalar transport models coupled with the Lattice Boltzmann method, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.01.017.

V.E. Küng et al. / Computers and Mathematics with Applications (

)

–

7

Fig. 6. Péclet numbers corresponding to the diffusion constants shown in Fig. 5.

Fig. 7. Order of accuracy for the four different solvers. The first and second order decays are indicated by dotted and dashed lines, respectively. Onishi’s method is represented without (downwards pointing triangles) and with (upwards pointing triangles) correction term. The results from Onishi with C = 1 and Succi coincide.

4. Order of accuracy and grid dependency 4.1. Order of accuracy The order of accuracy is evaluated from a Gaussian distribution in a fluid with constant velocity. As initial values the resolution of the Gaussian is chosen as σ0 = 5, box length N0 = 200 and u0 = 0.001. The deviation from the theoretical curve, defined as

[ δerr =

Nx −1 1 ∑

Nx

]−1/2 (c(xi , t) − cth (xi , t))2

(20)

i=0

with cth denoting the theoretical curve and Nx the current resolution, was calculated after t = 1000 time steps. While increasing the resolution the relation σx /Nx and the Reynolds number Re = ux Nx /ν were kept constant. The results are represented in Fig. 7 for each solver. The Lax–Wendroff scheme as well as the methods from Succi and Onishi is of second order accuracy, while Osmanlic’s method is only of first order. Nevertheless, the absolute error is significantly lower than the error of the Lax–Wendroff Please cite this article in press as: V.E. Küng, et al., Comparison of passive scalar transport models coupled with the Lattice Boltzmann method, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.01.017.

8

V.E. Küng et al. / Computers and Mathematics with Applications (

)

–

Table 1 Overview over different solvers. Diffusion

Oscillations

Order of accuracy

Grid dependency

Lax–Wendroff Succi

D = 0.5(1 − u2 ) No

No Yes

2nd 2nd

No (x, y)

Onishi

C

2nd

u

1st

u

Osmanlic

=0 D const. D = 0.5u(1 − u)

C

=1

No

C

=0

No –

C

=1

Yes

method and uncorrected Onishi. This behaviour is preserved for higher resolution, since the velocity has to be scaled with the resolution and cannot reach values larger than 1 without violating the Courant condition. 4.2. Grid dependency All previous investigations were carried out with 1D examples, however, dependency of results on the grid orientation is often an issue in LB simulations, due to its Cartesian grid. In order to detect grid dependencies of the different methods it is necessary to regard a two dimensional system with fluid velocity in off-grid direction. To do so we chose the advection of a two dimensional Gaussian

] [ (x − x0 )2 + (y − y0 )2 c(x, y, t = 0) = C0 exp − 2σ 2

(21)

(with amplitude C0 = 1 and centre (x0 , y0 ) = (100, 100)) through fluid with the constant velocity u = 0.01 · (2, 1). Not only is the fluid direction off the grid, but also ux and uy do not have the same value and grid dependency can become apparent easily. Fig. 8 displays the scalar field after t = 10 000 time steps for the different models. The Lax–Wendroff scheme exhibits, as we already know, large diffusion but no grid dependency. Similar behaviour shows Onishi’s method without diffusion correction (Fig. 8b), although with less pronounced dispersion. It can be further decreased by increasing the correction factor C up to values of 0.9 without introducing lattice dependencies (Fig. 8d), since the numerical diffusion still outweighs and smoothes out oscillations. However, when the diffusion term is fully corrected (Fig. 8f), oscillations begin to rise in −u direction. While this scheme was equivalent to Succi’s method in one dimension, in 2D now the differences between both methods become apparent: While for Onishi oscillations depend on the absolute fluid velocity and its direction, the oscillations introduced by Succi’s scheme are dependent on the velocities in grid direction ux and uy . Since ux is larger than uy , mainly the oscillation in x-direction is visible in Fig. 8e. Finally, the method of Osmanlic (Fig. 8c) shows asymmetric diffusion, spreading the Gaussian more along the velocity direction than in the perpendicular direction, due to the dependency of the diffusion constant on the fluid velocity. Yet, the diffusion constants are smaller than for the methods of Lax and Onishi with C ≤ 0.9. 5. Overview and conclusion Generally, all regarded methods either display numerical diffusion or oscillations (compare Table 1). The Lax–Wendroff scheme [37] is a very stable standard method to tackle advection, but the high numerical diffusion at small velocities poses a problem for coupling to the Lattice Boltzmann method. The extension made by Succi et al. [28] was able to eliminate this diffusion, but therefore introduced oscillations in the system. These oscillations are problematic for scalar fields restricted to positive values (e.g. relative concentrations) or for discontinuous and steep initial distributions. The diffusion correction also is grid dependent. Onishi et al. [29] took advantage of the Lattice Boltzmann distribution functions and introduced a scalable correction term, thus, tuning between diffusion and oscillations is possible. However, the suitable correction is velocity dependent and the Péclet number remains small for small velocities. By introducing a case distinction [30] the numerical diffusion was reduced significantly without an extra ‘anti-diffusion’ term and Péclet numbers are larger than 2, but it comes at the cost of obtaining only first order accuracy. As the focus of this work lies on the behaviour of the passive scalar advection schemes, independent of the fluid flow solver, a direct coupling between fluid flow and scalar transport is not regarded. However, these couplings can be found in the literature: Succi et al. compared their method to a standard finite volume based method in a Poiseuille velocity profile [28] and later applied it to the transport of chemical species of catalytic reactive flow for the same geometry [39], however the species distribution still acted as a passive scalar concerning the hydrodynamics. It was also discussed by van der Sman et al. in a comparison between different solvers for the convection–diffusion equation [40]. Onishi et al. validated their transport model for viscoelastic fluids with an oscillating flow [29]. This was also done by Osmanlic et al. for their own method as well as an investigation of the viscoelastic coupling for a planar contraction and a four rolls mill [30]. For future work it will be interesting to compare these four solvers with more complex Lattice Boltzmann models, for example the model from Ginzburg [26], as the ‘pure advection’ method shows good performance for large Péclet numbers. Please cite this article in press as: V.E. Küng, et al., Comparison of passive scalar transport models coupled with the Lattice Boltzmann method, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.01.017.

V.E. Küng et al. / Computers and Mathematics with Applications (

)

–

(a) Lax–Wendroff.

(b) Onishi C = 0.

(c) Osmanlic.

(d) Onishi C = 0.9.

(e) Succi.

(f) Onishi C = 1.

9

Fig. 8. Profile c(x, t) of a two dimensional Gaussian distribution with σ = 10 in a fluid with velocity u = (0.02, 0.01) after 10000 time steps. The results from the methods of Lax–Wendroff (8a), Osmanlic (8c) and Succi (8e) are shown on the right-hand side, while for Onishi’s method three different values for C are represented (8b, 8d, 8f). The dotted lines mark the theoretical form and position of the Gaussian. Please note that Figs. 8a and 8b are scaled differently than the other four.

In summary, finding a suitable scalar transport solver depends on the requirements of the observed scalar field. In principle a method without diffusion and with second order accuracy is preferable, especially if the regarded scalar contains absolute numbers, and has a continuous distribution where oscillations are small and bearable. However, scalar fields often represent relative concentrations [5,2] that are physically limited between 0 and 1 where oscillations below 0 and above 1 cannot be tolerated. Discontinuous distributions, for example Θ -distributions, are not uncommon as well [16,41]. In those cases accepting numerical diffusion is the preferable option. If second order accuracy is required a suitable parameter C , Please cite this article in press as: V.E. Küng, et al., Comparison of passive scalar transport models coupled with the Lattice Boltzmann method, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.01.017.

10

V.E. Küng et al. / Computers and Mathematics with Applications (

)

–

for scaling the correction, can be found. If the field is passive, however, first order accuracy is sufficient because the fluid solver will not be affected and a reasonable Péclet number is more important. For those cases the simple approach used by Osmanlic et al. [30] is a good alternative. The methods have different strengths and weaknesses therefore decisions must be made on a case by case basis. This work can be a guideline to find a suited advection model for a specific application. Acknowledgements The authors gratefully acknowledge funding of the German Research Foundation (DFG) within the project KO 1984/9-2 and the Collaborative Research Center 814 (CRC 814), sub-project B4. Appendix. Diffusion constant DOs We examine the numerical diffusion for Osmanlic’s method for a one dimensional fluid with velocity u in x-direction. The density and mass are constant and set to 1. For the 1D case the LB equilibrium distributions read eq

f1,2 =

1

(1 ± u) (22) 2 and the resulting mass exchange according to Eq. (14) is ∆m1 = −u = −∆m2 . In this case the numerical evolution equation with case distinction (15) becomes a simple forward in time, backward in space finite difference scheme1 : c(x, t + 1) = c(x, t) − u [c(x, t) − c(x − 1, t)] .

(23)

Applying the Taylor series expansion up to the second order

∂ c(x, t) 1 ∂ 2 c(x, t) + + ··· ∂t 2 ∂t2 2 ∂ c(x, t) 1 ∂ c(x, t) + + ··· c(x − 1, t) = c(x, t) − ∂x 2 ∂ x2 c(x, t + 1) = c(x, t) +

(24) (25)

results in the differential equation

∂c ∂c 1 ∂ 2c 1 ∂ 2c +u =− + u 2. 2 ∂t ∂x 2 ∂t 2 ∂x

(26)

Now the second order time derivative on the right hand side can be replaced with the second order space derivative according 2 2 to the advection equation (2) with ∂∂ t 2c = u2 ∂∂ x2c . As result we obtain a advection–diffusion equation with a velocity dependent diffusion constant:

∂ 2c ∂c ∂c 1 +u = u(1 − u) 2 . ∂t ∂ x 2   ∂ x

(27)

=DOs

In the two dimensional case with a D2Q9 grid the diffusion constant becomes a tensor with coefficients parallel and perpendicular to the velocity direction. References [1] P. Gillibrand, M. Herzfeld, A mass-conserving advection scheme for offline simulation of scalar transport in coastal ocean models, Ocean Modell. 101 (2016) 1–16. http://dx.doi.org/10.1016/j.ocemod.2016.02.008. [2] A. Tamburini, G.L. Barbera, A. Cipollina, G. Micale, M. Ciofalo, Cfd prediction of scalar transport in thin channels for reverse electrodialysis, Desalin. Water Treat. 55 (12) (2015) 3424–3445. http://dx.doi.org/10.1080/19443994.2014.959735. [3] V.I. Klyatskin, W.A. Woyczynski, D. Gurarie, Diffusing passive tracers in random incompressible flows: Statistical topography aspects, J. Stat. Phys. 84 (3) (1996) 797–836. http://dx.doi.org/10.1007/BF02179658. [4] A. Tiribocchi, N. Stella, G. Gonnella, A. Lamura, Hybrid lattice boltzmann model for binary fluid mixtures, Phys. Rev. E 80 (2009) 026701. http: //dx.doi.org/10.1103/PhysRevE.80.026701. [5] D.J. Allen, A.R. Douglass, R.B. Rood, P.D. Guthrie, Application of a monotonic upstream-biased transport scheme to three-dimensional constituent transport calculations, Mon. Weather Rev. 119 (10) (1991) 2456–2464. http://dx.doi.org/10.1175/1520-0493(1991)119<2456:AOAMUB>2.0.CO;2. [6] Y. Peng, J. Zhou, R. Burrows, Modelling solute transport in shallow water with the lattice boltzmann method, Comput. & Fluids 50 (1) (2011) 181–188. http://dx.doi.org/10.1016/j.compfluid.2011.07.008. [7] R. Zhang, H. Fan, H. Chen, A lattice boltzmann approach for solving scalar transport equations, Phil. Trans. R. Soc. A 369 (1944) (2011) 2264–2273. http://dx.doi.org/10.1098/rsta.2011.0019. [8] S. Chen, G.D. Doolen, Lattice boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1) (1998) 329–364. http://dx.doi.org/10.1146/annurev.fluid. 30.1.329. [9] E.G. Flekkøy, Lattice bhatnagar-gross-krook models for miscible fluids, Phys. Rev. E 47 (1993) 4247–4257. http://dx.doi.org/10.1103/PhysRevE.47.4247.

1 Forward in space if the velocity is negative. Please cite this article in press as: V.E. Küng, et al., Comparison of passive scalar transport models coupled with the Lattice Boltzmann method, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.01.017.

V.E. Küng et al. / Computers and Mathematics with Applications (

)

–

11

[10] M.R. Swift, E. Orlandini, W.R. Osborn, J.M. Yeomans, Lattice boltzmann simulations of liquid-gas and binary fluid systems, Phys. Rev. E 54 (1996) 5041–5052. http://dx.doi.org/10.1103/PhysRevE.54.5041. [11] X. Shan, H. Chen, Lattice boltzmann model for simulating flows with multiple phases and components, Phys. Rev. E 47 (1993) 1815–1819. http: //dx.doi.org/10.1103/PhysRevE.47.1815. [12] C. Körner, M. Thies, T. Hofmann, N. Thürey, U. Rüde, Lattice boltzmann model for free surface flow for modeling foaming, J. Stat. Phys. 121 (1) (2005) 179–196. http://dx.doi.org/10.1007/s10955-005-8879-8. [13] E. Attar, C. Körner, Lattice Boltzmann model for thermal free surface flows with liquid–solid phase transition, Int. J. Heat Fluid Flow 32 (1) (2011) 156–163. http://dx.doi.org/10.1016/j.ijheatfluidflow.2010.09.006. [14] C. Körner, E. Attar, P. Heinl, Mesoscopic simulation of selective beam melting processes, J. Mater Process. Technol. 211 (6) (2011) 978–987. http: //dx.doi.org/10.1016/j.jmatprotec.2010.12.016. [15] Z. Guo, T.S. Zhao, Lattice boltzmann model for incompressible flows through porous media, Phys. Rev. E 66 (2002) 036304. http://dx.doi.org/10.1103/ PhysRevE.66.036304. [16] W. Shyy, Computational Modeling for Fluid Flow and Interfacial Transport, Dover Publ., Mineola, N.Y., USA, 2006. [17] B. Leonard, A stable and accurate convective modelling procedure based on quadratic upstream interpolation, Comput. Methods Appl. Mech. Engrg. 19 (1979) 59–98. http://dx.doi.org/10.1016/0045-7825(79)90034-3. [18] W. Shyy, A study of finite difference approximations to steady-state, convection-dominated flow problems, J. Comput. Phys. 57 (3) (1985) 415–438. http://dx.doi.org/10.1016/0021-9991(85)90188-3. [19] B. Cockburn, C. Shu, The local discontinuous galerkin method for time-dependent convection–diffusion systems, SIAM J. Numer. Anal. 35 (6) (1998) 2440–2463. http://dx.doi.org/10.1137/S0036142997316712. [20] C. Dawson, S. Sun, M.F. Wheeler, Compatible algorithms for coupled flow and transport, Comput. Methods Appl. Mech. Engrg. 193 (23–26) (2004) 2565–2580. http://dx.doi.org/10.1016/j.cma.2003.12.059. [21] X. Shan, Simulation of rayleigh-bénard convection using a lattice boltzmann method, Phys. Rev. E 55 (1997) 2780–2788. http://dx.doi.org/10.1103/ PhysRevE.55.2780. [22] X. He, S. Chen, G. Doolen, A novel thermal model for the lattice boltzmann method in incompressible limit, J. Comput. Phys. 146 (1) (1998) 282–300. http://dx.doi.org/10.1006/jcph.1998.6057. [23] Y.Q. Zu, Y.Y. Yan, W.P. Shi, L.Q. Ren, Numerical method of lattice boltzmann simulation for flow past a rotating circular cylinder with heat transfer, Internat. J. Numer. Methods Heat Fluid Flow 18 (5–6) (2008) 766–782. http://dx.doi.org/10.1108/09615530810885560. [24] R. Merks, A. Hoekstra, P. Sloot, The moment propagation method for advection–diffusion in the lattice boltzmann method: Validation and péclet number limits, J. Comput. Phys. 183 (2) (2002) 563–576. http://dx.doi.org/10.1006/jcph.2002.7209. [25] D. Yu, S.S. Girimaji, A.J.C. Ladd, Revised moment propagation method for scalar transport, Phys. Rev. E 78 (2008) 056706. http://dx.doi.org/10.1103/ PhysRevE.78.056706. [26] I. Ginzburg, Equilibrium-type and link-type lattice boltzmann models for generic advection and anisotropic-dispersion equation, Adv. Water Resour. 28 (11) (2005) 1171–1195. http://dx.doi.org/10.1016/j.advwatres.2005.03.004. [27] I. Ginzburg, Generic boundary conditions for lattice boltzmann models and their application to advection and anisotropic dispersion equations, Adv. Water Resour. 28 (11) (2005) 1196–1216. http://dx.doi.org/10.1016/j.advwatres.2005.03.009. [28] S. Succi, H. Chen, C. Teixeira, G. Bella, A.D. Maio, K. Molvig, An integer lattice realization of a lax scheme for transport processes in multiple component fluid flows, J. Comput. Phys. 152 (2) (1999) 493–516. http://dx.doi.org/10.1006/jcph.1999.6242. [29] J. Onishi, Y. Chen, H. Ohashi, A lattice boltzmann model for polymeric liquids, Prog. Comput. Fluid Dyn. Int. J 5 (1) (2005) 75–84. http://dx.doi.org/10. 1504/PCFD.2005.005819. [30] F. Osmanlic, C. Körner, Lattice boltzmann method for oldroyd-b fluids, Comput. & Fluids 124 (2016) 190–196. http://dx.doi.org/10.1016/j.compfluid. 2015.08.004. special Issue for ICMMES-2014. [31] P. Lallemand, L.-S. Luo, Theory of the lattice boltzmann method: Dispersion, dissipation, isotropy, galilean invariance, and stability, Phys. Rev. E 61 (2000) 6546–6562. http://dx.doi.org/10.1103/PhysRevE.61.6546. [32] P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems, Phys. Rev. 94 (1954) 511–525. http://dx.doi.org/10.1103/PhysRev.94.511. [33] L.-S. Luo, Theory of the lattice boltzmann method: Lattice boltzmann models for nonideal gases, Phys. Rev. E 62 (2000) 4982–4996. http://dx.doi.org/ 10.1103/PhysRevE.62.4982. [34] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes 3rd Edition: The Art of Scientific Computing, third ed., Cambridge University Press, New York, NY, USA, 2007. [35] W.F. Ames, Numer. Methods Partial Differential Equations, Academic Press, New York, USA, 1977. [36] W.G. Collins, An accurate variation of the two-step lax-wendroff integration of horizontal advection, Q. J. R. Meteorol. Soc. 109 (459) (1983) 255–261. http://dx.doi.org/10.1002/qj.49710945913. [37] P. Lax, B. Wendroff, Systems of conservation laws, Comm. Pure Appl. Math. 13 (2) (1960) 217–237. http://dx.doi.org/10.1002/cpa.3160130205. [38] A. Klassen, V.E. Forster, C. Körner, A multi-component evaporation model for beam melting processes, Modelling Simul. Mater. Sci. Eng. 25 (2017). http://dx.doi.org/10.1088/1361-651X/aa5289. [39] S. Succi, M. Adamo, G. Bella, M. Bernaschi, Multi-representation techniques for multi-scale simulation: Reactive microflows in a catalytic converter, Mol. Simul. 25 (1–2) (2000) 13–26. http://dx.doi.org/10.1080/08927020008044109. [40] R. van der Sman, M. Ernst, Convection–diffusion lattice Boltzmann scheme for irregular lattices, J. Comput. Phys. 160 (2) (2000) 766–782. http: //dx.doi.org/10.1006/jcph.2000.6491. [41] C. Hirsch, Numerical Computation of Internal and External Flows, Wiley, New York, USA, 1990.

Please cite this article in press as: V.E. Küng, et al., Comparison of passive scalar transport models coupled with the Lattice Boltzmann method, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.01.017.