Rectangular lattice Boltzmann method using multiple relaxation time collision operator in two and three dimensions

Highlights

• A novel 2D and 3D rectangular grid, MRT lattice Boltzmann implementation.

• Equilibrium moments and matrices included for D2Q9 and D3Q27 lattice implementation.

• Linear stability analysis performed.

• Quasi-incompressible second order multi-scale expansion presented.

• Taylor-Green vortex and turbulent channel flow simulations performed.

Abstract

We present a lattice Boltzmann (lb) method using a rectangular, non-isotropic lattice based on d2q9 and d3q27 velocity sets in two and three dimensions. A second order multi-scale expansion ensures that the scheme correctly reproduces hydrodynamic behaviour. A novel set of basis vectors is introduced in order to allow independent adjustment of eigenvalues corresponding to second order moments as required in order to ensure correct hydrodynamic behaviour using the non-isotropic lattice. Errors are reduced compared to other rectangular grid implementations. Linear perturbation analysis indicates that our scheme has similar stability properties to the isotropic lb method. We investigate the error behaviour of our scheme by performing Taylor-Green vortex flow simulations and comparing our results to simulations using a square grid and also to analytical results. We demonstrate that our scheme is well suited to direct numerical simulation of wall bounded turbulent flows and compare to well known benchmark results.

Introduction

The lattice Boltzmann (lb) method has been increasing in popularity as a means of simulating fluid flow. The explicit and local nature of calculations allows highly efficient parallel implementations and the method has been successfully adapted to massively parallel computing architectures such as graphics processing units [1], [2], [3], [4] as well as achieving high performance on traditional parallel machines [5], [6], [7]. The kinetic nature of the lb method simplifies the implementation of additional physics and the method has enjoyed success in simulating flows involving multiple phases or components, reacting species and complex geometries such as porous media as well as particle laden flows, thermal flows and microfluidics. Thorough reviews of the method and its application to the range of problems mentioned, as well as incompressible and compressible fluid dynamics are presented in Annual Reviews [8], [9] and books [10], [11].

Despite these successes, one persistent drawback has been that the standard lattice Boltzmann method relies on a uniform, square or cube shaped grid coinciding with an isotropic velocity basis so that the advection part of the discrete velocity Boltzmann equation,

$${\partial }_t{f}_{\mathit{i}}(x_{\beta },t)+c_{\mathit{i}\alpha }{\partial }_{\alpha }{f}_{\mathit{i}}(x_{\beta },t)={\Omega }_{\mathit{i}}(x_{\beta },t),$$

may be solved by simple point to point transfer of particles. To date, most progress towards the use of non-uniform grids has been made by retaining the isotropic velocity basis and performing calculations on an unrelated, non-uniform computational grid. The advection of particles can be solved using finite volume methods [6], [12], [13], [14], finite element methods [15], [16], [17], [18], or finite difference methods [19], [20], [21]. Spectral schemes have also been used [22]. Decoupling the grid from the velocity basis even allows the use of velocity basis that does not tile space, such as the 13 velocity, icosahedral basis used by Tamura et al. [21]. Second order or higher difference schemes must be used in order to avoid introducing numerical viscosity, resulting in additional stability concerns compared to the unconditionally stable point to point advection of the normal lb method. Stability can be ensured for finite volume schemes by using flux limiters [6] while others have simply reported reduced stability and the presence of oscillations as drawbacks of their proposed schemes [21].

Standard lb methods use an upwind collision operator, where the collision term is calculated and applied at the beginning of the time step, followed by an independent streaming operation. This decoupling of advection and collision terms introduces an artificial negative viscosity which is related to the grid size and time step and is also referred to as the discrete lattice effect. Thus when the grid is non-uniform, the upwind collision operator may not be used since the discrete lattice effect would vary across the domain and the collision operator must then be solved using a higher order scheme, which must also take into account the interdependent effects of particle advection and collision. Ubertini, Succi and Bella [14] have noted issues in reaching a small viscosity while other studies have not tested the stability at small viscosities [12], [13], [19]. Advanced treatments of the collision operator have been proposed [15], [16], [17], [18], [20] and have demonstrated success with higher Reynolds number unsteady flows, although these methods have thus far only been demonstrated in two dimensions.

The interpolation supplemented lb method [23], [24], [25], [26], [27] is an alternative approach keeping simple point to point advection of particles and an upwind collision operator. Particle locations no longer coincide with grid points after advection so interpolation is used to reconstruct the populations at grid sites. Second order or higher interpolation must be used in order to maintain the second order accuracy of the lb method and results can be improved by using a least squares interpolation [26], [27]. These methods have had success at higher Reynolds numbers, however they require significant additional computational effort and introduce additional stability concerns.

Another alternative to these methods is to use local grid refinement, using smaller cube elements [5], [28], [29], [30]. Grid refinement does not allow the treatment of curved boundaries or the use of non-rectangular elements as the other mentioned methods do however it has been shown to successfully simulate fully turbulent three dimensional flows [5], [31].

All of these methods introduce significant additional computational effort and complexity and some of them suffer from reduced stability. Karlin, Succi and Orszag Karlin et al. [32] were the first to suggested a lattice Boltzmann method on a completely irregular grid but noted that there were significant problems to be solved before it would become usable. Lattice Boltzmann schemes using a non-isotropic but uniform rectangular grid and corresponding velocities basis such that particle populations transfer directly from one node to another, have since been implemented. The generalised lb method, also known as the multiple relaxation time (mrt) method, simplifies implementation of these schemes since it is sufficient to define equilibrium moments rather than an explicit formula for each equilibrium population.

The standard d2q9 mrt lattice can reproduce the correct hydrodynamic moments up to second order. Bouzidi et al. [33] introduced a two dimensional scheme based on the d2q9 lattice using this approach noting stability decreases as the aspect ratio of each cell was increased from one. Zhou [34], [35], Ren et. al. [36] and Hegeler et. al. [37] have also demonstrated similar schemes however the aforementioned models do not correctly model diffusive terms since the third order and higher moments are not hydrodynamically correct and hence viscosity is not isotropic. Adjusting the relaxation rates alone is not sufficient to recover hydrodynamic behaviour using the standard velocity basis.

Peng et. al. [38] have developed a mrt based lb method which correctly reproduces hydrodynamic behaviour. The errors introduced by incorrect rank three moments previously mentioned are corrected by using a quasi-equilibrium collision operator. Peng et. al. perform a multi-scale expansion based on each moment that is used as a basis in order to develop the coefficients that are used in the quasi-equilibrium terms. As a consequence, the analysis must be repeated in order to develop coefficients for each new lattice that may be developed in the future. Results are presented for the d2q9 velocity set using the He and Luo [39] incompressible approximation.

In this paper we present a novel rectangular lattice Boltzmann method including two and three dimensional lattices as shown in Fig. 1. Our method uses a different set of basis vectors compared to typical mrt methods that allows hydrodynamic behaviour to be restored by adjusting relaxation parameters alone without requiring any quasi-equilibrium collision step. Both our work and the work of Peng et. al. [38] correctly eliminates all diffusive errors in the second order multi-scale expansion however our work is also extended to three dimensions. Our theoretical derivation is based on the standard multi-scale expansion typically applied to lb methods rather than the perturbation analysis or moment based analysis used by Peng et. al., which means that our method is applicable to a larger range of lattices without requiring extensive recalculation of coefficients. Our analysis is also novel since it does not use the He and Luo [39] incompressible approximation during multi-scale expansion.

We perform a linear perturbation analysis and find that the maximum stable Mach number decreases if the aspect ratio is increased in the same direction as the mean flow. The scheme retains similar stability compared to uniform grid methods for moderate aspect ratios. This instability is sensitive to the flow being investigated and will be more problematic for configurations with a significant mean velocity. The d3q27 lattice used has a chequerboard invariant that can cause stability problems in some situations [40], [41], [42]. We have successfully used fractional propagation to eliminate this problem. Adjustment of relaxation rates corresponding to non-hydrodynamic moments has been reported as a means of stabilising lb simulations [41], [43] although we have not found it necessary to use this stabilization method in our turbulent calculations.

Our scheme is tested for convergence by comparing to the analytic solutions of Taylor-Green vortex flow with varying Mach number and grid resolution. We perform a direct numerical simulation of turbulent channel flow using both normal and non-isotropic grids and compare turbulence statistics to the results of Moser et al. [44], finding good agreement between lb simulations using isotropic and rectangular grids.

The use of a non-isotropic grid allows improvements in computational efficiency for flow configurations that have larger gradients in one direction. The advantage of the present scheme, compared to other methods capable of using a fully non-uniform grid, is in its simplicity. The scheme does not introduce any new calculation steps or theoretical assumptions compared to the well understood mrt method. Research on the use of non-isotropic lattices is also important because it will advance progress towards the use of a fully non-uniform lattice Boltzmann method [32].