Lattice Boltzmann algorithms without cubic defects in Galilean invariance on standard lattices

Abstract

The vast majority of lattice Boltzmann algorithms produce a non-Galilean invariant viscous stress. This defect arises from the absence of a term in the third moment, the equilibrium heat flow tensor, proportional to the cube of the fluid velocity. This moment cannot be specified independently of the lower moments on the standard lattices such as D2Q9, D3Q15, D3Q19 or D3Q27. A partial correction has recently been demonstrated that restores some of these missing cubic terms on the D2Q9 and D3Q27 tensor product lattices. This correction restores Galilean invariance for shear flows aligned with the coordinate axes, but flows inclined at arbitrary angles may show larger errors than before. These remaining errors are due to the diagonal terms of the equilibrium heat flow tensor, which cannot be corrected on standard lattices. However, the remaining errors may be largely absorbed by introducing a matrix collision operator with velocity-dependent collision rates for the diagonal components of the momentum flux tensor. This completely restores Galilean invariance for flows with uniform density, and in general reduces the magnitude of the defect in Galilean invariance from Mach number cubed to Mach number to the fifth power. The effectiveness of the resulting algorithm is demonstrated by comparisons with the standard and partially corrected lattice Boltzmann algorithms for two- and three-dimensional flows.

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 $f(\mathbf{x},\mathit{\xi },t)$ for a dilute monatomic gas as

$${\partial }_{t}f+\mathit{\xi }\cdot \mathrm{\nabla }f=\mathcal{C}[f,f].$$

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 $\mathbf{c}=\mathit{\xi }-\mathbf{u}$, 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 $\{{\mathit{\xi }}_i:i=0,\dots ,N-1\}$, and evolves the finite set of functions $f_i(\mathbf{x},t)$ for $i=0,\dots ,N-1$ rather than $f(\mathbf{x},\mathit{\xi },t)$. 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 $\mathcal{C}[f,f]$. However, the adoption of a finite velocity set $\{{\mathit{\xi }}_i\}$ 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 $\overline{i}$ such that ${\mathit{\xi }}_i+{\mathit{\xi }}_{\overline{i}}=0$, and includes a so-called rest particle with velocity ${\mathit{\xi }}_0=0$. 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]

$$f_i^{(0)}=w_i\rho (1+3\mathbf{u}\cdot {\mathit{\xi }}_i+\frac{9}{2}{(\mathbf{u}\cdot {\mathit{\xi }}_i)}^2-\frac{3}{2}{|\mathbf{u}|}^2),$$

where ρ is the fluid density, and the $w_i$ are weights attached to the different particle velocities ${\mathit{\xi }}_i$. 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 $\rho \mathbf{u}\mathbf{u}\mathbf{u}$ term due to the truncation of the discrete equilibria at $O({|\mathbf{u}|}^2)$. This leads to the viscous momentum flux [6]

$${\mathit{\Pi }}^{(1)}=-\tau \rho \theta [(\mathrm{\nabla }\mathbf{u})+{(\mathrm{\nabla }\mathbf{u})}^{\mathsf{T}}]+\tau \mathrm{\nabla }\cdot (\rho \mathbf{u}\mathbf{u}\mathbf{u}).$$

The first term is a Newtonian viscous stress with dynamic viscosity $\mu =\tau \rho \theta$, where θ is the temperature, but the second term $\tau \mathrm{\nabla }\cdot (\rho \mathbf{u}\mathbf{u}\mathbf{u})$ is an error term that breaks Galilean invariance. The ratio of the two terms scales as ${|\mathbf{u}|}^2/\theta ={\mathrm{Ma}}^2$, 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 $O({|\mathbf{u}|}^3)$ 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 ${\mathit{\xi }}_i$ 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 ${\xi }_{i\alpha }$ all lie in the set $\{-1,0,1\}$, so ${\xi }_{i\alpha }^3={\xi }_{i\alpha }$ and ${\sum }_i{\xi }_{i\alpha }^3{f}_i={\sum }_i{\xi }_{i\alpha }{f}_i$. 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 $\{-2,-1,0,1,2\}$, 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 ${\Pi }_{xy}$ component of the momentum flux, but errors remain in the diagonal components ${\Pi }_{xx}$ and ${\Pi }_{yy}$. 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 ${\mathsf{Q}}^{(0)}$. 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 $O({\mathrm{Ma}}^3)$ to $O({\mathrm{Ma}}^5)$ in the presence of density variations.