Optimizing free parameters in the D3Q19 Multiple-Relaxation lattice Boltzmann methods to simulate under-resolved turbulent ﬂows

We present a D3Q19 lattice scheme based in MRT with central moments (MRT-CM), where the free parameters of the model are optimized to dissipate under-resolved ﬂow structures with high wavenumbers. In Ch´avez-Modena et al. in Computers & Fluids 172:397-409, 2018 [1], we compared the BGK, MRT-RM and MRT-CM for the D2Q9 lattice scheme using von Neumann analyses and quantiﬁed their numerical properties. Based on this study, we proposed an optimized 2D MRT-CM scheme with enhanced stability for under-resolved ﬂows. Here, we extend this idea to the D3Q19 MRT-CM scheme. As before, we base our optimization for the free parameters, on the k-1% dispersion-error rule, that states that waves with dispersive errors above 1% should be dissipated since they pollute the solution and may cause instabilities. To this aim we increase dissipation in the scheme for waves with dispersive errors above 1%. The resulting optimized scheme is veriﬁed through a von Newmann analysis and validated for the three-dimensional Taylor-Green isotropic turbulent ﬂow. We show how the original D3Q19 MRT-CM (d’Humi`eres version) leads to unrealistic kinetic energy accumulation at high wave numbers, whilst our optimized MRT-CM provides the correct energy dissipated rate, avoiding energy build up

function. The increased complexity with respect to BGK resulted in an improved numerical stability [18].
Despite the enhancement in the stability of the MRT-RM with respect to the BGK collision operator, the MRT-RM can still show instabilities for small fluid viscosities [19]. Consequently, the MRT with central moments (MRT-CM) collision operator was introduced [20]. Premnath and Banerjee [21] showed that the cascaded LBM with central moments is consistent with the Navier-Stokes equations via a Chapman-Enskog expansion, and that in this approach, Galilean invariance is naturally enforced based on the relaxation central moments, see also [22]. Using central moments for the collision operator, it is possible to enhance the numerical stability by increasing the numerical dissipation at high wavenumbers [1].
Von Neumann stability analyses [23] enable the quantification of numerical errors in numerical schemes with periodic boundary conditions. Sterling and Chen [24] were the first to apply this analysis to the LB BGK approach. Then, Lallemand and Luo [18] used it to compare the enhanced stability of the MRT-RM over the BGK approach. Subsequently, Siebert et al. [25] included high order terms into the equilibrium distribution of the D2Q9 model to improve the linear stability of the scheme (as shown depicted Fig. 2 in [25]). Malaspinas [26] proposed a new version of the BGK with improved stability and based on recursive relations and regularization for the LB posed as Hermite series, which has been subsequently validated by Mattila et al. [27] and Coreixas et al. [28]. Later, in our previous work [1], we showed how the MRT-CM is more dissipative at higher wavenumbers compared with BGK and MRT-RM, providing better numerical stability. Also, we explored higher order terms in the fluid velocity, following Malaspinas' work [26] to observe that our proposed approach provides similar results, with comparable stability, and for similar range of Mach numbers.
In addition to quantifying the numerical stability, this technique can also be used to provide insights into dispersive and dissipative errors. Marié et al. [29] compared BGK and MRT-RM collision operators. Dubois et al. [30] studied the numerical stability of the relative velocity (MRT-RM and MRT-CM) D2Q9 schemes with two different set of moments, proposed by Lallemand and Luo [18] and Geier et al. [19]. They concluded that MRT-CM with Geier's moment basis [19] had better stability properties. Recently, Gauthier et al. analyzed the interaction between modes [31] to understand the reason of numerical instabilities.
As mentioned, both the MRT-RM and the MRT-CM relax each hydrodynamic moment with a different relaxation time but all combinations lead to the same macroscopic state through the Chapman-Enskog expansion [32]. The relaxation times that are not fixed by the physics of interest become free parameters that can be optimized to enhance particular numerical aspects. Lallemand and Luo [18] optimized these parameters (for the D2Q9 lattice scheme) maximizing the Galilean invariance of the scheme, while reducing numerical errors (i.e. dispersion and dissipation). Similarly, Xu and Sagaut [7,33] proposed an optimization to minimize dispersion/dissipation errors for the MRT-RM in D2Q9 scheme. A different approach is to adjust the high order relaxation parameters to enhance the stability for under resolved simulations. This approach was used by Ning et al. [34] to improve the stability of a 2D central moment LB method for a lid driven cavity flow problem. However, the authors of that work stated that their results were problem-dependent, and that additional work was required in 3D to optimize the relaxation parameters of the cascaded MRT LBM by means of a Fourier linear analysis. Adam et al. [22] proposed a cascaded LBM for the D3Q19 in the context of non-Newtonian flows. In their approach the free parameter were set to 1 (i.e., equilibrium), although it was already noticed that these could be adjusted independently to improve numerical stability by means of a Fourier linear stability analysis. These ideas were recently exploited by the authors [1], where the D2Q9 MRT-CM scheme was optimized to increase dissipation only for high under-resolved wavenumbers (above the k-1% dispersion-error), leaving low wavenumbers (i.e. well resolved scales) unchanged. This was tested successfully in the double periodic shear layer test.
In this work we extend our previous work to improve the MRT-CM for three-dimensional flows. The optimization strategy used previously for the D2Q9 lattice scheme [1] is now extended to the D3Q19 scheme. Again, instead of minimizing dissipation errors for all wavenumbers, we propose to maintain numerical dissipation for well resolved wavenumbers whilst increasing dissipation for under-resolved wavenumbers. The optimization is inspired in the rule of k-1% dispersion-error presented by Moura et al. [35] in the context of high order numerical methods. They suggested that waves are only accurately resolved if the dispersion error (difference between theoretical and numerical) is below 1%. The wavenumber at which the error becomes 1% was named "k-1% dispersion-error" and lead to the "1% rule". Following this rule, all waves above the k-1% should be dissipated since these are poorly resolved and may pollute the solution. We follow the idea of damping under-resolved waves and apply it to the LBM for the first time in a 3D lattice scheme.
To assess the D3Q19 optimized MRT-CM, we simulate the Taylor-Green vortex (TGV) case [36] that includes starting transitional flow followed by decaying homogeneous turbulence. This case enables the quantification of vortex stretching/pairing processes leading to energy cascading from large to small eddies, allowing the study of the dynamics of transition from laminar to turbulent flow and subsequent turbulent energy decay. This test-case has been widely used to study dissipation errors of numerical schemes, of high order type, e.g. [37,38] and for LB schemes, e.g. [39,40,41].
The remaining of this text is organized as follows. First, in Section 2, we describe the numerical methodology which is divided in two parts: first, the Lattice Boltzmann method with the different collision operators and second the optimization strategy. Then, in Section 3, the results of the optimized approach are tested for the turbulent Taylor Green Vortex case. Finally, in Section 4 conclusions are presented.

Methodology
In this section the numerical methodology used in this work is presented.
First, the Lattice Boltzmann method (LBM) is introduced. Special attention is paid to the definition of the collision operator. Secondly, an optimization method based on linear stability analysis is shown. The optimization aims at maximizing the robustness of the scheme for under-resolver simulations without penalizing its accuracy. The final objective is to improve the performance of the scheme for turbulent flows, therefore a three dimensional LBM scheme, in particular the D3Q19, is considered.

Generalities
The LBM is a numerical technique that provides numerical solutions of the continuous Boltzmann equation. The discrete Boltzmann equation reads: The discrete set of velocities is a vector of Q components represented by the symbol e i . As a consequence, the discrete probability distribution functions (PDFs), f i (x, e i , t), are stored at each lattice node for each time step, t. At each time step, the information stored in the discrete PDFs is streamed through the lattice and collided. The discrete collision operator, ⌦ i , is responsible of computing the post-collision state conserving mass and linear momentum.
The collision operator is of critical importance in the LBM, as it is responsible of the modelization of the physics. The first approach for the collision operator was proposed by McNamara [42], but it was still rather complicated. The high cost of its evaluation precluded its use until Higuera and Jimenez [13] simplified this operator by performing linearization (under the assumption that the discrete PDF, f i , is close to its equilibrium state). Expanding the discrete collision operator, ⌦ i , around the equilibrium state of the discrete PDF, f eq i , leads to a linearized operator: where K ji is known as collision matrix, and was further simplified to obtain the different collision operators. The collision matrix describes how the PDFs relax towards the equilibrium state. It is directly related with the viscosity and Mach number.
The discrete local equilibrium function f eq i is usually computed as a second order Taylor expansion of the Maxwell-Boltzmann distribution, where w i are weighting constants built to preserve isotropy and c s is the speed of sound. The particular values of the weighting constants w i , depend on the discrete set of velocities [43].

Single-relaxation time
The single-relaxation time based on the Bhatnagar-Gross-Krook (BGK) [15] approximation is the most popular approach for the collision operator. In this case, using K BGK ji = (1/⌧ ) ij (assuming ij Kronecker delta notation) the collision operator simplifies to: Only one relaxation time for all PDFs, ⌧ , is considered. The Chapman-Enskog expansion applied to the classical LBM equation (Eq. 1), with Eq. 4 as the discrete collision operator, establishes the relation between the relaxation time, the kinematic shear, ⌫, and bulk viscosities, ⌘, of the macroscopic fluid [44]: where D is the dimension space. The BGK collision operator, is still widely used but has several shortcomings. For example, the kinematic viscosity is conditioned by the relaxation parameter, ⌧ , (see Eq. 5) which causes numerical instabilities for values near ⌧ = 0.5 (with t = 1), therefore complicating the simulation of high Reynolds number flows [16]. It should also be noticed that the bulk viscosity, ⌘, cannot be freely chosen because it is constrained by the kinematic viscosity, ⌫.

Multiple-relaxation time
In the MRT, the relaxation matrix is computed as a product of three matrices, The matrix M accounts for the definition of the moments while the matrix S is defined as a diagonal matrix: where s is a vector with relaxation times for the different moments. The MRT-CM collision operator is obtained defining the collision matrix as K MRT CM ji = M jk (u) 1 S kl M li (u): This notation, first presented by Dubois et al. [45], proposes a more general formulation than Geier's [19], which allows working with the equilibrium PDF, f eq i defined in the BGK model. Under this framework, taking u = 0, one recovers the MRT with raw moments (MRT-RM) presented by d'Humières [17], whilst setting u as the fluid velocity provides Geier's method. Let us note that Dubois et al. formulation is equivalent to Geier's (for u set to the fluid velocity) but the latter uses equilibrium distributions directly derived from the Maxwellian distributions, unlike the former who considered arbitrary velocity fields. Additionally, note that the BGK collision operator can be recovered by setting all the diagonal elements S ii to 1/⌧ . The formulation introduced up to this point is general for any lattice distribution and spatial dimensions. In this work we aim at improving the performance of LBM for under-resolved turbulent flows, therefore a three dimensional LBM is considered. In particular we focus on the popular D3Q19, whose moment matrix, M , can be found in Appendix A.

Optimization method for the MRT-CM
In this section the optimization method to improve the performance of LBM for turbulent flows is introduced. This method, proposed by the authors in a previous work [1], modifies the relaxation times of the collision operator in order to maximize the dissipation of under-resolved wavelengths. The optimization method makes use of linear stability analysis to separate well-resolved and under-resolved wavelengths.

Theoretical modes
The theoretical modes are obtained through an analytic linearization of the NS equations [46]. The resulting hydrodynamic modes are known as shear mode, ! s t , and acoustic modes, ! ± t . The first is related to the kinematic shear viscosity. The second is related to both kinematic shear and bulk viscosities. The analytical expressions for the theoretical hydrodynamic modes read: where k is the wavenumber. For an illustrative purpose Figure 1 depicts dispersion (real part of !, Re(!)) and dissipation (imaginary part of !, Im(!)) behavior for a range of wavenumbers [0, ⇡], following Eq. 8 for three-dimensions The numerical modes can be obtained by means of von Neumann (VN) analysis [24]. Von Neumann analysis splits the PDF into an equilibrium state,f i ,

Numerical modes
The first term in the right hand side,f i , is the global PDF, which does not vary with time and space i.e., it depends only on the average density and velocity.
The second term in the right hand side, f i (x, t), accounts for the fluctuations from equilibrium. The fluctuation is assumed as a sinusoidal wave F i e i(kx !t) with an amplitude F i . Substituting Eq. 9 into Eq. 2, expanding the f eq i by means of Taylor series centred at the global distribution and after rearranging the expression, an eigenvalue problem is obtained, F i = F j G ij , where = e i! t are the eigenvalues of the amplification matrix, G ij , defined as: where the different collision matrices, K ij , are defined in Section 2.1.2. For more details on the VN analysis for the LBM, see [24,1].

Optimization process: Rule of k-1% dispersion-error
The optimization described here has been inspired by the rule of k-1% dispersionerror proposed by Moura et al. [35] in the context of hp-spectral methods [47,48,38]. Moura et al. suggested that wavenumbers with dispersion errors (difference between theoretical and numerical dispersion modes) higher than 1% should be dissipated, since high dispersion errors tend to pollute the solution. In other words, a high dissipative error at high wavenumbers is a positive feature of a numerical method. In a previous work [1], we proposed to optimize the values of the free relaxation times (those affecting moments of order higher than two), to maximize the dissipation of the shear mode in wavenumbers higher than the k-1% dispersion-error wavenumber for the D2Q9 MRT-CM.
The dissipation level is measured as the area defined by the theoretical shear mode and numerical shear mode curves for wavenumbers higher than k-1% dispersion-error wavenumber (green area in Figure 2). The error for the k-1% dispersion-error is calculated with the shear mode (theoretical and numerical) by means of the following expression: As the rule describes, k-1% dispersion-error wavenumber occurs when this error is 1%. The dissipation of acoustic modes has not been taken into account because they barely change with the free parameters.

Optimized MRT-CM
This section details the optimization process for D3Q19 MRT-CM lattice scheme. The optimization aims at maximizing the dissipation of the shear mode in wavenumbers higher than the k-1% dispersion-error wavenumber by modifying the values of the free relaxation times (s 10 18   The optimization starts by evaluating the 1% wavenumbers for d'Humières scheme, shown in Table 2, to define the wavenumber from which we hope to increase dissipation (and defines our objective function). Once determined, the optimal values of s 10 12 and s 13 15 are found by a brute force approach where the parametric space of the free parameters is discretized with s i = 0.01 (i.e. a tolerance of 1%). Note that only values of the relaxation times in the range [0, 2[ are considered. For each parameter combination (Ma, ⌫, s i ) the shear dissipation values above the k-1% dispersion-error wavenumber is obtained. Table   1 shows the optimal values of these parameters (the combination that maximizes the shear dissipation above the k-1%). These values are very different from the free parameters proposed by d'Humières [49] where s 10 12 = 1.98 and s 13 15 = 1.98.
It was previously stated that the group of parameters s 16 18 does not affect the shear dissipation mode, and therefore was kept constant (s 16 18 = 1.4) in the optimization process. In the D3Q19, the first, s 10 12 , and second, s 13 15 groups of parameters relax the third order moments, while the last group, s 16 18 , relaxes the 4th order moments. The relaxation of the 4th order moments affects the energy equation of the compressible Navier-Stokes equations, therefore we do not expect the relaxation times s 16 18 to have a big impact in the dissipation of the shear mode. This behaviour was already reported in our previous work [1], for the D2Q9 scheme.
To confirm this, we compute the shear mode dissipation surface for the different relaxation times. Figure 3 shows the evolution of the shear mode dissipation when s 16 18 and s 10 12 are modified, while the rest of the relaxation parameters remain fixed. As can be seen, the last group of parameters, s 16 18 , does not affect the shear mode dissipation. Note that the maximum dissipation corresponds to the optimized value (s 10 12 = 0.9) for Ma=0.2 and ⌫ = 10 4 . Finally, Table 2 summarizes the k-1% error for d'Humières and our optimized parameters for Mach numbers 0.1  Ma  0.3 and viscosities 10 5  ⌫  10 3 .
Our optimized values provide a slight increase in the wavenumber associated to the 1% dispersion-error, which suggests an increase in dissipation after the k-1%, and a reduction of overall dispersion error. Additionally, since our optimization ensures that the dissipation is maximized for wavenumbers above the k-1%, we can also ensure that above this threshold numerical dissipation will damp dispersive errors.   To illustrate the effect of the optimization, the numerical shear and acous- To provide a wider context, in Appendix B we include a comparison between the numerical modes of the BGK, MRT-RM and MRT-CM D3Q19 lattice schemes.
In addition, some insight on the effect of the relaxation parameters on numerical modes is given in Appendix C where d'Humières [49] and Lallemand's [18] relaxation times are compared for D3Q19 lattice scheme.

Numerical validation
In this section the new optimized scheme is compared with previous approaches found in the literature. The aim of this section is to show that the proposed scheme provides increased stability without penalizing the accuracy.
In order to study the effect of the present optimized MRT-CM on a threedimensional turbulent configuration, the decaying Taylor-Green vortex (TGV) [36] has been simulated. It is a fundamental test case used as prototype for vortex stretching and production of small-scale eddies and therefore allows the study of the dynamics of transition to turbulence. This test-case has been widely used to study the dissipation errors of numerical schemes [37]. In the LB context, the three-dimensional TGV problem has been used to validate different LB approaches in terms of turbulence dissipation [39,40].
In this work, the TGV is simulated with BGK, d'Humières' MRT-CM and

Simulation setup
The TGV problem can be run using a variety of flow and initial conditions.
where the reference density and length, ⇢ 0 and L, are set to one. The equation of state used is p = c 2 s ⇢, therefore the reference pressure is p 0 = 1/3.

Results
To analyze the capabilities of each collision operator the dynamics of a threedimensional decaying vortex has been computed and compared with a DNS simulation performed with a spectral method [37]. In the following, the dissipation rate and the spectrum of the kinetic energy is scrutinized for various under-resolved grids.

Reference results
The TGV evolution is characterized by three main steps visible in the time trace of kinetic energy dissipation rate, " (see Eq. 13). First, the initial laminar state is transitioning to turbulence until the stretched vortex tubes break down into small scales around t/t c = 5. Then the dissipation rate rises to a sharp peak near t/t c = 9 corresponding to the fully turbulent state which is then decaying similarly to an isotropic and homogeneous turbulence. Figure

Kinetic energy dissipation rate
As illustrated before, in order to study the evolution of the fluid, the temporal evolution of the kinetic energy dissipation rate " is calculated: As shown by Eq. 13, the kinetic energy dissipation rate is the temporal derivative of the kinetic energy, E k , which is estimated using the following formula: where the density is ⇢ = ⇢ 0 and V is the volume of the domain.
The results obtained with this definition are compared to the kinetic energy dissipation rate of the reference [37]. Figure 6  the results of the three approaches are overlapped.

Kinetic Energy spectra
The transfer of energy from large scales to small scales is referred to as in- space will be computed as: Subsequently, the energy spectrum function, E(||k||), which represents the contribution to the turbulent kinetic energy over the surfaces of spheres S(||k||) with radius ||k|| = q k 2 x + k 2 y + k 2 z , can be computed as: Continuing with the analysis, in order to understand how the optimiza- ison. The cut-off wavenumber, k cut of f , is the maximum wavenumber for resolved length scales. In particular for the TGV, where x = L N = 2⇧ N , the cut-off wavenumber, assuming a minimum wavelenght of 2 x by Nyquist criterion, is First, notice how the differences between both sets of parameters decrease when the grid resolution increases. Again, with the resolutions N = 128 and N = 256, both sets of parameters reproduce with a good accuracy the dynamics of transition to turbulence. Note that for N = 256 the results of the two approaches are overlapped.
Then, with N = 32 and 64 at t/t c = 6 and 8, before and after the dissipation peak (see Figure 6) the optimized values are slightly less dissipative enabling higher wavenumbers (smaller eddies). However, differences appear at high wavenumbers for t/t c = 14 and 20. We observe that the curves bifurcate after the k 1%opt in all cases. Let us remind the reader that these values approximate the wavenumber at which our scheme differs from d'Humières and adds dissipations for badly resolved waves. Our optimized scheme damps the extra energy build up that appear at high wavenumbers when using d 'Humières' scheme, and allows us to obtain energy slopes that are close to the theoretical value of k 5/3 as shown in Table 3. These numerical slopes, m n are computed at t/t c = 20 in the wavenumber range: k EI < k < k cut of f , where k EI is the demarcation wavenumber between the energy and inertial range [52]. Both schemes are compared with the theoretical slope value, m t = 5/3 ⇡ 1.667.
It may be concluded that our scheme is indeed useful for under-resolved simulations as it prevents energy build up and enables to retrieve the correct cascade of energy and underlying turbulent physics.  Table 3: Numerical slope, mn, and relative error, " = mn m t m t · 100, of the kinetic energy spectra computed for wavenumber range: k EI < k < k cut of f at t/tc = 20.
Finally, our criterion for designing a scheme suitable for LES is not based on physical arguments but rather on numerical aspects. Our approach is based on obtaining the best results with the limited resources available rather than trying to simulate turbulence up to a certain percentage from the Kolmogorov length

Conclusions
In this paper, the free parameter (

Appendix A. Moments matrix for D3Q19
Based on Dubois et al. [45], the following set of moments are proposed: Applying the Gram-Schmidt orthogonalization procedure to Eq. A.1, the matrix M (u) is obtained: As can be seen, the selection of moments has a physical foundation. ⇢ is the density, e is the kinetic energy, ✏ is related to the kinetic energy square, j = ⇢u components correspond to components of momentum, q components correspond to the internal energy components, p components correspond to the symmetric traceless viscous stress tensor, ⇡ is related with the kinetic energy and the viscous stress tensor and m components are the asymmetric third-order moments.
Hence the transformation matrix in raw moments is M (u = 0) with X = e x , Y = e y and Z = e z .
Finally, as it was explained before, S kl is a diagonal matrix. The entries of this diagonal matrix account for the relaxation times, s i , of the different moments. In particular, in the D3Q19 scheme this reads: It is important to notice that, having identified the relations between the kinematic and bulk viscosities with the relaxation times, the higher order moments relaxation times, s 10 18 , remain free parameters. Finally, it should be mentioned that the relaxation values, s i , related with the viscosities must lie between 0 and 2 since the kinematic and bulk viscosity cannot be negative (see Eq. A.6 with t = 1).

Appendix B. BGK, MRT-RM and MRT-CM collision operators
In this appendix, the numerical modes of three different collision operators (BGK, MRT-RM and MRT-CM) are compared. Following our previous work [1] regarding an optimization of the D2Q9 MRT-CM, the bulk viscosity is fixed in BGK and MRT approaches to the same value to permit the comparison between both approaches. So, the relation ⌘ BGK = 2 3 ⌫ is also applied to the MRT collision operator. In particular, the kinematic viscosity is fixed to 10 3 kg ms and consequently the bulk viscosity is fixed to 6.7 · 10 4 kg ms . That means a value of   in Eq. 8. Also as expected, the modes have a similar behaviour as in D2Q9 [1] regarding the dissipation of the shear mode. The MRT-CM approaches presents a higher dissipation rate at high wavenumbers.
Besides, the MRT-RM seems more unstable compared with the D2Q9 lattice scheme at the same fluid conditions (see Figure 3(b) in [1]). These differences between D2Q9 and D3Q19 lattice schemes could come from the different values used from the literature for the free parameters, Lallemand's [18] and d'Humières' [49] values for D2Q9 and D3Q19 respectively. This assumption is clarified in the Appendix C.

Appendix C. Comparison between Lallemand's and d'Humières' parameters:
As mentioned before, initially Lallemand and Luo [18] proposed values for the free parameters of the MRT collision operator for the D2Q9 lattice scheme.
Later, d'Humières' [49] suggested different values for the D3Q19 MRT relaxation times. Both sets of parameters were obtained through an optimization process based on linear stability analysis.