Lattice Boltzmann algorithms without cubic defects in Galilean invariance on standard lattices
Introduction
The Navier–Stokes equations may be derived from the Boltzmann equation that describes dilute monatomic gases. This derivation seeks slowly varying solutions using the Chapman–Enskog expansion [1], [2]. This expansion involves only the lowest few moments of the Maxwell–Boltzmann distribution. The lattice Boltzmann approach to computational fluid dynamics seeks to construct discrete kinetic theories that also yield the Navier–Stokes equations under the same limit. The derivations may proceed exactly in parallel if one replaces the integral moments of continuous kinetic theory with discrete moments or sums over a discrete velocity space.
The continuous Boltzmann equation expresses the evolution of the distribution function for a dilute monatomic gas as The distribution function gives the number density of particles at position x and time t moving with velocity ξ. The left-hand side is a Lagrangian derivative along particle trajectories, while the integral collision operator on the right hand side involves only relative velocities between pairs of particles. The Boltzmann equation is thus invariant under Galilean transformations. This is commonly expressed in kinetic theory by working with the peculiar velocity , the difference between ξ and the local fluid velocity u defined in (6) below.
The lattice Boltzmann approach achieves computational efficiency by restricting the particle velocity ξ to a finite set , and evolves the finite set of functions for rather than . The corresponding discrete Boltzmann equation is a set of N linear, constant coefficient, hyperbolic PDEs, such as (4) below, that is readily discretised in space and time. All nonlinearity is confined to an algebraic collision term, the discrete analogue of . However, the adoption of a finite velocity set inevitably breaks Galilean invariance at the level of the discrete Boltzmann equation. The set of particle velocities typically satisfies the symmetry property that for each i there is an such that , and includes a so-called rest particle with velocity . The frame in which these two properties hold defines a preferred frame.
A much more serious concern is loss of Galilean invariance in the slowly varying limit. The common lattice Boltzmann formulations contain too few degrees of freedom to correctly reproduce the viscous stress in the Navier–Stokes equations. The common formulations use discrete equilibrium distributions that are quadratic polynomials in the fluid velocity u [3], [4], [5] where ρ is the fluid density, and the are weights attached to the different particle velocities . The first three moments of these discrete equilibria coincide with the corresponding moments of the Maxwell–Boltzmann distribution, but the next moment is lacking a term due to the truncation of the discrete equilibria at . This leads to the viscous momentum flux [6] The first term is a Newtonian viscous stress with dynamic viscosity , where θ is the temperature, but the second term is an error term that breaks Galilean invariance. The ratio of the two terms scales as , so the error term is commonly treated as negligible on the assumption that the Mach number Ma is sufficiently small. However, it is an obstacle to the accuracy of the lattice Boltzmann approach at finite Ma, and we show in Section 8 below that it has a detrimental effect on numerical stability.
The natural remedy to this loss of Galilean invariance is to restore the missing terms to the discrete equilibria, but this is not possible when the particle velocities form one of the standard lattices known as D2Q9, D3Q15, D3Q19 or D3Q27 [5], [7]. The particle velocities for the three different 3D lattices are given in Table 1, while the velocities for the D2Q9 lattice arise from projecting any of the 3D lattices onto the xy plane. The dimensionless particle velocity components all lie in the set , so and . The third moment that appears in the calculation of the viscous stress thus cannot be adjusted independently of the first moment.
The next simplest approach is to adopt a larger lattice with more discrete particle velocities [8], [9], [10], although a relatively early attempt [11] to restore the cubic terms on a two-dimensional 17 velocity lattice used incorrect equilibria [12]. Besides adding additional particle velocities whose components lie in a larger integer set such as , one may use higher order Gauss–Hermite quadratures to motivate other particle velocities whose components are not integer multiples of each other [13], [14], [15]. The latter require interpolation of distribution functions between lattice points or some alternative space–time discretisation, while any addition of more degrees of freedom per lattice point increases computational cost, and introduces more scope for grid-scale instabilities of the type analysed in [16].
The free-energy approach of Swift et al. [17] simulates non-ideal fluids by adding terms proportional to the density gradient to the equilibrium momentum flux. These gradients are calculated from finite difference approximations on the lattice. Further finite difference approximations to velocity gradients may also be added to cancel the non-Galilean-invariant defects introduced by the non-ideal terms [18], [19], and also those due to the underlying ideal fluid [20]. However, this again increases complexity and scope for grid-scale instabilities. The standard lattice Boltzmann space/time discretisation (see Section 5) defines effective finite difference stencils that differ from those commonly used [21]. Any inconsistency between the two sets of finite difference stencils may create grid-scale artifacts.
More recently, it has been realised that a partial correction of the third moment is possible on the D2Q9 lattice [12], [22], [23], [24]. This restores Galilean invariance for shear flows aligned with the coordinate axes, which depend only upon the component of the momentum flux, but errors remain in the diagonal components and . In this paper we show that these remaining errors may be largely absorbed using a nonlinear matrix collision operator with relaxation times that depend upon the local fluid velocity u. Matrix collision operators originally appeared as linearisations of lattice gas binary collision operators [25], but have since been developed to improve the stability of lattice Boltzmann algorithms while preserving the isotropy of the viscous stress [26], [27], [28], [29]. In fact, the D3Q13 algorithm requires a matrix collision operator to produce an isotropic viscous stress [30]. However, they have also been used to model the strongly anisotropic stress-strain relation that holds in strongly magnetised plasmas [31]. In this paper we introduce small intentional anisotropies into the matrix collision operator to correct the anisotropies caused by the partial restoration of the cubic terms in the tensor . The resulting algorithm completely restores Galilean invariance at the Navier–Stokes level for flows with uniform density, and reduces the error in the viscous stress from to in the presence of density variations.
Section snippets
The viscous stress in lattice Boltzmann hydrodynamics
Lattice Boltzmann algorithms are space–time discretisations of discrete Boltzmann equations of the form Each is a particle distribution function that propagates with velocity , and the are the components of an collision matrix Ω. The equilibrium distributions are prescribed functions of the moments ρ and u defined in (6) below. For example, the common quadratic polynomial equilibria are given in (2).
The discrete Boltzmann equation (4)
Partially corrected third moment
The natural remedy is to restore the missing cubic terms to so that the discrete third moment coincides with (12) for the Maxwell–Boltzmann distribution. However, the commonly used lattices such as D2Q9, D3Q15, D3Q19, and D3Q27 do not have enough degrees of freedom to specify the whole tensor independently of the lower moments. The particle velocity components for these lattices so, for example, , which implies . It is possible to correct the
Further correction using a matrix collision operator
The viscous stress for a general discrete equilibrium third moment may be written as where is the difference between the third moment of the discrete equilibrium and the corresponding third moment (12) of the Maxwell–Boltzmann distribution.
The partial correction (14) described above leaves only the diagonal components of the tensor non-zero, while the off-diagonal terms and
Implementation
Further discretising the discrete Boltzmann equation (4) in space and time leads to the lattice Boltzmann equation for the transformed distribution functions and the discrete collision matrix [29], [35]. These expressions reduce to those derived by He et al. [36] for the single-relaxation-time collision operator . A recent rederivation of these results using
Linear shear flows at arbitrary angles
Measuring the decay rate of axis-aligned shear flows such as superimposed on a uniform background velocity U is the common benchmark to test loss of Galilean invariance [6], [12], [23], [24]. However, this benchmark only tests the correctness of the off-diagonal component of the viscous momentum flux. We test the correctness of the full viscous momentum flux tensor using shear flows of the form The small deviations to the uniform background velocity U are
Nonlinear simulations
For the more demanding test of an evolving flow with density variations we employ the initial conditions used by Minion and Brown [38] with parameters , in the doubly-periodic domain . The sinusoidal perturbation triggers a Kelvin–Helmholtz instability in the two anti-parallel shear layers, which subsequently roll up into spiral vortices as shown in Fig. 2.
In the absence of cubic and other errors, the lattice
Three dimensions
Keating et al. [24] investigated partial cubic corrections in three dimensions using the equilibria These equilibria give correct values for the off-diagonal components such as on the D3Q27 lattice with the weights given in Table 1. However, contrary to the assertion of Keating et al. [24], these equilibria give incorrect off-diagonal components when used with the D3Q15 or D3Q19 lattices. A direct evaluation using velocities and
Further correction with adjustable bulk viscosity
The more general form of the Navier–Stokes viscous stress is The λ term adds an additional bulk viscosity that enhances the dissipation of sound waves. A bulk viscosity much larger than the shear viscosity is commonly used to improve the stability of lattice Boltzmann algorithms in the weakly compressible regime [27], [28], [48], [49], but also enlarges the non-Galilean contribution to the viscous stress from the trace of the spurious term.
The off-diagonal
Conclusion
The derivation of hydrodynamic equations from discrete kinetic theory via a system of moment equations follows the same route as the derivation from the continuous Boltzmann equation. However, most common lattice Boltzmann algorithms use quadratic polynomials in u for their equilibrium distributions. The third equilibrium moment therefore omits the cubic term that is present in the corresponding equilibrium moment for continuous kinetic theory. This leads to an erroneous term
Acknowledgements
It is a pleasure to acknowledge the hospitality of the Isaac Newton Institute programme Partial Differential Equations in Kinetic Theories, and of the Beijing Computational Science Research Center, where parts of this work were carried out. The computations employed the University of Oxfordʼs Advanced Research Computing facilities, the Emerald GPU cluster of the e-Infrastructure South Centre for Innovation, and the fast Fourier transform library FFTW3 [51]. The authorʼs research was supported
References (51)
- et al.
Nonlinear reactions advected by a flow
Physica A
(1996) - et al.
A Galilean invariant model of the lattice Boltzmann method for multiphase fluid flows using free-energy approach
Comput. Phys. Commun.
(2000) - et al.
Investigation of Galilean invariance of multi-phase lattice Boltzmann methods
Physica A
(2006) - et al.
Simulation of two-dimensional decaying turbulence using the “incompressible” extensions of the lattice Boltzmann method
Comput. Fluids
(2006) Incompressible limits of lattice Boltzmann equations using multiple relaxation times
J. Comput. Phys.
(2003)Lattice Boltzmann formulation for Braginskii magnetohydrodynamics
Comput. Fluids
(2011)An interpretation and derivation of the lattice Boltzmann method using Strang splitting
Comput. Math. Appl.
(2013)- et al.
A novel thermal model of the lattice Boltzmann method in incompressible limit
J. Comput. Phys.
(1998) - et al.
Performance of under-resolved two-dimensional incompressible flow simulations, II
J. Comput. Phys.
(1997) - et al.
Implementation aspects of 3D lattice-BGK: boundaries, accuracy, and a new fast relaxation method
J. Comput. Phys.
(1999)
Lattice Boltzmann method for 3-D flows with curved boundary
J. Comput. Phys.
Direct numerical and large eddy simulation of longitudinal flow along triangular array of rods using the lattice Boltzmann method
Math. Comput. Simul.
Rotational invariance in the three-dimensional lattice Boltzmann method is dependent on the choice of lattice
J. Comput. Phys.
Turbulent jet computations based on MRT and Cascaded Lattice Boltzmann models
Comput. Math. Appl.
The effect of lattice models within the lattice Boltzmann method in the simulation of wall-bounded turbulent flows
J. Comput. Phys.
A numerical method for solving incompressible viscous flow problems
J. Comput. Phys.
The Mathematical Theory of Non-uniform Gases
The Boltzmann Equation and Its Applications
A simple lattice Boltzmann scheme for Navier–Stokes fluid flow
Europhys. Lett.
Lattice Boltzmann model for simulation of magnetohydrodynamics
Phys. Rev. Lett.
Lattice BGK models for the Navier–Stokes equation
Europhys. Lett.
Lattice BGK models for the Navier–Stokes equation: Nonlinear deviation in compressible regimes
Europhys. Lett.
Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation
Phys. Rev. E
Thermal lattice Bhatnagar–Gross–Krook model without nonlinear deviations in macrodynamic equations
Phys. Rev. E
Galilean invariance of lattice Boltzmann models
Europhys. Lett.
Cited by (95)
Lattice Boltzmann method with artificial bulk viscosity using a neural collision operator
2024, Computers and FluidsHPC performance study of different collision models using the Lattice Boltzmann solver Musubi
2023, Computers and FluidsSimulation of the FDA nozzle benchmark: A lattice Boltzmann study
2022, Computer Methods and Programs in BiomedicineA linear stability analysis of compressible hybrid lattice Boltzmann methods
2021, Journal of Computational PhysicsCitation Excerpt :In this context, Feng et al. [37] recently developed a stable ideal gas HLBM on a standard D2Q9 lattice. Using a second-order regularized collision operator [38], along with a body-force term addressing the standard lattice Galilean invariance [39], they were able to simulate natural convection phenomena with large temperature differences. This model was later applied to the simulation of low Mach reactive flows [40].
Recursive finite-difference Lattice Boltzmann schemes
2021, Computers and Mathematics with ApplicationsDirect aeroacoustic simulation with a cumulant Lattice-Boltzmann model
2021, Computers and Fluids