Fokker-Planck Central Moment Lattice Boltzmann Method for Effective Simulations of Fluid Dynamics

Abstract

We present a new formulation of the central moment lattice Boltzmann (LB) method based on a minimal continuous Fokker-Planck (FP) kinetic model, originally proposed for stochastic diffusive-drift processes (e.g., Brownian dynamics), by adapting it as a collision model for the continuous Boltzmann equation (CBE) for fluid dynamics. The FP collision model has several desirable properties, including its ability to preserve the quadratic nonlinearity of the CBE, unlike that based on the common Bhatnagar-Gross-Krook model. Rather than using an equivalent Langevin equation as a proxy, we construct our approach by directly matching the changes in different discrete central moments independently supported by the lattice under collision to those given by the CBE under the FP-guided collision model. This can be interpreted as a new path for the collision process in terms of the relaxation of the various central moments to “equilibria”, which we term as the Markovian central moment attractors that depend on the products of the adjacent lower order moments and a diffusion coefficient tensor, thereby involving of a chain of attractors; effectively, the latter are nonlinear functions of not only the hydrodynamic variables, but also the non-conserved moments; the relaxation rates are based on scaling the drift coefficient by the order of the moment involved. The construction of the method in terms of the relevant central moments rather than via the drift and diffusion of the distribution functions directly in the velocity space facilitates its numerical implementation and analysis. We show its consistency to the Navier-Stokes equations via a Chapman-Enskog analysis and elucidate the choice of the diffusion coefficient based on the second order moments in accurately representing flows at relatively low viscosities or high Reynolds numbers. We will demonstrate the accuracy and robustness of our new central moment FP-LB formulation, termed as the FPC-LBM, using the D3Q27 lattice for simulations of a variety of flows, including wall-bounded turbulent flows. We show that the FPC-LBM is more stable than other existing LB schemes based on central moments, while avoiding numerical hyperviscosity effects in flow simulations at relatively very low physical fluid viscosities through a refinement to a model founded on kinetic theory.

1. Introduction

The lattice Boltzmann Method (LBM) is a numerical scheme generally used for simulating fluid flows and various associated physical phenomena [1], and has been shown to deal with complex boundaries and interfacial dynamics, as well as multiscale flows such as turbulence and particulate suspension flows, quite well [2,3,4]. The key idea in the LBM is that simplified models of the continuous Boltzmann equation [5], are employed to track distributions of fluid particles which are restricted to collide and stream along specific velocity directions designated by lattice links. More generally, the collision process is modeled as to relax the distributions towards an equilibrium state designated by attractors, that are typically defined by the Maxwell distribution. Those particle distributions are then shifted, or streamed, back along the same lattice links. In turn, the hydrodynamic variables are then recovered by taking the relevant statistical moments of those distributions.

Moreover, modeling the collision process in the LBM plays an important role in tuning the physics of the flow characteristics as well as in the numerical stability of the scheme itself. To that end, many different collision models have been developed with the overall goal of improving its ability to represent multi-physics effects or to increase the limits of its applicability in flow simulations that are prone to numerical instabilities. In many practical applications, fluid flows occur with relatively large Reynolds numbers or low viscosities, and thus numerical simulation of those flows requires a robust numerical method capable of dealing with the associated numerical artifacts that arise when using a less robust model.

The first and simplest collision model, known as the single relaxation time (SRT) collision model [6], paved the way for treating the collision step as a relaxation process, but it was also found to have sever limitations in terms of numerical stability. More specifically, the SRT model breaks down when applied to problems that involve large Reynolds numbers or very small fluid viscosities. The primary issue being that the distributions are relaxed at the same rate, which causes some of them to incur large numerical errors. As a response to these issues, the multiple relaxation model was developed, which performs the collision step in moment space instead of velocity space, and allows different moments of the distributions to relax at different rates [7,8]. The MRT model was found to significantly improve the numerical stability when compared to the SRT collision model. Next, a natural extension of the MRT collision model was developed, which performs the relaxation in a new space known as central moment space that can be described as moment space with a moving reference frame. This new development is commonly known as the central moment collision model or as the cascaded collision model, and the attractors are obtained by matching the discrete central moment equilibria supported by a lattice with the corresponding central moments of the continuous Maxwell distribution [9]. This new model again showed great improvements in numerical stability when compared to the MRT collision model.

When simulating flows with relatively low fluid viscosities, some of such LBM collision models involve disparities between the relaxation rates of the second and higher order moments, and as such, some amendments are required to avoid numerical hyper-viscosities. One ad hoc solution in this regard for the central moment family of LB schemes involves constructing the attractors for the higher order central moments by exploiting the factorization property of the continuous Maxwellian and using the products of the post-collision central moments of lower orders to drive the relaxation process, without relying on the Maxwellian by itself [10]. This scheme is known as the factorized LBM. Furthermore, another solution that deals with the hyper-viscosity artifact, which is much more sophisticated with additional attendant transformations involves relaxing cumulants in favor of central moments under collision [11]. In this work, we propose and investigate a new and alternative formulation of the central moment method LB based on the relaxation of central moments towards attractors defined by taking the central moments of the Fokker-Planck model for the Boltzmann collision integral and eliminates these artifacts while achieving highly stable flow simulations through adopting a different perspective and refinements from kinetic theory.

The Fokker-Planck (FP) equation is an important formulation in statistical mechanics originally constructed for describing the evolution of the distribution of particles, whose motion stems from the effects of stochastic forces [12,13,14,15,16,17,18], the prototypical example being the Brownian dynamics arising from the stochastic diffusive and drift processes [19,20,21,22,23]. Thus, the continuous-time FP equation represents the changes in the distribution function due to drift and diffusion in the phase space, and is equivalent to the Langevin equation which is a classical example of the stochastic differential equation containing contributions from deterministic and random terms [24]. It is also referred to as the Kolmogorov forward equation in the mathematical literature and can be derived from the master equation, a fundamental approach for describing probabilistically-driven phenomena involving Markovian stochastic processes—a differential form of the so-called Chapman-Kolmogorov equation, via the Kramers-Moyal expansion. As such, the FP equation has been generalized further and applied to solve a large class of problems [25,26,27], which includes its recent use in constructing reduced order models for complex systems represented as stochastic cascade processes from empirical data [28,29].

In the context of fluid dynamics, based on stochastic molecular models pioneered by Kirkwood and collaborators [30,31,32,33], which are essentially Langevin equations and hence correspond to appropriate FP equations, the equations of fluid motion and their transport coefficients can be derived systematically [34]. On the other hand, of specific interest to this work, Lebowitz et al. [35] proposed the FP formulation, viz., the drift-and-diffusion driven changes of a distribution function in the phase space, as a model for the collision term in the continuous Boltzmann equation, where their drift and diffusion coefficients are related to the transport coefficients of the fluid. Moreover, Cercignani [36] demonstrated that such a FP model for collision represents a consistent and rigorous approximation (or a diffusion limit) of the Boltzmann’s collision integral term, especially for grazing collisions (see also [37,38,39,40,41]). Also, it has been shown by Pawula [42] that Boltzmann’s collision term leads to a differential operator under a finite truncation via the Kramers-Moyal expansion, involving only the first two terms corresponding to the FP equation. As such, the continuous Boltzmann equation itself can be derived in a few different ways, the most general approach among them is based on the so-called Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy [43]. Figure 1 shows schematically the different routes of how the celebrated Boltzmann equation can be derived under the repeated randomness or molecular chaos (stosszahlansatz) assumption [44] or via a reduction from the BBGKY hierarchy [43] and the modeling of its collision term via either the popular Bhatnagar-Gross-Krook (BGK) approximation [45] or the FP model under appropriate simplifying considerations; both these models reduces the complexity of using the Boltzmann equation in representing the dynamics of fluids (or plasmas). Clearly, Markov process considerations underlie the derivation of the Boltzmann equation or the development of the FP model, and is thus a crucial aspect in modeling the physical processes during collision.

In addition to the fact that the FP model directly arises as a simplification of the Boltzmann collision term, it has other advantages. The resulting FP collision model preserves the quadratic nonlinearity of the Boltzmann’s collision term, unlike the BGK model, which allows only a relaxation to the Maxwell distribution function (see [36]). Moreover, by exploiting the flexibility available in tuning the drift and diffusion coefficients appropriately and with necessary modifications, variety of complex fluid flows with attendant multiphysics effects can be modeled, such as flows involving heat transfer with adjustable Prandtl number, microscale gas flows, polyatomic gas flows, mixtures of gases, etc.

A stochastic numerical scheme for practical implementation, based on using an equivalent Langevin equation in place of the FP collision model in the Boltzmann equation was constructed by Jenny et al. in [46], which has been extended further and applied to various fluid dynamical applications more recently (see e.g., [47,48,49,50,51,52,53,54,55]). Generally, such algorithms serve as a computationally more efficient alternatives to the direct simulation Monte Carlo method for the solution of the Boltzmann equation for simulation of microscale gas flows at moderate Knudsen numbers. In the LB framework, Moroni et al. [56] presented a scheme for the solution of the FP equation whose discretization is effected by the Hermite projections in the velocity space. Such velocity space discretizations, which retain terms up to the third order in the fluid velocities for the changes in the distribution function under collision, was then extended for hydrodynamics in [57] but was found to have stability limitations, which will be addressed through a different approach in this work. More recently, the FP equation has been applied via its proxy, viz., an equivalent Langevin equation and then solved along with an LBM for simulating polymeric fluid flows (see e.g., [58,59,60]).

In this work, we present a new central moment LBM constructed from the simplified formulation of the continuous Boltzmann equation with the FP collision model via a matching principle, i.e., the changes in different discrete central moments independently supported by a chosen lattice of the LB algorithm under collision are equal to those given by the central moments of the continuous FP-based collision model of the Boltzmann equation for effective flow simulations. Thus, we propose to construct a LBM that implements the FP collision model in terms of their resulting system of central moments rather than relying on the solution of a stochastic proxy formulation based on Langevin-type equation, or trying to solve the FP equation directly as given in the velocity space. Such an approach is reminiscent of and has some analogy to formulating and solving the problem of dispersion of a tracer species under the combined action of molecular diffusion and the driving nonuniform velocity field—the so-called classical Taylor dispersion [61]; in such cases, it is known that solving the distribution of the tracer species directly depends on invoking various approximations and hence restrictive [61], while recasting the problem in terms of the evolution equations of the various moments of the distribution of the tracer species avoids the need to use restrictive assumptions and hence yields more general solutions as shown by Aris [62]. As such, in our case, the use of central moments avoids the need to approximate the cumbersome derivatives of the distribution functions with respect to the particle velocities present in the FP model, which can be represented exactly in terms of the central moments of lower orders, with an underlying local structure which is naturally amenable for highly stable LB implementations. Interestingly, as shown schematically in Figure 2, from a mathematical perspective, the Boltzmann’s collision term is an integral operator, while its modeling via a FP kinetic model leads to a differential operator, which when replaced with its central moments along with suitable refinements under the standard LB discretization leads a significantly simpler algebraic operator that greatly facilitates algorithmic implementation.

Our formulation for the LBM based on FP-guided collision results in a new path for the evolution of different central moments, especially at higher orders, under collision and can be re-interpreted in terms of the relaxation of the various central moments to “equilibria” that depend only on the adjacent, lower order post-collision moments, appearing in the form of recurrence relations. Since the collision process is fundamentally stochastic and Markovian in nature as noted above, we designate such equilibria identified under the equivalent relaxation-type interpretation, which also involves scaling the drift coefficient by the order of the participating moment to define the relaxation rates, as Markovian central moment attractors. This is distinct from the path of relaxation of central moments to the central moments of the Maxwellian under collision, since the Markovian central moment attractors depend on products of the lower order moments with components of a diffusion coefficient tensor, thereby allowing for considerations of their non-equilibrium effects. This is especially the case for evolving the fourth and higher order central moments under collision. As we will show in a later section, only in the special case where we assume the lower order moments appearing in the Markovian central moment attractors to be in equilibrium, will the latter coincide with those of the central moments of the Maxwellian. In addition, we will present and implement a method to include body forces in our formulation based directly on the Boltzmann’s acceleration term by invoking the matching principle, i.e., by setting the discrete central moment changes due to body forces to be equal to those given by the central moments of the acceleration term of the continuous Boltzmann equation.

We will construct central moment LBMs using FP-guided collision, referred to as the FPC-LBM, in two (2D) and three (3D) dimensions using the standard D2Q9 and D3Q27 lattices, respectively, and present the implementation details of the resulting algorithms. As such, the two-dimensional, nine velocity (D2Q9) and three-dimensional, twenty-seven velocity (D3Q27) lattices are the tensor product lattice sets in the respective dimensions and commonly used in LBM for flow simulations. In order to establish their consistency to the Navier-Stokes equations, we will present Chapman-Enskog analysis of the equivalent central moment system of the Boltzmann equation using the FP collision model in 2D and 3D. We will perform a detailed numerical study of our central moment LBM using FP-guided collision to demonstrate its accuracy and stability characteristics for various standard benchmark flow problems in 2D and 3D, and also present comparisons of its performance when compared to those using other collision models.

Finally, we note that for simulation of flows at relatively very low fluid viscosities, there is necessarily large disparities in the relaxation rates of the non-conserved second order and the higher order moments when the LB formulation is implemented using the multiple relaxation times, and if the central moments of the Maxwellian are used as equilibria, the contributions from such higher moments involves terms analogous to the non-equilibrium momentum fluxes (related to the strain rate tensor) and can dominate the corresponding physical contributions from the second-order non-equilibrium moments, which manifest as numerical hyperviscosities. On the other hand, this can be addressed by suitably prescribing the model parameters associated with the FP model, viz., the diffusion coefficient tensor in the Markovian central moment attractor. More specifically, the diffusion tensor parameter associated with the FP model is related to the variance of the distribution function, which corresponds to the different components of the second order central moments. When we consider the latter to be at equilibrium in evolving the conserved hydrodynamic variables via the relaxation of second order moments as required to recover the Navier-Stokes equations, while retaining them to be in general non-equilibrium states for evolving the higher moments, the resulting Markovian central moment attractors naturally become nonlinear functions of lower non-conserved moments, which, as shown later in this work, can effectively avoid such numerical hyperviscosities, while also greatly enhancing numerical stability in simulations of flows at higher Reynolds numbers. Furthermore, to demonstrate the capability of our approach for accurate turbulence simulations, we will presents results obtained using the FPC-LBM for turbulent channel flow and compare them with state-of-the-art direct numerical simulations (DNS) data obtained previously using a Navier-Stokes based solver.

In summary, our present work is motivated by the fact that various collision models of the LBM are prone to numerical instabilities or are subject to numerical hyperviscosities when simulating flows at relatively low fluid viscosities or high Reynolds numbers. We address both these issues by constructing FPC-LBM in 2D and 3D based on refining a model from kinetic theory that results in a more robust and accurate formulation for flow simulations in such cases.

This paper is organized as follows. In the next section (Section 2), we will present the derivation of a central moment formulation of the continuous FP collision operator of the Boltzmann equation in 2D, followed by the construction of the discrete Markovian central moment attractors using the D2Q9 lattice for the 2D FPC-LBM. Section 3 presents the corresponding 3D continuous central moment FP collision model and the respective discrete attractors using the D3Q27 lattice for the 3D FPC-LBM. Mathematical consistency analysis via the Chapman-Enskog expansions carried out directly using central moments in 2D and 3D are given in Appendix A and Appendix C, respectively; complete algorithmic implementation details of the 2D and 3D FPC-LBMs are discussed in Appendix B and Appendix D, respectively. Simulations using the FPC-LBMs for various standard benchmark flow problems that establish their accuracy, improvements in numerical stability when compared to using a variety of other existing collision models, and its ability in avoiding hyperviscosity effects are presented in Section 4. Finally, Section 5 briefly summarizes our main contributions and the key conclusions drawn from this work.

2. Construction and Analysis of 2D Central Moment Fokker-Planck Collision Model of Boltzmann Equation and 2D FPC-LBM

We start with the continuous Boltzmann equation given by (see e.g., [63])

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial f}{\partial t}}+\mathbf{\xi }·\mathbf{\nabla }f={\left({\displaystyle \frac{\delta f}{\delta t}}\right)}_{coll}^{Boltz}+{\left({\displaystyle \frac{\delta f}{\delta t}}\right)}_{forcing},\end{array}$$

(1)

where $f=f(\mathbf{x},\mathbf{\xi },t)$ is the single-particle density distribution function. In 2D, the position vector and the particle velocity are given by $\mathbf{x}=(x,y)$ and $\mathbf{\xi }=({\xi }_x,{\xi }_y)$, respectively, with $\mathbf{\nabla }=({\partial }_x,{\partial }_y)$. Considering athermal flows for simplicity, the fluid density $\rho$ and bulk velocity $\mathit{u}=(u_x,u_y)$ are obtained through the leading zeroth and first moments of f with respect to $\mathbf{\xi }$ given by

$$\rho =\int fd\mathbf{\xi },\rho \mathit{u}=\int f\mathbf{\xi }d\mathbf{\xi }.$$

(2)

Here, the left side of Equation (1) accounts for the variation in f due to particle advection or streaming, and if the fluid is subjected to an external force $\mathit{F}=(F_x,F_y)$, the particles undergo acceleration whose effect on f is represented through the term ${\left({\displaystyle \frac{\delta f}{\delta t}}\right)}_{forcing}$ in its right side given by

$${\left({\displaystyle \frac{\delta f}{\delta t}}\right)}_{forcing}=-{\displaystyle \frac{\mathit{F}}{\rho }}·{\mathbf{\nabla }}_{\xi }f,$$

(3)

where ${\mathbf{\nabla }}_{\xi }=({\partial }_{{\xi }_x},{\partial }_{{\xi }_y})$ is the gradient in the velocity space. The particles interact through collision that redistributes their velocities, which is accounted for on the rate of change in f via the Boltzmann’s collision term for binary collisions ${\left({\displaystyle \frac{\delta f}{\delta t}}\right)}_{coll}$ given by [63]

$$\begin{array}{c}\hfill {\left({\displaystyle \frac{\delta f}{\delta t}}\right)}_{coll}^{Boltz}=\int gb(f^{\prime }{f}_1^{\prime }-f{f}_1)\mathrm{d}b\mathrm{d}\epsilon \mathrm{d}{\xi }_i,\end{array}$$

(4)

where the distribution functions of particle pairs before and after collision are given by $f\equiv f(\mathbf{x},\mathbf{\xi },t)$, $f_1\equiv f(\mathbf{x},{\mathbf{\xi }}_1,t)$, and $f^{\prime }\equiv f(\mathbf{x},{\mathbf{\xi }}^{\prime },t)$, $f_1^{\prime }\equiv f(\mathbf{x},{\mathbf{\xi }}_1^{\prime },t)$. Here, $g=|\mathbf{\xi }-{\mathbf{\xi }}_1|$ is the relative speed, $\epsilon$ is the azimuthal collision angle and b is the impact parameter. Equation (4) is a multidimensional integral operator and represents one of the main difficulties with using the Boltzmann equation. Thus, it is often replaced by simpler models that represent certain essential features of the collision process. Varieties of collision models have been constructed over the years. One such popular model is the BGK model [45] given by

$$\begin{array}{c}\hfill {\left({\displaystyle \frac{\delta f}{\delta t}}\right)}_{coll}^{BGK}={\omega }_{BGK}(f^M-f),\end{array}$$

(5)

where $f^M$ is the local Maxwell distribution function given by

$$\begin{array}{c}\hfill f^M={\displaystyle \frac{\rho }{{\left(2\pi c_s^2\right)}^{\frac{D}{2}}}}exp\left[{\displaystyle \frac{-{(\mathbf{\xi }-\mathit{u})}^2}{2c_s^2}}\right],\end{array}$$

(6)

and ${\omega }_{BGK}$ is the collision frequency. In Equation (6), $\rho$ is the density, $\mathit{u}$ is the fluid velocity, $c_s$ is the speed of sound (generally related to temperature, but for athermal flows held as a constant), and D is the number of spatial dimensions. Equation (5) represents the effect of collisions as a relaxation toward an equilibrium state defined in Equation (6), and can be derived as a drastic approximation of the collision term (Equation (4)).

Another simplification of the Boltzmann’s collision integral (Equation (4)) is the Fokker-Planck (FP) model proposed in Lebowitz (1960) [35]. While originally used to represent the drift and diffusion processes of a Brownian particle, it can be shown to a rigorous approximation of the collision term (Equation (4)) in the Boltzmann equation under the assumptions of the small velocity changes under collision and the particles undergo grazing collision, which is often the case for scattering processes involving inverse power interaction potentials that result in large impact parameters (see Cercignani [36], Gombosi [40] and Liboff [41]). The Fokker-Planck collision model can be represented as [35]

$$\begin{array}{c}\hfill {\left({\displaystyle \frac{\delta f}{\delta t}}\right)}_{coll}^{FP}={\omega }_{FP}\left[{\displaystyle \frac{\partial }{\partial {\xi }_i}}\left(({\xi }_i-u_i)f\right)+D_{ij}^{\prime }{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_i\partial {\xi }_j}}\right],\end{array}$$

(7)

where ${\omega }_{FP}$ is a frequency characterizing collisions and $D_{ij}^{\prime }$ is the diffusion tensor parameter related to the variance of the distribution function or its second order central moment, which will be discussed in more detail later. The effective diffusion rate coefficient $D_{ij}$ itself is related to the diffusion tensor parameter $D_{ij}^{\prime }$ via the FP relaxation frequency as $D_{ij}=D_{ij}^{\prime }{\omega }_{FP}$. Equation (7) represents the effect of collision in terms of drift and diffusion in the velocity space. The drift term models a process analogous to the dynamical friction and reflects the idea that collisions tend to gradually eliminate all gradients in the particle velocities of the fluid; the diffusion term represents the effect of diffusion of f in the velocity space resulting in a broadening of the distribution function via molecular chaos [40]. It is interesting to note that unlike the BGK model (Equation (5)), the FP model (Equation (7)) preserves the quadratic nonlinearity of the Boltzmann’s collision integral (Equation (4)) since $u_i$ and $D_{ij}^{\prime }$ appearing in Equation (7) could themselves depend on f. We note here that Equation (7) is a generalization of the FP model developed in [35] by considering a tensor for the diffusion parameter $D_{ij}^{\prime }$ rather than a scalar $D^{\prime }$.

We will now formulate a central moment representation of the continuous FP model (Equation (7)), which will then serve as a basis for constructing a corresponding discrete central moment lattice Boltzmann (LB) method for flow simulations. We start by defining the inner products of any two objects a and b as the integration over velocity space in 2D as

$$\begin{array}{c}\hfill 〈a,b〉=\underset{-\infty }{\overset{\infty }{\int }}\underset{-\infty }{\overset{\infty }{\int }}ab\mathrm{d}{\xi }_x\mathrm{d}{\xi }_y.\end{array}$$

Now, we introduce the weights $W_{mn}$ of order ($m+n$) based on the integral powers of the components of the peculiar velocity ${\xi }_i-u_i$ (i.e., the particle velocity relative to the bulk fluid velocity) as

$$\begin{array}{c}\hfill W_{mn}={({\xi }_x-u_x)}^m{({\xi }_y-u_y)}^n.\end{array}$$

(8)

Then, we can define the central moment of the distribution function f of order ($m+n$) as the inner product of f and $W_{mn}$ as

$$\begin{array}{c}\hfill {\Pi }_{mn}=〈f,W_{mn}〉.\end{array}$$

For convenience, we rewrite the FP collision operator given in Equation (7) in a shorthand notation as

$$\begin{array}{c}\hfill {\left({\displaystyle \frac{\delta f}{\delta t}}\right)}_{coll}^{FP}={\omega }_{FP}\left({{\displaystyle \frac{\delta f}{\delta t}}}^{FP1}+{{\displaystyle \frac{\delta f}{\delta t}}}^{FP2}\right),\end{array}$$

(9)

where

$$\begin{array}{c}\hfill {{\displaystyle \frac{\delta f}{\delta t}}}^{FP1}={\displaystyle \frac{\partial }{\partial {\xi }_i}}\left(({\xi }_i-u_i)f\right),{{\displaystyle \frac{\delta f}{\delta t}}}^{FP2}=D_{ij}^{\prime }{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_i\partial {\xi }_j}}=D_{xx}^{\prime }{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_x^2}}+D_{yy}^{\prime }{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_y^2}}+2D_{xy}^{\prime }{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_x\partial {\xi }_y}}.\end{array}$$

(10)

Here, FP1 represents the `drift’ term of the FP collision operator while FP2 corresponds to the `diffusion term’.

We then define the rate of change of the central moment of order ($m+n$) under collision via the FP model, by taking the inner product of each term in the right hand side of Equation (9) along with the attendant terms in Equation (10), with the central moment weights of order $(m+n)$ given in Equation (8), as in

$$\begin{array}{c}\hfill {\left({\displaystyle \frac{\delta {\Pi }_{mn}}{\delta t}}\right)}_{coll}^{FP}={\omega }_{FP}\left[〈{{\displaystyle \frac{\delta f}{\delta t}}}^{FP1},W_{mn}〉+〈{{\displaystyle \frac{\delta f}{\delta t}}}^{FP2},W_{mn}〉\right]={\omega }_{FP}\left[{\displaystyle \frac{\delta {\Pi }_{mn}^{FP1}}{\delta t}}+{\displaystyle \frac{\delta {\Pi }_{mn}^{FP2}}{\delta t}}\right].\end{array}$$

(11)

Next, we take these inner products in Equation (11) term by term, starting with the first term

$$\begin{array}{c}\hfill {\displaystyle \frac{\delta {\Pi }_{mn}^{FP1}}{\delta t}}=〈{\displaystyle \frac{\partial }{\partial {\xi }_x}}\left(({\xi }_x-u_x)f\right),W_{mn}〉+〈{\displaystyle \frac{\partial }{\partial {\xi }_y}}\left(({\xi }_y-u_y)f\right),W_{mn}〉.\end{array}$$

(12)

To simplify this last term, consider the following chain rule ${\displaystyle \frac{\partial }{\partial {\xi }_i}}\left(({\xi }_i-u_i)f\right)=f+({\xi }_i-u_i){\displaystyle \frac{\partial f}{\partial {\xi }_i}}$ and then using this to rewrite Equation (12), and absorbing the $({\xi }_i-u_i)$ term into the respective central moment weight, $W_{mn}$, as

$$\begin{array}{c}\hfill {\displaystyle \frac{\delta {\Pi }_{mn}^{FP1}}{\delta t}}=〈{\displaystyle \frac{\partial f}{\partial {\xi }_x}},W_{m+1,n}〉+〈{\displaystyle \frac{\partial f}{\partial {\xi }_y}},W_{m,n+1}〉+2〈f,W_{mn}〉.\end{array}$$

(13)

We now further evaluate the right hand side of Equation (13) term by term, starting with the first term and rewriting it as

$$\begin{array}{c}\hfill 〈{\displaystyle \frac{\partial f}{\partial {\xi }_x}},W_{m+1,n,p}〉=\left(\underset{-\infty }{\overset{\infty }{\int }}{({\xi }_x-u_x)}^{m+1}{\displaystyle \frac{\partial f}{\partial {\xi }_x}}\mathrm{d}{\xi }_x\right)\left(\underset{-\infty }{\overset{\infty }{\int }}{({\xi }_y-u_y)}^n\mathrm{d}{\xi }_y\right),\end{array}$$

(14)

where the first term in the right hand side of Equation (14) can be evaluated via integration by parts as

$$\begin{array}{c}\hfill \underset{-\infty }{\overset{\infty }{\int }}{({\xi }_x-u_x)}^{m+1}{\displaystyle \frac{\partial f}{\partial {\xi }_x}}\mathrm{d}{\xi }_x={({\xi }_x-u_x)}^{m+1}f{|}_{-\infty }^{\infty }-\underset{-\infty }{\overset{\infty }{\int }}(m+1){({\xi }_x-u_x)}^{m}f\mathrm{d}{\xi }_x.\end{array}$$

Assuming that f decays faster than the increase due to any polynomials of ${\xi }_x-u_x$, i.e., ${({\xi }_x-u_x)}^{m+1}f=0$ as ${\xi }_x\to \pm \infty$, the first term on the right hand side of Equation (13) becomes

$$\begin{array}{c}\hfill 〈{\displaystyle \frac{\partial f}{\partial {\xi }_x}},W_{m+1,n}〉=-(m+1)\left(\underset{-\infty }{\overset{\infty }{\int }}{({\xi }_x-u_x)}^{m}f\mathrm{d}{\xi }_x\right)\left(\underset{-\infty }{\overset{\infty }{\int }}{({\xi }_y-u_y)}^n\mathrm{d}{\xi }_y\right),\end{array}$$

which is simplified to

$$\begin{array}{c}\hfill 〈{\displaystyle \frac{\partial f}{\partial {\xi }_x}},W_{m+1,n}〉=-(m+1)〈f,W_{mn}〉=-(m+1){\Pi }_{mn}.\end{array}$$

(15)

Similarly, the second term on the right hand side of Equation (13) is evaluated the same way as

$$\begin{array}{c}\hfill 〈{\displaystyle \frac{\partial f}{\partial {\xi }_y}},W_{m,n+1}〉=-(n+1){\Pi }_{mn}.\end{array}$$

(16)

Finally, using Equations (15) and (16) back into Equation (13) and simplifying, we obtain the following expression for the rate of change in the central moment of order $(m+n)$ under collision due to the drift term in the FP collision model as

$$\begin{array}{c}\hfill {\displaystyle \frac{\delta {\Pi }_{mn}^{FP1}}{\delta t}}=-(m+n){\Pi }_{mn}.\end{array}$$

(17)

Next, we consider the rate of change of central moment under collision due to the diffusion term in the FP collision model by evaluating the second term on the right hand side of Equation (11) and using Equation (10) as

$$\begin{array}{c}\hfill {\displaystyle \frac{\delta {\Pi }_{mn}^{FP2}}{\delta t}}=〈D_{xx}^{\prime }{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_x^2}},W_{mn}〉+〈D_{yy}^{\prime }{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_y^2}},W_{mn}〉+〈2D_{xy}^{\prime }{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_x\partial {\xi }_y}},W_{mn}〉.\end{array}$$

(18)

The first term on the right hand side of Equation (18) is evaluated as

$$\begin{array}{c}\hfill 〈D_{xx}^{\prime }{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_x^2}},W_{mn}〉=D_{xx}^{\prime }\left[\underset{-\infty }{\overset{\infty }{\int }}{({\xi }_x-u_x)}^m\mathrm{d}\left({\displaystyle \frac{\partial f}{\partial {\xi }_x}}\right)\right]\left(\underset{-\infty }{\overset{\infty }{\int }}{({\xi }_y-u_y)}^n\mathrm{d}{\xi }_y\right)\end{array}$$

(19)

where the term inside the square brackets on the right hand side of Equation (19) is readily evaluated via integration by parts as

$$\begin{array}{c}\hfill \underset{-\infty }{\overset{\infty }{\int }}{({\xi }_x-u_x)}^m\mathrm{d}\left({\displaystyle \frac{\partial f}{\partial {\xi }_x}}\right)={({\xi }_x-u_x)}^m{\displaystyle \frac{\partial f}{\partial {\xi }_x}}{|}_{-\infty }^{\infty }-\underset{-\infty }{\overset{\infty }{\int }}m{({\xi }_x-u_x)}^{m-1}\left({\displaystyle \frac{\partial f}{\partial {\xi }_x}}\right)\mathrm{d}{\xi }_x.\end{array}$$

As before, assuming ${({\xi }_x-u_x)}^m{\displaystyle \frac{\partial f}{\partial {\xi }_x}}=0$ as ${\xi }_x\to \pm \infty$, it follows that

$$\begin{array}{c}\hfill 〈D_{xx}^{\prime }{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_x^2}},W_{mn}〉=D_{xx}^{\prime }\left(-m\underset{-\infty }{\overset{\infty }{\int }}{({\xi }_x-u_x)}^{m-1}\left({\displaystyle \frac{\partial f}{\partial {\xi }_x}}\right)\mathrm{d}{\xi }_x\right)\left(\underset{-\infty }{\overset{\infty }{\int }}{({\xi }_y-u_y)}^n\mathrm{d}{\xi }_y\right),\end{array}$$

which then yields

$$\begin{array}{c}\hfill 〈D_{xx}^{\prime }{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_x^2}},W_{mn}〉=-m{D}_{xx}^{\prime }〈{\displaystyle \frac{\partial f}{\partial x}},W_{m-1,n}〉.\end{array}$$

(20)

Subsequently, invoking the identity in Equation (15), we get

$$\begin{array}{c}\hfill 〈{\displaystyle \frac{\partial f}{\partial {\xi }_x}},W_{m-1,n}〉=-(m-1)〈f,W_{m-2,n}〉=-(m-1){\Pi }_{m-2,n},\end{array}$$

and hence substituting this last equation in Equation (20) results in the following reduced expression for the latter

$$\begin{array}{c}\hfill 〈D_{xx}^{\prime }{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_x^2}},W_{mn}〉=m(m-1)D_{xx}^{\prime }{\Pi }_{m-2,n}.\end{array}$$

(21)

Similarly, the second and third terms in the right hand side of Equation (18) follow from the above considerations as

$$\begin{array}{c}\hfill 〈D_{yy}^{\prime }{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_y^2}},W_{mn}〉=n(n-1)D_{yy}^{\prime }{\Pi }_{m,n-2},〈D_{xy}^{\prime }{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_x\partial {\xi }_y}},W_{mn}〉=mn{D}_{xy}^{\prime }{\Pi }_{m-1,n-1}.\end{array}$$

(22)

Thus, the total rate of change of the central moment of order $(m+n)$ under collision due to the diffusion of the particle distribution function given in Equation (18) follows by combining the last two equations as

$$\begin{array}{c}\hfill {\displaystyle \frac{\delta {\Pi }_{mn}^{FP2}}{\delta t}}=m(m-1)D_{xx}^{\prime }{\Pi }_{m-2,n}+n(n-1)D_{yy}^{\prime }{\Pi }_{m,n-2}+2mn{D}_{xy}^{\prime }{\Pi }_{m-1,n-1}\end{array}$$

(23)

Finally, combining Equations (17) and (23), the net rate of change of central moment of order $(m+n)$ due to the drift and diffusion processes based on the Fokker-Planck collision model is given by

$$\begin{array}{ccc}\hfill {\left({\displaystyle \frac{\delta {\Pi }_{mn}}{\delta t}}\right)}_{coll}^{FP}& =& {\omega }_{FP}\left[-(m+n){\Pi }_{mn}+m(m-1)D_{xx}^{\prime }{\Pi }_{m-2,n}\right.\hfill \\ & & \left.+n(n-1)D_{yy}^{\prime }{\Pi }_{m,n-2}+2mn{D}_{xy}^{\prime }{\Pi }_{m-1,n-1}\right].\hfill \end{array}$$

Since ${\omega }_{FP}$ is a free parameter it can be rescaled by absorbing the prefactor $(m+n)$ of ${\Pi }_{mn}$ appearing inside the square bracket of the right side of the last equation and setting it as a renormalized relaxation frequency ${\omega }_{mn}$ given by ${\omega }_{mn}=(m+n){\omega }_{FP}$. In effect, we can then recast the central moment FP collision model as a relaxation of ${\Pi }_{mn}$ to its “equilibrium state" ${\Pi }_{mn}^{Mv}$, which we term as the Markovian central moment attractor of order $(m+n)$, at a rate ${\omega }_{mn}$. Such an attractor can be equivalently obtained from the stationary condition on the net effect of collision due to drift and diffusion given by

$$\begin{array}{ccc}\hfill \left[{\displaystyle \frac{\delta {\Pi }_{mn}^{FP1}}{\delta t}}+{\displaystyle \frac{\delta {\Pi }_{mn}^{FP2}}{\delta t}}\right]& =& 0when{\Pi }_{mn}={\Pi }_{mn}^{Mv}.\hfill \end{array}$$

(24)

Then, the refined central moment-based FP rate equation under collision can be finally expressed as

$$\begin{array}{c}\hfill {\left({\displaystyle \frac{\delta }{\delta t}}{\Pi }_{mn}\right)}_{coll}^{FP}={\omega }_{mn}\left[{\Pi }_{mn}^{Mv}-{\Pi }_{mn}\right],\end{array}$$

(25)

where

$$\begin{array}{c}\hfill {\Pi }_{mn}^{Mv}={\displaystyle \frac{m(m-1)}{m+n}}{D}_{xx}^{\prime }{\Pi }_{m-2,n}+{\displaystyle \frac{n(n-1)}{m+n}}{D}_{yy}^{\prime }{\Pi }_{m,n-2}+{\displaystyle \frac{2mn}{m+n}}{D}_{xy}^{\prime }{\Pi }_{m-1,n-1}.\end{array}$$

(26)

A few important remarks are in order here. Clearly, ${\Pi }_{mn}^{Mv}$ involves a recurrence relationship to evolve a higher central moment of order $(m+n)$ in terms of products of lower central moments of order $(m+n-2)$ and the diffusion tensor parameter $D_{ij}^{\prime }$. In other words, the attractor is a nonlinear function of not only the conserved hydrodynamic moments but also that of the various lower order non-conserved (kinetic) central moments since as discussed next, in general, the components of $D_{ij}^{\prime }$ themselves are related to the second order central moments. Moreover, as seen in Equation (25) and also emphasized earlier schematically in Figure 2, the mathematical character of the central moment formulation of the FP collision operator is algebraic and local rather than being differential or integral in nature, which lends itself to an effective numerical implementation using the LBM framework as discussed later in this section. Finally, since ${\omega }_{mn}$ can differ depending on the order of the moment, viz., $(m+n)$, the formulation given in Equations (25) and (26) naturally allows the use of multiple relaxation times or rates for relaxing different central moments to their respective Marvokian attractors. Such a generalization is analogous to adapting the original BGK model that uses a single relaxation time to using multiple relaxation times that relax different moments to their equilibria based on the Maxwell distribution as commonly adopted in previous LB approaches. However, the central moment FP collision model as developed here is conceptually different as it arises from modeling the drift-diffusion processes leading to a more general form for its “equilibria” as noted above.

2.1. Choice of Diffusion Tensor Parameter

The selection of the diffusion tensor parameter $D_{ij}^{\prime }$ is a key aspect in the numerical implementation of the central moment FP collision model. Before discussing it, let’s introduce the following convenient notation:

$$\begin{array}{c}\hfill D_{20}^{\prime }=D_{xx}^{\prime },D_{02}^{\prime }=D_{yy}^{\prime },D_{11}^{\prime }=D_{xy}^{\prime }.\end{array}$$

(27)

In general, $D_{ij}^{\prime }$ is related to the respective components of the variance (or the spread) of the distribution function f or its second order central moments. Now, if we evaluate the second order central moment attractors given in Equation (26) and use the notation in the last equation, we get ${\Pi }_{20}^{Mv}={\Pi }_{00}{D}_{20}^{\prime }$, ${\Pi }_{02}^{Mv}={\Pi }_{00}{D}_{02}^{\prime }$, and ${\Pi }_{11}^{Mv}={\Pi }_{00}{D}_{11}^{\prime }$, which upon inverting yield $D_{20}^{\prime }={\Pi }_{20}^{Mv}/{\Pi }_{00}$, $D_{02}^{\prime }={\Pi }_{02}^{Mv}/{\Pi }_{00}$, and $D_{11}^{\prime }={\Pi }_{11}^{Mv}/{\Pi }_{00}$. As shown in the Chapman-Enskog analysis of the continuous kinetic equation with the central moment FP collision model as developed above (see Appendix A) that for correctly recovering the conserved momentum fields, the leading components of second order moments ${\Pi }_{20}^{\left(0\right)}={\Pi }_{20}^{Mv}$, ${\Pi }_{02}^{\left(0\right)}={\Pi }_{02}^{Mv}$ and ${\Pi }_{11}^{\left(0\right)}={\Pi }_{11}^{Mv}$ should be isotropic with the diagonal parts related to the pressure field $P=\rho c_s^2$, i.e., ${\Pi }_{20}^{\left(0\right)}={\Pi }_{02}^{\left(0\right)}=\rho c_s^2$ and ${\Pi }_{11}^{\left(0\right)}=0$. Then, in view of the above, and noting that ${\Pi }_{00}=\rho$, it follows that for evolving the conserved modes, which in turn require the computation of the rate of change of second order moments, we require the diffusion coefficient tensor to be chosen based on the equilibrium second order central moments, i.e., $D_{20}^{\prime }={\Pi }_{20}^{\left(0\right)}/{\Pi }_{00}=c_s^2$, $D_{02}^{\prime }={\Pi }_{02}^{\left(0\right)}/{\Pi }_{00}=c_s^2$, and $D_{11}^{\prime }={\Pi }_{11}^{\left(0\right)}/{\Pi }_{00}=0$.

On the other hand, for the physically important and challenging case involving the simulation of high Reynolds number flows, which are associated with fluids with relatively low viscosities that depend on the relaxation frequencies of the second order moments (see Appendix A), we generally have $1/{\omega }_{II}\ll 1/{\omega }_h$, where ${\omega }_{II}$ is taken as the relaxation frequency of the second order moments (i.e., $(m+n)=2$) and ${\omega }_h$ correspond to that of higher order moments (i.e., $(m+n)>2$). In other words, the second order moments are fast modes while the higher order ones are slower to relax. In such situations, when evolving the higher moments, the fast second order moments would have already approached their post-collision states which are not in equilibrium but contain their non-equilibrium contributions as well, which should be accounted for in their associated Markovian attractors that depend on $D_{ij}^{\prime }$. Hence, in such situations, to relax the higher order moments with $(m+n)>2$ under collisions, we prescribe $D_{ij}^{\prime }$ to be equal to ${\Pi }_{ij}/{\Pi }_{00}$ or $D_{20}^{\prime }={\Pi }_{20}/\rho$, $D_{02}^{\prime }={\Pi }_{02}/\rho$, and $D_{11}^{\prime }={\Pi }_{11}/\rho$. One may interpret such modification to the diffusion tensor parameter for evolving the higher order moments when compared to the second order moments as a form of renormalization.

Then, in summary, the diffusion tensor parameter appearing in the Markovian central moment attractor in Equation (26) for evolving the moment ${\Pi }_{mn}$ of order $(m+n)$ under collision is selected as follows:

$$\begin{array}{ccc}& & D_{20}^{\prime }=c_s^2,D_{02}^{\prime }=c_s^2,D_{11}^{\prime }=0for(m+n)\le 2,\hfill \\ & & D_{20}^{\prime }={\displaystyle \frac{{\Pi }_{20}}{\rho }},D_{02}^{\prime }={\displaystyle \frac{{\Pi }_{02}}{\rho }},D_{11}^{\prime }={\displaystyle \frac{{\Pi }_{11}}{\rho }}for(m+n)>2.\hfill \end{array}$$

(28)

2.2. Selected Continuous Markovian Central Moment Attractors in 2D

From the above considerations, we can then evaluate the continuous Markovian central moment attractors given in Equation (26) for the components of moments that are independently supported by the D2Q9 lattice, which are required for the LB algorithm discussed in this section. Then, we get

$$\begin{array}{ccc}\hfill & & {\Pi }_{00}^{Mv}=\rho ,{\Pi }_{10}^{Mv}=0,{\Pi }_{01}^{Mv}=0,\hfill \\ & & {\Pi }_{20}^{Mv}=c_s^2\rho ,{\Pi }_{02}^{Mv}=c_s^2\rho ,\hfill \\ & & {\Pi }_{21}^{Mv}=0,{\Pi }_{12}^{Mv}=0,\hfill \\ & & {\Pi }_{22}^{Mv}={\displaystyle \frac{1}{\rho }}\left[{\Pi }_{20}{\Pi }_{02}+2{\Pi }_{11}^2\right].\hfill \end{array}$$

(29)

Interestingly, the attractor for the fourth order moment is now a nonlinear function of the second order central moments.

2.3. Central Moment of Boltzmann’s Acceleration Term Due to Body Force in 2D

In addition, where there is a body force, it results in a rate of change contribution to the distribution function through ${\left({\displaystyle \frac{\delta f}{\delta t}}\right)}_{forcing}$ in the continuous Boltzmann equation (Equation (1)). We define the central moment of the rate of change of the distribution function due to the body force of order ($m+n$) as

$$\begin{array}{c}\hfill {\Gamma }_{mn}=〈{\left({\displaystyle \frac{\delta f}{\delta t}}\right)}_{forcing},W_{mn}〉,\end{array}$$

(30)

which can be evaluated by using Equation (3) in Equation (30) to get

$$\begin{array}{c}\hfill {\Gamma }_{mn}=-{\displaystyle \frac{F_x}{\rho }}〈{\displaystyle \frac{\partial f}{\partial {\xi }_x}},W_{mn}〉-{\displaystyle \frac{F_y}{\rho }}〈{\displaystyle \frac{\partial f}{\partial {\xi }_y}},W_{mn}〉.\end{array}$$

(31)

Then, from Equations (15) and (16), it follows that $〈{\displaystyle \frac{\partial f}{\partial {\xi }_x}},W_{mn}〉=-m{\Pi }_{m-1,n}$ and $〈{\displaystyle \frac{\partial f}{\partial {\xi }_y}},W_{mn}〉=-n{\Pi }_{m,n-1}$. Hence, the rate of change of the central moment of order ($m+n$) due to the body force is given by

$$\begin{array}{c}\hfill {\Gamma }_{mn}=m{\displaystyle \frac{F_x}{\rho }}{\Pi }_{m-1,n}+n{\displaystyle \frac{F_y}{\rho }}{\Pi }_{m,n-1},\end{array}$$

(32)

which relates it to $\mathit{F}=(F_x,F_y)$ and to the central moment components of f of a lower order (reduced by a degree of 1) and is exact (see [64]).

A Chapman-Enskog analysis of the 2D central moment-based continuous Boltzmann equation with the above FP collision model that is directly based on an expansion in terms of the central moments to establish its consistency with the 2D Navier-Stokes equations is given in Appendix A.

2.4. Construction of 2D FPC-LBM

We will now discretize the Boltzmann equation (see Equation (1)) with using the FP guided collision operator given in Equation (25) along with the forcing term (see Equation (30)) as follows: first we replace the continuous particle velocity $\mathbf{\xi }$ with a discrete particle velocity set ${\mathbf{e}}_{\alpha }$, with each direction corresponding to a lattice link, such that we have $f_{\alpha }(\mathbf{x},t)=f(\mathbf{x},{\mathbf{e}}_{\alpha },t)$, and then integrate the resulting discrete velocity Boltzmann equation along the particle characteristics over a time step ${\delta }_t$ that exactly spans a distance equal to the magnitude of ${\mathbf{e}}_{\alpha }{\delta }_t$ which then further discretizes both the space and time (see He et al. [65,66]). In 2D, the standard D2Q9 lattice is used whose discrete particle velocities ${\mathbf{e}}_{\alpha }$ have the following Cartesian components:

$$\begin{array}{c}\hfill |{\mathbf{e}}_x\rangle ={(0,1,0,-1,0,1,-1,-1,1)}^{†},\\ \hfill |{\mathbf{e}}_y\rangle ={(0,0,1,0,-1,1,1,-1,-1)}^{†},\end{array}$$

where, henceforth, † denotes the transpose operator and we use the Dirac’s notation of $|·\rangle$ and $\langle ·|$ for denoting the column and row vectors, respectively. Then the resulting LBM can be generically represented as a two-step process involving a collision step which is then followed by a streaming step:

$$\begin{array}{}\mathrm{(33a)}& \hfill {\tilde{f}}_{\alpha }(\mathbf{x},t)& =& f_{\alpha }(\mathbf{x},t)+{\Omega }_{\alpha }(\mathbf{x},t),\alpha =0,1,2,\dots ,q\hfill \mathrm{(33b)}& f_{\alpha }(\mathbf{x},t+{\delta }_t)& =& {\tilde{f}}_{\alpha }(\mathbf{x}-{\mathbf{e}}_{\alpha }{\delta }_t,t),\hfill \end{array}$$

where $q=9$ and ${\Omega }_{\alpha }(\mathbf{x},t)$ represents the cumulative changes in the distribution function $f_{\alpha }$ due to collision and the effect of body force during a time step ${\delta }_t$. In practice, such changes will be obtained in terms of central moments via the FP-guided collision and source term updates derived for the continuous Boltzmann equation projected to the countable independent moments represented by the D2Q9 lattice, which will then be mapped back to those in terms of the distribution functions. To accomplish this, we first define the discrete central moments and discrete raw moments of order ($m+n$) of the distribution function $f_{\alpha }$, its corresponding Markovian attractor $f_{\alpha }^{Mv}$, and the source term $S_{\alpha }$ for the body force, respectively, as follows:

$$\begin{array}{}\mathrm{(34a)}& \hfill \left(\begin{array}{c}{\kappa }_{mn}\\ {\kappa }_{mn}^{Mv}\\ {\sigma }_{mn}\end{array}\right)& =& \sum _{\alpha =0}^8\left(\begin{array}{c}{f}_{\alpha }\\ f_{\alpha }^{Mv}\\ S_{\alpha }\end{array}\right){(e_{\alpha x}-u_x)}^m{(e_{\alpha y}-u_y)}^n,and\hfill \mathrm{(34b)}& \hfill \left(\begin{array}{c}{\kappa }_{mn}^{\prime }\\ {\kappa }_{mn}^{M{v}^{\prime }}\\ {\sigma }_{mn}^{\prime }\end{array}\right)& =& \sum _{\alpha =0}^8\left(\begin{array}{c}{f}_{\alpha }\\ f_{\alpha }^{Mv}\\ S_{\alpha }\end{array}\right)e_{\alpha x}^m{e}_{\alpha y}^n.\hfill \end{array}$$

Then, the countable independent central moments and raw moments for the D2Q9 lattice respectively are given by

$$\begin{array}{c}\hfill {\mathbf{m}}^{\mathbf{c}}=({\kappa }_{00},{\kappa }_{10},{\kappa }_{01},{\kappa }_{20},{\kappa }_{02},{\kappa }_{11},{\kappa }_{21},{\kappa }_{12},{\kappa }_{22}),\end{array}$$

(35a)

$$\begin{array}{c}\hfill \mathbf{m}=({\kappa }_{00}^{\prime },{\kappa }_{10}^{\prime },{\kappa }_{01}^{\prime },{\kappa }_{20}^{\prime },{\kappa }_{02}^{\prime },{\kappa }_{11}^{\prime },{\kappa }_{21}^{\prime },{\kappa }_{12}^{\prime },{\kappa }_{22}^{\prime }),\end{array}$$

(35b)

and, similarly, those for the attractors and sources can be listed. Let us also define a 9-dimensional vector $\mathbf{f}$ containing the distribution functions via

$$\begin{array}{c}\hfill \mathbf{f}={(f_0,f_1,f_2,\cdots ,f_8)}^{†},\end{array}$$

(36)

and one can similarly group the attractors $f_{\alpha }^{Mv}$ and sources $S_{\alpha }$ in respective vectors.

To perform the collision step in terms of relaxations of different central moments to their respective attractors along with source term updates, we need to first successively map the distribution functions $\mathbf{f}$ to raw moments $\mathbf{m}$ and then to central moments ${\mathbf{m}}^{\mathbf{c}}$; then the post-collision central moments need to be transformed back to raw moments and then to the distribution functions to complete the collision step. To achieve these, we introduce the mappings between the raw moments and the distribution functions as

$$\mathbf{m}=\mathsf{P}\mathbf{f},\mathbf{f}={\mathsf{P}}^{-1}\mathbf{m},$$

(37)

where $\mathsf{P}$ is a matrix dependent on the basis vectors for moments, i.e., via the components of the monomials $\left.〈{\mathbf{e}}_x^m{\mathbf{e}}_y^n\right|$ given by

$$\mathsf{P}=\left[\right|\mathbf{1}\rangle ,|{\mathbf{e}}_x\rangle ,|{\mathbf{e}}_y\rangle ,|{\mathbf{e}}_x^2\rangle ,|{\mathbf{e}}_y^2\rangle ,|{\mathbf{e}}_x{\mathbf{e}}_y\rangle ,|{\mathbf{e}}_x^2{\mathbf{e}}_y\rangle ,|{\mathbf{e}}_x{\mathbf{e}}_y^2\rangle ,|{\mathbf{e}}_x^2{\mathbf{e}}_y^2\rangle {]}^{†},$$

(38)

where $\left|\mathbf{1}\right.〉$ is a 9-dimensional unit vector given by $|\mathbf{1}\rangle =\left|{\mathbf{e}}_x^0{\mathbf{e}}_y^0\right.〉={(1,1,1,1,1,1,1,1,1)}^{†}$. In addition, we represent the transformations between central moments and raw moments as

$${\mathbf{m}}^{\mathbf{c}}=\mathsf{F}\mathbf{m},\mathbf{m}={\mathsf{F}}^{-1}{\mathbf{m}}^{\mathbf{c}},$$

(39)

where both the frame transformation matrix $\mathsf{F}$ and its inverse ${\mathsf{F}}^{-1}$ are lower triangular matrices dependent on the fluid velocity components $\mathit{u}=(u_x,u_y)$ and can be readily evaluated through collecting the coefficients of the binomial expansions and their inversions (see e.g., [67]).

Before we describe our 2D FPC-LBM, we first obtain the expressions for the discrete central moments of the Markovian attractor ${\kappa }_{mn}^{Mv}$ and the source term ${\sigma }_{mn}$ due to force field by matching the corresponding continuous versions as given in Equation (26) (or Equation (29)) and Equation (32), respectively, i.e., ${\kappa }_{mn}^{Mv}={\Pi }_{mn}^{Mv}$ and ${\sigma }_{mn}={\Gamma }_{mn}$. Thus, the discrete central moments of the Markovian attractors for the D2Q9 lattice read as

$$\begin{array}{c}{\kappa }_{00}^{Mv}=\rho ,{\kappa }_{10}^{Mv}=0,{\kappa }_{01}^{Mv}=0,\\ {\kappa }_{20}^{Mv}=c_s^2\rho ,{\kappa }_{02}^{Mv}=c_s^2\rho ,{\kappa }_{11}^{Mv}=0,\\ {\kappa }_{21}^{Mv}=0,{\kappa }_{12}^{Mv}=0,\\ \begin{array}{|c|}\hline {\kappa }_{22}^{Mv}={\displaystyle \frac{1}{\rho }}\left({\tilde{\kappa }}_{20}{\tilde{\kappa }}_{02}+2{\tilde{\kappa }}_{11}{\tilde{\kappa }}_{11}\right)\\ \hline\end{array}.\end{array}$$

(40)

It can be seen that the fourth order attractor depends on the sum of the products of the current states of the second order central moments. In practical implementations, we will consider the post-collision state of the central moments ${\kappa }_{20}$, ${\kappa }_{02}$, and ${\kappa }_{11}$ (designated by ${\tilde{\kappa }}_{20}$, ${\tilde{\kappa }}_{02}$, and ${\tilde{\kappa }}_{11}$) in evaluating the fourth order Markovian attractor. Moreover, the contributions of the body force to the various central moments can be listed as follows:

$$\begin{array}{c}{\sigma }_{00}=0,{\sigma }_{10}=F_x,{\sigma }_{01}=F_y,\\ {\sigma }_{20}=2{\displaystyle \frac{F_x}{\rho }}{\kappa }_{10}=0,{\sigma }_{02}=2{\displaystyle \frac{F_y}{\rho }}{\kappa }_{01}=0,{\sigma }_{11}={\displaystyle \frac{F_x}{\rho }}{\kappa }_{10}+{\displaystyle \frac{F_y}{\rho }}{\kappa }_{01}=0,\\ {\sigma }_{21}=2{\displaystyle \frac{F_x}{\rho }}{\kappa }_{11}+{\displaystyle \frac{F_y}{\rho }}{\kappa }_{20},{\sigma }_{12}={\displaystyle \frac{F_x}{\rho }}{\kappa }_{02}+{\displaystyle \frac{F_y}{\rho }}{\kappa }_{11},\\ {\sigma }_{22}=2\left({\displaystyle \frac{F_x}{\rho }}{\kappa }_{12}+{\displaystyle \frac{F_y}{\rho }}{\kappa }_{21}\right).\end{array}$$

(41)

The complete algorithmic details of the implementation of the 2D FPC-LBM using the D2Q9 lattice are given in Appendix B. The solution method can be briefly summarized as follows. During every time step, the discrete distribution functions are first mapped to a set of raw moments and then to central moments. Then, the collision step is performed by relaxing various central moments to their equilibria based on the Markovian attractors of the Fokker-Planck collision model discussed earlier. Once the post-collision central moments are computed, they are transformed back to corresponding raw moments and then to distribution functions. Subsequently, the streaming step is performed by lock-step advection of the post-collision distribution functions along the respective lattice directions. Finally, using these updated distribution functions, the hydrodynamic fields, including the local density and velocity fields are obtained via taking their kinetic moments, which completes one time step of the FPC-LBM.

3. Construction and Analysis of 3D Central Moment Fokker-Planck Collision Model of Boltzmann Equation and 3D FPC-LBM

We will now extend our formulation discussed in the previous section to 3D. Here, we consider the Boltzmann equation, its acceleration term, and the FP collision operator given in Equations (1), (3), and (7), respectively, by taking the phase space coordinates, their derivatives, and the body force as $\mathbf{x}=(x,y,z)$, $\mathbf{\xi }=({\xi }_x,{\xi }_y,{\xi }_z)$, $\mathbf{\nabla }=({\partial }_x,{\partial }_y,{\partial }_z)$, ${\mathbf{\nabla }}_{\xi }=({\partial }_{{\xi }_x},{\partial }_{{\xi }_y},{\partial }_{{\xi }_z})$, and $\mathit{F}=(F_x,F_y,F_z)$. Then, to obtain the rate of change of central moments of f due to collision guided by the Fokker-Planck model in 3D, we follow the same procedure as before, by taking the inner product of each term in the Fokker-Planck collision model with the central moment weights. Here, the inner product of two objects, a and b, is now defined in 3D as

$$\begin{array}{c}\hfill 〈a,b〉=\underset{-\infty }{\overset{\infty }{\int }}\underset{-\infty }{\overset{\infty }{\int }}\underset{-\infty }{\overset{\infty }{\int }}ab\mathrm{d}{\xi }_x\mathrm{d}{\xi }_y\mathrm{d}{\xi }_z\end{array}$$

(42)

and the corresponding central moment weights of order ($m+n+p$) as

$$\begin{array}{c}\hfill W_{mnp}={({\xi }_x-u_x)}^m{({\xi }_y-u_y)}^n{({\xi }_z-u_z)}^p.\end{array}$$

(43)

Then, the continuous central moment of order ($m+n+p$) is written as

$$\begin{array}{c}\hfill {\Pi }_{mnp}=〈f,W_{mnp}〉.\end{array}$$

(44)

As before, we rewrite the FP collision operator given in Equation (7) conveniently in a shorthand notation as shown in Equation (9), where ${{\displaystyle \frac{\delta f}{\delta t}}}^{FP1}$ is given in Equation (10) and ${{\displaystyle \frac{\delta f}{\delta t}}}^{FP2}$ now modifies in 3D to

$$\begin{array}{c}\hfill {{\displaystyle \frac{\delta f}{\delta t}}}^{FP2}=D_{xx}^{\prime }{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_x^2}}+D_{yy}^{\prime }{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_y^2}}+D_{zz}^{\prime }{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_z^2}}+2D_{xy}^{\prime }{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_x\partial {\xi }_y}}+2D_{xz}^{\prime }{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_x\partial {\xi }_z}}+2D_{yz}^{\prime }{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_y\partial {\xi }_z}},\end{array}$$

(45)

which expresses the total rate of change of f under collision due to its diffusion in different possible directions in the velocity space, and whose magnitude is determined by diffusion tensor parameter $D_{ij}^{\prime }$. For ease of presentation, we rewrite the latter by following a slightly more general and compact notation:

$$\begin{array}{c}\hfill D_{200}^{\prime }=D_{xx}^{\prime },D_{020}^{\prime }=D_{yy}^{\prime },D_{002}^{\prime }=D_{zz}^{\prime },D_{110}^{\prime }=D_{xy}^{\prime },D_{011}^{\prime }=D_{yz}^{\prime },D_{101}^{\prime }=D_{xz}^{\prime }.\end{array}$$

(46)

Then, as in Equation (11), we define the rate of change of the central moment of order ($m+n+p$) under collision via the FP model, by taking the inner product of each term in the right hand side of Equation (9) using the attendant terms ${{\displaystyle \frac{\delta f}{\delta t}}}^{FP1}$ in Equation (10) and ${{\displaystyle \frac{\delta f}{\delta t}}}^{FP2}$ in Equation (45) with the central moment weights of order $(m+n+p)$ given in Equation (43) as

$$\begin{array}{c}\hfill {\left({\displaystyle \frac{\delta {\Pi }_{mnp}}{\delta t}}\right)}_{coll}^{FP}={\omega }_{FP}\left[〈{{\displaystyle \frac{\delta f}{\delta t}}}^{FP1},W_{mnp}〉+〈{{\displaystyle \frac{\delta f}{\delta t}}}^{FP2},W_{mnp}〉\right]={\omega }_{FP}\left[{\displaystyle \frac{\delta {\Pi }_{mnp}^{FP1}}{\delta t}}+{\displaystyle \frac{\delta {\Pi }_{mnp}^{FP2}}{\delta t}}\right].\end{array}$$

(47)

Now, following the derivation given in the previous section, the central moments of the derivatives of the distribution function in different directions in the velocity space can be expressed as follows:

$$\begin{array}{c}\hfill 〈{\displaystyle \frac{\partial f}{\partial {\xi }_x}},W_{mnp}〉=-m{\Pi }_{m-1,n,p},〈{\displaystyle \frac{\partial f}{\partial {\xi }_y}},W_{mnp}〉=-n{\Pi }_{m,n-1,p},〈{\displaystyle \frac{\partial f}{\partial {\xi }_z}},W_{mnp}〉=-p{\Pi }_{m,n,p-1}.\end{array}$$

(48)

Using these successively in pairs, we get the central moments of the cross-diffusion terms in the velocity space as

$$\begin{array}{ccc}\hfill 〈{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_x{\xi }_y}},W_{mnp}〉& =& mn{\Pi }_{m-1,n-1,p},〈{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_y{\xi }_z}},W_{mnp}〉=np{\Pi }_{m,n-1,p-1},\hfill \\ \hfill 〈{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_x{\xi }_z}},W_{mnp}〉& =& mp{\Pi }_{m-1,n,p-1}.\hfill \end{array}$$

(49)

Also, the central moments of the diffusion terms in the principal directions in the velocity space read as

$$\begin{array}{ccc}\hfill 〈{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_x^2}},W_{mnp}〉& =& m(m-1){\Pi }_{m-2,n,p},〈{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_y^2}},W_{mnp}〉=n(n-1){\Pi }_{m,n-2,p},\hfill \\ \hfill 〈{\displaystyle \frac{{\partial }^{2}f}{\partial {\xi }_z^2}},W_{mnp}〉& =& p(p-1){\Pi }_{m,n,p-2},\hfill \end{array}$$

(50)

Using the above identities, the drift and diffusion contributions to the rate of change of central moment of order ($m+n+p$) under collision can be written as follows:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\delta {\Pi }_{mnp}^{FP1}}{\delta t}}& =& 〈{\displaystyle \frac{\partial f}{\partial {\xi }_x}},W_{m+1,n,p}〉+〈{\displaystyle \frac{\partial f}{\partial {\xi }_y}},W_{m,n+1,p}〉+〈{\displaystyle \frac{\partial f}{\partial {\xi }_z}},W_{m,n,p+1}〉+3〈f,W_{mnp}〉\hfill \\ & =& -(m+n+p){\Pi }_{mnp}.\hfill \end{array}$$

(51)

and

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\delta {\Pi }_{mnp}^{FP2}}{\delta t}}& =& m(m-1)D_{200}^{\prime }{\Pi }_{m-2,n,p}+n(n-1)D_{020}^{\prime }{\Pi }_{m,n-2,p}+p(p-1)D_{002}^{\prime }{\Pi }_{m,n,p-2}+\hfill \\ & & 2mn{D}_{110}^{\prime }{\Pi }_{m-1,n-1,p}+2np{D}_{011}^{\prime }{\Pi }_{m,n-1,p-1}+2mp{D}_{101}^{\prime }{\Pi }_{m-1,n,p-1}.\hfill \end{array}$$

(52)

As such, from these last two equations, it can be seen that the drift represents a depletion process, while the diffusion is a gain process under collision. Combining them and since the characteristic Fokker-Planck inverse time scale ${\omega }_{FP}$ is a free parameter, we can rescale it using $(m+n+p)$ as ${\omega }_{mnp}=(m+n+p){\omega }_{FP}$, such that it can then be used to effectively control the rate of relaxation of the central moment ${\Pi }_{mnp}$ to its attractor designated as ${\Pi }_{mnp}^{Mv}$. Then, the central moment FP collision operator can be compactly rewritten in 3D as

$${\left({\displaystyle \frac{\delta }{\delta t}}{\Pi }_{mnp}\right)}_{coll}^{FP}={\omega }_{mnp}\left[{\Pi }_{mnp}^{Mv}-{\Pi }_{mnp}\right],$$

(53)

where the attendant Markovian central moment attractor reads as

$$\begin{array}{ccc}\hfill {\Pi }_{mnp}^{Mv}& =& D_{200}^{\prime }{\displaystyle \frac{m(m-1)}{(m+n+p)}}{\Pi }_{m-2,n,p}+D_{020}^{\prime }{\displaystyle \frac{n(n-1)}{(m+n+p)}}{\Pi }_{m,n-2,p}+\hfill \\ & & D_{002}^{\prime }{\displaystyle \frac{p(p-1)}{(m+n+p)}}{\Pi }_{m,n,p-2}+2D_{110}^{\prime }{\displaystyle \frac{mn}{(m+n+p)}}{\Pi }_{m-1,n-1,p}+\hfill \\ & & 2D_{011}^{\prime }{\displaystyle \frac{np}{(m+n+p)}}{\Pi }_{m,n-1,p-1}+2D_{101}^{\prime }{\displaystyle \frac{mp}{(m+n+p)}}{\Pi }_{m-1,n,p-1}.\hfill \end{array}$$

(54)

We emphasize that all the key remarks made regarding the Markovian central moment attractor for the 2D case in the paragraph below Equation (26) is applicable for the above Equation (54) as well by replacing ($m+n$) with ($m+n+p$) along with Equation (53) by replacing ${\omega }_{mn}$ with ${\omega }_{mnp}$.

3.1. Choice of Diffusion Tensor Parameter

The choice of the diffusion tensor parameter $D_{ij}^{\prime }$, which may also be referred to as the second order Kramers-Moyal expansion coefficients [42], plays a crucial role in the performance of the numerical simulations that use the central moment FP collision formulation. A detailed discussion of the considerations involved in their selection in the 2D case is given earlier in Section 2.1, which also readily extends for the present 3D model given above. Following this, for evolving moments of second order (i.e., $(m+n+p)=2$) that is used to update the hydrodynamic fields, we require $D_{200}^{\prime }={\Pi }_{200}^{\left(0\right)}/{\Pi }_{000}^{\left(0\right)}$, $D_{020}^{\prime }={\Pi }_{020}^{\left(0\right)}/{\Pi }_{000}^{\left(0\right)}$, $D_{002}^{\prime }={\Pi }_{002}^{\left(0\right)}/{\Pi }_{000}^{\left(0\right)}$, $D_{110}^{\prime }={\Pi }_{110}^{\left(0\right)}/{\Pi }_{000}^{\left(0\right)}$, $D_{101}^{\prime }={\Pi }_{101}^{\left(0\right)}/{\Pi }_{000}^{\left(0\right)}$, and $D_{011}^{\prime }={\Pi }_{011}^{\left(0\right)}/{\Pi }_{000}^{\left(0\right)}$, where ${\Pi }_{mnp}^{\left(0\right)}={\Pi }_{mnp}^{Mv}=c_s^2\rho$ for $\left(mnp\right)=\left(200\right),\left(020\right)$ and $\left(002\right)$, and ${\Pi }_{mnp}^{\left(0\right)}={\Pi }_{mnp}^{Mv}=0$ for $\left(mnp\right)=\left(110\right),\left(101\right)$ and $\left(011\right)$; on the other hand, for evolving higher moments (i.e., $(m+n+p)>2$), we prescribe $D_{200}^{\prime }={\Pi }_{200}/{\Pi }_{000}$, $D_{020}^{\prime }={\Pi }_{020}/{\Pi }_{000}$, $D_{002}^{\prime }={\Pi }_{002}/{\Pi }_{000}$, $D_{110}^{\prime }={\Pi }_{110}/{\Pi }_{000}$, $D_{101}^{\prime }={\Pi }_{101}/{\Pi }_{000}$, and $D_{011}^{\prime }={\Pi }_{011}/{\Pi }_{000}$. Then, using ${\Pi }_{000}={\Pi }_{000}^{\left(0\right)}=\rho$, in summary, the diffusion tensor parameter appearing in the Markovian central moment attractor in Equation (54) for evolving the moment ${\Pi }_{mnp}$ of order $(m+n+p)$ under collision is selected as follows:

$$\begin{array}{ccc}\hfill & & D_{200}^{\prime }=D_{020}^{\prime }=D_{002}^{\prime }=c_s^2,\mathrm{for}(m+n+p)\le 2,\hfill \\ & & D_{110}^{\prime }=D_{101}^{\prime }=D_{011}^{\prime }=0,\hfill \\ & & D_{200}^{\prime }={\displaystyle \frac{{\Pi }_{200}}{\rho }},D_{020}^{\prime }={\displaystyle \frac{{\Pi }_{020}}{\rho }},D_{002}^{\prime }={\displaystyle \frac{{\Pi }_{002}}{\rho }},\mathrm{for}(m+n+p)>2,\hfill \\ & & D_{110}^{\prime }={\displaystyle \frac{{\Pi }_{110}}{\rho }},D_{101}^{\prime }={\displaystyle \frac{{\Pi }_{101}}{\rho }},D_{011}^{\prime }={\displaystyle \frac{{\Pi }_{011}}{\rho }}.\hfill \end{array}$$

(55)

3.2. Selected Continuous Markovian Central Moment Attractors in 3D

Based on the above selection of the diffusion tensor parameter, in anticipation of constructing the FPC-LBM in 3D later in this section we now evaluate the continuous Markovian central moment attractor given in Equation (54) for the components independently supported by the D3Q27 lattice. Up to the third moments, which directly appear in the hydrodynamics in the Chapman-Enskog analysis, they read as

$$\begin{array}{ccc}\hfill & & {\Pi }_{000}^{Mv}=\rho ,{\Pi }_{100}^{Mv}={\Pi }_{010}^{Mv}={\Pi }_{001}^{Mv}=0,\hfill \end{array}$$

(56)

$$\begin{array}{ccc}& & {\Pi }_{110}^{Mv}={\Pi }_{101}^{Mv}={\Pi }_{011}^{Mv}=0,{\Pi }_{200}^{Mv}={\Pi }_{020}^{Mv}={\Pi }_{002}^{Mv}=c_s^2\rho ,\hfill \end{array}$$

$$\begin{array}{ccc}& & {\Pi }_{120}^{Mv}={\Pi }_{102}^{Mv}={\Pi }_{210}^{Mv}={\Pi }_{012}^{Mv}={\Pi }_{201}^{Mv}={\Pi }_{021}^{Mv}={\Pi }_{111}^{Mv}=0.\hfill \end{array}$$

In addition, the higher order continuous central moment attractors from the fourth to the sixth order, which complete all the 26 independent moments that can be evolved when using the D3Q27 lattice, are given by

$$\begin{array}{c}\hfill {\Pi }_{220}^{Mv}={\displaystyle \frac{1}{\rho }}\left[{\Pi }_{200}{\Pi }_{020}+2{\Pi }_{110}^2\right],{\Pi }_{202}^{Mv}={\displaystyle \frac{1}{\rho }}\left[{\Pi }_{200}{\Pi }_{002}+2{\Pi }_{101}^2\right],\end{array}$$

$$\begin{array}{c}\hfill {\Pi }_{022}^{Mv}={\displaystyle \frac{1}{\rho }}\left[{\Pi }_{020}{\Pi }_{002}+2{\Pi }_{011}^2\right],{\Pi }_{211}^{Mv}={\displaystyle \frac{1}{\rho }}\left[{\Pi }_{200}{\Pi }_{011}+2{\Pi }_{110}{\Pi }_{101}\right],\end{array}$$

$$\begin{array}{c}\hfill {\Pi }_{121}^{Mv}={\displaystyle \frac{1}{\rho }}\left[{\Pi }_{020}{\Pi }_{101}+2{\Pi }_{110}{\Pi }_{011}\right],{\Pi }_{112}^{Mv}={\displaystyle \frac{1}{\rho }}\left[{\Pi }_{002}{\Pi }_{110}+2{\Pi }_{101}{\Pi }_{011}\right],\end{array}$$

$$\begin{array}{c}\hfill {\Pi }_{122}^{Mv}={\displaystyle \frac{2}{5\rho }}\left[{\Pi }_{020}{\Pi }_{102}+{\Pi }_{002}{\Pi }_{120}+4{\Pi }_{011}{\Pi }_{111}+2({\Pi }_{101}{\Pi }_{021}+{\Pi }_{110}{\Pi }_{012})\right],\end{array}$$

$$\begin{array}{c}\hfill {\Pi }_{212}^{Mv}={\displaystyle \frac{2}{5\rho }}[{\Pi }_{200}{\Pi }_{012}+{\Pi }_{002}{\Pi }_{210}+4{\Pi }_{101}{\Pi }_{111}+2({\Pi }_{110}{\Pi }_{102}+{\Pi }_{011}{\Pi }_{201})],\end{array}$$

$$\begin{array}{c}\hfill {\Pi }_{221}^{Mv}={\displaystyle \frac{2}{5\rho }}[{\Pi }_{200}{\Pi }_{021}+{\Pi }_{020}{\Pi }_{201}+4{\Pi }_{110}{\Pi }_{111}+2({\Pi }_{011}{\Pi }_{210}+{\Pi }_{101}{\Pi }_{120})],\end{array}$$

$$\begin{array}{c}\hfill {\Pi }_{222}^{Mv}={\displaystyle \frac{1}{3\rho }}[{\Pi }_{200}{\Pi }_{022}+{\Pi }_{020}{\Pi }_{202}+{\Pi }_{002}{\Pi }_{220}+4({\Pi }_{110}{\Pi }_{112}+{\Pi }_{101}{\Pi }_{121}+{\Pi }_{011}{\Pi }_{211})].\end{array}$$

(57)

It can be seen from the above set of equations that the Markovian central moment attractors of 4th and higher orders involve various combinations of the quadratic products of the lower order central moments, analogous, but differing in detail, to those appearing in the transformations involved between central moments and cumulants before performing collision using the latter in the cumulant LBM [11]. On the other hand, the factorized LBM [10] prescribes products of lower moments for higher moment attractors in a rather ad hoc manner with substantially fewer terms as the order of the moment increases. Thus, in terms of complexity, we may regard the central moment FP collision model to be intermediate between the factorized and cumulant formulations of the LBM while derivable from refining a collision model of the Boltzmann equation that also has various other applications in statistical mechanics and beyond as noted in the introduction.

Another interesting point to note here is that the representation of the collision step as a set of relaxations to quadratically nonlinear combinations of lower order central moments implies that it cannot be written in a matrix form, just like the cumulant LBM, even if the mappings needed prior to and following collision be expressed in that way; by contrast, in other existing LB models, typically a matrix formulation is contrived to represent the effect of collisions, which effectively implies considering relaxations to central moment equilibria that are linear functions of non-conserved lower order moments. As such, there is no a priori physical reason to express the effect of collisions in a matrix form noting that the Boltzmann’s collision integral term itself is quadratically nonlinear in terms of the distribution functions and the present central moment FP formulation broadly reflects such nonlinearity in the attractors.

3.3. Central Moment of Boltzmann’s Acceleration Term Due to Body Force in 3D

To account for the effect of any external body force $\mathit{F}=(F_x,F_y,F_z)$, as in the 2D case discussed earlier, we first introduce the central moment of the rate of change of the distribution function due to the body force of order ($m+n+p$) as

$$\begin{array}{c}\hfill {\Gamma }_{mnp}=〈{\left({\displaystyle \frac{\delta f}{\delta t}}\right)}_{forcing},W_{mnp}〉.\end{array}$$

(58)

Evaluating this last equation using Equation (3) along with the identities for the central moments of the derivatives of f in the velocity space given in Equation (48), we obtain the rate of change of the central moment of order ($m+n+p$) due to the body force given by

$$\begin{array}{c}\hfill {\Gamma }_{mnp}=m{\displaystyle \frac{F_x}{\rho }}{\Pi }_{m-1,n,p}+n{\displaystyle \frac{F_y}{\rho }}{\Pi }_{m,n-1,p}+p{\displaystyle \frac{F_z}{\rho }}{\Pi }_{m,n,p-1},\end{array}$$

(59)

which generalizes the expression derived earlier for the 2D case to 3D [64].

For completeness, consistency of the above central moment-based FP collision model to the 3D Navier-Stokes equations is shown using a Chapman-Enskog analysis in Appendix C.

3.4. Construction of 3D FPC-LBM

We will now construct a FPC-LBM in 3D on the D3Q27 lattice, with particle velocities ${\mathbf{e}}_{\alpha }$, where $\alpha =1,2,\cdots ,26$, whose cartesian components are given by

$$\begin{array}{ccc}\hfill |{\mathbf{e}}_x\rangle & =& {(0,1,-1,0,0,0,0,1,-1,1,-1,1,-1,1,-1,0,0,0,0,1,-1,1,-1,1,-1,1,-1)}^{†},\hfill \\ \hfill |{\mathbf{e}}_y\rangle & =& {(0,0,0,1,-1,0,0,1,1,-1,-1,0,0,0,0,1,-1,1,-1,1,1,-1,-1,1,1,-1,-1)}^{†},\hfill \\ \hfill |{\mathbf{e}}_z\rangle & =& {(0,0,0,0,0,1,-1,0,0,0,0,1,1,-1,-1,1,1,-1,-1,1,1,1,1,-1,-1,-1,-1)}^{†}.\hfill \end{array}$$

Following the discretization process as outlined in the 2D case, we arrive at a LB scheme with the standard collide-and-stream steps generically represented in Equations (33a) and (33b), where $q=27$. To express the collision step in terms of changes to various central moments via their relaxations to attractors as prescribed by the FP collision model along with those due to the body force for the D3Q27 lattice, we first define the discrete central moments and discrete raw moments of order ($m+n+p$) of the distribution function $f_{\alpha }$, its corresponding Markovian attractor $f_{\alpha }^{Mv}$, and the source term $S_{\alpha }$ for the body force, respectively, as follows:

$$\begin{array}{lll}\left(\begin{array}{c}{\kappa }_{mnp}\\ {\kappa }_{mnp}^{Mv}\\ {\sigma }_{mnp}\end{array}\right)& =& \sum _{\alpha =0}^{26}\left(\begin{array}{c}{f}_{\alpha }\\ f_{\alpha }^{Mv}\\ S_{\alpha }\end{array}\right){(e_{\alpha x}-u_x)}^m{(e_{\alpha y}-u_y)}^n{(e_{\alpha z}-u_z)}^p,\mathrm{and}\end{array}$$

(60a)

$$\begin{array}{lll}\left(\begin{array}{c}{\kappa }_{mnp}^{\prime }\\ {\kappa }_{mnp}^{M{v}^{\prime }}\\ {\sigma }_{mnp}^{\prime }\end{array}\right)& =& \sum _{\alpha =0}^{26}\left(\begin{array}{c}{f}_{\alpha }\\ f_{\alpha }^{Mv}\\ S_{\alpha }\end{array}\right)e_{\alpha x}^m{e}_{\alpha y}^n{e}_{\alpha z}^p.\end{array}$$

(60b)

Then, the countable independent central moments and raw moments for the D3Q27 lattice are, respectively, given by

$$\begin{array}{cc}\hfill {\mathbf{m}}^{\mathbf{c}}=(& {\kappa }_{000},{\kappa }_{100},{\kappa }_{010},{\kappa }_{001},{\kappa }_{110},{\kappa }_{101},{\kappa }_{011},{\kappa }_{200},{\kappa }_{020},{\kappa }_{002},{\kappa }_{120},{\kappa }_{102},{\kappa }_{210},{\kappa }_{012},\hfill \\ & {\kappa }_{201},{\kappa }_{021},{\kappa }_{111},{\kappa }_{220},{\kappa }_{202},{\kappa }_{022},{\kappa }_{211},{\kappa }_{121},{\kappa }_{112},{\kappa }_{122},{\kappa }_{212},{\kappa }_{221},{\kappa }_{222}),\hfill \end{array}$$

(61a)

and

$$\begin{array}{cc}\hfill \mathbf{m}=(& {\kappa }_{000}^{\prime },{\kappa }_{100}^{\prime },{\kappa }_{010}^{\prime },{\kappa }_{001}^{\prime },{\kappa }_{110}^{\prime },{\kappa }_{101}^{\prime },{\kappa }_{011}^{\prime },{\kappa }_{200}^{\prime },{\kappa }_{020}^{\prime },{\kappa }_{002}^{\prime },{\kappa }_{120}^{\prime },{\kappa }_{102}^{\prime },{\kappa }_{210}^{\prime },{\kappa }_{012}^{\prime },\hfill \\ & {\kappa }_{201}^{\prime },{\kappa }_{021}^{\prime },{\kappa }_{111}^{\prime },{\kappa }_{220}^{\prime },{\kappa }_{202}^{\prime },{\kappa }_{022}^{\prime },{\kappa }_{211}^{\prime },{\kappa }_{121}^{\prime },{\kappa }_{112}^{\prime },{\kappa }_{122}^{\prime },{\kappa }_{212}^{\prime },{\kappa }_{221}^{\prime },{\kappa }_{222}^{\prime }),\hfill \end{array}$$

(61b)

and, similarly, one can write those for the attractors and sources. In addition, we can also express a 27-dimensional vector $\mathbf{f}$ containing the distribution functions via

$$\mathbf{f}=(f_0,f_1,f_2,\cdots ,f_{26}).$$

(62)

To express the collision step in terms of relaxations of different central moments to their respective attractors along with source term updates, we have previously introduced the mappings between distribution functions $\mathbf{f}$ to raw moments $\mathbf{m}$ (see Equation (37)) and those between the raw moments $\mathbf{m}$ and central moments ${\mathbf{m}}^{\mathbf{c}}$(see Equation (39)). The projection matrix $\mathsf{P}$ appearing in the former case for the D3Q27 lattice can be written as

$$\begin{array}{c}\hfill \mathsf{P}=\left[\right|\mathbf{1}\rangle ,|{\mathbf{e}}_x\rangle ,|{\mathbf{e}}_y\rangle ,|{\mathbf{e}}_z\rangle ,|{\mathbf{e}}_x{\mathbf{e}}_y\rangle ,|{\mathbf{e}}_x{\mathbf{e}}_z\rangle ,|{\mathbf{e}}_y{\mathbf{e}}_z\rangle ,|{\mathbf{e}}_x^2\rangle ,|{\mathbf{e}}_y^2\rangle ,|{\mathbf{e}}_z^2\rangle ,|{\mathbf{e}}_x{\mathbf{e}}_y^2\rangle ,|{\mathbf{e}}_x{\mathbf{e}}_z^2\rangle ,\\ \hfill |{\mathbf{e}}_x^2{\mathbf{e}}_y\rangle ,|{\mathbf{e}}_y{\mathbf{e}}_z^2\rangle ,|{\mathbf{e}}_x^2{\mathbf{e}}_z\rangle ,|{\mathbf{e}}_y^2{\mathbf{e}}_z\rangle ,|{\mathbf{e}}_x{\mathbf{e}}_y{\mathbf{e}}_z\rangle ,|{\mathbf{e}}_x^2{\mathbf{e}}_y^2\rangle ,|{\mathbf{e}}_x^2{\mathbf{e}}_z^2\rangle ,|{\mathbf{e}}_y^2{\mathbf{e}}_z^2\rangle ,|{\mathbf{e}}_x^2{\mathbf{e}}_y{\mathbf{e}}_z\rangle ,\\ \hfill |{\mathbf{e}}_x{\mathbf{e}}_y^2{\mathbf{e}}_z\rangle ,|{\mathbf{e}}_x{\mathbf{e}}_y{\mathbf{e}}_z^2\rangle ,|{\mathbf{e}}_x{\mathbf{e}}_y^2{\mathbf{e}}_z^2\rangle ,|{\mathbf{e}}_x^2{\mathbf{e}}_y{\mathbf{e}}_z^2\rangle ,|{\mathbf{e}}_x^2{\mathbf{e}}_y^2{\mathbf{e}}_z\rangle ,|{\mathbf{e}}_x^2{\mathbf{e}}_y^2{\mathbf{e}}_z^2\rangle {]}^{†}\end{array}$$

in which $\left|\mathbf{1}\right.〉$ is a 27-dimensional vector with unit elements forming a basis for the zeroth moment given by $\left|\mathbf{1}\right.〉={(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)}^{†}$; the frame transformation matrix $\mathsf{F}$ (as well as its inverse ${\mathsf{F}}^{-1}$) appearing in the latter case depend on the local fluid velocity components $\mathit{u}=(u_x,u_y,u_z)$ and are low triangular matrices in structure that can be readily obtained via extracting the coefficients of the binomial expansions of central moments in terms of raw moments (see e.g., [68]).

The discrete Markovian central moment attractors ${\kappa }_{mnp}^{Mv}$ needed in the construction of the 3D FPC-LBM are obtained, as in the 2D case, by matching the corresponding continuous versions given in Equations (56) and (57), i.e., ${\kappa }_{mnp}^{Mv}={\Pi }_{mnp}^{Mv}$ for all the independent countable moments supported by the D3Q27 lattice. Thus, we get

$${\kappa }_{000}^{Mv}=\rho ,{\kappa }_{100}^{Mv}=0,{\kappa }_{010}^{Mv}=0,{\kappa }_{001}^{Mv}=0,$$

$${\kappa }_{110}^{Mv}=0,{\kappa }_{101}^{Mv}=0,{\kappa }_{011}^{Mv}=0,{\kappa }_{200}^{Mv}=\rho c_s^2,{\kappa }_{020}^{Mv}=\rho c_s^2,{\kappa }_{002}^{Mv}=\rho c_s^2,$$

$${\kappa }_{120}^{Mv}=0,{\kappa }_{102}^{Mv}=0,{\kappa }_{210}^{Mv}=0,{\kappa }_{012}^{Mv}=0,{\kappa }_{201}^{Mv}=0,{\kappa }_{021}^{Mv}=0,{\kappa }_{111}^{Mv}=0,$$

$$\begin{array}{|c|}\hline {\kappa }_{220}^{Mv}={\displaystyle \frac{1}{\rho }}({\tilde{\kappa }}_{200}{\tilde{\kappa }}_{020}+2{\tilde{\kappa }}_{110}{\tilde{\kappa }}_{110})\\ \hline\end{array},\begin{array}{|c|}\hline {\kappa }_{202}^{Mv}={\displaystyle \frac{1}{\rho }}({\tilde{\kappa }}_{200}{\tilde{\kappa }}_{002}+2{\tilde{\kappa }}_{101}{\tilde{\kappa }}_{101})\\ \hline\end{array},$$

$$\begin{array}{|c|}\hline {\kappa }_{022}^{Mv}={\displaystyle \frac{1}{\rho }}({\tilde{\kappa }}_{020}{\tilde{\kappa }}_{002}+2{\tilde{\kappa }}_{011}{\tilde{\kappa }}_{011})\\ \hline\end{array},$$

$$\begin{array}{|c|}\hline {\kappa }_{211}^{Mv}={\displaystyle \frac{1}{\rho }}({\tilde{\kappa }}_{200}{\tilde{\kappa }}_{011}+2{\tilde{\kappa }}_{110}{\tilde{\kappa }}_{101})\\ \hline\end{array},\begin{array}{|c|}\hline {\kappa }_{121}^{Mv}={\displaystyle \frac{1}{\rho }}({\tilde{\kappa }}_{020}{\tilde{\kappa }}_{101}+2{\tilde{\kappa }}_{110}{\tilde{\kappa }}_{011})\\ \hline\end{array},$$

$$\begin{array}{|c|}\hline {\kappa }_{112}^{Mv}={\displaystyle \frac{1}{\rho }}({\tilde{\kappa }}_{002}{\tilde{\kappa }}_{110}+2{\tilde{\kappa }}_{011}{\tilde{\kappa }}_{101})\\ \hline\end{array},$$

$$\begin{array}{|c|}\hline {\kappa }_{122}^{Mv}={\displaystyle \frac{2}{5\rho }}({\tilde{\kappa }}_{020}{\tilde{\kappa }}_{102}+{\tilde{\kappa }}_{002}{\tilde{\kappa }}_{120}+4{\tilde{\kappa }}_{011}{\tilde{\kappa }}_{111}+2({\tilde{\kappa }}_{101}{\tilde{\kappa }}_{021}+{\tilde{\kappa }}_{011}{\tilde{\kappa }}_{012}))\\ \hline\end{array},$$

$$\begin{array}{|c|}\hline {\kappa }_{212}^{Mv}={\displaystyle \frac{2}{5\rho }}({\tilde{\kappa }}_{200}{\tilde{\kappa }}_{012}+{\tilde{\kappa }}_{002}{\tilde{\kappa }}_{210}+4{\tilde{\kappa }}_{101}{\tilde{\kappa }}_{111}+2({\tilde{\kappa }}_{110}{\tilde{\kappa }}_{102}+{\tilde{\kappa }}_{011}{\tilde{\kappa }}_{201}))\\ \hline\end{array},$$

$$\begin{array}{|c|}\hline {\kappa }_{221}^{Mv}={\displaystyle \frac{2}{5\rho }}({\tilde{\kappa }}_{200}{\tilde{\kappa }}_{021}+{\tilde{\kappa }}_{020}{\tilde{\kappa }}_{201}+4{\tilde{\kappa }}_{110}{\tilde{\kappa }}_{111}+2({\tilde{\kappa }}_{011}{\tilde{\kappa }}_{210}+{\tilde{\kappa }}_{101}{\tilde{\kappa }}_{120}))\\ \hline\end{array},$$

$$\begin{array}{|c|}\hline {\kappa }_{222}^{Mv}={\displaystyle \frac{1}{3\rho }}({\tilde{\kappa }}_{200}{\tilde{\kappa }}_{022}+{\tilde{\kappa }}_{020}{\tilde{\kappa }}_{202}+{\tilde{\kappa }}_{002}{\tilde{\kappa }}_{220}+4({\tilde{\kappa }}_{110}{\tilde{\kappa }}_{112}+{\tilde{\kappa }}_{101}{\tilde{\kappa }}_{121}+{\tilde{\kappa }}_{011}{\tilde{\kappa }}_{211}))\\ \hline\end{array}.$$

Thus, the above form a chain of attractors, especially for the fourth and higher central moments, each of which depends on certain combinations of the most recent (current) state of the adjacent and relevant lower moment, differing by a degree of two. In practical implementations, such lower central moments needed in the attractors are already in the post-collision states and the updated values will be used to prescribe the attractors for the higher central moments when they relax under collision. In addition, the source terms to the various central moments expressing the contributions from the presence of any body force is expressed from matching with their continuous counterparts given in Equation (59), which can be compactly written as

$$\begin{array}{c}\hfill {\sigma }_{mnp}=m{\displaystyle \frac{F_x}{\rho }}{\kappa }_{m-1,n,p}+n{\displaystyle \frac{F_y}{\rho }}{\kappa }_{m,n-1,p}+p{\displaystyle \frac{F_z}{\rho }}{\kappa }_{m,n,p-1}.\end{array}$$

(63)

The implementation details of the 3D FPC-LBM using the D3Q27 lattice are given in Appendix D.

4. Results and Discussion

We now demonstrate numerically the accuracy and stability properties of the new LB schemes based on central moment Fokker-Planck guided collision model, viz., the FPC-LBMs developed in the previous sections through a series of simulations of benchmark flow problems in two and three dimensions respectively. In two dimensions using the D2Q9 lattice, we show the accuracy of the FPC-LBM using a common test case, viz., the lid-driven square cavity flow, by comparing simulation results with those found in literature. Next, the stability of the FPC-LBM is demonstrated in two dimensions with the double periodic shear layer simulations, which give a visual representation of stability. Specifically, in this test case, two shear layers with an initial perturbation begin to roll up on themselves. When the grid resolution is relatively too coarse, secondary vortices begin to form at the thinnest parts of each layer for some of the LBM formulations, especially the SRT and MRT formulations. On the other hand, the central moment based LBM collision models are found to more robustly deal with such instabilities and this is shown visually when the secondary vortices mentioned above do not form for the same simulation parameters as those used with a less stable method. Furthermore, we show a direct comparison between two central moment collision models, the FPC-LBM and the MCM-LBM, which uses Maxwellian central moments for equilibria, by simulating this flow problem at relatively high Mach numbers as well as beyond the usual limits of the relaxation parameter for the bulk viscosity when the simulation is generally accurate, and demonstrate the superior numerical performance of the former compared to the latter.

Then, in three dimensions, the lid-driven cubic cavity flow benchmark is used to indicate both accuracy and stability characteristics of FPC-LBM using the D3Q27 lattice. The orthogonal crossing shear wave test is then used to show how the FPC-LBM can deal with numerical hyperviscosity effects that arise in simulations of flows with extremely small physical fluid viscosities. For these cases, when appropriate, we show comparisons of the FPC-LBM with various other existing collision models in LBM. These include the single relaxation time (SRT)-LBM, the non-orthogonal raw moments-based multiple relaxation time (MRT)-LBM, the cascaded collision model with non-orthogonal moment basis using Maxwellian central moments for equilibria referred to as the MCM-LBM, the factorized LBM, and the cumulant LBM. Lastly, we show the capability of FPC-LBM in accurately and robustly resolving complex flows by simulating the canonical turbulent channel flow and comparing the resulting statistics against a recent direct numerical simulations (DNS) data available in the literature.

4.1. Two-Dimensional Lid-Driven Square Cavity Flow: An Accuracy Study

The benchmark known as lid-driven cavity flow consists of a square cavity with a viscous fluid inside; the cavity has a lid moving in the horizontal or x direction with constant velocity U. In turn, the moving lid acts to shear the fluid which then forms vortical structures due to interaction with the surrounding stationary walls of the cavity. Furthermore, as the Reynolds number is increased, the characteristic velocity of the flow is then increased which causes counter-rotating vortices of varying sizes to form in the corners of the cavity. Here, we perform an accuracy study by comparing the results obtained using the FPC-LBM with those of Ghia et al. (1982) [69] for the cases of Reynolds numbers of $Re=1000,3200,5000$, and 7500. Here, the Reynolds number is defined as $Re=U_o{L}_o/\nu$, where $U_o$ is the characteristic plate velocity or U, $L_o$ is the characteristic length of one side of the square cavity, and $\nu$ is the kinematic viscosity of the fluid. The relaxation parameters related to the kinematic viscosity of the fluid, which depends on the choice of $Re$, are set according to Equation (A22) (or Equation (A54) in the 3D case), and, unless otherwise stated, that related to the bulk viscosity $\zeta$ and the others for relaxing third and higher moments are set to unity; the speed of sound $c_s$ is taken to be $1/\sqrt{3}$ in all the simulations reported in this work. We used a grid resolution of $500\times 500$ and performed all our simulations at various Reynolds numbers noted above until they reach steady state by ensuring that the 2-norm of the residual error is less than $1\times {10}^{-15}$.

Shown in Figure 3 are the streamlines that come from the results using the FPC-LBM, and are found to be consistent with those found in literature [69]. However, to make a direct comparison, we find the location of the ($x,y$) coordinates of the center of each of the corner vortices and compare them with the reference data reported by Ghia et al. (1982) [69] for all the cases of $Re$ noted above. These results are tabulated in

Table 1. Comparison of the location of the center of various vortices in a 2D lid-driven square cavity flow for different Reynolds numbers. Computed results are obtained using the FPC-LBM with a grid resolution of $500\times 500$ and are compared with those of Ghia et al. (1982) [69]. PV—Primary Vortex, BR1—Bottom Right 1, BR2—Bottom Right 2, BL—Bottom Lefft, TL—Top Left.

Vortex	Model	$\mathbf{Re}=1000$	$\mathbf{Re}=3200$	$\mathbf{Re}=5000$	$\mathbf{Re}=7500$
PV	FPC-LBM	$(0.5306,0.5650)$	$(0.5185,0.5395)$	$(0.5158,0.5344)$	$(0.5138,0.5318)$
PV	Ghia et al.	$(0.5313,0.5625)$	$(0.5165,0.5469)$	$(0.5117,0.5352)$	$(0.5117,0.5322)$
BR1	FPC-LBM	$(0.8646,0.1115)$	$(0.8249,0.0834)$	$(0.8054,0.0720)$	$(0.7913,0.0643)$
BR1	Ghia et al.	$(0.8594,0.1094)$	$(0.8125,0.0859)$	$(0.8086,0.0742)$	$(0.7813,0.0625)$
BR2	FPC-LBM	$(0.9929,0.0057)$	$(0.9903,0.0083)$	$(0.9795,0.0172)$	$(0.9533,0.0401)$
BR2	Ghia et al.	$(0.9922,0.0078)$	$(0.9844,0.0078)$	$(0.9805,0.0195)$	$(0.9492,0.0430)$
BL	FPC-LBM	$(0.0830,0.0771)$	$(0.0810,0.1191)$	$(0.0736,0.1344)$	$(0.0649,0.1510)$
BL	Ghia et al.	$(0.0859,0.0781)$	$(0.0859,0.1094)$	$(0.0703,0.1367)$	$(0.0645,0.1504)$
TL	FPC-LBM	$N/A$	$(0.0574,0.8974)$	$(0.0635,0.9089)$	$(0.0669,0.9102)$
TL	Ghia et al.	$N/A$	$(0.0547,0.8984)$	$(0.0625,0.9102)$	$(0.0664,0.9141)$

, side by side with our comparison case, and are found to be in relatively very good quantitative agreement. Moreover, the centerline horizontal and vertical components of the velocity profiles computed using the FPC-LBM for all the above four representative $Re$ are in excellent agreement with the reference results [69] (see Figure 4 and Figure 5).

4.2. Doubly Periodic Shear Layers: Numerical Performance Study

The benchmark commonly known as doubly periodic shear layers [70] is employed here to indicate a level of stability of the FPC-LBM as compared to other collision models such as MCM-LBM, which uses the Maxwellian based equilibria, i.e., ${\kappa }_{mn}^{eq}={\Pi }_{mn}^{\left(0\right)}=〈f^M,W_{mn}〉$, where $f^M$ is given in Equation (6). Since the SRT-LBM and MRT-LBM are already known to be less stable compared to central moment-based LB formulations, our focus will be on comparing the different possible central moment-based LB scheme for simulating this flow problem. The key idea is that a pair of shear layers in a periodic box will roll up on themselves if given an initial perturbation in the direction perpendicular to the shear layers. Then, looking at the resulting vorticity field at a dimensionless time equal to one, we observe if the flow field is fully resolved. When using a relatively fine enough grid resolution, the layers will each be shown to roll up in a single vortex at the location of the initial perturbation for any LBM collision operator being used, which is considered to be a fully resolved flow field. However, if the discretized grid is made to be relatively coarse, then spurious secondary vortices can form at the thinnest parts of each layer due to the presence of numerical artifacts associated with the numerical method being used, which could ultimately destabilize the simulation. Furthermore, the Mach number is used here to enforce the characteristic speed of the shear layers, which can trigger numerical instabilities when made progressively larger for a chosen grid resolution, and one of the goals is to determine which one of the central moment LB formulation sustains stable simulations relative to the other in such cases.

The following equation is used to initialize the flow field $\mathit{u}(t=0)=(u_x,u_y)$ by non-dimensionalizing the spatial coordinates with the side length of the periodic square domain, $L_o$:

$$\begin{array}{ccc}& u_x& =\left\{\begin{array}{c}{U}_{o}tanh(4(y-0.25)/w)y\le 0.5,\hfill \\ U_{o}tanh(4(0.75-y)/w)y>0.5,\hfill \end{array}\right.\hfill \\ & u_y& =\delta sin\left(2\pi \right(x+0.25\left)\right),\hfill \end{array}$$

(64)

where $\delta$ is related to the magnitude of the initial perturbation and w is related to the thickness of the shear layers; here, we use $\delta =0.05$ and $w=0.05$. In all of the following simulations, we set the Reynolds number to $Re=$ 30,000, where the Reynolds number for this problem is defined by $Re=U_o{L}_o/\nu$. The Mach number associated in the simulations is defined as $Ma=U_o/c_s$, where $c_s=1/\sqrt{3}$, and the characteristic time scale reads as $T=L_o/U_o$. Here, we compare the different collision models by running simulations over a range of grid resolutions as well as Mach numbers, and compare the resulting flow fields at the dimensionless time t, obtained by scaling time using T, of 1.0.

In what follows, we make a direct comparison of the resulting vorticity fields between the FPC-LBM and the MCM-LBM. The results shown in Figure 6 and Figure 7 consist of baseline cases using the MCM-LBM and FPC-LBM, respectively, using the Mach numbers of $Ma=0.05,0.2$, and $0.3$ along with grid resolutions of $L^2={64}^2,{128}^2$, and ${256}^2$. As such, the Mach numbers are generally in the lower range and hence all of the cases were found to be numerically stable without the formation of any spurious secondary vortices, and the results from both the collision models are almost indistinguishable. This generally shows the robustness of central moment collision models and also indicates that we must consider more extreme parametric conditions to see significant differences between them.

In this regard, we perform two more comparison studies that attempt this by first using larger Mach numbers and second by using a relatively large increase in bulk viscosity that exceeds the point of having positive effects. The results from the larger Mach number cases of $Ma=0.4,0.5$, and $0.6$, and using the same three different grid resolutions as before are presented in Figure 8 using MCM-LBM and Figure 9 using FPC-LBM. Evidently, from Figure 8, the MCM-LBM simulations show the formation of spurious secondary vortices for all cases of these higher $Ma$ for the coarsest grid resolution of ${64}^2$, and these artifacts progressively become more prominent as $Ma$ increases. On the other hand, the results from the FPC-LBM (see Figure 9) show no such spurious secondary vortices for those same cases, which indicates that the FPC-LBM can sustain stable simulations without any noticeable artifacts even at large Mach numbers and coarse grid resolutions.

Lastly, we consider cases in the lower Mach number ranges again ($Ma=0.05$, $0.2$ and $0.3$) and with the same sets of grid resolutions as before, but at an extremely large bulk viscosity by setting the relaxation parameter associated with the bulk viscosity to ${\omega }_3=0.35$. Figure 10 shows the vorticity contours at $t=1$ computed using MCM-LBM. In general, increasing the bulk viscosity dampens any spurious pressure waves; however, by increasing the relaxation time ${\tau }_3=1/{\omega }_3$ associated with the bulk viscosity to $2.857$ from being ∼1.0 or smaller, the range of beneficial limits is exceeded, and instead the simulations begin to numerically destabilize with MCM-LBM with the incipient spurious secondary vortices at the coarsest resolution even with the usage of lower range of Mach numbers (see Figure 10); by contrast, as shown in Figure 11 using the FPC-LBM does not result in spurious vortices for these cases and the simulations remain stable. This again indicates that the FPC-LBM is a more robust central moment LB scheme than the MCM-LBM.

4.3. Three-Dimensional Lid-Driven Cubic Cavity Flow: An Accuracy Study

Next, let’s validate the 3D FPC-LBM by performing simulations of flow inside a cubic cavity driven by the motion of one of its lids and compare the results with the benchmark numerical data available in the literature [71,72]. This problem is set up with a cubic grid domain of $L_x\times L_y\times L_z=150\times 150\times 150$, with a lid, located at $y=L_y$, that moves with a constant velocity, $U_o$, in the x direction. The characteristic plate velocity is set by specifying the Mach number to be relatively small, with $Ma=U_o/c_s=0.1$, so that the flow is maintained in the weakly compressible range; the choice of the Reynolds number, defined by $Re=U_o{L}_o/\nu$ with $L_o=L_x=L_y=L_z$, sets up the fluid viscosity $\nu$. All of the other walls are stationary, and the standard half-way bounce back boundary condition is applied to all boundaries in order to enforce the no slip condition on each wall, including a momentum correction to the moving wall.

Figure 11. Vorticity contours of doubly periodic shear layers that roll up due to an applied perturbation at $t=1$ for different sets of lower Mach numbers of $0.05$, $0.2$ and $0.3$ (along rows) and at grid resolutions of ${64}^2$, ${128}^2$, and ${256}^2$ (along columns) at an extremely large bulk viscosity by setting the relaxation parameter associated with the bulk viscosity to ${\omega }_3=0.35$ and computed using the Fokker-Planck equilibria based FPC-LBM. It is seen that the FPC-LBM remains stable even for the cases with an extreme increase in bulk viscosity shown previously, where the MCM-LBM did not remain stable (see Figure 10). This indicates that the FPC-LBM is numerically more stable when compared to the MCM-LBM in such cases as well.

Figure 11. Vorticity contours of doubly periodic shear layers that roll up due to an applied perturbation at $t=1$ for different sets of lower Mach numbers of $0.05$, $0.2$ and $0.3$ (along rows) and at grid resolutions of ${64}^2$, ${128}^2$, and ${256}^2$ (along columns) at an extremely large bulk viscosity by setting the relaxation parameter associated with the bulk viscosity to ${\omega }_3=0.35$ and computed using the Fokker-Planck equilibria based FPC-LBM. It is seen that the FPC-LBM remains stable even for the cases with an extreme increase in bulk viscosity shown previously, where the MCM-LBM did not remain stable (see Figure 10). This indicates that the FPC-LBM is numerically more stable when compared to the MCM-LBM in such cases as well.

Then, we simulate the three-dimensional lid-driven cavity flow for Reynolds numbers of $Re=100,400,1000$. Figure 12 shows the computed streamlines at these Reynolds numbers along three different centerplanes. As $Re$ increases, the flow patterns become progressively more complex with additional vortical structures near the corner edges, and the observed patterns are consistent with those reported in [71,72]. In addition, comparisons of the horizontal and vertical velocity profiles computed using the FPC-LBM along different directions in the centerplanes for the above three $Re$ with the prior numerical results obtained from Navier-Stokes solvers [71,72] are presented in Figure 13. Very good quantitative agreement is seen for all the cases considered. These results give an indication that the 3D FPC-LBM accurately simulates flow fields with rather complex structures.

One of the major advantages of the LBM is its natural parallelization property. We have parallelized our implementation of the 3D FPC-LBM using the Message Passing Interface (MPI) library via a domain decomposition approach. We then performed simulations of the cubic cavity flow resolved with $150\times 150\times 150$ using $p=16,32,64,128,256$, 512 processors in our in-house computer cluster (see https://ccm-docs.readthedocs.io/en/latest/alderaan (accessed on 10 October 2024)). Figure 14 shows the parallel performance of our 3D FPC-LBM. Evidently, near-linear speed-up is seen, which is ideal for simulating large domain size problems effectively. We will exploit this feature for turbulent channel flow simulations discussed later.

4.4. Three-Dimensional Lid-Driven Cubic Cavity Flow: A Stability Study

In order to further illustrate the numerical properties of the 3D FPC-LBM in perspective with other existing LB collision models, next we investigate their stability characteristics in the simulation of three-dimensional, lid-driven cubic cavity flow. If ${\omega }_j$, where $j=110,101,011,2d1$, and $2d2$ represent the relaxation parameters of the second order moments that determine the fluid viscosity $\nu$ (see Equation (A54)), for brevity in discussion, we let them to be related to a relaxation time $\tau$ via $\tau =1/{\omega }_j$. It is well known that as $\tau$ approaches $1/2$ or equivalently, as the viscosity becomes relatively small, the LB schemes are susceptible to numerical instabilities that can cause any errors in simulations to grow exponentially large, and different collision models can exhibit such behavior to a different degree depending on when a certain threshold viscosity or the Reynolds number is reached before they become unstable. Here, we leverage this feature to form a description for the numerical stability properties of the FPC-LBM as compared to other collision models.

We use the three-dimensional lid driven cavity flow simulation as a prototypical example in this regard (see e.g., Ref. [73] for a recent such simulation study using the cascaded LBM and other collision models). In this regard, we start with the case that we know is numerically unstable, namely $\tau =0.5$, which has a corresponding Reynolds number, $Re=\infty$, when the characteristic velocity and length scales are fixed, as the relaxation time $\tau$ is related to the fluid viscosity, $\nu$. We then systematically increase the relaxation time until the simulation can remain stable, and report the resulting threshold Reynolds number associated with that relaxation time as being the maximum possible Reynolds number for which the simulation is stable. Specifically, we increase the relaxation time as ${\tau }_{new}={\tau }_{old}+\delta$, where $\delta$ is the amount of increase for each case after the previous simulation has become unstable. Here, we use $\delta =5\times {10}^{-5}$, and we consider any simulation that reaches 500,000 time steps to be stable. Furthermore, after we find the minimum value for the relaxation time that allows the simulation to remain stable for 500,000 time steps, we then validate that every value above this minimum also remains stable for at least ten more increments in it. At this point, we define the simulation as having crossed the threshold from unstable to stable, and we note the Reynolds number at the edge of that threshold. We perform this study using four different grid resolutions and compare the results from five different collision models. The grid resolutions we consider for the cubic cavity of dimension $L_x\times L_y\times L_z=L^3$, are $L=48,64,80$, and 96 and the plate velocity $U_0$ is based on the Mach number via $U_o=Ma\times c_s$ by fixing $Ma=0.2$ in all cases. We then compare the results of these cases across the MRT-LBM, MCM-LBM, factorized LBM, cumulant LBM, and the FPC-LBM approaches, but do not consider the SRT-LBM as it is already known to be the least stable approach in LBM in general. As a result, this study gives an indication of the numerical stability features of the FPC-LBM as compared to other collision models. Figure 15 shows the maximum possible Reynolds number achieved using different collision models at various resolutions. In general, the central moment-based schemes (MCM-LBM, factorized LBM, and FPC-LBM) are more stable compared to the raw moment-based MRT-LBM. Among the central moment-based schemes, MCM-LBM is the least stable and the FPC-LBM consistently outperforms it by a factor greater than 2 for the maximum threshold $Re$, and the latter is more stable than the ad hoc factorized LBM as well. As such, we find that the FPC-LBM has stability characteristics that are very close to those observed for the cumulant LBM (see Figure 15).

A note regarding the computational cost of various LB schemes is in order here. In general, depending on the computing platform used, the stream-and-collide procedure of the LB methods are known to be memory-bound rather than compute-bound and hence the computational efforts associated with different collision models do not vary significantly. On a per-iteration basis, for example, the raw moment-based MRT-LBM requires about 20–30% more computational time than the simplest of the collision models, viz., the SRT-LBM. Our FPC-LBM, which is based on central moments and using Markovian attractors, requires similar computational effort as the Maxwellian-based central moment (MCM)-LBM and involves fewer operations than the state-of-the-art cumulant LBM. As such, while our FPC-LBM requires about 50% additional overhead when compared to the SRT-LBM, it is offset by the several numerical advantages associated with the Fokker-Planck formulation using central moments as discussed in this work, including its dramatically improved numerical stability for simulations of fluid flows at lower viscosities or higher Reynolds numbers much more effectively.

4.5. Orthogonal Crossing Shear Waves: A Numerical Hyperviscosity Study

This benchmark is applied in effort to asses the presence of a numerical artifact known as hyperviscosity that has been observed in some of the LBM formulations, especially those that relax different moments at different rates during collision [10,11]. More specifically, when computing flows with relatively very small fluid viscosities, there can exist large disparities between the relaxation rates of the second order moments compared to the higher order moments. In turn, depending on the choice of equilibria in such cases, the contributions from such higher order moments involve terms similar to the non-equilibrium momentum fluxes (related to the strain rate tensor) which can dominate the corresponding physical contributions from the second order non-equilibrium moments, and manifest as numerical hyperviscosities [11]. Such a numerical artifact can be investigated by studying the decay rate of three-dimensional crossing shear waves [10].

If we consider the fluid velocities and the pressure gradients to be small, and in the absence of body forces or boundaries, for a pair of three-dimensional crossing shear waves involving orthogonal wave vectors, the Navier-Stokes equations can be simplified to obtain the following analytical solution based on the product solution of the respective shear waves:

$$\begin{array}{c}\hfill u_y\left(t\right)=U_{y}cos\left(k_{x}x\right)cos\left(k_{z}z\right)exp(-\nu k_x^2)exp(-\nu k_z^2),\end{array}$$

(65)

where $k_x=2\pi /L_x$ and $k_z=2\pi /L_z$ are the wavenumbers in the two orthogonal coordinate directions x and z, respectively, $U_y$ is the initial amplitude of the wave perturbed along the y direction, and $\nu$ is the fluid viscosity. Thus, starting from the initial condition in a periodic box

$$\begin{array}{c}\hfill u_x=0,u_z=0,u_y=U_{y}cos\left(k_{x}x\right)cos\left(k_{z}z\right),\end{array}$$

(66)

according to Equation (65), the crossing shear waves decay at a rate that depends on the wavenumbers and the viscosity of the fluid. Furthermore, given that the decay rate depends on the fluid viscosity, the presence of any numerical hyperviscosities can significantly alter the amplitude decay rate from the predictions above. Here, we discretize such a periodic box with a relatively coarse grid resolution and apply the initial conditions above, so that we can study different LBM formulations and investigate their ability to predict the decay rate.

As in [10], we consider an extreme case of relatively small fluid viscosity of $\nu =1\times {10}^{-7}$, with wavenumbers defined using $L_x=L_z=30$ and an initial amplitude of $U_y=1\times {10}^{-5}$, and we study the decay rate for simulations using various LBM collision formulations relative to the analytically predicted decay rate shown above. Namely, we compare the decay rates as produced by the single relaxation time collision model or the SRT-LBM, the central moment collision model which uses the Maxwellian central moments for equilibria or the MCM-LBM, the central moment collision model with the Fokker-Planck guided collisions or the FPC-LBM, and the cumulant LBM involving relaxations based on cumulants during collision.

The results in Figure 16a indicate that the decay rate for the MCM-LBM is much higher than the analytically derived decay rate because the velocity amplitude decays rapidly to zero, and thus is not able to effectively overcome the hyperviscosity artifact. On the other hand the SRT-LBM, FPC-LBM, and the cumulant LBM have apparently similar decay rates as compared to what is predicted analytically. Figure 16b shows a closer view of these results, which more clearly illuminates the key features of the other three collision models shown here. To begin with, the SRT-LBM produces a decay rate for the crossing shear waves that is consistent with what is expected; however it has other features which cause noise in the decay rate, and this is consistent with observations shown in [10]. On the other hand, we observed that the FPC-LBM and cumulant LBM both produce nearly identical decay rates in this benchmark, which are in agreement with those predicted by the analytical solution. As such, it underscores the importance of equilibria in reducing the hyperviscosity effects dramatically. Specifically, with the FPC-LBM, the choice of the central moment equilibria based on the Markovian attractor given in Equation (54) and its attendant tensor diffusion parameter as presented in Equation (55) (which controls the rate of diffusion of the distribution function in different directions in the velocity space) play a crucial role in this regard.

4.6. Fully Developed Turbulent Channel Flow: Demonstration Case Study of FPC-LBM for Turbulence Simulations

In our final case study, we demonstrate the capability of the FPC-LBM to perform turbulence simulations using a prototypical example involving the turbulent channel flow and compare the computed turbulence statistics against a recent DNS data of Lee and Moser (2015) [74]. To reduce the computational effort, we perform simulations within a domain that encompasses only the half-channel height H by utilizing a specular refection boundary condition on the symmetry plane [75]. For resolving the turbulent flow structures adequately, the computational domain is taken to have an aspect ratio of $2\pi H\times \pi H\times H$. The flow is driven by a body force along the periodic streamwise direction given by $F_x=-{\displaystyle \frac{dP}{dx}}={\displaystyle \frac{{\tau }_w}{H}}={\displaystyle \frac{\rho u_{*}^2}{H}}$, where ${\tau }_w$ is the wall shear stress and $u_{*}$ is the shear or friction velocity, which satisfies $u_{*}={({\tau }_w/\rho )}^{1/2}$. For performing turbulence simulations, we consider a shear Reynolds number ${Re}_{*}=180$, where ${Re}_{*}=u_{*}H/\nu$ with $\nu$ being the molecular viscosity of the fluid, by resolving the computational domain using the above aspect ratio by taking $H=100$ grid nodes in the wall normal direction. We choose a similar approach as that of [75], where we perform large eddy simulations (LES) by utilizing a common subgrid scale model known as the Smagorinsky model along with the van Driest wall dampening function to represent the variations in the subgrid effects near the wall. The characteristic length scale near the wall is the viscous length scale ${\delta }_{\nu }=\nu /u_{*}$, which can be used to express the characteristic time scale of eddy as $T^{*}=H/u_{*}$. The initial run for the simulations are performed using the 3D FPC-LBM for a duration of $50T^{*}$ until stationary turbulence statistics (such as the invariance in the Reynolds stress profiles) is achieved; then an additional run for a period of $30T^{*}$ is carried out to collect the turbulence statistics by averaging the flow field and turbulence fluctuations in time and in space along the homogeneous directions (i.e., the horizontal planes). Large eddy simulations performed using parallel computations with 512 processors of our in-house cluster took less than 24 h of wall-clock time to complete and to collect the turbulence statistics. Figure 17 shows comparisons of the mean streamwise velocity, root-mean-square (rms) velocity fluctuations, and the Reynolds stress as a function of the distance from the wall in wall units, i.e., $z^{+}=z/{\delta }_{\nu }$ computed using the FPC-LBM with the DNS data of Lee and Moser (2015) [74]. For completeness, also shown is the related experimental data of Kreplin and Eckelmann (1979) [76]. Evidently, good agreement is seen, which demonstrates that the FPC-LBM is well-suited for simulating turbulent flows accurately and effectively.

5. Summary and Conclusions

Collision models play a critical role in determining the numerical accuracy and stability of simulations using LBM. Fokker-Planck (FP) kinetic equation, which represents stochastic processes such as the prototypical Brownian motion, can be used as a model for the collision integral of the Boltzmann equation and involves variations in the distribution function under collision due to its drift and diffusion in the phase space. In this paper, we have derived a new approach based on a FP-guided central moment collision operator for the LBM. It effectively involves relaxation of different central moments to their respective attractors or “equilibria” that depend on the products of lower order central moments and the components of the diffusion coefficient tensor; the latter is related to the variance of the distribution function or its second order central moments. We designate such attractors as the Markovian central moment attractors reflecting the repeated randomness nature of the collision process.

We have constructed central moment LB algorithms based on the FP collision model, which is referred to as the FPC-LBM, in 2D and 3D using D2Q9 and D3Q27 lattices, respectively, through a matching principle, viz., by matching the changes in different discrete central moments independently supported by the lattice under collision to those given by the continuous Boltzmann equation with its collision term expressed in terms of the FP model. We have shown the consistency of our approach to the Navier-Stokes equations via a Chapman-Enskog analysis based directly on expansions of central moments about their attractors. The development of the continuous central moment based FP collision model and their use in the construction of novel LB schemes based on discrete velocities highlighted the important role of the rates of diffusion of the distribution function along different directions in the velocity space, or the second order central moments, in determining the evolution of still higher order moments, which can in turn influence their overall numerical properties.

Simulations of a variety of benchmark problems in both 2D and 3D established its accuracy and demonstrated its superior numerical properties such as enhanced numerical stability and an ability to avoid numerical hyperviscosities in simulations with extremely low physical viscosities of the fluid. It is shown to have cumulant LBM-like behavior in terms of stability while being simpler, and does not involve ad hoc factorized collision formulation to reduce hyperviscosity effects in central moment LBMs as it is rooted and derived directly from a well-founded kinetic model based on the Fokker-Planck formulation. The FPC-LBM is also shown to be effective in computing wall-bounded turbulent flows as it is able to predict turbulence statistics accurately that is in very good agreement when compared to state-of-the-art direct numerical simulations (DNS) data based on a NS-based solver.