Mathematical Processes 4.1d 4.1.f Number Operation 4.2.f Answers 2019

1. Introduction

Properties such as simplicity, full explicitness and locality of the discrete time-evolution equations have made the lattice Boltzmann method (LBM) a popular tool and serious alternative to classical approaches for the simulation of flows in the limit of vanishing Mach numbers [1,2]. In practice, the classical lattice Boltzmann (LB) formulation is a solver for the isothermal compressible Navier–Stokes (NS) system intended as an approximation to the incompressible NS equations. A proof of this assertion is the presence of acoustic modes propagating at finite speeds, while the incompressible flow assumption, with a divergence-free flow field, enforced through the Poisson equation supposes instant propagation of normal modes and therefore infinite sound speed. As such, the compressible nature of the LBM allows the equation of state to be kept local. Furthermore, boundary closures such, as the bounce-back family of boundary conditions provide for a simple and efficient treatment of solid walls, as one does not have to explicitly define the pressure and/or density there. These properties have allowed this numerical method to be successfully applied to a wide variety of complex weakly compressible flows [3–8].

In recent years, considerable effort has been made in developing suitable LB-based solutions to model fully compressible flows. To extend the LBM to fully compressible flows, models relying on larger stencils were first proposed. The larger stencils allow for a larger number of constraints to be satisfied, and hence a larger number of moments of the discrete distribution function matching their continuous counterparts. Given the larger number of modes (both physical and ghost) present in such models and hence more frequent modal interactions [9], stability issues in the most basic collision models were more pronounced compared to standard low-order stencil solvers. To corroborate this statement one can refer to the usable temperature ranges for different higher-order stencils reported in [10]. In spite of these issues, a number of more advanced collision models such as the entropic [11] and thermal recursive regularized collision models [12] were shown to effectively allow for stable simulations using those extended stencils.

In contrast with previous studies, a number of authors have opted for what one might call a decoupled–recoupled formulation for compressible flows. In such formulations, a thermal LB model is used to solve the compressible NS equations while the balance equation for the second-order conserved variable (energy) is modelled using another solver (either LB-based [13–19] or classical [20–23]). The reduced number of constraints on the flow solver allows for LB models based on standard first-neighbour stencils. However, just as for higher-order stencils, these models are restricted to a limited range of temperatures and velocities due to stability issues and/or errors from higher-order moments. Owing to the insufficient order of the quadrature, first-neighbour stencils fail in recovering the correct third-order moments and therefore lead to errors in the stress tensor at the NS level [14,17,19,22–24]. The Galilean invariance issue associated to those moments, manifesting itself in the form of velocity-dependent dissipation rates for the acoustic (and shear) modes, can be dealt with, to some extent, through appropriate correction terms.

A number of researchers have argued that the reference state (used in the classical approach to construct the discrete equilibrium state) is a major obstacle to the extension of LB solvers to high Mach number flows and flows involving large temperature gradients [25]. They claimed that the classical approach of deriving an equilibrium state via entropy minimization at the continuum level followed by a discretization of this state would result in a discrete equilibrium that inevitably violates the second law of thermodynamics. As a solution to this shortcoming and in line with the general philosophy of the entropic model, the authors proposed a solver relying on a discrete equilibrium with exponential form by solving the discrete entropy minimization problem and finding the corresponding Lagrange multipliers at each node [25]. Such an approach to the construction of stable LB solvers for compressible flows was further extended in [26] by using all of Grad's 13 moments (raw moments instead of temperature-scaled central) in the numerical reconstruction of the discrete equilibrium state. Another approach to guarantee stability of LB solvers based on polynomial equilibria for extreme configurations is to minimize deviations from the reference state by adapting it to the local state of the flow. This rational is based on the fact that the reference state is the only state in which such discrete equilibria exactly satisfy entropy minimization, and hence all deviations in higher-order moments go to zero. In practice a LBM-compatible formulation of this concept was first introduced through shifted stencils in [27]. Later on, a linear stability study along with theoretical analysis of the phase-space discretization procedure confirmed the Galilean invariance of the discrete equilibrium distribution functions (EDFs) and that, when using shifted stencils (with the shifts closely matching the fluid velocity), one could achieve unrestricted Mach numbers [28]. The concept of shifted stencils has also been formulated within the context of a more general framework with temperature- and velocity-dependent discrete stencils [29–31] or through the particle on demand formulation [32]. It is worth noting that such an approach in extending the operation range of discrete kinetic models and reducing higher-order moments errors had previously been used in the context of Eulerian solvers (as opposed to Lagrangian solvers such as LBM). Interested readers are referred to [33–35] for additional details.

Following up on this research, and considering the growing interest for numerical artefacts tied to the reference state in the velocity space discretization process, the aim of the present work is to study the effects of deviation of the fluid temperature from the reference temperature for polynomial discrete equilibria. After a brief introduction of the formalism leading to the thermal discrete equilibrium function, the effect of fluid temperature on errors in the third-order moment at the continuum level are evaluated. Furthermore, the appropriate form of the correction term for the space- and time-discretized LB equations is derived through a Chapman–Enskog (CE) analysis for different orders of the EDF. The von Neumann method is then used to evaluate the stability boundaries of the classical first-neighbour LB method (in terms of fluid temperature). The effects of the order of the EDF, collision model, correction term and discretization order of that term are also studied. Finally, to corroborate the results from the spectral analysis and as a proof of concept, a number of canonical test-cases are considered. These include the decaying shear and acoustic wave cases, freely propagating pressure front, supersonic flow around an NACA0012 aerofoil and the isothermal convected vortex on shifted stencils.

2. Formalism of thermal lattice Boltzmann on standard stencils

Before going into the von Neumann linear analysis and studying the properties of the different models, the mathematical formulations are introduced following these three axes: (a) discretization of the continuous Boltzmann equation and derivation of discrete EDFs on standard stencils, (b) derivation of appropriate correction terms for the third-order moments tensor, and (c) introduction of the different collision models studied here.

(a) Discretization of Boltzmann equation

Phase-space discretization and derivation of appropriate discrete attractors (EDFs) are the most important steps in the development of an LB solver. As such, two different approaches for discretizing phase-space will be detailed here: (a) Gauss–Hermite quadrature-based discretization, and (b) direct moment matching method [36].

(i) Discretization of phase-space: Hermite expansion and Gauss–Hermite quadrature

While the LBM in its initial formulation was developed based on the assumption of an isothermal compressible flow (via a second-order Taylor–McLaurin expansion of the EDF around Ma=0 [37–39]) so as to converge to the incompressible formulation in the limit of vanishing Mach numbers, it has since been extended to thermal flows [40]. The extension to thermal flows, among other methods, can be achieved using the Hermite polynomial expansion approach to construct a discrete kinetic scheme [12,41–43]. This systematic approach to the construction of the discrete equilibrium state also opened the door for higher order EDFs [12,44]. While already documented in a number of articles, the major steps of this derivation will be briefly reviewed here for the sake of clarity. Starting with the Maxwell–Boltzmann equilibrium in the continuum limit in a non-dimensionalized form [45]:

$$f_{\text{MB}}^{(\text{eq})}({\mathit{\xi }}^{\ast },\rho ,{\mathit{u}}^{\ast },{\theta }^{\ast })=\rho {(2\pi {\theta }^{\ast })}^{-D/2}\exp -\frac{{({\mathit{\xi }}^{\ast }-{\mathit{u}}^{\ast })}^2}{2{\theta }^{\ast }},$$ 2.1

where ξ* and u * are the dimensionless particle and fluid velocities, ρ the local density, θ* the non-dimensional temperature and D the physical space dimension, it can be expanded using a class of orthonormal functions, the Hermite polynomials in this case, as [42,43]

$$f_{\text{MB}}^{(\text{eq})}({\mathit{\xi }}^{\ast },\rho ,{\mathit{u}}^{\ast },{\theta }^{\ast })=w({\mathit{\xi }}^{\ast })\sum _{n=0}^{\mathrm{\infty }}\frac{1}{n!}{\mathit{a}}_n^{(\text{eq})}(\rho ,{\mathit{u}}^{\ast },{\theta }^{\ast }):{\mathit{H}}_n({\mathit{\xi }}^{\ast }),$$ 2.2

where ':' is the Frobenius inner product, ${\mathit{H}}_n({\mathit{\xi }}^{\ast })$ is the Hermite polynomial tensor of order n and ${\mathit{a}}_n^{(\text{eq})}$ is the equilibrium coefficient tensor of the corresponding order. This specific choice of orthonormal functions is motivated by the upcoming Gauss–Hermite quadrature (for phase-space discretization) [46] and the fact that the Hermite coefficients are tied to physically meaningful moments of the distribution function [42]. The variables ξ*, u * and θ* are non-dimensionalized using a reference sound speed (based on a reference temperature and molecular mass) $\sqrt{k_B{T}_0/m_0}$ . While not necessary for the expansion (and then later on truncation and quadrature) process, the non-dimensionalization is essential in deriving a stable solver, as it controls errors coming from higher-order components (not correctly recovered by the approximation due to truncation and limited degrees of freedom in the quadrature).

The Hermite polynomials are defined as [47–49]

$${\mathit{H}}_n({\mathit{\xi }}^{\ast })={(-1)}^n\frac{1}{w({\mathit{\xi }}^{\ast })}{\mathbf{\nabla }}_{{\mathit{\xi }}^{\ast }}^{n}w({\mathit{\xi }}^{\ast }),$$ 2.3

where the weight function , w(ξ*), is defined as

$$w({\mathit{\xi }}^{\ast })={(2\pi )}^{-D/2}\exp -\frac{{{\mathit{\xi }}^{\ast }}^2}{2}.$$ 2.4

The coefficients appearing in equation (2.2) are computed as

$${\mathit{a}}_n^{(\text{eq})}(\rho ,{\mathit{u}}^{\ast },{\theta }^{\ast })={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{\mathit{H}}_n({\mathit{\xi }}^{\ast })f_{\text{MB}}^{(\text{eq})}({\mathit{\xi }}^{\ast },\rho ,{\mathit{u}}^{\ast },{\theta }^{\ast })\text{d}{\mathit{\xi }}^{\ast }.$$ 2.5

The first few Hermite polynomials obtained using equation (2.3) are

$$\begin{array}{rl}& {\mathcal{H}}_0=1,\end{array}$$ 2.6a

$$\begin{array}{rl}& {\mathcal{H}}_{i_1}={\xi }_{i_1}^{\ast },\end{array}$$ 2.6b

$$\begin{array}{rl}& {\mathcal{H}}_{i_1{i}_2}={\xi }_{i_1}^{\ast }{\xi }_{i_2}^{\ast }-{\delta }_{i_1{i}_2},\end{array}$$ 2.6c

$$\begin{array}{rl}& {\mathcal{H}}_{i_1{i}_2{i}_3}={\xi }_{i_1}^{\ast }{\xi }_{i_2}^{\ast }{\xi }_{i_3}^{\ast }-{[{\xi }_{i_1}^{\ast }{\delta }_{i_2{i}_3}]}_{\text{cyc}},\end{array}$$ 2.6d

$$\begin{array}{rl}\text{and}& {\mathcal{H}}_{i_1{i}_2{i}_3{i}_4}={\xi }_{i_1}^{\ast }{\xi }_{i_2}^{\ast }{\xi }_{i_3}^{\ast }{\xi }_{i_4}^{\ast }+{[{\delta }_{i_1{i}_2}{\delta }_{i_3{i}_4}]}_{\text{cyc}}-{[{\xi }_{i_3}^{\ast }{\xi }_{i_4}^{\ast }{\delta }_{i_1{i}_2}]}_{\text{cyc}},\end{array}$$ 2.6e

where indexes i ₁, i ₂, … go over all directions in space (i.e. x,y,z), ${\delta }_{i_1{i}_2}$ is the Kronecker delta function and []_cyc designates cyclic permutations without repetition over the indices, while the corresponding thermal coefficients are shown to be (using equation (2.5))

$$\begin{array}{rl}& a_0^{(\text{eq})}=\rho ,\end{array}$$ 2.7a

$$\begin{array}{rl}& a_{i_1}^{(\text{eq})}=\rho u_{i_1}^{\ast },\end{array}$$ 2.7b

$$\begin{array}{rl}& a_{i_1{i}_2}^{(\text{eq})}=\rho u_{i_1}^{\ast }{u}_{i_2}^{\ast }+\rho ({\theta }^{\ast }-1){\delta }_{i_1{i}_2},\end{array}$$ 2.7c

$$\begin{array}{rl}& a_{i_1{i}_2{i}_3}^{(\text{eq})}=\rho u_{i_1}^{\ast }{u}_{i_2}^{\ast }{u}_{i_3}^{\ast }+\rho ({\theta }^{\ast }-1){[u_{i_1}^{\ast }{\delta }_{i_2{i}_3}]}_{\text{cyc}},\end{array}$$ 2.7d

$$\begin{array}{rl}\text{and}& a_{i_1{i}_2{i}_3{i}_4}^{(\text{eq})}=\rho u_{i_1}^{\ast }{u}_{i_2}^{\ast }{u}_{i_3}^{\ast }{u}_{i_4}^{\ast }+\rho {({\theta }^{\ast }-1)}^2{[{\delta }_{i_1{i}_2}{\delta }_{i_3{i}_4}]}_{\text{cyc}}+\rho ({\theta }^{\ast }-1){[u_{i_1}^{\ast }{u}_{i_2}^{\ast }{\delta }_{i_3{i}_4}]}_{\text{cyc}}.\end{array}$$ 2.7e

Being interested in the low-order dynamics of the distribution function (the Euler and NS levels) in the context of the LB method and limited by the order of the quadrature, a truncated approximation of the Hermite expansion can be used [42,43]:

$$f^{(\text{eq},N)}({\mathit{\xi }}^{\ast },\rho ,{\mathit{u}}^{\ast },{\theta }^{\ast })=w({\mathit{\xi }}^{\ast })\sum _{n=0}^N\frac{1}{n!}{\mathit{a}}_n^{(\text{eq})}(\rho ,{\mathit{u}}^{\ast },{\theta }^{\ast }):{\mathit{H}}_n({\mathit{\xi }}^{\ast }).$$ 2.8

An a priori argument for the order of truncation N can be provided via a simple Chapman–Enskog (CE) analysis of the continuous Boltzmann equation (to find highest order moments of the EDF involved at the NS scale). For example, to correctly recover the mass and momentum balance equations at the NS level, one must set N to three. This rational is further backed by the fact that Hermite polynomials are mutually orthonormal. In the remainder of this article it will be shown that a posteriori considerations such as stability enhancement through error control of higher order moments can affect the choice of this truncation [50].

The truncated Hermite polynomial expansion of the distribution function can then be used in a Gauss–Hermite quadrature to discretize phase space (ξ*):

$$\int P^M({\mathit{\xi }}^{\ast })w({\mathit{\xi }}^{\ast })d{\mathit{\xi }}^{\ast }\cong \sum _{\alpha =0}^Q{w}_{\alpha }{P}^M({\mathit{c}}_{\mathit{\alpha }}^{\ast }),$$ 2.9

where c_α * are the discrete abscissae and $w_{\alpha }$ the corresponding weights. An optimal choice of the discrete abscissae (as roots of the corresponding Hermite polynomial) guarantees that equation (2.9) holds exactly for polynomials P(ξ*) up to order 2Q − 1 [41]. In the context of the present derivation, P ^M (ξ*) is defined as

$$P^M({\mathit{\xi }}^{\ast })=\int \frac{{\mathcal{H}}_{\mathit{I}}({\mathit{\xi }}^{\ast })f^N({\mathit{\xi }}^{\ast },\rho ,{\mathit{u}}^{\ast },{\theta }^{\ast })}{w({\mathit{\xi }}^{\ast })}\text{d}{\mathit{\xi }}^{\ast },$$ 2.10

where I +N =M, I being the order of the highest moment to be exactly computed by the quadrature and N the truncation order of the EDF expansion. Following up the rational previously put forward on the order of the truncation, to correctly recover NS level dynamics one would need a minimum quadrature order of Q ≥ (N +M + 1)/2 = 7/2. While a fourth-order quadrature looks like the best choice, it results in the following discrete abscissae (in a one-dimensional configuration) $c_{\alpha ,x}^{\ast }=\pm \sqrt{3\pm \sqrt{6}}$ . These abscissae are not adapted to the LB method as they do not lead to on-lattice propagation of populations. Such stencils can be used in the context of LB models with off-lattice propagation [16,51]. That is why a third-order quadrature (not correctly recovering the third-order moment) is usually employed instead. This third-order quadrature in one dimension leads to abscissae $c_{\alpha ,x}^{\ast }\in \{-\sqrt{3},0,\sqrt{3}\}$ and weights w _α,x  ∈ {1/6, 2/3, 1/6}. Successive application of the previously underlined procedure results in the following tensor product-based EDF in discrete phase-space:

$$f_{\alpha }^{(\text{eq},N)}({\mathit{c}}_{\mathit{\alpha }}^{\mathbf{\ast }},\rho ,{\mathit{u}}^{\ast },{\theta }^{\ast })=w_{\alpha }\sum _{n=0}^N\frac{1}{n!}{\mathit{a}}_n^{(\text{eq})}(\rho ,{\mathit{u}}^{\ast },{\theta }^{\ast }):{\mathcal{H}}_{\alpha ,n}({\mathit{c}}_{\mathit{\alpha }}^{\mathbf{\ast }}),$$ 2.11

where $w_{\alpha }$ is defined as

$$w_{\alpha }=\prod _{i=x,y,z}{w}_{\alpha ,i}.$$ 2.12

The last step in the derivation of the LB scheme is an integration of the partial differential equations (PDEs) along characteristic lines, followed by a change of variable to make the scheme fully explicit, which results in the following discrete time-evolution equations [21,52,53]:

$$f_{\alpha }(\mathit{x}+{\mathit{c}}_{\mathit{\alpha }}{\delta }_t,t+{\delta }_t)=f_{\alpha }(\mathit{x},t)+{\delta }_t{\mathrm{\Omega }}_{\alpha }(\mathit{x},t),$$ 2.13

where ${\mathrm{\Omega }}_{\alpha }(\mathit{x},t)$ is the collision operator detailed in the next section and c_α is non-dimensionalized with δ _x /δ _t . The relaxation coefficient is defined as

$$\tau =\frac{\mu }{p}+\frac{{\delta }_t}{2},$$ 2.14

where μ is the dynamic viscosity. It should be noted that the term δ _t /2 comes from the change of variable in the space-time discretization process. The pressure p is defined as

$$p=\rho \frac{k_{B}T}{m},$$ 2.15

and

$$\frac{{\delta }_t}{{\delta }_x}=\sqrt{\frac{3k_B{T}_0}{m_0}}.$$ 2.16

It should be noted that the factor three appearing inside the square root is specific to the third-order quadrature-based stencil and its tensor products, i.e. D1Q3, D2Q9 and D3Q27, and the D3Q19 stencil. The EDFs and Hermite polynomials and coefficients after discretization in space and time, with standard LBM notation, can be found in appendix A.

(ii) Discretization of phase-space: direct moment-matching

Another approach to discretize phase-space into a set of finite discrete particle velocities is through direct matching of moments [36,54]. One usually starts from the moments of the EDF appearing in the targeted system of macroscopic equations as constraints. Then admissible one-dimensional stencils with enough degrees of freedom and satisfying previously set constraints are selected. For example, for a discrete solver intended for the mass and momentum balance equations up to the NS level, the conditions are the following:

$$\begin{array}{rl}& \sum _{\alpha }{f}_{\alpha }^{(\text{eq})}=\rho ,\end{array}$$ 2.17a

$$\begin{array}{rl}& \sum _{\alpha }{c}_{\alpha }{f}_{\alpha }^{(\text{eq})}=\rho u,\end{array}$$ 2.17b

$$\begin{array}{rl}& \sum _{\alpha }{c}_{\alpha }^2{f}_{\alpha }^{(\text{eq})}=\rho (u^2+\theta ),\end{array}$$ 2.17c

$$\begin{array}{rl}\text{and}& \sum _{\alpha }{c}_{\alpha }^3{f}_{\alpha }^{(\text{eq})}=\rho u(u^2+3\theta ),\end{array}$$ 2.17d

where, for the sake of clarity, both u and $c_{\alpha }$ are non-dimensionalized by δ _x /δ _t . Given the number of free parameters in the D1Q3 stencil, only the first three conditions can be satisfied by the stencil. Solving the following system of linear equations:

$$\mathit{M}{\mathit{f}}_{\mathit{\alpha }}^{\mathbf{(}\text{eq}\mathbf{)}}={\mathbf{\Pi }}^{(\text{eq})},$$ 2.18

where Π ^(eq) are the moments appearing on the right-hand side (r.h.s.) of equations (2.17) and the components of the transform matrix M correspond to the coefficients appearing on the left-hand side (l.h.s.) of equations (2.17), one obtains the following discrete EDFs:

$$f_0^{(\text{eq})}=\rho (1-\frac{{\mathrm{\Pi }}_2^{(\text{eq})}}{\rho })$$ 2.19a

and

$$f_{\pm 1}^{(\text{eq})}=\rho \frac{1}{2}({\mathrm{\Pi }}_2^{(\text{eq})}/\rho \pm u).$$ 2.19b

Using a full tensor product of the D1Q3 stencil, one gets the D2Q9 stencil in two-dimensional with the following discrete EDFs:

$$f_{\alpha }^{(\text{eq})}=\rho \prod _{i=x,y}{(1-{\mathrm{\Pi }}_{ii}^{(\text{eq})}/\rho )}^{1-\mid c_{\alpha ,i}\mid }{\left(\frac{{\mathrm{\Pi }}_{ii}^{(\text{eq})}/\rho +c_{\alpha ,i}{u}_i}{2}\right)}^{\mid c_{\alpha ,i}\mid }.$$ 2.20

Such a discrete equilibrium was used in [14,17,18], among other sources. As for the previous phase-space discretization approach, one can get the discrete (in space in time) time-evolution equation by integrating the system of PDEs along characteristic lines. It is interesting to note that applying the closure relation and matching condition, as in [54], to the obtained discrete EDFs one would recover 1/3 as the non-dimensional (i.e. $k_{B}T/m/({\delta }_x^2/{\delta }_t^2)$ ) reference temperature, and weights corresponding to those derived using the Gausse–Hermite quadrature. Expanding the EDF of equation (2.20), it can be readily shown that it is equivalent to the fourth-order Hermite expansion of the previous subsection [55]. As previously mentioned, the use of the third-order quadrature results in errors at the NS level, regardless of the phase-discretization approach. The different a posteriori corrections proposed in the literature to account for this error will be derived in the next subsection.

(b) Correction of third-order diagonal components of EDF

Given that the EDF of equation (2.20) is equivalent to the fourth-order Hermite-expanded EDF, only the Hermite expansion-based EDFs will be considered for the theoretical analysis of this subsection. A simple CE analysis shows that at the NS level, to match the viscous (non-equilibrium) stress tensor for the continuous Boltzmann equation, moments of orders two and three of the EDF must be exactly recovered. Integrating in phase space, the following continuous second- and third-order moments are in fact recovered [42]:

$${\mathrm{\Pi }}_{i_1{i}_2}^{\text{MB}}=\rho u_{i_1}{u}_{i_2}+\rho c_s^2{\delta }_{i_1,i_2}{\theta }^{\ast }$$ 2.21a

and

$${\mathrm{\Pi }}_{i_1{i}_2{i}_3}^{\text{MB}}=\rho u_{i_1}{u}_{i_2}{u}_{i_3}+\rho c_s^2{[u_{i_1}{\delta }_{i_2,i_3}{\theta }^{\ast }]}_{\text{cyc}},$$ 2.21b

while one gets the following moments with the second- and third-order discrete EDFs:

$${\mathrm{\Pi }}_{i_1{i}_2}^{(N=2)}=\rho u_{i_1}{u}_{i_2}+\rho c_s^2{\delta }_{i_1,i_2}{\theta }^{\ast }$$ 2.22a

and

$${\mathrm{\Pi }}_{i_1{i}_2{i}_3}^{(N=2)}=\rho c_s^2{[u_{i_1}{\delta }_{i_2,i_3}]}_{\text{cyc}},$$ 2.22b

and

$${\mathrm{\Pi }}_{i_1{i}_2}^{(N=3)}=\rho u_{i_1}{u}_{i_2}+\rho c_s^2{\delta }_{i_1,i_2}{\theta }^{\ast }$$ 2.23a

and

$${\mathrm{\Pi }}_{i_1{i}_2{i}_3}^{(N=3)}=\rho c_s^2{\delta }_{i_1{i}_2{i}_3}{[u_{i_1}{\delta }_{i_2,i_3}]}_{\text{cyc}}+\rho (1-{\delta }_{i_1{i}_2{i}_3})\{u_{i_1}{u}_{i_2}{u}_{i_3}+c_s^2{[u_{i_1}{\delta }_{i_2,i_3}{\theta }^{\ast }]}_{\text{cyc}}\}.$$ 2.23b

The discrete EDF of equation (2.20) recovers the same second- and third-order moments as the third-order Hermite expansion. To further illustrate the errors in the moments of the discrete EDFs, they are plotted in figure 1 as a function of the non-dimensional temperature. Given the focus of the present manuscript on thermal effects, the velocity-dependence of these errors is not shown here. Interested readers are referred to [28] for a detailed study of this specific aspect. As observed in these figures, the error goes to zero at θ* = 1 for the Hermite-expanded EDFs, while deviations from this state result in very pronounced errors in third-order moments. The EDF of equation (2.20) exhibits the same behaviour and the error goes to zero at θ = 1/3, which is equivalent to θ* = 1.

Figure 1. Error in the diagonal components of the third-order moments tensor using (a) second-order, shown in red, and third-order, shown in blue, Hermite-expanded EDFs and (b) EDF of equation (2.20). (Online version in colour.)

• Download figure
• Open in new tab
• Download PowerPoint

To better put forward the shortcoming of the first-order stencil in recovering the NS level terms, let us perform now a brief CE analysis. Introducing a second-order Taylor–McLaurin expansion, the discrete time-evolution equation changes into

$$D_t{f}_{\alpha }+\frac{1}{2}{D_t}^2{f}_{\alpha }=\frac{{\delta }_t}{\tau }(f_{\alpha }^{(\text{eq})}-f_{\alpha })+{\mathrm{\Psi }}_{\alpha }+\mathcal{O}({\delta }_t^3),$$ 2.24

where $D_t={\mathrm{\partial }}_t+{\mathit{c}}_{\mathit{\alpha }}\cdot \mathbf{\nabla }$ and ${\mathrm{\Psi }}_{\alpha }$ is the correction term for the third-order moments of the EDF, to be determined via the CE analysis of this subsection. Two different approaches for correcting the third-order moments tensor have been proposed in the literature. Both will be derived here. As usual with such developments a multi-scale expansion is introduced:

$$\begin{array}{rl}& f_{\alpha }=f_{\alpha }^{(0)}+\epsilon f_{\alpha }^{(1)}+{\epsilon }^2{f}_{\alpha }^{(2)}+\mathcal{O}({\epsilon }^3),\end{array}$$ 2.25a

$$\begin{array}{rl}& {\mathrm{\partial }}_t=\epsilon {\mathrm{\partial }}_t^{(1)}+{\epsilon }^2{\mathrm{\partial }}_t^{(2)}+\mathcal{O}({\epsilon }^3),\end{array}$$ 2.25b

$$\begin{array}{rl}& \mathbf{\nabla }=\epsilon {\mathbf{\nabla }}^{(1)}\end{array}$$ 2.25c

$$\begin{array}{rl}\text{and}& {\mathrm{\Psi }}_{\alpha }=\epsilon {\mathrm{\Psi }}_{\alpha }^{(1)}.\end{array}$$ 2.25d

Sorting terms of different orders in ε into separate equations and taking moments of orders zero and one to get mass and momentum conservation, at order two in ε (NS level) the following momentum conservation equation is recovered:

$$\begin{array}{rl}& {\mathrm{\partial }}_t^{(2)}\rho \mathit{u}+{\mathbf{\nabla }}^{(1)}\cdot (\frac{1}{2}-\frac{\tau }{{\delta }_t})\{{\mathrm{\partial }}_t^{(1)}{\mathbf{\Pi }}_2^{(\text{eq},N)}+{\mathbf{\nabla }}^{(1)}\cdot {\mathbf{\Pi }}_3^{(\text{eq},N)}\}\\ & +{\mathbf{\nabla }}^{(1)}\cdot \left(\frac{\tau }{{\delta }_t}\right)(\sum _{\alpha }{\mathit{c}}_{\mathit{\alpha }}\otimes {\mathit{c}}_{\mathit{\alpha }}{\mathrm{\Psi }}_{\alpha }^{(1)})=0.\end{array}$$ 2.26

As such, for the stress tensor to be correctly recovered at this scale one must have

$${\mathrm{\Psi }}_{\alpha }=(1-\frac{{\delta }_t}{2\tau })\frac{w_{\alpha }}{2c_s^4}(\mathbf{\nabla }\otimes {\mathcal{H}}_{\alpha ,2}):{\mathit{\delta }\mathbf{\Pi }}_3^{(\text{eq},N)},$$ 2.27

where ${\mathit{\delta }\mathbf{\Pi }}_3^{(\text{eq},N)}$ designates the deviation of the discrete EDF moment from its continuous counterpart. The form of the correction, as seen here, involves a first-order space-derivative of a combination of macroscopic variables and as such provides an a posteriori justification to the multi-scale expansion introduced in equation (2.25). An interesting point to note here is the presence of a coefficient (1 − (δ _t /2τ)) similar to Guo's forcing scheme [56], in front of the correction term (it does not appear in the space-continuous development) [19,22,23]. Omission of this factor results in pronounced deviations regarding acoustic modes dispersion, especially as the relaxation coefficient τ/δ _t gets away from 1.

It is also worth noting that the treatment for third- and fourth-order EDFs differs from that at second-order. Regarding third- and fourth-order EDFs, one obtains

$${\mathrm{\Psi }}_{\alpha }^{(N>2)}=(1-\frac{{\delta }_t}{2\tau })\frac{w_{\alpha }}{2c_s^4}[{\mathcal{H}}_{\alpha ,xx}{\mathrm{\partial }}_x\delta {\mathrm{\Pi }}_{xxx}^{(\text{eq},N)}+{\mathcal{H}}_{\alpha ,yy}{\mathrm{\partial }}_y\delta {\mathrm{\Pi }}_{yyy}^{(\text{eq},N)}].$$ 2.28

Instead, for the second-order EDF, additional correction terms are required, or one fails to recover correctly deviatoric components of the third-order moments:

$$\begin{array}{rl}{\mathrm{\Psi }}_{\alpha }^{(N=2)}& ={\mathrm{\Psi }}_{\alpha }^{(N>2)}\\ & +(1-\frac{{\delta }_t}{2\tau })\frac{w_{\alpha }}{c_s^4}{\mathcal{H}}_{\alpha ,xy}[{\mathrm{\partial }}_x(\delta {\mathrm{\Pi }}_{xxy}^{(\text{eq},N)}+\delta {\mathrm{\Pi }}_{xyy}^{(\text{eq},N)})+{\mathrm{\partial }}_y(\delta {\mathrm{\Pi }}_{xxy}^{(\text{eq},N)}+\delta {\mathrm{\Pi }}_{xyy}^{(\text{eq},N)})].\end{array}$$ 2.29

As opposed to the approach taken previously, one can also directly introduce the correction term at order ε² as proposed in [17,18]. In practice, this means that the correction would involve a Laplacian and as such be expanded as

$${\mathrm{\Psi }}_{\alpha }^{\mathrm{\prime }}={\epsilon }^2{{\mathrm{\Psi }}_{\alpha }^{\mathrm{\prime }}}^{(2)}.$$ 2.30

Re-writing the momentum conservation equation at the NS level with this new correction term:

$${\mathrm{\partial }}_t^{(2)}\rho \mathit{u}+{\mathbf{\nabla }}^{(1)}\cdot (\frac{1}{2}-\frac{\tau }{{\delta }_t})\{{\mathrm{\partial }}_t^{(1)}{\mathbf{\Pi }}_2^{(\text{eq},N)}+{\mathbf{\nabla }}^{(1)}\cdot {\mathbf{\Pi }}_3^{(\text{eq},N)}\}-\sum _{\alpha }{\mathit{c}}_{\mathit{\alpha }}{{\mathrm{\Psi }}_{\alpha }^{\mathrm{\prime }}}^{(2)}=0,$$ 2.31

one gets the following restrictions on the correction term:

$$\sum _{\alpha }{{\mathrm{\Psi }}_{\alpha }^{\mathrm{\prime }}}^{(2)}=0$$ 2.32

and

$$\sum _{\alpha }{\mathit{c}}_{\alpha }{{\mathrm{\Psi }}_{\alpha }^{\mathrm{\prime }}}^{(2)}={\mathbf{\nabla }}^{(1)}\cdot (\frac{\mu }{p}{\mathbf{\nabla }}^{(1)}\cdot \delta {\mathbf{\Pi }}_3^{(\text{eq},N)}).$$ 2.33

The correction term using the second approach can therefore be defined as

$${\mathrm{\Psi }}_{\alpha }^{\mathrm{\prime }}=\frac{w_{\alpha }}{c_s^2}{\mathit{c}}_{\alpha }\cdot \mathbf{\nabla }\cdot (\frac{\mu }{p}\mathbf{\nabla }\cdot \delta {\mathbf{\Pi }}_3^{(\text{eq},N)}).$$ 2.34

The effectiveness of both of these approaches in restoring the Galilean invariance of the shear and acoustic modes dissipation rates will be analysed in the next sections.

(c) Collision operators

The ability of four different classes of collision models to handle deviations from the reference temperature will be studied in this manuscript: (i) single relaxation time (SRT), (ii) regularized, (iii) multiple relaxation time (MRT) in central Hermite moments space, and (iv) MRT in temperature-scaled central Hermite moments space. As such, they will be briefly introduced.

(i) Single relaxation time

The SRT is the most basic form of the Bhatnagar–Gross–Krook (BGK) collision operator [57], relaxing all components of the distribution function at the same rate. It is defined as

$${\mathrm{\Omega }}_{\alpha }^{(\text{SRT})}=\frac{1}{\tau }(f_{\alpha }^{(\text{eq},N)}-f_{\alpha }).$$ 2.35

This collision operator, while efficient and easy to implement becomes practically unusable when the non-dimensional relaxation coefficient, τ/δ _t , gets close to 0.5 [2]. As such, to overcome this issue and extend the usability range of the solver a number of more complex collision models have been developed.

(ii) Regularized

As readily shown through a CE analysis, at the NS level moments of the EDF up to order three are involved, while the non-equilibrium part only affects the dynamics through its second-order moment. One could, in principle, filter out all higher-order contributions to get better accuracy and stability. The original regularized collision model was developed and proposed with this argument in mind [58] even though other rationals could lead to the same operator, e.g. [59,60]. The regularized time evolution equation is

$$f_{\alpha }(\mathit{x}+{\mathit{c}}_{\mathit{\alpha }}{\delta }_t,t+{\delta }_t)=f_{\alpha }^{(\text{eq},N)}(\mathit{x},t)+(1-\frac{{\delta }_t}{\tau })f_{\alpha }^{(1,N)}(\mathit{x},t),$$ 2.36

where $f_{\alpha }^{(1)}$ is the first-order contribution (in the CE expansion) of the non-equilibrium part, which similarly to the EDF can be expanded in terms of Hermite polynomials:

$$f_{\alpha }^{(1,N)}=w_{\alpha }\sum _{n=2}^N\frac{1}{n!c_s^{2n}}{\mathit{a}}_{\mathit{n}}^{\mathbf{(}\mathbf{1}\mathbf{)}}:{\mathcal{H}}_{\mathit{\alpha }\mathbf{,}\mathit{n}}.$$ 2.37

In the original regularized model [58], N = 2, and ${\mathit{a}}_{\mathbf{2}}^{\mathbf{(}\mathbf{1}\mathbf{)}}$ is approximated by projection as

$${\mathit{a}}_{\mathbf{2}}^{\mathbf{(}\mathbf{1}\mathbf{)}}\approx {\mathit{a}}_{\mathbf{2}}^{\mathbf{(}\mathbf{n}\mathbf{e}\mathbf{q}\mathbf{)}}=\sum _{\alpha }{\mathcal{H}}_{\alpha ,n}(f_{\alpha }-f_{\alpha }^{(\text{eq},N)}).$$ 2.38

This approach is limited to second-order EDFs, and has limited effects on the stability domain for vanishing viscosities [28,50]. Higher-order regularization in combination with closures extracted from a CE analysis have been shown to drastically improve the performances of the regularized collision operator [12,61,62]. Using the CE expansion to reconstruct the first-order non-equilibrium contribution [12,63]:

$$\begin{array}{rl}& a_{xyy}^{(1)}=u_x{a}_{yy}^{(1)}+2u_y{a}_{xy}^{(1)}-\tau \rho c_s^4{\theta }^{\ast }{\mathrm{\partial }}_x{\theta }^{\ast },\end{array}$$ 2.39a

$$\begin{array}{rl}& a_{xxy}^{(1)}=u_y{a}_{xx}^{(1)}+2u_x{a}_{xy}^{(1)}-\tau \rho c_s^4{\theta }^{\ast }{\mathrm{\partial }}_y{\theta }^{\ast },\\ \text{and}& a_{xxyy}^{(1)}=u_x^2{a}_{yy}^{(1)}+u_y^2{a}_{xx}^{(1)}+4u_x{u}_y{a}_{xy}^{(1)}+c_s^2({\theta }^{\ast }-1)(a_{yy}^{(1)}+a_{xx}^{(1)})\end{array}$$ 2.39b

$$\begin{array}{rl}& -\tau \rho c_s^4{\theta }^{\ast }(u_x{\mathrm{\partial }}_x{\theta }^{\ast }+u_y{\mathrm{\partial }}_y{\theta }^{\ast }).\end{array}$$ 2.39c

It must be noted that all terms involving space derivatives of the temperature will be omitted in the von Neumann analysis as the temperature is supposed uniform in space. All three orders, i.e. two, three and four, of the regularization will be considered in the present study. At order two, the reconstruction is based on projection, while at orders three and four the recursive regularization approach is used.

The regularized family of collision operators can also be seen as pertaining to the MRT family of collision operators, as they rely on relaxing physical and ghost modes at different rates. The projection-based regularization approach can be readily shown to pertain to the MRT family of operators in raw Hermite moments space [50,64]. The recursive regularization approach, on the other hand, for isothermal flows has also been shown to be an MRT model in central Hermite moments space [65]. To better clarify this point, and further corroborate this equivalence for thermal flows, two classes of MRT collision models will be considered in this paper, and therefore detailed in the next part.

(iii) Multiple relaxation time collision model in central moments and temperature-scaled central

moments space)

As an extension to the generalized BGK collision model proposed in [66,67], the cascaded operator was first introduced in [68]. While relying on the same general philosophy of applying the collision operator in moments space, the central moments space was introduced as a way to further de-alias relaxation of different modes. The collision operator is written as [69]

$${\mathrm{\Omega }}_{\alpha }^{\text{(CM-MRT)}}={\mathcal{T}}^{-1}\mathit{S}\mathcal{T}(f_{\alpha }^{(\text{eq})}-f_{\alpha }),$$ 2.40

where $\mathcal{T}$ is the central moments transform matrix defined as

$$\tilde{\mathbf{\Pi }}=\mathcal{T}{\mathit{f}}_{\mathit{\alpha }},$$ 2.41

and $\tilde{\mathbf{\Pi }}$ are moments of the distribution function in central moments space:

$${\tilde{\mathrm{\Pi }}}_{\mathit{i}}=\sum _{\alpha }\prod _{j\in \mathit{i}}(c_{\alpha ,j}-u_j)f_{\alpha },$$ 2.42

with i  = {i ₁, i ₂…i _N }. While one can chose among a wide variety of moments bases, for the sake of clarity and to better illustrate the relationship between the regularized collision operators and the MRT collision model, central Hermite polynomials will be chosen as the moments basis here:

$$\begin{array}{rl}& {\mathcal{T}}_{1,\beta }={\tilde{\mathcal{H}}}_{\beta ,0}=1,\end{array}$$ 2.43a

$$\begin{array}{rl}& {\mathcal{T}}_{2,\beta }={\tilde{\mathcal{H}}}_{\beta ,x}=c_{\beta ,x}-u_x,\end{array}$$ 2.43b

$$\begin{array}{rl}& {\mathcal{T}}_{3,\beta }={\tilde{\mathcal{H}}}_{\beta ,y}=c_{\beta ,y}-u_y,\end{array}$$ 2.43c

$$\begin{array}{rl}& {\mathcal{T}}_{4,\beta }={\tilde{\mathcal{H}}}_{\beta ,xy}=(c_{\beta ,x}-u_x)(c_{\beta ,y}-u_y),\end{array}$$ 2.43d

$$\begin{array}{rl}& {\mathcal{T}}_{5,\beta }={\tilde{\mathcal{H}}}_{\beta ,xx}={(c_{\beta ,x}-u_x)}^2-c_s^2,\end{array}$$ 2.43e

$$\begin{array}{rl}& {\mathcal{T}}_{6,\beta }={\tilde{\mathcal{H}}}_{\beta ,yy}={(c_{\beta ,y}-u_y)}^2-c_s^2,\end{array}$$ 2.43f

$$\begin{array}{rl}& {\mathcal{T}}_{7,\beta }={\tilde{\mathcal{H}}}_{\beta ,xyy}=(c_{\beta ,x}-u_x)[{(c_{\beta ,y}-u_y)}^2-c_s^2],\end{array}$$ 2.43g

$$\begin{array}{rl}& {\mathcal{T}}_{8,\beta }={\tilde{\mathcal{H}}}_{\beta ,xxy}=(c_{\beta ,y}-u_y)[{(c_{\beta ,x}-u_x)}^2-c_s^2]\end{array}$$ 2.43h

$$\begin{array}{rl}\text{and}& {\mathcal{T}}_{9,\beta }={\tilde{\mathcal{H}}}_{\beta ,xxyy}=[{(c_{\beta ,x}-u_x)}^2-c_s^2][{(c_{\beta ,y}-u_y)}^2-c_s^2].\end{array}$$ 2.43i

This choice of moments basis results in the following equilibrium moments:

$${\tilde{\mathrm{\Pi }}}^{(\text{eq})}=\{\rho ,0,0,0,\rho c_s^2({\theta }^{\ast }-1),\rho c_s^2({\theta }^{\ast }-1),0,0,\rho c_s^4{({\theta }^{\ast }-1)}^2\}.$$ 2.44

It should be noted that the full fourth-order Hermite expansion has been used for the EDF here.

To adapt this collision model to thermal flows, it can be further extended to temperature-scaled central moments [70] as originally proposed by Grad [47]. This is achieved via the following collision operator:

$${\mathrm{\Omega }}_{\alpha }^{\text{(sCM-MRT)}}={\mathcal{T}}_{{\theta }^{\ast }}^{-1}\mathit{S}{\mathcal{T}}_{{\theta }^{\ast }}(f_{\alpha }^{(\text{eq})}-f_{\alpha }),$$ 2.45

with

$${\tilde{\mathbf{\Pi }}}_{{\theta }^{\ast }}={\mathcal{T}}_{{\theta }^{\ast }}{\mathit{f}}_{\mathit{\alpha }},$$ 2.46

and

$${\tilde{\mathrm{\Pi }}}_{{\theta }^{\ast },\mathit{i}}=\sum _{\alpha }\prod _{j\in \mathit{i}}\left(\frac{c_{\alpha ,j}-u_j}{\sqrt{{\theta }^{\ast }}}\right)f_{\alpha }.$$ 2.47

As for the previous model, choosing temperature-scaled central Hermite polynomials one gets the following equilibrium moments:

$${\tilde{\mathrm{\Pi }}}^{(\text{eq})}=\{\rho ,0,0,0,0,0,0,0,0\}.$$ 2.48

The performances of both these collision operators will be studied through the von Neumann analysis in the next section.

3. Linear von Neumann analysis

To better understand the effect of temperature deviations from the reference state on stability, dispersion and dissipation the von Neumann method is used.

(a) Formalism

For nonlinear systems of equations, this approach consists of introducing a perturbation into the linearized version of the considered system of equations in a fully periodic domain, then following its time evolution [71]. For nonlinear systems, under the assumption of a linear regime (weak perturbation) the field can be expanded as a combination of a uniform reference state and a first-order deviation:

$$f_{\alpha }\approx {\overline{f}}_{\alpha }+f_{\alpha }^{\mathrm{\prime }}.$$ 3.1

Introducing this first-order Taylor–McLaurin expansion into the discrete-time evolution equation one can recover the discrete linearized equations for the perturbation. Given that the topic has been thoroughly treated in the literature [9,28,50,72,73], details and derivation of the final equations will be omitted. For example for the SRT collision operator the following equation is recovered:

$$f_{\alpha }^{\mathrm{\prime }}(\mathit{x}+{\mathit{c}}_{\mathit{\alpha }}{\delta }_t,t+{\delta }_t)={\delta }_{\alpha \beta }(1-\frac{{\delta }_t}{\tau })+\frac{{\delta }_t}{\tau }{J}_{\alpha \beta }(f_{\alpha }^{(\text{eq},N)}),$$ 3.2

while the linearized time evolution equation for the regularized model is expressed as

$$f_{\alpha }^{\mathrm{\prime }}(\mathit{x}+{\mathit{c}}_{\mathit{\alpha }}{\delta }_t,t+{\delta }_t)=J_{\alpha \beta }(f_{\alpha }^{(\text{eq},N)})+(1-\frac{{\delta }_t}{\tau })J_{\alpha \beta }(f_{\alpha }^{(1,N)}),$$ 3.3

where $J_{\alpha \beta }(f_{\alpha }^{(\text{eq},N)})$ and $J_{\alpha \beta }(f_{\alpha }^{(1,N)})$ are the Jacobians of the equilibrium and first-order non-equilibrium parts of the distribution function. These Jacobians can be computed as

$$J_{\alpha \beta }(f_{\alpha }^{(\text{eq},N)})=w_{\alpha }\sum _{n=0}^N\frac{1}{n!c_s^{2n}}{\mathcal{H}}_{\mathit{\alpha }\mathbf{,}\mathit{n}}:{\mathrm{\partial }}_{f_{\beta }}{\mathit{a}}_{\mathit{n}}^{\mathbf{(}\text{eq}\mathbf{,}\mathit{N}\mathbf{)}}{\mid }_{{\overline{f}}_{\beta }}$$ 3.4

and

$$J_{\alpha \beta }(f_{\alpha }^{(1,N)})=w_{\alpha }\sum _{n=2}^N\frac{1}{n!c_s^{2n}}{\mathcal{H}}_{\mathit{\alpha }\mathbf{,}\mathit{n}}:{\mathrm{\partial }}_{f_{\beta }}{\mathit{a}}_{\mathit{n}}^{\mathbf{(}\mathbf{1}\mathbf{,}\mathit{N}\mathbf{)}}{\mid }_{{\overline{f}}_{\beta }}.$$ 3.5

The correction terms are also linearized via their Jacobians as

$$J_{\alpha \beta }({\mathrm{\Psi }}_{\alpha })=(1-\frac{{\delta }_t}{2\tau })\frac{w_{\alpha }}{2c_s^4}({\mathcal{H}}_{\alpha ,2}\otimes \mathbf{\nabla }):{\mathrm{\partial }}_{f_{\beta }}{\mathit{\delta }\mathbf{\Pi }}_3^{(\text{eq},N)}{\mid }_{{\overline{f}}_{\beta }}$$ 3.6

and

$$J_{\alpha \beta }({\mathrm{\Psi }}_{\alpha }^{\mathrm{\prime }})=\frac{w_{\alpha }}{c_s^2}{\mathit{c}}_{\alpha }\cdot (\mathbf{\nabla }\cdot \frac{\mu }{p}\mathbf{\nabla })\cdot {\mathrm{\partial }}_{f_{\beta }}{\mathit{\delta }\mathbf{\Pi }}_3^{(\text{eq},N)}{\mid }_{{\overline{f}}_{\beta }}.$$ 3.7

For a detailed derivation of these Jacobians, interested readers are referred, e.g. to [28,50].

Introducing the standing waves into the linearized discrete-time evolution equation:

$$f_{\alpha }^{\mathrm{\prime }}=F_{\alpha }^{\mathrm{\prime }}\exp (i\omega t-\mathit{k}\cdot \mathit{x}),$$ 3.8

one obtains a system of equations, through which the time-amplification factor of the perturbation ω can be obtained for different wavenumbers k . In order for the system to be linearly stable for the chosen set of parameters, i.e. velocity, temperature and non-dimensional viscosity, the real component of ω must remain negative for all possible values of k .

For all stability domains presented in this section, the two-dimensional wavenumber space k  ∈ [ −ππ] × [0π] is discretized using a resolution of 100 points in each direction, i.e. $\delta k_x,\delta k_y=\frac{\pi }{100}$ . Furthermore, the orientation of the velocity vector was also taken into account by considering different angles by reference to the x-axis, i.e. ϕ ∈ [ −π/2, π/2] and δϕ =π/20. The stability of the different schemes has been evaluated for a set of parameters $({\theta }^{\ast },\nu {\delta }_t/{{\delta }_x}^2)$ with ${\theta }^{\ast }\in \{0.001,0.05,0.1,0.2,\dots ,2.3\}$ and $\nu {\delta }_t/{{\delta }_x}^2=p\times {10}^{-q}$ with p ∈ {1, 5} and q ∈ {0, 1, 2, 3, 4, 5, 6}.

(b) Effect of EDF order

As shown in previous sections, in practice and on standard first-order stencils, only the addition of off-diagonal components of the third-order Hermite coefficients affects the NS-level dynamics. The addition of fourth-order components of the EDF, however, while not affecting the NS-level terms can improve the stability domain of the scheme [50,62]. To that end, the linear stability domains of the SRT collision operator with different orders of the EDF were evaluated. The results (maximum achievable non-dimensional velocities and Mach numbers) are shown in figure 2. As expected from the LB development and continuum error analysis, at θ* = 1, the solver exhibits a much wider usability domain (in terms of non-dimensional viscosity). It is also interesting to note that the maximum achievable Mach number/non-dimensional velocity is reached at 0.65 <θ < 0.85. This can in part be justified by the fact that for θ* < 1, temperature deviation affects the behaviour of the numerical solver in the form of hyperviscosity, while ${\theta }^{\ast }>1$ results in reduced dissipation at higher wavenumbers. Furthermore, lowering the non-dimensional temperature at constant non-dimensional velocity results in increased Mach numbers.

Figure 2. Stability domain (in the ${\theta }^{\ast }-\nu {\delta }_t/{\delta }_x^2$ space) of the SRT collision operator. (a) Maximum achievable Mach number for (from left to right) orders two, three and four (the colourbar shows the Mach number). (b) corresponding maximum non-dimensional velocities. Dashed lines are iso-contours with δu = 0.05 and $\delta \text{Ma}=0.1$ . (Online version in colour.)

• Download figure
• Open in new tab
• Download PowerPoint

While the addition of third-order components does not have any major impact on the stability domain (especially for small non-dimensional viscosities), the fourth-order EDF extends it noticeably, making the solver usable for non-dimensional viscosities as small as $\nu {\delta }_t/{\delta }_x^2={10}^{-3}$ at θ = 0.1. This observation is in agreement with results reported on the effect of the fourth order-term on the stability domain at θ* = 1 [50,62]. It is also worth noting that the use of the EDF of equation (2.20) resulted in the same stability domain as the fourth-order EDF.

Apart from the equilibrium state, the equilibration path is a feature of the collision operator that can have non-negligible effects on the stability domain. As such, the next section will focus on the effects of the different collision models on the linear stability domain.

(c) Advanced collision models: regularized and MRT

The use of the regularized family of collision models is also observed to further extend the stability domain. The obtained results are illustrated in figure 3. Although the fourth-order recursive regularized scheme, as expected, allows the highest non-dimensional velocities to be achieved for θ* ∈ [0.5, 2.3], it is outperformed by the classical second-order regularized collision operator for θ* < 0.5. Although showing good linear stability properties at low non-dimensional temperatures, it must be noted that the second-order expansion used in the classical regularized collision model results in over-dissipation of all physical modes (both shear and acoustic) even in the continuum limit, for θ* < 1. This over-dissipation is also the reason behind the wider linear stability domain of this collision operator at low non-dimensional temperatures. To better illustrate this assertion the spectral dissipation of physical modes using this collision model at three different temperatures are shown in figure 4. It can clearly be observed that lower temperatures results in pronounced over-dissipation of all physical modes, even in the continuum limit.

Figure 3. Stability domain (in the ${\theta }^{\ast }-\nu {\delta }_t/{\delta }_x^2$ space) of the regularized collision operator. (a) Maximum achievable Mach number for (from left to right) orders two, three and four; corresponding maximum non-dimensional velocities are shown in the (b). Plain lines are iso-contours with δu = 0.05 and $\delta \text{Ma}=0.1$ . (Online version in colour.)

• Download figure
• Open in new tab
• Download PowerPoint

Figure 4. Spectral dissipation of the second-order regularized collision operator for three different temperatures: (black symbols) θ* = 0.2, (red symbols) θ* = 0.5 and (blue symbols) θ* = 1. The velocity and non-dimensional viscosities are set to u = 0.02 and $\nu {\delta }_t/{\delta }_x^2={10}^{-4}$ . Only physical modes are displayed. (Online version in colour.)

• Download figure
• Open in new tab
• Download PowerPoint

As mentioned earlier, and detailed in [65], the isothermal fourth-order recursive regularized collision operator is equivalent to an MRT collision model in central Hermite moments space with all non-physical relaxation coefficients set to 1. However, on the same stencil, taking into account thermal effects in the regularization, i.e. temperature-dependent terms in equation (2.39), one would have to use a temperature-scaled central Hermite moments space to recover the same dynamics. Furthermore, the central moments formulation has the additional advantage of allowing for freely tunable parameters for an even wider stability domain and/or enhanced numerical properties (although MRT formulations of the regularized collision operators can also allow for such additional degrees of freedom [74]). As such, to corroborate the previous statements the linear stability domains of both operators were computed. The stability domains for three different choices of ghost relaxation coefficients are shown in figure 5. Looking at these stability domains three observations can be made.

Figure 5. Stability domain (in the ${\theta }^{\ast }-\nu {\delta }_t/{\delta }_x^2$ space) of the MRT collision operator. (a) Central Hermite moments space with the ghost relaxation coefficients set to (from left to right) 0.6, 1, 1.4; (b) temperature-scaled central Hermite moments space with the ghost relaxation coefficients set to (from left to right) 0.6, 1, 1.4. Plain lines are iso-contours with δu = 0.05 and $\delta \text{Ma}=0.1$ . (Online version in colour.)

• Download figure
• Open in new tab
• Download PowerPoint

Regardless of the value of the relaxation coefficient of ghost moments, the temperature-scaled operator systematically outperformed the central moments MRT collision model, especially away from the stencil reference temperature. As such, the temperature-scaled central moments MRT collision operator is a more suitable choice for thermal flows. This observation is in agreement with simulation results reported in [70].

Setting the relaxation coefficients of ghost moments to one, the temperature-scaled central Hermite moments MRT recovers exactly the same linear stability domain as the SRT formulation of the thermal fourth-order recursive regularized operator (on the D2Q9 stencil). For a three-dimensional stencil (such as the D3Q27) one would have to reconstruct the non-equilibrium contributions up to order six to recover the corresponding temperature-scaled central Hermite moments operator. The equivalent MRT model, i.e. with all ghost relaxation coefficients set to one, drastically reduces the number of operations.

Changing the relaxation rate of ghost moments affects the linear stability domain of the operators (both central and temperature-scaled central moments). As such appropriate closures for these free parameters can enhance stability and numerical properties.

Now that the effect of the EDF and collision model on the stability domain have been analysed, the next subsection will focus on validating previously derived correction terms through spectral analyses.

(d) Effect of correction term

To see whether the correction terms found in the literature and derived in previous sections through a CE analysis correct the shortcomings of the standard stencil, spectral dispersion–dissipation curves of the discrete solver were analysed both with and without these terms. The correction terms are evaluated using second-order central difference approximations, unless stated otherwise. The obtained results for the second-order EDF are displayed in figures 6 and 7. In both figures, the flow velocity vector is in the x-direction and only results on the k _x axis are shown. To better showcase the effect of this correction term, figure 6 considers a configuration where θ* = 1 and u _x  = 0.3. It can be seen that the correction terms remove the velocity-dependence of the effective viscosity of shear and acoustic modes in the limit of vanishing wavenumbers (continuum limit). The next plot (figure 7) illustrates this same effect for deviations from the reference temperature. While, as shown there, one observes over-dissipation of the acoustic modes at θ* = 0.5, the introduction of the correction terms restores the correct dissipation rate (independent of temperature) at the continuum limit.

Figure 6. Effect of the correction terms on the spectral behaviour of the shear and acoustic modes with a second-order EDF: using (a) ${\mathrm{\Psi }}_{\alpha }$ and (b) ${\mathrm{\Psi }}_{\alpha }^{\mathrm{\prime }}$ . Blue symbols represent the spectral dissipations without the correction terms while corrected ones are shown in red. u _x  = 0.3, θ* = 1 and $\nu {\delta }_t/{\delta }_x^2=1\times {10}^{-4}$ . Plain black lines are the reference spectral dissipation curves. (Online version in colour.)

• Download figure
• Open in new tab
• Download PowerPoint

Figure 7. Effect of the correction terms on the spectral behaviour of the shear and acoustic modes with a second-order EDF: using (a) ${\mathrm{\Psi }}_{\alpha }$ and (b) ${\mathrm{\Psi }}_{\alpha }^{\mathrm{\prime }}$ . Blue symbols represent the spectral dissipations without the correction terms while corrected ones are shown in red. u _x  = 0.05, θ* = 0.5 and $\nu {\delta }_t/{\delta }_x^2=1\times {10}^{-4}$ . Plain black lines are the reference spectral dissipation curves. (Online version in colour.)

• Download figure
• Open in new tab
• Download PowerPoint

Given the effect of the correction term on the dissipation of acoustic modes at low wavenumbers, it might also be interesting to see whether higher-order evaluations of the correction terms would lead to enhanced results. This effect is illustrated in figures 8 and 9, where the correction terms are computed using second-, fourth- and sixth-order central differences for two different fluid temperatures, i.e. θ* = 0.5 and 1.5. The coefficients of the corresponding finite difference approximation for ${\mathrm{\Psi }}_{\alpha }$ are given in table 1 for the sake of clarity.

Figure 8. Spectral dissipation of the acoustic modes using a fourth-order EDF and the correction terms, i.e. (a) ${\mathrm{\Psi }}_{\alpha }$ and (b) ${\mathrm{\Psi }}_{\alpha }^{\mathrm{\prime }}$ , evaluated using (in black) second-, (in red) fourth-, and (in blue) sixth-order central differences for θ* = 0.5. Plain black lines are the reference spectral dissipation curves. (Online version in colour.)

• Download figure
• Open in new tab
• Download PowerPoint

Figure 9. Spectral dissipation of the acoustic modes using a fourth-order EDF and the correction terms, i.e. (a) ${\mathrm{\Psi }}_{\alpha }$ and (b) ${\mathrm{\Psi }}_{\alpha }^{\mathrm{\prime }}$ , evaluated using (in black) second-, (in red) fourth-, and (in blue) sixth-order central differences for θ* = 1.5. Plain black lines are the reference spectral dissipation curves. (Online version in colour.)

• Download figure
• Open in new tab
• Download PowerPoint

Table 1. Coefficients for finite difference approximations of first-order derivative.

Collapse

order	x − 3δ x	x − 2δ x	x −δ x	x +δ x	x + 2δ x	x + 3δ x
2	0	0	−1/2	1/2	0	0
4	0	1/12	−2/3	2/3	−1/12	0
6	−1/60	3/20	−3/4	3/4	−3/20	1/60

The coefficients for ${\mathrm{\Psi }}_{\alpha }^{\mathrm{\prime }}$ corresponding to those of a second-order derivative are also given in table 2. As observed in figures 8 and 9, while the use of higher-order approximations reduces the deviations from the physical dissipation rate at higher wavenumbers in the first configuration to some extent, the effect is less pronounced in the second case. Furthermore, the first correction, i.e. using ${\mathrm{\Psi }}_{\alpha }$ is more effective in reducing the errors at non-vanishing wavenumbers. It is also worth noting that the second approach to correcting the third-order moments, i.e. using ${\mathrm{\Psi }}_{\alpha }^{\mathrm{\prime }}$ is more complicated than the first approach in practice, as it involves evaluating a second-order derivative with a variable coefficient, i.e. μ/p. The spectral analyses performed in this subsection, confirmed that both correction terms were able to restore (in the continuum limit) the Galilean invariance of the dissipation rate of acoustic and shear modes caused by errors in the third-order moments tensor tied to deviations of both temperature and velocity from the reference state of the stencil. Finally, the correction terms, using central difference approximations, were observed to have little effect on the linear stability domains. However, it is interesting to note that using a first-order upwind approximation to evaluate the correction term ${\mathrm{\Psi }}_{\alpha }$ allows for supersonic simulations for a range of θ* and non-dimensional viscosities. These effects are illustrated in figure 10. To further establish these observations, a number of numerical simulations were performed. Corresponding results are detailed in the next section.

Figure 10. Stability domain (in the ${\theta }^{\ast }-\nu {\delta }_t/{\delta }_x^2$ space) of the MRT collision operator in temperature-scaled central Hermite moments space with the correction term ${\mathrm{\Psi }}_{\alpha }$ evaluated using (a) second-order central difference and (b) first-order upwind approximations. The ghost relaxation coefficients are set to 1. The plain red contours on the r.h.s. figure show the supersonic regions. Colourbars show the Mach number. The spacing between the black dashed contours is δMa = 0.1. (Online version in colour.)

• Download figure
• Open in new tab
• Download PowerPoint

Table 2. Coefficients for finite difference approximations of second-order derivative.

Collapse

order	x − 3δ x	x − 2δ x	x −δ x	x	x +δ x	x + 2δ x	x + 3δ x
2	0	0	1	−2	1	0	0
4	0	−1/12	4/3	−5/2	4/3	−1/12	0
6	1/90	−3/10	3/2	−49/18	3/2	−3/10	1/90

4. Numerical application

(a) Physical properties: sound speed and dissipation rates of shear and acoustic waves

In order to further corroborate observations made in previous sections, three simple pseudo one-dimensional cases are studied here: (a) a freely travelling pressure front, (b) decaying shear waves, and (c) decaying acoustic waves. While the first cases will be used to show that all thermal models studied here can recover the correct temperature-dependent speed of sound, the last two aim at validating the correction terms derived earlier in the manuscript. It is also worth nothing that for all simulations involving energy transfer, a finite-difference solver for the corresponding balance equation is coupled to the LB flow solver [21].

To validate the temperature-dependence of the speed of sound recovered by the LB models, a pseudo one-day domain is separated into two regions with the same initial temperature and velocity and a non-dimensional pressure difference of Δp = 10⁻⁴ (between them). The system is left to evolve and the shock front is tracked over time to compute the speed of sound. The simulations are performed for three different specific heat capacity ratios, i.e. γ = 1.4, 2 and 3, and non-dimensional temperatures θ*. Given that the third-order tensor is not expected to affect the dispersion of acoustic modes (as confirmed by spectral analyses), only the fourth-order Hermite-expanded EDF and those of equation (2.20) are considered here. Both rely on an SRT collision operator. The obtained results are displayed in figure 11 and compared to the theoretical isentropic sound speed $c_s(T)=\sqrt{\gamma k_{B}T/m}$ . As observed in this figure both EDFs correctly capture the sound speed. Furthermore, they exhibit exactly the same behaviour. This was expected as these EDFs are equivalent to each other.

Figure 11. Speed of sound as obtained through the travelling pressure front using the moment matching and Hermite polynomial (fourth-order) EDFs at different temperatures. Circles represent results using the Hermite polynomials EDF, crosses those using the moment matching EDF and plain lines theoretical isentropic sound speeds. (Online version in colour.)

• Download figure
• Open in new tab
• Download PowerPoint

In the second case, a decaying shear wave is modelled. By monitoring the time evolution of the maximum velocity of this wave and fitting an exponential function to it, the effective numerical kinetic viscosity can be measured. A domain of size 2 ×N _y (with N _y  = 200) is used to perform the simulations. The initial conditions are defined as [14]

$$\begin{array}{rl}& \rho (x,y)=1,\end{array}$$ 4.1a

$$\begin{array}{rl}& u_x(x,y)=1\times {10}^{-4}\sin \left(\frac{2\pi y}{N_y}\right),\end{array}$$ 4.1b

$$\begin{array}{rl}& u_y(x,y)=0,\end{array}$$ 4.1c

$$\begin{array}{rl}\text{and}& {\theta }^{\ast }={\theta }_0^{\ast },\end{array}$$ 4.1d

while the non-dimensional kinematic viscosity is set to $\nu {\delta }_t/{\delta }_x^2=0.1$ . The simulations are performed over 50 000 time steps using second- (both with and without the corrections for off-diagonal components) and third-order EDFs for θ* ∈ [0.1 − 1.2]. The obtained results are displayed in figure 12. As observed in figure 12, the second-order EDF is unable to recover the correct dissipation rate for the shear mode. At lower temperatures, the deviations from the reference temperature result in pronounced over-dissipation of this mode, explaining the unexpected stability of the second-order regularized collision operator at very low temperatures. It can also be observed that the addition of the off-diagonal components corrections (both approaches) restores the appropriate dissipation rate for the second-order EDF.

Figure 12. Effective dissipation rate of the shear mode for the third-order EDF and the second-order EDF with (both using ${\mathrm{\Psi }}_{\alpha }$ and ${\mathrm{\Psi }}_{\alpha }^{\mathrm{\prime }}$ ) and without the correction terms obtained from the decaying shear wave test-case. The reference viscosity is shown with a plain line. (Online version in colour.)

• Download figure
• Open in new tab
• Download PowerPoint

To assess the accuracy of the acoustic modes dissipation, the test-case presented in [75] is used. In the limit of small velocity and density variations, the linear regime holds and acoustic waves can be modelled as (only considering propagation in the x-direction) [75]:

$${\mathrm{\partial }}_t^{2}u=c_s^2{\mathrm{\partial }}_x^{2}u+(\frac{4}{3}\nu +\eta ){\mathrm{\partial }}_t{\mathrm{\partial }}_x^{2}u,$$ 4.2

where η is the bulk viscosity. To evaluate the effective viscosity of the acoustic modes, the density is initialized as

$$\rho (x)={\rho }_0+\delta \rho \sin \left(\frac{2\pi y}{N_y}\right),$$ 4.3

and the initial velocity is set to zero. The time evolution of the energy in the domain, evaluated as E =u ² +θ*(ρ −ρ ₀), is monitored over 50 000 steps. Given that in the context of the linear regime, wave decay is an exponential function of the effective viscosity, the effective viscosity is then extracted via a simple fitting process. In the context of this study, N _y  = 200, ρ ₀ = 1 and δρ = 10⁻⁶ are used. As with the previous case, different temperatures are considered for the SRT collision operator with a third-order equilibrium both with (using both approaches) and without the correction term. The results are shown in figure 13. The introduction of the correction, either one, effectively removes the temperature-dependence of the acoustic modes dissipation rate in the limit of continuity (fully resolved features).

Figure 13. Effective dissipation rate of the acoustic modes for the third-order EDF with (both ${\mathrm{\Psi }}_{\alpha }$ and ${\mathrm{\Psi }}_{\alpha }^{\mathrm{\prime }}$ ) and without the correction term obtained from the decaying acoustic wave test-case. The reference viscosity is shown with a plain line. (Online version in colour.)

• Download figure
• Open in new tab
• Download PowerPoint

(b) Supersonic flow around the NACA0012 aerofoil

To illustrate the possibilities within the context of a compressible LBM formulation on standard stencils and its limitations, the temperature-scaled central moments MRT collision model along with the correction term ${\mathrm{\Psi }}_{\alpha }$ have been coupled to a finite-difference solver for the energy balance equation [21] to model the flow past the NACA0012 aerofoil in the supersonic regime. The configuration modelled in this subsection follows that used in [76,77]. The results shown here being merely a proof of concept, at the difference of previous studies where the Reynolds number was set to 10⁴, here the simulation is performed at Re = 2000. Given the restrictive stability domain at θ* ≠ 1, a simulation at Re = 10⁴ would have resulted in a large domain.

The simulation domain consists of a rectangular box of size 12 l _c  × 6 l _c , where l _c is the aerofoil chord length. It is bound by constant velocity and second-order zero-gradient boundary conditions on the left- and right-hand sides, respectively [78]. To emulate far-field free flow conditions, the inlet velocity is also enforced at the top and bottom boundaries. The solid boundaries on the aerofoil surface are taken into account through the half-way bounce-back scheme [79]. The Mach number at the inlet is set to 1.5. To achieve stability at such a Mach number, the temperature and velocity are set to θ* = 0.13 and u ₀ = 0.3125, respectively. Furthermore, based on results from the stability analysis, the correction term ${\mathrm{\Psi }}_{\alpha }$ was evaluated using a first-order upwind approximation. The chord length was resolved by 160 grid-points. To accelerate convergence, the fluid velocity in the domain was initialized with the inlet velocity. The simulation, e.g. shock front position, was observed to converge after approximately 15 000 steps. The obtained Mach number and density fields are displayed in figure 14. The flow structure, as observed here, is able to capture the main features of the flow, namely a shock upstream the leading edge and a weaker front in the vicinity of the trailing edge. To further validate the results obtained in the present work, the distribution of the pressure coefficient, C _p around the aerofoil is compared to results reported in [76]. The results along with pressure coefficient field are shown in figure 15.

Figure 14. (a) Mach number and (b) density fields for the NACA0012 supersonic case. (Online version in colour.)

• Download figure
• Open in new tab
• Download PowerPoint

Figure 15. (a) Pressure coefficient field and (b) distribution around the aerofoil. The sampling line used for the plot in (b) is shown with black/white dashed lines in (a). Red plain lines in (b) are reference data from [76] while black symbols are from the present study. (Online version in colour.)

• Download figure
• Open in new tab
• Download PowerPoint

The pressure coefficient is computed as $C_p=2(p-p_0)/{\rho }_0{u}_0^2$ , with $p=\rho c_s^2{\theta }^{\ast }$ and where ρ ₀ and p ₀ are the inlet density and pressure. As observed in figure 15, although relying on a relatively coarse grid (as compared to [76] where the chord is resolved with 800 points), the maximum pressure coefficient at the shock front, its position and the shock wave thickness are in good agreement with reference solutions.

While able to model this supersonic flow, it is interesting to note that the solver was unstable for higher Reynolds and/or Mach numbers. Furthermore, other collision operators such as SRT or MRT in central moments space were unstable for this configuration. While one can further extend the stability domain via a full use of the free parameters space or adaptive dissipation mechanisms such as the entropic family of collision models, it can also be used with the concept of shifted stencils to reach higher Reynolds and Mach numbers at a relatively low cost. As a proof of concept, the next subsection will focus on modelling the isothermal vortex convection case at different Mach numbers using shifted stencils.

(c) Isothermal convected vortex on shifted stencils

The in-viscid convected vortex is a case that shows the ability of the scheme to deal with vortical structures. The velocity field is initialized as

$$\begin{array}{rl}& u_x=u_0-\delta u_0\frac{y-y_c}{r_0}\exp -\frac{{\mathit{r}}^2}{2r_0^2},\end{array}$$ 4.4a

$$\begin{array}{rl}& u_y=\delta u_0\frac{x-x_c}{r_0}\exp -\frac{{\mathit{r}}^2}{2r_0^2}\end{array}$$ 4.4b

$$\begin{array}{rl}\text{and}& \rho =1-\frac{{\delta u_0}^2}{2{\theta }^{\ast }}\exp -\frac{{\mathit{r}}^2}{r_0^2},\end{array}$$ 4.4c

where u ₀ is the convection velocity, δu ₀ is the vortex strength and r ₀ is the vortex radius. In a previous publication, shifted stencils were used to attain very large convective Mach numbers [28]. However, it was shown that for inviscid flows, simulations were limited to Mach numbers corresponding to integer multiples of the stencil Mach number, i.e. $Ma=k\sqrt{3}$ with $k\in \mathbb{Z}$ . For θ* = 1, the only way to go beyond this restriction is through off-lattice propagation supplemented with higher-order reconstruction on lattices. For θ* ≠ 1, the effective Mach number is $\text{Ma}=k\sqrt{3/{\theta }^{\ast }}$ . As such deviations from the reference temperature can be used to overcome the previously cited issue on shifted stencils. For a shift of U _x  =δx/δt one has Ma = 1.73. To showcase the added value of shifted stencils accompanied with a flexible reference temperature simulations over a wide range of Mach numbers are performed using shifts of U _x  =δx/δt and U _x  = 2δx/δt. The results are shown in figure 16. The simulations were performed on a domain of size 128 × 128 with δu ₀ = 10⁻³. The results presented here show that changing the non-dimensional temperature θ* and using shifted stencils to match the background flow velocity one can perform inviscid simulations of the isothermal convected vortex at Mach numbers others than those imposed by the shift in the stencil. The concept of shifted stencils can therefore be used to perform supersonic simulations with higher Reynolds number at low computational costs.

Figure 16. Comparison of (red dashed) initial and (black plain) final pressure iso-contours for different Mach numbers: (a) 1.73 and θ* = 1, (b) 2 and θ* = 0.75, (c) 2.3 and θ* = 0.56, (d) 2.6 and θ* = 0.44, (e) 3.1 and θ* = 1.27, (f ) 3.5 and θ* = 0.98, (g) 3.8 and θ* = 0.83 and (h) 4.1 and θ* = 0.71. In the upper row a shift of $U_x=\frac{\delta x}{\delta t}$ is used while the lower row uses U _x  = 2(δx/δt). (Online version in colour.)

• Download figure
• Open in new tab
• Download PowerPoint

5. Conclusion

The aim of the present work was to systematically study the effect of local temperature on the behaviour of the lattice Boltzmann solver and its numerical properties for polynomials EDFs. With the help of the von Neumann formalism, the linear stability domain for different orders of the EDF and more advanced collision models such as the regularized and MRT collision operators in different moments bases were established. Furthermore, this spectral analysis established the effect of the different correction terms on shear and acoustic modes dissipation and the stability domain.

A comparative study of the properties of Hermite expansion-based EDFs of different orders and those obtained through moment matching showed that the latter is equivalent to the former with a fourth-order expansion. The latter is in practice more interesting as it is a (computationally) more efficient way of writing the EDF. Among the different collision operators considered in this study the MRT model in temperature-scaled Hermite polynomial space showed the widest stability domain. It was shown that setting the free parameters in the model to one, the SRT thermal fourth-order recursive regularized operator would be recovered. The addition of the different correction terms for the third-order moments tensor was also shown to restore the Galilean invariance of the dissipation rate of physical modes at the NS level. It was also observed that using a first-order upwind approximation (with the first approach, ${\mathrm{\Psi }}_{\alpha }$ ) could allow for supersonic flow simulations. This observation was further confirmed through a simulation of a supersonic flow around the NACA0012 aerofoil.

While with the addition of the correction term and careful consideration of stability domains one can perform thermal simulations on standard stencils, as readily shown in the literature, all results presented here show that the reference state is the global attractor of the classical LB scheme. Therefore, as shown through the isothermal vortex convection case, adapting the global attractor of the model, i.e. reference speed and temperature, to the local fluid state is a promising approach for high Mach/Reynolds number applications [29–32].

Data accessibility

This article has no additional data.

Authors' contributions

S.A.H. carried out all mathematical developments and simulations and drafted the manuscript. N.D. and D.T. coordinated the project and contributed to the writing of the manuscript. All authors read and approved the manuscript.

Competing interests

The authors declare that they have no conflict of interest.

Funding

The financial support of the German research foundation (DFG), Germany within the graduate college for 'Micro-Macro-Interactions in Structured Media and Particle Systems' (GRK 1554) is gratefully acknowledged.

Acknowledgements

S.A.H. would like to acknowledge fruitful discussions on the regularized collision models and von Neumann analysis with C. Coreixas and G. Wissocq.

Appendix A. Details on Hermite polynomials and coefficients

For the sake of clarity, the Hermite polynomials and coefficients after discretization in space and time are given here using the notation adopted by the LBM community:

$$\begin{array}{rl}& {\mathcal{H}}_{\alpha ,0}=1,\end{array}$$ A 1a

$$\begin{array}{rl}& {\mathcal{H}}_{\alpha ,i_1}=c_{\alpha ,i_1},\end{array}$$ A 1b

$$\begin{array}{rl}& {\mathcal{H}}_{i_1{i}_2}=c_{\alpha ,i_1}{c}_{\alpha ,i_2}-c_s^2{\delta }_{i_1{i}_2},\end{array}$$ A 1c

$$\begin{array}{rl}& {\mathcal{H}}_{\alpha ,i_1{i}_2{i}_3}=c_{\alpha ,i_1}{c}_{\alpha ,i_2}{c}_{\alpha ,i_3}-c_s^2{[c_{\alpha ,i_1}{\delta }_{i_2{i}_3}]}_{\text{cyc}},\end{array}$$ A 1d

$$\begin{array}{rl}\text{and}& {\mathcal{H}}_{\alpha ,i_1{i}_2{i}_3{i}_4}=c_{\alpha ,i_1}{c}_{\alpha ,i_2}{c}_{\alpha ,i_3}{c}_{\alpha ,i_4}+c_s^4{[{\delta }_{i_1{i}_2}{\delta }_{i_3{i}_4}]}_{\text{cyc}}-c_s^2{[c_{\alpha ,i_3}{c}_{\alpha ,i_4}{\delta }_{i_1{i}_2}]}_{\text{cyc}},\end{array}$$ A 1e

and

$$\begin{array}{rl}& a_0^{(\text{eq})}=\rho ,\end{array}$$ A 2a

$$\begin{array}{rl}& a_{i_1}^{(\text{eq})}=\rho u_{i_1},\end{array}$$ A 2b

$$\begin{array}{rl}& a_{i_1{i}_2}^{(\text{eq})}=\rho u_{i_1}{u}_{i_2}+\rho c_s^2({\theta }^{\ast }-1),\end{array}$$ A 2c

$$\begin{array}{rl}& a_{i_1{i}_2{i}_3}^{(\text{eq})}=\rho u_{i_1}{u}_{i_2}{u}_{i_3}+\rho c_s^2({\theta }^{\ast }-1){[u_{i_1}{\delta }_{i_2{i}_3}]}_{\text{cyc}},\end{array}$$ A 2d

$$\begin{array}{rl}\text{and}& a_{i_1{i}_2{i}_3{i}_4}^{(\text{eq})}=\rho u_{i_1}{u}_{i_2}{u}_{i_3}{u}_{i_4}+\rho c_s^4{({\theta }^{\ast }-1)}^2{[{\delta }_{i_1{i}_2}{\delta }_{i_3{i}_4}]}_{\text{cyc}}+\rho c_s^2({\theta }^{\ast }-1){[u_{i_1}{u}_{i_2}{\delta }_{i_3{i}_4}]}_{\text{cyc}},\end{array}$$ A 2e

where c _s is the sound speed at the reference state non-dimensionalized with δ _x /δ _t , i.e. $c_s^2=1/3$ . Furthermore, c _α,i  ∈ {0, ± 1}. The EDFs are then expanded as

$$f_{\alpha }^{(\text{eq},N)}=w_{\alpha }\sum _{n=0}^N\frac{1}{n!c_s^{2n}}{\mathit{a}}_n^{(\text{eq})}:{\mathcal{H}}_{\alpha ,n}.$$ A 3

Appendix B. Jacobian of the moment-matching discrete EDF

Given that detailed derivation and corresponding expressions for the Jacobians of the collision operators based on the Hermite-expanded discrete EDF have already been given in [28], only the Jacobian of the discrete EDF resulting from the moment-matching procedure will be given here. Let us first rewrite the EDF of equation (2.20) as

$$f_{\alpha }^{(\text{eq})}=\rho \prod _{i=x,y}{\mathrm{\Lambda }}_{\alpha ,i}^{1-\mid c_{\alpha ,i}\mid }{\mathrm{\Gamma }}_{\alpha ,i}^{\mid c_{\alpha ,i}\mid },$$ B 1

where:

$${\mathrm{\Lambda }}_{\alpha ,i}=(1-u_i^2-\theta ),$$ B 2

and

$${\mathrm{\Gamma }}_{\alpha ,i}=\left(\frac{u_i^2+\theta +c_{\alpha ,i}{u}_i}{2}\right).$$ B 3

With this notation, the Jacobian can be developed as

$${\mathrm{\partial }}_{f_{\beta }}{f}_{\alpha }^{(\text{eq})}=\prod _{i=x,y}{\mathrm{\Lambda }}_{\alpha ,i}^{1-\mid c_{\alpha ,i}\mid }{\mathrm{\Gamma }}_{\alpha ,i}^{\mid c_{\alpha ,i}\mid }(1+\rho \sum _{j=x,y}\frac{1-\mid c_{\alpha ,i}\mid }{{\mathrm{\Lambda }}_{\alpha ,j}^{1-\mid c_{\alpha ,j}\mid }}{\mathrm{\partial }}_{f_{\beta }}{\mathrm{\Lambda }}_{\alpha ,j}+\frac{\mid c_{\alpha ,i}\mid }{{\mathrm{\Gamma }}_{\alpha ,j}^{\mid c_{\alpha ,j}\mid }}{\mathrm{\partial }}_{f_{\beta }}{\mathrm{\Gamma }}_{\alpha ,j}),$$ B 4

where the Jacobians of Λ _α,i and Γ _α,i can be readily computed as

$${\mathrm{\partial }}_{f_{\beta }}{\mathrm{\Lambda }}_{\alpha ,i}=\frac{2}{\rho }(c_{\beta ,i}{u}_i-u_i^2),$$ B 5

and

$${\mathrm{\partial }}_{f_{\beta }}{\mathrm{\Gamma }}_{\alpha ,i}=\frac{1}{\rho }(c_{\beta ,i}{u}_i-u_i^2)+\frac{1}{2\rho }(c_{\alpha ,i}{c}_{\beta ,i}-c_{\alpha ,i}{u}_i),$$ B 6

Footnotes

One contribution of 15 to a theme issue 'Fluid dynamics, soft matter and complex systems: recent results and new methods'.

© 2020 The Author(s)

Published by the Royal Society. All rights reserved.

References

Reference

• 1.
  Succi S

  . 2015 Lattice Boltzmann 2038. Europhys. Lett. 109 , 50001. (doi:10.1209/0295-5075/109/50001) Crossref, Google Scholar

• 2.
  Succi S

  . 2001 The lattice Boltzmann equation: for fluid dynamics and beyond . Oxford, UK: Oxford University Press. Google Scholar

• 3.
  Montessori A, Falcucci G

  . 2018 Lattice Boltzmann modeling of complex flows for engineering applications . San Rafael, CA: Morgan & Claypool Publishers. Google Scholar

• 4.
  Falcucci G, Ubertini S, Bella G, Succi S

  . 2013 Lattice Boltzmann simulation of cavitating flows. Commun. Comput. Phys. 13 , 685–695. (doi:10.4208/cicp.291011.270112s) Crossref, ISI, Google Scholar

• 5.
  Falcucci G, Chibbaro S, Succi S, Shan X, Chen H

  . 2008 Lattice Boltzmann spray-like fluids. Europhys. Lett. 82 , 24005. (doi:10.1209/0295-5075/82/24005) Crossref, Google Scholar

• 6.
  Chopard B, Ouared R, Ruefenacht DA, Yilmaz H

  . 2007 Lattice Boltzmann modeling of thrombosis in giant aneurysms. Int. J. Mod. Phys. C 18 , 712–721. (doi:10.1142/S0129183107010978) Crossref, ISI, Google Scholar

• 7.
  Dorschner B, Bösch F, Chikatamarla SS, Boulouchos K, Karlin IV

  . 2016 Entropic multi-relaxation time lattice Boltzmann model for complex flows. J. Fluid Mech. 801 , 623–651. (doi:10.1017/jfm.2016.448) Crossref, ISI, Google Scholar

• 8.
  Chikatamarla SS, Karlin IV

  . 2015 Entropic lattice Boltzmann method for multiphase flows. Phys. Rev. Lett. 114 , 174502. (doi:10.1103/PhysRevLett.114.174502) Crossref, PubMed, ISI, Google Scholar

• 9.
  Wissocq G, Sagaut P, Boussuge JF

  . 2019 An extended spectral analysis of the lattice Boltzmann method: modal interactions and stability issues. J. Comput. Phys. 380 , 311–333. (doi:10.1016/j.jcp.2018.12.015) Crossref, ISI, Google Scholar

• 10.
  Frapolli N

  . 2017 Entropic lattice Boltzmann models for thermal and compressible flows. PhD thesis, ETH Zurich. Google Scholar

• 11.
  Frapolli N, Chikatamarla SS, Karlin IV

  . 2014 Multispeed entropic lattice Boltzmann model for thermal flows. Phys. Rev. E 90 , 043306. (doi:10.1103/PhysRevE.90.043306) Crossref, ISI, Google Scholar

• 12.
  Coreixas C, Wissocq G, Puigt G, Boussuge JF, Sagaut P

  . 2017 Recursive regularization step for high-order lattice Boltzmann methods. Phys. Rev. E 96 , 033306. (doi:10.1103/PhysRevE.96.033306) Crossref, PubMed, ISI, Google Scholar

• 13.
  Li Q, He YL, Wang Y, Tao WQ

  . 2007 Coupled double-distribution-function lattice Boltzmann method for the compressible Navier–Stokes equations. Phys. Rev. E 76 , 056705. (doi:10.1103/PhysRevE.76.056705) Crossref, ISI, Google Scholar

• 14.
  Saadat MH, Bösch F, Karlin IV

  . 2019 Lattice Boltzmann model for compressible flows on standard lattices: variable Prandtl number and adiabatic exponent. Phys. Rev. E 99 , 013306. (doi:10.1103/PhysRevE.99.013306) Crossref, PubMed, ISI, Google Scholar

• 15.
  Saadat MH, Bösch F, Karlin IV

  . 2020 Semi-Lagrangian lattice Boltzmann model for compressible flows on unstructured meshes. Phys. Rev. E E101, 023311. Crossref, ISI, Google Scholar

• 16.
  Wilde D, Krämer A, Reith D, Foysi H

  . 2019Semi-Lagrangian lattice Boltzmann method for compressible flows. (http://arxiv.org/abs/1910.13918) Google Scholar

• 17.
  Prasianakis NI, Karlin IV

  . 2007 Lattice Boltzmann method for thermal flow simulation on standard lattices. Phys. Rev. E 76 , 016702. (doi:10.1103/PhysRevE.76.016702) Crossref, ISI, Google Scholar

• 18.
  Prasianakis NI, Karlin IV

  . 2008 Lattice Boltzmann method for simulation of compressible flows on standard lattices. Phys. Rev. E 78 , 016704. (doi:10.1103/PhysRevE.78.016704) Crossref, ISI, Google Scholar

• 19.
  Feng YL, Sagaut P, Tao WQ

  . 2015 A three dimensional lattice model for thermal compressible flow on standard lattices. J. Comput. Phys. 303 , 514–529. (doi:10.1016/j.jcp.2015.09.011) Crossref, ISI, Google Scholar

• 20.
  Safari H, Krafczyk M, Geier M

  . In press. A Lattice Boltzmann model for thermal compressible flows at low Mach numbers beyond the Boussinesq approximation. Comput. Fluids . (doi:10.1016/j.compfluid.2018.04.016) Google Scholar

• 21.
  Hosseini SA, Safari H, Darabiha N, Thévenin D, Krafczyk M

  . 2019 Hybrid lattice Boltzmann-finite difference model for low Mach number combustion simulation. Combust. Flame 209 , 394–404. (doi:10.1016/j.combustflame.2019.07.041) Crossref, ISI, Google Scholar

• 22.
  Feng YL, Guo SL, Tao WQ, Sagaut P

  . 2018 Regularized thermal lattice Boltzmann method for natural convection with large temperature differences. Int. J. Heat Mass Transfer 125 , 1379–1391. (doi:10.1016/j.ijheatmasstransfer.2018.05.051) Crossref, ISI, Google Scholar

• 23.
  Feng YL, Boivin P, Jacob J, Sagaut P

  . 2019 Hybrid recursive regularized thermal lattice Boltzmann model for high subsonic compressible flows. J. Comput. Phys. 394 , 82–99. (doi:10.1016/j.jcp.2019.05.031) Crossref, ISI, Google Scholar

• 24.
  Dellar PJ

  . 2014 Lattice Boltzmann algorithms without cubic defects in Galilean invariance on standard lattices. J. Comput. Phys. 259 , 270–283. (doi:10.1016/j.jcp.2013.11.021) Crossref, ISI, Google Scholar

• 25.
  Frapolli NN, Chikatamarla SS, Karlin IV

  . 2015 Entropic lattice Boltzmann model for compressible flows. Phys. Rev. E 92 , 061301. (doi:10.1103/PhysRevE.92.061301) Crossref, ISI, Google Scholar

• 26.
  Latt J, Coreixas C, Beny J, Parmigiani A

  . 2020 Efficient supersonic flow simulations using lattice Boltzmann methods based on numerical equilibria. Phil. Trans. R. Soc. A 378 , 20190559. (doi:10.1098/rsta.2019.0559) Link, Google Scholar

• 27.
  Frapolli N, Chikatamarla SS, Karlin IV

  . 2016 Lattice kinetic theory in a comoving Galilean reference frame. Phys. Rev. Lett. 117 , 010604. (doi:10.1103/PhysRevLett.117.010604) Crossref, PubMed, ISI, Google Scholar

• 28.
  Hosseini SA, Coreixas C, Darabiha N, Thévenin D

  . 2019 Extensive analysis of the lattice Boltzmann method on shifted stencils. Phys. Rev. E 100 , 063301. (doi:10.1103/PhysRevE.100.063301) Crossref, PubMed, ISI, Google Scholar

• 29.
  Sun C

  . 1998 Lattice Boltzmann models for high speed flows. Phys. Rev. E 58 , 7283–7287. (doi:10.1103/PhysRevE.58.7283) Crossref, ISI, Google Scholar

• 30.
  Sun C

  . 2000 Adaptive lattice Boltzmann model for compressible flows: viscous and conductive properties. Phys. Rev. E 61 , 2645–2653. (doi:10.1103/PhysRevE.61.2645) Crossref, ISI, Google Scholar

• 31.
  Sun C, Hsu AT

  . 2003 Three-dimensional lattice Boltzmann model for compressible flows. Phys. Rev. E 68 , 016303. (doi:10.1103/PhysRevE.68.016303) Crossref, ISI, Google Scholar

• 32.
  Dorschner B, Bösch F, Karlin IV

  . 2018 Particles on demand for kinetic theory. Phys. Rev. Lett. 121 , 130602. (doi:10.1103/PhysRevLett.121.130602) Crossref, PubMed, ISI, Google Scholar

• 33.
  Nadiga BT, Pullin DI

  . 1994 A method for near-equilibrium discrete-velocity gas flows. J. Comput. Phys. 112 , 162–172. (doi:10.1006/jcph.1994.1089) Crossref, ISI, Google Scholar

• 34.
  Nadiga BT

  . 1995 An Euler solver based on locally adaptive discrete velocities. J. Stat. Phys. 81 , 129–146. (doi:10.1007/BF02179972) Crossref, ISI, Google Scholar

• 35.
  Huang J, Xu F, Vallières M, Feng DH, Qian YH, Fryxell B, Strayer MR

  . 1997 A thermal LBGK model for large density and temperature differences. Int. J. Mod. Phys. C 8 , 827–841. (doi:10.1142/S0129183197000710) Crossref, ISI, Google Scholar

• 36.
  Guo Z, Shu C

  . 2013 Lattice Boltzmann method and its applications in engineering , vol. 3. Singapore: World Scientific. Crossref, Google Scholar

• 37.
  He X, Luo LS

  . 1997 A priori derivation of the lattice Boltzmann equation. Phys. Rev. E 55 , R6333. (doi:10.1103/PhysRevE.55.R6333) Crossref, ISI, Google Scholar

• 38.
  He X, Luo LS

  . 1997 Theory of the lattice Boltzmann method: from the Boltzmann equation to the lattice Boltzmann equation. Phys. Rev. E 56 , 6811. (doi:10.1103/PhysRevE.56.6811) Crossref, ISI, Google Scholar

• 39.
  Shim JW, Gatignol R

  . 2011 Thermal lattice Boltzmann method based on a theoretically simple derivation of the Taylor expansion. Phys. Rev. E 83 , 046710. (doi:10.1103/PhysRevE.83.046710) Crossref, ISI, Google Scholar

• 40.
  He X, Shan X, Doolen GD

  . 1998 Discrete Boltzmann equation model for nonideal gases. Phys. Rev. E 57 , R13. (doi:10.1103/PhysRevE.57.R13) Crossref, ISI, Google Scholar

• 41.
  Shan X, He X

  . 1998 Discretization of the velocity space in the solution of the Boltzmann equation. Phys. Rev. Lett. 80 , 65. (doi:10.1103/PhysRevLett.80.65) Crossref, ISI, Google Scholar

• 42.
  Shan X, Yuan XF, Chen H

  . 2006 Kinetic theory representation of hydrodynamics: a way beyond the Navier–Stokes equation. J. Fluid Mech. 550 , 413–441. (doi:10.1017/S0022112005008153) Crossref, ISI, Google Scholar

• 43.
  Philippi PC, Hegele LA, Surmas R, Siebert DN, Dos Santos LOE

  . 2007 From the Boltzmann to the lattice Boltzmann equation: beyond BGK collision models. Int. J. Mod. Phys. C 18 , 556–565. (doi:10.1142/S0129183107010796) Crossref, ISI, Google Scholar

• 44.
  Philippi PC, Hegele LA, Dos Santos LOE, Surmas R

  . 2006 From the continuous to the lattice Boltzmann equation: the discretization problem and thermal models. Phys. Rev. E 73 , 056702. (doi:10.1103/PhysRevE.73.056702) Crossref, ISI, Google Scholar

• 45.
  Struchtrup H

  . 2005Macroscopic transport equations for rarefied gas flows. In Macroscopic transport equations for rarefied gas flows, pp. 145–160. Berlin, Germany: Springer. Google Scholar

• 46.
  Piaud B, Blanco S, Fournier R, Ambruş VE, Sofonea V

  . 2014 Gauss quadratures–the keystone of lattice Boltzmann models. Int. J. Mod. Phys. C 25 , 1340016. (doi:10.1142/S0129183113400160) Crossref, ISI, Google Scholar

• 47.
  Grad H

  . 1949 On the kinetic theory of rarefied gases. Commun. Pure Appl. Math. 2 , 331–407. (doi:10.1002/cpa.3160020403) Crossref, ISI, Google Scholar

• 48.
  Grad H

  . 1958Principles of the kinetic theory of gases. In Thermodynamics of gases (ed S FlŁügge), pp. 205–294. Berlin, Germany: Springer. Google Scholar

• 49.
  Grad H

  . 1949 Note on N-dimensional Hermite polynomials. Commun. Pure Appl. Math. 2 , 325–330. (doi:10.1002/cpa.3160020402) Crossref, ISI, Google Scholar

• 50.
  Hosseini SA, Coreixas C, Darabiha N, Thévenin D

  . 2019 Stability of the lattice kinetic scheme and choice of the free relaxation parameter. Phys. Rev. E 99 , 063305. (doi:10.1103/PhysRevE.99.063305) Crossref, PubMed, ISI, Google Scholar

• 51.
  Krämer A, Küllmer K, Reith D, Joppich W, Foysi H

  . 2017 Semi-Lagrangian off-lattice Boltzmann method for weakly compressible flows. Phys. Rev. E 95 , 023305. (doi:10.1103/PhysRevE.95.023305) Crossref, PubMed, ISI, Google Scholar

• 52.
  He X, Chen S, Doolen GD

  . 1998 A novel thermal model for the lattice Boltzmann method in incompressible limit. J. Comput. Phys. 146 , 282–300. (doi:10.1006/jcph.1998.6057) Crossref, ISI, Google Scholar

• 53.
  Dellar PJ

  . 2013 An interpretation and derivation of the lattice Boltzmann method using Strang splitting. Comput. Math. Appl. 65 , 129–141. (doi:10.1016/j.camwa.2011.08.047) Crossref, ISI, Google Scholar

• 54.
  Karlin IV, Asinari P

  . 2010 Factorization symmetry in the lattice Boltzmann method. Physica A 389 , 1530–1548. (doi:10.1016/j.physa.2009.12.032) Crossref, ISI, Google Scholar

• 55.
  Shim JW, Gatignol R

  . 2013 How to obtain higher-order multivariate Hermite expansion of Maxwell–Boltzmann distribution by using Taylor expansion? Z Angew. Math. Phys. 64 , 473–482. (doi:10.1007/s00033-012-0265-1) Crossref, ISI, Google Scholar

• 56.
  Guo ZL, Zheng C, Shi B

  . 2002 Discrete lattice effects on the forcing term in the lattice Boltzmann method. Phys. Rev. E 65 , 046308. (doi:10.1103/PhysRevE.65.046308) Crossref, PubMed, ISI, Google Scholar

• 57.
  Bhatnagar PL, Gross EP, Krook M

  . 1954 A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys. Rev. 94 , 511. (doi:10.1103/PhysRev.94.511) Crossref, ISI, Google Scholar

• 58.
  Latt J, Chopard B

  . 2006 Lattice Boltzmann method with regularized pre-collision distribution functions. Math. Comput. Simul. 72 , 165–168. (doi:10.1016/j.matcom.2006.05.017) Crossref, ISI, Google Scholar

• 59.
  Ladd AJC, Verberg R

  . 2001 Lattice-Boltzmann simulations of particle-fluid suspensions. J. Stat. Phys. 104 , 1191–1251. (doi:10.1023/A:1010414013942) Crossref, ISI, Google Scholar

• 60.
  Chen H, Zhang R, Staroselsky I, Jhon M

  . 2006 Recovery of full rotational invariance in lattice Boltzmann formulations for high Knudsen number flows. Physica A 362 , 125–131. (doi:10.1016/j.physa.2005.09.008) Crossref, ISI, Google Scholar

• 61.
  Malaspinas O

  . 2015Increasing stability and accuracy of the lattice Boltzmann scheme: recursivity and regularization. arXiv:1505.06900 [physics.u-dyn]. Google Scholar

• 62.
  Coreixas C

  . 2018 High-order extension of the recursive regularized lattice Boltzmann method. PhD thesis, Université de Toulouse. Google Scholar

• 63.
  Mattila KK, Philippi PC, Hegele LA

  . 2017 High-order regularization in lattice-Boltzmann equations. Phys. Fluids 29 , 046103. (doi:10.1063/1.4981227) Crossref, ISI, Google Scholar

• 64.
  Hosseini SA, Darabiha N, Thévenin D

  . 2019 Theoretical and numerical analysis of the lattice kinetic scheme for complex flow simulations. Phys. Rev. E 99 , 023305. (doi:10.1103/PhysRevE.99.023305) Crossref, PubMed, ISI, Google Scholar

• 65.
  Coreixas C, Chopard B, Latt J

  . 2019 Comprehensive comparison of collision models in the lattice Boltzmann framework: theoretical investigations. Phys. Rev. E 100 , 033305. (doi:10.1103/PhysRevE.100.033305) Crossref, PubMed, ISI, Google Scholar

• 66.
  d'Humières D

  . 1992Generalized lattice-Boltzmann equations. In Rarefied gas dynamics: theory and simulation, pp. 450–458. Progress in Astronautics and Aeronautics, vol. 159. Washington, DC: AIAA. Google Scholar

• 67.
  d'Humières D

  . 2002 Multiple–relaxation–time lattice Boltzmann models in three dimensions. Phil. Trans. R. Soc. Lond. A 360 , 437–451. (doi:10.1098/rsta.2001.0955) Link, ISI, Google Scholar

• 68.
  Geier M, Greiner A, Korvink JG

  . 2006 Cascaded digital lattice Boltzmann automata for high Reynolds number flow. Phys. Rev. E 73 , 066705. (doi:10.1103/PhysRevE.73.066705) Crossref, ISI, Google Scholar

• 69.
  Dubois F, Fevrier T, Graille B

  . 2015 Lattice Boltzmann schemes with relative velocities. Commun. Comput. Phys. 17 , 1088–1112. (doi:10.4208/cicp.2014.m394) Crossref, ISI, Google Scholar

• 70.
  Li X, Shi Y, Shan X

  . 2019 Temperature-scaled collision process for the high-order lattice Boltzmann model. Phys. Rev. E 100 , 013301. (doi:10.1103/PhysRevE.100.013301) Crossref, PubMed, ISI, Google Scholar

• 71.
  Hirsch C

  . 2007 Numerical computation of internal and external flows: the fundamentals of computational fluid dynamics . Amsterdam, The Netherlands: Elsevier. Google Scholar

• 72.
  Wissocq G

  . 2019Etudes des méthodes lattice Boltzmann pour les simulations de systèmes d'air secondaires de turbomachines. PhD thesis, Aix-Marseille Universite, France. Google Scholar

• 73.
  Hosseini SA, Darabiha N, Thévenin D, Eshghinejadfard A

  . 2017 Stability limits of the single relaxation-time advection–diffusion lattice Boltzmann scheme. Int. J. Mod. Phys. C 28 , 1750141. (doi:10.1142/S0129183117501418) Crossref, ISI, Google Scholar

• 74.
  Shan X, Chen H

  . 2007 A general multiple-relaxation-time Boltzmann collision model. Int. J. Mod. Phys. C 18 , 635–643. (doi:10.1142/S0129183107010887) Crossref, ISI, Google Scholar

• 75.
  Dellar PJ

  . 2001 Bulk and shear viscosities in lattice Boltzmann equations. Phys. Rev. E 64 , 031203. (doi:10.1103/PhysRevE.64.031203) Crossref, ISI, Google Scholar

• 76.
  Frapolli N, Chikatamarla SS, Karlin IV

  . 2016 Entropic lattice Boltzmann model for gas dynamics: theory, boundary conditions, and implementation. Phys. Rev. E 93 , 063302. (doi:10.1103/PhysRevE.93.063302) Crossref, PubMed, ISI, Google Scholar

• 77.
  Hafez M, Wahba E

  . 2007 Simulations of viscous transonic flows over lifting airfoils and wings. Comput. Fluids 36 , 39–52. (doi:10.1016/j.compfluid.2005.07.002) Crossref, ISI, Google Scholar

• 78.
  Mohamad AA

  . 2011 Lattice Boltzmann method . Berlin, Germany: Springer. Crossref, Google Scholar

• 79.
  Krüger T, Kusumaatmaja H, Kuzmin A, Shardt O, Silva G, Viggen EM

  . 2017 The lattice Boltzmann method . Berlin, Germany: Springer International Publishing. Crossref, Google Scholar

Mathematical Processes 4.1d 4.1.f Number Operation 4.2.f Answers 2019

Source: https://royalsocietypublishing.org/doi/10.1098/rsta.2019.0399

Posted by: rountreeandestalmoss.blogspot.com

Share this post