Free surface Neumann boundary condition for the advection–diffusion lattice Boltzmann method

Journal of Computational Physics 301 (2015) 230–246

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

Free surface Neumann boundary condition for the advection–diffusion lattice Boltzmann method Matthias Markl ∗ , Carolin Körner Chair of Metals Science and Technology, Friedrich-Alexander-Universität Erlangen-Nürnberg, Martensstr. 5, 91058 Erlangen, Germany

a r t i c l e

i n f o

Article history: Received 17 June 2015 Accepted 15 August 2015 Available online 28 August 2015 Keywords: Neumann boundary condition Lattice Boltzmann method Free surface Advection–diffusion equation

a b s t r a c t The main objective of this paper is the derivation and validation of a free surface Neumann boundary condition for the advection–diffusion lattice Boltzmann method. Most literature boundary conditions are applied on straight walls and sometimes on curved geometries or fixed free surfaces, but dynamic free surfaces, especially with fluid motion in normal direction, are hardly addressed. A Chapman–Enskog Expansion is the basis for the derivation of the advection–diffusion equation using the advection–diffusion lattice Boltzmann method and the BGK collision operator. For this numerical scheme, a free surface Neumann boundary condition with no flux in normal direction to the free surface is derived. Finally, the boundary condition is validated in different static and dynamic test scenarios, including a detailed view on the conservation of the diffusive scalar, the normal and tangential flux components to the free surface and the accuracy. The validation scenarios reveal the superiority of the new approach to the compared literature schemes, especially for arbitrary fluid motion. © 2015 Elsevier Inc. All rights reserved.

1. Introduction The lattice Boltzmann (LB) method [1] is commonly applied to hydrodynamic fluid flow [2–4] described by the Navier– Stokes equations. Additionally, the LB method is suitable to solve the diffusion and advection–diffusion equation. Possible applications are thermal flows [5,6], multi-component flows [7] or solute transport in porous media [8]. Most of the LB implementations use the Bhatnagar–Gross–Krook (BGK) collision operator [9], which limits its applications to isotropic diffusion phenomena. By replacing the BGK with a multi-relaxation time (MRT) operator, an anisotropic diffusion behavior is possible [10]. A Neumann boundary condition specifies the solution of a first derivative at a boundary. Regarding the advection– diffusion LB method, the boundary condition is defined on the flux, which is the spatial derivative of the diffusive scalar. This paper focuses on a zero flux condition in normal direction to the free surface. Additionally, there are no physical restrictions in tangential direction. Neumann and other boundary conditions are usually applied on straight walls [8]. Towards curved or free surface boundaries there are only few approaches in literature. The easiest way to impose a boundary condition is the bounce-back scheme [11], resulting in a zero flux in normal direction. Unfortunately, the tangential parts also vanish. Modified bounceback schemes, altered by the boundary direction [12,13], still affect the tangential flux. A Neumann boundary condition is

*

Corresponding author. E-mail address: [email protected] (M. Markl).

http://dx.doi.org/10.1016/j.jcp.2015.08.033 0021-9991/© 2015 Elsevier Inc. All rights reserved.

M. Markl, C. Körner / Journal of Computational Physics 301 (2015) 230–246

231

presented in [14], which interpolates the macroscopic quantities at the boundary, and demonstrates an appropriate reconstruction of the tangential flux. However, if interpolation is not possible, the fall-back boundary condition is the bounce-back rule. A similar interpolation approach basing on the bounce-back approach and using the particle distribution functions is presented in [15]. Fixed curved boundaries are examined and a first order accuracy for the Neumann boundary condition with arbitrary fluxes is achieved. Most of these schemes are derived for static boundary problems, like simulations of porous media or curved channels. But considering free surface flows, especially movements in normal direction occur. All mentioned schemes are not appropriate to model these kinds of moving boundary problems [14], because the first order moment of the LB method is not conserved resulting in an unphysical concentration of the diffusive scalar at the boundary depending on the velocity. The free surface LB method described in [5] provides a suitable boundary condition for arbitrary free surface flows. Nevertheless, some drawbacks are identified during the revision of this approach. Therefore, after a short introduction into the free surface advection–diffusion LB method, a new approach is derived, where these drawbacks are eliminated or minimized. Finally, the boundary condition is validated in different static and dynamic test scenarios, including a detailed view on conservation of the diffusive scalar, the normal and tangential flux components to the free surface and the accuracy. 2. Free surface advection–diffusion lattice Boltzmann method The single phase-continuum conservation equations to simulate incompressible fluid transport are the incompressible Navier–Stokes equations. Suppose an additional diffusive scalar  , e.g., representing the thermal energy density or the concentration of an additive. Due to diffusion and the dynamics of the fluid the diffusive scalar is distributed described by the advection–diffusion equation

      ∂t  + ∇ · u = ∇ · a  ∇  ,

(1)

where t is the time and u the divergence free macroscopic velocity and a the diffusion coefficient, which can either depend on  or is a constant value.   LB models base on particle distribution functions (pdf) f x, v, t , describing the probability of finding a particle with microscopic velocity v at position x at time t. The basic idea for the derivation of the LB method [2,4,16,17] is to solve the linear transport equation for pdfs in the physical momentum space. The D3Q19 stencil [18], a finite set of 19 discrete velocities c i with lattice weights ωi , discretizes the microscopic phase space. The domain is covered by a regular lattice consisting of cubic cells with side length δx . In each cell the diffusive scalar is modeled by the discretized pdfs f i and the macroscopic quantity is computed by

=



fi.

(2)

i

For all cells the pdf values collide and stream to the neighbor cells using the lattice displacement vector e i = c i δt and the time resolution δt . The so-called BGK collision operator [9] uses a single relaxation time τ to relax towards the local equilibrium, which reads in the discretized form









f i x + e i , t + δt = f i x, t + eq

where f i

eq 

fi

δt 

τ

eq 

fi







x, t − f i x, t ,

is the Maxwell equilibrium distribution [16]





x, t = ωi  1 +

√

ci · u c 2s

 (4)

,

where c s = δx /( 3δt ) is the lattice speed of sound. The relaxation time

a

(3)

= c 2s (

τ − 0.5δt ).

τ is related to the diffusion coefficient a by (5)

The Mach number describes the ratio of the characteristic fluid velocity to the lattice speed of sound. In the incompressible flow limit, i.e., for small Mach numbers Ma < 0.1, the advection–diffusion equation (1) is derived by a Chapman–Enskog expansion in Appendix A. The validation scenarios require a Dirichlet boundary condition for straight walls, which ensures a constant value  w by a non-equilibrium bounce-back approach [12]





eq 

f i x, t = f i











 w + f ieq  w − f i x, t ,

(6)

where i denotes the inverse direction of i. A free surface advection–diffusion LB method is necessary for the simulation of moving interfaces between immiscible gas and liquid fluids. It has to be guaranteed that the gas phase is separated by a closed interface layer from the fluid phase [19]. Thus each lattice cell has a certain state: gas, interface or liquid. To ensure mass conservation a volume of

232

M. Markl, C. Körner / Journal of Computational Physics 301 (2015) 230–246

fluid approach is used, where a fill level ϕ ∈ [0; 1] is defined for each lattice cell. The normal n of the free surface then is computed by a modified gradient of the fill levels [20]. Due to the large density difference between the gas and the fluid, the gas phase is neglected, i.e. there are no pdfs available in gas cells. Because of this fact missing pdfs after the streaming step in interface cells have to be reconstructed by the free surface boundary condition. Due to the dynamics of the liquid, cell conversions between the cell states are executed. All cells are able to change their state but a direct change from a gas to a liquid cell and vice versa are not allowed [19]. In this paper, the fluid motion is achieved by direct modifications of the fill levels without hydrodynamic computations. 3. Free surface Neumann boundary condition In this section the most suitable literature approach is analyzed and the drawbacks are identified. Subsequently, a modified Neumann boundary condition is derived on the basis of results of a Chapman–Enskog expansion. The algorithm is explained, the zero flux condition is verified and the section closes with a summary including expected impacts of the method improvements. Without loss of generality, the dependency of  the  diffusion coefficient on the diffusive scalar () is skipped, and the dependency on the current position and time step x, t is omitted for a better readability. 3.1. Literature review The thermal free surface LB method described in [5] provides a Neumann boundary condition for arbitrary free surface flows, where some drawbacks are reviewed. The total flux q in the moving frame of reference for an arbitrary cell is









computed by the in- and outstreaming pdfs, which are denoted by f i x, t and f i x + e i , t , respectively. Converting the microscopic view of the pdfs into the macroscopic view using Equations (3), (5), (A.19), (A.20) and (A.23), the flux is approximated by [5]





q x, t =



 



c i − u f i x, t +

i

≈ −2a∇  − 2

a

∂t c 2s





 



ci − u f i x + ei , t .

(7)

i



u .

(8)

The right hand side contains, besides the diffusive part, an additional term representing the time derivative of the advective flow. After identifying the flux in an interface cell, the normal part is inverted and applied to the free surface to ensure the zero flux boundary condition. Similar to bounce-back approaches, this scheme reduces the tangential flux. The first application of this model are metal foams, where thin walls with a thickness of only a few cells are simulated [19]. Therefore, during the computation of the inverse flux, this boundary condition explicitly modifies the tangential part, which establishes the heat conduction through these thin walls. However, the approximation of the tangential flux at the free surface is the drawback of this computation, because the approximation introduces errors into the system. Additionally, the recommendation in [12] to omit restrictions on the flux in tangential direction at least to leading orders is ignored. Unlike to hydrodynamics, a prescribed value at the boundary is in general not compatible to solutions of the advection–diffusion equation [12]. This formulation assumes a non-computable tangential flux at the free surface. The compensating source term F i , which causes the flux to vanish at the free surface is given by [5]



Fi =

2 c 2s





c i ωi · Q +

4 c 4s









i , c i −u ·n<0 c i

ωi · Q







c i ωi · u ,

0,





c i − u · n < 0,

(9)

else,

where Q is the inverse flux. The authors state a second order accuracy of the boundary condition with respect to the fluid velocity, which is reviewed by

⎛







F i ci − u =

i



2 c 2s





i , c i −u ·n<0





c i c i ωi · Q −

⎞

2⎜ ⎝  c 2s

 





⎟





c i ωi · Q ⎠u ≈ Q − O u 2 .

(10)

i , c i −u ·n<0

This computation assumes the reduction of the third lattice symmetry in Equation (A.8) for the half of the lattice direc- tions to c 2s /2 · I . This is the case for a symmetric split, where each opposite split velocities of the lattice stencil, c i − u





and c i − u , are divided into one pointing towards the gas phase and one pointing towards the fluid phase. However, this assumption is not always satisfied, depending on the orientations of the fluid velocity and free surface normal, as illustrated in Fig. 1. The right configuration demonstrates an asymmetric split, because both horizontal split velocities are pointing towards the gas phase. The split criteria operates in the moving reference frame (see condition in Equation (9)), which is visible by the modified lattice velocities. The increasing fluid velocity from the left to the right sketch results in an asymmetric instead of a symmetric split. Due to the insufficient symmetry an additional error term is added to Equation (10),

M. Markl, C. Körner / Journal of Computational Physics 301 (2015) 230–246

233

Fig. 1. Split criteria for free surface boundary condition of [5]. Symmetric split of opposite directions due to slow and advantageously directed fluid velocity (left). The increasing and disadvantageously directed fluid velocity causes an asymmetric split, because the split velocity pointing to the right crosses the free surface (right).

which only depends on the lattice configuration. In three dimensions this drawback is aggravated, because the asymmetric split of opposite directions is strengthened. Therefore, a new approach basing on the same idea of applying the inverse flux is presented, where these drawbacks are eliminated or minimized. 3.2. Derivation The condition for the flux in the moving frame of reference in normal direction n to the free surface has to vanish



!



0 = n · −τ c 2s ∇  .

(11)

This condition is well chosen to omit the pdfs like in Equation (7). It exploits the flux term identified in  outstreaming  Equation (8) by the instreaming pdfs f i x, t . By means of a Chapman–Enskog expansion, the first order moment in the moving reference system is given by combining Equations (A.19), (A.20) and (A.23)





c i − u f i = −τ c 2s ∇  − τ ∂t





u .

(12)

i

In the next step the flux is isolated from the additional term on the right hand side of Equation (12). The time derivative is approximated by a Taylor-Series expansion

  1         ∂t u ≈  x, t u x, t −  x, t − δt u x, t − δt , δt

(13)

where a first order approximation is sufficient to achieve second order accuracy [21]. Inserting Equation (13) into Equation (12) and replacing the unknown diffusive scalar at the current time step results in

−τ c 2s ∇  =





ci − u f i +

i

=





ci − u f i −

i

    τ      x, t u x, t −  x, t − δt u x, t − δt

δt

   τ   x, t − δt u x, t − δt ,

δt



u = 1 −

τ δt

 u.

(14)

Finally, the zero flux Neumann boundary condition is reformulated by inserting Equation (14) in Equation (11)



!

0=n·





   τ  c i − u f i −  x, t − δt u x, t − δt .

i

 

δt

(15)

There are two main differences to the literature approach. The first is the consideration of only instreaming pdfs in Equation (12) in contrast to the application of outstreaming pdfs in Equation (7), which solely modifies the corresponding prefactor compared to Equation (8). Secondly, the first time derivative of the advective flow is taken into account, which is not considered in the literature model. 3.3. Algorithm As described in the previous section, the numerical model omits the gas phase, where all gas cells summarized in the set Cg , and no pdfs are streaming into neighboring interface cells. Therefore, the zero flux condition is applied by an appropriate reconstruction of these missing pdfs. The algorithm consists of the same three steps as described in [5], whereby the single computations of the second and third step are modified:

234

M. Markl, C. Körner / Journal of Computational Physics 301 (2015) 230–246

1. Estimation of unknown pdfs In this step, all available pdfs from neighboring interface as  well as fluid cells are taken into account and only the unknown pdfs from the gas phase are approximated, where x − e i ∈ Cg . Assuming small deviations per time step it is appropriate to approximate them by their outstreaming values resulting in the estimation f i 













f i x, t =



f i x + ei , t ,



f i x, t ,



x − e i ∈ Cg ,

(16)

else.

2. Computation of flux in normal direction Subsequently, the resulting flux in normal direction to the free surface qn is computed by inserting Equation (16) in Equation (15)





qn = n ·



   τ  c i − u f i −  x, t − δt u x, t − δt . 





(17)

δt

i

3. Application of inverse flux In the last step, the resulting flux in normal direction is inverted and applied to the approximated pdfs to ensure the zero flux condition



fi =

  f i − f i− , c i − u · n < 0, f i ,

(18)

else,

where f i− is the corresponding compensating source term to F i in Equation (9). The aim is to find a term f i− , which prevents additional error terms during the reconstruction process, in contrast to Equation (10). Therefore, the assumption of the reduction of the third lattice symmetry for the half stencil to c 2s /2 · I is dropped. Additionally, the intermediate terms of the partial sum over c i ωi in Equation (10) needs to be omitted, because they achieve many different values depending on the split condition. These considerations result in the idea to exploit the third Binomial formula











ωi c i − u  c i + u  = ωi c i c i − u  u  .

(19)

A similar ansatz as in the previous section to review the literature approach is used to find an appropriate formulation for the compensating source term



!

nqn =









c i − u  f i− .

(20)

i , c i −u ·n<0

The compensating   source term is mathematically constructed using the residual terms of Equation (19) on the left hand side ωi c i + u  and an inverse matrix M −1 basing on the right hand side, which reduces to the identity when summing over all reconstructed pdfs f i by Equation (15):

 





f i− = ωi c i + u  · M −1 · n qn ,

M=









ωi c i c i − u  u  .



(21)

i , c i −u ·n<0

Consequently, the resulting flux applied on the pdfs is exactly restored, which is shown in the further course of this section. Applying the inverse of the reconstruction matrix M circumvents the lattice symmetry problem, which occurs in the approach of [5]. For the new approach there are no assumptions on the consideration of opposite lattice directions. The numerical model requires an invertible reconstruction matrix M, which is verified by its regularity demonstrated in Appendix B. 3.4. Zero flux condition Equation (15) defines the zero flux condition, which is satisfied by inserting Equations (18) and (21)

⎛

!

⎞

⎜ 

=n·







c i − u  f i− −



i , c i −u ·n<0

i

=n·





c i − u  f i −

0=n·⎝

 



   τ  c i − u f i − nqn −  x, t − δt u x, t − δt

 



i



I − nn ·



  ⎟ τ   x, t − δt u x, t − δt ⎠

δt



δt



 i





   τ  c i − u f i −  x, t − δt u x, t − δt 



δt

.

M. Markl, C. Körner / Journal of Computational Physics 301 (2015) 230–246



For any normalized vector n the factor n · I − nn

235

results in a zero vector. This method achieves the zero flux condition

without an additional error term during reconstruction. Therefore, the new approach in principle reaches the LB method accuracy in normal direction. However, a small error is introduced during the approximation of the unknown pdfs, which is not revertible by the reconstruction algorithm. In tangential direction no further assumptions are made, but the tangential part is also modified by the approximation procedure, which lowers the accuracy in tangential direction. 3.5. Summary The thermal free surface Neumann boundary condition varies from the literature approach of [5]. The following enumeration summarizes all differences and discusses their expected impact: 1. Computation of flux   The derived scheme of this section applies only the instreaming pdfs of the current position and time step f i x, t   in Equation (12), whereby the literature approach relies additionally on the outstreaming pdfs f i x + e i , t . During the revision of the literature approach, only the assumption of small differences per time step allows the derivation of the flux in Equation (8). Besides this difference, the new scheme also takes the first time derivative of the advective flow in Equation (12) into account. However, comparing this term to the flux, it acts on a smaller scale, because the velocity is restricted to small values by the low Mach number assumption of the LB method. Solely modifying these aspects and applying them on the validation examples of the next section cause minor effects on the final behavior of the literature approach. Therefore, these modifications rather clarify and reduce the necessary assumptions of the derivation. 2. Consideration of tangential flux Hardly any thermal Neumann boundary condition in literature takes the tangential aspects of the flux into account, because it should not be restricted [12]. However, the approach of [5] determines the tangential flux and applies it into the reconstruction algorithm of the missing pdfs. This treatment has a major influence on the thermal behavior at the free surface and is separately studied in the next validation section. 3. Application of compensating source term The literature review reveals a second order error term during the application of the compensating source term in Equation (10). Additionally, the reconstruction assumes a symmetric split of the lattice directions, which is not generally satisfied and results in a further error term. In contrast, the Neumann approach is independent of the split and achieves a perfect application of the source term without an error term as described in the previous paragraph. Regarding the original approach of [5], the parts of the compensating source term in tangential direction are large. Therefore, an asymmetric split has a large expected impact on the boundary behavior. If the tangential aspects are neglected, an impact is only expected for dynamic scenarios, because the additional error term depends on the fluid velocity. 4. Validation In this section the thermal free surface Neumann boundary condition (NB) is validated with static and dynamic scenarios. All scenarios are compared to three reference implementations. The first approach (A0) is the original method derived in [5]. The second scheme (A1) uses the same algorithm, except of neglecting all tangential parts during the reconstruction, due to the annotations in [12]. This constellation demonstrates the single effects of the tangential and the additional error term of the first reference boundary condition. The third reference implementation (BB) is the bounce-back approach. Furthermore, the approaches are compared to a one dimensional solution without a free surface (1D), if there is an appropriate setup. The scenarios include static or dynamic, partially with cell conversions, and lattice aligned or inclined setups in order to verify the universal accuracy of the boundary conditions. First, two static diffusion scenarios without a hydrodynamic velocity are examined. In the beginning an isolated system within a sphere is studied, where different diffusive scalar distributions are applied. Subsequently, the Dirichlet boundary condition on straight walls is added by a diffusion through a cylinder between two walls. Afterwards, two dynamic diffusion scenarios with a velocity tangential and normal to the free surface are studied. In the first scenario a cylinder is moved in tangential direction and in the second scenario a plate is skewed by a normal velocity. The setups per scenario comprise the accuracy, which is the power of the decrease of the absolute error with the spatial resolution. At best, the boundary condition does not decrease the accuracy of the LB method, which is a second order accurate scheme. In all numerical simulations a zero flux is applied on the free surface, δt scales with δx2 and the analytical velocity is coupled into the thermodynamics, i.e., no hydrodynamics are simulated. 4.1. Spherical scenario Two spherical diffusion scenarios are used in [14] to validate their Neumann boundary conditions and to compare the results with the bounce-back scheme (G2). A non-conserving (G0) and a conserving (G1) scheme are derived, which both base on an interpolation scheme at a surface point. In this section, these scenarios are repeated and the numerical results are compared. The bounce-back case G2 is recomputed, while the results of G0 and G1 are directly taken from literature.

236

M. Markl, C. Körner / Journal of Computational Physics 301 (2015) 230–246

Fig. 2. Setup of spherical scenario, where a sphere with radius rs and midpoint xs is applied in the center of the simulation domain and is initialized with different diffusive scalar distributions from red to blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

A sphere with radius rs = 0.4 and midpoint xs is located in the center of a cubic domain with side length L = 1. Depending on the free surface boundary condition, a volume- or center-based approach to initialize the cells is used. The schemes A0, A1, NB and BB use the volume-based approach, where all cells crossing the sphere are computed. The boundary conditions G0 to G2 apply the center approach, where only cell center points contained by the sphere are considered. The analytical solution is given by







 

 x, t = exp −aα 2t  x .

(22)

The parameter α is adjusted to ensure a zero total diffusive scalar inside the sphere [14]. The diffusion coefficient a is computed by Equation (5), the corresponding relaxation time of the scenario and δt = δx2 . Fig. 2 shows two different initial conditions for the diffusive scalar distribution. In the radial symmetric scenario this condition is

 

1



 x =  

2

πα

r x

sin



 





 

r x = x − xs  ,

αr x ,

α = 4.493 41/rs .

(23)

In contrast to the radial symmetric scenario, where only the normal parts of the flux are present, the plane symmetric scenario involves tangential parts with

 

1



 x =   r x



2



1

 

  sin α r x

πα αr x

     θ x = arccos x2  /r x , 



eq

 



  

− cos α r x

    φ x = arctan x1  , x0  ,

The pdfs are initialized with

f i x, t = 0 = f i





 x, t , −τ c 2s ∇  x, t



  

sin θ x

  

cos φ x ,

x = x − xs ,

α = 2.081 58/rs .

(24) (25)

(26)

in accordance to Equation (A.23) [14]. The analytical solution is compared with the numerical simulation by the absolute error

  2        1  2 L a (t ) =  f i x, t −  x, t · , x∈D

i

nc

(27)

where nc is the number of cells. In Fig. 3 the evolution of the absolute error L 2a depending on the grid resolution at t = 0.15 is plotted. The relaxation time of 1.1 δt is examined in order to compare the numerical results with [14]. The dashed lines indicate a grid resolution of δx = 0.05. The upper and lower dashed–dotted lines indicate first and second order accuracy, respectively, and Table 1 summarizes the accuracies. At t = 0.15 the absolute error in the radial symmetric scenario is almost constant. The BB approach has the highest errors and results in an accuracy of 1.28, whereby it perfectly conserves the diffusive scalar in the system. The accuracy

M. Markl, C. Körner / Journal of Computational Physics 301 (2015) 230–246

237

Fig. 3. Grid resolution dependent absolute error L 2a for τ /δt = 1.1 at t = 0.15. All schemes show an accuracy between first and second order (indicated by the dash-dotted lines), except for G0, G1 and G2 in (a). Three different error levels, where BB is worst, error level of NB, A0 and A1 is one fifth of BB, and errors of G0, G1 and G2 are one order of magnitude lower (a). Because L 2a is not constant at t = 0.15 in (b), only the accuracies are compared.

Table 1 Accuracy of spherical scenario for accuracies around 1.5.

τ /δt = 1.1 at t = 0.15. Extreme accuracies for G0, G1 and G2 (radial), whereby all other simulations reach acceptable

Symmetry

NB

A0

A1

BB

G0 [14]

G1 [14]

G2

Radial Plane

1.50 1.68

1.49 1.41

1.50 1.68

1.28 1.37

0.93 1.53

2.27 1.53

3.03 1.57

of NB and the two reference implementations A0 and A1 is approximately 1.5 and the corresponding errors are significantly smaller than for the BB approach, although the total diffusive scalar in the system is slightly increased. Nevertheless, the lowest errors are achieved with the center-based initialization routine, which is used for all reference conditions G0 to G2. The conserving schemes G1 and G2 achieve accuracies larger than 2, whereby both graphs are strongly dependent on the grid resolution. Finally, the non-conserving scheme G0 results in a first order accuracy. Concerning the absolute error value, the latter three schemes in general lower the errors up to a factor of 1/5 compared to the former three approaches. In the plane symmetric scenario, all approaches reach an accuracy of approximately 1.5. This is due to the more complex initial diffusive scalar distribution in the system. Nevertheless, similar to the radial symmetric setup, there are three error levels, where the BB scheme is worst and solely G0 is significantly better than the other approaches. The results of the boundary conditions G0 to G2, cause similar accuracies but smaller errors than the Neumann boundary condition NB, due to the more accurate initialization routine. The scenario also demonstrates the similar behavior of A1 and NB, which cause slightly better accuracies than A0. Especially for NB a second order accuracy is expected, because the LB method is a second order accurate scheme, but the reconstruction at the wrong surface point degrades the result. The second order accuracy is only achievable, if the boundary condition is applied at second order accurate surface points instead of the cell centers. 4.2. Cylindrical scenario In this section the evolution of a linear profile through a cylinder between two walls with different energy densities is studied. The first setup validates to an axis aligned cylinder, while the second setup applies two arbitrary inclination angles to the lattice grid. Both setups are compared to an analytical solution and the first configuration is additionally examined by a one dimensional simulation, where the center line of the cylinder is simulated. The setup of the cylinder scenario is sketched in Fig. 4. A cylinder with radius rc = 0.4 and length L = 1 is located in the center xc of a cubic domain with side length L + 1.5 rc . The cylinder is initialized with 0 = 1 and equilibrium pdfs. The cylinder caps are closed with two solid walls with diffusive scalar offsets from 0 of e = 1 and zero for the left and right wall, respectively. The cylinder and the walls are rotated by the inclination angles φ and θ . The domain is periodically

238

M. Markl, C. Körner / Journal of Computational Physics 301 (2015) 230–246

Fig. 4. Setup of cylindrical scenario between two walls with φ = 18.0◦ , θ = 4.5◦ . Sketch shows slice through the center of the x2 axis with the initial setup of the inclined cylinder surrounded by gas and the left heated wall at t = 0 (left). Three dimensional view with diffusive scalar distribution from 0 + e (red) to 0 (blue) in wireframed wall and filled cylinder cells at t = 0.5 (right). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

connected in all dimensions and the domain size is large enough to have the complete setup surrounded by at least one layer of gas cells. The normal vector pointing from the left towards the right wall is

1 xv =   v , v 

v = ( 1, tan φ, tan θ ).

(28)

The diffusion coefficient a is computed by Equation (5), the corresponding relaxation time of the scenario and δt = δx2 . The analytical solution converging to a linear profile is given by [22]





 x, t =  

2

π 

e

∞  (−1)n n =1

dist x = x − xc +

 

n L 2

 sin nπ

 

dist x L

 exp −n2 π 2

at L2

 

 +

dist x L

e + 0 ,

(29)

 xv

· xv ,

(30)

where dist x is the distance from the current position to the heated wall. The solution is compared to the numerical simulation by the absolute error computed with Equation (27). First, the absolute error L 2a is studied for arbitrary inclination angles for τ /δt = 0.8 at t = 4. The interpolated numerical results for φ and θ in {0.0◦ , 4.5◦ , . . . , 45.0◦ } are illustrated in Fig. 5. As expected, the numerical results are symmetric to the dashed line through the origin. All approaches cause similar errors for any configuration, except for the reference implementation A0. These results deviate in general about one order of magnitude. Additionally, there exist some configurations, where the lattice dependent effect on the error is visible, e.g., there is a negative effect at (45.0◦ , 0.0◦ ) or (4.5◦ , 0◦ ) and a positive effect at (45.0◦ , 45.0◦ ). The circles indicate the configurations analyzed in the further course of this section. The first setup at (0.0◦ , 0.0◦ ) is chosen because it is additionally comparable to a 1D solution with no free surface. The second configuration at (18.0◦ , 4.5◦ ) represents a setup for any arbitrary inclination angles. In Fig. 6 the absolute error L 2a depending on the spatial resolution is plotted. The dashed lines indicate the connection to the previous results at δx = 0.05. The dashed–dotted lines indicate first (top) and second order accuracy (bottom). Table 2 summarizes the corresponding accuracies. As expected, the one dimensional simulation 1D reaches second order accuracy. The absolute error value is also more than one order of magnitude lower than for all setups with a free surface. Depending on the inclination angle all other schemes nearly achieve first order, except the original reference implementation A0. In the inclined configuration the finer grid resolution nearly has no effect. Although, there are local deviations, the total diffusive scalar in the system is reached in the steady state. Summarizing, the tangential parts of the flux should not be modified in the free surface reconstruction algorithm, because the resulting absolute errors are increased, and the flux profile in the core of the cylinder is disturbed. Like in the spherical scenario, the reference implementation A1 and the Neumann boundary condition NB cause similar errors. The next scenarios introduce a fluid velocity and leave the steady state.

M. Markl, C. Körner / Journal of Computational Physics 301 (2015) 230–246

239

Fig. 5. Inclination angle dependent absolute error L 2a for τ /δt = 0.8 and δx = 0.05 at t = 4. Similar errors for NB, A1 and BB. Strong lattice dependency with higher errors for A0. Dashed line and circles indicate symmetry axis and configurations for further analysis. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Fig. 6. Grid resolution dependent absolute error L 2a for τ /δt = 0.8 at t = 4. Errors of 1D are one order of magnitude lower compared to all other schemes (a). First order accuracy and similar errors for all approaches except for A0 with a significant lower accuracy (b).

240

M. Markl, C. Körner / Journal of Computational Physics 301 (2015) 230–246

Table 2 Accuracy of cylindrical scenario for τ /δt = 0.8 at t = 4. Expected second order accuracy for the 1D simulation. All other schemes are around first order accuracy, except A0 in the inclined configuration, where almost no error reduction is achieved with higher grid resolutions.

φ

θ

NB

A0

A1

BB

1D

0.0◦ 18.4◦

0.0◦ 4.5◦

1.29 0.76

0.78 0.08

1.29 0.76

0.86 1.43

2.00

Fig. 7. Setup of parallel movement scenario at t = 0. Sketch shows slice through the center of the x2 axis with the initial setup of the cylinder surrounded by gas (left). The diffusive scalar distribution from 20 (red) to zero (blue) and the velocities (white arrows in horizontal planes) are shown (right). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

4.3. Parallel movement scenario The next two sections switch from steady state to dynamic scenarios. First, an axis aligned cylinder in a periodic domain with a velocity parallel to the free surface is studied. This scenario refers to [21], where the one dimensional advection– diffusion LB scheme is validated. These results are compared with the three dimensional scenarios. The setup of the parallel movement scenario is illustrated in Fig. 7. An axis aligned cylinder of radius rc = 0.4 and length L = 1 is located in the center of a periodic domain with side length L. The diffusive scalar is initialized with [21]

 

 x = 0 + 0 cos (α x0 ), with

α=

2π L

(31)

,

0 = 1 and the temporal evolution of the velocity is   t u 0,0 , u 0 (t ) = u 0,0 cos 2π Ma = , tp

(32)

cs

with the maximum velocity u 0,0 and the period time t p = 1. In the one dimensional setup, only the center line along the cylinder axis is simulated. The time dependent analytical solution for the advection–diffusion equation is [21]

















 x, t = 0 + 0 exp −α 2at cos (α x0 ) cos α x0  (t ) + sin (α x0 ) sin α x0  (t ) , x0  (t ) =

u 0 (t )t p 2π



sin 2π

t tp

(33)



(34)

,

where the diffusion coefficient a is computed by Equation (5) with the corresponding relaxation time of the scenario and δt = 3.2 δx2 . The absolute error per period is evaluated by [21]

   2 nt        1  2 La =  f i x, nδt −  x, nδt · , n=0 x∈D

i

nt nc

(35)

where nt and nc are the number of time steps and cells, respectively. The grid resolution is modified in order to determine the model accuracy. The numerical results are plotted in Fig. 8, where the dashed dotted lines indicate first and second order accuracy. In Table 3 the accuracies are summarized.

M. Markl, C. Körner / Journal of Computational Physics 301 (2015) 230–246

241

2

Fig. 8. Relaxation time and lattice spacing dependent absolute error per period La for u 0,0 = 0.1, t s = 1 and δx = 0.05 (a) or τ /δt = 0.8 (b). Similar behavior with minimum at 1.0 for NB and A1 schemes, A0 reaches higher errors at low relaxation times, and BB already attains its minimum at 0.7 (a). Nearly second order accuracy of 1D, accuracy level between 1 and 2 for NB and A1, and other approaches do not reach first order accuracy (b).

Table 3 Accuracy of parallel movement scenario. 1D achieves expected second order accuracy of the LB method. NB and A1 reach an accuracy of 1.6. Modifying the tangential flux deteriorates the accuracy to first order, which is even worse compared to the BB scheme. NB

A0

A1

BB

1D

1.60

0.90

1.60

1.06

1.99

The one dimensional simulation 1D nearly reaches the second order accuracy of the LB scheme. NB and A1 achieve an accuracy of 1.6. This is a good result, regarding the high velocity gradient of the scenario and the approximation of the unknown pdfs by the values of the previous time step. As in the previous sections, the BB and the original reference implementation A0 only reach first order accuracy. Summarizing, this scenario demonstrates the disadvantage of modifying the tangential parts of the thermal flux in the free surface reconstruction algorithm. Although a velocity is added to the scenario, the Neumann boundary condition NB and the modified reference implementation A1 cause similar errors. This indicates the improvement of the compensating source term application only for velocities in normal direction to the free surface, which is finally studied in the last scenario. 4.4. Normal movement scenario In this section an infinite plate is skewed by a velocity normal to the free surface. This scenario demonstrates the difference between the Neumann boundary condition NB and the other approaches. The setup of the normal movement scenario is sketched in Fig. 9. A plate with a side length L 0 = 1 and a height of H = 0.79 is located in the center of a cubic domain with equal side length L 0 and a height of L 2 = 1.2. The height is chosen to have the free surface exactly located in the center of a cell with a spatial resolution of δx = 0.01. All domain boundaries are periodic, which mimics an infinite elongated plate. The plate is initialized with 0 = 1 and equilibrium pdfs. A velocity in normal direction u 2 is skewing the plate







u 2 x, t = u 2,0 sin (α x0 ) sin (α x1 ) cos 2π u 2,0 = 2π

h 2t p

,

t tp



,

α= Ma =

2π L0

u 2,0 cs

(36)

, ,

(37)

where h is the free surface amplitude of the block and t p the period time of 0.5, where δt = δx2 . The first simulation setup is taken at an amplitude of 0.01, where no cell conversions occur. The second setup is chosen near the Mach number limit of 0.1, represented by an amplitude h = 0.20. The absolute error per period is computed by Equation (35) in combination with the analytical diffusive scalar  x, t = 0 . The grid resolution is modified in order to determine the model accuracy. The numerical results are shown in Fig. 10, where the dashed dotted lines indicate first and second order accuracy. The accuracies are presented in Table 4.

242

M. Markl, C. Körner / Journal of Computational Physics 301 (2015) 230–246

Fig. 9. Setup of normal movement scenario. Sketches show a slice of the domain through the x1 axis at L 0 /4 at t = 0 and after a quarter period t = t p /4 including the velocity and skewed block (left and center). Isosurfaces on the fill level indicate free surfaces of the skewed block, colored slices and arrows illustrate velocities in top and bottom direction at t = t p /8 (right). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

2

Fig. 10. Grid resolution dependent absolute error per period L a for τ /δt = 0.8, u 2,0 = 0.1 and t s = 0.5. Error of BB approach increases with lattice resolution due to increased number of cell conversions. A0 and A1 show quite similar results at approximately second order accuracy. NB reaches lowest errors, while cell conversions destroy advantage and cause the convergence to A0/A1 (b).

Table 4 Accuracy of normal movement scenario. Second order accuracy of NB and A0 is destroyed by cell conversions and tangential modifications, respectively, in contrast to A1. BB scheme is insufficient due to negative accuracy. h

NB

A0

A1

BB

0.01 0.20

2.65 1.43

1.57 1.94

2.14 2.09

−0.02 −0.19

The BB approach results in a negative accuracy. These large errors occur because of the diffusive scalar concentration in interface cells, where the fill level increases, while the diffusive scalar vanishes in cells, where the fill level lowers. Unless there are no cell conversions this effect is symmetric and vanishes over time. But, when cell conversions are necessary, always cells with a low value are deleted and new cells with a high value are generated. Thus, the total diffusive scalar in the system is constantly increasing. This behavior obviously justifies the insufficiency of this approach in the scenario of moving free surfaces. In the original approach A0 the tangential parts destroy the second order accuracy, which is achieved by the modified scheme A1, in the case of no cell conversions. Regarding converted cells, the reset of the pdfs to the equilibrium values

M. Markl, C. Körner / Journal of Computational Physics 301 (2015) 230–246

243

increase the absolute error, but then both schemes achieve a second order accuracy. Finally, the best results are reached by NB with up to one order of magnitude lower error values and second order accuracy in the case of no cell conversions. Unfortunately, the increasing amount of cell conversions with lower lattice spacings degrades the result and the scheme achieves an accuracy of approximately 1.4. This scenario demonstrates the insufficiency of the bounce-back approach for arbitrary dynamic diffusion scenarios. The normal movement to the free surface highlights the advantages of the Neumann boundary condition NB. In this scenario, the improved application of the compensating source term during the NB scheme has the most impact. With increasing fluid velocities, the improvements to the literature approach A0 [5] as well as the modified scheme A1 are obvious. 5. Conclusion In this paper, a free surface Neumann boundary condition for the advection–diffusion LB method is derived after the analysis of a similar literature approach [5]. The validation scenarios demonstrate the advantages of the new approach. The modification of the flux computation has a minor effect in all scenarios. With the spherical diffusion scenarios, comparable results to other literature boundary conditions are achieved. The inclined cylinder scenario demonstrates the disadvantage of modifying the tangential parts of the heat flux density. The results of the tangential movement scenario agree with one dimensional reference simulations. The most impressive results are achieved by the normal movement scenario, where on the one hand, the insufficiency of the bounce-back approach and, on the other hand, the superiority of the Neumann boundary condition to all other approaches is illustrated. In this setup, the most important aspect is the improvement of the application of the compensating source term. Both boundary conditions assume, that the surface point is directly located on the lattice point, thus analytical solutions depending on a certain shape are hardly achieved, like it is shown in the spherical tests. Therefore, the application of this boundary condition at the free surface point may improve the behavior. The hydrodynamic free surface boundary condition proposed in [5] follow the same reconstruction procedure of estimation, recalculation and application of a compensating term on the unknown pdfs. Thus, the improvements of this Neumann boundary condition are also applicable to a hydrodynamic free surface boundary condition. Acknowledgements We gratefully thank the German Research Foundation (DFG) for funding our research in the Collaborative Research Centre 814, CRC814, SFB814, project B4. Appendix A. Derivation of the advection–diffusion lattice Boltzmann method In this appendix the advection–diffusion equation (1) is derived by a Chapman–Enskog expansion on the discretized BGK collision operator inspired by [23]. A small parameter is used, which scales with δt . It expresses the physical effects on different time scales. With this parameter all quantities have the same order of magnitude and the different scales are introduced by k . To expand the collision operator i a Taylor-Series is adapted to approximate the collision operator by

   

i = f i x + e i , t + δt − f i x, t =

∞ n   n =1

n!

∂t + ∇ · c i

n 

f i x, t



     2  2 ∂t + 2∂t ∇ · c i + ∇∇ : c i c i f i x, t . ≈ ∂t + ∇ · c i + 

(A.1)

2





During the course of this section f i always refers to the current position in time x, t . Now each component of the approximation is expanded, whereby the corresponding limits are the 0 scale for ∇ and the 1 scale for ∂t , f i and i . The results are applied to Equation (A.1) and collected by equal order of to 0

(i ) = 0, 

(1 ) (0) ( 0)

i = ∂t + ∇ (0) · c i f i ,   

1  (0) 2 1 (2 (0 (1 (1 (0 (0

i ) = ∂t ) + ∇ (0) · c i f i ) + ∂t ) + + ∂t ) ∇ (0) · c i + ∇ (0) ∇ (0) : c i c i f i ) . ∂t 2

2

(A.2) (A.3) (A.4)

The macroscopic values of the diffusive scalar  and the first and second order moment tensors b and B, are gained by summation of the pdfs f i . Subsequently, they are expanded like f i on order k to

(k) =

 i

(k)

fi ,

b(k) =

 i

(k)

ci f i ,

B (k) =

 i

(k)

ci ci f i .

(A.5)

244

M. Markl, C. Körner / Journal of Computational Physics 301 (2015) 230–246

The LB approach operates on a lattice stencil, whereby the D3Q19 stencil is applied

1

(12 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1)T , ⎛ ⎞T δx ⎝ 0 1 −1 0 0 0 0 1 1 1 1 −1 −1 −1 −1 0 0 0 0 ⎠ 0 0 0 1 −1 0 0 1 −1 0 0 1 −1 0 0 1 1 −1 −1 , c= δt 0 0 0 0 0 1 −1 0 0 1 −1 0 0 1 −1 1 −1 1 −1

ωi =

which provides the following lattice symmetries to ensure isotropy



(A.6)

36



ω i = 1,

i



ω i c i = 0,

i

ωi c i c i = c 2s I ,



i

(A.7)

ω i c i c i c i = 0.

(A.8)

i

These symmetries are used to take moments from the equilibrium distribution of Equation (4). The requirement for zeroth order moment conservation or diffusive scalar conservation is an invariant collision operator



!

i = 0.

(A.9)

i

Due to the scale separation it has to be satisfied for each part of the expanded collision operator. The requirement is already fulfilled by Equation (A.2) on the 0 scale. The 1 scale of Equation (A.3) leads to !

0=



1

0)

(i ) = ∂t(



i

( 0)

fi

+ ∇ (0)

i



(0)

ci f i

0

= ∂t( ) (0) + ∇ (0) · b(0) .

(A.10)

i

2 scale Equation (A.4) results in !  (2) 0=

i

On the

i

=

 

(1)

+

∂t

i

(1) (0)

= ∂t 

1  (0) 2 1 (0 + ∂t ) ∇ (0) · c i + ∇ (0) ∇ (0) : c i c i ∂ 2 t 2



(0)

fi

 

(0 (1 + ∂t ) + ∇ (0) · c i f i )

1  (0) 2 (0) 1  + ∂t(0) ∇ (0) · b(0) + ∇ (0) ∇ (0) : B (0) + ∂t(0) (1) + ∇ (0) · b(1) . ∂ 2 t 2

+

(A.11) (0)

Equation (A.10) is inserted into Equation (A.11) in order to eliminate the second order time derivative ∂t !

(1) (0)

0 = ∂t



1

1

0

0

+ ∂t( ) ∇ (0) · b(0) + ∇ (0) ∇ (0) : B (0) + ∂t( ) (1) + ∇ (0) · b(1) . 2

Finally, Equations (A.10) and (A.12) are recombined using the corresponding



!

(0) (0)

0 = 0 ∂t





+

1

(A.12)

2

+ ∇ (0) · b(0)

(1) (0)

∂t 



scale

1 0 1 0 + ∂t( ) ∇ (0) · b(0) + ∇ (0) ∇ (0) : B (0) + ∂t( ) (1) + ∇ (0) · b(1) 2 2 1

(0

1

(0

= ∂t (0) + ∇ · b(0) + ∂t ) ∇ · b(0) + ∇∇ : B (0) + ∂t ) (1) + ∇ · b(1) . 2

2



(A.13)

The macroscopic advection–diffusion equation (1) describes the evolution of the diffusive scalar  under advection by the fluid velocity u. This equation is recovered from the zeroth order moment conservation law, the equilibrium distribution of Equation (4) and the discretized BGK collision operator [9]

i = −

δt 

τ

eq 

fi − fi

(A.14)

.

Without loss of generality, the dependency of the thermal diffusivity and the relaxation time on the temperature are skipped for a better readability. The collision operator given in Equation (A.9) ensures the zeroth order moment conservation. The idea for the reconstruction is to take Equation (A.13) and replace all unknown terms. Therefore, the collision operator of Equation (A.14) is expanded to ( 0)

=−

δt 

(0)

eq

− fi τ 

δt (1) (1 )

i = − fi .

i

τ

fi

,

(A.15) (A.16)

M. Markl, C. Körner / Journal of Computational Physics 301 (2015) 230–246

245

Simplifying the expansion by using Equations (4), (A.2) and (A.3) leads to ( 0)

fi

(1)

fi

eq

= fi , =−

τ

(1)

δt

=−

i

τ

(0)

+ ∇ ( 0) · c i ω i  1 +

∂t

δt





(0)

Taking moments of the equilibrium function (4) for f i

 ( 0) =



eq 

fi





, u = 

i

b ( 0) =



eq 

ci f i



, u = 





i eq 

ci ci f i



, u = 

ci · u

ωi 1 +

i

i

B ( 0) =







i

(A.18)

ωi c i c i 1 +

(A.19)

 = u ,

c 2s

i

(A.17)

results in

= ,

ci · u



c 2s

 .



c 2s

ωi c i 1 +

ci · u

ci · u

(A.20)

 = c 2s  I .

c 2s

(A.21)

The missing macroscopic terms are gained by Equations (A.10), (A.18) and (A.20)

 (1 ) =



(1)

=−

fi

i

b(1) =



(1)

ci f i

τ δt

=−

i

0

∂t( )  −

τ δt

τ δt

  ∇ ( 0) ·  u = 0,

c 2s ∇ (0)  −

τ δt

0) 

∂t(

(A.22)



u .

(A.23)

By inserting Equations (A.19) to (A.23) into Equation (A.13)

1

0 = ∂t (0) + ∇ · b(0) +

2

1

∂t(0) ∇ · b(0) + ∇∇ : B (0) + 2 ∂t(0) (1) + ∇ · b(1) 2

  

1 = ∂t  + ∇ · u + ∂t ∇ · u + ∇ c 2s ∇  − ∇ · τ c 2s ∇  + τ ∂t u 2 2         1 1 2 = ∂t  + ∇ · u − ∇ c s τ − δt ∇  − τ − δt ∇ · ∂t u , 



1



(0)



2

2

in combination with the definition of the thermal diffusivity a in Equation (5) and by neglecting the error term ∇ · ∂t the advection–diffusion equation is recovered

  ∂t  + ∇ · u = ∇ a∇ .



u



(A.24)

Appendix B. Invertible reconstruction matrix The numerical model of the thermal free surface Neumann boundary condition requires an invertible reconstruction matrix M, which is verified by its regularity. The regularity is achieved if all eigenvalues are unequal to zero. First, a symmetric split is assumed, i.e., all opposite split velocities are separated in two different sets pointing to the gas or material, respectively. Consequently, the corresponding eigenvalues are

λ1,2 =

c 2s 2

λ3 =

,

c 2s 2

−

1 3

 1−

τ δt

2





u·u ,

(B.1)

using the lattice symmetry in Equation (A.8) and the weights of the D3Q19 stencil given in Equation (A.7). The first and second eigenvalues are always unequal to zero. Setting the third eigenvalue to zero and computing the root of the equation for τ results in



τ1,2 = 1 ∓



3 cs



  δt .

(B.2)

2 u 

The LB method is designed for low Mach numbers Ma < 0.1, where the maximum fluid velocity is Ma · c s . The relaxation time τ has to be greater than δt /2 due to stability issues of the LB method. Therefore, the relaxation time is always larger than the first root

1



τ > δt > τ1 > 1 − 2



3 1 2 Ma



δt .

(B.3)

246

M. Markl, C. Körner / Journal of Computational Physics 301 (2015) 230–246

The second root limits the maximum relaxation time to





τ < τ2 < 1 +

3 1



2 Ma

(B.4)

δt .

Choosing the relaxation time in accordance to Equation (B.4) ensures the regularity of the reconstruction matrix. Finally, the symmetric split assumption is dropped, which results in modified factors of the eigenvalues, which are all larger than zero. The worst case constellation of lattice directions, where only the five directions pointing in positive direction of an axis are taken into account, modifies the third eigenvalue to

λ3 =

c 2s 6

−

1 6

 1−

τ

2

δt





u·u ,

(B.5)

which further reduces the upper bound of the relaxation time to



τ < τ2 < 1 +

1 Ma

 δt .

(B.6)

Because all eigenvalues are unequal to zero, if τ is chosen according to the worst case of Equation (B.6), the reconstruction matrix is always invertible for low Mach numbers. References [1] S. Succi, The Lattice Boltzmann Equation – For Fluid Dynamics and Beyond, Clarendon Press, 2001. [2] F. Higuera, J. Jimenez, Boltzmann approach to lattice gas simulations, Europhys. Lett. 9 (1989), http://dx.doi.org/10.1209/0295-5075/9/7/009. [3] H. Chen, S. Chen, W.H. Matthaeus, Recovery of the Navier–Stokes equations using a lattice-gas Boltzmann method, Phys. Rev. A 45 (1992), http://dx.doi.org/10.1103/PhysRevA.45.R5339. [4] S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329–364, http://dx.doi.org/10.1146/ annurev.fluid.30.1.329. [5] 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. [6] R. Ammer, M. Markl, U. Ljungblad, C. Körner, U. Rüde, Simulating fast electron beam melting with a parallel thermal free surface lattice Boltzmann method, Comput. Math. Appl. 67 (2) (2014) 318–330, http://dx.doi.org/10.1016/j.camwa.2013.10.001. [7] X. Shan, G. Doolen, Diffusion in a multicomponent lattice Boltzmann equation model, Phys. Rev. E, Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics 54 (4 Suppl. A) (1996) 3614–3620, http://dx.doi.org/10.1103/PhysRevE.54.3614. [8] Q. Kang, D. Zhang, S. Chen, X. He, Lattice Boltzmann simulation of chemical dissolution in porous media, Phys. Rev. E, Stat. Nonlinear Soft Matter Phys. 65 (3) (2002) 036318, http://dx.doi.org/10.1103/PhysRevE.65.036318. [9] 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. A 94 (3) (1954) 511–525, http://dx.doi.org/10.1103/PhysRev.94.511. [10] 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. [11] H.-B. Huang, X.-Y. Lu, M.C. Sukop, Numerical study of lattice Boltzmann methods for a convection–diffusion equation coupled with Navier–Stokes equations, J. Phys. A, Math. Theor. 44 (5) (2011), http://dx.doi.org/10.1088/1751-8113/44/5/055001. [12] 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. [13] H. Yoshida, M. Nagaoka, Multiple-relaxation-time lattice Boltzmann model for the convection and anisotropic diffusion equation, J. Comput. Phys. 229 (20) (2010) 7774–7795, http://dx.doi.org/10.1016/j.jcp.2010.06.037. [14] T. Gebäck, A. Heintz, A lattice Boltzmann method for the advection–diffusion equation with Neumann boundary conditions, Commun. Comput. Phys. 15 (2) (2014) 487–505, http://dx.doi.org/10.4208/cicp.161112.230713a. [15] L. Li, R. Mei, J.F. Klausner, Boundary conditions for thermal lattice Boltzmann equation method, J. Comput. Phys. 237 (0) (2013) 366–395, http://dx.doi.org/10.1016/j.jcp.2012.11.027. [16] X. He, L.-S. Luo, A priori derivation of the lattice Boltzmann equation, Phys. Rev. E 55 (1997) R6333–R6336, http://dx.doi.org/10.1103/ PhysRevE.55.R6333. [17] X. He, L.-S. Luo, Theory of the lattice Boltzmann method: from the Boltzmann equation to the lattice Boltzmann equation, Phys. Rev. E 56 (1997), http://dx.doi.org/10.1103/PhysRevE.56.6811. [18] Y.H. Qian, D. d’Humières, P. Lallemand, Lattice BGK models for Navier–Stokes equation, Europhys. Lett. 17 (6) (1992) 479–484, http://dx.doi.org/10.1209/ 0295-5075/17/6/001. [19] 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–2) (2005) 179–196, http://dx.doi.org/10.1007/s10955-005-8879-8. [20] B.J. Parker, D.L. Youngs, Two and three dimensional Eulerian simulation of fluid flow with material interfaces, Atomic Weapons Establishment, 1992, AWE preprint. [21] B. Chopard, J.-L. Falcone, J. Latt, The lattice Boltzmann advection–diffusion model revisited, Eur. Phys. J. Spec. Top. 171 (1) (2009) 245–249, http:// dx.doi.org/10.1140/epjst/e2009-01035-5. [22] R. Ammer, M. Markl, V. Jüchter, C. Körner, U. Rüde, Validation experiments for LBM simulations of electron beam melting, Int. J. Mod. Phys. C 25 (2) (2014), http://dx.doi.org/10.1142/S0129183114410095. [23] J. Lätt, Hydrodynamic limits of lattice Boltzmann equations, Ph.D. thesis, Université de Genève, 2007.